text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Solids (Modia3D/Solids/_module.jl)
#
# Content:
#
# Utility functions for SolidFileMesh.jl
# function getObjInfos read file and returns: centroid, longestEdge and all objPoints

function getObjInfos(filename::AbstractString, scaleFactor::MVector{3,Float64})
  objPoints = []
  facesIndizes = []
  areTriangles = true
  centroid = ModiaMath.ZeroVector3D
  longestEdge = 0
  if filename[end-3:end] == ".obj"
    open(filename,"r") do file
    i = 0
    x_max::Float64 = 0.0
    x_min::Float64 = 0.0
    y_max::Float64 = 0.0
    y_min::Float64 = 0.0
    z_max::Float64 = 0.0
    z_min::Float64 = 0.0
    sum = ModiaMath.ZeroVector3D

    for line in eachline(file; keep=true)
      if isequal(line[1],'v') && isequal(line[2],' ')
        i += 1
        tmp = rsplit(line,' ')
        push!(objPoints,[parse(Float64,tmp[end-2])*scaleFactor[1],parse(Float64,tmp[end-1])*scaleFactor[2],parse(Float64,tmp[end])*scaleFactor[3]])
        if i == 1
          x_max = objPoints[i][1]
          x_min = objPoints[i][1]
          y_max = objPoints[i][2]
          y_min = objPoints[i][2]
          z_max = objPoints[i][3]
          z_min = objPoints[i][3]
        else
          (x_min, x_max) = check_MinMax(x_min, x_max, objPoints[i][1])
          (y_min, y_max) = check_MinMax(y_min, y_max, objPoints[i][2])
          (z_min, z_max) = check_MinMax(z_min, z_max, objPoints[i][3])
        end
        sum += objPoints[i]
      end
      if isequal(line[1],'f') && isequal(line[2],' ')
        if areTriangles
            tmp = rsplit(line,' ')
            if length(tmp) == 4 || ( length(tmp) == 5 && isequal(tmp[5],"\r\n") )
                push!(facesIndizes,[parse(Int64,rsplit(tmp[2],"/")[1]), parse(Int64,rsplit(tmp[3],"/")[1]), parse(Int64,rsplit(tmp[4],"/")[1])])
            else
                areTriangles = false
    end; end; end; end

    if length(objPoints) != 0
      centroid = sum / length(objPoints)
    else
      centroid = ModiaMath.ZeroVector3D
    end
    if !areTriangles
      empty!(facesIndizes)
    end
    longestEdge = max((x_max-x_min), (y_max-y_min), (z_max-z_min))
    end
    return (centroid, longestEdge, objPoints, facesIndizes)
  else
    println("Just .obj files are supported.")
  end
end

getObjInfos(filename::AbstractString, scaleFactor::AbstractVector) = getObjInfos(filename, MVector{3,Float64}(scaleFactor))

# utility function for getObjInfos
function check_MinMax(act_min::Float64, act_max::Float64, value)
  if act_min >= value
    act_min = value
  end
  if act_max <= value
    act_max = value
  end
  return (act_min, act_max)
end

"""
VHACD(filename::String, outputDirectory::String, outputBasename::String;
              showLog=false, resolution=100000, depth=20, concavity=0.0025, planeDownsampling=4, convexhullDownsampling=4,
              alpha=0.05, beta=0.05, gamma=0.00125, delta=0.05, pca=false, mode=false, convexhullApproximation=true,
              maxNumVerticesPerCH=64, minVolumePerCH=0.0001)

Decomposes a 3D surface into a set of near convex solids.
Calls external executable V-HACD based on https://github.com/kmammou/v-hacd

execute this function once before usage

filename ... Input file, OBJ Files (.obj)
outputDirectory ... Output directory
outputBasename ... Base name for output file (a number and file extension will be added automatically)
showLog ... show log file upon completion
resolution ... Maximum number of voxels generated during the voxelization stage (default=100,000, range=10,000-16,000,000)
depth ... Maximum number of clipping stages. During each split stage, solids with a concavity higher than the user defined threshold are clipped according the 'best' clipping plane (default=20, range=1-32)
concavity ... Maximum allowed concavity (default=0.0025, range=0.0-1.0)
planeDownsampling ... Controls the granularity of the search for the 'best' clipping plane (default=4, range=1-16)
convexhullDownsampling ... Controls the precision of the convex-hull generation process during the clipping plane selection stage (default=4, range=1-16)
alpha ... Controls the bias toward clipping along symmetry planes (default=0.05, range=0.0-1.0)
beta ... Controls the bias toward clipping along revolution axes (default=0.05, range=0.0-1.0)
gamma ... Controls the maximum allowed concavity during the merge stage (default=0.00125, range=0.0-1.0)
delta ... Controls the bias toward maximaxing local concavity (default=0.05, range=0.0-1.0)
pca ... Enable/disable normalizing the mesh before applying the convex decomposition (default=false)
mode ... false: voxel-based approximate convex decomposition, true: tetrahedron-based approximate convex decomposition
maxNumVerticesPerCH ... Controls the maximum number of triangles per convex-hull (default=64, range=4-1024)
minVolumePerCH ... Controls the adaptive sampling of the generated convex-hulls (default=0.0001, range=0.0-0.01)
convexhullApproximation ... Enable/disable approximation when computing convex-hulls
"""
function VHACD(filename::String, concaveProperties::ConcaveProperties)
  @assert(isfile(filename))
  @assert(filename[end-3:end] == ".obj")
  @assert(isinteger(concaveProperties.resolution))
  @assert(concaveProperties.resolution>=10000 && concaveProperties.resolution<=16000000)
  @assert(isinteger(concaveProperties.depth))
  @assert(concaveProperties.depth>=1 && concaveProperties.depth<=32)
  @assert(concaveProperties.concavity>=0.0 && concaveProperties.concavity<=1.0)
  @assert(isinteger(concaveProperties.planeDownsampling))
  @assert(concaveProperties.planeDownsampling>=1 && concaveProperties.planeDownsampling<=16)
  @assert(isinteger(concaveProperties.convexhullDownsampling))
  @assert(concaveProperties.convexhullDownsampling>=1 && concaveProperties.convexhullDownsampling<=16)
  @assert(concaveProperties.alpha>=0.0 && concaveProperties.alpha<=1.0)
  @assert(concaveProperties.beta>=0.0 && concaveProperties.beta<=1.0)
  @assert(concaveProperties.gamma>=0.0 && concaveProperties.gamma<=1.0)
  @assert(concaveProperties.delta>=0.0 && concaveProperties.delta<=1.0)
  @assert(isinteger(concaveProperties.maxNumVerticesPerCH))
  @assert(concaveProperties.maxNumVerticesPerCH>=4 && concaveProperties.maxNumVerticesPerCH<=1024)
  @assert(concaveProperties.minVolumePerCH>=0.0 && concaveProperties.minVolumePerCH<=0.01)

  logfilename = tempname()

  input = "--input \"" * filename * "\""
  output = "--output " * concaveProperties.outputBasename
  logfilename = "--log " * logfilename
  resolution = "--resolution " * string(concaveProperties.resolution)
  depth = "--depth " * string(concaveProperties.depth)
  concavity = "--concavity " * string(concaveProperties.concavity)
  planeDownsampling = "--planeDownsampling " * string(concaveProperties.planeDownsampling)
  convexhullDownsampling = "--convexhullDownsampling " * string(concaveProperties.convexhullDownsampling)
  alpha = "--alpha " * string(concaveProperties.alpha)
  beta = "--beta " * string(concaveProperties.beta)
  gamma = "--gamma " * string(concaveProperties.gamma)
  delta = "--detlta " * string(concaveProperties.delta)
  pca = "--pca " * (concaveProperties.pca ? "1" : "0")
  mode = "--mode " * (concaveProperties.mode ? "1" : "0")
  maxNumVerticesPerCH = "--maxNumVerticesPerCH " *string(concaveProperties.maxNumVerticesPerCH)
  minVolumePerCH = "--minVolumePerCH " * string(concaveProperties.minVolumePerCH)
  convexhullApproximation = "--convexhullApproximation " * (concaveProperties.convexhullApproximation ? "1" : "0")

  parameters = input * " " * output * " " * logfilename * " " * resolution * " " * depth * " " * concavity * " " * planeDownsampling * " " * convexhullDownsampling * " " * alpha * " " * beta * " " * gamma * " " * delta * " " * pca * " " * mode * " " * maxNumVerticesPerCH * " " * minVolumePerCH * " " * convexhullApproximation
  #println(parameters)

  if isdir(concaveProperties.outputDirectory) != true
    mkdir(concaveProperties.outputDirectory)
  end

  if isempty(readdir(concaveProperties.outputDirectory))   # outputDirectory is empty
    println("Dieser Ordner ist leer")
  else          # outputDirectory isn't empty
  #  println("The output directory: ", shape.concaveProperties.outputDirectory ," is not empty, please check if the concave object: ", splitdir(shape.filename)[2]," is already fragmented into disjoint convex objects, otherwise delete the content of the output directory and retry.")
  end
end


text:
using LoopVectorization
using LinearAlgebra: BlasInt, BlasFloat, LU, UnitLowerTriangular, ldiv!, checknonsingular, BLAS, LinearAlgebra

function lu(A::AbstractMatrix, pivot::Union{Val{false}, Val{true}} = Val(true); kwargs...)
    return lu!(copy(A), pivot; kwargs...)
end

function lu!(A, pivot::Union{Val{false}, Val{true}} = Val(true); check=true, kwargs...)
    m, n  = size(A)
    minmn = min(m, n)
    F = if minmn < 10 # avx introduces small performance degradation
        LinearAlgebra.generic_lufact!(A, pivot; check=check)
    else
        lu!(A, Vector{BlasInt}(undef, minmn), pivot; check=check, kwargs...)
    end
    return F
end

const RECURSION_THRESHOLD = Ref(-1)

# AVX512 needs a smaller recursion limit
function pick_threshold()
    blasvendor = BLAS.vendor()
    RECURSION_THRESHOLD[] >= 0 && return RECURSION_THRESHOLD[]
    if blasvendor === :openblas || blasvendor === :openblas64
        LoopVectorization.VectorizationBase.AVX512F ? 110 : 72
    else
        LoopVectorization.VectorizationBase.AVX512F ? 48 : 72
    end
end

function lu!(A::AbstractMatrix{T}, ipiv::AbstractVector{<:Integer},
             pivot::Union{Val{false}, Val{true}} = Val(true);
             check::Bool=true,
             # the performance is not sensitive wrt blocksize, and 16 is a good default
             blocksize::Integer=16,
             threshold::Integer=pick_threshold()) where T
    info = zero(BlasInt)
    m, n = size(A)
    mnmin = min(m, n)
    if A isa StridedArray && mnmin > threshold
        info = reckernel!(A, pivot, m, mnmin, ipiv, info, blocksize)
        if m < n # fat matrix
            # [AL AR]
            AL = @view A[:, 1:m]
            AR = @view A[:, m+1:n]
            apply_permutation!(ipiv, AR)
            ldiv!(UnitLowerTriangular(AL), AR)
        end
    else # generic fallback
        info = _generic_lufact!(A, pivot, ipiv, info)
    end
    check && checknonsingular(info)
    LU{T, typeof(A)}(A, ipiv, info)
end

function nsplit(::Type{T}, n) where T
    k = 512 ÷ (isbitstype(T) ? sizeof(T) : 8)
    k_2 = k ÷ 2
    return n >= k ? ((n + k_2) ÷ k) * k_2 : n ÷ 2
end

Base.@propagate_inbounds function apply_permutation!(P, A)
    for i in axes(P, 1)
        i′ = P[i]
        i′ == i && continue
        @simd for j in axes(A, 2)
            tmp = A[i, j]
            A[i, j] = A[i′, j]
            A[i′, j] = tmp
        end
    end
    nothing
end

function reckernel!(A::AbstractMatrix{T}, pivot::Val{Pivot}, m, n, ipiv, info, blocksize)::BlasInt where {T,Pivot}
    @inbounds begin
        if n <= max(blocksize, 1)
            info = _generic_lufact!(A, pivot, ipiv, info)
            return info
        end
        n1 = nsplit(T, n)
        n2 = n - n1
        m2 = m - n1

        # ======================================== #
        # Now, our LU process looks like this
        # [ P1 ] [ A11 A21 ]   [ L11 0 ] [ U11 U12  ]
        # [    ] [         ] = [       ] [          ]
        # [ P2 ] [ A21 A22 ]   [ L21 I ] [ 0   A′22 ]
        # ======================================== #

        # ======================================== #
        # Partition the matrix A
        # [AL AR]
        AL = @view A[:, 1:n1]
        AR = @view A[:, n1+1:n]
        #  AL  AR
        # [A11 A12]
        # [A21 A22]
        A11 = @view A[1:n1, 1:n1]
        A12 = @view A[1:n1, n1+1:n]
        A21 = @view A[n1+1:m, 1:n1]
        A22 = @view A[n1+1:m, n1+1:n]
        # [P1]
        # [P2]
        P1 = @view ipiv[1:n1]
        P2 = @view ipiv[n1+1:n]
        # ========================================

        #   [ A11 ]   [ L11 ]
        # P [     ] = [     ] U11
        #   [ A21 ]   [ L21 ]
        info = reckernel!(AL, pivot, m, n1, P1, info, blocksize)
        # [ A12 ]    [ P1 ] [ A12 ]
        # [     ] <- [    ] [     ]
        # [ A22 ]    [ 0  ] [ A22 ]
        Pivot && apply_permutation!(P1, AR)
        # A12 = L11 U12  =>  U12 = L11 \ A12
        ldiv!(UnitLowerTriangular(A11), A12)
        # Schur complement:
        # We have A22 = L21 U12 + A′22, hence
        # A′22 = A22 - L21 U12
        #mul!(A22, A21, A12, -one(T), one(T))
        schur_complement!(A22, A21, A12)
        # record info
        previnfo = info
        # P2 A22 = L22 U22
        info = reckernel!(A22, pivot, m2, n2, P2, info, blocksize)
        # A21 <- P2 A21
        Pivot && apply_permutation!(P2, A21)

        info != previnfo && (info += n1)
        @avx for i in 1:n2
            P2[i] += n1
        end
        return info
    end # inbounds
end

function schur_complement!(𝐂, 𝐀, 𝐁)
    @avx for m ∈ 1:size(𝐀,1), n ∈ 1:size(𝐁,2)
        𝐂ₘₙ = zero(eltype(𝐂))
        for k ∈ 1:size(𝐀,2)
            𝐂ₘₙ -= 𝐀[m,k] * 𝐁[k,n]
        end
        𝐂[m,n] = 𝐂ₘₙ + 𝐂[m,n]
    end
end

#=
    Modified from https://github.com/JuliaLang/julia/blob/b56a9f07948255dfbe804eef25bdbada06ec2a57/stdlib/LinearAlgebra/src/lu.jl
    License is MIT: https://julialang.org/license
=#
function _generic_lufact!(A, ::Val{Pivot}, ipiv, info) where Pivot
    m, n = size(A)
    minmn = length(ipiv)
    @inbounds begin
        for k = 1:minmn
            # find index max
            kp = k
            if Pivot
              amax = abs(zero(eltype(A)))
              for i = k:m
                  absi = abs(A[i,k])
                  if absi > amax
                      kp = i
                      amax = absi
                  end
              end
            end
            ipiv[k] = kp
            if !iszero(A[kp,k])
                if k != kp
                    # Interchange
                    @simd for i = 1:n
                        tmp = A[k,i]
                        A[k,i] = A[kp,i]
                        A[kp,i] = tmp
                    end
                end
                # Scale first column
                Akkinv = inv(A[k,k])
                @avx for i = k+1:m
                    A[i,k] *= Akkinv
                end
            elseif info == 0
                info = k
            end
            # Update the rest
            @avx for j = k+1:n
                for i = k+1:m
                    A[i,j] -= A[i,k]*A[k,j]
                end
            end
        end
    end
    return info
end

ext: jl


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

__precompile__(true)

module ValueShapes

using ArraysOfArrays
using ElasticArrays
using FillArrays
using Requires

import TypedTables

include("value_shape.jl")
include("value_accessor.jl")
include("scalar_shape.jl")
include("array_shape.jl")
include("const_value_shape.jl")
include("named_tuple_shape.jl")

function __init__()
    @require Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" include("distributions_support.jl")
end

end # module


text:
using Documenter
using Neet

makedocs(
    sitename = "Neet",
    format = :html,
    modules = [Neet],
    pages = [
        "index.md",
        "core.md"
    ]
)

deploydocs(
    repo = "github.com/dglmoore/Neet.jl.git"
)


text:
function parentind(ind::Symbol, op::Operation)
    for (id,opp) ∈ enumerate(parents(op))
        name(opp) === ind && return id
    end
    -1
end
function symbolind(ind::Symbol, op::Operation, td::UnrollArgs)
    id = parentind(ind, op)
    id == -1 && return ind, op
    @unpack u₁, u₁loopsym, u₂loopsym, vloopsym, u₂max, suffix = td
    parent = parents(op)[id]
    pvar, u₁op, u₂op = variable_name_and_unrolled(parent, u₁loopsym, u₂loopsym, vloopsym, suffix)
    Symbol(pvar, '_', Core.ifelse(u₁op, u₁, 1)), parent
end

staticexpr(x::Int) = Expr(:call, Expr(:curly, lv(:Static), x))
staticexpr(x::Union{Symbol,Expr}) = Expr(:call, lv(:Static), x)
maybestatic(x::Int) = staticexpr(x)
maybestatic(x::Union{Symbol,Expr}) = x

_MMind(ind::Union{Symbol,Expr}, str::Int) = Expr(:call, lv(:MM), VECTORWIDTHSYMBOL, ind, staticexpr(str))
_MMind(ind::Int, str::Int) = Expr(:call, lv(:MM), VECTORWIDTHSYMBOL, staticexpr(ind), staticexpr(str))

_MMind(ind::Union{Int,Symbol,Expr}, str::Union{Symbol,Expr}) = addexpr(mulexpr(_MMind(0,1),str),ind)
_MMind(ind, str::MaybeKnown) = isknown(str) ? _MMind(ind, gethint(str)) : _MMind(ind, getsym(str))

_iszero(::Symbol) = false
_iszero(::Expr) = false
_iszero(x) = iszero(x)
_isone(::Symbol) = false
_isone(::Expr) = false
_isone(x) = isone(x)
# To make this code type stable, despite potentially pushing `Int`s, `Symbol`s, or `Expr`s into `ret`,
# it takes the form of nested calls, each branching and popping off an argument.
function addoffset!(ret::Expr, stride, ind) # 3 args
    if _iszero(stride)
        pushexpr!(ret, ind)
    else
        pushexpr!(ret, _MMind(ind, stride))
    end
    nothing
end
# dropping `calcbypointeroffset` has to be delayed until after multiplying the `indexstride` by `index.
function addoffset!(ret::Expr, stride, ind, offset, calcbypointeroffset::Bool) # 5 -> 3 args
    if calcbypointeroffset
        addoffset!(ret, stride, offset)
    elseif _iszero(offset)
        addoffset!(ret, stride, ind)
    else
        addoffset!(ret, stride, addexpr(ind, offset))
    end
end
function _addoffset!(ret::Expr, vloopstride, indexstride, index, offset, calcbypointeroffset::Bool) # 6 -> 5 args
    if _isone(indexstride)
        addoffset!(ret, vloopstride, index, offset, calcbypointeroffset)
    elseif isknown(vloopstride) & isknown(indexstride)
        addoffset!(ret, gethint(vloopstride)*gethint(indexstride), index, offset, calcbypointeroffset)
    else
        addoffset!(ret, mulexpr(vloopstride,indexstride), index, offset, calcbypointeroffset)
    end
end
# multiply `index` by `indexstride`
function addoffset!(ret::Expr, vloopstride, indexstride, index, offset, calcbypointeroffset::Bool) # 6 -> (5 or 6) args
    if _isone(indexstride)
        addoffset!(ret, vloopstride, index, offset, calcbypointeroffset) # 5
    elseif calcbypointeroffset # `ind` is getting dropped, no need to allocate via `mulexpr`
        _addoffset!(ret, vloopstride, indexstride, index, offset, calcbypointeroffset) # 6
    else # multiply index by stride
        _addoffset!(ret, vloopstride, indexstride, mulexpr(index, indexstride), offset, calcbypointeroffset) # 6
    end
end

function addoffset!(ret::Expr, indvectorized::Bool, vloopstride, indexstride, index, offset, calcbypointeroffset::Bool) # 7 -> (5 or 6) args
    if indvectorized
        addoffset!(ret, vloopstride, indexstride, index, offset, calcbypointeroffset)
    elseif _isone(indexstride)
        addoffset!(ret, 0, index, offset, calcbypointeroffset)
    else
        addoffset!(ret, 0, lazymulexpr(index, indexstride), offset, calcbypointeroffset)
    end
end

function addvectoroffset!(ret::Expr, indvectorized::Bool, unrolledsteps, vloopstride, indexstride, index, offset, calcbypointeroffset::Bool) # 8 -> 7 args
    # if _iszero(unrolledsteps) # if no steps, pass through; should be unreachable
    #     addoffset!(ret, indvectorized, vloopstride, indexstride, index, offset, calcbypointeroffset)
    # else
    if calcbypointeroffset # if we previously would've dropped the index, we now replace it with the `VECTORWIDTH` step
        if _isone(unrolledsteps)
            addoffset!(ret, indvectorized, vloopstride, indexstride, VECTORWIDTHSYMBOL, offset, false)
        else
            addoffset!(ret, indvectorized, vloopstride, indexstride, mulexpr(VECTORWIDTHSYMBOL,unrolledsteps), offset, false)
        end
    elseif _isone(unrolledsteps) # add the step to the index
        addoffset!(ret, indvectorized, vloopstride, indexstride, addexpr(VECTORWIDTHSYMBOL,index), offset, false)
    else
        addoffset!(ret, indvectorized, vloopstride, indexstride, addexpr(mulexpr(VECTORWIDTHSYMBOL,unrolledsteps),index), offset, false)
    end
end
# unrolledloopstride is a stride multiple on `unrolledsteps`
function addvectoroffset!(
    ret::Expr, mm::Bool, unrolledsteps::Int, unrolledloopstride, vloopstride, indexstride::Integer, index, offset::Integer, calcbypointeroffset::Bool, indvectorized::Bool
) # 10 -> (7 or 8) args
    if unrolledsteps == 0 # neither unrolledloopstride or indexstride can be 0
        addoffset!(ret, mm, vloopstride, indexstride, index, offset, calcbypointeroffset) # 7 arg
    elseif indvectorized
        unrolledsteps *= indexstride
        if isknown(unrolledloopstride)
            addvectoroffset!(ret, mm, gethint(unrolledloopstride)*unrolledsteps, vloopstride, indexstride, index, offset, calcbypointeroffset) # 8 arg
        elseif unrolledsteps == 1
            addvectoroffset!(ret, mm, unrolledloopstride, vloopstride, indexstride, index, offset, calcbypointeroffset) # 8 arg
        else
            addvectoroffset!(ret, mm, mulexpr(unrolledloopstride,unrolledsteps), vloopstride, indexstride, index, offset, calcbypointeroffset) # 8 arg
        end
    else
        addoffset!(ret, mm, vloopstride, indexstride, index, offset + unrolledsteps, calcbypointeroffset) # 7 arg
    end
end

"""
unrolled loads are calculated as offsets with respect to an initial gesp. This has proven important to helping LLVM generate efficient code in some cases.
Therefore, unrolled === true results in inds being ignored.
_mm means to insert `mm`s.
"""
function mem_offset(op::Operation, td::UnrollArgs, inds_calc_by_ptr_offset::Vector{Bool}, _mm::Bool)
    # @assert accesses_memory(op) "Computing memory offset only makes sense for operations that access memory."
    ret = Expr(:tuple)
    indices = getindicesonly(op)
    offsets = getoffsets(op)
    strides = getstrides(op)
    loopedindex = op.ref.loopedindex
    # inds_calc_by_ptr_offset = indices_calculated_by_pointer_offsets(ls, op.ref)
    @unpack vloopsym = td
    for (n,ind) ∈ enumerate(indices)
        indvectorized = _mm & (ind === vloopsym)
        offset = offsets[n] % Int
        stride = strides[n] % Int
        ind_by_offset = inds_calc_by_ptr_offset[n] | (ind === CONSTANTZEROINDEX)
        if !ind_by_offset
            offset += (stride - 1)
        end
        @unpack vstep = td
        if loopedindex[n]
            addoffset!(ret, indvectorized, vstep, stride, ind, offset, ind_by_offset) # 7 arg
        else
            offset -= 1
            newname, parent = symbolind(ind, op, td)
            # _mmi = indvectorized && parent !== op && (!isvectorized(parent))
            # addoffset!(ret, newname, stride, offset, _mmi)
            _mmi = indvectorized && parent !== op && (!isvectorized(parent))
            @assert !_mmi "Please file an issue with an example of how you got this."
            if isu₁unrolled(parent) & (td.u₁ > 1)
                gf = GlobalRef(Core,:getfield)
                firstnew = Expr(:call, gf, Expr(:call, gf, newname, 1), 1, false)
                if isvectorized(parent) & (!_mm)
                    firstnew = Expr(:call, lv(:unmm), firstnew)
                end
                addoffset!(ret, 0, firstnew, offset, false)
            else
                if isvectorized(parent) & (!_mm)
                    addoffset!(ret, 0, Expr(:call, lv(:unmm), newname), offset, false)
                else
                    addoffset!(ret, 0, newname, offset, false)
                end
            end
        end
    end
    ret
end
# function unrolled_curly(op::Operation, u₁::Int, u₁loopsym::Symbol, vectorized::Symbol, mask::Bool)

# interleave: `0` means `false`, positive means literal, negative means multiplier
function unrolled_curly(op::Operation, u₁::Int, u₁loop::Loop, vloop::Loop, mask::Bool, interleave::Int=0)
    u₁loopsym = u₁loop.itersymbol
    vloopsym = vloop.itersymbol
    indices = getindicesonly(op)
    vstep = step(vloop)
    li = op.ref.loopedindex
    # @assert all(loopedindex)
    # @unpack u₁, u₁loopsym, vloopsym = td
    # @show vptr(op), inds_calc_by_ptr_offset
    # isone(u₁) && return mem_offset_u(op, td, inds_calc_by_ptr_offset, true)
    AV = AU = -1
    for (n,ind) ∈ enumerate(indices)
        # @show AU, op, n, ind, vloopsym, u₁loopsym
        if li[n]
            if ind === vloopsym
                @assert AV == -1 # FIXME: these asserts should be replaced with checks that prevent using `unrolled_curly` in these cases (also to be reflected in cost modeling, to avoid those)
                AV = n
            end
            if ind === u₁loopsym
                @assert AU == -1
                AU = n
            end
        else
            opp = findop(parents(op), ind)
            # @show opp
            if isvectorized(opp)
                @assert AV == -1
                AV = n
            end
            if (u₁loopsym === CONSTANTZEROINDEX) ? (CONSTANTZEROINDEX ∈ loopdependencies(opp)) : (isu₁unrolled(opp))
                @assert AU == -1
                AU = n
            end
        end
    end
    # if AU == -1
    #     return mem_offset_u(op, td, inds_calc_by_ptr_offset, true)
    # end
    # ind = mem_offset_u(op, td, inds_calc_by_ptr_offset, false)
    vecnotunrolled = AU != AV
    conditional_memory_op = isconditionalmemop(op)
    if mask || conditional_memory_op
        M = one(UInt)
        # `isu₁unrolled(last(parents(op)))` === is condop unrolled?
        # isu₁unrolled(last(parents(op)))
        if vecnotunrolled || conditional_memory_op || (interleave > 0) # mask all
            M = (M << u₁) - M
        else # mask last
            M <<= (u₁ - 1)
        end
    else
        M = zero(UInt)
    end
    @assert isknown(step(u₁loop)) "Unrolled loops must have known steps to use `Unroll` type; this is a bug, shouldn't have reached here"
    if AV > 0
        @assert isknown(step(vloop)) "Vectorized loops must have known steps to use `Unroll` type; this is a bug, shouldn't have reached here."
        X = convert(Int, getstrides(op)[AV])
        X *= gethint(step(vloop))
        intvecsym = :(Int($VECTORWIDTHSYMBOL))
        if interleave > 0
            Expr(:curly, lv(:Unroll), AU, interleave, u₁, AV, intvecsym, M, X)
        elseif interleave < 0
            unrollstepexpr = :(Int($(mulexpr(VECTORWIDTHSYMBOL, -interleave))))
            Expr(:curly, lv(:Unroll), AU, unrollstepexpr, u₁, AV, intvecsym, M, X)
        else
            if vecnotunrolled
                # Expr(:call, Expr(:curly, lv(:Unroll), AU, 1, u₁, AV, intvecsym, M, 1), ind)
                Expr(:curly, lv(:Unroll), AU, gethint(step(u₁loop)), u₁, AV, intvecsym, M, X)
            else
                if isone(X)
                    Expr(:curly, lv(:Unroll), AU, intvecsym, u₁, AV, intvecsym, M, X)
                else
                    unrollstepexpr = :(Int($(mulexpr(VECTORWIDTHSYMBOL, X))))
                    Expr(:curly, lv(:Unroll), AU, unrollstepexpr, u₁, AV, intvecsym, M, X)
                end
            end
        end
    else
        Expr(:curly, lv(:Unroll), AU, gethint(step(u₁loop)), u₁, 0, 1, M, 1)
    end
end
function unrolledindex(op::Operation, td::UnrollArgs, mask::Bool, inds_calc_by_ptr_offset::Vector{Bool})
    @unpack u₁, u₁loopsym, u₁loop, vloop = td
    isone(u₁) && return mem_offset_u(op, td, inds_calc_by_ptr_offset, true)
    any(==(u₁loopsym), getindicesonly(op)) || return mem_offset_u(op, td, inds_calc_by_ptr_offset, true)

    unrollcurl = unrolled_curly(op, u₁, u₁loop, vloop, mask)
    ind = mem_offset_u(op, td, inds_calc_by_ptr_offset, false)
    Expr(:call, unrollcurl, ind)
end

function mem_offset_u(op::Operation, td::UnrollArgs, inds_calc_by_ptr_offset::Vector{Bool}, _mm::Bool, incr₁::Int = 0)
    @assert accesses_memory(op) "Computing memory offset only makes sense for operations that access memory."
    @unpack u₁loopsym, u₂loopsym, vloopsym, u₁step, u₂step, vstep, suffix = td

    #u₁
    incr₂ = max(suffix, 0)#(suffix == -1) ? 0 : suffix
    ret = Expr(:tuple)
    indices = getindicesonly(op)
    offsets = getoffsets(op)
    strides = getstrides(op)
    # allbasezero = all(inds_calc_by_ptr_offset) && all(iszero, offsets)
    loopedindex = op.ref.loopedindex
    if iszero(incr₁) & iszero(incr₂)
        return mem_offset(op, td, inds_calc_by_ptr_offset, _mm)
        # append_inds!(ret, indices, loopedindex)
    else
        for (n,ind) ∈ enumerate(indices)
            ind_by_offset = inds_calc_by_ptr_offset[n] | (ind === CONSTANTZEROINDEX)
            offset = convert(Int, offsets[n])
            stride = convert(Int, strides[n])
            indvectorized = ind === vloopsym
            indvectorizedmm = _mm & indvectorized
            if !ind_by_offset
                offset += (stride - 1)
            end
            if ind === u₁loopsym
                addvectoroffset!(ret, indvectorizedmm, incr₁, u₁step, vstep, stride, ind, offset, ind_by_offset, indvectorized) # 9 arg
            elseif ind === u₂loopsym
                # if isstore(op)
                #     @show indvectorized, ind === vloopsym, u₂loopsym, incr₂
                # end
                addvectoroffset!(ret, indvectorizedmm, incr₂, u₂step, vstep, stride, ind, offset, ind_by_offset, indvectorized) # 9 arg
            elseif loopedindex[n]
                addoffset!(ret, indvectorizedmm, vstep, stride, ind, offset, ind_by_offset) # 7 arg
            else
                offset -= 1
                newname, parent = symbolind(ind, op, td)
                _mmi = indvectorizedmm && parent !== op && (!isvectorized(parent))
                #                              addoffset!(ret, newname, 1, offset, _mmi)
                @assert !_mmi "Please file an issue with an example of how you got this."
                if isvectorized(parent) & (!_mm)
                    if isu₁unrolled(parent) & (td.u₁ > 1)
                        gf = GlobalRef(Core,:getfield)
                        newname_unmm = Expr(:call, lv(:unmm), Expr(:call, gf, Expr(:call, gf, newname, 1), incr₁+1, false))
                    else
                        newname_unmm = Expr(:call, lv(:unmm), newname)
                    end
                    if stride ≠ 1
                        newname_unmm = mulexpr(newname_unmm,stride)
                    end
                    addoffset!(ret, 0, newname_unmm, offset, false)
                elseif isu₁unrolled(parent) & (td.u₁ > 1)
                    gf = GlobalRef(Core,:getfield)
                    firstnew = Expr(:call, gf, Expr(:call, gf, newname, 1), incr₁+1, false)
                    if stride ≠ 1
                        firstnew = mulexpr(firstnew,stride)
                    end
                    addoffset!(ret, 0, firstnew, offset, false)
                elseif stride == 1
                    addoffset!(ret, 0, newname, offset, false)
                else
                    addoffset!(ret, 0, mulexpr(newname,stride), offset, false)
                end
            end
        end
    end
    ret
end

@inline and_last(a, b) = a & b
@generated function and_last(v::VecUnroll{N}, m) where {N}
    q = Expr(:block, Expr(:meta,:inline), :(vd = data(v)))
    t = Expr(:call, lv(:promote))
    gf = GlobalRef(Core, :getfield)
    for n ∈ 1:N
        push!(t.args, :($gf(vd, $n, false)))
    end
    push!(t.args, :($gf(vd, $(N+1), false) & m))
    push!(q.args, Expr(:call, lv(:VecUnroll), t))
    q
end

isconditionalmemop(op::Operation) = (instruction(op).instr === :conditionalload) || (instruction(op).instr === :conditionalstore!)
function add_memory_mask!(memopexpr::Expr, op::Operation, td::UnrollArgs, mask::Bool)
    @unpack u₁, u₁loopsym, u₂loopsym, vloopsym, u₂max, suffix = td
    if isconditionalmemop(op)
        condop = last(parents(op))
        opu₂ = (suffix ≠ -1) && isu₂unrolled(op)
        condvar, condu₁unrolled = condvarname_and_unroll(condop, u₁loopsym, u₂loopsym, vloopsym, suffix, opu₂)
        # if it isn't unrolled, then `m`
        u = condu₁unrolled ? u₁ : 1
        # u = isu₁unrolled(condop) ? u₁ : 1
        condvar = Symbol(condvar, '_', u)
        # If we need to apply `MASKSYMBOL` and the condvar
        # 2 condvar possibilities:
        #   `VecUnroll` applied everywhere
        #    single mask "broadcast"
        # 2 mask possibilities
        #    u₁loopsym ≠  vloopsym, and we mask all
        #    u₁loopsym == vloopsym, and we mask last
        # broadcast both, so can do so implicitly
        # this is true whether or not `condbroadcast`
        if !mask || (!isvectorized(op))
            push!(memopexpr.args, condvar)
        elseif (u₁loopsym ≢ vloopsym) | (u₁ == 1) # mask all equivalenetly
            push!(memopexpr.args, Expr(:call, lv(:&), condvar, MASKSYMBOL))
            # if the condition `(u₁loopsym ≢ vloopsym) | (u₁ == 1)` failed, we need to apply `MASKSYMBOL` only to last unroll.
        elseif !condu₁unrolled && isu₁unrolled(op) # condbroadcast
            # explicitly broadcast `condvar`, and apply `MASKSYMBOL` to end
            t = Expr(:call, lv(:promote))
            for um ∈ 1:u₁-1
                push!(t.args, condvar)
            end
            push!(t.args, Expr(:call, lv(:&), condvar, MASKSYMBOL))
            push!(memopexpr.args, Expr(:call, lv(:VecUnroll), t))
        else# !condbroadcast && !vecunrolled
            push!(memopexpr.args, Expr(:call, lv(:and_last), condvar, MASKSYMBOL))
        end
    elseif mask && isvectorized(op)
        push!(memopexpr.args, MASKSYMBOL)
    end
    nothing
end

varassignname(var::Symbol, u::Int, isunrolled::Bool) = isunrolled ? Symbol(var, u) : var
# name_memoffset only gets called when vectorized
function name_memoffset(var::Symbol, op::Operation, td::UnrollArgs, u₁unrolled::Bool, inds_calc_by_ptr_offset::Vector{Bool})
    @unpack u₁, u₁loopsym, u₂loopsym, suffix = td
    if (suffix == -1) && u₁ < 0 # u₁ == -1 sentinel value meaning not unrolled
        name = var
        mo = mem_offset(op, td, inds_calc_by_ptr_offset, true)
    else
        name = u₁unrolled ? Symbol(var, u₁) : var
        mo = mem_offset_u(op, td, inds_calc_by_ptr_offset)
    end
    name, mo
end

function condvarname_and_unroll(cond::Operation, u₁loop::Symbol, u₂loop::Symbol, vloop::Symbol, suffix::Int, opu₂::Bool)
    condvar, condu₁, condu₂ = variable_name_and_unrolled(cond, u₁loop, u₂loop, vloop, Core.ifelse(opu₂, suffix, -1))
    condvar, condu₁
end


text:
using DocumenterTravis
using Test

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


text:
"""
    struct GenericMesh{N,T} <: AbstractMesh{N,T}

Data structure representing a generic mesh in an ambient space of dimension `N`,
with data of type `T`.
"""
Base.@kwdef struct GenericMesh{N,T} <: AbstractMesh{N,T}
    nodes::Vector{SVector{N,T}} = Vector{SVector{N,T}}()
    # for each element type (key), store the data required to reconstruct the
    # elements (value).
    elements::Dict{DataType,Any} = Dict{DataType,Any}()
    # mapping from entity to (etype,tags)
    ent2tags::Dict{AbstractEntity,Dict{DataType,Vector{Int}}} = Dict{AbstractEntity,Dict{DataType,Vector{Int}}}()
end

nodes(m::GenericMesh)    = m.nodes

elements(m::GenericMesh) = m.elements

ent2tags(m::GenericMesh) = m.ent2tags

Base.keys(m::GenericMesh) = keys(elements(m))

Geometry.entities(m::GenericMesh)  = keys(ent2tags(m)) |> collect

domain(m::GenericMesh)    = entities(m) |> Domain

# implement the interface for ElementIterator of lagrange elements on a generic mesh
function Base.size(iter::ElementIterator{<:LagrangeElement,<:GenericMesh})
    msh               = mesh(iter)
    E                 = eltype(iter)
    tags::Matrix{Int} = msh.elements[E]
    _, Nel           = size(tags)
    return (Nel,)
end

function Base.getindex(iter::ElementIterator{<:LagrangeElement,<:GenericMesh},i::Int)
    E                   = eltype(iter)
    M                   = mesh(iter)
    tags::Matrix{Int}   = M.elements[E]
    node_tags           = view(tags,:,i)
    vtx                 = view(M.nodes,node_tags)
    el                  = E(vtx)
    return el
end

function Base.iterate(iter::ElementIterator{<:LagrangeElement,<:GenericMesh}, state=1)
    state > length(iter) && (return nothing)
    iter[state], state + 1
end

"""
    dom2elt(m::GenericMesh,Ω,E)

Compute the element indices `idxs` of the elements of type `E` composing `Ω`, so that `m[E][idxs]` gives all
the elements of type `E` meshing `Ω`.
"""
function dom2elt(m::GenericMesh, Ω, E::DataType)
    idxs = Int[]
    for ent in entities(Ω)
        tags = get(m.ent2tags[ent], E, Int[])
        append!(idxs, tags)
    end
    return idxs
end

"""
    dom2elt(m::GenericMesh,Ω)

Return a `Dict` with keys being the element types of `m`, and values being the
indices of the elements in `Ω` of that type.
"""
function dom2elt(m::GenericMesh, Ω)
    dict = Dict{DataType,Vector{Int}}()
    for E in keys(m)
        tags = dom2elt(m, Ω, E)
        if !isempty(tags)
            dict[E] = tags
        end
    end
    return dict
end

# convert a mesh to 2d by ignoring third component. Note that this also requires
# converting various element types to their 2d counterpart. These are needed
# because some meshers like gmsh always create three-dimensional objects, so we
# must convert after importing the mesh
function convert_to_2d(mesh::GenericMesh{3})
    @assert all(E -> geometric_dimension(domain(E)) < 3, keys(mesh))
    T = primitive_type(mesh)
    # create new dictionaries for elements and ent2tags with 2d elements as keys
    els  = empty(elements(mesh))
    e2t  = empty(ent2tags(mesh))
    for (E, tags) in elements(mesh)
        E2d = convert_to_2d(E)
        els[E2d] = tags
    end
    for (ent, dict) in ent2tags(mesh)
        new_dict = empty(dict)
        for (E, tags) in dict
            E2d = convert_to_2d(E)
            new_dict[E2d] = tags
        end
        e2t[ent] = new_dict
    end
    # construct new 2d mesh
    GenericMesh{2,T}(;
        nodes=[x[1:2] for x in nodes(mesh)],
        elements=els,
        ent2tags=e2t
    )
end

function convert_to_2d(::Type{LagrangeElement{R,N,SVector{3,T}}}) where {R,N,T}
    LagrangeElement{R,N,SVector{2,T}}
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
using BenchmarkingEconomicEfficiency
using Test

@testset "BenchmarkingEconomicEfficiency.jl" begin

    include("deacostddf.jl")
    include("dearevenueddf.jl")
    include("deacostadd.jl")
    include("dearevenueadd.jl")
    include("deaprofitadd.jl")
    include("deacostrussell.jl")
    include("dearevenuerussell.jl")
    include("deaprofitrussell.jl")
    include("deaprofiterg.jl")
    include("deaprofitmddf.jl")
    include("deacostholder.jl")
    include("dearevenueholder.jl")
    include("deaprofitholder.jl")
    include("deacostrddf.jl")
    include("dearevenuerddf.jl")
    include("deaprofitrddf.jl")
    include("deacostgda.jl")
    include("dearevenuegda.jl")
    include("deaprofitgda.jl")

end


text:
# Belief Propagation                #
# --------------------------------- #
@rule MvNormalMeanCovariance(:μ, Marginalisation) (m_out::PointMass, m_Σ::PointMass) = MvNormalMeanCovariance(mean(m_out), mean(m_Σ))

@rule MvNormalMeanCovariance(:μ, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_Σ::PointMass) = begin
    m_out_mean, m_out_cov = mean_cov(m_out)
    return MvNormalMeanCovariance(m_out_mean, m_out_cov + mean(m_Σ))
end

# Variational                       #
# --------------------------------- #
@rule MvNormalMeanCovariance(:μ, Marginalisation) (q_out::Any, q_Σ::Any) = MvNormalMeanCovariance(mean(q_out), mean(q_Σ))

@rule MvNormalMeanCovariance(:μ, Marginalisation) (m_out::PointMass, q_Σ::Any) = MvNormalMeanCovariance(mean(m_out), mean(q_Σ))

@rule MvNormalMeanCovariance(:μ, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, q_Σ::Any) = begin
    m_out_mean, m_out_cov = mean_cov(m_out)
    return MvNormalMeanCovariance(m_out_mean, m_out_cov + mean(q_Σ))
end


text:
"""
    struct AngularBasis <: BasisSet

Defines a set of states with angular momentum `j`. The dimensionality if `2j+1`.

# Interface

* `m(b::AngularBasis, i)`: returns the ``m`` quantum number of the `i`th state.
"""
struct AngularBasis <: BasisSet
    j :: HalfInteger
    function AngularBasis(j)
        @assert j >= 0
        new(j)
    end
end
Base.length(b::AngularBasis) = convert(Int, 2 * b.j + 1)
m(b::AngularBasis, idx) = -b.j + (idx - 1)

"""
    struct JOperator{K} <: LinearOperator{AngularBasis}

Type to represent the angular momentum operators ``J_z``, ``J_+`` and ``J_-``, acting on
states with defined total angular momentum ([`AngularBasis`](@ref)).
"""
struct JOperator{K} <: LinearOperator{AngularBasis}
    b :: AngularBasis
end
function JOperator(b::AngularBasis, kind::Symbol)
    kind in [:z, :+, :-] || throw(ArgumentError("Bad kind: $kind"))
    JOperator{kind}(b)
end
basis(op::JOperator) = op.b
function apply(i::Integer, op::JOperator{K}, j::Integer) :: ComplexF64 where K
    angJ = convert(ComplexF64, op.b.j)
    if K == :z
        (i == j) ? m(op.b, j) : 0.0
    elseif K == :+
        (i == j + 1) ? sqrt((angJ - m(op.b, j)) * (angJ + m(op.b, j) + 1)) : 0.0
    elseif K == :-
        (i == j - 1) ? sqrt((angJ + m(op.b, j)) * (angJ - m(op.b, j) + 1)) : 0.0
    else
        error("Unreachable reached.")
    end
end

struct JOperatorSet{B <: BasisSet, T <: LinearOperator{B}, U <: LinearOperator{B}, V <: LinearOperator{B}}
    J₊ :: T
    J₋ :: U
    Jz :: V
end
function JOperatorSet(JOp, b::BasisSet)
    Jz = JOp(b, :z)
    J₊, J₋ = JOp(b, :+), JOp(b, :-)
    JOperatorSet{typeof(b),typeof(J₊),typeof(J₋),typeof(Jz)}(J₊, J₋, Jz)
end
JOperatorSet(b::BasisSet) = JOperatorSet(JOperator, b)

function JOperatorSet(f::Function, js::Vararg{JOperatorSet})
    length(js) > 0 || throw(ArgumentError("Must provide at least one J operator"))
    Angular.JOperatorSet(
        f((j.J₊ for j in js)...),
        f((j.J₋ for j in js)...),
        f((j.Jz for j in js)...)
    )
end

basis(op::JOperatorSet) = basis(op.Jz)

Jx(op::JOperatorSet) = (op.J₊ + op.J₋) / 2
Jy(op::JOperatorSet) = (op.J₊ - op.J₋) / (2im)
J2(op::JOperatorSet) = J2Operator(op)

"""
    findj(λ) -> HalfInteger

Convert and ``J^2`` operator eigenvalue into the corresponding half-integer ``j`` value.
"""
findj(λ) = HalfInteger(isqrt(1 + round(Int, real(4λ))) - 1, 2)

struct J2Operator{B <: BasisSet, JS <: JOperatorSet{B}, EOP <: LinearOperator{B}} <: LinearOperator{B}
    J :: JS
    J2 :: EOP

    function J2Operator(jops :: JOperatorSet)
        b = basis(jops)
        J2 = (jops.J₊*jops.J₋ + jops.J₋*jops.J₊)/2 + jops.Jz^2
        new{typeof(b),typeof(jops),typeof(J2)}(jops, J2)
    end
end
basis(op::J2Operator) = basis(op.J)
apply(i::Integer, op::J2Operator, j::Integer) = apply(i, op.J2, j)
function matrix(j2op::J2Operator)
    J₊ = matrix(j2op.J.J₊)
    J₋ = matrix(j2op.J.J₋)
    Jz = matrix(j2op.J.Jz)
    (J₊*J₋ .+ J₋*J₊)./2 .+ Jz^2
end


text:
using LazySets, MathematicalPredicates, SparseArrays

function heat_specification()
    # initial set: xᵢ ∈ [0.6, 0.625] for i ≤ 2 and xᵢ = 0 for i > 2
    X0 = Hyperrectangle(low=[fill(0.6, 2); zeros(198)],
                        high=[fill(0.625, 2); zeros(198)])

    # safety property: x133 ≤ 0.1
    property = is_contained_in(HalfSpace(sparsevec([133], [1.0], 200), 0.1))

    # time horizon: 20 time units
    time_horizon = 20.0

    # specification
    O = Dict(:T => time_horizon, :property => property)

    return X0, O
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
using Test

ProjDir = joinpath(dirname(@__FILE__), "..", "examples", "BernoulliInitTheta")
cd(ProjDir) do

  isdir("tmp") &&
    rm("tmp", recursive=true);

    bernoullimodel = "
    data {
      int<lower=1> N;
      int<lower=0,upper=1> y[N];
    }
    parameters {
      real<lower=0,upper=1> theta;
    }
    model {
      theta ~ beta(1,1);
      y ~ bernoulli(theta);
    }
    "

    bernoullidata = [
      Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1]),
      Dict("N" => 10, "y" => [0, 1, 0, 0, 0, 0, 1, 0, 0, 1]),
      Dict("N" => 10, "y" => [0, 0, 0, 0, 0, 0, 1, 0, 1, 1]),
      Dict("N" => 10, "y" => [0, 0, 0, 1, 0, 0, 0, 1, 0, 1])
    ]

    theta_init = [0.1, 0.4, 0.5, 0.9]
    bernoulliinit = Dict[
      Dict("theta" => theta_init[1]),
      Dict("theta" => theta_init[2]),
      Dict("theta" => theta_init[3]),
      Dict("theta" => theta_init[4]),
    ]

    monitor = ["theta", "lp__", "accept_stat__"]

    stanmodel = Stanmodel(name="bernoulli",
      model=bernoullimodel,
      init=Stan.Init(init=bernoulliinit),
      adapt=1);

      sim1 = stan(stanmodel, bernoullidata, CmdStanDir=CMDSTAN_HOME)
      global sim1_values = round.(sim1.value[1, end, :], 1)

  isdir("tmp") &&
    rm("tmp", recursive=true);

  sim1_values |> display

  @assert sim1_values == theta_init "Not close to theta input values."
end # cd


text:
using JuMP
using Ipopt
using Hiop
using DelimitedFiles
using LinearAlgebra
using Test

include("opfdata.jl")

# export solve, model, initialPt_IPOPT, outputAll, computeAdmitances
function solve(opfmodel, opf_data)
    optimize!(opfmodel)
    status = termination_status(opfmodel)
    if status != MOI.LOCALLY_SOLVED
        println("Could not solve the model to optimality.")
    end
    return opfmodel,status
end

function model(opf_data; max_iter=100, solver="Ipopt", hiopalgebra = :Dense, scaling = 1.0)
    Pg0, Qg0, Vm0, Va0 = initialPt_IPOPT(opf_data)
    lines = opf_data.lines; buses = opf_data.buses; generators = opf_data.generators; baseMVA = opf_data.baseMVA
    busIdx = opf_data.BusIdx; FromLines = opf_data.FromLines; ToLines = opf_data.ToLines; BusGeners = opf_data.BusGenerators;

    nbus  = length(buses); nline = length(lines); ngen  = length(generators)

    #branch admitances
    YffR,YffI,YttR,YttI,YftR,YftI,YtfR,YtfI,YshR,YshI = computeAdmitances(lines, buses, baseMVA)

    #
    # JuMP model now
    #
    if solver == "Hiop"
        opfmodel = Model(optimizer_with_attributes(Hiop.Optimizer, "algebra" => hiopalgebra))
    else
        opfmodel = Model(optimizer_with_attributes(Ipopt.Optimizer, "max_iter" => max_iter))
    end

    @variable(opfmodel, generators[i].Pmin <= Pg[i=1:ngen] <= generators[i].Pmax, start = Pg0[i])
    @variable(opfmodel, generators[i].Qmin <= Qg[i=1:ngen] <= generators[i].Qmax, start = Qg0[i])

    @variable(opfmodel, Va[i=1:nbus], start = Va0[i])
    @variable(opfmodel, buses[i].Vmin <= Vm[i=1:nbus] <= buses[i].Vmax, start = Vm0[i])
    # fix the voltage angle at the reference bus
    set_lower_bound(Va[opf_data.bus_ref], buses[opf_data.bus_ref].Va)
    set_upper_bound(Va[opf_data.bus_ref], buses[opf_data.bus_ref].Va)

    # minimize active power
    coeff0 = Vector{Float64}(undef, ngen)
    for (i,v) in enumerate(generators)
        coeff0[i] = v.coeff[v.n]
    end
    coeff1 = Vector{Float64}(undef, ngen)
    for (i,v) in enumerate(generators)
        coeff1[i] = v.coeff[v.n - 1]
    end
    coeff2 = Vector{Float64}(undef, ngen)
    for (i,v) in enumerate(generators)
        coeff2[i] = v.coeff[v.n - 2]
    end

    @NLobjective(opfmodel, Min, scaling*sum( coeff2[i]*(baseMVA*Pg[i])^2
                                    + coeff1[i]*(baseMVA*Pg[i])
                                    + coeff0[i] for i=1:ngen))

    #
    # power flow balance
    #

    for b in 1:nbus
        # real part
        @NLconstraint(
        opfmodel,
        ( sum( YffR[l] for l in FromLines[b]) + sum( YttR[l] for l in ToLines[b]) + YshR[b] ) * Vm[b]^2
        + sum( Vm[b]*Vm[busIdx[lines[l].to]]  *( YftR[l]*cos(Va[b]-Va[busIdx[lines[l].to]]  ) + YftI[l]*sin(Va[b]-Va[busIdx[lines[l].to]]  )) for l in FromLines[b] )
        + sum( Vm[b]*Vm[busIdx[lines[l].from]]*( YtfR[l]*cos(Va[b]-Va[busIdx[lines[l].from]]) + YtfI[l]*sin(Va[b]-Va[busIdx[lines[l].from]])) for l in ToLines[b]   )
        - ( sum(baseMVA*Pg[g] for g in BusGeners[b]) - buses[b].Pd ) / baseMVA      # Sbus part
        ==0)
        # imaginary part
    end
    for b in 1:nbus
        @NLconstraint(
        opfmodel,
        ( sum(-YffI[l] for l in FromLines[b]) + sum(-YttI[l] for l in ToLines[b]) - YshI[b] ) * Vm[b]^2
        + sum( Vm[b]*Vm[busIdx[lines[l].to]]  *(-YftI[l]*cos(Va[b]-Va[busIdx[lines[l].to]]  ) + YftR[l]*sin(Va[b]-Va[busIdx[lines[l].to]]  )) for l in FromLines[b] )
        + sum( Vm[b]*Vm[busIdx[lines[l].from]]*(-YtfI[l]*cos(Va[b]-Va[busIdx[lines[l].from]]) + YtfR[l]*sin(Va[b]-Va[busIdx[lines[l].from]])) for l in ToLines[b]   )
        - ( sum(baseMVA*Qg[g] for g in BusGeners[b]) - buses[b].Qd ) / baseMVA      #Sbus part
        ==0)
    end
    #
    # branch/lines flow limits
    #
    nlinelim=0
    for l in 1:nline
        if lines[l].rateA!=0 && lines[l].rateA<1.0e10
        nlinelim += 1
        flowmax=(lines[l].rateA/baseMVA)^2

        #branch apparent power limits (from bus)
        Yff_abs2=YffR[l]^2+YffI[l]^2;        Yft_abs2=YftR[l]^2+YftI[l]^2
        Yre=YffR[l]*YftR[l]+YffI[l]*YftI[l]; Yim=-YffR[l]*YftI[l]+YffI[l]*YftR[l]
        @NLconstraint(
            opfmodel,
            Vm[busIdx[lines[l].from]]^2 *
            ( Yff_abs2*Vm[busIdx[lines[l].from]]^2 + Yft_abs2*Vm[busIdx[lines[l].to]]^2
            + 2*Vm[busIdx[lines[l].from]]*Vm[busIdx[lines[l].to]]
                *(Yre*cos(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])
                -Yim*sin(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])
                )
            )
            - flowmax <=0)

        #branch apparent power limits (to bus)
        Ytf_abs2=YtfR[l]^2+YtfI[l]^2; Ytt_abs2=YttR[l]^2+YttI[l]^2
        Yre=YtfR[l]*YttR[l]+YtfI[l]*YttI[l]; Yim=-YtfR[l]*YttI[l]+YtfI[l]*YttR[l]
        @NLconstraint(
            opfmodel,
            Vm[busIdx[lines[l].to]]^2 *
            ( Ytf_abs2*Vm[busIdx[lines[l].from]]^2 + Ytt_abs2*Vm[busIdx[lines[l].to]]^2
            + 2*Vm[busIdx[lines[l].from]]*Vm[busIdx[lines[l].to]]
                *(Yre*cos(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])
                -Yim*sin(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])
                )
            )
            - flowmax <=0)
        end
    end
    println("Buses: $nbus  Lines: $nline  Generators: $ngen")
    println("Lines with limits: $nlinelim")

    return opfmodel, Pg, Qg, Va, Vm
end
# Compute initial point for IPOPT based on the values provided in the case data
function initialPt_IPOPT(opfdata)
    Pg=zeros(length(opfdata.generators)); Qg=zeros(length(opfdata.generators)); i=1
    for g in opfdata.generators
        # set the power levels in in between the bounds as suggested by matpower
        # (case data also contains initial values in .Pg and .Qg - not used with IPOPT)
        Pg[i]=0.5*(g.Pmax+g.Pmin)
        Qg[i]=0.5*(g.Qmax+g.Qmin)
        i=i+1
    end
    @assert i-1==length(opfdata.generators)

    Vm=zeros(length(opfdata.buses)); i=1;
    for b in opfdata.buses
        # set the ini val for voltage magnitude in between the bounds
        # (case data contains initials values in Vm - not used with IPOPT)
        Vm[i]=0.5*(b.Vmax+b.Vmin);
        i=i+1
    end
    @assert i-1==length(opfdata.buses)

    # set all angles to the angle of the reference bus
    Va = opfdata.buses[opfdata.bus_ref].Va * ones(length(opfdata.buses))

    return Pg,Qg,Vm,Va
end

# cases with objective scaling
# Dense Hiop algebra
for case in [[0.0283286,"case9"],[0.00357143, "case118"],[0.00645161, "case300"]]
    casepath = joinpath(dirname(@__FILE__), "data", case[2])
    max_iter=100
    opfdata = opf_loaddata(casepath)
    Pg0, Qg0, Vm0, Va0 = initialPt_IPOPT(opfdata)
    opfmodel_ref, Pg_ref, Qg_ref, Va_ref, Vm_ref = model(opfdata; solver="Ipopt", scaling = case[1])
    opfmodel, Pg, Qg, Va, Vm = model(opfdata; solver="Hiop", hiopalgebra = :Dense, scaling = case[1])
    opfmodel_ref ,status = solve(opfmodel_ref,opfdata)
    opfmodel,status = solve(opfmodel,opfdata)
    @test value.(Pg) ≈ value.(Pg)
    @test value.(Qg) ≈ value.(Qg)
    @test value.(Vm) ≈ value.(Vm)
    @test value.(Va) ≈ value.(Va)
    @test objective_value(opfmodel) ≈ objective_value(opfmodel_ref)
end

# Sparse Hiop algebra
# This works with sparse algebra
for case in [[0.0283286,"case9"],[0.00357143, "case118"]]
    casepath = joinpath(dirname(@__FILE__), "data", case[2])
    max_iter=100
    opfdata = opf_loaddata(casepath)
    Pg0, Qg0, Vm0, Va0 = initialPt_IPOPT(opfdata)
    opfmodel_ref, Pg_ref, Qg_ref, Va_ref, Vm_ref = model(opfdata; solver="Ipopt", scaling = case[1])
    opfmodel, Pg, Qg, Va, Vm = model(opfdata; solver="Hiop", hiopalgebra = :Sparse, scaling = case[1])
    opfmodel_ref ,status = solve(opfmodel_ref,opfdata)
    opfmodel,status = solve(opfmodel,opfdata)
    @test value.(Pg) ≈ value.(Pg)
    @test value.(Qg) ≈ value.(Qg)
    @test value.(Vm) ≈ value.(Vm)
    @test value.(Va) ≈ value.(Va)
    @test objective_value(opfmodel) ≈ objective_value(opfmodel_ref)
end

# Crashes
# for case in [[0.00645161, "case300"]]
#     casepath = joinpath(dirname(@__FILE__), "data", case[2])
#     max_iter=100
#     opfdata = opf_loaddata(casepath)
#     Pg0, Qg0, Vm0, Va0 = initialPt_IPOPT(opfdata)
#     opfmodel_ref, Pg_ref, Qg_ref, Va_ref, Vm_ref = model(opfdata; solver="Ipopt", scaling = case[1])
#     opfmodel, Pg, Qg, Va, Vm = model(opfdata; solver="Hiop", hiopalgebra = :Sparse, scaling = case[1])
#     opfmodel_ref ,status = solve(opfmodel_ref,opfdata)
#     opfmodel,status = solve(opfmodel,opfdata)
#     @test_broken value.(Pg) ≈ value.(Pg)
#     @test_broken value.(Qg) ≈ value.(Qg)
#     @test_broken value.(Vm) ≈ value.(Vm)
#     @test_broken value.(Va) ≈ value.(Va)
#     @test_broken objective_value(opfmodel) ≈ objective_value(opfmodel_ref)
# end


text:
import MaximumLikelihoodProblems

import LinearAlgebra
import Random
import Statistics
import Test

_this_test_filename = @__FILE__
_test_directory = dirname(_this_test_filename)
root_directory = dirname(_test_directory)
examples_directory = joinpath(root_directory, "examples")
test_directory = joinpath(root_directory, "test")

Test.@testset "MaximumLikelihoodProblems.jl" begin
    Test.@testset "linear_algebra.jl" begin
        Test.@testset "_pseudoinverse" begin
            Test.@test MaximumLikelihoodProblems._pseudoinverse([1 2; 3 4]) == LinearAlgebra.inv([1 2; 3 4])
            Test.@test MaximumLikelihoodProblems._pseudoinverse([0 0; 0 0]) == [0 0; 0 0]
        end
    end

    Test.@testset "examples" begin
        include("examples.jl")
    end
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

# TODO: radical

const SMARTS_ATOM_COND_SYMBOL = Dict(
    'X' => :connectivity,
    'D' => :nodedegree,
    'v' => :valence,
    'H' => :hcount,
    'r' => :atom_sssrsizes,
    'R' => :atom_sssrcount
)

const SMARTS_CHARGE_SIGN = Dict(
    '+' => 1,
    '-' => -1
)

function atomsymbol!(state::SmartsParserState)
    """ Atomsymbol <- Br / Cl / [AaBCcNnOoPpSsFI*]
    """
    if read(state) == 'C' && lookahead(state, 1) == 'l'
        forward!(state, 2)
        return :and => (:atomsymbol => :Cl, :isaromatic => false)
    elseif read(state) == 'B' && lookahead(state, 1) == 'r'
        forward!(state, 2)
        return :and => (:atomsymbol => :Br, :isaromatic => false)
    elseif read(state) in "BCNOPSFI"
        sym = Symbol(read(state))
        forward!(state)
        return :and => (:atomsymbol => sym, :isaromatic => false)
    elseif read(state) in "cnops"
        sym = Symbol(uppercase(read(state)))
        forward!(state)
        return :and => (:atomsymbol => sym, :isaromatic => true)
    elseif read(state) == 'A'
        forward!(state)
        return :isaromatic => false
    elseif read(state) == 'a'
        forward!(state)
        return :isaromatic => true
    elseif read(state) == '*'
        forward!(state)
        return :any => true
    end
end

function atom!(state::SmilesParser)
    """ Atom <- '[' AtomProp+ ']' / AtomSymbol
    """
    c = read(state)
    if c == '['
        forward!(state)
        a = lghighand!(state, atomprop!)
        ca = read(state)
        @assert ca == ']' "(atom!) unexpected token: $(ca)"
        forward!(state)
        prop = collectand(a)
        @assert haskey(prop, :atomsymbol) || prop[:hcount] >= 1
        atoms = [SmilesAtom(
            get(prop, :atomsymbol, :H),
            get(prop, :charge, 0),
            1,
            get(prop, :mass, nothing),
            get(prop, :isaromatic, false),
            get(prop, :stereo, :unspecified)
        )]
        hcnt = get(prop, :hcount, 0)
        if !haskey(prop, :atomsymbol)
            hcnt -= 1
        end
        for h in 1:hcnt
            push!(atoms, SmilesAtom(:H))
        end

        return atoms
    else
        a = atomsymbol!(state)
        if a === nothing
            return SmilesAtom[]
        else
            sym = a.second[1].second
            arom = a.second[2].second
            return [SmilesAtom(sym, 0, 1, nothing, arom, :unspecified)]
        end
    end
end

function atom!(state::SmartsParser)
    """ Atom <- '[' (AtomProp / LogicalOperator)+ ']' / AtomSymbol
    """
    c = read(state)
    if c == '['
        forward!(state)
        q = lglowand!(state, atomprop!)
        cq = read(state)
        @assert cq == ']' "(atomquery!) unexpected token: $(cq)"
        forward!(state)
        return [SmartsAtom(q)]
    else
        a = atomsymbol!(state)
        if a === nothing
            return SmartsAtom[]
        else
            return [SmartsAtom(a)]
        end
    end
end

function atomprop!(state::SmartsParserState)
    """ AtomProp <- '\$(' RecursiveQuery ')' / Mass / Symbol / AtomNum /
        Stereo / CHG / [DHRrvX]
    """
    c = read(state)
    c2 = lookahead(state, 1)
    if haskey(ATOMSYMBOLMAP, string(uppercase(c), c2))
        # Two letter atoms
        forward!(state, 2)
        if string(uppercase(c), c2) in ("As", "Se")
            return :and => (
                :atomsymbol => Symbol(uppercase(c), c2),
                :isaromatic => islowercase(c)
            )
        end
        @assert isuppercase(c)
        return :atomsymbol => Symbol(c, c2)
    elseif haskey(SMARTS_ATOM_COND_SYMBOL, c)
        # Neighbor and ring conditions
        forward!(state)
        if isdigit(c2)
            num = parse(Int, c2)
            forward!(state)
        else
            c in (:r, :R) && return :not => (:atom_sssrcount => 0)
            num = 1
        end
        return SMARTS_ATOM_COND_SYMBOL[c] => num
    elseif c in ('A', 'a', '*')
        forward!(state)
        c == 'A' && return :isaromatic => false
        c == 'a' && return :isaromatic => true
        c == '*' && return :any => true
    elseif haskey(ATOMSYMBOLMAP, string(uppercase(c)))
        # Single letter atoms
        forward!(state)
        if uppercase(c) in "BCNOPS"
            return :and => (
                :atomsymbol => Symbol(uppercase(c)),
                :isaromatic => islowercase(c)
            )
        end
        @assert isuppercase(c)
        return :atomsymbol => Symbol(c)
    elseif c == '#'
        # Atomic number
        forward!(state)
        start = state.pos
        while isdigit(lookahead(state, 1))
            forward!(state)
        end
        num = parse(Int, SubString(state.input, start, state.pos))
        forward!(state)
        return :atomsymbol => atomsymbol(num)
    elseif c in keys(SMARTS_CHARGE_SIGN)
        # Charge
        forward!(state)
        c2 = read(state)
        if isdigit(c2)
            chg = parse(Int, c2)
            forward!(state)
        else
            chg = 1
            while read(state) == c
                forward!(state)
                chg += 1
            end
        end
        return :charge => chg * SMARTS_CHARGE_SIGN[c]
    elseif c == '@'
        # Stereo
        # @ -> anticlockwise, @@ -> clockwise, ? -> or not specified
        forward!(state)
        cw = false
        if read(state) == '@'
            forward!(state)
            cw = true
        end
        if read(state) == '?'
            forward!(state)
            return :not => (:stereo => cw ? :anticlockwise : :clockwise)
        else
            return :stereo => cw ? :clockwise : :anticlockwise
        end
    elseif isdigit(c)
        # Isotope
        start = state.pos
        while isdigit(lookahead(state, 1))
            forward!(state)
        end
        num = SubString(state.input, start, state.pos)
        forward!(state)
        return :mass => parse(Int, num)
    elseif c == '$' && lookahead(state, 1) == '('
        # Recursive
        forward!(state, 2)
        start = state.pos
        while read(state) != ')'
            forward!(state)
        end
        q = SubString(state.input, start, state.pos - 1)
        forward!(state)
        return :recursive => q
    end
end


text:
using PyCall
using SparseArrays
using Images
using Colors
import Base: length

struct FluxDiagram{T<:Real}
    ts::Array{T,1}
    outputdirectory::String
end
export FluxDiagram

length(p::FluxDiagram) = 1
export length

Broadcast.broadcastable(p::FluxDiagram) = Ref(p)
export broadcastable

"""
return an array of all png images associated with a FluxDiagram object
"""
getimages(fd::FluxDiagram) = getdiagram.(fd,1:length(fd.ts))
export getimages

"""
return the png image associated with the index ind
"""
getdiagram(fd::FluxDiagram,ind::Int64) = load(string(joinpath(fd.outputdirectory,"flux_diagram_"),ind,".png"))
export getdiagram

"""
generate a png representing spc at location path
"""
function drawspc(spc::Species,path::String=".")
    name = spc.name
    fname = string(name,".png")
    if !(path in readdir("."))
        mkdir(path)
    else
        if fname in readdir(path)
            return
        end
    end
    if spc.inchi != ""
        mol = molecule.Molecule().fromInChI(spc.inchi)
    elseif spc.smiles != ""
        mol = molecule.Molecule().fromSMILES(spc.smiles)
    else
        throw(error("no smiles or inchi for molecule $name"))
    end
    mol.draw(joinpath(path,fname))
end
export drawspc

"""
generate pngs for all species in phase and store them in the "species" folder
"""
function drawspecies(phase::T) where {T<:AbstractPhase}
    for spc in phase.species
        drawspc(spc,"species")
    end
end
export drawspecies

"""
generates and returns the image of a single flux diagram at the given time point
all PyPlot colorscheme names are valid inputs for colorscheme
"""
function getfluxdiagram(bsol,t;centralspecieslist=Array{String,1}(),superimpose=false,
    maximumnodecount=50, maximumedgecount=50, concentrationtol=1e-6, speciesratetolerance=1e-6,
    maximumnodepenwidth=10.0,maximumedgepenwidth=10.0,radius=1,centralreactioncount=-1,outputdirectory="fluxdiagrams",
    colorscheme="viridis")

    fd = makefluxdiagrams(bsol,[t]; centralspecieslist=centralspecieslist,superimpose=superimpose,
        maximumnodecount=maximumnodecount, maximumedgecount=maximumedgecount, concentrationtol=concentrationtol,
        speciesratetolerance=speciesratetolerance,maximumnodepenwidth=maximumnodepenwidth,
        maximumedgepenwidth=maximumedgepenwidth,radius=radius,centralreactioncount=centralreactioncount,
        outputdirectory=outputdirectory,colorscheme=colorscheme)

    return getdiagram(fd,1)
end
export getfluxdiagram

"""
generates a series of flux diagrams at the time points indicated
each flux diagram will have the same nodes and edges as determined by the options
all PyPlot colorscheme names are valid inputs for colorscheme
"""
function makefluxdiagrams(bsol,ts;centralspecieslist=Array{String,1}(),superimpose=false,
    maximumnodecount=50, maximumedgecount=50, concentrationtol=1e-6, speciesratetolerance=1e-6,
    maximumnodepenwidth=10.0,maximumedgepenwidth=10.0,radius=1,centralreactioncount=-1,outputdirectory="fluxdiagrams",
    colorscheme="viridis")

    specieslist = bsol.domain.phase.species
    speciesnamelist = getfield.(specieslist,:name)
    numspecies = length(specieslist)
    reactionlist = bsol.domain.phase.reactions
    phase = bsol.domain.phase
    if !isdir(outputdirectory)
        mkdir(outputdirectory)
    end

    concentrations = hcat([bsol.sol(t)[1:numspecies]./getV(bsol,t) for t in ts]...)

    reactionrates = zeros(length(reactionlist),length(ts))
    for (i,t) in enumerate(ts)
        cs,kfs,krevs = calcthermo(bsol.domain,bsol.sol(t),t)[[2,9,10]]
        reactionrates[:,i] = [getrate(rxn,cs,kfs,krevs) for rxn in reactionlist]
    end

    drawspecies(bsol.domain.phase)
    speciesdirectory = joinpath(pwd(),"species")

    #find central species
    centralspeciesindices = Array{Int64,1}()
    if length(centralspecieslist) != 0
        for (j,centralspecies) in enumerate(centralspecieslist)
            for (i,species) in enumerate(specieslist)
                if species.name == centralspecies
                    push!(centralspeciesindices,i)
                    break
                end
            end
            if length(centralspeciesindices) != j
                throw(error("Central species $centralspecies could not be found"))
            end
        end
    end

    speciesrates = zeros(numspecies,numspecies,length(ts))
    for (index,reaction) in enumerate(reactionlist)
        rate = reactionrates[index]
        if length(reaction.pairs[1]) > 1
            pairs = reaction.pairs
        else
            pairs = getpairs(reaction)
        end
        for (reactant,product) in pairs
            reactantindex = findfirst(y->y==reactant,speciesnamelist)
            productindex = findfirst(y->y==product,speciesnamelist)
            speciesrates[reactantindex,productindex,:] .+= rate
            speciesrates[productindex,reactantindex,:] .-= rate
        end
    end

    maxconcentrations = maximum(concentrations,dims=2)
    maxconcentration = maximum(maxconcentrations)

    maxreactionrates = maximum(abs.(reactionrates),dims=2)

    maxspeciesrates = maximum(abs.(speciesrates),dims=3)
    maxspeciesrate = maximum(maxspeciesrates)

    speciesindex = sortperm(reshape(maxspeciesrates,numspecies^2),rev=true)

    nodes = Array{Int64,1}()
    edges = []
    if !superimpose && length(centralspecieslist) != 0
        for centralspeciesindex in centralspeciesindices
            push!(nodes,centralspeciesindex)
            addadjacentnodes!(centralspeciesindex,nodes,edges,phase,
                maxreactionrates,maxspeciesrates,centralreactioncount,radius,Array{Int64,1}(),speciesnamelist)
        end
    else
        for i = 1:numspecies^2
            productindex = div(speciesindex[i],numspecies)+1
            reactantindex = rem(speciesindex[i],numspecies)
            if reactantindex == 0
                reactantindex = numspecies
                productindex -= 1
            end
            if reactantindex > productindex
                continue
            end
            if maxspeciesrates[reactantindex,productindex] == 0
                break
            end
            if !(reactantindex in nodes) && length(nodes) < maximumnodecount
                push!(nodes,reactantindex)
            end
            if !(productindex in nodes) && length(nodes) < maximumnodecount
                push!(nodes,productindex)
            end
            if !((reactantindex,productindex) in edges) && !((productindex,reactantindex) in edges)
                push!(edges,(reactantindex,productindex))
            end
            if length(nodes) > maximumnodecount
                break
            end
            if length(edges) >= maximumedgecount
                break
            end
        end
        if superimpose && length(centralspecieslist) > 0
            nodescopy = nodes[:]
            for centralspeciesindex in centralspeciesindices
                if !(centralspeciesindex in nodes)
                    push!(nodes,centralspeciesindex)
                    addadjacentnodes!(centralspeciesindex,nodes,edges,phase,
                        maxreactionrates,maxspeciesrates,centralreactioncount,-1,nodescopy,speciesnamelist)
                end
            end
        end
    end

    graph = pydot.Dot("flux_diagram",graph_type="digraph",overlap="false")
    graph.set_rankdir("LR")
    graph.set_fontname("sans")
    graph.set_fontsize("10")

    for index in nodes
        species = specieslist[index]
        node = pydot.Node(name=species.name)
        node.set_penwidth(maximumnodepenwidth)
        graph.add_node(node)

        speciesindex = string(species.name,".png")
        imagepath = ""
        if !isdir(speciesdirectory)
            continue
        end
        for (root,dirs,files) in walkdir(speciesdirectory)
            for f in files
                if f == speciesindex
                    imagepath = joinpath(root,f)
                    break
                end
            end
        end
        if isfile(imagepath)
            node.set_image(imagepath)
            node.set_label(" ")
        end
    end

    for (reactantindex,productindex) in edges
        if reactantindex in nodes && productindex in nodes
            reactant = specieslist[reactantindex]
            product = specieslist[productindex]
            edge = pydot.Edge(reactant.name,product.name)
            edge.set_penwidth(maximumedgepenwidth)
            graph.add_edge(edge)
        end
    end

    graph = pydot.graph_from_dot_data(graph.create_dot(prog="dot"))[1]

    for t in 1:length(ts)
        slope = -maximumnodepenwidth / log10(concentrationtol)
        for index in nodes
            species = specieslist[index]
            if occursin(r"^[a-zA-Z0-9_]*$",species.name)
                species_string = species.name
            else
                species_string = string("\"",species.name,"\"")
            end

            node = graph.get_node(species_string)[1]
            concentration = concentrations[index,t] / maxconcentration
            if concentration < concentrationtol
                penwidth = 0.0
            else
                penwidth = round((slope*log10(concentration)+maximumnodepenwidth)*1.0e3)/1.0e3
            end

            node.set_penwidth(penwidth)
        end

        slope = -maximumedgepenwidth / log10(speciesratetolerance)
        minspeciesrate = Inf
        for index in 1:length(edges)
            reactantindex,productindex = edges[index]
            sprate = abs(speciesrates[reactantindex,productindex,t])
            if sprate > speciesratetolerance && minspeciesrate > sprate
                minspeciesrate = sprate
            end
        end

        for index in 1:length(edges)
            reactantindex,productindex = edges[index]
            if reactantindex in nodes && productindex in nodes
                reactant = specieslist[reactantindex]
                product = specieslist[productindex]

                if occursin(r"^[a-zA-Z0-9_]*$",reactant.name)
                    reactant_string = reactant.name
                else
                    reactant_string = string("\"",reactant.name,"\"")
                end

                if occursin(r"^[a-zA-Z0-9_]*$",product.name)
                    product_string = product.name
                else
                    product_string = string("\"",product.name,"\"")
                end

                edge = graph.get_edge(reactant_string,product_string)[1]

                speciesrate = speciesrates[reactantindex,productindex,t] / maxspeciesrate
                if speciesrate < 0
                    edge.set_dir("back")
                    speciesrate = -speciesrate
                else
                    edge.set_dir("forward")
                end

                if speciesrate < speciesratetolerance
                    penwidth = 0.0
                    edge.set_dir("none")
                else
                    penwidth = round((slope*log10(speciesrate) + maximumedgepenwidth)*1.0e3)/1.0e3
                end

                edge.set_penwidth(penwidth)
                edge.set_color(getcolor(speciesrates[reactantindex,productindex,t],maxspeciesrate,minspeciesrate,colorscheme))

            end
        end

        if ts[t] == 0.0
            label = "t = 0 s"
        else
            tval = log10(ts[t])
            label = "t = 10^$tval s"
        end

        graph.set_label(label)
        graph.write_dot(joinpath(outputdirectory,"flux_diagram_$t.dot"))
        graph.write_png(joinpath(outputdirectory,"flux_diagram_$t.png"))
        graph.write_svg(joinpath(outputdirectory,"flux_diagram_$t.svg"))
    end
    return FluxDiagram(ts,outputdirectory)
end

export makefluxdiagrams

function addadjacentnodes!(targetnodeindex,nodes,edges,phase,maxreactionrates,maxspeciesrates,reactioncount,rad,mainnodes,speciesnamelist)
    if rad == 0
        return
    elseif rad < 0 && targetnodeindex in mainnodes
        return
    else
        lt(x::Int64,y::Int64) = maxreactionrates[x] < maxreactionrates[y]
        targetreactionindices = Array{Int64,1}()
        for reaction in phase.reactions
            if targetnodeindex in reaction.reactantinds || targetnodeindex in reaction.productinds
                push!(targetreactionindices,reaction.index)
            end
        end
        sort!(targetreactionindices,lt=lt,rev=true)

        if reactioncount == -1
            targetreactionlist = phase.reactions[targetreactionindices]
        else
            targetreactionlist = phase.reactions[targetreactionindices[1:reactioncount]]
        end
        for reaction in targetreactionlist
            if length(reaction.pairs[1]) > 1
                pairs = reaction.pairs
            else
                pairs = getpairs(reaction)
            end
            for (reactant,product) in pairs
                rindex = findfirst(y->y==reactant,speciesnamelist)
                pindex = findfirst(y->y==product,speciesnamelist)
                if rindex == targetnodeindex
                    if !(pindex in nodes)
                        push!(nodes,pindex)
                        addadjacentnodes!(pindex,nodes,edges,phase,maxreactionrates,maxspeciesrates,reactioncount,rad-1,mainnodes,speciesnamelist)
                    end
                    if !((rindex,pindex) in edges) && !((pindex,rindex) in edges)
                        push!(edges,(rindex,pindex))
                    end
                end
                if pindex == targetnodeindex
                    if !(rindex in nodes)
                        push!(nodes,rindex)
                        addadjacentnodes!(rindex,nodes,edges,phase,maxreactionrates,maxspeciesrates,reactioncount,rad-1,mainnodes,speciesnamelist)
                    end
                    if !((rindex,pindex) in edges) && !((pindex,rindex) in edges)
                        push!(edges,(rindex,pindex))
                    end
                end
            end
        end
    end
end

function getcolor(speciesrate,maxspeciesrate,minspeciesrate,colorscheme="viridis")
    """
    gives the color corresponding to the scaled log species rate
    for a given PyPlot color scheme
    """
    scale = log(maxspeciesrate)-log(minspeciesrate)
    value = (log(abs(speciesrate))-log(minspeciesrate))/scale
    out = PyPlot.get_cmap(colorscheme)(value)[1:3]
    return "#"*hex(RGB(out...))
end


text:
module TestMatrixOfConstraints

using Test

import VecMathOptInterface
import SparseArrays

const MOI = VecMathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

MOI.Utilities.@product_of_sets(
    ScalarSets,
    MOI.EqualTo{T},
    MOI.GreaterThan{T},
    MOI.LessThan{T},
    MOI.Interval{T},
)

function _new_ScalarSets()
    return MOI.Utilities.GenericOptimizer{
        Float64,
        MOI.Utilities.ObjectiveContainer{Float64},
        MOI.Utilities.VariablesContainer{Float64},
        MOI.Utilities.MatrixOfConstraints{
            Float64,
            MOI.Utilities.MutableSparseMatrixCSC{
                Float64,
                Int,
                MOI.Utilities.ZeroBasedIndexing,
            },
            MOI.Utilities.Hyperrectangle{Float64},
            ScalarSets{Float64},
        },
    }()
end

MOI.Utilities.@product_of_sets(
    VectorSets,
    MOI.Nonnegatives,
    MOI.Nonpositives,
    MOI.SecondOrderCone,
)

function _new_VectorSets()
    return MOI.Utilities.GenericOptimizer{
        Int,
        MOI.Utilities.ObjectiveContainer{Int},
        MOI.Utilities.VariablesContainer{Int},
        MOI.Utilities.MatrixOfConstraints{
            Int,
            MOI.Utilities.MutableSparseMatrixCSC{
                Int,
                Int,
                MOI.Utilities.OneBasedIndexing,
            },
            Vector{Int},
            VectorSets{Int},
        },
    }()
end

function test_ScalarSets_basic()
    model = _new_ScalarSets()
    @test MOI.is_empty(model)
    src = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(
        src,
        """
variables: x, y
minobjective: x + y
x >= 1.0
y == 2.1
c: x + 2.0 * y <= 3.0
""",
    )
    index_map = MOI.copy_to(model, src)
    @test MOI.is_empty(model) == false
    for (k, v) in index_map
        @test MOI.is_valid(src, k)
        @test MOI.is_valid(model, v)
    end
    @test length(MOI.get(src, MOI.ListOfConstraintTypesPresent())) == 3
    MOI.empty!(model)
    @test MOI.is_empty(model)
    return
end

function test_VectorSets_basic()
    model = _new_VectorSets()
    @test MOI.is_empty(model)
    src = MOI.Utilities.Model{Int}()
    x = MOI.add_variables(src, 2)
    c = MOI.add_constraint(
        src,
        MOI.VectorAffineFunction{Int}(
            MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.(1, x)),
            [1, 3],
        ),
        MOI.SecondOrderCone(2),
    )
    index_map = MOI.copy_to(model, src)
    @test MOI.is_empty(model) == false
    for (k, v) in index_map
        @test MOI.is_valid(src, k)
        @test MOI.is_valid(model, v)
    end
    @test length(MOI.get(src, MOI.ListOfConstraintTypesPresent())) == 1
    MOI.empty!(model)
    @test MOI.is_empty(model)
    return
end

function test_ScalarSets_supports_constraint()
    model = _new_ScalarSets()
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Float64},
        MOI.LessThan{Float64},
    )
    @test !MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{Float64},
        MOI.Nonnegatives,
    )
    return
end

function test_VectorSets_supports_constraint()
    model = _new_VectorSets()
    @test !MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Float64},
        MOI.LessThan{Float64},
    )
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{Int},
        MOI.Nonnegatives,
    )
    return
end

function test_delete()
    model = _new_ScalarSets()
    x = MOI.add_variable(model)
    @test_throws MOI.DeleteNotAllowed MOI.delete(model, x)
    @test_throws MOI.DeleteNotAllowed MOI.delete(model, [x])
    return
end

function test_ScalarSets_get_ConstraintFunction()
    model = _new_ScalarSets()
    src = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(
        src,
        """
variables: x, y
minobjective: x + y
c: x + 2.0 * y <= 3.0
""",
    )
    index_map = MOI.copy_to(model, src)
    c = MOI.get(model, MOI.ConstraintIndex, "c")
    c_f = MOI.get(model, MOI.ConstraintFunction(), c)
    @test MOI.Utilities.map_indices(c_f) do x
        return index_map[x]
    end ≈ MOI.get(src, MOI.ConstraintFunction(), index_map[c])
    return
end

function test_VectorSets_get_ConstraintFunction()
    model = _new_VectorSets()
    src = MOI.Utilities.Model{Int}()
    x = MOI.add_variables(src, 2)
    c = MOI.add_constraint(
        src,
        MOI.VectorAffineFunction{Int}(
            MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.(1, x)),
            [1, 3],
        ),
        MOI.SecondOrderCone(2),
    )
    index_map = MOI.copy_to(model, src)
    c_f = MOI.get(model, MOI.ConstraintFunction(), c)
    @test MOI.Utilities.map_indices(c_f) do x
        return index_map[x]
    end ≈ MOI.get(src, MOI.ConstraintFunction(), index_map[c])
    return
end

function test_ScalarSets_get_ConstraintSet()
    model = _new_ScalarSets()
    src = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(
        src,
        """
variables: x, y
minobjective: x + y
c: 2.0 * x + y >= 4.0
""",
    )
    index_map = MOI.copy_to(model, src)
    c = MOI.get(model, MOI.ConstraintIndex, "c")
    c_set = MOI.get(model, MOI.ConstraintSet(), c)
    @test c_set == MOI.get(src, MOI.ConstraintSet(), index_map[c])
    return
end

function test_VectorSets_get_ConstraintSet()
    model = _new_VectorSets()
    src = MOI.Utilities.Model{Int}()
    x = MOI.add_variables(src, 2)
    c = MOI.add_constraint(
        src,
        MOI.VectorAffineFunction{Int}(
            MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.(1, x)),
            [1, 3],
        ),
        MOI.SecondOrderCone(2),
    )
    index_map = MOI.copy_to(model, src)
    c_set = MOI.get(model, MOI.ConstraintSet(), c)
    @test c_set == MOI.get(src, MOI.ConstraintSet(), index_map[c])
    return
end

function test_add_after_final_touch()
    model = _new_ScalarSets()
    @test MOI.is_empty(model)
    src = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(
        src,
        """
variables: x, y
minobjective: x + y
x >= 1.0
y == 2.1
c: x + 2.0 * y <= 3.0
""",
    )
    index_map = MOI.copy_to(model, src)
    x = MOI.get(src, MOI.VariableIndex, "x")
    @test_throws(
        MOI.AddConstraintNotAllowed,
        MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(
                [MOI.ScalarAffineTerm(1.0, index_map[x])],
                0.0,
            ),
            MOI.LessThan(2.0),
        ),
    )
    return
end

function test_UnsupportedConstraint()
    model = _new_ScalarSets()
    src = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(
        src,
        """
variables: x, y
c: [x, y] in Nonnegatives(2)
""",
    )
    @test_throws MOI.UnsupportedConstraint MOI.copy_to(model, src)
    return
end

function test_ScalarSets()
    optimizer = MOI.Utilities.GenericOptimizer{
        Float64,
        MOI.Utilities.ObjectiveContainer{Float64},
        MOI.Utilities.VariablesContainer{Float64},
        MOI.Utilities.MatrixOfConstraints{
            Float64,
            MOI.Utilities.MutableSparseMatrixCSC{
                Float64,
                Int,
                MOI.Utilities.ZeroBasedIndexing,
            },
            MOI.Utilities.Hyperrectangle{Float64},
            ScalarSets{Float64},
        },
    }()
    MOI.Test.runtests(
        optimizer,
        MOI.Test.Config(exclude = Any[MOI.optimize!]),
        include = ["test_model_", "test_objective_", "test_variable_"],
        exclude = [
            "test_model_ScalarFunctionConstantNotZero",
            "test_model_default_DualStatus",
            "test_model_default_PrimalStatus",
            "test_model_default_TerminationStatus",
        ],
    )
    return
end

function test_modify()
    model = MOI.Utilities.GenericOptimizer{
        Int,
        MOI.Utilities.ObjectiveContainer{Int},
        MOI.Utilities.VariablesContainer{Int},
        MOI.Utilities.MatrixOfConstraints{
            Int,
            MOI.Utilities.MutableSparseMatrixCSC{
                Int,
                Int,
                MOI.Utilities.OneBasedIndexing,
            },
            MOI.Utilities.Hyperrectangle{Int},
            ScalarSets{Int},
        },
    }()
    x = MOI.add_variable(model)
    func = 2x
    set = MOI.EqualTo(1)
    c = MOI.add_constraint(model, func, set)
    MOI.Utilities.final_touch(model, nothing)
    @test_throws MOI.DeleteNotAllowed(c) MOI.delete(model, c)
    err = MOI.AddConstraintNotAllowed{typeof(func),typeof(set)}(
        MOI.Utilities._MATRIXOFCONSTRAINTS_MODIFY_NOT_ALLOWED_ERROR_MESSAGE,
    )
    @test_throws err MOI.add_constraint(model, func, set)
    return
end

MOI.Utilities.@struct_of_constraints_by_set_types(
    ZerosOrNot,
    MOI.Zeros,
    Union{MOI.Nonnegatives,MOI.Nonpositives},
)

MOI.Utilities.@product_of_sets(Zeros, MOI.Zeros)

function test_multicone()
    T = Int
    Indexing = MOI.Utilities.OneBasedIndexing
    model = MOI.Utilities.GenericOptimizer{
        T,
        MOI.Utilities.ObjectiveContainer{T},
        MOI.Utilities.VariablesContainer{T},
        ZerosOrNot{T}{
            MOI.Utilities.MatrixOfConstraints{
                T,
                MOI.Utilities.MutableSparseMatrixCSC{T,Int,Indexing},
                Vector{T},
                Zeros{T},
            },
            MOI.Utilities.MatrixOfConstraints{
                T,
                MOI.Utilities.MutableSparseMatrixCSC{T,Int,Indexing},
                Vector{T},
                VectorSets{T},
            },
        },
    }()
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    MOI.add_constraint(
        model,
        MOI.Utilities.vectorize([T(5) * x + T(2)]),
        MOI.Zeros(1),
    )
    MOI.add_constraint(
        model,
        MOI.Utilities.vectorize([T(3) * y + T(1)]),
        MOI.Nonnegatives(1),
    )
    MOI.add_constraint(
        model,
        MOI.Utilities.vectorize([T(6), T(7) * x, T(4)]),
        MOI.Nonpositives(1),
    )
    MOI.Utilities.final_touch(model, nothing)
    @test convert(
        SparseArrays.SparseMatrixCSC{T,Int},
        model.constraints.moi_zeros.coefficients,
    ) == [T(5) zero(T)]
    @test model.constraints.moi_zeros.constants == T[2]
    @test convert(
        SparseArrays.SparseMatrixCSC{T,Int},
        model.constraints.moi_nonnegatives.coefficients,
    ) == [zero(T) T(3); zero(T) zero(T); T(7) zero(T); zero(T) zero(T)]
    @test model.constraints.moi_nonnegatives.constants == T[1, 6, 0, 4]
    return
end

MOI.Utilities.@product_of_sets(
    MatrixSets,
    MOI.PositiveSemidefiniteConeTriangle,
    MOI.PositiveSemidefiniteConeSquare,
    MOI.LogDetConeTriangle,
    MOI.LogDetConeSquare,
    MOI.RootDetConeTriangle,
    MOI.RootDetConeSquare,
)

function test_matrix_sets(::Type{T} = Int) where {T}
    model = MOI.Utilities.GenericOptimizer{
        T,
        MOI.Utilities.ObjectiveContainer{T},
        MOI.Utilities.VariablesContainer{T},
        MOI.Utilities.MatrixOfConstraints{
            T,
            MOI.Utilities.MutableSparseMatrixCSC{
                T,
                T,
                MOI.Utilities.OneBasedIndexing,
            },
            Vector{T},
            MatrixSets{T},
        },
    }()
    function test_set(set)
        MOI.empty!(model)
        func = MOI.Utilities.zero_with_output_dimension(
            MOI.VectorAffineFunction{T},
            MOI.dimension(set),
        )
        ci = MOI.add_constraint(model, func, set)
        MOI.Utilities.final_touch(model, nothing)
        @test MOI.get(model, MOI.ConstraintSet(), ci) == set
    end
    test_set(MOI.PositiveSemidefiniteConeTriangle(2))
    test_set(MOI.PositiveSemidefiniteConeSquare(2))
    test_set(MOI.LogDetConeTriangle(2))
    test_set(MOI.LogDetConeSquare(2))
    test_set(MOI.RootDetConeTriangle(2))
    test_set(MOI.RootDetConeSquare(2))
    return
end

end

TestMatrixOfConstraints.runtests()


text:
module ASTGen
export to_ast
using MLStyle
using JSON

@nospecialize
@inline function ID(x)
    x
end

count = 0

IS_BREAK = Symbol("____Py2Jl____UNSAFE_MANGLE_IS_BREAK")
MANGLE_BASE = Symbol("____Py2Jl____UNSAFE_MANGLE_TMP")

mangle() = begin
    global count
    let name = Symbol(MANGLE_BASE, string(count))
        count = count + 1
        name
    end
end

struct Record end

"""
define a pattern `Record(a, b, c=>c, ...)`
identical to Dict `Dict(a=a, b=b, c=>, ...)`
"""
function MLStyle.pattern_uncall(
    ::Type{Record},
    self,
    tparams::AbstractArray,
    targs::AbstractArray,
    args::AbstractArray,
)
    isempty(tparams) || error("A Record pattern requires no type params.")
    isempty(targs) || error("A Record pattern requires no type arguments.")
    isempty(tparams) || error("A Record pattern requires no type parameters.")
    args = ((arg isa Symbol ?
                let key = QuoteNode(arg)
                    :($key => $arg)
                end : arg)
            for arg in args)
    pat = :($Dict($(args...)))
    self(pat)
end

macro not_implemented_yet()
    :(throw("notimplemented yet"))
end

empty_block = Expr(:block)

function ret_nil(node)
    Expr(:block, node, nothing)
end

"""
Cannot generate python annotations for
the semantics is not the same as Julia's annoations.
"""
function annotate(sym, ty)
    @not_implemented_yet
end

function assign(target, value)
    Expr(:(=), target, value)
end

function (<|)(f, arg)
    f(arg)
end

function gather(args)
    isempty(args) ? nothing :
    length(args) === 1 && args[1] isa Expr ? args[1] :
    Expr(:block, args...)

end

function for_iter(f :: Function, iter_arg, seq, body)
    basic = Expr(:for, assign(iter_arg, seq), body)
    token = mangle()
    result = mangle()
    check_break = Expr(
        :block,
        assign(token, Expr(:call, Ref, false)),
        assign(IS_BREAK, token),
        basic,
        f(token),
    )
end

function for_iter(iter_arg, seq, body)
    for_iter((_) -> nothing, iter_arg, seq, body)
end

function while_loop(f :: Function, cond, body)
    basic = Expr(:while, cond, body)
    token = mangle()
    result = mangle()
    check_break = Expr(:block,
        assign(token, Expr(:call, Ref, false)),
        assign(IS_BREAK, token),
        basic,
        f(token),
    )
end

function while_loop(cond, body)
    while_loop((_) -> nothing, cond, body)
end

function call(fn, args...)
    Expr(:call, fn, args...)
end

function as_global(names...)
    Expr(:global, names...)
end

function get_attr(expr, attr :: Symbol)
    Expr(:., expr, QuoteNode(attr))
end

function break!()
    Expr(:block, assign(get_attr(IS_BREAK, :x), true), Expr(:break))
end

function continue!()
    Expr(:continue)
end

function ifelse(cond, then, else_ :: Nothing)
    Expr(:if, cond, then)
end

function ifelse(cond, then, else_)
    Expr(:if, cond, then, else_)
end

function ifelse(cond, then)
    Expr(:if, cond, then)
end

function isinstance(inst, typs :: Union{Tuple, Vector})
    foldr(typs, init=true) do typ, last
        Expr(:||, isinstance(inst, typ), last)
    end
end

function isinstance(inst, typ)
    Expr(:call, isa, inst, typ)
end

struct Context
    filename
    imports
end

filename(x) = x.filename

function tag_loc(ctx, x)
    _filename = filename(ctx)
    @λ begin
        (if _filename === nothing end &&
         Record(lineno, colno)) -> LineNumberNode(lineno)

        Record(lineno, colno) -> LineNumberNode(lineno, _filename)

        _ -> nothing
    end
end

function trans_block(ctx, seq)
    res = []
    for each in seq
        loc = tag_loc(ctx, each)
        if loc !== nothing
            push!(res, loc)
        end
        push!(res, apply(ctx, each))
    end
    res
end

function fundef(ctx, class, body, kwarg, args, kw_defaults, kwonlyargs, defaults, vararg, fn_ast)
    if kwarg === nothing && isempty(kwonlyargs)
        if isempty(kw_defaults) && isempty(defaults)
            f = @λ Expr(:no_eval, arg, annotation) -> (arg, annotation)
            arg_anno = map(x -> f(apply(ctx, x)), args)
            args = map(first, arg_anno)
            annos = [annotate(arg, anno) for (arg, anno) in arg_anno if anno !== nothing]

            if class == "Lambda"
                Expr(:function, Expr(:tuple, args...), Expr(:block, annos..., apply(ctx, body)))
            else

                body = gather <| trans_block(ctx, body)
                decorator_list = fn_ast[:decorator_list]
                fn_name = Symbol(fn_ast[:name])
                init = Expr(:function, Expr(:call, fn_name, args...), Expr(:block, annos..., body))
                reduce(decorator_list, init=init) do last, decorator
                    decorator = apply(ctx, decorator)
                    wrapped = call(decorator, last)
                    Expr(:(=), fn_name, Expr(:let, empty_block, wrapped))
                end |> ret_nil
            end
        else
            @not_implemented_yet
        end
    else
        @not_implemented_yet
    end
end

function apply(ctx, x)
    _apply = @λ begin
        (num :: Number)  -> num
        (str :: String)  -> str
        (:: Nothing)     -> nothing

        Record(:class => "Module", body) ->
        let body = gather <| trans_block(ctx, body)

            filename === nothing ?
                body                 :
                Expr(:module, true, Symbol(filename), body)
        end

        Record(:class => "Name", id) ->
        Symbol(id)
        Record(:class => "Num", n) -> n

        Record(:class => "List", elts) ->
        Expr(:vect, map(apply, elts)...)

        Record(:class => "Tuple", elts) ->
        Expr(:tuple, map(apply, elts)..., )

        Record(:class => "Return", value) -> Expr(:return, apply(ctx, value))

        Record(:class => "arg", annotation, arg) -> Expr(:no_eval, Symbol(arg), apply(ctx, annotation))

        # FunctionDef
        (Record(class,
                body,
                :args => Record(
                    kwarg,
                    args,
                    kw_defaults,
                    kwonlyargs,
                    defaults,
                    vararg
                )) && fn_ast) -> fundef(ctx, class, body, kwarg, args, kw_defaults, kwonlyargs, defaults, vararg, fn_ast)

        Record(:class   => "Assign", targets, value) ->
        (reduce(targets, init = apply(ctx, value)) do last, target
         Expr(:(=), apply(ctx, target), last)
         end |> ret_nil)

        Record(:class => "AugAssign", target, op, value) -> @not_implemented_yet

        Record(:class => "AnnAssign", target, annotation, value) ->
        (annotate(apply(ctx, target), apply(ctx, annotation)) |>
         target -> assign(target, value)            |>
         ret_nil)

        Record(:class => "For", target, iter, body, orelse) ->
        let target = apply(ctx, target),
            iter = apply(ctx, iter),
            body = gather <| trans_block(ctx, body),
            orelse = gather <| trans_block(ctx, orelse)

            for_iter(target, iter, body) do token
                ifelse(Expr(:call, !, get_attr(token, :x)), orelse)
            end

        end

        Record(:class => "While", test, body, orelse) ->
        let cond = apply(ctx, test),
            body = gather <| trans_block(ctx, body),
            orelse = gather <| trans_block(ctx, orelse)

            while_loop(cond, body) do token
                ifelse(Expr(:call, !, get_attr(token, :x)), orelse)
            end
        end

        Record(:class => "With") -> @not_implemented_yet

        Record(:class => "ClassDef") -> @not_implemented_yet

        Record(:class => "Try", body, handlers, orelse, finalbody) ->
        if !isempty(finalbody) || !isempty(orelse)
            @not_implemented_yet
        else
            ret = Expr(:try, gather <| trans_block(ctx, body))
            if isempty(handlers)
                return ret
            end
            except = mangle()
            args = ret.args
            push!(args, except)
            init = call(throw, except)
            foldr(handlers, init = init) do handler, last
                @match handler begin
                    Record(:type => exc, name, body) =>
                        let exc = apply(ctx, exc),

                            body = gather <| trans_block(ctx, body),

                            tc = isinstance(except, exc),

                            case = name === nothing ? body : gather([
                                assign(name, except),
                                body
                            ])

                            ifelse(tc, body, last)
                        end
                    _ => @error "Unknown python ast."
                end
            end |> it -> push!(args, it)
            ret
        end

        # runtime error
        Record(:class => "Raise", :exc => nothing, :cause => nothing)   -> @not_implemented_yet

        Record(:class => "Raise", exc, :cause => nothing) -> call(throw, apply(ctx, exc))

        Record(:class => "Raise", :exc => _, :cause => _) -> @not_implemented_yet

        Record(:class => "Import") -> @not_implemented_yet

        Record(:class => "ImportFrom") -> @not_implemented_yet

        Record(:class => "Global", names) -> Expr(:global, map(Symbol, names))

        Record(:class => "Pass") -> nothing

        Record(:class => "Break") -> break!()

        Record(:class => "Continue") -> continue!()

        Record(:class => "If", test, body, orelse) ->
        let cond = apply(ctx, test),
            body = gather <| trans_block(ctx, body),
            orelse = gather <| trans_block(ctx, orelse)

            ret_nil <| ifelse(cond, body, orelse)
        end

        Record(:class => "Expr", value) -> ret_nil <| apply(ctx, value)

        Record(:class => "Starred", value) -> Expr(:..., apply(ctx, value))

        Record(:class => "Call", func, args, keywords) ->
        begin
            func = apply(ctx, func)
            args = map(apply, args)
            keywords = [(it[:arg], apply(it[:value])) for it in keywords]
            kw_unpack = [Expr(:..., snd) for (fst, snd) in keywords if fst === nothing]
            kw_args = [Expr(:kw, fst, snd) for (fst, snd) in keywords if fst !== nothing]
            Expr(:call, func, Expr(:parameters, kw_unpack...), args..., kw_args...)
        end
        Record(:class => "BinOp", op, left, right) ->
        let op =  @match op[:class] begin
            # TODO, binary operator in Python cannot be mapped to Julia directly.
            # We should implement Python.(+), Python.(-)...
            "Add"     => (+)
            "Sub"     => (-)
            "Mult"    => (*)
            "Div"     => (/)
            "MatMult" => @not_implemented_yet
            "Mod"     => (%)
            "Pow"     => (^)
            "LShift"  => (<<)
            "RShift"  => (>>)
            "BitOr"   => (|)
            "BitXor"  => xor
            "BitAnd"  => (&)
            "FloorDiv"=> floor ∘ (/)
        end
            call(op, apply(ctx, left), apply(ctx, right))
        end

        # for Python 3.8+
        (Record(:class => "Constant", value)
         || # for Python 3.7- and 3.7
         Record(:class => "NameConstant", value)
         ) -> value

        Record(:class=> "Compare", left, ops, comparators) ->
        foldl(zip(ops, comparators) |> collect, init=apply(ctx, left)) do last, (op, comparator)
            let comparator = apply(ctx, comparator)
                f = @match op[:class] begin
                    #  Eq | NotEq | Lt | LtE | Gt | GtE | Is | IsNot | In | NotIn
                    "Eq" => (==)
                    "NotEq" => (!=)
                    "Lt" => (<)
                    "LtE" => (<=)
                    "Gt"  => (>)
                    "GtE" => (>=)
                    "Is"  => (===)
                    "IsNot" => (!==)
                    "In"    => (in)
                    "NotIn" => (!in)
                end
                :($f($last, $comparator))
            end
        end
        this ->
        let msg = "class: $(this[:class]), attributes: $(keys(this))."
            @match this begin
                Dict(:class=>"Module") => println(:aaa)
                ::T where T => println(T)
            end
            throw(msg)
        end
    end
    _apply(x)
end

"""
Check https://github.com/python/cpython/blob/master/Parser/Python.asdl
for more implementation details.
"""
function to_ast(ctx, python :: Dict)
    apply(ctx, python)
end

@specialize

end