text:
######################################################################
# Distributions Package MultivariateDistribution
######################################################################

#################### Dirichlet ####################

Dirichlet(alpha::AbstractVector{T}) where {T<:Real} =
  Dirichlet(convert(Vector{Float64}, alpha))

#################### Multinomial ####################

Multinomial(n::Real, p::AbstractVector{T}) where {T<:Real} =
  Multinomial(convert(Int, n), convert(Vector{Float64}, p))

#################### MvNormal ####################

MvNormal(μ::AbstractVector{T}, Σ::AbstractMatrix{U}) where {T<:Real, U<:Real} =
  MvNormal(convert(Vector{Float64}, μ), convert(Matrix{Float64}, Σ))

MvNormal(μ::AbstractVector{T}, σ::AbstractVector{U}) where {T<:Real, U<:Real} =
  MvNormal(convert(Vector{Float64}, μ), convert(Vector{Float64}, σ))

MvNormal(Σ::AbstractMatrix{T}) where {T<:Real} =
  MvNormal(convert(Matrix{Float64}, Σ))

MvNormal(σ::AbstractVector{T}) where {T<:Real} =
  MvNormal(convert(Vector{Float64}, σ))

#################### MvNormalCanon ####################

MvNormalCanon(h::AbstractVector{T}, J::U) where {T<:Real, U<:AbstractMatrix} =
  MvNormalCanon(convert(Vector{Float64}, h), convert(Matrix{Float64}, J))

MvNormalCanon(h::AbstractVector{T}, prec::U) where {T<:Real, U<:AbstractVector} =
  MvNormalCanon(convert(Vector{Float64}, h), convert(Vector{Float64}, prec))

MvNormalCanon(J::AbstractMatrix{T}) where {T<:Real} =
  MvNormalCanon(convert(Matrix{Float64}, J))

MvNormalCanon(prec::AbstractVector{T}) where {T<:Real} =
  MvNormalCanon(convert(Vector{Float64}, prec))

MvNormalCanon(d::Real, prec::T) where {T<:Real} =
  MvNormalCanon(convert(Int, d), convert(Float64, prec))

#################### MvTDist ####################

MvTDist(df::Real, μ::AbstractVector{T}, C::PDMat) where {T<:Real} =
  MvTDist(convert(Float64, df), convert(Vector{Float64}, μ), C)

MvTDist(df::Real, μ::AbstractVector{T}, Σ::AbstractMatrix{U}) where {T<:Real, U<:Real} =
  MvTDist(convert(Float64, df), convert(Vector{Float64}, μ),
          convert(Matrix{Float64}, Σ))

MvTDist(df::Real, Σ::AbstractMatrix{T}) where {T<:Real} =
  MvTDist(convert(Float64, df), convert(Matrix{Float64}, Σ))

#################### VonMisesFisher ####################

VonMisesFisher(μ::AbstractVector{T}, κ::Real) where {T<:Real} =
  VonMisesFisher(convert(Vector{Float64}, μ), convert(Float64, κ))

######################################################################
# Distributions Package MatrixDistribution
######################################################################

#################### InverseWishart ####################

InverseWishart(df::Real, Ψ::AbstractMatrix{T}) where {T<:Real} =
  InverseWishart(convert(Float64, df), convert(Matrix{Float64}, Ψ))

#################### Wishart ####################

Wishart(df::Real, S::AbstractMatrix{T}) where {T<:Real} =
  Wishart(convert(Float64, df), convert(Matrix{Float64}, S))


text:
using OrdinaryDiffEq
using Trixi
using LinearAlgebra

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 1.4
equations = CompressibleEulerEquations3D(gamma)

# Initial condition for an idealized baroclinic instability test
# https://doi.org/10.1002/qj.2241, Section 3.2 and Appendix A
function initial_condition_pressure_gradient(x, t, equations::CompressibleEulerEquations3D)
  r = x[3]

  radius_earth = 6.371229e6
  # Make sure that r is not smaller than radius_earth
  z = max(r - radius_earth, 0.0)

  # Unperturbed basic state
  rho, u, p = basic_state_baroclinic_instability_longitudinal_velocity(0.0, 0.0, z)

  # No velocity
  v1 = 0.0
  v2 = 0.0
  v3 = 0.0

  return prim2cons(SVector(rho, v1, v2, v3, p), equations)
end

# Unperturbed balanced steady-state.
# Returns primitive variables with only the velocity in longitudinal direction (rho, u, p).
# The other velocity components are zero.
function basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z)
  # Parameters from Table 1 in the paper
  # Corresponding names in the paper are commented
  radius_earth                   = 6.371229e6  # a
  half_width_parameter           = 2           # b
  gravitational_acceleration     = 9.80616     # g
  k                              = 3           # k
  surface_pressure               = 1e5         # p₀
  gas_constant                   = 287         # R
  surface_equatorial_temperature = 310.0       # T₀ᴱ
  surface_polar_temperature      = 240.0       # T₀ᴾ
  lapse_rate                     = 0.005       # Γ
  angular_velocity               = 7.29212e-5  # Ω

  # Distance to the center of the Earth
  r = z + radius_earth

  # In the paper: T₀
  temperature0 = 0.5 * (surface_equatorial_temperature + surface_polar_temperature)
  # In the paper: A, B, C, H
  const_a = 1 / lapse_rate
  const_b = (temperature0 - surface_polar_temperature) /
            (temperature0 * surface_polar_temperature)
  const_c = 0.5 * (k + 2) * (surface_equatorial_temperature - surface_polar_temperature) /
                            (surface_equatorial_temperature * surface_polar_temperature)
  const_h = gas_constant * temperature0 / gravitational_acceleration

  # In the paper: (r - a) / bH
  scaled_z = z / (half_width_parameter * const_h)

  # Temporary variables
  temp1 = exp(lapse_rate/temperature0 * z)
  temp2 = exp(-scaled_z^2)

  # In the paper: ̃τ₁, ̃τ₂
  tau1 = const_a * lapse_rate / temperature0 * temp1 + const_b * (1 - 2 * scaled_z^2) * temp2
  tau2 = const_c * (1 - 2 * scaled_z^2) * temp2

  # In the paper: ∫τ₁(r') dr', ∫τ₂(r') dr'
  inttau1 = const_a * (temp1 - 1) + const_b * z * temp2
  inttau2 = const_c * z * temp2

  # Temporary variables
  temp3 = r/radius_earth * cos(lat)
  temp4 = temp3^k - k/(k + 2) * temp3^(k+2)

  # In the paper: T
  temperature = 1 / ((r/radius_earth)^2 * (tau1 - tau2 * temp4))

  # In the paper: U, u (zonal wind, first component of spherical velocity)
  big_u = gravitational_acceleration/radius_earth * k * temperature * inttau2 * (temp3^(k-1) - temp3^(k+1))
  temp5 = radius_earth * cos(lat)
  u = -angular_velocity * temp5 + sqrt(angular_velocity^2 * temp5^2 + temp5 * big_u)

  # Hydrostatic pressure
  p = surface_pressure * exp(-gravitational_acceleration/gas_constant * (inttau1 - inttau2 * temp4))

  # Density (via ideal gas law)
  rho = p / (gas_constant * temperature)

  return rho, u, p
end

@inline function source_terms_gravity(u, x, t, equations::CompressibleEulerEquations3D)
  radius_earth               = 6.371229e6  # a
  gravitational_acceleration = 9.80616     # g
  angular_velocity           = 7.29212e-5  # Ω

  r = x[3]
  # Make sure that r is not smaller than radius_earth
  z = max(r - radius_earth, 0.0)
  r = z + radius_earth

  du1 = zero(eltype(u))

  # Gravity term
  temp = -gravitational_acceleration * radius_earth^2 / r^3
  du2 = zero(eltype(u))
  du3 = zero(eltype(u))
  du4 = temp * u[1] * x[3]
  du5 = temp * u[4] * x[3]

  return SVector(du1, du2, du3, du4, du5)
end

function indicator_test(u::AbstractArray{<:Any,5},
                        mesh, equations, dg::DGSEM, cache;
                        kwargs...)
  alpha = zeros(Int, nelements(dg, cache))
  radius_earth = 6.371229e6

  for element in eachelement(dg, cache)
    # Calculate coordinates at Gauss-Lobatto nodes
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, i, j, k, element)
      lambda = x[1] / radius_earth

      if pi/16 + 1e-3 < lambda < 3*pi/16 - 1e-3
        alpha[element] = 1
      end
    end
  end

  return alpha
end

function Trixi.get_element_variables!(element_variables, indicator::typeof(indicator_test), ::AMRCallback)
  return nothing
end

initial_condition = initial_condition_pressure_gradient

boundary_conditions = Dict(
  :z_neg  => boundary_condition_slip_wall,
  :z_pos  => boundary_condition_slip_wall
)

# This is a good estimate for the speed of sound in this example.
# Other values between 300 and 400 should work as well.
surface_flux = FluxLMARS(340)
volume_flux  = flux_kennedy_gruber
solver = DGSEM(polydeg=3, surface_flux=surface_flux, volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

radius_earth = 6.371229e6
circumference_earth = 2 * pi * radius_earth

trees_per_dimension = (4, 1, 4)
mesh = P4estMesh(trees_per_dimension, polydeg=1,
                 coordinates_min=(0.0, -1/64 * circumference_earth, radius_earth),
                 coordinates_max=(circumference_earth / 8, 1/64 * circumference_earth, radius_earth + 30000.0),
                 periodicity=(true, true, false), initial_refinement_level=0)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    source_terms=source_terms_gravity,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 5.0e4)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=1000,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

# Abuse the variable `initial_refinement_level` for this because this will be increased
# during convergence tests
initial_refinement_level = 0
amr_controller = ControllerThreeLevel(semi, indicator_test,
                                      base_level=initial_refinement_level,
                                      max_level=initial_refinement_level + 1, max_threshold=0.6)
amr_callback = AMRCallback(semi, amr_controller,
                           interval=0,
                           adapt_initial_condition=true,
                           adapt_initial_condition_only_refine=true)

# This variable can be changed by passing the keyword argument conforming=true
# to trixi_include or convergence_test
conforming=false
if conforming
  callbacks = CallbackSet(summary_callback,
                          analysis_callback,
                          save_solution,
                          alive_callback)
else
  callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        amr_callback,
                        save_solution,
                        alive_callback)
end

###############################################################################
# run the simulation

# Use a Runge-Kutta method with automatic (error based) time step size control
# Enable threading of the RK method for better performance on multiple threads
sol = solve(ode, RDPK3SpFSAL49(thread=OrdinaryDiffEq.True()), abstol=1.0e-6, reltol=1.0e-6,
            save_everystep=false, callback=callbacks);

summary_callback() # print the timer summary

ext: jl


text:
using SIMD
using Base.Cartesian: @nexprs

struct Block{W,MR,NR,G,TA,TB,TAB,TC}
    Ac::TA
    Bc::TB
    AB::TAB
    C::TC
    mc::Int
    kc::Int
    nc::Int
    inc1A::Int
    inc2A::Int
    inc1B::Int
    inc2B::Int
    inc1C::Int
    inc2C::Int
end

const Ac = Vector{UInt8}(undef, 110592)
const Bc = Vector{UInt8}(undef, 6266880)
const AB = Vector{UInt8}(undef, 1 << 12)

function Block(A::X, B::Y, C::Z, ::Val{Width}, ::Val{MR}, ::Val{NR}, ::Val{Generic}) where {X,Y,Z,Width,MR,NR,Generic}
    global Ac, Bc, AB
    m, n = size(C)
    mc = 72
    kc = 192
    nc = 4080
    T = promote_type(eltype(X), eltype(Y), eltype(Z))
    siz = sizeof(T)
    _Ac = unsafe_wrap(Array, Ptr{T}(pointer(Ac)), length(Ac)÷siz)
    _Bc = unsafe_wrap(Array, Ptr{T}(pointer(Bc)), length(Bc)÷siz)
    _AB = unsafe_wrap(Array, Ptr{T}(pointer(AB)), length(AB)÷siz)
    Block{Width,MR,NR,Generic,typeof(_Ac),typeof(_Bc),typeof(_AB),typeof(C)}(_Ac, _Bc, _AB, C, mc, kc, nc, strides(A)..., strides(B)..., strides(C)...,)
end
getmr(blk::Block{W,MR,NR}) where {W,MR,NR} = MR
getnr(blk::Block{W,MR,NR}) where {W,MR,NR} = NR

"""
    addmul!(C, A, B, blk::Block=Block(A, B, C, false)) -> C

`addmul!` computes ``C = AB + C``, where ``A``, ``B``, and ``C`` are matrices.
"""
function addmul!(C, A, B; width=Val(4), mr=Val(8), nr=Val(6), generic=Val(false))
    blk = Block(A, B, C, width, mr, nr, generic)
    m,  k = size(A); _k, n = size(B)
    @assert k == _k
    _m, _n = size(C)
    @assert m == _m && n == _n
    mb, _mc = cld(m, blk.mc), m % blk.mc
    nb, _nc = cld(n, blk.nc), n % blk.nc
    kb, _kc = cld(k, blk.kc), k % blk.kc
    for j in 1:nb # Loop 5
        nc = (j!=nb || _nc==0) ? blk.nc : _nc
        for l in 1:kb # Loop 4
            kc = (l!=kb || _kc==0) ? blk.kc : _kc
            #_β = l==1 ? β : 1.0
            offsetB = blk.kc*(l-1)*blk.inc1B + blk.nc*(j-1)*blk.inc2B
            pack_B!(blk, B, kc, nc, offsetB)
            for i in 1:mb # Loop 3
                mc = (i!=mb || _mc==0) ? blk.mc : _mc
                offsetA = blk.mc*(i-1)*blk.inc1A + blk.kc*(l-1)*blk.inc2A
                offsetC = blk.mc*(i-1)*blk.inc1C + blk.nc*(j-1)*blk.inc2C
                pack_A!(blk, A, mc, kc, offsetA)
                macro_ker!(blk, C, mc, nc, kc, offsetC)
            end # Loop 3
        end # Loop 4
    end # Loop 5
    C
end

@inline function pack_MRxK!(blk::Block, A, k::Int, offsetA::Int, offsetAc::Int)
    mr = getmr(blk)
    @inbounds for j in 1:k
        for i in 1:mr
            blk.Ac[offsetAc+i] = A[offsetA + (i-1)*blk.inc1A + 1]
        end
        offsetAc += mr
        offsetA  += blk.inc2A
    end
    return nothing
end

function pack_A!(blk::Block, A, mc::Int, kc::Int, offsetA::Int)
    mr = getmr(blk)
    mp, _mr = divrem(mc, mr)
    offsetAc = 0
    for i in 1:mp
        pack_MRxK!(blk, A, kc, offsetA, offsetAc)
        offsetAc += kc*mr
        offsetA  += mr*blk.inc1A
    end
    if _mr > 0
        @inbounds for j in 1:kc
            for i in 1:_mr
                blk.Ac[offsetAc+i] = A[offsetA + (i-1)*blk.inc1A + 1]
            end
            for i in _mr+1:mr
                blk.Ac[offsetAc+i] = zero(eltype(A))
            end
            offsetAc += mr
            offsetA  += blk.inc2A
        end
    end
    return nothing
end

@inline function pack_KxNR!(blk::Block, B, k::Int, offsetB::Int, offsetBc::Int)
    nr = getnr(blk)
    @inbounds for i = 1:k
        for j = 1:nr
            blk.Bc[offsetBc+j] = B[offsetB + (j-1)*blk.inc2B + 1]
        end
        offsetBc += nr
        offsetB  += blk.inc1B
    end
    return nothing
end

function pack_B!(blk::Block, B, kc::Int, nc::Int, offsetB::Int)
    nr = getnr(blk)
    np, _nr = divrem(nc, nr)
    offsetBc = 0
    for j in 1:np
        pack_KxNR!(blk, B, kc, offsetB, offsetBc)
        offsetBc += kc*nr
        offsetB  += nr*blk.inc2B
    end
    if _nr > 0
        @inbounds for i in 1:kc
            for j in 1:_nr
                blk.Bc[offsetBc+j] = B[offsetB + (j-1)*blk.inc2B + 1]
            end
            for j in _nr+1:nr
                blk.Bc[offsetBc+j] = zero(eltype(B))
            end
            offsetBc += nr
            offsetB  += blk.inc1B
        end
    end
    return nothing
end

@inline function macro_ker!(blk::Block, C, mc::Int, nc::Int, kc::Int, offsetC::Int)
    mr′, nr′ = getmr(blk), getnr(blk)
    mp, _mr = cld(mc, mr′), mc % mr′
    np, _nr = cld(nc, nr′), nc % nr′
    for j in 1:np
        nr = (j!=np || _nr==0) ? nr′ : _nr
        for i in 1:mp
            mr = (i!=mp || _mr==0) ? mr′ : _mr
            offsetA = (i-1)*kc*mr′
            offsetB = (j-1)*kc*nr′
            if mr == mr′ && nr == nr′
                micro_ker!(blk, kc, offsetA, offsetB, offsetC+(i-1)*mr′*blk.inc1C + (j-1)*nr′*blk.inc2C, Val(true))
            else
                micro_ker!(blk, kc, offsetA, offsetB, 0, Val(false))
                _axpy!(C, 1, blk.AB, mr, nr, offsetC+(i-1)*mr′*blk.inc1C + (j-1)*nr′*blk.inc2C,
                       0, 1, mr′)
            end
        end
    end
    return nothing
end

# generic
@inline function micro_ker!(blk::Block{<:Any,<:Any,<:Any,true}, kc::Int, offsetA::Int, offsetB::Int, offsetC::Int, ::Val{loadC}) where loadC
    mr, nr = getmr(blk), getnr(blk)
    @inbounds begin
        fill!(blk.AB, zero(eltype(blk.AB)))
        for k in 1:kc
            for j in 1:nr, i in 1:mr
                blk.AB[i + (j-1)*mr] += blk.Ac[offsetA+i] * blk.Bc[offsetB+j]
            end
            offsetA += mr
            offsetB += nr
        end
        if loadC
            for j in 1:nr, i in 1:mr
                blk.C[offsetC+(i-1)*blk.inc1C+(j-1)*blk.inc2C+1] += blk.AB[i + (j-1)*mr]
            end
        end
    end
    return nothing
end

# SIMD
@inline @generated function micro_ker!(blk::Block{Width,MR,NR,false,T1}, kc::Int, # TODO: mixed-precision
                            offsetA::Int, offsetB::Int, offsetC::Int,
                            ::Val{loadC}) where {Width,MR,NR,Generic,T1,loadC}
    mr = MR÷Width
    quote
        @inbounds begin
            pA, pAB = pointer(blk.Ac), pointer(blk.AB)
            pB = pointer(blk.Bc)
            T = eltype(T1)
            siz = sizeof(T)
            VT = Vec{$Width, T}
            pC = pointer(blk.C)
            @nexprs $NR i ->
                @nexprs $mr j ->
                    ab_i_j = loadC ? vload(VT, pC + (offsetC+(i-1)*blk.inc2C+(j-1)*$Width)siz) :
                                     zero(VT)
            for k in 1:kc
                # TODO: prefetching
                #_prefetch_r(pA + (offsetA+$MR*(k+2))siz)
                #_prefetch_r(pB + (offsetB+(k+2)*$NR)siz)
                #_prefetch_r(pA + (offsetA+$MR*(k-1))siz + 512)
                #_prefetch_r(pB + (offsetB+$NR*(k-1))siz + 512)
                @nexprs $mr j -> a_j = vload(VT, pA + (offsetA+$MR*(k-1)+(j-1)*$Width)siz)
                @nexprs $NR i -> begin
                    b_i = VT(blk.Bc[offsetB+(k-1)*$NR+i])
                    @nexprs $mr j -> ab_i_j = muladd(a_j, b_i, ab_i_j)
                end
            end
            #loadC && _prefetch_r(pC + offsetC*siz + 512)
            @nexprs $NR i ->
                @nexprs $mr j -> begin
                    ptr = loadC ? pC + (offsetC+(i-1)*blk.inc2C+(j-1)*$Width)siz :
                                  pAB + ((i-1)*$MR + (j-1)*$Width)*siz
                    vstore(ab_i_j, ptr)
                end
        end
        return nothing
    end #quote
end

@inline function _axpy!(Y, α, X, m::Int, n::Int,
                        offsetY::Int, offsetX::Int, inc1X::Int, inc2X::Int)
    inc1Y, inc2Y = stride(Y, 1), stride(Y, 2)
    @inbounds for j in 1:n, i in 1:m
        Y[offsetY+(i-1)*inc1Y+(j-1)*inc2Y+1] += α*X[offsetX+(i-1)*inc1X+(j-1)*inc2X+1]
    end
    return nothing
end


text:
function test_types(::Type{T}) where T <: CasadiSymbolicObject
    @testset "$( string("Console output for ", T, "                            ") )" begin
        @test ( @capture_out show( T(1) ) ) == "1"
        @test ( @capture_out show( T([1 ; 2.1]) ) ) == "[1, 2.1]"
        @test ( @capture_out show( T([1 0 ; 2.1 -3]) ) ) == "\n[[1, 0], \n [2.1, -3]]"
        @test ( @capture_out show( T("x") ) ) == "x"
    end
end


text:
module ScientificTypesBase

# Type exports
export Convention

# re-export-able types and methods
export Scientific, Found, Unknown, Known, Finite, Infinite,
    OrderedFactor, Multiclass, Count, Continuous,
    ScientificTimeType, ScientificDate, ScientificDateTime,
    ScientificTime,
    Textual, Binary,
    ColorImage, GrayImage, Image, Table,
    Density, Sampleable,
    ManifoldPoint,
    Annotated, AnnotationOf, Multiset, Iterator,
    Compositional

export elscitype, nonmissing
# -------------------------------------------------------------------
# Scientific Types

abstract type Found          end
abstract type Known <: Found end
abstract type Unknown <: Found end

abstract type Annotated{S} <: Known end
abstract type AnnotationOf{S} <: Known end
abstract type Multiset{S} <: Known end

# for iterators over objects with scitype Ω that do not have some
# AbstractVector scitype:
abstract type Iterator{Ω} end

abstract type Infinite <: Known end
abstract type Finite{N} <: Known end
abstract type Image{W,H} <: Known end
abstract type ScientificTimeType <: Known end
abstract type Textual <: Known end
abstract type Table{K} <: Known end

abstract type Continuous <: Infinite end
abstract type Count <: Infinite end

abstract type Multiclass{N} <: Finite{N} end
abstract type OrderedFactor{N} <: Finite{N} end

abstract type ScientificDate <: ScientificTimeType end
abstract type ScientificTime <: ScientificTimeType end
abstract type ScientificDateTime <: ScientificTimeType end

abstract type GrayImage{W,H} <: Image{W,H} end
abstract type ColorImage{W,H} <: Image{W,H} end

# when sampled, objects with these scitypes return objects of scitype Ω:
abstract type Sampleable{Ω} end
abstract type Density{Ω} <: Sampleable{Ω} end

abstract type ManifoldPoint{M} <: Known end

# compositional data with D components, see CoDa.jl
abstract type Compositional{D} <: Known end

# aliases:
const Binary = Finite{2}
const Scientific = Union{Missing,Found} # deprecated (no longer publicized)

# for internal use
const Arr = AbstractArray
const TUPLE_SPECIALIZATION_THRESHOLD = 20
# ------------------------------------------------------------------
# Convention

abstract type Convention end

# -----------------------------------------------------------------
# nonmissing

if VERSION < v"1.3"
    # see also discourse.julialang.org/t/get-non-missing-type-in-the-case-of-parametric-type/29109
    """
        nonmissingtype(TT)

    Return the type `T` if `TT = Union{Missing,T}` for some `T` and return `TT`
    otherwise.
    """
    function nonmissingtype(::Type{T}) where T
        return T isa Union ? ifelse(T.a == Missing, T.b, T.a) : T
    end
end

const nonmissing = nonmissingtype

# -----------------------------------------------------------------
# Constructor for table scientific type

"""
Table(...)

Constructor for the `Table` scientific type with:

```
Table(S1, S2, ..., Sn) <: Table
```

where  `S1, ..., Sn` are the scientific type of the table's columns which
are expected to be represented by abstract vectors.
"""
Table(Ts::Type...) = Table{<:Union{(Arr{<:T,1} for T in Ts)...}}

# -----------------------------------------------------------------
# scitype

include("scitype.jl")

end # module


text:
################################################################################
#
# Full batch training using the nlopt.jl package
#
################################################################################

# Outer wrapper into NLopt - This is the function called by the user
function train_nlopt(net::NeuralNet; maxiter=100, algorithm = :LD_LBFGS)

    # Create vector of initial parameters from net
    params = zeros(Float64, num_params(net))
    unroll_param_vec!(params, net)

    # Train the neural net
    train_nlopt(net, params, maxiter=maxiter, algorithm=algorithm)
end

# Wrapper into NLopt
function train_nlopt(net::NeuralNet, params::Vector; maxiter=100, algorithm = :LD_LBFGS)

    # Create "short hand" for the cost function to be minimized
    function nlopt_cost_function(params, grad)
        cost_and_grad!(net, params, grad)
    end

    # Options
    options = Opt(algorithm, length(params))
    maxeval!(options, maxiter)

    # Objective
    min_objective!(options, nlopt_cost_function)

    # Train the neural net
    @printf "Training Loss:\n"
    cost, params, _ = optimize(options, params)
end

# Single function to return cost and gradient - input into convex optimiser
function cost_and_grad!{F<:AbstractFloat}(net::NeuralNet, params::Vector{F}, grad::Vector{F})

    # Copy params into net
    input_param_vec!(params, net)

    # Perform forward propogation through the layers
    fwd_prop!(net)

    # Cost / Loss / Error
    L = net.layers[end]
    J = L.cost(L.A, L.y) / L.m + net.reg_cost(net) + sparse_cost(net)
    @printf " -> %.3f" J

    # Perform back propogation to calculate gradients
    bwd_prop!(net)

    # Unroll and set gradient vector in-place
    unroll_grad_vec!(grad, net)

    # Return cost
    J
end

# Roll single vector of parameters into the network
function input_param_vec!(params, net)
    @assert length(params) == num_params(net)
    start = 1
    for ℓ in net.layers[2:end]

        # No parameters in pooling layers
        typeof(ℓ) <: PoolLayer && continue

        # Place bias
        idx = length(ℓ.b)
        stop = start + idx - 1
        ℓ.b[:] = view(params, start:stop)
        start += idx

        # Place weight
        idx = length(ℓ.W)
        stop = start + idx - 1
        ℓ.W[:] = view(params, start:stop)
        start += idx
    end
end

# Unroll parameters from the network into a single vector
function unroll_param_vec!(params, net)
    @assert length(params) == num_params(net)
    start = 1
    for ℓ in net.layers[2:end]

        # No parameters in pooling layers
        typeof(ℓ) <: PoolLayer && continue

        # Place bias
        idx = length(ℓ.b)
        stop = start + idx - 1
        params[start:stop] = ℓ.b[:]
        start += idx

        # Place weight
        idx = length(ℓ.W)
        stop = start + idx - 1
        params[start:stop] = ℓ.W[:]
        start += idx
    end
end

# Unroll gradient parameters from the network into a single vector
function unroll_grad_vec!(grad, net)
    @assert length(grad) == num_params(net)
    start = 1
    for ℓ in net.layers[2:end]

        # No parameters in pooling layers
        typeof(ℓ) <: PoolLayer && continue

        # Place bias
        idx = length(ℓ.∇b)
        stop = start + idx - 1
        grad[start:stop] = ℓ.∇b[:]
        start += idx

        # Place weight
        idx = length(ℓ.∇W)
        stop = start + idx - 1
        grad[start:stop] = ℓ.∇W[:]
        start += idx
    end
end


text:
# Test collection_domains function
@testset "collection_domains" begin
    domain1 = IntervalDomain(0, 1)
    domain2 = UniDistributionDomain(Uniform())
    domain = CollectionDomain([domain1, domain2])
    @test_throws ErrorException collection_domains(domain1)
    @test collection_domains(domain) == [domain1, domain2]
    @test_deprecated collection_sets(domain)
end

# Test length
@testset "Base.length" begin
    sdomain1 = IntervalDomain(0, 1)
    sdomain2 = UniDistributionDomain(Uniform())
    domain1 = CollectionDomain([sdomain1, sdomain2])
    domain2 = MultiDistributionDomain(MvNormal(ones(4)))
    @test length(sdomain1) == 1
    @test length(sdomain2) == 1
    @test length(domain1) == 2
    @test length(domain2) == 4
end

# Test supports_in_domain functions
@testset "supports_in_domain" begin
    # supports_in_domain (IntervalDomain)
    @testset "IntervalDomain" begin
        domain = IntervalDomain(0, 1)
        @test supports_in_domain(0, domain)
        @test !supports_in_domain(-1, domain)
        @test !supports_in_domain(2, domain)
    end
    # supports_in_domain (UnivariateDistribution domain)
    @testset "Univariate Distribution" begin
        domain = UniDistributionDomain(Uniform())
        @test supports_in_domain(0, domain)
        @test !supports_in_domain(-1, domain)
        @test !supports_in_domain(2, domain)
    end
    # supports_in_domain (MultivariateDistribution domain)
    @testset "Multivariate Distribution" begin
        domain = MultiDistributionDomain(MvNormal(ones(2)))
        @test supports_in_domain(ones(2), domain)
        @test supports_in_domain(ones(2, 10), domain)
        bad_supports = [1 1; 2 2; 3 3];
        @test_throws ErrorException supports_in_domain(bad_supports, domain)
    end
    # supports_in_domain (CollectionDomain)
    @testset "CollectionDomain" begin
        domain1 = IntervalDomain(0, 1)
        domain2 = UniDistributionDomain(Uniform())
        domain = CollectionDomain([domain1, domain2])
        bad_supports = [1 1; 2 2; 3 3];
        @test_throws ErrorException supports_in_domain(bad_supports, domain)
        supports = [1 2; 3 4];
        @test !supports_in_domain(supports, domain)
        supports = [0 1; 0 1];
        @test supports_in_domain(supports, domain)
    end
    # supports_in_domain (Fallback)
    @testset "Fallback" begin
        @test supports_in_domain(0, BadDomain())
        @test_deprecated supports_in_set(0, BadDomain())
    end
end

# Test lower bound functions
@testset "Lower Bound" begin
    domain1 = IntervalDomain(0, 1)
    domain2 = UniDistributionDomain(Uniform())
    domain3 = MultiDistributionDomain(MvNormal([0., 0.], [1. 0.; 0. 2.]))
    domain4 = CollectionDomain([domain1, domain2])
    domain5 = BadDomain()
    # JuMP.has_lower_bound
    @testset "JuMP.has_lower_bound" begin
        @test has_lower_bound(domain1)
        @test has_lower_bound(domain2)
        @test !has_lower_bound(domain3)
        @test has_lower_bound(domain4)
        @test !has_lower_bound(domain5)
    end
    # JuMP.lower_bound
    @testset "JuMP.lower_bound" begin
        @test lower_bound(domain1) == 0.0
        @test lower_bound(domain2) == 0.0
        @test_throws ErrorException lower_bound(domain3)
        @test lower_bound(domain4) == [0.0, 0.0]
        @test_throws ErrorException lower_bound(domain5)
    end
    # JuMP.set_lower_bound
    @testset "JuMP.set_lower_bound" begin
        @test set_lower_bound(domain1, 0.5) == IntervalDomain(0.5, 1)
        @test_throws ErrorException set_lower_bound(domain2, 2)
        @test_throws ErrorException set_lower_bound(domain3, 2)
        @test_throws ErrorException set_lower_bound(domain4, [0.1, 0.1])
        @test_throws ErrorException set_lower_bound(domain5, 2)
        domain6 = CollectionDomain([domain1, domain1])
        @test set_lower_bound(domain6, [0.1, 0.1]).domains ==
                                    [IntervalDomain(0.1, 1.), IntervalDomain(0.1, 1.)]
    end
end

# Test upper bound functions
@testset "Upper Bound" begin
    domain1 = IntervalDomain(0, 1)
    domain2 = UniDistributionDomain(Uniform())
    domain3 = MultiDistributionDomain(MvNormal([0., 0.], [1. 0.; 0. 2.]))
    domain4 = CollectionDomain([domain1, domain2])
    domain5 = BadDomain()
    # JuMP.has_upper_bound
    @testset "JuMP.has_upper_bound" begin
        @test has_upper_bound(domain1)
        @test has_upper_bound(domain2)
        @test !has_upper_bound(domain3)
        @test has_upper_bound(domain4)
        @test !has_upper_bound(domain5)
    end
    # JuMP.upper_bound
    @testset "JuMP.upper_bound" begin
        @test upper_bound(domain1) == 1.0
        @test upper_bound(domain2) == 1.0
        @test_throws ErrorException upper_bound(domain3)
        @test upper_bound(domain4) == [1.0, 1.0]
        @test_throws ErrorException upper_bound(domain5)
    end
    # JuMP.set_upper_bound
    @testset "JuMP.set_upper_bound" begin
        @test set_upper_bound(domain1, 0.5) == IntervalDomain(0, 0.5)
        @test_throws ErrorException set_upper_bound(domain2, 2)
        @test_throws ErrorException set_upper_bound(domain3, 2)
        @test_throws ErrorException set_upper_bound(domain4, [0.1, 0.1])
        @test_throws ErrorException set_upper_bound(domain5, 2)
        domain6 = CollectionDomain([domain1, domain1])
        @test set_upper_bound(domain6, [0.9, 0.9]).domains ==
                                    [IntervalDomain(0., 0.9), IntervalDomain(0., 0.9)]
    end
end

# Test support labels
@testset "Support Labels" begin
    @testset "DataTypes" begin
        @test AbstractSupportLabel isa DataType
        @test All <: AbstractSupportLabel
        @test PublicLabel <: AbstractSupportLabel
        @test UserDefined <: PublicLabel
        @test UniformGrid <: PublicLabel
        @test SampleLabel <: PublicLabel
        @test MCSample <: SampleLabel
        @test WeightedSample <: SampleLabel
        @test Mixture <: PublicLabel
        @test UniqueMeasure isa UnionAll
        @test InternalLabel <: AbstractSupportLabel
        @test InfiniteOpt._NoLabel <: AbstractSupportLabel
    end
    @testset "generate_unique_label" begin
        @test InfiniteOpt.generate_unique_label() <: UniqueMeasure
    end
end

# Test support generation
@testset "generate_support_values" begin
    @testset "IntervalDomain" begin
        domain = IntervalDomain(0., 1.)
        @test generate_support_values(domain, num_supports = 10, sig_digits = 3)[1] isa Vector{<:Number}
        @test generate_support_values(domain, num_supports = 10, sig_digits = 3)[2] == UniformGrid
        @test generate_support_values(domain, num_supports = 10, sig_digits = 3)[1][2] == 0.111
        @test generate_support_values(domain, num_supports = 10, sig_digits = 3)[1][2] != 1/11
        @test length(generate_support_values(domain, num_supports = 10, sig_digits = 3)[1]) == 10
        @test generate_support_values(domain, MCSample, num_supports = 10)[1] isa Vector{<:Number}
        @test generate_support_values(domain, MCSample, num_supports = 10)[2] == MCSample
        @test_throws ErrorException generate_support_values(domain, Val(:a))
        @test generate_support_values(domain, All)[2] == UniformGrid
    end
    @testset "Distribution Domains" begin
        dist1 = Normal(0., 1.)
        dist2 = MvNormal([0.; 0.], [1. 0.; 0. 2.])
        domain1 = UniDistributionDomain(dist1)
        domain2 = MultiDistributionDomain(dist2)
        @test generate_support_values(domain1, num_supports = 10)[1] isa Vector{<:Number}
        @test generate_support_values(domain2, num_supports = 10)[1] isa Array{<:Number, 2}
        @test generate_support_values(domain2, num_supports = 10)[2] == WeightedSample
        @test length(generate_support_values(domain1, num_supports = 10)[1]) == 10
        @test size(generate_support_values(domain2, num_supports = 10)[1]) == (2, 10)
        @test_throws ErrorException generate_support_values(domain1, Val(:a))
        @test_throws ErrorException generate_support_values(domain2, Val(:a))
        @test generate_support_values(domain1, All)[2] == WeightedSample
        @test generate_support_values(domain2, All)[2] == WeightedSample
        @test generate_support_values(domain1, MCSample)[2] == MCSample
    end
    @testset "Matrix Distribution Domains" begin
        dist = MatrixBeta(2, 2, 2)
        domain = MultiDistributionDomain(dist)
        @test generate_support_values(domain, num_supports = 10)[1] isa Array{<:Number, 2}
        @test generate_support_values(domain, num_supports = 10)[2] == WeightedSample
        @test size(generate_support_values(domain, num_supports = 10)[1]) == (4, 10)
        @test_throws ErrorException generate_support_values(domain, Val(:a))
        @test generate_support_values(domain, All)[2] == WeightedSample
    end
    @testset "_generate_collection_supports" begin
        domain1 = IntervalDomain(0., 1.)
        domain2 = IntervalDomain(0., 1.)
        domain = CollectionDomain([domain1, domain2])
        @test InfiniteOpt._generate_collection_supports(domain, 10, 3) isa Array{<:Number, 2}
        @test InfiniteOpt._generate_collection_supports(domain, 10, 3)[2, 2] == 0.111
        @test InfiniteOpt._generate_collection_supports(domain, 10, 3)[2, 2] != 1/11
        @test size(InfiniteOpt._generate_collection_supports(domain, 10, 3)) == (2, 10)
    end
    @testset "CollectionDomain (IntervalDomains)" begin
        domain1 = IntervalDomain(0., 1.)
        domain2 = IntervalDomain(0., 1.)
        domain = CollectionDomain([domain1, domain2])
        @test generate_support_values(domain, num_supports = 10, sig_digits = 3)[1] isa Array{<:Number, 2}
        @test generate_support_values(domain, num_supports = 10, sig_digits = 3)[2] == UniformGrid
        @test generate_support_values(domain, num_supports = 10, sig_digits = 3)[1][2, 2] == 0.111
        @test generate_support_values(domain, num_supports = 10, sig_digits = 3)[1][2, 2] != 1/11
        @test size(generate_support_values(domain, num_supports = 10, sig_digits = 3)[1]) == (2, 10)
        @test_throws ErrorException generate_support_values(domain, Val(:a))
        @test generate_support_values(domain, All)[2] == UniformGrid
        @test generate_support_values(domain, MCSample)[2] == MCSample
        @test generate_support_values(domain, MCSample)[1] isa Array{<:Number, 2}
    end
    @testset "CollectionDomain (UniDistributionDomains)" begin
        domain1 = UniDistributionDomain(Normal())
        domain2 = UniDistributionDomain(Normal())
        domain = CollectionDomain([domain1, domain2])
        @test generate_support_values(domain, num_supports = 10, sig_digits = 3)[1] isa Array{<:Number, 2}
        @test generate_support_values(domain, num_supports = 10, sig_digits = 3)[2] == WeightedSample
        @test size(generate_support_values(domain, num_supports = 10, sig_digits = 3)[1]) == (2, 10)
        @test_throws ErrorException generate_support_values(domain, Val(:a))
        @test generate_support_values(domain, All)[2] == WeightedSample
    end
    @testset "CollectionDomain (InfiniteScalarDomains)" begin
        domain1 = UniDistributionDomain(Normal())
        domain2 = IntervalDomain(0., 1.)
        domain = CollectionDomain([domain1, domain2])
        @test generate_support_values(domain, num_supports = 10, sig_digits = 3)[1] isa Array{<:Number, 2}
        @test generate_support_values(domain, num_supports = 10, sig_digits = 3)[2] == Mixture
        @test generate_support_values(CollectionDomain([domain2, domain2]), MCSample, num_supports = 10, sig_digits = 3)[2] == MCSample
        @test size(generate_support_values(domain, num_supports = 10, sig_digits = 3)[1]) == (2, 10)
        @test_throws ErrorException generate_support_values(domain, Val(:a))
        @test generate_support_values(domain, All)[2] == Mixture
        @test generate_support_values(domain, MCSample)[2] == MCSample
        @test generate_support_values(domain, MCSample)[1] isa Array{<:Number, 2}
    end
    @testset "Fallback" begin
        @test_throws ErrorException generate_support_values(BadDomain())
    end
    @testset "User Interface" begin
        domain = IntervalDomain(0., 1.)
        @test generate_supports(domain)[1] isa Vector{<:Number}
        @test generate_supports(domain, MCSample)[1] isa Vector{<:Number}
    end
end


text:
module EIIndex

using LightGraphs
using LightGraphs.SimpleGraphs
using Plots: histogram, vline!
import Plots: plot

export ei_global, ei_group, ei_indiv, cugtest_edges, plot

include("ei_index.jl")
include("cugtest.jl")

end#module

max_stars_repo_name: jfb-h/EIIndex

max_issues_repo_name: jfb-h/EIIndex

max_forks_repo_name: jfb-h/EIIndex


text:
module SparseCCSA
    using CatViews
    using LinearAlgebra
    export CCSAState
    export optimize
    include("ccsa.jl")
end


text:
export NestedEinsum

"""
    parse_nested(s::AbstractString (, iy = []))

return a contraction-tree consisting of `NestedEinsumConstructor`
objects as nodes and `IndexGroup` as leaves.
`s` should only contain the `ixs` part
"""

function parse_nested(s::AbstractString, iy = [])
    count(==('('),s) == count(==(')'),s) || throw(
        ArgumentError("Parentheses don't pair up in $s"))

    _, out = parse_parens(s, firstindex(s), 1)
    append!(out.iy, iy)
    filliys!(out)
    return construct(out)
end

"""
    parse_parens(s::AbstractString, i, narg)

parse one level of parens starting at index `i` where `narg` counts which tensor the
current group of indices, e.g. "ijk", belongs to.
Recursively calls itself for each new opening paren that's opened.
"""
function parse_parens(s::AbstractString, i, narg)
    out = NestedEinsumConstructor{Char}([], [], [])
    g = IndexGroup{Char}([],narg)

    while i <= lastindex(s)
        c = s[i]
        j = nextind(s,i)
        if c === '('
            # opening a parens means that next a new contraction is parsed
            j, out2, narg = parse_parens(s, j, narg)
            push!(out, out2)
            union!(out.inds, out2.inds)
        elseif c === ')' || c === ','
            # either closing a parens or a comma means that a tensors indices
            # are complete and can be added to the contraction-object
            if !isempty(g)
                push!(out, g)
                union!(out.inds, g.inds)
            end
            if  c === ','
                # comma implies that a new tensor is parsed for the current contraction
                narg += 1
                g = IndexGroup{Char}([], narg)
            else
                # parens implies that the current contraction is complete
                return j, out, narg
            end
        elseif isletter(c) # could allow isnumeric() too
            push!(g,c)
        else
            throw(ArgumentError("parsing $s failed, $c is not a valid entry"))
        end
        i = j
    end
    if !isempty(g)
        push!(out, g)
        union!(out.inds, g.inds)
    end
    return i, out
end

"""
    filliys!(neinsum::NestedEinsumConstructor)

goes through all `NestedEinsumConstructor` objects in the tree and saves the correct `iy` in them.
"""
function filliys!(neinsum)
    iy = neinsum.iy
    args = neinsum.args
    for i in 1:length(neinsum.args)
        arg = args[i]
        arg isa IndexGroup && continue
        union!(arg.iy, intersect(arg.inds, iy))
        for j in 1:length(neinsum.args)
            i === j && continue
            union!(arg.iy, intersect(arg.inds, args[j].inds))
        end
        filliys!(arg)
    end
    return neinsum
end

"""
    IndexGroup

Leaf in a contractiontree, contains the indices and the number of the tensor it
describes, e.g. in "ij,jk -> ik", indices "ik" belong to tensor `1`, so
would be described by IndexGroup(['i','k'], 1).
"""
struct IndexGroup{T}
    inds::Vector{T}
    n::Int
end

Base.push!(ig::IndexGroup, c) = (push!(ig.inds,c); ig)
Base.isempty(ig::IndexGroup) = isempty(ig.inds)

"""
    NestedEinsumConstructor

describes a (potentially) nested einsum. Important fields:
- `args`, vector of all inputs, either `IndexGroup` objects corresponding to tensors or `NestedEinsumConstructor`
- `iy`, indices of output
"""
struct NestedEinsumConstructor{T}
    args::Vector{Union{NestedEinsumConstructor{T}, IndexGroup{T}}}
    inds::Vector{T}
    iy::Vector{T}
end

Base.push!(neinsum::NestedEinsumConstructor, x) = (push!(neinsum.args,x); neinsum)

using MacroTools

function _nested_ein_macro(ex; einsum=:einsum)
    @capture(ex, (left_ := right_)) || throw(ArgumentError("expected A[] := B[]... "))
    @capture(left, Z_[leftind__] | [leftind__] ) || throw(
        ArgumentError("can't understand LHS, expected A[i,j] etc."))
    Z === nothing && @gensym Z
    primefix!(leftind)
    allinds = unique(leftind)

    MacroTools.postwalk(right) do x
        @capture(x, A_[inds__]) && union!(allinds, inds)
        x
    end
    primefix!(allinds)

    tensors = Symbol[]
    nein = parse_nested_expr(right, tensors, allinds)
    append!(nein.iy, indexin(leftind,allinds))
    filliys!(nein)
    snein = construct(nein)

    tensornames = map(esc, tensors)
    :($(esc(Z)) = $snein(($(tensornames...),)...))
end

function parse_nested_expr(expr, tensors, allinds)
    if @capture(expr, *(args__))
        einargs = map(x -> parse_nested_expr(x,tensors, allinds), args)
        intinds = union(mapreduce(x -> x.inds, vcat, einargs))
        return NestedEinsumConstructor{Int}(einargs, intinds, Int[])
    elseif @capture(expr, A_[inds__])
        push!(tensors,A)
        return IndexGroup{Int}(indexin(primefix!(inds), allinds), length(tensors))
    end
end

struct NestedEinsum{ET}
    args::Vector{NestedEinsum{ET}}
    tensorindex::Int  # -1 if not leaf
    eins::ET

    NestedEinsum(args::Vector{NestedEinsum{DynamicEinCode{LT}}}, eins::DynamicEinCode{LT}) where LT = new{DynamicEinCode{LT}}(args, -1, eins)
    NestedEinsum{DynamicEinCode{LT}}(arg::Int) where LT = new(NestedEinsum{DynamicEinCode{LT}}[], arg)
    function NestedEinsum(args::Tuple, eins::DynamicEinCode{LT}) where LT
        @assert length(args) == length(getixs(eins))
        new{DynamicEinCode{LT}}([arg isa Int ? NestedEinsum{DynamicEinCode{LT}}(arg) : arg for arg in args], -1, eins)
    end

    NestedEinsum(args::Vector{NestedEinsum{StaticEinCode}}, eins::StaticEinCode) = new{StaticEinCode}(args, -1, eins)
    NestedEinsum{StaticEinCode}(arg::Int) = new{StaticEinCode}(NestedEinsum{StaticEinCode}[], arg)
    function NestedEinsum(args::Tuple, eins::StaticEinCode)
        @assert length(args) == length(getixs(eins))
        new{StaticEinCode}([arg isa Int ? NestedEinsum{StaticEinCode}(arg) : arg for arg in args], -1, eins)
    end
end

isleaf(ne::NestedEinsum) = ne.tensorindex != -1

function Base.:(==)(a::NestedEinsum, b::NestedEinsum)
    ex = a.args == b.args && a.tensorindex == b.tensorindex
    if isdefined(a, :eins) && isdefined(b, :eins)
        return ex && a.eins == b.eins
    elseif !(isdefined(a, :eins)) && !(isdefined(b, :eins))
        return ex
    else
        return false
    end
end

function construct(nein::NestedEinsumConstructor{T}) where T
    ixs = Tuple(map(extractixs, nein.args))
    iy = Tuple(nein.iy)
    eins = StaticEinCode{ixs,iy}()
    args = Tuple(map(x -> x isa NestedEinsumConstructor ? construct(x) : x.n,nein.args))
    return NestedEinsum(args, eins)
end
extractixs(x::IndexGroup) = Tuple(x.inds)
extractixs(x::NestedEinsumConstructor) = Tuple(x.iy)

function (neinsum::NestedEinsum)(@nospecialize(xs::AbstractArray...); size_info = nothing)
    size_dict = size_info===nothing ? Dict{labeltype(neinsum.eins),Int}() : copy(size_info)
    get_size_dict!(neinsum, xs, size_dict)
    return einsum(neinsum, xs, size_dict)
end

function get_size_dict!(ne::NestedEinsum, @nospecialize(xs), size_info::Dict{LT}) where LT
    d = collect_ixs!(ne, Dict{Int,Vector{LT}}())
    ks = sort!(collect(keys(d)))
    ixs = [d[i] for i in ks]
    return get_size_dict_!(ixs, [collect(Int, size(xs[i])) for i in ks], size_info)
end

collect_ixs(ne::EinCode) = [_collect(ix) for ix in getixs(ne)]
function collect_ixs(ne::NestedEinsum)
    d = OMEinsum.collect_ixs!(ne, Dict{Int,Vector{OMEinsum.labeltype(ne.eins)}}())
    ks = sort!(collect(keys(d)))
    return @inbounds [d[i] for i in ks]
end
function collect_ixs!(ne::NestedEinsum, d::Dict{Int,Vector{LT}}) where LT
    @inbounds for i=1:length(ne.args)
        arg = ne.args[i]
        if isleaf(arg)
            d[arg.tensorindex] = _collect(LT, OMEinsum.getixs(ne.eins)[i])
        else
            collect_ixs!(arg, d)
        end
    end
    return d
end

function einsum(neinsum::NestedEinsum, @nospecialize(xs::NTuple{N,AbstractArray} where N), size_dict::Dict)
    # do not use map because the static overhead is too large
    # do not use `setindex!` because we need to make the AD work
    mxs = Vector{AbstractArray}(undef, length(neinsum.args))
    for (i, arg) in enumerate(neinsum.args)
        mxs = _safe_set(mxs, i, isleaf(arg) ? xs[arg.tensorindex] : einsum(arg, xs, size_dict))
    end
    return einsum(neinsum.eins, (mxs...,), size_dict)
end

_safe_set(lst, i, y) = (lst[i] = y; lst)

# Better printing
using AbstractTrees

function AbstractTrees.children(ne::NestedEinsum)
    d = Dict()
    for (k,item) in enumerate(ne.args)
        d[k] = isleaf(item) ? _join(OMEinsum.getixs(ne.eins)[k]) : item
    end
    d
end

function AbstractTrees.printnode(io::IO, x::String)
    print(io, x)
end
function AbstractTrees.printnode(io::IO, x::NestedEinsum)
    isleaf(x) ? print(io, x.tensorindex) : print(io, x.eins)
end

function Base.show(io::IO, e::EinCode)
    s = join([_join(ix) for ix in getixs(e)], ", ") * " -> " * _join(getiy(e))
    print(io, s)
end
function Base.show(io::IO, e::NestedEinsum)
    print_tree(io, e)
end
Base.show(io::IO, ::MIME"text/plain", e::NestedEinsum) = show(io, e)
Base.show(io::IO, ::MIME"text/plain", e::EinCode) = show(io, e)
_join(ix) = isempty(ix) ? "" : join(ix, connector(eltype(ix)))
connector(::Type{Char}) = ""
connector(::Type{Int}) = "∘"
connector(::Type) = "-"

# flattten nested einsum
function _flatten(code::NestedEinsum, iy=nothing)
    isleaf(code) && return [code.tensorindex=>iy]
    ixs = []
    for i=1:length(code.args)
        append!(ixs, _flatten(code.args[i], OMEinsum.getixsv(code.eins)[i]))
    end
    return ixs
end

flatten(code::EinCode) = code
function flatten(code::NestedEinsum)
    ixd = Dict(_flatten(code))
    if code.eins isa DynamicEinCode
        DynamicEinCode([ixd[i] for i=1:length(ixd)], collect(OMEinsum.getiy(code.eins)))
    else
        StaticEinCode{ntuple(i->(ixd[i]...,), length(ixd)), OMEinsum.getiy(code.eins)}()
    end
end


text:
module test_millboard_unicode

using Millboard
using Test

Millboard.set_table_mode(:markdown)

board = ["한글"]
@test sprint(show, table(board)) == """
|   | 1    |
|---|------|
| 1 | 한글 |"""

board = ["한\n글"; "abc"]
@test sprint(show, table(board)) == """
|   | 1   |
|---|-----|
|   | 한  |
| 1 | 글  |
| 2 | abc |"""

end # module test_millboard_unicode


text:
@testset "Velocity Constraint" begin

    # UnicycleGame
    T = Float64
    N = 10
    dt = 0.1
    model = UnicycleGame(p=3)
    probsize = ProblemSize(N,model)
    pdtraj = PrimalDualTraj(probsize, dt, f=ones, amplitude=0.1)
    game_con = GameConstraintValues(probsize)

    v_max =  ones(model.p)
    v_min = -ones(model.p)
    add_velocity_bound!(model, game_con, v_max, v_min)

    cval = game_con.state_conval[1]
    @test length(game_con.state_conval) == model.p
    @test length(cval) == model.p

    # BicycleGame
    model = BicycleGame(p=3)
    probsize = ProblemSize(N,model)
    game_con = GameConstraintValues(probsize)

    v_max =  [1,+Inf,+Inf]
    v_min =  [1,-1,-Inf]
    add_velocity_bound!(model, game_con, v_max, v_min)

    cval = game_con.state_conval[1]
    @test length(game_con.state_conval) == model.p
    @test length(cval) == model.p-1

    # Empty constraints
    game_con = GameConstraintValues(probsize)

    v_max =  Inf*ones(model.p)
    v_min = -Inf*ones(model.p)
    add_velocity_bound!(model, game_con, v_max, v_min)

    cval = game_con.state_conval[1]
    @test length(game_con.state_conval) == model.p
    @test length(cval) == 0

end


text:
function spiral_matrix(n)

end


text:
#---------------------------------------------------------------------#
#This function computes the Grid and Elements in 1D.
#Written by F.X. Giraldo on April 19, 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function create_grid(Np,Npoin_cg,Npoin_dg,Ne,x,DFloat)

    coord_cg=zeros(DFloat,Npoin_cg)
    coord_dg=zeros(DFloat,Npoin_dg)
    intma_cg=zeros(Int64,Np,Ne)
    intma_dg=zeros(Int64,Np,Ne)
    periodicity_cg=zeros(Int64,Npoin_cg)
    periodicity_dg=zeros(Int64,Npoin_dg)

    #Set some constants
    xmin=-1
    xmax=+1
    dx=(xmax-xmin)/Ne

    #Generate Grid Points
    I=1
    coord_cg[1]=xmin
    for e=1:Ne
        x0=xmin + (e-1)*dx
        intma_cg[1,e]=I
        for i=2:Np
            I+=1
            coord_cg[I]=( x[i]+1 )*dx/2 + x0
            intma_cg[i,e]=I
        end
    end

    #Generate Periodicity_CG array
    for I=1:Npoin_cg
        periodicity_cg[I]=I
    end
    periodicity_cg[Npoin_cg]=1

    #Generate INTMA_DG
    I=0
    for e=1:Ne, i=1:Np
        I+=1
        intma_dg[i,e]=I
    end

    #Generate COORD_DG
    for e=1:Ne, i=1:Np
        Icg=intma_cg[i,e]
        Idg=intma_dg[i,e]
        coord_dg[Idg]=coord_cg[Icg]
    end

    #Generate Periodicity_DG array
    for I=1:Npoin_dg
        periodicity_dg[I]=I
    end

    return (coord_cg, coord_dg, intma_cg, intma_dg, periodicity_cg, periodicity_dg)
end


text:
using SortingAlgorithms

raw_rel = readdlm("C:\\Users\\bernhard.konig\\Documents\\ASync\\irobot\\data_tmp\\eraw_rel.csv");
orig = deepcopy(raw_rel[:]);

r = deepcopy(orig);i = zeros(Int, length(r)); @benchmark mysortperm!(i, r)
r2 = deepcopy(orig);i2 = zeros(Int, length(r2)); @benchmark sortperm!(i2, r2)
@assert isequal(r2[i2], r[i])

isequal(i,i2)

r = deepcopy(orig); a = sortperm(r)
r = deepcopy(orig); a1, b1 = mysortperm(r)

isequal(a,a1)
isequal(raw_rel[a],b1)

r = deepcopy(orig); @btime sort!(r);

r = deepcopy(orig); @btime sort!(r);
r = deepcopy(orig); @btime sort!(r, alg=QuickSort);
r = deepcopy(orig); @btime sort!(r, alg=RadixSort);

r = deepcopy(orig); @btime sortperm(r);
r = deepcopy(orig); @btime mysortperm(r);
r = deepcopy(orig); @btime sortperm(r, alg=QuickSort);
r = deepcopy(orig); @benchmark mysortperm(r)
r = deepcopy(orig); @benchmark sortperm(r, alg=QuickSort)

r = deepcopy(orig);i = zeros(Int, length(r)); @btime sortperm!(i, r, alg=QuickSort);

size(orig)

xx = rand(length(orig));@btime sortperm(xx);
xx = rand(length(orig));@btime sort!(xx);

export mysortperm

function mysortperm!(ii, A)
    n = length(A)
   # ii = Array(Int,n)

    for i = 1:n
        ii[i] = i
    end

    B = copy(A)
    quicksort!(B, ii)

   # return ii, B
    return nothing
end # function mysortperm!

function mysortperm(A)
    n = length(A)
    ii = Array(Int, n)

    for i = 1:n
        ii[i] = i
    end

    B = copy(A)
    quicksort!(B, ii)

    return ii, B
end # function mysortperm

function quicksort!(A, order, i=1, j=length(A))
# modified from:
# http://rosettacode.org/wiki/Sorting_algorithms/Quicksort#Julia
    if j > i
        if j - i <= 50
              # Insertion sort for small groups is faster than Quicksort
            InsertionSort!(A, order, i, j)
            return A
        end

        # pivot = A[rand(i:j)] # random element of A
        pivot = A[ div(i + j, 2) ]
        left, right = i, j
        while left <= right
            while A[left] < pivot
                left += 1
            end
            while A[right] > pivot
                right -= 1
            end
            if left <= right
                A[left], A[right] = A[right], A[left]
                order[left], order[right] = order[right], order[left]

                left += 1
                right -= 1
            end
        end  # left <= right

        quicksort!(A, order, i,   right)
        quicksort!(A, order, left, j)
    end  # j > i

    return A
end # function quicksort!

function InsertionSort!(A, order, ii=1, jj=length(A))

    for i = ii + 1:jj
        j = i - 1
        temp  = A[i]
        itemp = order[i]

        while true
            if j == ii - 1
                break
            end
            if A[j] <= temp
                break
            end
            A[j + 1] = A[j]
            order[j + 1] = order[j]
            j -= 1
        end

        A[j + 1] = temp
        order[j + 1] = itemp
    end  # i

    return
end # function InsertionSort!


text:
using Reexport

@reexport module UncertainValues
    using IntervalArithmetic
    import IntervalArithmetic: interval
    using Distributions
    using RecipesBase
    using StaticArrays
    using Measurements
    import KernelDensity: KernelDensity, UnivariateKDE, default_bandwidth, kde
    import StatsBase: StatsBase,
        ProbabilityWeights, pweights,
        FrequencyWeights, fweights,
        Weights, weights,
        AnalyticWeights, aweights

    ##########################################
    # Abstract types
    ##########################################
    include("abstract_types/AbstractUncertainValue.jl")

    include("abstract_types/AbstractTheoreticalDistributionScalarValue.jl")

    include("abstract_types/AbstractPopulation.jl")
    include("abstract_types/AbstractScalarPopulation.jl")

    include("abstract_types/AbstractEmpirical.jl")

    include("abstract_types/AbstractUncertainScalarKDE.jl")

    ##########################################
    # Distributions
    ##########################################
    include("distributions/assign_dist.jl")

    ##########################################
    # Fitted distribution type
    ##########################################
    include("distributions/fitted_distribution.jl")

    ##########################################
    # Composite uncertain scalar types
    ##########################################
    # Theoretical distributions with known parameters
    include("UncertainScalarsTheoretical.jl")

    # Theoretical distributions with fitted parameters
    include("UncertainScalarsTheoreticalFitted.jl")

    # Kernel density estimated distributions
    include("UncertainScalarsKDE.jl")

    # Populations with weighted probabilities
    include("UncertainScalarPopulation.jl")

    # Certain values (i.e. values without uncertainty)
    include("CertainValue.jl")

    ##########################################
    # Composite uncertain vector types
    ##########################################
    include("UncertainVectorsTheoretical.jl")
    include("UncertainVectorsTheoreticalFitted.jl")
    include("UncertainVectorsKDE.jl")

    # Define common constructor, so a similar syntax may be used to construct
    # all types of uncertain values.
    include("UncertainValue.jl")

    ##########################################
    # Comparisons
    ##########################################
    include("operations/comparisons.jl")

    # Conversions
    include("convert.jl")

    export KernelDensity, UnivariateKDE, default_bandwidth, kde
        ProbabilityWeights, pweights,
        FrequencyWeights, fweights,
        Weights, weights,
        AnalyticWeights, aweights
end #module

"""
	UncertainValues

A module defining uncertain value types.
"""
UncertainValues


text:
using PyCall

## Translate Symata to SymPy and vice versa.
##
## Symata -> SymPy is in general much faster than the reverse, 10x or more

## We don't need to define all of these.
const sympysyms1 = [:Symbol,:Number,:Integral, :Sum, :StrictLessThan, :StrictGreaterThan,
                    :LessThan, :GreaterThan, :Unequality,
                    :Integer, :Rational, :Float]

## These are buried in arbitrary class hierarchies
const sympysyms2 = [:True,:False]

const sympy  = PyCall.PyNULL()
const mpmath = PyCall.PyNULL()

for s in (sympysyms1..., sympysyms2...)
    @eval const $(Symbol("sympy_", s)) = PyCall.PyNULL()
end

import Base: isless

const PyObject = PyCall.PyObject

# Notes for evenutally extracting mpz and mpf from python mpmath objects
#  julia> ex[:_mpf_][2]
#  PyObject mpz(5400274185508743)
#  julia> ex
#  PyObject mpf('0.29977543700203352')

# Francesco Bonazzi contributed code for an early version of this file.
# FIXME: We may want to SymPy integers to BigInts in Symata when bigint_input is set.
# SymPy integers are already BigInt, but we are currently converting them to Int (I think)
# Eg.
# symata> BigIntInput(True)
# syamta> Head(Cos(0))
#     Int64

function import_sympy()
    copy!(sympy, PyCall.pyimport_conda("sympy", "sympy"))
    copy!(mpmath, PyCall.pyimport_conda("mpmath", "mpmath"))
end

# PYDEBUGLEVEL and SJDEBUGLEVEL defined in debug.jl
macro pydebug(level, a...)
    if level <= PYDEBUGLEVEL
        :((println("pydeb: ", $(a...));println()))
    else
        nothing
    end
end

macro sjdebug(level, a...)
    if level <= SJDEBUGLEVEL
        :((symprintln("sjdeb: ", $(a...));println()))
    else
        nothing
    end
end

######################################################
##   Translation Dicts
##   Following is somewhat out of data.

# SYMPY_TO_SYMATA_FUNCTIONS
# 1) keys are sympy function names. values are the Symata Head
# The sympy functions are found via sympy.key and are stored
# in the dict py_to_mx_dict, with values being the Symata heads.
# Some other special cases,
# eg Add Mul are also stored in py_to_mx_dict.
# The Symata heads are looked up only at the beginning
# of the main _pytosj routine.
# An Mxpr is created by looking up the head and mapping
# _pytosj over the args.
#
# 2) A reversed dict SYMATA_TO_SYMPY_FUNCTIONS is created from SYMPY_TO_SYMATA_FUNCTIONS.
# This dict is used for two things a) to automatically lookup sympy docstrings associated
# with Symata heads in doc.jl. b) to populate mx_to_py_dict. This dict stores callable references
# to sympy functions with keys being the Symata heads. A few other sympy functions are put into
# mx_to_py_dict "by hand" in populate_mx_to_py_dict().(does it need to be done this way ?)

const SYMPY_TO_SYMATA_FUNCTIONS = Dict{Symbol,Symbol}()
const SYMATA_TO_SYMPY_FUNCTIONS = Dict{Symbol,Symbol}()

# Dict containing symbols created via sympy.Symbol("...")
# This maps Sympy symbols to Symata (Julia) Symbols
const SYMPY_USER_SYMBOLS = Dict{Symbol,Any}()

# This maps in the other direction, Julia to Sympy symbols.
# But, we put more than user symbols in the following
const SYMPY_USER_SYMBOLS_REVERSE = Dict{PyCall.PyObject,Any}()  # type must be Any, we want more than symbols

function set_pytosj(py, sj)
    spy = Symbol(py)
    ssj = Symbol(sj)
    if haskey(SYMPY_TO_SYMATA_FUNCTIONS, spy)
        @warn("*** set_pytosj ", spy, " already has value ", SYMPY_TO_SYMATA_FUNCTIONS[spy], " can't set it to ", ssj)
        return
    end
    SYMPY_TO_SYMATA_FUNCTIONS[spy] = ssj
end

function set_sjtopy(sj,py)
    spy = Symbol(py)
    ssj = Symbol(sj)
    if haskey(SYMATA_TO_SYMPY_FUNCTIONS,ssj)
        @symwarn("!!! set_sjtopy ", sj, " already has value ", SYMATA_TO_SYMPY_FUNCTIONS[ssj], " can't set it to ", py)
        return
    end
    SYMATA_TO_SYMPY_FUNCTIONS[ssj] = spy
end

get_sjtopy(sj) = SYMATA_TO_SYMPY_FUNCTIONS[sj]

#######################################################

# described above
const mx_to_py_dict = Dict()   # Can we specify types for these Dicts ?

# This does not refer the sympy 'rewrite' capability.
# This refers to some kind of rewriting of functions or arguments
# during sympy to symata translation.
# This dict is populated below.
const py_to_mx_rewrite_function_dict = Dict{Any,Any}()

# End translation Dicts
######################################################

## Digits of precision for Sympy
const MPMATH_DPS = Int[0]

"""
    get_mpmath_dps()

Get the number of digits of precision used by mpmath.
"""
get_mpmath_dps() = mpmath[:mp][:dps]

"""
    set_mpmath_dps(n)

Set the number of digits of precision used by mpmath.
"""
function set_mpmath_dps(n)
    push!(MPMATH_DPS, mpmath[:mp][:dps])
    mpmath[:mp][:dps] = n
end

"""
    restore_mpmath_dps()

Pop the most recent number of digits of precision used by mpmath
from a stack, and set the number of DPS.
"""
restore_mpmath_dps() = (mpmath[:mp][:dps] = pop!(MPMATH_DPS))

# This translates key value pairs that we use to describe the mapping from sj to py.
function get_sympy_math(x)
    if length(x) == 1
        jf = x[1]
        st = string(jf)
        sjf = uppercasefirst(st)
    elseif length(x) == 2
        (jf,sjf) = x
    else
        error("get_sympy_math: wrong number of args in $x")
    end
    return jf,sjf
end

# These are used for rewriting in both directions and for calling functions via sympy.symbol
# Functions whose rules are written by hand, ie. that are not handled by
# the lists single_arg_float_complex, etc. must be entered here.
function make_sympy_to_symata()
    symbolic_misc = [ (:Order, :Order), (:LaplaceTransform, :laplace_transform),
                      ( :InverseLaplaceTransform, :inverse_laplace_transform ),
                      (:InverseFourierTransform, :inverse_fourier_transform ),
                      (:FourierTransform, :fourier_transform),
                      (:Log, :log), ( :Sqrt, :sqrt),
                      (:ProductLog, :LambertW),
                      (:Exp, :exp), (:Abs, :Abs), (:MeijerG, :meijerg), (:PolarLift, :polar_lift),
                      (:ExpPolar, :exp_polar), (:LowerGamma, :lowergamma), (:Sign,:sign),
                      (:PeriodicArgument, :periodic_argument), (:Max, :Max), (:Min, :Min),
                      (:DirichletEta, :dirichlet_eta),(:PolyLog, :polylog),
                      (:Conjugate, :conjugate), (:Factorial, :factorial)
                      ]

##  Cannot put Equality here. It is not a symbol. It is not a function. It is a something else.
## (:Equal,:Equality)

    for funclist in (single_arg_float_complex, single_arg_float_int_complex, single_arg_float,
                     single_arg_float_int, single_arg_int, two_arg_int,
                     two_arg_float_and_float_or_complex, two_arg_float,
                     one_or_two_args1
                     )
        for x in funclist
            if length(x) != 3 continue end
            (julia_func, symata_func, sympy_func) = x
            set_pytosj(sympy_func, symata_func)
            set_sjtopy(symata_func, sympy_func)
        end
    end

    for funclist in (symbolic_misc, no_julia_function, no_julia_function_one_arg,
                     no_julia_function_one_or_two_int,
                     no_julia_function_two_args, no_julia_function_two_or_three_args,
                     no_julia_function_three_args, no_julia_function_four_args)
        for x in funclist
            symata_func, sympy_func = get_sympy_math(x)
            set_pytosj(sympy_func, symata_func)
            set_sjtopy(symata_func, sympy_func)
        end
    end
    set_pytosj(:InverseLaplaceTransform, :InverseLaplaceTransform) # FIXME. Do this elsewhere ?
end

function register_sjfunc_pyfunc(sj::SymString, py::SymString)
    set_pytosj(Symbol(py), Symbol(sj))
    set_sjtopy(Symbol(sj), Symbol(py))
end

function register_only_pyfunc_to_sjfunc(sj::SymString, py::SymString)
    set_pytosj(py,sj)
end

## These two functions are only used in doc.jl to look up the
# docstring corresponding to the symp function called by an Symata head.
function have_pyfunc_symbol(sjsym)
    haskey(SYMATA_TO_SYMPY_FUNCTIONS, sjsym)
end
function lookup_pyfunc_symbol(sjsym)
    get_sjtopy(sjsym)
end

### Convert SymPy to Mxpr
#  Note: To convert, say a hypergeometric function to a float, do this:  f1[:evalf]()
function populate_py_to_mx_dict()
    # If this is Dict{Any,Any}, then nothing is translated until the
    # catchall branch at the end of sympyt2mxpr. Why ?
    eval(Meta.parse("const py_to_mx_dict = Dict{PyCall.PyObject,Symbol}()"))
    for onepair in (
                    (sympy[:Add], :Plus),
                    (sympy[:Mul], :Times),
                    (sympy[:Pow] ,:Power),
                    (sympy[:Derivative], :D),  ## We need to implement or handle Derivative
                    (sympy[:containers][:Tuple], :List),
                    (sympy[:oo], :Infinity),
                    (sympy[:zoo],:ComplexInfinity),
                    (sympy[:Eq], :Equal),
                    (sympy[:functions][:elementary][:piecewise][:ExprCondPair], :ConditionalExpression))
        py_to_mx_dict[onepair[1]] = onepair[2]
        if haskey(sympy[:functions][:special][:hyper], :TupleArg)  # This is missing in older versions of SymPy. (But so are many other symbols)
            py_to_mx_dict[sympy[:functions][:special][:hyper][:TupleArg]] = :List
        end
    end
end

# These functions are also contained in sympy.C (but sympy.C has been deprecated)
function mk_py_to_mx_funcs()
    for (pysym,sjsym) in SYMPY_TO_SYMATA_FUNCTIONS
        pystr = string(pysym)
        sjstr = string(sjsym)
        if haskey(sympy[:functions], pystr)
            py_to_mx_dict[sympy[:functions][pystr]] =  Symbol(sjstr)
        end
    end
end

function rewrite_function_sympy_to_julia(expr, pytype)
    func = py_to_mx_rewrite_function_dict[pytype]
    return func(expr)
end

have_function_sympy_to_symata_translation(pytype::PyCall.PyObject) = haskey(py_to_mx_dict, pytype)
get_function_sympy_to_symata_translation(pytype::PyCall.PyObject) = py_to_mx_dict[pytype]
have_rewrite_function_sympy_to_julia(pytype::PyCall.PyObject) = haskey(py_to_mx_rewrite_function_dict, pytype)

## We put more than user symbols here
function populate_user_symbol_dict()
    for onepair in (
                    (sympy[:numbers]["ComplexInfinity"](), :ComplexInfinity),
                    (sympy[:numbers]["Pi"](), :Pi),
                    (sympy[:numbers]["EulerGamma"](), :EulerGamma),
                    (sympy[:numbers]["Exp1"](),  :E),
                    (sympy[:numbers]["ImaginaryUnit"](), complex(0,1)),
                    (sympy[:numbers]["NegativeInfinity"](), MinusInfinity))
        SYMPY_USER_SYMBOLS_REVERSE[onepair[1]] = onepair[2]
    end
end

sympy_to_mxpr_symbol(s::Symbol) = get(py_to_mx_symbol_dict, s, s)
sympy_to_mxpr_symbol(s::String) = sympy_to_mxpr_symbol(Symbol(s))

function pytosj(x)
    getkerneloptions(:return_sympy) && return x  # The user can disable conversion. Eg. for debugging.
    maybetrace_pytosj(x)
end

mutable struct SympyTrace
    trace::Bool
end
const SYMPYTRACE = SympyTrace(false)
const SYMPYTRACENUM = Int[0]
get_pytosj_count() = SYMPYTRACENUM[1]
increment_pytosj_count() = SYMPYTRACENUM[1] += 1
decrement_pytosj_count() = SYMPYTRACENUM[1] -= 1
is_pytosj_trace() = SYMPYTRACE.trace

## This probably does not work, because we call the inner function _pytosj when doing recursion
##
function maybetrace_pytosj(expr)
    increment_pytosj_count()
    if is_pytosj_trace()
        ind = " " ^ (get_pytosj_count() - 1)
        println(ind,">>", get_pytosj_count(), " " , expr)
    end
    res = _pytosj(expr)
    if is_pytosj_trace()
        ind = " " ^ (get_pytosj_count() - 1)
        println(ind,"<<", get_pytosj_count(), " " , res)
    end
    decrement_pytosj_count()
    mergeargs(res)  # this may not be enough. only looks at one level, I think
    res
end

_pytosj(x) = x

function pytosj_map(head,args)
    nargs = newargs()
    for a in args
        push!(nargs,pytosj(a))
#        push!(nargs,_pytosj(a))  ## Using underscore here breaks some things. But pytosj_map is only called when _pytosj should be safe.
    end
    mxpra(head,nargs)
end

function pytosj_map2(head,args)
    nargs = newargs()
    for a in args
        push!(nargs,_pytosj(a))  ## Using underscore here
    end
    mxpra(head,nargs)
end

function pytosj_Function(pyexpr,pytype)
    head = Symbol(pytype[:__name__])
    targs = pyexpr[:args]
    if  head == :_context  return Qsym(pytosj(targs[1]),pytosj(targs[2])) end
    if targs[1] == dummy_arg  # sympy does not allow functions without args, so we pass a dummy arg.
        return mxprcf(head, [])
    else
        return pytosj_map2(head, targs)
    end
end

## TODO: we probably want to get rid of all of these.
## We are currently not translating !=, etc. correctly. At the moment it does not trigger any bugs.
##  Do not use Comparison.
##    Sympy        Symata
##    ==           SameQ     '==' is probably just python '=='. I can't convert this to sympy
##   Eq,Equality   Equal     Eq and Equality are aliases. In sympy Eq neither as symbol, nor a function. It is its own class
##                           This makes Eq more like Add and not Cos.
function populate_rewrite_dict()
    py_to_mx_rewrite_function_dict[sympy_True] = x -> true
    py_to_mx_rewrite_function_dict[sympy_False] = x -> false
    py_to_mx_rewrite_function_dict[sympy_Sum] = pyexpr -> deepsetfixed(mxpr(:Sum, mapmargs(pytosj, pyexpr[:args])...))
    py_to_mx_rewrite_function_dict[sympy_Integral] = pyexpr -> deepsetfixed(mxpr(:Integrate, mapmargs(pytosj, pyexpr[:args])...))

    for (fname,pyftype, sjsymbolstr) in (("greater_than_equal", :sympy_GreaterThan, ">="),
                                         ("less_than_equal", :sympy_LessThan, "<="),
                                         ("less_than", :sympy_StrictLessThan, "<"),
                                         ("greater_than", :sympy_StrictGreaterThan, ">"),
#                                         ("equality", :sympy_Equality, "=="),
                                         ("unequality", :sympy_Unequality, "!="))

        sfname = Symbol("pytosj_" * fname)
        sjsymbol = QuoteNode(Symbol(sjsymbolstr))
        @eval begin
            function ($sfname)(pyexpr)
                args = pyexpr[:args]
                return mxpr(:Comparison, _pytosj(args[1]), $sjsymbol, _pytosj(args[2]))
            end
            py_to_mx_rewrite_function_dict[$(eval(pyftype))] = $sfname
        end
    end
end

####
####   Main _pytosj method
####
## expr[:__class__] is a bit faster and does fewer allocations
function _pytosj(expr::T) where T <: PyCall.PyObject
    res = get(SYMPY_USER_SYMBOLS_REVERSE, expr,false) ## This is much faster, but may not be significant
    res !== false && return res
    pytype = expr[:__class__]
    if pytype == sympy_Symbol
        @pydebug(3, "pytype Symbol trans. ", expr)
        return Symbol(expr[:name])
    end
    @pydebug(3, "Entering with ", expr)
#    if pytype == sympy_Integer  ## This does not work. 1, 1/2, etc, are in their own classes.
    ## Many things apparently cannot be checked via the 'pytype' or class.
    ## There are paraphyletic groups whose membership is tested with functions (well, class member functions)
    if expr[:is_Integer]
        #  returns Int or BigInt. It just returns a member of a python object. It would be nice to get this faster.
        # `_to_mpath` in evalf.py checks again expr.is_Integer.
#        num = expr[:_to_mpmath](-1)  #  precision '-1' is ignored for integers.
        num = expr[:p]  ## This is 5 or so times faster and does the same thing.
        return num
    end
#    if pytype == sympy_Rational
    if expr[:is_Rational]
        p = expr[:p]
        q = expr[:q]
        return Rational(expr[:p],expr[:q])  # These are Int64's. Don't know when or if they are BigInts
    end
#    if pytype == sympy_Float
    if expr[:is_Float]
        return convert(AbstractFloat, expr) # Need to check for BigFloat
    end
    if pyisinstance(expr, sympy[:Number]) && (! expr[:is_finite])
            # we should check name. there are several infinities
        return Infinity  # need to see if this is maybe -Infinity
    end
    if have_function_sympy_to_symata_translation(pytype)
        @pydebug(3, "function lookup trans. ", expr)
        return pytosj_map(get_function_sympy_to_symata_translation(pytype), expr[:args])
    end
    if expr[:is_Function] ## User defined functions whose pytpe is the name of the function
        @pydebug(3, "is_Function trans. ", expr)
        return pytosj_Function(expr,pytype)
    end   # perhaps a user defined function
    if have_rewrite_function_sympy_to_julia(pytype)
        @pydebug(3, "rewrite trans. ", expr)
        return rewrite_function_sympy_to_julia(expr,pytype)
    end
    head = Symbol(expr[:func][:__name__]) # this is the default if not handled above
    @pydebug(3, "default trans. ", expr, " new head ", head)
    @pydebug(3, "typeof head ", typeof(head))
    @pydebug(3, "args  ", expr[:args])
    mxout =  pytosj_map(head,  expr[:args])
    @pydebug(3, "mxout ", mxout)
    mxout
end

# By default, Dict goes to Dict
function _pytosj(expr::Dict)
    ndict = Dict()
    for (k,v) in expr
        ndict[pytosj(k)] = pytosj(v)
    end
    return ndict
end

function _pytosj(expr::Array{T,1}) where T
    return pytosj_map(:List, expr)
end

function _pytosj(expr::Tuple)
    return pytosj_map(:List, expr)
end

#### Convert Mxpr to SymPy

## We use instances here and classes to go the other direction.
## This is because the name cannot be obtained from an instance.
## TODO: I think the above is not true. We can merge this with another dict
function populate_mx_to_py_dict()
    for onepair in (
         (:Plus,sympy[:Add]),
         (:Times, sympy[:Mul]),
         (:Power, sympy[:Pow]),
         (:E, sympy[:E]),
         (:EulerGamma, sympy[:EulerGamma]),
         (:I, sympy[:I]),
         (:Pi,  sympy[:pi]),
         (:Log, sympy[:log]),
         (:Infinity, sympy[:oo]),
         (:ComplexInfinity, sympy[:zoo]))
        mx_to_py_dict[onepair[1]] = onepair[2]
    end
end

function mk_mx_to_py_funcs()
    for (sjsym,pysym) in SYMATA_TO_SYMPY_FUNCTIONS
        pystr = string(pysym)
        sjstr = string(sjsym)
        local obj
        if length(pystr) > 6 && pystr[1:6] == "sympy_"  # These are wrapper function around sympy functions
            try
                obj = eval(pysym)
            catch
                continue
            end
        else
            try
                obj = eval(Meta.parse("sympy[:" * pystr * "]"))   # These call the sympy functions directly
            catch
                continue
            end
        end
        mx_to_py_dict[Symbol(sjstr)] = obj
    end
end

#######################################################
##     sjtopy
#######################################################

# NB wrapper to call _sjtopy
function sjtopy(args...)
    if length(args) == 1
        res = _sjtopy(args[1])
        return res
    end
    res = _sjtopy(args)
    (res...,)
end

function _sjtopy(z::Complex)
    @sjdebug(3,"complex ", z)
    if real(z) == 0
        res = mxpr(:Times, :I, imag(z))
    else
        res = mxpr(:Plus, real(z), mxpr(:Times, :I, imag(z)))
    end
    return _sjtopy(res)
end

## FIXME: this can probably be simplified, but we we would need to test it.
function _sjtopy(mx::Mxpr{:List})
    @sjdebug(3,"List ", mx)
    a = Array{Any}(undef, 0)
    rs = map(_sjtopy, mx.args)
    for el in rs
        push!(a,el)
    end
    return a
#    return [map(_sjtopy, mx.args)...]  This fails for some versions of Julia
end

# For now all infinities are mapped to one of two infinities
function _sjtopy(mx::Mxpr{:DirectedInfinity})
    @sjdebug(3,"Infinity ", mx)
    if mx == ComplexInfinity
        sympy[:zoo]
    elseif mx == MinusInfinity
        SymPyMinusInfinity
    elseif mx[1] == I
        sympy[:Mul](sympy[:I],sympy[:oo])
    else
        sympy[:oo]
    end
end

### Equal

function _sjtopy(mx::Mxpr{:Equal})
    if length(mx) == 2
        sympy[:Eq](_sjtopy(mx[1]),_sjtopy(mx[2]))
    else
        @symwarn("Unimplemented translation to Sympy", mx)
    end
end

### Part  (defined in parts.jl)

function get_part_one_ind(pyobj::PyObject, ind)
    if pyhead(pyobj) == :Tuple
        return pyobj[ind]
    end
    symerror(string("Can't find part number $ind of ", pyobj))
end

### Args, This function is defined in lists.jl

@doap Args(pyobj::PyCall.PyObject) = pyobj[:args]

### PyHead

function pyhead(pyobj::PyCall.PyObject)
    return Symbol(pyobj[:func][:__name__])
end

@mkapprule PyHead :nargs => 1

@sjdoc PyHead """
    PyHead(pyobj)

return the sympy "Head" of `pyobj` as a Symata (Julia) symbol.
"""

@doap PyHead(pyobj::PyCall.PyObject) = pyhead(pyobj)

### PyClass

@sjdoc PyClass """
    PyClass(pyobj)

return the python class of the python object `pyobj`.
"""

@mkapprule PyClass :nargs => 1
@doap PyClass(pyobj::PyCall.PyObject) = pyobj[:__class__]

### Rewrite

const RewriteDict = Dict( :Gamma => :gamma,
                          :HypergeometricPFQ => :hyper,
                          :Binomial => :binomial,
                          :Sin => :sin)

@mkapprule Rewrite :nargs => 2

@sjdoc Rewrite """
    Rewrite(expr,form)

rewrite `expr` in term of `form`. This is
implemented for some functions that call SymPy. The second argument is
translated, but you may give the SymPy symbol, as well

## Examples

```
Rewrite(CatalanNumber(n), Gamma)
Rewrite(CatalanNumber(1/2), Gamma)
Rewrite(CatalanNumber(n), HypergeometricPFQ)
```
"""

@doap function Rewrite(expr, form::Symbol)
    arg1 = _sjtopy(mx[1])
    transform = get(RewriteDict, form, form)
    pyres = arg1[:rewrite](transform)
    pytosj(pyres)
end

### HypergeometricPFQ

do_HypergeometricPFQ(mx::Mxpr{:HypergeometricPFQ}, p::Mxpr{:List}, q::Mxpr{:List}, z::W) where {W<:AbstractFloat} =
    eval_hypergeometric(mx,p,q,z)

do_HypergeometricPFQ(mx::Mxpr{:HypergeometricPFQ}, p::Mxpr{:List}, q::Mxpr{:List}, z::Complex{W}) where {W<:AbstractFloat} =
    eval_hypergeometric(mx,p,q,z)

function eval_hypergeometric(mx, p, q, z)
    result = sjtopy(mx)
    fresult =
        try
            result[:evalf]()
        catch
            result
        end
    pytosj(fresult)
end

#### MeijerG

@mkapprule MeijerG

@sjdoc MeijerG """
    MeijerG

Meijer G function.
"""

function do_MeijerG(mx::Mxpr{:MeijerG}, p::Mxpr{:List}, q::Mxpr{:List}, z)
    mxc = copy(mx)
    mxc[1] = (p[1],p[2])
    mxc[2] = (q[1],q[2])
    try
        zpy = sjtopy(z)
        ppy = (_sjtopy(p[1]), _sjtopy(p[2]))
        qpy = (_sjtopy(q[1]), _sjtopy(q[2]))
        pyres = sympy[:meijerg](ppy,qpy,zpy)
        sjres = pytosj(pyres)
    catch
        mx
    end
end

function _sjtopy(mx::Mxpr{:MeijerG})
    @sjdebug(3,"MeijerG ", mx)
    pyhead = mx_to_py_dict[mhead(mx)]
    p = mx[1]
    q = mx[2]
    z = mx[3]
    pyhead((_sjtopy(p[1]), _sjtopy(p[2])), (_sjtopy(q[1]), _sjtopy(q[2])), _sjtopy(z))
end

do_MeijerG(mx::Mxpr{:MeijerG}, p::Mxpr{:List}, q::Mxpr{:List}, z::W) where {W<:AbstractFloat} =
    eval_meijerg(mx,p,q,z)

do_MeijerG(mx::Mxpr{:MeijerG}, p::Mxpr{:List}, q::Mxpr{:List}, z::Complex{W}) where {W<:AbstractFloat} =
    eval_meijerg(mx,p,q,z)

function eval_meijerg(mx, p, q, z)
    mxc = copy(mx)
    mxc[1] = (p[1],p[2])
    mxc[2] = (q[1],q[2])
    result = sjtopy(mxc)
    fresult =
        try
            result[:evalf]()
        catch
            result
        end
    pytosj(fresult)
end

#### Tuple

function _sjtopy(t::Tuple)
    @sjdebug(3,"Tuple ", mx)
    (map(_sjtopy,t)...,)
end

## Sympy functions are called here.
function _sjtopy(mx::Mxpr)
    @sjdebug(3,"Mxpr ", mx)
    if mhead(mx) in keys(mx_to_py_dict)
        @sjdebug(1,"In mx_to_py_dict ", mx, " pyhead ", mx_to_py_dict[mhead(mx)] )
        return mx_to_py_dict[mhead(mx)](map(_sjtopy, mx.args)...) # calling a function in our dictionary
    end
    @sjdebug(1,"Make function ", mx)
    pyfunc = sympy[:Function](string(mhead(mx)))  # Don't recognize the head, so make it a user function
    mxargs = margs(mx)
    if length(mxargs) == 0
        return pyfunc(dummy_arg)  # sympy functions must have at least one argument
    else
        return pyfunc(map(_sjtopy, mxargs)...)
    end
end

#### Symbol

# mx is a julia Symbol.
# We create a sympy symbol, unless it has already been created
function _sjtopy(mx::Symbol)
    @sjdebug(3,"Symbol ", mx)
    if haskey(mx_to_py_dict,mx)
        @sjdebug(1,"In mx_to_py_dixt ", mx)
        return mx_to_py_dict[mx]
    end
    if haskey(SYMPY_USER_SYMBOLS,mx)
        sym = SYMPY_USER_SYMBOLS[mx]
        return sym
    else
        sym = sympy[:Symbol](mx)
        SYMPY_USER_SYMBOLS[mx] = sym
        SYMPY_USER_SYMBOLS_REVERSE[sym] = mx
        return sym
    end
end

#### Assume

@mkapprule Assume

@sjdoc Assume """
    Assume(sym,prop1,prop2,...)

sets properties `prop1, prop2,...` for symbols `sym`.

These properties are used in simplifying relations, `Refine`, integral transforms, etc.

Properties are: commutative, complex, imaginary, real, integer, odd,
even, prime, composite, zero, nonzero, rational, algebraic,
transcendental, irrational, finite, infinite, negative, nonnegative,
positive, nonpositive, hermitian, antihermitian.

These property symbols may contain uppercase characters. For example
`Positive` and `positive` are the same. Use of properties with
lowercase initial is deprecated.
"""

# TODO. We implemented this to use, e.g. the inequality solvers
# Can't find a better way to create the symbol
# prop is not evaluated, rather the symbol 'prop' is set to true
# This uses the 'old' SymPy assumptions system
@doap function Assume(s::Symbol, inprops...)
    ss = string(s)
    props = [ Symbol(lowercase(string(x))) for x in inprops]
    propstrs = map( (x) -> " $x=true ", props)
    propstr = join(propstrs, ", ")
    evalstr = "sympy[:Symbol](\"$ss\", $propstr)"
    sym = eval(Meta.parse(evalstr))
    SYMPY_USER_SYMBOLS[s] = sym
    s
end

#### Max

# TODO Max(x,y) == Max(y,x)
# TODO Do numerical arguments in Julia
# TODO Flatten lists
@mkapprule Max :nodefault => true
@doap function Max(args...)
    args = margs(flatten_recursive(mxpr(:List,args...)))
    return pytosj(sympy[:Max](mapmargs(sjtopy, args)...))
end

@doap Max() = MinusInfinity

#### Min

@mkapprule Min :nodefault => true
@doap function Min(args...)
    args = margs(flatten_recursive(mxpr(:List,args...)))
    return pytosj(sympy[:Min](mapmargs(sjtopy, args)...))
end

@doap Min() = Infinity

_sjtopy(mx::Rational{T}) where {T<:Integer} = sympy[:Rational](numerator(mx),denominator(mx))

_sjtopy(mx::Number) = mx

function _sjtopy(s::Qsym)
    f = sympy[:Function]("_context")
    f(_sjtopy(s.context),_sjtopy(s.name))
end

# For our LaplaceTransform code, (etc.)
_sjtopy(a::Array{T,1}) where {T} =  map(_sjtopy, a)
_sjtopy(expr::PyCall.PyObject) = expr
_sjtopy(x::AbstractString) =  x

# Default conversion. Don't error, but only warn. Then return x so that we
# can capture and inspect it.
function _sjtopy(x)
    @symwarn("sjtopy: Unable to convert $x from Symata to SymPy")
    return x
end

"""
    init_sympy()

Initialize SymPy. SymPy must be loaded at runtime. `init_sympy()` is called
in the Symata module ` __init__` function.
"""
function init_sympy()
    import_sympy()
    eval(Meta.parse("const dummy_arg = sympy[:Symbol](\"DUMMY\")"))
    eval(Meta.parse("const SymPyMinusInfinity = sympy[:Mul](-1 , sympy[:oo])"))
    make_sympy_to_symata()
    populate_py_to_mx_dict()
    mk_py_to_mx_funcs()
    populate_user_symbol_dict()
    populate_mx_to_py_dict()
    mk_mx_to_py_funcs()
    for s in sympysyms1
        @eval copy!($(Symbol("sympy_", s)) , sympy[$(QuoteNode(s))])
    end
    copy!(sympy_True,sympy[:boolalg][:BooleanTrue])
    copy!(sympy_False,sympy[:boolalg][:BooleanFalse])
    populate_rewrite_dict()
    init_isympy()
end

const symatapydir = joinpath(dirname(@__FILE__), "..", "pysrc")

const symatapy = PyCall.PyNULL()

function init_isympy()
    try
        push!(pyimport("sys")["path"], symatapydir)
        copy!(symatapy, pyimport("isympy"))
    catch
        nothing
    end
end

"""
    isympy()

enters the ipython shell with sympy imported.
"""
function isympy()
    try
        symatapy[:main]()
    catch
        nothing
    end
end

## These are not used anywhere
# typename(x::PyCall.PyObject) = pytypeof(x)[:__name__]
# name(x::PyCall.PyObject) = x[:__name__]

#  Try the sympy function 'pycall'. If there is an error, give warning
#  'errstr' and return (from surrounding function body)
#  'return_val_err' Store the error message in the kernel state.
#  On success, return the result of the function call.
macro try_sympyfunc(pycall, errstr, return_val_err)
    qpycall = QuoteNode(pycall)
    return :(
             begin
             (sflag, _pyres) =
                 try
                   (true, $(esc(pycall)))
                 catch pyerr
                  (false,pyerr)
                 end
                 if sflag == false
                 @warn($errstr)
                   setkerneloptions(:sympy_error, _pyres)
                   return $(esc(return_val_err))
                 end
                 _pyres
              end
           )
end

#### PyDoc

@mkapprule PyDoc :nargs => 1

function do_PyDoc(mx::Mxpr{:PyDoc},sym)
    try eval(Meta.parse("println(sympy[:$(string(sym))][:__doc__])"))
    catch
        Null
    end
end

@sjdoc PyDoc """
    PyDoc(sym)

prints the documentation for the symbol `sym`, if available.
"""

#### Arithmetic

## Define methods for Sympy objects.

## These break the ability to include sympy objects in Symata expressions with
## heads such as :Times, via `*`. It's not clear which, if either is more useful

## PyCall.jl has now implemented these.
## I comment out the definitions here, so they don't conflict. GJL 11MAY2018
# for sympair in ( (:*, :Mul), (:+, :Add), (:^, :Pow) )
#     jop = sympair[1]
#     qj = QuoteNode(sympair[1])
#     qsp = QuoteNode(sympair[2])
#     @eval Base.$(jop)(a::PyCall.PyObject, b::PyCall.PyObject) = sympy[$qsp](a,b)
#     @eval Base.$(jop)(a::PyCall.PyObject, b) = sympy[$qsp](a,b)
#     @eval Base.$(jop)(a, b::PyCall.PyObject) = sympy[$qsp](a,b)
#     @eval Base.$(jop)(a::PyCall.PyObject, b::Number) = sympy[$qsp](a,b)
#     @eval Base.$(jop)(a::PyCall.PyObject, b::Integer) = sympy[$qsp](a,b)
# end

# Base.:-(a::PyCall.PyObject) = -1 * a
# Base.:-(a::PyCall.PyObject, b::PyCall.PyObject) = a + (-b)
# Base.:-(a, b::PyCall.PyObject) = a + (-b)
# Base.:-(a::PyCall.PyObject, b) = a + (-b)


text:
#
# This file is part of module
#   ModiaMedia (ModiaMedia/src/ModiaMedia.jl)
#

"""
    data = SimpleMediumData(;cp_const=nothing, cv_const=nothing, d_const=nothing, eta_const=nothing,
                             lambda_const=nothing, a_const=nothing, T0=nothing, MM_const=nothing)

Generate a `SimpleMediumData` object containing the data
for a SimpleMedium medium.
"""
mutable struct SimpleMediumData
    cp_const::Float64       # "Constant specific heat capacity at constant pressure";
    cv_const::Float64       # "Constant specific heat capacity at constant volume";
    d_const::Float64        # "Constant density";
    eta_const::Float64      # "Constant dynamic viscosity";
    lambda_const::Float64   # "Constant thermal conductivity";
    a_const::Float64        # "Constant velocity of sound";
    T_min::Float64          # Minimum temperature valid for medium model
    T_max::Float64          # Maximum temperature valid for medium model
    T0::Float64             # "Zero enthalpy temperature";
    MM_const                # "Molar mass";

    SimpleMediumData(;cp_const=nothing, cv_const=nothing, d_const=nothing, eta_const=nothing,
                      lambda_const=nothing, a_const=nothing, T_min=nothing, T_max=nothing,
                      T0=nothing, MM_const=nothing) =
         new(cp_const, cv_const, d_const, eta_const, lambda_const, a_const, T_min, T_max, T0, MM_const)
end

"""
    medium = SimpleMedium(; mediumName     = nothing,
                            reference_p    = 101325.0,
                            reference_T    = 298.15,
                            p_default      = 101325.0,
                            T_default      = 293.15,
                            fluidConstants = nothing,
                            data           = nothing)

Generate a `SimpleMedium <: PureSubstance` medium object.
"""
mutable struct SimpleMedium <: PureSubstance
    infos::FluidInfos
    fluidConstants::SVector{1,BasicFluidConstants}
    fluidLimits::FluidLimits
    data::SimpleMediumData
end

"""
    state = SimpleMediumState(Medium, p, T)

Generate a `SimpleMediumState <: ThermodynamicState` object containing
pressure `p` [Pa] and temperature `T` [K] as thermodynamic states.
"""
mutable struct SimpleMediumState <: ThermodynamicState
    Medium::SimpleMedium
    p::Float64
    T::Float64
end

function SimpleMedium(; mediumName=nothing,
                        reference_p=101325.0,
                        reference_T=298.15,
                        p_default=101325.0,
                        T_default=293.15,
                        fluidConstants=nothing,
                        data=nothing)

    infos = FluidInfos(mediumName           = mediumName,
                       substanceNames       = [mediumName],
                       extraPropertiesNames = fill("",0),
                       ThermoStates         = IndependentVariables_pT,
                       singleState          = true,
                       reducedX             = true,
                       fixedX               = false,
                       reference_p          = reference_p,
                       reference_T          = reference_T,
                       reference_X          = fill(1.0,1),
                       p_default            = p_default,
                       T_default            = T_default)

    fluidLimits = FluidLimits(TMIN = data.T_min,
                              TMAX = data.T_max)

    Medium = SimpleMedium(infos, fill(fluidConstants,1), fluidLimits, data)

    infos.h_default  = specificEnthalpy(data, SimpleMediumState(Medium,p_default,T_default))
    fluidLimits.HMIN = specificEnthalpy(data, SimpleMediumState(Medium,p_default,data.T_min))
    fluidLimits.HMAX = specificEnthalpy(data, SimpleMediumState(Medium,p_default,data.T_max))

    return Medium
end

setState_pTX(m::SimpleMedium,p,T,X) = SimpleMediumState(m,p,T)
setState_phX(m::SimpleMedium,p,h,X) = SimpleMediumState(m,p,m.data.T0+h/m.data.cp_const)
setState_psX(m::SimpleMedium,p,s,X) = SimpleMediumState(m,p,exp(s/m.data.cp_const + log(m.infos.reference_T)))
setState_dTX(m::SimpleMedium,d,T,X) = error("From setState_dTX: Pressure cannot be computed from temperature and density for the incompressible fluid $(m.infos.mediumName)!")
isenthalpicState(m::SimpleMedium, state::SimpleMediumState, dp::Float64) = SimpleMediumState(m,state.p+dp, state.T)

setState_pTX!(state::SimpleMediumState,p,T,X) = begin state.p=p; state.T=T; nothing end
setState_phX!(state::SimpleMediumState,p,h,X) = begin state.p=p; state.T=state.Medium.data.T0+h/state.Medium.data.cp_const; nothing end
setState_psX!(state::SimpleMediumState,p,s,X) = begin state.p=p; state.T=exp(s/state.Medium.data.cp_const + log(state.Medium.infos.reference_T)); nothing end
setState_dTX!(state::SimpleMediumState,d,T,X) = error("From setState_dTX!: Pressure cannot be computed from temperature and density for the incompressible fluid $(state.Mediumm.infos.mediumName)!")
isenthalpicState!(state_b::SimpleMediumState, state_a::SimpleMediumState, dp::Float64) = begin state_b.p = state_a.p+dp; state_b.T = state_a.T; nothing end

specificEnthalpy(data::SimpleMediumData, state::SimpleMediumState)::Float64 = data.cp_const*(state.T - data.T0)
pressure(               m::SimpleMedium, state::SimpleMediumState)::Float64 = state.p
temperature(            m::SimpleMedium, state::SimpleMediumState)::Float64 = state.T
density(                m::SimpleMedium, state::SimpleMediumState)::Float64 = m.data.d_const
specificEnthalpy(       m::SimpleMedium, state::SimpleMediumState)::Float64 = specificEnthalpy(m.data,state)
specificInternalEnergy( m::SimpleMedium, state::SimpleMediumState)::Float64 = m.data.cv_const*(state.T - m.data.T0)
specificHeatCapacityCp( m::SimpleMedium, state::SimpleMediumState)::Float64 = m.data.cp_const
dynamicViscosity(       m::SimpleMedium, state::SimpleMediumState)::Float64 = m.data.eta_const

specificEnthalpy_T(m::SimpleMedium, T)::Float64 = m.data.cp_const*(T - m.data.T0)

### Functions for mass and energy balance

density(      m::SimpleMedium)::Float64 = m.data.d_const
density_der_1(m::SimpleMedium)::Float64 = 0.0

specificInternalEnergy_T(      m::SimpleMedium,T)::Float64 = m.data.cv_const*(T - m.data.T0)
specificInternalEnergy_T_der_1(m::SimpleMedium,T)::Float64 = 0.0
specificInternalEnergy_T_der_2(m::SimpleMedium,T)::Float64 = m.data.cv_const

function standardCharacteristics(m::SimpleMedium)::Dict{AbstractString,Any}
    p = m.infos.reference_p
    T = collect( range(m.fluidLimits.TMIN, stop=m.fluidLimits.TMAX, length=101) )
    h = zeros(length(T))
    u = zeros(length(T))
    state = setState_pT(m,p,T[1])

    for i in 1:length(T)
        setState_pT!(state,p,T[i])
        h[i] = specificEnthalpy(state)
        u[i] = specificInternalEnergy(state)
    end

    mediumDict = Dict{AbstractString,Any}()
    mediumDict["T"]   = uconvert.(u"°C", T*1u"K")
    mediumDict["h"]   = h*u"J/kg"
    mediumDict["u"]   = u*u"J/kg"
    mediumDict["d"]   = to_DensityDisplayUnit(density(state))*1u"g/cm^3"

    return mediumDict
end

function standardPlot(m::SimpleMedium; figure=1)
    mediumDict = standardCharacteristics(m)
    ModiaMath.plot(mediumDict, [("h", "u"), "d"], xAxis="T", heading=m.infos.mediumName, figure=figure)
end


text:
using Test
using WaveProp.Geometry

@testset "HyperRectangle tests" begin
    low_corner  = (0.0,0.0)
    high_corner = (1.0,2.0)
    mid         = (low_corner .+ high_corner) ./ 2 |> SVector
    rec = HyperRectangle(low_corner,high_corner)
    @test mid == center(rec)
    @test (mid ∈ rec)  == true
    @test !in(high_corner .+ (1,1),rec)
    @test diameter(rec) == sqrt(1^2 + 2^2)
    @test radius(rec) == sqrt(1^2 + 2^2)/2
    # bbox
    pts = NTuple{2,Float64}[]
    for x=-1:0.1:1
        for y=-1:0.1:1
            push!(pts,(x,y))
        end
    end
    @test bounding_box(pts) == HyperRectangle((-1.,-1),(1,1.))
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
using FEM_P2
using PyPlot

bc = (Dirichlet, Neumann)
L = 2.0
γ0 = 1.0
γL = -1.5
a = 1.0
c = 0.0
f(x) = 2.0*a
exact_u(x) = γ0 + (γL+2L)*x - x^2

mesh = Mesh(L .* [0.0, 0.25, 0.6, 1.0], bc)
A = assemble_matrix(mesh, deriv_times_deriv!, a)
F = assemble_vector(mesh, f)
G = Neumann_bc_vector(mesh, [γL])
F .+= G
ufix = [γ0]
ndof = mesh.ndof
F .-= A[:,ndof+1:end] * ufix

ufree = A[:,1:ndof] \ F
Und = [ ufree; ufix ]
x, u = evaluate_from_ndvals(mesh, Und)
xnd, und = evaluate_from_ndvals(mesh, Und, 2)

figure(1)
plot(x, u, "y-" , xnd, und, "o", x, exact_u.(x), "r:")
legend(("FEM", "Nodal values", "Exact solution"),
       loc="lower center")
grid(true)
xlabel(L"x")
title("Quadratic Finite Elements")

max_stars_repo_name: billmclean/ComputationalMathsNotes

max_issues_repo_name: billmclean/ComputationalMathsNotes

max_forks_repo_name: billmclean/ComputationalMathsNotes


text:
using Inform
using Base.Test

@test Inform._libinform_found
@test Inform._libinform != C_NULL
@test isa(Inform._libinform, Ptr{Void})

include("activeinfo.jl")


text:
function revers(side :: HorizonSide)
    side = HorizonSide((Int(side)+2)%4)
    return side
end

function movements!(r,sides)

    last_side = nothing;
    pl_side = nothing;

    while (isborder(r, sides[1]) ==  false) && (isborder(r, sides[2]) ==  false)
        for side in sides
            if isborder(r,side) == false
                move!(r,side)
                pl_side = last_side
                last_side = side
            end
        end
        putmarker!(r)
    end
    putmarker!(r)
    return revers(last_side), revers(pl_side)
end

function movements_mark(r, revers_sides)

    while ismarker(r) == true
        for side in revers_sides
            move!(r,side)
        end
    end

end

function krest(r)

    for sides in [ (Nord, Ost), (Sud, Ost), (Sud, West), (Nord, West) ]
        revers_sides = movements!(r,sides)
        print(revers_sides)
        movements_mark(r,revers_sides)
    end

    putmarker!(r)
    show(r)
end

krest(r)
show(r)

max_stars_repo_name: AntonDudDm/mirea-progs

max_issues_repo_name: AntonDudDm/mirea-progs

max_forks_repo_name: AntonDudDm/mirea-progs


text:
@testset "isconvex" begin
    m = JuMP.Model()

    JuMP.@variables m begin
       x
       y
       z
   end

   # AffExpr
   @test MultilinearOpt.isconvex(x + y)
   @test MultilinearOpt.isconvex(x - z - 3)

   # QuadExpr
   @test MultilinearOpt.isconvex(x^2)
   @test MultilinearOpt.isconvex(x^2 + 0 * z^2) # test positive semidefinite gramian
   @test MultilinearOpt.isconvex(x^2 + 3 * y - z + 5)
   @test !MultilinearOpt.isconvex(-x^2)
   @test !MultilinearOpt.isconvex(x + y - z ^2 + x^2 + y^2)
   G = rand(3, 3)
   G = G * G'
   vars = [x, y, z]
   expr = dot(vars, G * vars)
   G_back, vars_back = MultilinearOpt.gramian(expr)
   @test vars_back == vars
   @test G_back ≈ G atol=1e-10
   @test MultilinearOpt.isconvex(expr)

   # LinearConstraint
   @test MultilinearOpt.isconvex(JuMP.constraint_object(JuMP.@constraint(m, x + 3 * y == 0)))
   @test MultilinearOpt.isconvex(JuMP.constraint_object(JuMP.@constraint(m, 2 * x - y >= z)))
   @test MultilinearOpt.isconvex(JuMP.constraint_object(JuMP.@constraint(m, 2 * x - y <= z)))

   # QuadConstr
   @test !MultilinearOpt.isconvex(JuMP.constraint_object(JuMP.@constraint(m, x == y * z)))

   @test MultilinearOpt.isconvex(JuMP.constraint_object(JuMP.@constraint(m, x^2 + y^2 <= z)))

   @test !MultilinearOpt.isconvex(JuMP.constraint_object(JuMP.@constraint(m, x^2 + y^2 <= z^2)))

   @test MultilinearOpt.isconvex(JuMP.constraint_object(JuMP.@constraint(m, -x^2 - y^2 >= -z)))
end