text:
using ConservationLawsParticles
using SpecialFunctions, RecursiveArrayTools, DifferentialEquations, Plots

# model
V(t, x) = 0.0 # 2 + 0.5sin(2x) + 0.3sin(2sqrt(2)*x) + 0.2cos(2sqrt(7)*x)
W(t, x) = 0.0 # x^2/2 - abs(x)
W′(t, x) = 0.0 # x - sign(x)
mob(ρ) = 0.0 # max((ρ-1)^2 * (2ρ+1), 0)
smodel = DiffusiveSampledModel((V,), ((W′,),), (mob,), (Diffusion(1),))
imodel = DiffusiveIntegratedModel((V,), ((W,),), (mob,), (Diffusion(1),))

# initial condition
n = 20
x0 = ArrayPartition(vcat(
    range(-3, -2; length=2n),
    range(-1, 0; length=n),
    range(1, 2; length=3n)
    ))

boxsol(t, x; a, b) = (erf((x-a)/2sqrt(t)) - erf((x-b)/2sqrt(t))) / (2*(b-a))
refsol(t, x) = (
        boxsol(t, x; a=-3, b=-2) * (2n-1) +
        boxsol(t, x; a=-1, b=0) * (n-1) +
        boxsol(t, x; a=1, b=2) * (3n-1) +
        boxsol(t, x; a=-2, b=-1) + boxsol(t, x; a=0, b=1)
    ) / (6n-1)

# time span
tspan = (0., 1.0)

# ODE system for the particles
sprob = ODEProblem(velocities_diff!, x0, tspan, smodel)
# iprob = ODEProblem(velocities_diff!, x0, tspan, imodel)

abstol = reltol = 2e-8

@time ssol = solve(sprob, BS5(); abstol=abstol, reltol=reltol);

# plot the particle trajectories
plot(legend=false, title="Trajectories", xlabel="time", ylabel="position")
plot!(ssol, vars=eachindex(x0)[1:5:end]; color=:blue)
savefig("diffusion3.png")

# plot an animation of the density
anim = @animate for t in range(tspan...; length=101)
    plot(title="Density", xlabel="position", ylabel="density",
        xlims=(-7,7), ylims=(0,1/sqrt(pi)))
    plot!(x -> refsol(t, x); color=:blue, label="analytic")
    densityplot!(ssol(t); color=:red, label="particles")
end
gif(anim, "diffusion3.gif")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct HttpMessage <: SwaggerModel
    content::Any # spec type: Union{ Nothing, Any } # spec name: content

    function HttpMessage(;content=nothing)
        o = new()
        validate_property(HttpMessage, Symbol("content"), content)
        setfield!(o, Symbol("content"), content)
        o
    end
end # type HttpMessage

const _property_map_HttpMessage = Dict{Symbol,Symbol}(Symbol("content")=>Symbol("content"))
const _property_types_HttpMessage = Dict{Symbol,String}(Symbol("content")=>"Any")
Base.propertynames(::Type{ HttpMessage }) = collect(keys(_property_map_HttpMessage))
Swagger.property_type(::Type{ HttpMessage }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_HttpMessage[name]))}
Swagger.field_name(::Type{ HttpMessage }, property_name::Symbol) =  _property_map_HttpMessage[property_name]

function check_required(o::HttpMessage)
    true
end

function validate_property(::Type{ HttpMessage }, name::Symbol, val)
end

ext: jl


text:
#createmarketsapi.jl
# Use the newapp method to create the markets template
using Genie
Genie.newapp_webservice("markets")


text:
####################
# Helper Constants #
####################

# Intermidiate calculations calculations
const ECC2 = WGS84_f * (2.0 - WGS84_f) # Square of eccentricisty

##############
# Geocentric #
##############

export sGEOCtoECEF
"""
Convert geocentric position to equivalent Earth-fixed position.

Arguments:
- `geoc::Array{<:Real, 1}`: Geocentric coordinates (lon, lat, altitude) [rad] / [deg]
- `use_degrees:Bool`: If `true` interpret input as being in degrees.

Returns:
- `ecef::Array{<:Real, 1}`: Earth-fixed coordinates [m]
"""
function sGEOCtoECEF(geoc::Array{<:Real, 1} ; use_degrees::Bool=false)
    # Extract lat and lon
    lon = geoc[1]
    lat = geoc[2]

    # Handle non-explict use-degrees
    if length(geoc) == 3
        alt = geoc[3]
    else
        alt = 0.0
    end

    # Convert input to radians
    if use_degrees
        lat = lat*pi/180.0
        lon = lon*pi/180.0
    end

    # Check validity of input
    if lat < -pi/2 || lat > pi/2
        throw(ArgumentError("Lattiude, $lat, out of range. Must be between -90 and 90 degrees."))
    end

    # Compute Earth fixed coordinates
    r       = WGS84_a + alt
    x = r*cos(lat)*cos(lon)
    y = r*cos(lat)*sin(lon)
    z = r*sin(lat)

    return [x, y, z]
end

export sECEFtoGEOC
"""
Convert Earth-fixed position to geocentric location.

Arguments:
- `ecef::Array{<:Real, 1}`: Earth-fixed coordinated [m]
- `use_degrees:Bool`: If `true` returns result in units of degrees

Returns:
- `geoc`: Geocentric coordinates (lon, lat, altitude) [rad] / [deg]
"""
function sECEFtoGEOC(ecef::Array{<:Real, 1} ; use_degrees::Bool=false)
    # Expand ECEF coordinates
    x, y, z = ecef

    # Compute geocentric coordinates
    lat = atan(z, sqrt(x*x + y*y))
    lon = atan(y, x)
    alt = sqrt(x*x + y*y + z*z) - WGS84_a

    # Convert output to degrees
    if use_degrees
        lat = lat*180.0/pi
        lon = lon*180.0/pi
    end

    return [lon, lat, alt]
end

########################
# Geodetic Convertions #
########################

export sGEODtoECEF
"""
Convert geodetic position to equivalent Earth-fixed position.

Arguments:
- `geod::Array{<:Real, 1}`: Geodetic coordinates (lon, lat, altitude) [rad] / [deg]
- `use_degrees:Bool`: If `true` interpret input as being in degrees.

Returns:
- `ecef::Array{<:Real, 1}`: Earth-fixed coordinates [m]
"""
function sGEODtoECEF(geod::Array{<:Real, 1} ; use_degrees::Bool=false)
    # Extract lat and lon
    lon = geod[1]
    lat = geod[2]

    # Handle non-explict use-degrees
    if length(geod) == 3
        alt = geod[3]
    else
        alt = 0.0
    end

    # Convert input to radians
    if use_degrees
        lat = lat*pi/180.0
        lon = lon*pi/180.0
    end

    # Check validity of input
    if lat < -pi/2 || lat > pi/2
        throw(ArgumentError("Lattiude, $lat, out of range. Must be between -90 and 90 degrees."))
    end

    # Compute Earth-fixed position vector
    N = WGS84_a / sqrt(1.0 - ECC2*sin(lat)^2)
    x =           (N+alt)*cos(lat)*cos(lon)
    y =           (N+alt)*cos(lat)*sin(lon)
    z =  ((1.0-ECC2)*N+alt)*sin(lat)

    return [x, y, z]
end

export sECEFtoGEOD
"""
Convert geodetic coordinaties to Earth-fixed position

Arguments:
- `ecef::Array{<:Real, 1}`: Earth-fixed position [m]
- `use_degrees:Bool`: If `true` returns result in units of degrees

Returns:
- `geod::Array{<:Real, 1}`: Geocentric coordinates (lon, lat, altitude) [rad] / [deg]
"""
function sECEFtoGEOD(ecef::Array{<:Real, 1} ; use_degrees::Bool=false)
    # Expand ECEF coordinates
    x, y, z = ecef

    # Compute intermediate quantities
    epsilon  = eps(Float64) * 1.0e3 * WGS84_a # Convergence requirement as function of machine precision
    rho2 = x^2 + y^2                      # Square of the distance from the z-axis
    dz   = ECC2 * z
    N    = 0.0

    # Iteratively compute refine coordinates
    while true
        zdz    = z + dz
        Nh     = sqrt(rho2 + zdz^2)
        sinphi = zdz / Nh
        N      = WGS84_a / sqrt(1.0 - ECC2 * sinphi^2)
        dz_new = N * ECC2 * sinphi

        # Check convergence requirement
        if abs(dz - dz_new) < epsilon
            break
        end

        dz = dz_new
    end

    # Extract geodetic coordinates
    zdz = z + dz
    lat = atan(zdz, sqrt(rho2))
    lon = atan(y, x)
    alt = sqrt(rho2 + zdz^2) - N

    # Convert output to degrees
    if use_degrees
        lat = lat*180.0/pi
        lon = lon*180.0/pi
    end

    return [lon, lat, alt]
end

#######
# ENZ #
#######

export rECEFtoENZ
"""
Compute the rotation matrix from the Earth-fixed to the East-North-Up
coorindate basis.

Arguments:
- `station_ecef::Array{<:Real, 1}`: Earth-fixed cartesian station coordinates
- `conversion::Bool`: Conversion type to use. Can be "geocentric" or "geodetic"

Returns:
- `E::Array{Real, 2}`: Topocentric rotation matrix
"""
function rECEFtoENZ(ecef::Array{<:Real, 1} ; conversion::String="geodetic")
    if length(ecef) < 3
        throw(ArgumentError("Input coordinates must be length 3."))
    end

    # Compute Station Lat-Lon-Altitude
    if conversion == "geodetic"
        lon, lat, = sECEFtoGEOD(ecef, use_degrees=false)
    elseif conversion == "geocentric"
        lon, lat, = sECEFtoGEOC(ecef, use_degrees=false)
    else
        throw(ArgumentError("Unknown conversion method: $conversion"))
    end

    # Compute ENZ basis vectors
    eE = [-sin(lon) ; cos(lon) ; 0]
    eN = [-sin(lat)*cos(lon) ; -sin(lat)*sin(lon) ; cos(lat)]
    eZ = [cos(lat)*cos(lon) ; cos(lat)*sin(lon) ; sin(lat)]

    # Construct Rotation matrix
    E = hcat(eE, eN, eZ)'

    # Return Result
    return E
end

export rENZtoECEF
"""
Compute the rotation matrix from the Earth-fixed to the South-East-Zenith
coorindate basis.

Arguments:
- `station_ecef::Array{<:Real, 1}`: Earth-fixed cartesian station coordinates
- `conversion::Bool`: Conversion type to use. Can be "geocentric" or "geodetic"

Returns:
- `E::Array{Float64, 2}`: Topocentric rotation matrix
"""
function rENZtoECEF(ecef::Array{<:Real, 1} ; conversion::String="geodetic")
    # Check input coordinates
    if length(ecef) < 3
        throw(ArgumentError("Input coordinates must be length 3."))
    end

    return rECEFtoENZ(ecef, conversion=conversion)'
end

export sECEFtoENZ
"""
Compute the coordinates of an object in the topocentric frame of an
Earth-fixed frame

Arguments:
- `station_ecef::Array{<:Real, 1}`: Earth-fixed cartesian station coordinates
- `ecef::Array{<:Real, 1}`: Coordinates of the object in Earth-fixed station
- `conversion::Bool`: Conversion type to use. Can be "geocentric" or "geodetic"

Returns:
- `E::Array{Float64, 2}`: Topocentric rotation matrix
"""
function sECEFtoENZ(station_ecef::Array{<:Real, 1}, ecef::Array{<:Real, 1} ; conversion::String="geodetic")
    # Check input sizes
    if length(ecef) < 3
        throw(ArgumentError("Input ecef state must be at least length 3."))
    end

    if length(station_ecef) < 3
        throw(ArgumentError("Input station coordinates must be length 3."))
    end

    # Compute ENZ Rotation matrix
    E = rECEFtoENZ(station_ecef, conversion=conversion)

    # Transform range
    range_ecef = ecef[1:3] - station_ecef
    range_enz  = E * range_ecef

    # Transform range-rate (if necessary)
    if length(ecef) == 6
        range_rate_ecef = ecef[4:6]
        range_rate_enz  = E * range_rate_ecef
    end

    # Return
    if length(ecef) == 6
        sat_enz = vcat(range_enz, range_rate_enz)
    else
        sat_enz = range_enz
    end

    return sat_enz
end

export sENZtoECEF
"""
Compute the coordinates of an object in the topocentric frame of an
Earth-fixed frame

Arguments:
- `station_ecef::Array{<:Real, 1}`: Earth-fixed cartesian station coordinates
- `sez::Array{<:Real, 1}`: SEZ coordinates of the object
- `conversion::Bool`: Conversion type to use. Can be "geocentric" or "geodetic"

Returns:
- `E::Array{Float64, 2}`: Topocentric rotation matrix
"""
function sENZtoECEF(station_ecef::Array{<:Real, 1}, enz::Array{<:Real, 1} ; conversion::String="geodetic")
    # Check input sizes
    if length(enz) < 3
        throw(ArgumentError("Input ENZ state must be at least length 3."))
    end

    if length(station_ecef) < 3
        throw(ArgumentError("Input station coordinates must be length 3."))
    end

    # Compute ENZ Rotation matrix
    E = rENZtoECEF(station_ecef, conversion=conversion)

    # Transform range
    range_enz  = enz[1:3]
    range_ecef = E * range_enz

    # Transform range-rate (if necessary)
    if length(enz) == 6
        range_rate_enz  = enz[4:6]
        range_rate_ecef = E * range_rate_enz
    end

    # Return
    if length(enz) == 6
        sat_ecef = vcat(range_ecef + station_ecef, range_rate_ecef)
    else
        sat_ecef = range_ecef + station_ecef
    end

    return sat_ecef
end

#######
# SEZ #
#######

export rECEFtoSEZ
"""
Compute the rotation matrix from the Earth-fixed to the South-East-Zenith
coorindate basis.

Arguments:
- `station_ecef::Array{<:Real, 1}`: Earth-fixed cartesian station coordinates
- `conversion::Bool`: Conversion type to use. Can be "geocentric" or "geodetic"

Returns:
- `E::Array{Float64, 2}`: Topocentric rotation matrix
"""
function rECEFtoSEZ(ecef::Array{<:Real, 1} ; conversion::String="geodetic")
    if length(ecef) < 3
        throw(ArgumentError("Input coordinates must be length 3."))
    end

    # Compute Station Lat-Lon-Altitude
    if conversion == "geodetic"
        lon, lat, = sECEFtoGEOD(ecef, use_degrees=false)
    elseif conversion == "geocentric"
        lon, lat, = sECEFtoGEOC(ecef, use_degrees=false)
    else
        throw(ArgumentError("Unknown conversion method: $conversion"))
    end

    # Compute SEZ basis vectors
    eS = [sin(lat)*cos(lon) ; sin(lat)*sin(lon) ; -cos(lat)]
    eE = [-sin(lon) ; cos(lon) ; 0]
    eZ = [cos(lat)*cos(lon) ; cos(lat)*sin(lon) ; sin(lat)]

    # Construct Rotation matrix
    E = hcat(eS, eE, eZ)'

    # Return Result
    return E
end

export rSEZtoECEF
"""
Compute the rotation matrix from the Earth-fixed to the South-East-Zenith
coorindate basis.

Arguments:
- `station_ecef::Array{<:Real, 1}`: Earth-fixed cartesian station coordinates
- `conversion::Bool`: Conversion type to use. Can be "geocentric" or "geodetic"

Returns:
- `E::Array{Float64, 2}`: Topocentric rotation matrix
"""
function rSEZtoECEF(ecef::Array{<:Real, 1} ; conversion::String="geodetic")
    # Check input coordinates
    if length(ecef) < 3
        throw(ArgumentError("Input coordinates must be length 3."))
    end

    return rECEFtoSEZ(ecef, conversion=conversion)'
end

export sECEFtoSEZ
"""
Compute the coordinates of an object in the topocentric frame of an
Earth-fixed frame

Arguments:
- `station_ecef::Array{<:Real, 1}`: Earth-fixed cartesian station coordinates
- `ecef::Array{<:Real, 1}`: Coordinates of the object in Earth-fixed station
- `conversion::Bool`: Conversion type to use. Can be "geocentric" or "geodetic"

Returns:
- `E::Array{Float64, 2}`: Topocentric rotation matrix
"""
function sECEFtoSEZ(station_ecef::Array{<:Real, 1}, ecef::Array{<:Real, 1} ; conversion::String="geodetic")
    # Check input sizes
    if length(ecef) < 3
        throw(ArgumentError("Input ecef state must be at least length 3."))
    end

    if length(station_ecef) < 3
        throw(ArgumentError("Input station coordinates must be length 3."))
    end

    # Construct SEZ Rotation matrix
    E = rECEFtoSEZ(station_ecef, conversion=conversion)

    # Transform range
    range_ecef = ecef[1:3] - station_ecef
    range_sez  = E * range_ecef

    # Transform range-rate (if necessary)
    if length(ecef) == 6
        range_rate_ecef = ecef[4:6]
        range_rate_sez  = E * range_rate_ecef
    end

    # Return
    if length(ecef) == 6
        sez = vcat(range_sez, range_rate_sez)
    else
        sez = range_sez
    end

    return sez
end

export sSEZtoECEF
"""
Compute the coordinates of an object in the topocentric frame of an
Earth-fixed frame

Arguments:
- `station_ecef::Array{<:Real, 1}`: Earth-fixed cartesian station coordinates
- `sez::Array{<:Real, 1}`: SEZ coordinates of the object
- `conversion::Bool`: Conversion type to use. Can be "geocentric" or "geodetic"

Returns:
- `E::Array{Float64, 2}`: Topocentric rotation matrix
"""
function sSEZtoECEF(station_ecef::Array{<:Real, 1}, sez::Array{<:Real, 1} ; conversion::String="geodetic")
    # Check input sizes
    if length(sez) < 3
        throw(ArgumentError("Input SEZ state must be at least length 3."))
    end

    if length(station_ecef) < 3
        throw(ArgumentError("Input station coordinates must be length 3."))
    end

    # Compute ENZ Rotation matrix
    E = rSEZtoECEF(station_ecef, conversion=conversion)

    # Transform range
    range_sez  = sez[1:3]
    range_ecef = E * range_sez

    # Transform range-rate (if necessary)
    if length(sez) >= 6
        range_rate_sez  = sez[4:6]
        range_rate_ecef = E * range_rate_sez
    end

    # Return
    if length(sez) >= 6
        sat_ecef = vcat(range_ecef + station_ecef, range_rate_ecef)
    else
        sat_ecef = range_ecef + station_ecef
    end

    return sat_ecef
end

###############
# Topocentric #
###############

export sENZtoAZEL
"""
Convert East-North-Zenith topocentric state to azimuth, elevation, and range.

Arguments:
- `x::Array{<:Real, 1}`: East-North-Up state
- `use_degrees:Bool`: If `true` returns result in units of degrees

Returns:
- `azel::Array{<:Real, 1}`: Azimuth, elevation and range [rad; rad; m]
"""
function sENZtoAZEL(x::Array{<:Real, 1} ; use_degrees::Bool=false)
    # Check inputs
    if !(length(x) == 3 || length(x) == 6)
        throw(ArgumentError("Input ENZ state must be length 3 or 6."))
    end

    # Expand values
    rE, rN, rZ = x[1], x[2], x[3]

    # Range
    rho = norm(x[1:3])

    # Elevation
    el = atan(rZ, sqrt(rE^2 + rN^2))

    # Azimuth
    az = 0.0
    if el != pi/2 # Non-singular azimuth
        az = atan(rE, rN)
        if az < 0
            az += 2*pi
        end
    else # Azimuth may be singular for 90 deg elevation
        if length(x) != 6
            az = 0.0
            # @warn "Could not resolve singularity calculating azimuth."
        else
            # Use rate information to get azimuth if there is a singularity
            # in the position
            az = atan(x[4], x[5])
        end
    end

    # Output
    azel = [az ; el ; rho]

    if use_degrees
        azel[1] *= 180.0/pi
        azel[2] *= 180.0/pi
    end

    # Process Rate information
    if length(x) == 6
        rdE, rdN, rdZ = x[4], x[5], x[6]

        # Range-rate
        rhod = dot(x[1:3], x[4:6])/rho

        # Elevation-rate
        eld = (rdZ - norm(x[4:6])*sin(el))/sqrt(rE^2 + rN^2)

        # Azimuth-rate
        azd = (rdE*rN - rdN*rE)/(rE^2 + rN^2)

        # Output
        azel_rate = [azd ; eld ; rhod]
        if use_degrees
            azel_rate[1] *= 180/pi
            azel_rate[2] *= 180/pi
        end
    end

    # Return
    if length(x) == 6
        return vcat(azel, azel_rate)
    else
        return azel
    end
end

export sSEZtoAZEL
"""
Convert South-East-Zenith topocentric state to azimuth, elevation, and range.

Arguments:
- `x::Array{<:Real, 1}`: South-East-Zenith state
- `use_degrees:Bool`: If `true` returns result in units of degrees

Returns:
- `azel::Array{<:Real, 1}`: Azimuth, elevation and range [rad; rad; m]
"""
function sSEZtoAZEL(x::Array{<:Real, 1} ; use_degrees::Bool=false)
    # Check inputs
    if !(length(x) == 3 || length(x) == 6)
        throw(ArgumentError("Input rECEFtoSEZ state must be length 3 or 6."))
    end

    # Expand values
    rS, rE, rZ = x[1], x[2], x[3]

    # Range
    rho = norm(x[1:3])

    # Elevation
    el = atan(rZ, sqrt(rS^2 + rE^2))

    # Azimuth
    az = 0.0
    if el != pi/2
        az = atan(rE, -rS)
        if az < 0
            az += 2*pi
        end
    else
        if length(x) != 6
            az = 0.0
            # @warn "Could not resolve singularity calculating azimuth."
        else
            # Use rate information to get azimuth if there is a singularity
            # in the position
            az = atan(x[6], -x[5])
        end
    end

    # Output
    azel = [az ; el ; rho]

    if use_degrees
        azel[1] *= 180.0/pi
        azel[2] *= 180.0/pi
    end

    # Process Rate information
    if length(x) == 6
        rdS, rdE, rdZ = x[4], x[5], x[6]

        # Range-rate
        rhod = dot(x[1:3], x[4:6])/rho

        # Elevation-rate
        eld = (rdZ - norm(x[4:6])*sin(el))/sqrt(rS^2 + rE^2)

        # Azimuth-rate
        azd = (rdS*rE - rdE*rS)/(rS^2 + rE^2)

        # Output
        azel_rate = [azd ; eld ; rhod]
        if use_degrees
            azel_rate[1] *= 180/pi
            azel_rate[2] *= 180/pi
        end
    end

    # Return
    if length(x) == 6
        return vcat(azel, azel_rate)
    else
        return azel
    end
end


text:
#=
THis file contains all my mesh-manadgement methods.
=#

# mesh related functions
include("utils.jl")

# MESH MANAGEMENT
mutable struct mesh
    vtx::Array{Float64, 2} # vertices
    elt::Array{Array{Int64, 1}, 1} # elt → vtx
    edg::Array{Array{Int64, 1}, 1} # elt → edg
    freeedg::Array{Bool, 1} # elt → free edge
    freenod::Array{Bool, 1} # elt → free node
    bnd::Array{Int64, 1} # list of boundary nodes
    nedg::Int64 # number of edges
end

function searchCommonElt(vtx2elt::Array{Array{Int64,1},1}, n1::Int64,
    n2::Int64, egs::Int64)
    #= Search for the common element. This function may be shortened by making
        better use of the language… =#
    nv2e1 = length(vtx2elt[n1])
    nv2e2 = length(vtx2elt[n2])
    egs0 = 0
    egs1 = 0
    egs2 = 0
    ind0 = 0
    j    = 1
    while ind0 == 0 && j <= nv2e1
        egs1 = vtx2elt[n1][j] # j-th elt of node n1
        if egs1 != egs
            ind = 0; k = 1
            while ind == 0 && k <= nv2e2
                egs2 = vtx2elt[n2][k]  # the candidate
                if egs1 == egs2
                    ind = 1
                else
                    k += 1
                end
            end
            if ind == 1
                ind0 = 1
            else
                j += 1
            end
        else
            j += 1
        end
    end
    #
    if ind0 == 1
        if egs1 == egs2 && egs1 != egs
            egs0 = egs1
        else
            egs0 = 0
        end
    else
        egs0 = 0
    end
    #
    return egs0
end

function getNeighbours(vtx::Array{Float64, 2}, elt::Array{Array{Int64,1},1})
    #= This function returns the elements on each side of an edge. The value 0
        means that there is only one element ≡ the edge is free =#
    nvtx = size(vtx, 2)
    nelt = length(elt)
    # 1st: for each node, determine all neighbour elements
    vtx2elt = [Int64[] for i=1:nvtx]
    for ie=1:nelt # loop over all elements
        for ix in elt[ie] # loop over all nodes of each element
            push!(vtx2elt[ix], ie) # add the element to each of these nodes
        end
    end
    # 2nd: set the correspondance between the edges and the elements
    e2neighs = [Int64[] for i=1:nelt]
    for ie=1:nelt
        neltie = length(elt[ie])
        for k=1:neltie
            kp  = mod_wrap(k+1, neltie)
            nk  = elt[ie][k]
            nkp = elt[ie][kp]
            egs0 = searchCommonElt(vtx2elt, nk, nkp, ie)
            push!(e2neighs[ie], egs0)
        end
    end
    #
    return e2neighs, vtx2elt
end

function getNbEdges(e2neighs::Array{Array{Int64,1},1})
    #= Compute the total number of edges =#
    nelt = length(e2neighs)
    flag = zeros(Bool, nelt)
    nedg = Int64(0)
    for e=1:nelt
        neloc = length(e2neighs[e])
        for eloc=1:neloc
            tmp = e2neighs[e][eloc]
            if tmp != 0
                if flag[tmp] != true
                    nedg += 1
                end
            else
                nedg += 1
            end
        end
        flag[e] = true
    end
    #
    return nedg
end

function getMeshEdges(vtx::Array{Float64, 2}, elt::Array{Array{Int64, 1}, 1},
    nedg::Int64, e2neighs::Array{Array{Int64,1},1})
    #= This function assigns numbers to the edges of the mesh and return the
        correspondance between elements and edges with a sign =#
    nnod = size(vtx, 2)
    nelt = length(elt)
    #
    freeedg = zeros(Bool, nedg)
    freenod = zeros(Bool, nnod)
    nedg = 0 # reinitialize
    e2edg = [Int64[] for i=1:nelt]
    for e=1:nelt
        nodloc   = elt[e]
        e2edg[e] = zeros(Int64,length(nodloc))
    end
    #
    for e=1:nelt
        nodloc   = elt[e]
        xt       = vtx[:, nodloc]
        for eloc = 1:length(nodloc)
            if e2edg[e][eloc] != 0
                continue
            end
            if e2neighs[e][eloc] != 0
                elocp = mod_wrap(eloc+1, length(nodloc))
                xmt   = (xt[:,eloc] + xt[:, elocp])/2
                nodlocneigh = elt[e2neighs[e][eloc]]
                xl    = vtx[:, nodlocneigh]
                for u=1:length(nodlocneigh)
                    up  = mod_wrap(u+1, length(nodlocneigh))
                    xml = (xl[:,u] + xl[:, up])/2
                    if norm(xmt - xml) < 1e-8
                        nedg += 1
                        e2edg[e][eloc] = nedg
                        e2edg[e2neighs[e][eloc]][u] = nedg
                        break
                    end
                end
            else
                nedg += 1
                e2edg[e][eloc] = nedg
                freeedg[nedg]  = true
            end
        end
    end
    # free nodes
    for e=1:nelt
        neloc = length(elt[e])
        for eloc=1:neloc
            if freeedg[e2edg[e][eloc]]
                freenod[elt[e][eloc]] = true
                freenod[elt[e][mod_wrap(eloc+1, neloc)]] = true
            end
        end
    end
    return e2edg, freeedg, freenod
end

function readMesh(fname::String)
    println("Reading mesh file: $fname")
    vtx = zeros(Float64, 1, 1)
    elt = Array{Int64, 1}[];
    bnd = zeros(Int64, 1)
    open(fname, "r") do fid
        ln = readline(fid)
        println(" -> Reading vertices ...")
        # vertices
        nvtx = parse(Int64, ln)
        vtx = zeros(Float64, 2, nvtx)
        for i=1:nvtx
            ln = readline(fid);
            vtx[:, i] .= parse.(Float64, split(ln, " ", keepempty=false))
        end
        #elements
        println(" -> Reading elements ...")
        ln = readline(fid)
        nelt = parse(Int64, ln)
        for i=1:nelt
            ln = readline(fid)
            ln_s = split(ln, " ", keepempty=false)
            push!(elt, parse.(Int64, ln_s))
        end
        # boundary
        println(" -> Reading boundary nodes ...")
        ln = readline(fid)
        nbnd = parse(Int64, ln)
        bnd = zeros(Int64, nbnd)
        for i=1:nbnd
            ln = readline(fid)
            bnd[i] = parse(Int64, ln)
        end
    end
    #
    e2neigh, _ = getNeighbours(vtx, elt)
    nedg    = getNbEdges(e2neigh)
    e2edg, freeedg, freenod = getMeshEdges(vtx, elt, nedg, e2neigh)
    return mesh(vtx, elt, e2edg, freeedg, freenod, bnd, nedg)
end

# compute connectivity table
function connectivityTableOrde0(msh::mesh)
    #= This function returns the connectivity table corresponding to a problem
        of order k = 0. =#
end


text:
# Sample postprocessing script for customized figure.
#
# Reading 2D simulation data, zooming-in to a region of interest,
# adding streamlines and contour at specific levels on top of colored mesh.
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, PyPlot
using Vlasiator: RE

file = "bulk.0000501.vlsv"
nameρ = "rho"
nameV = "rho_v"

meta = load(file)

fig, ax = plt.subplots(1,1; figsize=(8,6), constrained_layout=true)

pcolormesh(meta, nameρ, ax)
streamplot(meta, nameV, ax; comp="xy", color="w", density=2.0)

cbar = ax.collections[end].colorbar
boxcoords = [0, 20, -15, 15]
ax.set_xlim([boxcoords[1],boxcoords[2]])
ax.set_ylim([boxcoords[3],boxcoords[4]])

# Contour line at a specific level
pArgs = Vlasiator.set_args(meta, nameρ, EARTH)
x, y = Vlasiator.get_axis(pArgs)
data = Vlasiator.prep2d(meta, nameρ)'
CS = plt.contour(x, y, data, levels = [1e7],
                 colors=("k",),linestyles=("-",),linewidths=(0.5,))

# Add a rectangular box region
boxrange = (250:299, 200:249)

x = LinRange(meta.coordmin[1], meta.coordmax[1], meta.ncells[1])
y = LinRange(meta.coordmin[2], meta.coordmax[2], meta.ncells[2])
# Construct the patch
rect = matplotlib.patches.Rectangle((x[boxrange[1][1]]/RE, y[boxrange[2][1]]/RE),
   (x[boxrange[1][end]] - x[boxrange[1][1]])/RE,
   (y[boxrange[2][end]] - y[boxrange[2][1]])/RE,
   linewidth=1, edgecolor="r", facecolor="none")

ax.add_patch(rect)


text:
module CarlemanLinearization

# dependencies
include("init.jl")

# methods
include("kronecker.jl")
include("linearization.jl")
include("error_bounds.jl")
include("conversion.jl")

# exported methods and types
include("exports.jl")

end # module


text:
__precompile__(true)

"""
Main module of `Expokit` -- Julia implementations of some routines contained
in [EXPOKIT](https://www.maths.uq.edu.au/expokit/).

Those routines allow an efficient calculation of the action of matrix exponentials
on vectors for large sparse matrices. For more details about the methods see
R.B. Sidje, ACM Trans. Math. Softw., 24(1):130-156, 1998 (or its
[preprint](https://www.maths.uq.edu.au/expokit/paper.pdf)).
"""
module Expokit

using LinearAlgebra

const axpy! = LinearAlgebra.axpy!
const exp! = LinearAlgebra.exp!

include("expmv.jl")
include("padm.jl")
include("chbv.jl")
include("phimv.jl")

end # module


text:
abstract type AbstractValEdge{V <: Integer, E_VALS} <: AbstractEdge{V} end

"""
    ValEdge{V, E_VALS} <: AbstractValEdge{V}

A datastructure representing an undirected edge with multiple values.
"""
struct ValEdge{V<:Integer, E_VALS} <: AbstractValEdge{V, E_VALS}
    src::V
    dst::V
    vals::E_VALS

    function ValEdge(src::V, dst::V, vals::E_VALS) where {V, E_VALS}
        src, dst = minmax(src, dst) # TODO maybe use a branchless operator
        return new{V, E_VALS}(src, dst, vals)
    end
end

# TODO update docstring and ValDiEdge
"""
    ValEdge(s, d, v)
Create a `ValEdge` with source `s`, destination `d` and value `v`.
# Examples
```
julia> e = ValEdge(4, 2, 'A')
Edge 4 => 2 with value 'A'
```
"""
function ValEdge end

"""
    ValDiEdge{V, E_VALS} <: AbstractValEdge{V}

A datastructure representing a directed edge with values.
"""
struct ValDiEdge{V<:Integer, E_VALS} <: AbstractValEdge{V, E_VALS}
    src::V
    dst::V
    vals::E_VALS

    function ValDiEdge(src::V, dst::V, vals::E_VALS) where {V, E_VALS}
        return new{V, E_VALS}(src, dst, vals)
    end
end

src(e::AbstractValEdge) = e.src
dst(e::AbstractValEdge) = e.dst

"""
    vals(e::ValEdge)
Returns the value attached to the edge `e`.
# Examples
```
julia> g = EdgeValGraph(3, String)

julia> add_edge!(g, 2, 3, ("xyz",))

julia> first(edges(g)) |> vals
"xyz"
```
"""
vals(e::AbstractValEdge) = e.vals

# TODO Maybe that should be declared somewhere else
const SingleValTuple{T} = Union{Tuple{T}, NamedTuple{S, Tuple{T}} where {S}}

val(e::AbstractValEdge; key::Union{Integer, Symbol, NoKey}=nokey) =
    _val(e, key)

_val(e::AbstractValEdge{V, E_VAL}, ::NoKey) where {V, E_VAL <: SingleValTuple} =
    e.vals[1]

_val(e::AbstractValEdge, key::Union{Integer, Symbol}) = e.vals[key]

reverse(e::ValEdge) = e
reverse(e::ValDiEdge) = ValDiEdge(dst(e), src(e), vals(e))

is_directed(::Type{<:ValEdge}) = false
is_directed(::Type{<:ValDiEdge}) = true
is_directed(e::AbstractValEdge) = is_directed(typeof(e))

function show(io::IO, e::AbstractValEdge)
    isdir = is_directed(e)
    e_keys = keys(vals(e))
    has_symbol_keys = eltype(e_keys) === Symbol

    print(io, isdir ? "ValDiEdge" : "ValEdge")
    arrow = isdir ? "->" : "--"
    print(io, " $(src(e)) $arrow $(dst(e))")

    if length(e_keys) == 1
        print(io, " with value " * (has_symbol_keys ? "$(e_keys[1]) = $(val(e))" :
                                                      "$(val(e))"))
    elseif length(e_keys) > 1
        print(io, " with values $(vals(e))")
    end

    println(io)
end


text:
abstract type AbstractImageNet224 <: Name end
struct ImageNet224MMAP <: AbstractImageNet224 end
struct ImageNet224 <: AbstractImageNet224 end

const IMAGENET_PATH = Ref{String}(joinpath(DATADIR, "ImageNet224.h5"))
const IMAGENET_FILE = Ref{HDF5.File}(HDF5.File(-1, "nnn"))

function set_imagenet_path!(path)
    IMAGENET_PATH[] = path
end

function open_imagenet()
    if !isfile(IMAGENET_PATH[])
        error("File $(IMAGENET_PATH[]) does not exists. Use DatasetProvider.set_imagenet_path! function to set the proper path.")
    end
    fid = IMAGENET_FILE[]
    if !isopen(fid)
        IMAGENET_FILE[] = fid = h5open(IMAGENET_PATH[], "r")
    end
    return fid
end

close_imagenet() = close(IMAGENET_FILE[])

# download options
task(::Type{<:AbstractImageNet224}) = MultiClass
format(::Type{<:AbstractImageNet224}) = ColorImages
source(::Type{<:AbstractImageNet224}) = "https://image-net.org/"
version(::Type{<:AbstractImageNet224}) = "ILSVRC2012"
make_datadep(::Type{<:AbstractImageNet224}) = nothing

nclasses(::Type{<:AbstractImageNet224}) = 1000
classes(::Type{<:AbstractImageNet224}) = collect(1:1000)
nattributes(::Type{<:AbstractImageNet224}) = (224, 224, 3)
ninstances(::Type{<:AbstractImageNet224}) = (1281167, 0, 50000)

# dataset description
name(::Type{<:AbstractImageNet224}) = "ImageNet Large Scale Visual Recognition Challenge 2012"

function author(::Type{<:AbstractImageNet224})
    return ["Olga Russakovsky", "Jia Deng", "Hao Su", "Jonathan Krause", "Sanjeev Satheesh", "Sean Ma", "Zhiheng Huang", "Andrej Karpathy", "Aditya Khosla", "Michael Bernstein", "Alexander C. Berg", "Li Fei-Fei"]
end

function citation(::Type{<:AbstractImageNet224})
    return """
    @article{ILSVRC15,
        Author = {Olga Russakovsky and Jia Deng and Hao Su and Jonathan Krause and Sanjeev Satheesh and Sean Ma and Zhiheng Huang and Andrej Karpathy and Aditya Khosla and Michael Bernstein and Alexander C. Berg and Li Fei-Fei},
        Title = {{ImageNet Large Scale Visual Recognition Challenge}},
        Year = {2015},
        journal   = {International Journal of Computer Vision (IJCV)},
        doi = {10.1007/s11263-015-0816-y},
        volume={115},
        number={3},
        pages={211-252}
    }
    """
end

function loadraw(N::Type{<:AbstractImageNet224}, type)
    hassubset(N, type)
    split = type == :test ? "val" : "train"

    fid = open_imagenet()
    if N <: ImageNet224MMAP
        data = HDF5.readmmap(fid["$(split)_data"]) |> HDF5Array
    else
        @warn "due to the size of the ImageNet dataset, it may take several minutes to load data into memory"
        data = HDF5.read(fid["$(split)_data"])
    end
    targets = fid["$(split)_targets"][:]
    close(fid)
    return data, targets
end


text:
"""
$(DocStringExtensions.README)
"""
module NeuralPDE

using DocStringExtensions
using Reexport, Statistics
@reexport using DiffEqBase

using Flux, Zygote, DiffEqSensitivity, ForwardDiff, Random, Distributions
using DiffEqFlux, Adapt, DiffEqNoiseProcess, CUDA, StochasticDiffEq
using ModelingToolkit
using GalacticOptim
using Quadrature,Cubature
using QuasiMonteCarlo
using RuntimeGeneratedFunctions
using SciMLBase
using Statistics
using ArrayInterface
import Tracker, Optim
using DomainSets
using Symbolics
import ModelingToolkit: value, nameof, toexpr, build_expr, expand_derivatives
import DomainSets: Domain, ClosedInterval
import ModelingToolkit: Interval, infimum, supremum #,Ball
import SciMLBase: @add_kwonly
using Flux: @nograd
abstract type NeuralPDEAlgorithm <: DiffEqBase.AbstractODEAlgorithm end
"""
    TerminalPDEProblem(g, f, μ, σ, x0, tspan)
A semilinear parabolic PDE problem with a terminal condition.
Consider `du/dt = l(u) + f(u)`; where l is the nonlinear Lipschitz function
# Arguments
* `g` : The terminal condition for the equation.
* `f` : The function f(u)
* `μ` : The drift function of X from Ito's Lemma
* `μ` : The noise function of X from Ito's Lemma
* `x0`: The initial X for the problem.
* `tspan`: The timespan of the problem.
"""
struct TerminalPDEProblem{G,F,Mu,Sigma,X,T,P,A,UD,K} <: SciMLBase.SciMLProblem
    g::G
    f::F
    μ::Mu
    σ::Sigma
    X0::X
    tspan::Tuple{T,T}
    p::P
    A::A
    u_domain::UD
    kwargs::K
    TerminalPDEProblem(g,f,μ,σ,X0,tspan,p=nothing;A=nothing,u_domain=nothing,kwargs...) = new{typeof(g),typeof(f),
                                                         typeof(μ),typeof(σ),
                                                         typeof(X0),eltype(tspan),
                                                         typeof(p),typeof(A),typeof(u_domain),typeof(kwargs)}(
                                                         g,f,μ,σ,X0,tspan,p,A,u_domain,kwargs)
end

Base.summary(prob::TerminalPDEProblem) = string(nameof(typeof(prob)))

function Base.show(io::IO, A::TerminalPDEProblem)
  println(io,summary(A))
  print(io,"timespan: ")
  show(io,A.tspan)
end

"""
    KolmogorovPDEProblem(f , g, phi , xspan , tspan, d, noise_rate_prototype = none)
A standard Kolmogorov PDE Problem.
# Arguments
* `f` : The drift function from Feynman-Kac representation of the PDE.
* `g` : The noise function from Feynman-Kac representation of the PDE.
* `phi` : The terminal condition for the PDE.
* `tspan`: The timespan for the problem.
* `xspan`: The xspan for the problem.
* `d`: The dimensions of the input x.
* `noise_rate_prototype`: A prototype type instance for the noise rates, that is the output g.
"""
struct KolmogorovPDEProblem{ F, G, Phi, X , T , D ,P,U0, ND} <: DiffEqBase.DEProblem
    f::F
    g::G
    phi::Phi
    xspan::Tuple{X,X}
    tspan::Tuple{T,T}
    d::D
    p::P
    u0::U0
    noise_rate_prototype::ND
    KolmogorovPDEProblem( f, g, phi , xspan , tspan , d, p=nothing, u0=0 , noise_rate_prototype= nothing) = new{typeof(f),typeof(g),typeof(phi),eltype(tspan),eltype(xspan),typeof(d),typeof(p),typeof(u0),typeof(noise_rate_prototype)}(f,g,phi,xspan,tspan,d,p,u0,noise_rate_prototype)
end

Base.summary(prob::KolmogorovPDEProblem) = string(nameof(typeof(prob)))
function Base.show(io::IO, A::KolmogorovPDEProblem)
  println(io,summary(A))
  print(io,"timespan: ")
  show(io,A.tspan)
  print(io,"xspan: ")
  show(io,A.xspan)
  println(io , "μ")
  show(io , A.f)
  println(io,"Sigma")
  show(io , A.g)
end

abstract type ParametersDomain end
struct KolmogorovParamDomain{T} <:ParametersDomain
           sigma::Tuple{T,T}
           mu::Tuple{T,T}
           phi::Tuple{T,T}
end

struct ParamKolmogorovPDEProblem{ F, G, Phi, X , T , D, YD , P,U0 , YSP , YMP ,YPH, NP} <: DiffEqBase.DEProblem
    f::F
    g::G
    phi::Phi
    xspan::Tuple{X,X}
    tspan::Tuple{T,T}
    d::D
    Y_domain::YD
    p::P
    u0::U0
    Y_sigma_prototype::YSP
    Y_mu_prototype::YMP
    Y_phi_prototype::YPH
    noise_rate_prototype::NP
    ParamKolmogorovPDEProblem( f, g, phi , xspan , tspan , d, Y_domain, p=nothing, u0=0 ; Y_sigma_prototype=nothing , Y_mu_prototype=nothing , Y_phi_prototype=nothing , noise_rate_prototype= nothing) = new{typeof(f),typeof(g),typeof(phi),eltype(tspan),eltype(xspan),typeof(d),typeof(Y_domain),typeof(p),typeof(u0),typeof(Y_sigma_prototype),typeof(Y_mu_prototype),typeof(Y_phi_prototype),typeof(noise_rate_prototype)}(f,g,phi,xspan,tspan,d,Y_domain,p,u0 , Y_sigma_prototype,Y_mu_prototype,Y_phi_prototype,noise_rate_prototype)
end

Base.summary(prob::ParamKolmogorovPDEProblem) = string(nameof(typeof(prob)))
function Base.show(io::IO, A::ParamKolmogorovPDEProblem)
  println(io,summary(A))
  print(io,"timespan: ")
  show(io,A.tspan)
  print(io,"xspan: ")
  show(io,A.xspan)
  println(io , "μ")
  show(io , A.f)
  println(io,"Sigma")
  show(io , A.g)
end

"""
Algorithm for solving differential equation using a neural network.
Arguments:
* `chain`: A Chain neural network
* `opt`: The optimizer to train the neural network. Defaults to `BFGS()`.
* `initθ`: The initial parameter of the neural network.
* `autodiff`: The switch between automatic and numerical differentiation for
              the PDE operators. The reverse mode of the loss function is always AD.
"""
struct NNODE{C,O,P,K} <: NeuralPDEAlgorithm
    chain::C
    opt::O
    initθ::P
    autodiff::Bool
    kwargs::K
end
function NNODE(chain,opt=Optim.BFGS(),init_params = nothing;autodiff=false,kwargs...)
    if init_params === nothing
        if chain isa FastChain
            initθ = DiffEqFlux.initial_params(chain)
        else
            initθ,re  = Flux.destructure(chain)
        end
    else
        initθ = init_params
    end
    NNODE(chain,opt,initθ,autodiff,kwargs)
end

include("ode_solve.jl")
include("pde_solve.jl")
include("pde_solve_ns.jl")
include("kolmogorov_solve.jl")
include("rode_solve.jl")
include("stopping_solve.jl")
include("transform_inf_integral.jl")
include("pinns_pde_solve.jl")
include("neural_adapter.jl")
include("param_kolmogorov_solve.jl")

export NNODE, TerminalPDEProblem, NNPDEHan, NNPDENS, NNRODE,
       KolmogorovPDEProblem, NNKolmogorov, NNStopping,ParamKolmogorovPDEProblem,KolmogorovParamDomain, NNParamKolmogorov,
       PhysicsInformedNN, discretize,
       GridTraining, StochasticTraining, QuadratureTraining, QuasiRandomTraining,
       build_loss_function, get_loss_function,
       generate_training_sets, get_variables, get_argument, get_bounds,
       get_phi, get_numeric_derivative, get_numeric_integral,
       build_symbolic_equation, build_symbolic_loss_function, symbolic_discretize,
       AbstractAdaptiveLoss, NonAdaptiveLoss, GradientScaleAdaptiveLoss, MiniMaxAdaptiveLoss,
       LogOptions

end # module


text:
using Test

@testset "preimages" begin

# Testing skip_beginning / skip_end

a1 = LinRange(0,1,12)
a2 = [0, 0..0.1, 0..0.2, 0.1..0.3, 0.2..0.4, 0.3..0.5, 0.3..0.6, 0.6..0.8, 0.8..1, 1..1] #tricky weakly sorted vector

x1 = 0.25..0.25
x2 = 0.1..0.8
x3 = 2
x4 = -1
x5 = 0
x6 = 0.3

for a = (a1, a2)
    for x = (x1, x2, x3, x4, x5,  x6)
        i = InvariantMeasures.first_overlapping(a, x)
        @test i==0 || Interval(a[i]).hi <= Interval(x).lo
        @test i==length(a) || !(Interval(a[i+1]).hi <= Interval(x).lo)

        j = InvariantMeasures.last_overlapping(a, x)
        @test j==0 || !(Interval(a[j]).lo >= Interval(x).hi)
        @test j==length(a) || Interval(a[j+1]).lo >= Interval(x).hi
    end
end

# checking for equality with looser tolerance
function approxintervals(a, b)
    return a.lo ≈ b.lo && a.hi ≈ b.hi
end

for y = (a1, a2)
    for f in (x->x/2, x->1+x/2, x-> 1-x/2)
        b = InvariantMeasures.Branch(f, (@interval(0), @interval(1)))
        ylabel = 1:length(y)
        x, xlabel = preimages(y, b, ylabel)
        @test x[1] == b.X[1]
        @test x[end] != b.X[2] # to make sure the last entry isn't there
        @test length(x) == length(xlabel)
        y1 = filter(x->!isempty(x),intersect.(map(Interval,y), hull(b.Y[1], b.Y[2]+1e-15))) # the 1e-15 is there to simulate "disjoint intervals"
        y2 = f.(x)
        if !b.increasing
            y2 = reverse(y2)
        end
        @test all(approxintervals.(y1, y2))
    end
end

D = mod1_dynamic(x -> 2x)
DD = ∘(D, D, D, D)
p = [0, 0.2, 0.4, 0.6, 0.8]

x, xlabel = preimages(p, DD)
@test x ≈ 0:1/80:79/80
@test xlabel == repeat([1,2,3,4,5],16)

# test composing two different functions, to check that composition is handled in the correct order

f = x-> 2x
g = x -> (x+x^2)/2

D1 = mod1_dynamic(f)
D2 = mod1_dynamic(g)

y = 0:0.2:1
x, xlabel = preimages(y, D1 ∘ D2)

@test f.(g.(x)) ≈ 0:0.2:1.8
@test xlabel ≈ repeat(1:5, 2)

D1 = PwMap([x->2x, x->6x-3, x->3x-2], [0, 0.5, @interval(2/3), 1], [0 1; 0 1; 0 1])
D2 = PwMap([x->2x, x->4x-2, x->4x-3], [0, 0.5, 0.75, 1], [0 1; 0 1; 0 1])

z = 0:0.3:1
y, ylabel, y′ = InvariantMeasures.preimages_and_derivatives(z, D1)

@test y ≈ [0, 0.15, 0.3, 0.45, 0.5, 0.5+0.05, 0.5+0.1, 0.5+0.15, 2/3, 2/3+0.1, 2/3+0.2, 2/3+0.3]
@test ylabel == repeat(1:4, 3)
@test y′ ≈ [2,2,2,2,6,6,6,6,3,3,3,3]

x, xlabel, x′ = InvariantMeasures.preimages_and_derivatives(z, D1∘D2)
@test all(x .≈ vcat(y/2, 0.5 .+ y/4, 0.75 .+ y/4))
@test x′ == kron([4,12,6,8,24,12,8,24,12], [1,1,1,1])

end #testset


text:
# DiagTree = GWKT.DiagTree
# Parquet = GWKT.Parquet
# Gtype, Wtype = 1, 2

MUL, ADD = 1, 2
"""
Build tree with KinL = KoutL = [1, 0, 0, 0], KinR = KoutR = [0, 1, 0]
"""
function build(chan, legK, kidx, spin, irreducible, hasBubble, Gsym, Wsym)
    diag = DiagTree.Diagrams{Float64}()
    dir, ex = [], []
    if 1 in chan
        Tdir, Tex = buildTU(diag, legK, kidx, spin, irreducible, hasBubble, Gsym, Wsym, true)
        append!(dir, Tdir)
        append!(ex, Tex)
    end

    if 2 in chan
        Udir, Uex = buildTU(diag, legK, kidx, spin, irreducible, hasBubble, Gsym, Wsym, false)
        append!(dir, Udir)
        append!(ex, Uex)
    end

    if 3 in chan
        Sdir, Sex = buildS(diag, legK, kidx, Gsym, Wsym)
        append!(dir, Sdir)
        append!(ex, Sex)
    end

    rootDir = DiagTree.addNode!(diag, ADD, 1.0, [], dir, isRoot = true)
    rootEx = DiagTree.addNode!(diag, ADD, 1.0, [], ex, isRoot = true)
    return diag, rootDir, rootEx
end

function addG(diag, K, T::Vector{Tuple{Int,Int}}, Gsym)
    return [DiagTree.addPropagator!(diag, 1, 0, K[i], T[i], Gsym)[1] for i = 1:length(T)] # G order is 0
end

function addV(diag, K, T::Tuple{Int,Int}, Wsym) #direct and exchange V
    factor = [1.0, -1.0]
    return [DiagTree.addPropagator!(diag, 2, 1, K[i], T, Wsym, factor[i])[1] for i = 1:2] # V order is 1
end

function addW(diag, K, T::Tuple{Int,Int}, Wsym) #direct and exchange W, both have the same tau variables
    factor = [1.0, -1.0]
    return [DiagTree.addPropagator!(diag, 3, 1, K[i], T, Wsym, factor[i])[1] for i = 1:2] # W order is 1
end

function Tpair(Tidx, isBare, isDirect)
    if isBare
        return (Tidx, Tidx, Tidx, Tidx)
    else
        if isDirect
            return (Tidx, Tidx, Tidx + 1, Tidx + 1)
        else
            return (Tidx, Tidx + 1, Tidx + 1, Tidx)
        end
    end
end

function buildTU(diag, legK, kidx, spin, irreducible, hasBubble, Gsym, Wsym, isT)
    """
        k1-q                      k2+q
        |                         |
    t1.L ↑     t1.L       t2.L     ↑ t2.L
        |-------------->----------|
        |       |    k4    |      |
        |   v   |          |  v   |
        |       |  k4-q    |      |
        |--------------<----------|
    t1.L ↑    t1.L        t2.L     ↑ t2.L
        |                         |
        k1                        k2
    """
    Sym = isT ? -1.0 : 1.0
    INL, OUTL, INR, OUTR = 1, 2, 3, 4
    D, E = 1, 2
    KinL, KoutL, KinR, KoutR = legK

    if isT == false
        KoutL, KoutR = KoutR, KoutL
    end

    qd = KinL - KoutL
    K = zero(KinL)
    K[kidx] = 1
    Kc = K - qd
    Td, Te = [], []

    function map(isLbare, isRbare, isLdirect, isRdirect)
        Lt = Tpair(1, isLbare, isLdirect)
        Rt = Tpair(3, isRbare, isRdirect)

        extT = isT ? (Lt[INL], Lt[OUTL], Rt[INR], Rt[OUTR]) : (Lt[INL], Rt[OUTR], Rt[INR], Lt[OUTL])
        # construct tau table for Green's functions, e.g, (1, 3) means G(t3-t1)
        gT = [(Lt[OUTR], Rt[INL]), (Rt[OUTL], Lt[INR])]
        gK = [K, Kc]
        g = addG(diag, gK, gT, Gsym)

        #construct the Momentum table, momentum configurations are independent of tau
        LwK, RwK = [qd, KinL - K], [qd, KinR - Kc]
        Lw = isLbare ? addV(diag, LwK, (1, 2), Wsym) : addW(diag, LwK, (1, 2), Wsym)
        Rw = isRbare ? addV(diag, RwK, (3, 4), Wsym) : addW(diag, RwK, (3, 4), Wsym)
        return g, Lw, Rw, extT
    end

    function makeTree!(isLbare, isRbare)
        if irreducible == false || isT == false #U channel diagram always has those terms
            if hasBubble
                g, Lw, Rw, extT = map(isLbare, isRbare, true, true)
                push!(Td, DiagTree.addNode!(diag, MUL, spin * Sym, [g[1], g[2], Lw[D], Rw[D]], [], extT = extT))
            end
            g, Lw, Rw, extT = map(isLbare, isRbare, true, false)
            vde = DiagTree.addNode!(diag, MUL, Sym, [g[1], g[2], Lw[D], Rw[E]], [], extT = extT)
            g, Lw, Rw, extT = map(isLbare, isRbare, false, true)
            ved = DiagTree.addNode!(diag, MUL, Sym, [g[1], g[2], Lw[E], Rw[D]], [], extT = extT)
            append!(Td, [vde, ved])
            # push!(Td, DiagTree.addNode!(diag, ADD, 1.0, [], [vdd, vde, ved]))
        end
        g, Lw, Rw, extT = map(isLbare, isRbare, false, false)
        push!(Te, DiagTree.addNode!(diag, MUL, Sym, [g[1], g[2], Lw[E], Rw[E]], []; extT = extT))
    end

    # construct the propagator table
    makeTree!(true, true) #vxv
    makeTree!(true, false) #vxw
    makeTree!(false, true) #wxv
    makeTree!(false, false) #wxw
    if isT
        return Td, Te
    else
        return Te, Td
    end
end

# function buildTUC(diag, legK, kidx, spin, irreducible, hasBubble, Gsym, Wsym, isT)
#     """
#         k1-q                      k2+q
#         |                         |
#     t1.L ↑     t1.L       t2.L     ↑ t2.L
#         |-------------->----------|
#         |       |    k4    |      |
#         |   v   |          |  v   |
#         |       |  k4-q    |      |
#         |--------------<----------|
#     t1.L ↑    t1.L        t2.L     ↑ t2.L
#         |                         |
#         k1                        k2
#     """
#     Sym = isT ? -1.0 : 1.0
#     INL, OUTL, INR, OUTR = 1, 2, 3, 4
#     D, E = 1, 2
#     KinL, KoutL, KinR, KoutR = legK

#     if isT == false
#         KoutL, KoutR = KoutR, KoutL
#     end

#     qd = KinL - KoutL
#     K = zero(KinL)
#     K[kidx] = 1
#     Kc = K - qd
#     Td, Te = [], []

#     function map(isLbare, isRbare, isLdirect, isRdirect)
#         Lt = Tpair(1, isLbare, isLdirect)
#         Rt = Tpair(3, isRbare, isRdirect)

#         extT = isT ? (Lt[INL], Lt[INL], Rt[INR], Rt[INR]) : (Lt[INL], Rt[INR], Rt[INR], Lt[INL])
#         # construct tau table for Green's functions, e.g, (1, 3) means G(t3-t1)
#         gT = [(Lt[OUTR], Rt[INL]), (Rt[OUTL], Lt[INR])]
#         gK = [K, K]
#         g = addG(diag, gK, gT, Gsym)

#         #construct the Momentum table, momentum configurations are independent of tau
#         LwK, RwK = [qd, KinL - K], [qd, KinR - Kc]
#         Lw = isLbare ? addV(diag, LwK, (1, 2), Wsym) : addW(diag, LwK, (1, 2), Wsym)
#         Rw = isRbare ? addV(diag, RwK, (3, 4), Wsym) : addW(diag, RwK, (3, 4), Wsym)
#         return g, Lw, Rw, extT
#     end

#     function makeTree!(isLbare, isRbare)
#         if irreducible == false || isT == false #U channel diagram always has those terms
#             if hasBubble
#                 g, Lw, Rw, extT = map(isLbare, isRbare, true, true)
#                 push!(Td, DiagTree.addNode!(diag, MUL, spin * Sym, [g[1], g[2], Lw[D], Rw[D]], [], extT = extT))
#             end
#             g, Lw, Rw, extT = map(isLbare, isRbare, true, false)
#             vde = DiagTree.addNode!(diag, MUL, Sym, [g[1], g[2], Lw[D], Rw[E]], [], extT = extT)
#             g, Lw, Rw, extT = map(isLbare, isRbare, false, true)
#             ved = DiagTree.addNode!(diag, MUL, Sym, [g[1], g[2], Lw[E], Rw[D]], [], extT = extT)
#             append!(Td, [vde, ved])
#             # push!(Td, DiagTree.addNode!(diag, ADD, 1.0, [], [vdd, vde, ved]))
#         end
#         g, Lw, Rw, extT = map(isLbare, isRbare, false, false)
#         push!(Te, DiagTree.addNode!(diag, MUL, Sym, [g[1], g[2], Lw[E], Rw[E]], []; extT = extT))
#     end

#     # construct the propagator table
#     makeTree!(true, true) #vxv
#     makeTree!(true, false) #vxw
#     makeTree!(false, true) #wxv
#     makeTree!(false, false) #wxw
#     if isT
#         return Td, Te
#     else
#         return Te, Td
#     end
# end

function buildS(diag, legK, kidx, Gsym, Wsym)
    # Sym = -0.5
    Sym = -1.0
    INL, OUTL, INR, OUTR = 1, 2, 3, 4
    D, E = 1, 2
    KinL, KoutL, KinR, KoutR = legK

    qd = KinL - KoutL
    K = zero(KinL)
    K[kidx] = 1
    Ks = KinL + KinR - K
    Sd, Se = [], []

    function map(isLbare, isRbare, isLdirect, isRdirect)
        Lt = Tpair(1, isLbare, isLdirect)
        Rt = Tpair(3, isRbare, isRdirect)

        extT = (Lt[INL], Rt[OUTL], Lt[INR], Rt[OUTR])
        # construct tau table for Green's functions, e.g, (1, 3) means G(t3-t1)
        gT = [(Lt[OUTR], Rt[INL]), (Lt[OUTL], Rt[INR])]
        gK = [K, Ks]
        g = addG(diag, gK, gT, Gsym)

        #construct the Momentum table, momentum configurations are independent of tau
        # LLegK = [KinL, Ks, KinR, K]
        # RLegK = [K, KoutL, Ks, KoutR]
        LwK, RwK = [KinL - Ks, KinL - K], [K - KoutL, K - KoutR]
        Lw = isLbare ? addV(diag, LwK, (1, 2), Wsym) : addW(diag, LwK, (1, 2), Wsym)
        Rw = isRbare ? addV(diag, RwK, (3, 4), Wsym) : addW(diag, RwK, (3, 4), Wsym)
        return g, Lw, Rw, extT
    end

    function makeTree!(isLbare, isRbare)
        # g, Lw, Rw, extT = map(isLbare, isRbare, true, false)
        # push!(Sd, DiagTree.addNode!(diag, MUL, Sym, [g[1], g[2], Lw[D], Rw[E]], [], extT = extT))
        g, Lw, Rw, extT = map(isLbare, isRbare, false, true)
        push!(Sd, DiagTree.addNode!(diag, MUL, Sym, [g[1], g[2], Lw[E], Rw[D]], [], extT = extT))

        g, Lw, Rw, extT = map(isLbare, isRbare, true, true)
        push!(Se, DiagTree.addNode!(diag, MUL, Sym, [g[1], g[2], Lw[D], Rw[D]], [], extT = extT))
        # g, Lw, Rw, extT = map(isLbare, isRbare, false, false)
        # push!(Se, DiagTree.addNode!(diag, MUL, Sym, [g[1], g[2], Lw[E], Rw[E]], [], extT = extT))
    end
    # construct the propagator table
    makeTree!(true, true) #vxv
    makeTree!(true, false) #vxw
    makeTree!(false, true) #wxv
    makeTree!(false, false) #wxw

    return Sd, Se
end

if abspath(PROGRAM_FILE) == @__FILE__
    KinL = KoutL = [1, 0, 0, 0]
    KinR = KoutR = [0, 1, 0, 0]
    legK = [KinL, KoutL, KinR, KoutR]
    Gsym = [:mirror]
    Wsym = [:mirror, :timereversal]

    diag = DiagTree.Diagrams{Float64}()
    Td, Te = buildT(diag, legK, 2, 4, false, Gsym, Wsym)
    for t in Td
        DiagTree.showTree(diag, t)
    end
    for t in Te
        DiagTree.showTree(diag, t)
    end
end


text:
    begin
    local n,a,asym,d,v
    n = 10
    areal  = randn(n,n)
    acmplx = complex(randn(n,n), randn(n,n))

    for elty in (Float32, Float64, Complex64, Complex128)
        if elty == Complex64 || elty == Complex128
            a = acmplx
        else
            a = areal
        end
        a     = convert(Matrix{elty}, a)
        asym  = a' + a                  # symmetric indefinite
        apd   = a'*a                    # symmetric positive-definite

	(d,v) = eigs(a, nev=3)
	@test_approx_eq a*v[:,2] d[2]*v[:,2]

	(d,v) = eigs(asym, nev=3)
	@test_approx_eq asym*v[:,1] d[1]*v[:,1]
#        @test_approx_eq eigs(asym; nev=1, sigma=d[3])[1][1] d[3]

	(d,v) = eigs(apd, nev=3)
	@test_approx_eq apd*v[:,3] d[3]*v[:,3]
#        @test_approx_eq eigs(apd; nev=1, sigma=d[3])[1][1] d[3]

    # test (shift-and-)invert mode
    (d,v) = eigs(apd, nev=3, sigma=0)
    @test_approx_eq apd*v[:,3] d[3]*v[:,3]

    end
end

# Example from Quantum Information Theory
import Base: size, issym, ishermitian

type CPM{T<:Base.LinAlg.BlasFloat}<:AbstractMatrix{T} # completely positive map
	kraus::Array{T,3} # kraus operator representation
end

size(Phi::CPM)=(size(Phi.kraus,1)^2,size(Phi.kraus,3)^2)
issym(Phi::CPM)=false
ishermitian(Phi::CPM)=false

function *{T<:Base.LinAlg.BlasFloat}(Phi::CPM{T},rho::Vector{T})
	rho=reshape(rho,(size(Phi.kraus,3),size(Phi.kraus,3)))
	rho2=zeros(T,(size(Phi.kraus,1),size(Phi.kraus,1)))
	for s=1:size(Phi.kraus,2)
		As=slice(Phi.kraus,:,s,:)
		rho2+=As*rho*As'
	end
	return reshape(rho2,(size(Phi.kraus,1)^2,))
end
# Generate random isometry
(Q,R)=qr(randn(100,50))
Q=reshape(Q,(50,2,50))
# Construct trace-preserving completely positive map from this
Phi=CPM(Q)
(d,v,nconv,numiter,numop,resid) = eigs(Phi,nev=1,which="LM")
# Properties: largest eigenvalue should be 1, largest eigenvector, when reshaped as matrix
# should be a Hermitian positive definite matrix (up to an arbitrary phase)

@test_approx_eq d[1] 1. # largest eigenvalue should be 1.
v=reshape(v,(50,50)) # reshape to matrix
v/=trace(v) # factor out arbitrary phase
@test isapprox(vecnorm(imag(v)),0.) # it should be real
v=real(v)
# @test isapprox(vecnorm(v-v')/2,0.) # it should be Hermitian
# Since this fails sometimes (numerical precision error),this test is commented out
v=(v+v')/2
@test isposdef(v)

# Repeat with starting vector
(d2,v2,nconv2,numiter2,numop2,resid2) = eigs(Phi,nev=1,which="LM",v0=reshape(v,(2500,)))
v2=reshape(v2,(50,50))
v2/=trace(v2)
@test numiter2<numiter
@test_approx_eq v v2

@test_approx_eq eigs(speye(50), nev=10)[1] ones(10) #Issue 4246

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia


text:
# ===================
# Inhomogeneous case
# ===================

function reach_inhomog_LGG09!(F::Vector{RT},
                              dirs::TN,
                              Ω₀::LazySet{N},
                              Φ::AbstractMatrix{N},
                              NSTEPS::Integer,
                              δ::N,
                              X::Universe,
                              U::LazySet,
                              Δt0::TimeInterval,
                              cache,
                              threaded) where {N, VN, TN, SN, RT<:TemplateReachSet{N, VN, TN, SN}}

    # transpose coefficients matrix
    Φᵀ = copy(transpose(Φ))

    # preallocate output sequence
    ρℓ = Matrix{N}(undef, length(dirs), NSTEPS)

    _reach_inhomog_dir_LGG09!(ρℓ, Ω₀, Φᵀ, U, dirs, NSTEPS, cache, threaded)

    # fill template reach-set sequence
    Δt = (zero(N) .. δ) + Δt0
    @inbounds for k in 1:NSTEPS
        F[k] = TemplateReachSet(dirs, view(ρℓ, :, k), Δt)
        Δt += δ
    end

    return ρℓ
end

function _reach_inhomog_dir_LGG09!(ρℓ, Ω₀, Φᵀ, U, dirs, NSTEPS, cache, threaded::Val{true})
    ℓ = _collect(dirs)
    Threads.@threads for j in 1:length(ℓ)
        reach_inhomog_dir_LGG09!(ρℓ, j, Ω₀, Φᵀ, U, ℓ[j], NSTEPS, cache)
    end
end

function _reach_inhomog_dir_LGG09!(ρℓ, Ω₀, Φᵀ, U, dirs, NSTEPS, cache, threaded::Val{false})
    #for each direction, compute NSTEPS iterations
    @inbounds for (j, ℓ) in enumerate(dirs)
        reach_inhomog_dir_LGG09!(ρℓ, j, Ω₀, Φᵀ, U, ℓ, NSTEPS, cache)
    end
end

function reach_inhomog_dir_LGG09!(ρvec_ℓ::AbstractMatrix{N}, j, Ω₀, Φᵀ, U, ℓ::AbstractVector{N}, NSTEPS, cache::Val{true}) where {N}
    rᵢ = copy(ℓ)
    rᵢ₊₁ = similar(rᵢ)
    sᵢ = zero(N)

    @inbounds for i in 1:NSTEPS
        ρvec_ℓ[j, i] = ρ(rᵢ, Ω₀) + sᵢ
        sᵢ += ρ(rᵢ, U)

        # update cache for the next iteration
        mul!(rᵢ₊₁, Φᵀ, rᵢ)
        copy!(rᵢ, rᵢ₊₁)   # rᵢ .= rᵢ₊₁
    end
    return ρvec_ℓ
end

function reach_inhomog_dir_LGG09!(ρvec_ℓ::AbstractMatrix{N}, j, Ω₀, Φᵀ, U, ℓ::AbstractVector{N}, NSTEPS, cache::Val{false}) where {N}
    rᵢ = copy(ℓ)
    sᵢ = zero(N)

    @inbounds for i in 1:NSTEPS
        ρvec_ℓ[j, i] = ρ(rᵢ, Ω₀) + sᵢ
        sᵢ += ρ(rᵢ, U)

        # update cache for the next iteration
        rᵢ = Φᵀ * rᵢ
    end
    return ρvec_ℓ
end

function reach_inhomog_dir_LGG09!(ρvec_ℓ::AbstractVector{N}, Ω₀, Φᵀ, U, ℓ::AbstractVector{N}, NSTEPS, cache::Val{true}) where {N}
    rᵢ = copy(ℓ)
    rᵢ₊₁ = similar(rᵢ)
    sᵢ = zero(N)

    @inbounds for i in 1:NSTEPS
        ρvec_ℓ[i] = ρ(rᵢ, Ω₀) + sᵢ
        sᵢ += ρ(rᵢ, U)

        # update cache for the next iteration
        mul!(rᵢ₊₁, Φᵀ, rᵢ)
        copy!(rᵢ, rᵢ₊₁)   # rᵢ .= rᵢ₊₁
    end
    return ρvec_ℓ
end

function reach_inhomog_dir_LGG09!(ρvec_ℓ::AbstractVector{N}, Ω₀, Φᵀ, U, ℓ::AbstractVector{N}, NSTEPS, cache::Val{false}) where {N}
    rᵢ = copy(ℓ)
    sᵢ = zero(N)

    @inbounds for i in 1:NSTEPS
        ρvec_ℓ[i] = ρ(rᵢ, Ω₀) + sᵢ
        sᵢ += ρ(rᵢ, U)

        # update cache for the next iteration
        rᵢ = Φᵀ * rᵢ
    end
    return ρvec_ℓ
end

#= ---- version using vector-of-vectors
    ρℓ = [Vector{N}(undef, NSTEPS) for _ in 1:length(dirs)]

    # for each direction, compute NSTEPS iterations
    @inbounds for (i, ℓ) in enumerate(dirs)
        reach_inhomog_dir_LGG09!(ρℓ[i], Ω₀, Φᵀ, U, ℓ, NSTEPS)
    end

    # fill template reach-set sequence
    Δt = (zero(N) .. δ) + time_shift
    @inbounds for k in 1:NSTEPS
        sf = [ρℓ[i][k] for i in 1:length(dirs)]
        F[k] = TemplateReachSet(dirs, sf, Δt)
        Δt += δ
    end
---- =#

#=
# TODO: specialization for static vector / static matrix

# compute NSTEPS iterations support function along direction ℓ
function reach_inhomog_dir_LGG09!(ρvec_ℓ, Φᵀ, ℓ, Ω₀, NSTEPS, W, time_shift)
    sᵢ = zero(T)
    rᵢ = copy(ℓ)
    rᵢ₊₁ = similar(rᵢ)

    @inbounds for i in 0:NSTEPS-2
        mul!(rᵢ₊₁, Φᵀ, rᵢ)
        sᵢ₊₁ = sᵢ + ρ(rᵢ, W)
        ρvec_ℓ[i+1] = ρ(rᵢ₊₁, Ω₀) + sᵢ₊₁

        #copy!(rᵢ, rᵢ₊₁) # update cache for the next iteration
        rᵢ .= rᵢ₊₁
        sᵢ = sᵢ₊₁
    end
    return ρvec_ℓ
end
=#

# ===================================
# Inhomogeneous case with invariant
# ===================================

# this version uses one matrix-matrix product for each time-step
function reach_inhomog_LGG09!(F::Vector{RT},
                              dirs::TN,
                              Ω₀::LazySet{N},
                              Φ::AbstractMatrix{N},
                              NSTEPS::Integer,
                              δ::N,
                              X::LazySet,
                              U::LazySet,
                              Δt0::TimeInterval,
                              cache,
                              threaded) where {N, VN, TN, SN, RT<:TemplateReachSet{N, VN, TN, SN}}

    # transpose coefficients matrix
    Φᵀ = copy(transpose(Φ))

    # preallocate output sequence
    ndirs = length(dirs)
    ρmat = Matrix{N}(undef, ndirs, NSTEPS)

    rᵢ = reduce(hcat, dirs) # one column per template direction
    rᵢ₊₁ = similar(rᵢ)
    sᵢ = zeros(N, ndirs)

    # fill template reach-set sequence
    Δt = (zero(N) .. δ) + Δt0
    k = 1
    @inbounds while k <= NSTEPS

        # loop dispatch in threaded
        _reach_inhomog_LGG09_invariant!(Ω₀, U, rᵢ, ρmat, sᵢ, k, ndirs, threaded)

        # update cache for the next iteration
        mul!(rᵢ₊₁, Φᵀ, rᵢ)
        copy!(rᵢ, rᵢ₊₁)

        F[k] = TemplateReachSet(dirs, view(ρmat, :, k), Δt)
        _is_intersection_empty(X, set(F[k])) && break  # TODO pass disjointness method
        Δt += δ
        k += 1
    end
    if k < NSTEPS
        resize!(F, k-1)
    end
    return ρmat
end

function _reach_inhomog_LGG09_invariant!(Ω₀, U, rᵢ, ρmat, sᵢ, k, ndirs, threaded::Val{false})
    for j in 1:ndirs
        d = view(rᵢ, :, j)
        ρmat[j, k] = ρ(d, Ω₀) + sᵢ[j]
        sᵢ[j] += ρ(d, U)
    end
end

function _reach_inhomog_LGG09_invariant!(Ω₀, U, rᵢ, ρmat, sᵢ, k, ndirs, threaded::Val{true})
    Threads.@threads for j in 1:ndirs
        d = view(rᵢ, :, j)
        ρmat[j, k] = ρ(d, Ω₀) + sᵢ[j]
        sᵢ[j] += ρ(d, U)
    end
end

#=
# in this version we use several matrix-vector products for each direction
function reach_inhomog_LGG09!(F::Vector{RT},
                              dirs::TN,
                              Ω₀::LazySet{N},
                              Φ::AbstractMatrix{N},
                              NSTEPS::Integer,
                              δ::N,
                              X::LazySet,
                              U::LazySet,
                              time_shift::N,
                              cache) where {N, VN, TN, SN, RT<:TemplateReachSet{N, VN, TN, SN}}

    # transpose coefficients matrix
    Φᵀ = copy(transpose(Φ))

    # preallocate output sequence
    ndirs = length(dirs)
    ρmat = Matrix{N}(undef, ndirs, NSTEPS)

    rᵢ = [copy(ℓ) for ℓ in dirs]
    rᵢ₊₁ = [similar(first(rᵢ)) for _ in 1:ndirs]
    sᵢ = zeros(N, ndirs)

    # fill template reach-set sequence
    Δt = (zero(N) .. δ) + time_shift
    k = 1
    @inbounds while k <= NSTEPS
        for (j, ℓ) in enumerate(dirs)
            ρmat[j, k] = ρ(rᵢ[j], Ω₀) + sᵢ[j]
            sᵢ[j] += ρ(rᵢ[j], U)

            # update cache for the next iteration
            mul!(rᵢ₊₁[j], Φᵀ, rᵢ[j])
            copy!(rᵢ[j], rᵢ₊₁[j])
        end

        F[k] = TemplateReachSet(dirs, view(ρmat, :, k), Δt)
        _is_intersection_empty(X, set(F[k])) && break  # TODO pass disjointness method
        Δt += δ
        k += 1
    end
    if k < NSTEPS
        resize!(F, k-1)
    end
    return ρmat
end
=#

#=
# in this version all the reach-sets are computed until NSTEPS, then we intersect
# with the invariants afterwards
function reach_inhomog_LGG09!(F::Vector{RT},
                              dirs::TN,
                              Ω₀::LazySet{N},
                              Φ::AbstractMatrix{N},
                              NSTEPS::Integer,
                              δ::N,
                              X::LazySet,
                              U::LazySet,
                              time_shift::N,
                              cache) where {N, VN, TN, SN, RT<:TemplateReachSet{N, VN, TN, SN}}

    # transpose coefficients matrix
    Φᵀ = copy(transpose(Φ))

    # preallocate output sequence
    ρℓ = Matrix{N}(undef, length(dirs), NSTEPS)

    # for each direction, compute NSTEPS iterations
    @inbounds for (j, ℓ) in enumerate(dirs)
        reach_inhomog_dir_LGG09!(ρℓ, j, Ω₀, Φᵀ, U, ℓ, NSTEPS, cache)
    end

    # fill template reach-set sequence
    Δt = (zero(N) .. δ) + time_shift
    k = 1
    @inbounds while k <= NSTEPS
        F[k] = TemplateReachSet(dirs, view(ρℓ, :, k), Δt)
        _is_intersection_empty(X, set(F[k])) && break  # TODO pass disjointness method
        Δt += δ
        k += 1
    end
    if k < NSTEPS
        resize!(F, k-1)
    end
    return ρℓ
end
=#

# ------------------------------------------------------------
# Functionality that requires ExponentialUtilities.jl
# ------------------------------------------------------------

function load_krylov_LGG09_inhomog()
return quote

# Compute the sequence with constant input sets:
#
# ρ(ℓ, Ω₀), ρ(ℓ, Φ Ω₀ ⊕ V), ρ(ℓ, Φ^2 Ω₀ ⊕ Φ V ⊕ V), ρ(ℓ, Φ^3 Ω₀ ⊕ Φ^2 V ⊕ Φ V ⊕ V), ...
#
# Using Krylov subspace approximations to compute the action of Φ := exp(Aδ) over
# the direction ℓ.
#
# Method (see [1]):
#
# out[1] <- ρ(ℓ, Ω₀)
#
# out[2] <- ρ(ℓ, Φ Ω₀ ⊕ V) = ρ(ℓ, Φ Ω0) + ρ(ℓ, V) = ρ(Φᵀ ℓ, Ω₀) + ρ(ℓ, V)
#
# out[3] <- ρ(ℓ, Φ^2 Ω₀ ⊕ Φ V ⊕ V) = ρ((Φᵀ)^2 ℓ, Ω₀) + ρ(Φᵀ ℓ, V) + ρ(ℓ, V)
#
# out[4] <- ρ(ℓ, Φ^3 Ω₀ ⊕ Φ^2 V ⊕ Φ V ⊕ V) = ρ((Φᵀ)^3 ℓ, Ω₀) + ρ((Φᵀ)^2 ℓ, V) + ρ(Φᵀ ℓ, V) + ρ(ℓ, V)
#
# and so on.
#
# [1] Reach Set Approximation through Decomposition with Low-dimensional Sets and
#     High-dimensional Matrices. Sergiy Bogomolov, Marcelo Forets, Goran Frehse,
#     Frédéric Viry, Andreas Podelski and Christian Schilling (2018) HSCC'18
#     Proceedings of the 21st International Conference on Hybrid Systems: Computation
#     and Control: 41–50.
#
function reach_inhomog_krylov_LGG09!(out, Ω₀::LazySet, V::LazySet, Aᵀδ::AbstractMatrix,
                                     ℓ::AbstractVector, NSTEPS;
                                     hermitian=false, m=min(30, size(Aᵀδ, 1)), tol=1e-7)

    # initialization of the krylov subspace
    TA, Tb = eltype(Aᵀδ), eltype(ℓ)
    T = promote_type(TA, Tb)
    Ks = KrylovSubspace{T, real(T)}(length(ℓ), m)
    arnoldi!(Ks, Aᵀδ, ℓ; m=m, ishermitian=hermitian, tol=tol)

    # rᵢ stores is the cache for each vector: (Φᵀ)^i ℓ
    rᵢ = deepcopy(ℓ)
    out[1] = ρ(ℓ, Ω₀)

    # accumulated support vector sum due to the inputs
    s = zero(T)

    @inbounds for i in 1:NSTEPS-1
        s += ρ(rᵢ, V)
        expv!(rᵢ, i*1.0, Ks)
        out[i+1] = ρ(rᵢ, Ω₀) + s
    end
    return out
end

# Compute the sequence with time-varying input sets:
#
# ρ(ℓ, Ω₀), ρ(ℓ, Φ Ω₀ ⊕ V₀), ρ(ℓ, Φ^2 Ω₀ ⊕ Φ V₀ ⊕ V₁), ρ(ℓ, Φ^3 Ω₀ ⊕ Φ^2 V₀ ⊕ Φ V₁ ⊕ V₂), ...
#
# Using Krylov subspace approximations to compute the action of Φ := exp(Aδ) over
# the direction ℓ.
#
# Method (see [1]):
#
# out[1] <- ρ(ℓ, Ω₀)
#
# out[2] <- ρ(ℓ, Φ Ω₀ ⊕ V₀) = ρ(ℓ, Φ Ω0) + ρ(ℓ, V₀) = ρ(Φᵀ ℓ, Ω₀) + ρ(ℓ, V₀)
#
# out[3] <- ρ(ℓ, Φ^2 Ω₀ ⊕ Φ V₀ ⊕ V₁) = ρ((Φᵀ)^2 ℓ, Ω₀) + ρ(Φᵀ ℓ, V₀) + ρ(ℓ, V₁)
#
# out[4] <- ρ(ℓ, Φ^3 Ω₀ ⊕ Φ^2 V₀ ⊕ Φ V₁ ⊕ V₂) = ρ((Φᵀ)^3 ℓ, Ω₀) + ρ((Φᵀ)^2 ℓ, V₀) + ρ(Φᵀ ℓ, V₁) + ρ(ℓ, V₂)
#
# and so on.
#
# [1] Reach Set Approximation through Decomposition with Low-dimensional Sets and
#     High-dimensional Matrices. Sergiy Bogomolov, Marcelo Forets, Goran Frehse,
#     Frédéric Viry, Andreas Podelski and Christian Schilling (2018) HSCC'18
#     Proceedings of the 21st International Conference on Hybrid Systems: Computation
#     and Control: 41–50.
#
function reach_inhomog_krylov_LGG09!(out, Ω₀::LazySet, V::Vector{<:LazySet}, Aᵀδ::AbstractMatrix,
                                     ℓ::AbstractVector, NSTEPS;
                                     hermitian=false, m=min(30, size(Aᵀδ, 1)), tol=1e-7)

    # initialization of the krylov subspace
    TA, Tb = eltype(Aᵀδ), eltype(ℓ)
    T = promote_type(TA, Tb)
    Ks = KrylovSubspace{T, real(T)}(length(ℓ), m)
    arnoldi!(Ks, Aᵀδ, ℓ; m=m, ishermitian=hermitian, tol=tol)

    # rᵢ stores is the cache for each vector: (Φᵀ)^i ℓ
    r = [similar(ℓ) for _ in NSTEPS]
    r[1] = deepcopy(ℓ)
    out[1] = ρ(ℓ, Ω₀)

    @inbounds for i in 1:NSTEPS-1
        s = zero(T)
        for j in 1:i
            s += ρ(r[i-j+1], V[j])
        end
        expv!(r[i+1], i*1.0, Ks)
        out[i+1] = ρ(r[i+1], Ω₀) + s
    end
    return out
end

end end  # quote / load_krylov_LGG09_inhomog()

# ------------------------------------------------------------
# Methods using eigenvalues of the transition matrix
# ------------------------------------------------------------

#= FIXME needs review <<

# it is assumed that (λ, d) is an eigenvalue-eigenvector pair of the matrix Φᵀ
function reach_inhomog_dir_eig_LGG09!(out::AbstractVector{N}, X₀, V, d::AbstractVector{N}, λ::N, NSTEPS) where {N}
    if iszero(λ)
        _reach_inhomog_dir_eig_LGG09_zero!(out, V, NSTEPS)
    elseif λ > zero(N)
        _reach_inhomog_dir_eig_LGG09_positive!(out, X₀, V, d, λ, NSTEPS)
    else
        _reach_inhomog_dir_eig_LGG09_negative!(out, X₀, V, d, λ, NSTEPS)
    end
    return out
end

function _reach_inhomog_dir_eig_LGG09_zero!(out::AbstractVector{N}, V::LazySet, NSTEPS) where {N}
    ρ_d_V = ρ(d, V)
    @inbounds for i in 1:NSTEPS
        out[i] = ρ_d_V
    end
end

function _reach_inhomog_dir_eig_LGG09_zero!(out::AbstractVector{N}, V::AbstractVector{<:LazySet}, NSTEPS) where {N}
    @inbounds for i in 1:NSTEPS
        out[i] = ρ(d, V[i])
    end
end

function _reach_inhomog_dir_eig_LGG09_positive!(out::AbstractVector{N}, X₀, V::LazySet, d, λ, NSTEPS) where {N}
    ρ_d_X₀ = ρ(d, X₀)
    ρ_d_V = ρ(d, V)
    @inbounds begin
        out[1] = ρ_d_X₀
        λⁱ = λ
        λacc = one(N)
        for i in 2:NSTEPS
            out[i] = λⁱ * ρ_d_X₀ + λacc * ρ_d_V
            λacc = λacc + λⁱ
            λⁱ = λⁱ * λ
        end
    end
end

function _reach_inhomog_dir_eig_LGG09_positive!(out::AbstractVector{N}, X₀, V::AbstractVector{<:LazySet}, d, λ, NSTEPS) where {N}
    ρ_d_X₀ = ρ(d, X₀)
    @inbounds begin
        out[1] = ρ_d_X₀
        ρ_d_V_acc = ρ(d, V[1])
        λⁱ = λ
        for i in 2:NSTEPS
            out[i] = λⁱ * ρ_d_X₀ + ρ_d_V_acc
            ρ_d_V_acc = λ * ρ_d_V_acc + ρ(d, V[i])
            λⁱ = λⁱ * λ
        end
    end
end

function _reach_inhomog_dir_eig_LGG09_negative!(out::AbstractVector{N}, X₀, V::LazySet, d, λ, NSTEPS) where {N}
    ρ_d_X₀ = ρ(d, X₀)
    ρ_md_X₀ = ρ(-d, X₀)
    ρ_d_V = ρ(d, V)
    ρ_md_V = ρ(-d, V)
    λⁱ = λ
    @inbounds begin
        out[1] = ρ_d_X₀
        for i in 2:NSTEPS
            if iseven(i)
                out[i] = -λⁱ * ρ_md_X₀
            else
                out[i] = λⁱ * ρ_d_X₀
            end
            λⁱ = λⁱ * λ
        end
    end
end

function _reach_inhomog_dir_eig_LGG09_negative!(out::AbstractVector{N}, X₀, V::AbstractVector{<:LazySet}, d, λ, NSTEPS) where {N}
    ρ_d_X₀ = ρ(d, X₀)
    ρ_md_X₀ = ρ(-d, X₀)
    λⁱ = λ
    @inbounds begin
        out[1] = ρ_d_X₀
        for i in 2:NSTEPS
            if iseven(i)
                out[i] = -λⁱ * ρ_md_X₀
            else
                out[i] = λⁱ * ρ_d_X₀
            end
            λⁱ = λⁱ * λ
        end
    end
end
=#


text:
module GenOneDTMazeDpi

using Curiosity
using Curiosity: DPI, HordeDPI
using Curiosity: OpenWorld
using MinimalRLCore
include("../experiments/1d-tmaze.jl")
const ODTMU = Curiosity.OneDTMazeUtils
using Random
using StatsBase
using GVFHordes
using Statistics
using JLD2

using ProgressMeter

StatsBase.sample(p::GVFHordes.GVFParamFuncs.FunctionalPolicy, s, actions) =
    sample(Weights([p.func(;state_t = s, action_t = a) for a in actions]))

StatsBase.sample(rng, p::GVFHordes.GVFParamFuncs.AbstractPolicy, s, actions) =
    sample(Weights([get(p;state_t = s, action_t = a) for a in actions]))

function gen_dpi(num_states=100000, seed=1; exploring_starts="beg", h_kwargs...)
    Random.seed!(seed)

    parsed = OneDTmazeExperiment.default_args()
    parsed["exploring_starts"] = exploring_starts
    parsed["demon_learner"] = "Q"

    dpis = DPI[]

    state_filter_func =  Curiosity.FeatureSubset(identity, 1:2)

    for gvf_i ∈ 1:4
        cumulant_schedule = ODTMU.get_cumulant_schedule(parsed)
        env = OneDTMaze(cumulant_schedule, parsed["exploring_starts"], parsed["env_step_penalty"])

        ss = Vector{Float64}[]
        as = Int[]
        policy = ODTMU.GoalPolicy(gvf_i)

        while length(ss) < num_states
            s = MinimalRLCore.start!(env)
            push!(ss, state_filter_func(s))

            a = policy(s)
            term = false
            while term == false
                s, r, term = MinimalRLCore.step!(env, a)
                a = policy(s)

                push!(ss, state_filter_func(s))
                if term
                    break
                end
            end
        end

        push!(dpis, DPI(ss, policy; h_kwargs...))
    end

    return HordeDPI(dpis, state_filter_func)
end

function save_dpi(num_states, seed; exploring_starts="beg", h_kwargs...)
    hdpi = gen_dpi(num_states, seed; exploring_starts=exploring_starts, h_kwargs...)
    @save "./src/data/dpi/oned_tmaze_$(exploring_starts).jld2" hdpi
end

end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
"""
```
default_settings!(m::AbstractModel)
```

Default Settings are constructed, initialized and added to `m.settings`.
"""
function default_settings!(m::AbstractModel)

    # I/O File locations
    saveroot = normpath(joinpath(dirname(@__FILE__), "..","save"))
    datapath = normpath(joinpath(dirname(@__FILE__), "..","save","input_data"))

    m <= Setting(:saveroot, saveroot, "Root of data directory structure")
    m <= Setting(:dataroot, datapath, "Input data directory path")

    # Data settings for released and conditional data. Default behavior is to set vintage
    # of data to today's date.
    vint = Dates.format(now(), DSGE_DATE_FORMAT)
    m <= Setting(:data_vintage, vint, true, "vint", "Vintage of data")
    m <= Setting(:cond_vintage, vint, "Vintage of conditional data")
    m <= Setting(:cond_id, "0000", "Identifier of conditional dataset")
    m <= Setting(:use_population_forecast, false, "Whether to use population forecasts as data")
    m <= Setting(:population_mnemonic, :CNP16OV, "Mnemonic of FRED data series for computing per-capita values")

    # Dates
    m <= Setting(:date_presample_start, quartertodate("1959-Q3"), "Start date of pre-sample")
    m <= Setting(:date_mainsample_start, quartertodate("1960-Q1"), "Start date of main sample")
    m <= Setting(:date_zlbregime_start, quartertodate("2008-Q4"), "Start date of zero lower bound regime")
    m <= Setting(:date_mainsample_end, Dates.lastdayofquarter(Dates.today()-Dates.Month(3)), "End date of main sample")
    m <= Setting(:date_forecast_start, Dates.lastdayofquarter(Dates.today()), "Start date of forecast period")
    m <= Setting(:date_forecast_end, Dates.lastdayofquarter(Dates.today()+Dates.Month(60*3)), "End date of forecast period")

    # Anticipated shocks
    m <= Setting(:n_anticipated_shocks, 0, "Number of anticipated policy shocks")
    m <= Setting(:n_anticipated_shocks_padding, 20, "Padding for anticipated policy shocks")

    # General computation
    m <= Setting(:use_parallel_workers, true, "Use available parallel workers in computations")

    # Estimation
    m <= Setting(:reoptimize, true, "Optimize the posterior mode. If false, reads in mode from a file.")
    m <= Setting(:calculate_hessian, true, "Calculate the hessian at the mode")
    m <= Setting(:n_hessian_test_params, typemax(Int), "Max number of free params for which to calculate Hessian")

    # Metropolis-Hastings
    m <= Setting(:n_mh_simulations, 5000, "Number of draws per block in Metropolis-Hastings")
    m <= Setting(:n_mh_blocks, 22, "Number of blocks for Metropolis-Hastings")
    m <= Setting(:n_mh_burn, 2, "Number of blocks to use as burn-in in Metropolis-Hastings")
    m <= Setting(:mh_thin, 5, "Metropolis-Hastings thinning step")

end


text:
import StatsFuns.gammacdf

seed!(12345)

function testGamma(α, β, low, high, h, N)
  function cdf(x::Float64)
    return gammacdf(α, 1/β, x)
  end
  function sampler()
    return sampleGamma(α, β)
  end
  return testGOFContinuous(cdf, sampler, low, high, h, N)
end

@test testGamma(3.0, 2.0, 0.01, 6.0, 0.01, 2^21) > 0.01
@test testGamma(0.1, 1.0, 0.00, 4.0, 0.01, 2^21) > 0.01
@test testGamma(0.8, 1.0, 0.00, 7.0, 0.01, 2^21) > 0.01
@test testGamma(10.0, 1.0, 3.0, 23.0, 0.01, 2^21) > 0.01


text:
function etree{Tv,Ti}(A::Symmetric{Tv,SparseMatrixCSC{Tv,Ti}},postorder::Bool=false)
    etree(A.uplo == 'U' ? A.data : A.data', postorder)
end

function Base.symperm{Tv,Ti}(A::Symmetric{Tv,SparseMatrixCSC{Tv,Ti}},invperm)
    Symmetric(symperm(A.uplo == 'U' ? A.data : A.data',invperm),:U)
end

function Metis.metisform{Tv,Ti}(S::Symmetric{Tv,SparseMatrixCSC{Tv,Ti}})
    up = S.uplo == 'U'
    Sd = S.data
    n = size(Sd,2)
    I = sizehint!(Cint[],2*nnz(Sd))
    J = sizehint!(Cint[],2*nnz(Sd))
    Si = rowvals(Sd)
    for j in 1:n, i in sub(Si,nzrange(Sd,j))
        if (up ? (i < j) : (i > j))
            push!(I,i)
            push!(J,j)
            push!(I,j)
            push!(J,i)
        end
    end
    sp = sparse(I,J,1.)
    convert(Int32,n),sp.colptr .- one(Int32),rowvals(sp) .- one(Int32)
end


text:
import Base.isequal
import Base.hash
import Base.show

"""
abstract type ASTNumNode
    represents the AST nodes which evaluate to a number, when supplied with
    all the data (mds, fluxes)
"""
abstract type ASTNumNode end

struct ASTMatrix
    m::Array{ASTNumNode,2}
end

"""
ASTSymmOp
has fields: op::Char , args::Vector{ASTNumNode}
represents addition and multiplication operations
"""
struct ASTSymmOp <: ASTNumNode
    op::Char
    args::Vector{ASTNumNode}
end

"""
ast_add(a::ASTNumNode,b::ASTNumNode)
creates sum of two ASTNumNode objects
"""
function ast_add(a::ASTNumNode,b::ASTNumNode)
    return ASTSymmOp('+',[a,b])
end

function ast_add(a::ASTSymmOp,b::ASTNumNode)
    if(a.op=='+')
        push!( a.args , b )
        return a
    else
        return ASTSymmOp('+',[a,b])
    end
end
function ast_add(a::ASTNumNode,b::ASTSymmOp)
    ast_add(b,a)
end

"""
ASTRealValue has field val::Float64
and represents a hard-coded real value
"""
struct ASTRealValue <: ASTNumNode
    val::Float64
end

"""
ASTEMUMassRatio has fields
  semu::SpeciesEMU , mass:Float64

and represents a mass ratio (which has to be supplied to compute the FW-Sim)
"""
struct ASTEMUMassRatio <: ASTNumNode
    semu::SpeciesEMU
    mass::Float64 # here we could also have int, basically..
end
function Base.hash(emr::ASTEMUMassRatio)
    mass = round(emr.mass)
    return xor( Base.hash(emr.semu) , Base.hash(mass) )
end
function Base.isequal(a::ASTEMUMassRatio,b::ASTEMUMassRatio)
    return Base.isequal(a.semu,b.semu) && isequal(round(a.mass),round(b.mass))
end
function Base.isless(a::ASTEMUMassRatio,b::ASTEMUMassRatio)
    if( isless( a.semu , b.semu ) )
        return true
    end
    if( isequal(a.semu,b.semu) )
        return round(a.mass) < round( b.mass )
    end
    return false
end
function Base.show(io::IO,emr::ASTEMUMassRatio)
    print(io,"ASTEMR:$(emr.semu):M=$(round(emr.mass))")
end

"""
ASTFluxRate has fields
  flux::String

and represents a flux rate (which has to be supplied to compute the FW-Sim)
"""
struct ASTFluxRate <: ASTNumNode
    flux::String
end
function Base.hash(f::ASTFluxRate)
    return Base.hash("flux$(f.flux)")
end
function Base.isequal(a::ASTFluxRate,b::ASTFluxRate)
    return isequal(a.flux,b.flux)
end
function Base.isless(a::ASTFluxRate,b::ASTFluxRate)
    return a.flux < b.flux
end

"""
ASTEMUSimStep has fields
  A::ASTMatrix , X::ASSTMatrix , B::ASTMatrix , Y::ASSTMatrix

and represents a step of the FW-Sim. Notation is according to Antoniewicz2006,
i.e. new SEMU masses are in matrix X, and given ones are in matrix Y,
i.e we have to compute the step X = inv(A) * B * Y
"""
struct ASTEMUSimStep
    A::ASTMatrix
    X::ASTMatrix
    B::ASTMatrix
    Y::ASTMatrix
end

"""
ASTEMUSim has fields
  sim_steps::Dict{Integer,Vector{ASTEMUSimStep}}
  configured_semus::Vector{SpeciesEMU}

  and represents all information required to compute the FW-Sim, for given fluxes
  and mds for the semus in .configured_semus
  the simulation steps are sorted by size of EMUs
"""
struct ASTEMUSim
    sim_steps::Dict{Integer,Vector{ASTEMUSimStep}}
    configured_semus::Vector{SpeciesEMU}
end

"""
create_ast_from_emu_graph( eg::Dict{Integer,Vector{AVGraph}} )
returns ASTEMUSim

computes the AST which represents the forward simulation
"""
function create_ast_from_emu_graph( eg::Dict{Integer,Vector{AVGraph}} )
    required_semus::Set{SpeciesEMU} = Set{SpeciesEMU}()
    determined_semus::Set{SpeciesEMU} = Set{SpeciesEMU}()
    for si in sort( collect( keys(eg) ) ) # we have to go from small to large, as in the emu decomp..
        for gi in eg[si]
            ( roots::Vector{SpeciesEMU} , non_roots::Vector{SpeciesEMU} ) = get_roots_of_av_graph(gi)
            req_new = setdiff( Set( roots ) , determined_semus )
            union!( required_semus , req_new )
            union!( determined_semus , Set(non_roots ) )
        end
    end
    print("\nidentified $(length(required_semus)) required semus:\n")
    for rsi in required_semus; print("  $(rsi)\n"); end
    # ok, unroll them into single sas..
    required_sas = reduce( (x,y) -> union(x,y) , required_semus )
    #required_sas = sort(required_sas)
    print("\nSAS Set: $(length(required_sas)) required sas:\n")
    for sasi in sort(collect(required_sas))
        print(" $(sasi) ")
    end
    print("\n")

    return create_ast_from_emu_graph( eg , required_sas)
end

# configured_sas : all SAS which are given by the substrate configuration
function create_ast_from_emu_graph( eg::Dict{Integer,Vector{AVGraph}} , configured_sas::Set{SAS})

    debug_level = 1;

    sizes = sort( collect( keys( eg ) ) )

    # for sanity we keep track of all determined semu on the way..
    all_determined_semus::Set{SpeciesEMU} = Set{SpeciesEMU}()

    ast_emu_sim = ASTEMUSim(Dict{Integer,Vector{ASTEMUSimStep}}() , Vector{SpeciesEMU}())

    for si in sizes
        ast_emu_sim.sim_steps[si] = Vector{ASTEMUSimStep}()
        ni_cnt = 0
        for ni in eg[si]
            ni_cnt += 1
            # ok, lets build the matrix..
            # Step 1: partition semus into determined / target
            # Explanation: determined are:
            #    1. all non-max-size EMU ,
            #    2. EMUs containing only configured SAS
            # OR: easier: undetermined EMUs are all the product EMUs
            s_det::Set{SpeciesEMU}    = Set{SpeciesEMU}()
            s_nondet::Set{SpeciesEMU} = Set{SpeciesEMU}()

            emu_rxns  = get_all_emurxns(ni)
            semu_rxns = map( x1 -> SpeciesEMURxn(x1) , emu_rxns )
            semus     = get_all_semus(ni)

            for sri in semus
                if( length(sri) < si )
                    push!( s_det , sri )
                else
                    if( issubset( sri , configured_sas ) )
                        push!( s_det , sri )
                    else
                        push!( s_nondet , sri )
                    end
                end
            end

            if(debug_level>=1)
                print(  "--------------------------------------------------------------\n")
                print(  "NETWORK OF SIZE: $(si) , $(ni_cnt)/$(length(eg[si])) \n")
                print(  "--------------------------------------------------------------\n")
                print(  "ALL EMU RXNs: ------------------------------------------------\n")
                for rxi in semu_rxns; print("  $(rxi)\n"); end
                print(  "SEMUs (determined): ------------------------------------------\n")
                #for si in semus; print("  $(si)\n"); end
                for sid in s_det; print("  $(sid)\n"); end
                print(  "SEMUs (non-determined): --------------------------------------\n")
                #for si in semus; print("  $(si)\n"); end
                for sind in s_nondet; print("  $(sind)\n"); end
                print(  "--------------------------------------------------------------\n\n")
            end

            ## NOW: Build the matrices:
            #  X : Rows correspond to the non-det EMUs
            #      Cols correspond to the weight
            #  Y : Rows correspond to the determined EMUs
            #      Cols correspond to the weight
            #  A : Rows correspond to the non-det EMUs
            #      Cols correspond to the non-det EMUs
            #  B : Rows correspond to the non-det EMUs
            #      Cols correspond to the det EMUs

            # Order of rows / columns determined by:
            nondet_semu_order = sort(collect(s_nondet))
            det_semu_order    = sort(collect(s_det))

            ndspos = (x1::SpeciesEMU) -> find( xi -> xi==x1 , nondet_semu_order )
            #dspos  = (x1::SpeciesEMU) -> find( xi -> xi==x1 , det_semu_order )

            # HOW TO BUILD:
            # For each EMU we consider the influxes.
            # The influxes times the EMU mixture has to be equal
            # to the combined influxes of the EMU mixtures which feed the EMU

            # Build A, that's easy;
            A::Array{ASTNumNode,2} = Array{ASTNumNode,2}(length(nondet_semu_order),length(nondet_semu_order))
            for ns in nondet_semu_order
                #build this row:

                # find all emu rxns with product ns
                emu_rxns_ns = filter( er -> isequal( pemu2semu(er.product_emu) , ns ) , emu_rxns )
                # create the "entry on the diagonal"
                diag_entry::ASTNumNode = ASTSymmOp('+',Vector{ASTNumNode}())
                for eri in emu_rxns_ns
                    push!( diag_entry.args , ASTSymmOp('*',[ASTRealValue(-1),ASTFluxRate(eri.flux)] ) )
                end
                A[ndspos(ns),ndspos(ns)] = diag_entry
                # create the other entries:
                for eri in emu_rxns_ns
                    # skip for emu_rxns with more than one remu (becuase those cannot have undetermined remus..)
                    if( length( eri.reactand_emus ) > 1 ); continue; end
                    rpos  = ndspos( remu2semu( eri.reactand_emus[1] ) )
                    if(!isempty(rpos))
                        ndpos = ndspos(ns) ; @assert length(ndpos)==1 ; ndpos = ndpos[1] ;
                        @assert length(rpos)==1 ; rpos = rpos[1] ;
                        if( ! isassigned(A,ndpos,rpos))
                            A[ndpos,rpos] = ASTFluxRate(eri.flux)
                        else
                            A[ndpos,rpos] = ast_add( A[ndpos,rpos] , ASTFluxRate(eri.flux) )
                        end
                    end
                end
            end
            if(debug_level>=2) ; print(A) ; end

            # Build X (trivial)
            #
            if(debug_level>=2); println("\n\nsi : $(si)\n\n") ; end
            X::Array{ASTNumNode,2} = Array{ASTNumNode,2}(length(nondet_semu_order),si+1)
            for zxi in 0:si
                for ns in nondet_semu_order
                    X[ ndspos(ns) ,zxi+1] = ASTEMUMassRatio(ns,zxi)
                end
            end

            # Build B, similar to A, just with cols are the emu rxns with det.
            # reactand emus. NOTE: if one reactand is det. then all are det. ( as all are non-maximal size)

            det_rxns_1 = filter( x1 -> isempty(setdiff( x1.reactand_emus , s_det )) , semu_rxns )
            det_rxns_2 = filter( x1 -> !isempty(intersect( x1.reactand_emus , s_det )) , semu_rxns )
            @assert length(det_rxns_1)==length(det_rxns_2)

            det_rxns      = det_rxns_1
            det_rxn_order = sort(det_rxns) # this order defines the rows of B and the cols of Y
            drpos         = (x1::SpeciesEMURxn) -> find( xi->xi==x1 , det_rxn_order )[1]

            B::Array{ASTNumNode,2} = Array{ASTNumNode,2}(length(nondet_semu_order),length(det_rxns))
            Y::Array{ASTNumNode,2} = Array{ASTNumNode,2}(length(det_rxns),si+1)
            # Build B
            for ns in nondet_semu_order
                for rd in det_rxn_order
                    if( isequal( rd.product_emu , ns ) )
                        # !! MINUS
                        B[ndspos(ns),drpos(rd)] = ASTSymmOp('*',[ASTRealValue(-1),ASTFluxRate(rd.flux)] ) #ASTFluxRate(rd.flux)
                    end
                end
            end
            # Build Y , colums are masses 0:si
            for rd in det_rxn_order
                for mi in 0:si
                    # now collect all combinations of masses to get mi
                    rd_semus = rd.reactand_emus
                    rd_semus_sizes = map( rxi -> length(rxi) ,  rd_semus)
                    mmcs = matching_mass_combinations( mi , rd_semus_sizes)
                    yi::ASTNumNode = ASTSymmOp('+',Vector{ASTNumNode}(0))
                    for mmcsi in mmcs
                        yi_p::ASTNumNode = ASTSymmOp('*',Vector{ASTNumNode}(0))
                        wmicnt = 0
                        for wmi in mmcsi
                            wmicnt += 1
                            push!( yi_p.args , ASTEMUMassRatio( rd_semus[wmicnt] , wmi ) )
                        end
                        push!(yi.args,yi_p)
                    end
                    Y[drpos(rd),mi+1] = yi
                end
            end

            emusimstep_i = ASTEMUSimStep(ASTMatrix(A),ASTMatrix(X),ASTMatrix(B),ASTMatrix(Y))
            push!( ast_emu_sim.sim_steps[si] , emusimstep_i )
        end
    end

    if(debug_level>=1) ; println("\ndone..\n") ; end
    return ast_emu_sim
end

function construct_all_mass_ratios_for_measurement(fd::FluxiData , measurement::String )
    construct_all_mass_ratios_for_measurement(fd.measurements[measurement])
end

function construct_all_mass_ratios_for_measurement(vmd::Vector{MassDistributionMeasurement})
    smrs::Set{ASTEMUMassRatio} = Set{ASTEMUMassRatio}()
    for mdi in vmd
        #print(mdi)
        semu = mdi.sass
        for mi=0:length(semu)
            push!(smrs,ASTEMUMassRatio(semu,mi))
        end
    end
    return sort(collect(smrs))
end

function construct_all_mass_ratios_for_measurement(semus::Vector{SpeciesEMU})
    smrs::Set{ASTEMUMassRatio} = Set{ASTEMUMassRatio}()
    for si in semus
        for mi=0:length(si)
            push!(smrs,ASTEMUMassRatio(si,mi))
        end
    end
    return sort(collect(smrs))
end

function construct_mds(fd::FluxiData,substrate_conf::String,mds::Vector{ASTEMUMassRatio})
    construct_mds( fd.substrate_conf[substrate_conf] , mds )
end

"""
construct_mds(vmd::Vector{ConfiguredSEMUMixture},mds::Vector{ASTEMUMassRatio})
returns Vector{Float64}
computes the mds which are required to run the fw-sim from a ConfiguredSEMUMixture
object.

# Note: Required to init the fw simulation
"""
function construct_mds(vmd::Vector{ConfiguredSEMUMixture},mds::Vector{ASTEMUMassRatio})
    md_values = Vector{Float64}(length(mds))
    # preprocess the vmd..
    d_cs = Dict{ConfiguredSEMU,Float64}()
    for vmdi in vmd; d_cs = merge(d_cs,vmdi); end

    cnt_md = 0
    for md in mds
        cnt_md +=1
        @assert ( isinsinglespecies( md.semu ) )
        atom_positions = sort( map( s1 -> s1.pos , collect(md.semu) ) )

        # find the entries in d_cs which cover this semu..
        csemus = filter( x1 -> isempty( setdiff( md.semu , x1.semu ) ) , collect(keys(d_cs)) )
        #println("\nTotal values: $(sum( map( x1 -> d_cs[x1] , csemus ) ))\n")
        @assert sum( map( x1 -> d_cs[x1] , csemus ) ) > 0.999
        @assert sum( map( x1 -> d_cs[x1] , csemus ) ) < 1.001

        # now just select all which have md.mass labeled atoms and sum up..
        m_semus = filter( x1 -> round( sum(x1.config[atom_positions]))==round(md.mass) , csemus )
        if(isempty(m_semus))
            md_values[cnt_md] = 0.0;
        else
            md_values[cnt_md] = sum( map( x -> d_cs[x] , m_semus ) )
        end
    end
    return md_values
end

function get_mass_distribution_measurement_value( fd::FluxiData , mname::String , emu::SpeciesEMU , mass::Float64 )
    return get_mass_distribution_measurement_value( fd , mname , ASTEMUMassRatio(emu,mass) )
end

function get_mass_distribution_measurement_value( fd::FluxiData , mname::String , edr::ASTEMUMassRatio )
    mds::Array{MassDistributionMeasurement,1} = fd.measurements[mname]
    # 1. find the right measurement:
    mdsi = filter( x -> x.sass == edr.semu , mds)
    if(length(mdsi)==0);
        println("\nNo measurement found for semu $(edr.semu)")
        return NaN;
    end
    if(length(mdsi)>1)
        println("\n!!MULTIPLE MEASUREMENTS FOUND!! for semu $(edr.semu)")
        data = Array{Float64,2}(3,length(mdsi))
        for zi in 1:length(mdsi)
            data[:,zi] = mdsi[zi].data[edr.mass]
        end
    else
        data = Array{Float64,2}(3,1)
        data[:,1] = mdsi[1].data[edr.mass]
    end
    return data
end

# AST iterator:
function ast_traverse_df!( n::ASTEMUSim , v::Vector{Any})
    for si in sort(collect(keys(n.sim_steps)))
        for ssi in n.sim_steps[si]
            ast_traverse_df!( ssi.A , v )
            ast_traverse_df!( ssi.X , v )
            ast_traverse_df!( ssi.B , v )
            ast_traverse_df!( ssi.Y , v )
            push!( v , ssi )
        end
    end
    push!( v , n )
end
function ast_traverse_df!( n::ASTMatrix , v::Vector{Any})
    for zi=1:size(n.m,2)
        for zj=1:size(n.m,1)
            if(isassigned(n.m,zj,zi))
                ast_traverse_df!( n.m[zj,zi] , v )
            end
        end
    end
    push!( v , n )
end
function ast_traverse_df!( n::ASTSymmOp , v::Vector{Any})
    for ai in n.args
        ast_traverse_df!( ai , v )
    end
    push!( v , n )
end
function ast_traverse_df!( n::ASTNumNode , v::Vector{Any})
    push!( v , n )
end


text:
# """
# Turn sheet in a two-sided mesh, i.e. a mesh that contains two copies of the
# same geometric entity, one part of the front of the mesh, the other part of
# the back.
# """
# function twosided(edges::AbstractMesh, faces::AbstractMesh)

#     @assert dimension(edges) == 1
#     @assert vertices(edges) == vertices(faces)

#     back = -deepcopy(edges)
#     bnd = boundary(faces)
#     int = submesh(!in(bnd), back)

#     @show length(int)

#     return union(edges, int)
# end

function twosided(sheet::Mesh)

    Mesh(
        sheet.vertices,
        [sheet.faces; flip.(sheet.faces)])
end

function skeleton_twosided(mesh::AbstractMesh, dim::Int)

    @assert dimension(mesh) == 2
    @assert dim == 1 "only construction of edges supported for now"
    meshdim = dimension(mesh)
    @assert 0 <= dim <= meshdim

    if dim == meshdim
        return mesh
    end

    nc = length(mesh)
    T = SVector{dim+1,Int}
    simplices = zeros(T, nc*binomial(meshdim+1,dim+1))

    n = 1
    for (c,cell) in enumerate(mesh)
        for k in 1:3
            simplex = T(cell[k], cell[mod1(k+1,3)])
            if simplex == sort(simplex)
                simplices[n] = simplex
                n += 1
            end
        end
    end

    resize!(simplices,n-1)

    for simplex in simplices
        pos = findall(==(simplex), simplices)
        length(pos) == 2 || continue
        simplices[pos[2]] = flip(simplices[pos[2]])
    end

    Mesh(vertices(mesh), simplices)

end


text:
import .ExpressionExplorer: UsingsImports

"A container for the result of parsing the cell code, with some extra metadata."
Base.@kwdef struct ExprAnalysisCache
    code::String=""
    parsedcode::Expr=Expr(:toplevel, LineNumberNode(1), Expr(:block))
	module_usings_imports::UsingsImports = UsingsImports()
    function_wrapped::Bool=false
end

ExprAnalysisCache(notebook, cell::Cell) = let
    parsedcode=parse_custom(notebook, cell)
    ExprAnalysisCache(
        code=cell.code,
        parsedcode=parsedcode,
        module_usings_imports=ExpressionExplorer.compute_usings_imports(parsedcode),
        function_wrapped=ExpressionExplorer.can_be_function_wrapped(parsedcode),
    )
end

struct DefaultDict{K,V} <: AbstractDict{K,V}
    default::Union{Function,DataType}
    container::Dict{K,V}
end

"The (information needed to create the) dependency graph of a notebook. Cells are linked by the names of globals that they define and reference. 🕸"
Base.@kwdef struct NotebookTopology
    nodes::DefaultDict{Cell,ReactiveNode} = DefaultDict{Cell,ReactiveNode}(ReactiveNode)
    codes::DefaultDict{Cell,ExprAnalysisCache}=DefaultDict{Cell,ExprAnalysisCache}(ExprAnalysisCache)
end

DefaultDict{K,V}(default::Union{Function,DataType}) where {K,V} = DefaultDict{K,V}(default, Dict{K,V}())

function Base.getindex(aid::DefaultDict{K,V}, key::K)::V where {K,V}
    get!(aid.default, aid.container, key)
end
function Base.merge(a1::DefaultDict{K,V}, a2::DefaultDict{K,V}) where {K,V}
    DefaultDict{K,V}(a1.default, merge(a1.container, a2.container))
end
function Base.merge(a1::DefaultDict{K,V}, a2::AbstractDict) where {K,V}
    DefaultDict{K,V}(a1.default, merge(a1.container, a2))
end
Base.setindex!(aid::DefaultDict, args...) = Base.setindex!(aid.container, args...)
Base.delete!(aid::DefaultDict, args...) = Base.delete!(aid.container, args...)
Base.keys(aid::DefaultDict) = Base.keys(aid.container)
Base.values(aid::DefaultDict) = Base.values(aid.container)
Base.iterate(aid::DefaultDict, args...) = Base.iterate(aid.container, args...)


text:
using RealInterface
using Base.Test
import SpecialFunctions, NaNMath

for f in RealInterface.UNARY_ARITHMETIC
    @test isa(eval(Base, f), Function)
end

for f in RealInterface.BINARY_ARITHMETIC
    @test isa(eval(Base, f), Function)
end

for f in RealInterface.UNARY_MATH
    @test isa(eval(Base, f), Function)
end

for f in RealInterface.BINARY_MATH
    @test isa(eval(Base, f), Function)
end

for f in RealInterface.MISC_FUNCTIONS
    @test isa(eval(Base, f), Function)
end

for f in RealInterface.UNARY_PREDICATES
    @test isa(eval(Base, f), Function)
end

for f in RealInterface.BINARY_PREDICATES
    @test isa(eval(Base, f), Function)
end

for f in RealInterface.UNARY_SPECIAL_MATH
    @test isa(eval(SpecialFunctions, f), Function)
end

for f in RealInterface.BINARY_SPECIAL_MATH
    @test isa(eval(SpecialFunctions, f), Function)
end

for f in RealInterface.TERNARY_SPECIAL_MATH
    @test isa(eval(SpecialFunctions, f), Function)
end

for f in RealInterface.UNARY_NAN_MATH
    @test isa(eval(NaNMath, f), Function)
end

for f in RealInterface.BINARY_NAN_MATH
    @test isa(eval(NaNMath, f), Function)
end