text:
# Inspired from `SetProg/test/Test/square.jl`

using LinearAlgebra, Test

using SwitchOnSafety
using Polyhedra
const Sets = SetProg.Sets
using MultivariatePolynomials
using SemialgebraicSets

function ci_square_test(
    factory, variable::SetProg.AbstractVariable, set_test, interval_test; kws...)
    □ = polyhedron(HalfSpace([1, 0], 1.0) ∩ HalfSpace([-1, 0], 1) ∩ HalfSpace([0, 1], 1) ∩ HalfSpace([0, -1], 1))
    interval = polyhedron(HalfSpace([1], 1.0) ∩ HalfSpace([-1], 1))
    Δt = 0.5
    control_system = ConstrainedLinearControlDiscreteSystem(reshape([1.0], 1, 1), reshape([Δt], 1, 1), interval, interval)
    @testset "$(typeof(system))" for system in [
        ConstrainedLinearAlgebraicDiscreteSystem([1.0 Δt], [1.0 0.0], □),
        ConstrainedLinearControlDiscreteSystem([1.0 Δt; 0.0 0.0], reshape([0.0 1.0], 2, 1), □, FullSpace()),
        algebraiclift(control_system)
    ]
        set = invariant_set(system, factory, variable; verbose=0, kws...)
        set_test(set)
    end
    @testset "$(typeof(control_system))" begin
        set = invariant_set(control_system, factory, variable; verbose=0, kws...)
        interval_test(set)
    end
end

const atol = 1e-4
const rtol = 1e-4

@testset "Ellipsoid homogeneous with $factory" for factory in sdp_factories
    ci_square_test(
        factory, Ellipsoid(symmetric=true, dimension=2),
        ◯ -> begin
            @test ◯ isa Sets.Polar{Float64, Sets.EllipsoidAtOrigin{Float64}}
            @test Sets.polar(◯).Q ≈ Symmetric([1.0 -1/4; -1/4 1.0]) atol=atol rtol=rtol
        end,
        interval -> begin
            @test interval isa Sets.Polar{Float64,SetProg.Sets.EllipsoidAtOrigin{Float64}}
            @test Sets.polar(interval).Q ≈ Symmetric(ones(1, 1)) atol=atol rtol=rtol
        end,
        volume_heuristic = nth_root)
end

@testset "Ellipsoid non-homogeneous with $factory" for factory in sdp_factories
    ci_square_test(
        factory, Ellipsoid(point=SetProg.InteriorPoint([0.0, 0.0])),
        ◯ -> begin
            @test ◯ isa Sets.PerspectiveDual{Float64, Sets.Householder{Float64, Sets.ShiftedEllipsoid{Float64}, Float64}}
            z = Sets.perspective_variable(◯)
            x, y = Sets.space_variables(◯)
            ◯_dual = Sets.perspective_dual(◯)
            @test ◯_dual.p ≈ -z^2 + x^2 - x*y/2 + y^2 atol=atol rtol=rtol
            @test ◯_dual.set.Q ≈ Symmetric([1.0 -1/4; -1/4 1.0]) atol=atol rtol=rtol
            @test ◯_dual.set.b ≈ [0.0, 0.0] atol=atol rtol=rtol
            @test ◯_dual.set.β ≈ -1.0 atol=atol rtol=rtol
            @test Sets._householder(◯_dual.h) ≈ [-1.0 0.0 0.0
                                                  0.0 1.0 0.0
                                                  0.0 0.0 1.0] atol=atol rtol=rtol
        end,
        interval -> begin
            @test interval isa Sets.Translation{Sets.Polar{Float64,Sets.EllipsoidAtOrigin{Float64}},Float64,Vector{Float64}}
            @test Sets.polar(interval.set).Q ≈ Symmetric(ones(1, 1)) atol=atol rtol=rtol
            @test interval.c ≈ [0.0] atol=atol rtol=rtol
        end,
        volume_heuristic = nth_root, λ = 1.0)
end

@testset "Quadratic non-homogeneous with $factory" for factory in sdp_factories
    ci_square_test(
        factory, PolySet(degree=2, convex=true, point=SetProg.InteriorPoint([0.0, 0.0])),
        ◯ -> begin
            @test ◯ isa Sets.PerspectiveDual{Float64, Sets.Householder{Float64, Sets.ConvexPolynomialSet{Float64}, Float64}}
            z = Sets.perspective_variable(◯)
            x, y = Sets.space_variables(◯)
            ◯_dual = Sets.perspective_dual(◯)
            # The coefficient of `x*y` does not influence the volume
            # and with the values of the other parameters, it should
            # simply be in the interval [-2, -0.5].
            α = coefficient(◯_dual.p, x*y)
            @test α ≥ -2 - 2atol - rtol
            @test α ≤ -0.5 + 0.5atol + rtol
            @test ◯_dual.p ≈ -z^2 + x^2 + α*x*y + y^2 atol=atol rtol=rtol
        end,
        interval -> begin
            @test interval isa Sets.PerspectiveDual{Float64, Sets.Householder{Float64, Sets.UnknownSet{Float64}, Float64}}
            z = Sets.perspective_variable(interval)
            x, = Sets.space_variables(interval)
            interval_dual = Sets.perspective_dual(interval)
            @test interval_dual.p ≈ -z^2 + x^2 atol=atol rtol=rtol
        end,
        volume_heuristic = set -> L1_heuristic(set, [1.0, 1.0]), λ = 1.0)
end

@testset "Quartic homogeneous with $factory" for factory in sdp_factories
    ci_square_test(
        factory, PolySet(symmetric=true, degree=4, convex=true),
        ◯ -> begin
            @test ◯ isa Sets.Polar{Float64, Sets.ConvexPolynomialSublevelSetAtOrigin{Float64}}
            @test Sets.polar(◯).degree == 4
            x, y = variables(Sets.polar(◯).p)
            α = coefficient(Sets.polar(◯).p, x^3*y) / 2
            q = x^4 + 2α*x^3*y + 6x^2*y^2 + 2α*x*y^3 + y^4
            @test all(eigvals(Matrix(Sets.polar(◯).p.Q)) .≥ -atol)
            @test polynomial(Sets.polar(◯).p) ≈ q atol=atol rtol=rtol
            convexity_proof = Sets.convexity_proof(◯)
            @test convexity_proof.n == 4

            hess = 6 * [2.0, α, 2.0, α, 2.0,
                        2.0, 2.0, α, α, 2.0]
            Hess = SetProg.SumOfSquares.MultivariateMoments.SymMatrix(hess, 4)
            @test all(eigvals(Matrix(Hess)) .≥ -atol)
            @test convexity_proof.Q ≈ hess atol=atol rtol=rtol
        end,
        interval -> begin
            @test interval isa Sets.Polar{Float64, Sets.ConvexPolynomialSublevelSetAtOrigin{Float64}}
            @test Matrix(Sets.polar(interval).p.Q) ≈ ones(1, 1) atol=atol rtol=rtol
        end,
        volume_heuristic = set -> L1_heuristic(set, [1.0, 1.0]))
end


text:
import Pkg
Pkg.add("HTTP")
Pkg.add("JSON")

# Author https://github.com/andrew000

using HTTP
using JSON

TOKEN = ""  # Enter token from https://t.me/BotFather
URL = "https://api.telegram.org/bot$TOKEN"
OFFSET = 0
LIMIT = 5

function get_updates(offset, limit)
	params = Dict("offset" => offset, "limit" => limit)

    response = HTTP.request(
        "GET",
        "$URL/getUpdates",
        ["Content-Type" => "application/json"],
        JSON.json(params))

    return JSON.parse(String(response.body))["result"]
end

function send_message(update)
	if !("message" in keys(update)) || !("text" in keys(update["message"]))
		return
    end

    chat_id = update["message"]["chat"]["id"]
    echo_text = update["message"]["text"]

	params = Dict("chat_id" => chat_id, "text" => echo_text)

    HTTP.request(
        "POST",
        "$URL/sendMessage",
        ["Content-Type" => "application/json"],
		JSON.json(params))
end

while true
    updates = get_updates(OFFSET, LIMIT)

	if length(updates) != 0
		global OFFSET = updates[length(updates)]["update_id"] + 1

        for update in updates
		    send_message(update)
        end
	end
	Base.sleep(0.5)
end

ext: jl


text:
# ---
# title: 928. Minimize Malware Spread II
# id: problem928
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Depth-first Search, Union Find, Graph
# link: <https://leetcode.com/problems/minimize-malware-spread-ii/description/>
# hidden: true
# ---
#
# (This problem is the same as _Minimize Malware Spread_ , with the differences
# bolded.)
#
# In a network of nodes, each node `i` is directly connected to another node `j`
# if and only if `graph[i][j] = 1`.
#
# Some nodes `initial` are initially infected by malware.  Whenever two nodes
# are directly connected and at least one of those two nodes is infected by
# malware, both nodes will be infected by malware.  This spread of malware will
# continue until no more nodes can be infected in this manner.
#
# Suppose `M(initial)` is the final number of nodes infected with malware in the
# entire network, after the spread of malware stops.
#
# We will remove one node from the initial list, **completely removing it and
# any connections from this node to any other node**.  Return the node that if
# removed, would minimize `M(initial)`.  If multiple nodes could be removed to
# minimize `M(initial)`, return such a node with the smallest index.
#
#
#
# **Example 1:**
#
#
#
#     Input: graph = [[1,1,0],[1,1,0],[0,0,1]], initial = [0,1]
#     Output: 0
#
#
# **Example 2:**
#
#
#
#     Input: graph = [[1,1,0],[1,1,1],[0,1,1]], initial = [0,1]
#     Output: 1
#
#
# **Example 3:**
#
#
#
#     Input: graph = [[1,1,0,0],[1,1,1,0],[0,1,1,1],[0,0,1,1]], initial = [0,1]
#     Output: 1
#
#
#
#
# **Note:**
#
#   1. `1 < graph.length = graph[0].length <= 300`
#   2. `0 <= graph[i][j] == graph[j][i] <= 1`
#   3. `graph[i][i] = 1`
#   4. `1 <= initial.length < graph.length`
#   5. `0 <= initial[i] < graph.length`
#
#
## @lc code=start
using LeetCode

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


text:
"""
    Network

Tools for creating axial (i.e., 1d) thermoacoustic network models.
The formulation is Riemann-based:
p=F exp(+iω*l/c)+G exp(-iω*l/c)
Au=A/ρc[ F exp(+iω*l/c)-G exp(-iω*l/c)]
"""
module Network
using WavesAndEigenvalues.NLEVP
export discretize
include("NLEVP_exports.jl")

"""
    out=duct(l,c,A,ρ=1.4*101325/c^2)

Create local matrix coefficients for duct element.

# Arguments
- `l`: length of the duct
- `c`: speed of sound
- `A` cross-sectional area of the duct
- `ρ` density
"""
function duct(l,c,A,ρ=1.4*101325/c^2)
    M=[-1      -1;        #p_previous-p_current=0
        -A/(ρ*c) A/(ρ*c); #A u_previous-A u_current=0
        0          0;     #p_current-p_next=0
        0          0]     #A u_current-A u_next=0

    M31=[0        0;
        0         0;
        1         0;
        0         0]

    M32= [0        0;
        0         0;
        0         1;
        0         0]
    M41= [0        0;
        0        0;
        0         0;
        1         0]
    M42= [0        0;
        0         0;
        0         0;
        0         -1]
    out=[[M,(),()],
         [M31, (generate_exp_az(1im*l/c),), ((:ω,),),],
         [M32, (generate_exp_az(-1im*l/c),), ((:ω,),),],
         [M41*A/(ρ*c), (generate_exp_az(1im*l/c),), ((:ω,),),],
         [M42*A/(ρ*c), (generate_exp_az(-1im*l/c),), ((:ω,),),],
         ]

    return out
end

"""
    out=terminal(R,c,A,ρ=1.4*101325/c^2;init=true)

Create local matrix coefficients for terminal element.

# Arguments
- `R`: reflection coefficient
- `c`: speed of sound
- `A`: cross-sectional area of the terminal
- `ρ`: density
- `init`: (optional) if true its the left end else it is the right one.
"""
function terminal(R,c,A,ρ=1.4*101325/c^2;init=true)
    #if hasmethod(R,(ComplexF64,Int))
    if typeof(R)<:Number
    else
        println("Error, unknown impedance type $(typeof(Z))")
    end

    if init
         M=[+R      -1;
           +1       +1;
           1*A/(ρ*c)  -1*A/(ρ*c)]

    else
        M=[ -1       -1;
            -1*A/(ρ*c)  +1*A/(ρ*c);
          -1     +R]
    end
    return [[M,(),()],]
end

"""
    out=flame(c1,c2,A,ρ=1.4*101325/c1^2)

Create local matrix coefficients for flame element with n-tau dynamics.

# Arguments
- `c1`: speed of sound on the unburnt side
- `c2`: speed of sound on the burnt side
- `A`: cross-sectional area of the terminal
- `ρ`: density
"""
function flame(c1,c2,A,ρ=1.4*101325/c1^2)
    M= [0        0;
        0        0;
        0         0;
        1         -1]

    out=duct(0,c1,A,ρ)
    push!(out,[M*(c2^2/c1^2-1)*A/(ρ*c1), (pow1,exp_delay), ((:n,),(:ω,:τ),),])
    return out
end

"""
    out=helmholtz(V,l_n,d_n,c,A,ρ=1.4*101325/c1^2)

Create local matrix coefficients for a sidewall Helmholtz damper element.

# Arguments
- `V`: Volume of the damper
- `l_n`: length of the damper neck
- `d_n`: diameter of the damper neck
- `c`: speed of sound
- `A`: cross-sectional area of the terminal
- `ρ`: density

# Notes:

The model is build on  the transfer matrix from [1]

# References

[1] F. P. Mechel, Formulas of Acoustic, Springer, 2004, p. 728
"""
function helmholtz(V,l_n,d_n,c,A,ρ=1.4*101325/c^2)
    #=====

    p_u=p_d
    u_u=1/Z p_d+u_d <==> A u_u=A/Z p_d+A u_d>

    TODO:
    <==> Z*A*u_u=A*p_d+Z*Au_d (roe easier representation)

    where

    Z(ω)=ρ[ω^2/(π*c)*[2-r_n/r_u]+0.425*M*c*/S_n+1im[ω*l/S_n-c^2/(ω*V)]]

    where
    ρ: density
    c: speed of sound
    r_n: neck radius of the Helmholtz damper
    r_u: radius of the main duct
    M: Mach number
    S_n = π*r_n2 : crossectional area of the neck
    V = : volume of the helmholtz damper
    and
    l=l_n+t_w+0.85 r_n*(2-r_n/r_u)

    with tw: wall thickness?

    ## adaption to Riemann invariants

    p_u-A_+*exp(im*ω*l/c)+A_-*exp(-im*ω*l/c)=0
    u_u-A_+/ρc*exp(im*ω*l/c)+A_+/ρc*exp(-im*ω*l/c)=0
    The last equation is to be multiplied by A.

    ## notes
    this element does not resolve what is between the usptream and the downstream
    site. It directly describes the jump!
    =====#

    r_n=d_n/2
    r_u=sqrt(A/pi)
    S_n=pi*r_n^2
    t_w=0
    l=l_n+t_w+0.85*r_n*(2-r_n/r_u)
    M=0
    M21= [0        0;
        -1       -1;
        0         0;
        0         0]
    #M22= [0        0;
    #    0         -1;
    #    0         0;
    #    0         0]

    function impedance(ω::ComplexF64,k::Int)::ComplexF64
        let ρ=ρ, r_n=r_n, r_u=r_u, M=M, c=c, l=l, V=V ,S_n=S_n
            return ρ*(pow2(ω,k)/(π*c)*(2-r_n/r_u)+pow0(ω,k)*0.425*M*c/S_n
            +1.0im*pow1(ω,k)*l/S_n-1.0im*c^2/V*pow(ω,k,-1))
        end
    end

    function admittance(ω::ComplexF64,k::Int)::ComplexF64
        let Z=impedance
            if k==0
                return 1/Z(ω,0)
            elseif k==1
                return -Z(ω,1)/Z(ω,0)^2 #TODO: implement arbitray order chain rule
            else
                return NaN+NaN*im
            end
        end
    end
    function admittance(ω::Symbol)::String
        return "1/Z($ω)"
    end

    out=duct(0,c,A,ρ)
    #push!(out,[M21, (admittance,), ((:ω,),),])
    push!(out,[-M21/ρ, (admittance,), ((:ω,),),])
    return out
end

"""
    sidewallimp(imp,c,A,ρ=1.4*101325/c^2)

Use function predefined  function `imp` to prescripe a frequency-dependent
side wall impedance.
"""
function sidewallimp(imp,c,A,ρ=1.4*101325/c1^2)
    M21= [0        0;
        -1       -1;
        0         0;
        0         0]
    function admittance(ω::ComplexF64,k::Int)::ComplexF64
        let Z=imp
            if k==0
                return 1/Z(ω,0)
            elseif k==1
                return -Z(ω,1)/Z(ω,0)^2 #TODO: implement arbitray order chain rule
            else
                return NaN+NaN*im
            end
        end
    end
    function admittance(ω::Symbol)::String
        return "1/Z($ω)"
    end

    out=duct(0,c,A,ρ)
    #push!(out,[M21*A*(ρ*c), (admittance,), ((:ω,),),])
    push!(out,[M21, (admittance,), ((:ω,),),])
end
"""
    out=lhr(V,l_n,d_n,c,A,ρ=1.4*101325/c^2)

Create linear Helmholtz model according to [1].

#Reference
[1] Nonlinear effects in acoustic metamaterial based on a cylindrical pipe with
    ordered Helmholtz resonators,J. Lan, Y. Li, H. Yu, B. Li, and X. Liu, Phys.
    Rev. Lett. A., 2017, [doi:10.1016/j.physleta.2017.01.036](https://doi.org/10.1016/j.physleta.2017.01.036)
"""
function lhr(V,l_n,d_n,c,A,ρ=1.4*101325/c^2)
    r_n=d_n/2
    S_n=pi*r_n^2
    B0=ρ*c^2
    η=1.5E-5 #dynamic viscosity or air in m^2/s
    R_vis=ρ*l_n/r_n*sqrt(η/2)*S_n# *ω^(1/2) #!!!! wrong in Lan it must be divided not miultiplied by r_n
    R_rad=1/4*ρ*r_n^2/c*S_n#*ω^2
    l=l_n+1.7*r_n #length correction
    Cm=V/(ρ*c^2*S_n^2)
    Mm=ρ*l*S_n
    #δ=Rm/Mm
    ω0=1/(sqrt(Cm*Mm))
    println("M: $Mm, C: $Cm, freq: $ω0")
    C=B0*S_n/(1im*ω0^2*V)/(S_n)  #last division to convert between impedance and flow impedance
    function impedance(ω::ComplexF64,k::Int)::ComplexF64
        let C=C, R_vis=R_vis, R_rad=R_rad, ω0=ω0
            return C*pow1(ω,k)-C*ω0^2*pow(ω,k,-1)-C*1im*R_vis/Mm*pow(ω,k,1/2)-C*1im*R_rad/Mm*pow2(ω,k)
        end
    end

    return sidewallimp(impedance,c,A,ρ)
end

"""
    L=discretize(network)

Discretize network model.

# Arguments
- `network`

# Returns
- `L::LinearOperatorFamily`

# The network is specified of a list of network elements.
elements together with there options. Available elements are.
- `:duct`: duct element with options `(l,c,A,ρ)`
- `:flame`: n-tau-flame element with options `(c1,c2,A,ρ)`
- `:helmholtz`:Helmholtz damper element with options `(V,l_n,d_n,c,A,ρ)`
- `:unode`: sound hard boundary, i.e. a velocity node with options `(c,A,ρ)`
- `:pnode`: sound soft boundary, i.e., a pressure node with options `(c,A,ρ)`
- `:anechoic`: anechoic boundary condition with option `(c,A,ρ)`

For all elements the density `ρ` is an optional parameter which is computed
as `ρ=1.4*101325/c^2` if unspecified (air at atmospheric pressure).

# Example
This is an example for the discretization of a Rijke tube, with a Helmholtz
damper, a velocity node at the inlet and a pressure node at the outlet:
    network=[(:unode,(347.0, 0.01), ),
            (:duct,(0.25, 347.0, 0.01)),
            (:flame,(347,347*2,0.01)),
            (:duct,(0.25, 347.0*2, 0.01)),
            (:helmholtz,(0.02^3,0.01,.005,347*2,0.01)),
            (:duct,(0.25, 347.0*2, 0.01)),
            (:pnode,(347.0*2, 0.01), )
            ]

    L=discretize(network)
    L.params[n]=1
    L.params[τ]=0.001
"""
function discretize(network)
    i,j=1,1
    L=LinearOperatorFamily((:ω,))
    N=length(network)
    A=zeros(2*N,2*N)
    for (idx,(elmt, data)) in enumerate(network)
        if elmt==:duct
            terms=duct(data...)
        elseif elmt==:unode
            R=1
            if idx==1
                init=true
            elseif idx==length(network)
                init=false
            else
                println("Error, terminal element at intermediate position $idx!")
            end
            #terms=velocity_node(data...)
            terms=terminal(R,data...,init=init)
        elseif elmt==:pnode
            R=-1
            if idx==1
                init=true
            elseif idx==length(network)
                init=false
            else
                println("Error, terminal element at intermediate position $idx!")
            end
            terms=terminal(R,data...,init=init)
        elseif elmt==:anechoic
            R=0
            if idx==1
                init=true
            elseif idx==length(network)
                init=false
            else
                println("Error, terminal element at intermediate position $idx!")
            end
            terms=terminal(R,data...,init=init)
        elseif elmt==:pnodeend
            terms=pressure_node_end(data...)
        elseif elmt==:vnodestart
            terms=velocity_node_start(data...)
        elseif elmt==:flame
            terms=flame(data...)
        elseif elmt==:helmholtz
            terms=helmholtz(data...)
        elseif elmt==:helmholtz2
            terms=lhr(data...)
        else
            println("Warning: unknown element $elmt")
            continue
        end
        I,J=size(terms[1][1])

        for (coeff,func,arg) in terms
            M=deepcopy(A)
            M[i:i+I-1,j:j+J-1]=coeff
            push!(L,Term(M,func,arg,"M"))
        end
        i+=I-2
        j+=2
    end
    return L
end
end #module network


text:
@with_kw mutable struct PDPWSolver <: AbstractPOMCPSolver
    rng::AbstractRNG        = Base.GLOBAL_RNG
    max_depth::Int          = 20
    c::Float64              = 1.0
    tree_queries::Int       = 1000
    max_time::Float64       = Inf
    k_action::Float64       = 10.0
    alpha_action::Float64   = 1/2
    k_observation::Float64  = 10.0
    alpha_observation::Float64 = 1/2
    enable_action_pw::Bool      = false
    check_repeat_obs::Bool      = false
    check_repeat_act::Bool      = false
    next_action::Any        = RandomActionGenerator(rng)
    default_action::Any     = ExceptionRethrow()
    estimate_value::Any     = RolloutEstimator(RandomSolver(rng))
end

struct PDPWTree{A,O,S}
    # for each observation-terminated history
    total_n::Vector{Int}
    children::Vector{Vector{Int}}
    o_labels::Vector{O}
    B::Vector{Vector{S}}

    # o_lookup::Dict{Tuple{Int, O}, Int}

    # for each action-terminated history
    n::Vector{Int}
    v::Vector{Float64}
    ha_children::Vector{Vector{Int}}
    a_labels::Vector{A}
end

function PDPWTree(pomdp::POMDP, b::ParticleCollection, enable_apw::Bool, sz::Int=1000)
    A = action_type(pomdp)
    O = obs_type(pomdp)
    S = state_type(pomdp)
    sz = min(100_000, sz)
    if enable_apw
        return PDPWTree{A,O,S}(sizehint!(Int[0], sz),
                               sizehint!(Vector{Int}[Int[]], sz),
                               sizehint!(Vector{O}(1), sz),
                               sizehint!([particles(b)], sz),

                               sizehint!(Int[], sz),
                               sizehint!(Float64[], sz),
                               sizehint!(Vector{Int}[], sz),
                               sizehint!(A[], sz)
                             )
    else
        acts = collect(iterator(actions(pomdp)))
        return PDPWTree{A,O,S}(sizehint!(Int[0], sz),
                              sizehint!(Vector{Int}[collect(1:length(acts))], sz),
                              sizehint!(Vector{O}(1), sz),
                              sizehint!([particles(b)], sz),

                              sizehint!(zeros(Int, length(acts)), sz),
                              sizehint!(zeros(Float64, length(acts)), sz),
                              sizehint!(fill(Int[], length(acts)), sz),
                              sizehint!(acts, sz)
                             )
    end
end

function insert_obs_node!(t::PDPWTree, pomdp::POMDP, ha::Int, o, s, enable_apw::Bool)
    push!(t.total_n, 0)
    if enable_apw
        push!(t.children, Int[])
    else
        push!(t.children, sizehint!(Int[], n_actions(pomdp)))
    end
    push!(t.B, [s])
    push!(t.o_labels, o)
    hao = length(t.total_n)
    # t.o_lookup[(ha, o)] = hao
    if !enable_apw
        for a in iterator(actions(pomdp))
            n = insert_action_node!(t, hao, a)
            push!(t.children[hao], n)
        end
    end
    return hao
end

function insert_action_node!(t::PDPWTree, h::Int, a)
    push!(t.n, 0)
    push!(t.v, 0.0)
    push!(t.ha_children, Int[])
    push!(t.a_labels, a)
    return length(t.n)
end

struct PDPWObsNode{A,O} <: BasicPOMCP.BeliefNode
    tree::PDPWTree{A,O}
    node::Int
end

MCTS.isroot(n::PDPWObsNode) = n.node == 1
MCTS.n_children(n::PDPWObsNode) = length(n.tree.children[n.node])

mutable struct PDPWPlanner{P, SE, NA, RNG} <: Policy
    solver::PDPWSolver
    problem::P
    solved_estimator::SE
    next_action::NA
    rng::RNG
    _best_node_mem::Vector{Int}
    _tree::Nullable
end

function PDPWPlanner(solver::PDPWSolver, pomdp::POMDP)
    se = BasicPOMCP.convert_estimator(solver.estimate_value, solver, pomdp)
    @assert solver.check_repeat_obs == false
    @assert solver.check_repeat_act == false
    return PDPWPlanner(solver, pomdp, se, solver.next_action, solver.rng, Int[], Nullable())
end

solve(solver::PDPWSolver, pomdp::POMDP) = PDPWPlanner(solver, pomdp)

Base.srand(p::PDPWPlanner, seed) = srand(p.rng, seed)

function action(p::PDPWPlanner, b)
    local a::action_type(p.problem)
    try
        tree = PDPWTree(p.problem, b, p.solver.enable_action_pw, p.solver.tree_queries)
        a = search(p, b, tree)
        p._tree = Nullable(tree)
    catch ex
        # Note: this might not be type stable, but it shouldn't matter too much here
        a = convert(action_type(p.problem), BasicPOMCP.default_action(p.solver.default_action, p.problem, b, ex))
    end
    return a
end

function search(p::PDPWPlanner, b, t::PDPWTree)
    all_terminal = true
    start_us = CPUtime_us()
    for i in 1:p.solver.tree_queries
        if CPUtime_us() - start_us >= 1e6*p.solver.max_time
            break
        end
        s = rand(p.rng, b)
        if !POMDPs.isterminal(p.problem, s)
            simulate(p, s, PDPWObsNode(t, 1), p.solver.max_depth)
            all_terminal = false
        end
    end

    if all_terminal
        throw(AllSamplesTerminal(b))
    end

    h = 1
    best_node = first(t.children[h])
    best_v = t.v[best_node]
    @assert !isnan(best_v)
    for node in t.children[h][2:end]
        if t.v[node] >= best_v
            best_v = t.v[node]
            best_node = node
        end
    end

    return t.a_labels[best_node]
end

function simulate(p::PDPWPlanner, s, hnode::PDPWObsNode, steps::Int)
    if steps == 0 || isterminal(p.problem, s)
        return 0.0
    end

    t = hnode.tree
    h = hnode.node

    ka = p.solver.k_action
    aa = p.solver.alpha_action
    if p.solver.enable_action_pw && length(t.children[h]) <= ka*t.total_n[h]^aa
        a = next_action(p.next_action, p.problem, ParticleCollection(t.B[h]), hnode)
        n = insert_action_node!(t, h, a)
        push!(t.children[h], n)
    end

    ltn = log(t.total_n[h])
    best_nodes = empty!(p._best_node_mem)
    best_criterion_val = -Inf
    for node in t.children[h]
        n = t.n[node]
        if n == 0 && ltn <= 0.0
            criterion_value = t.v[node]
        elseif n == 0 && t.v[node] == -Inf
            criterion_value = Inf
        else
            criterion_value = t.v[node] + p.solver.c*sqrt(ltn/n)
        end
        if criterion_value > best_criterion_val
            best_criterion_val = criterion_value
            empty!(best_nodes)
            push!(best_nodes, node)
        elseif criterion_value == best_criterion_val
            push!(best_nodes, node)
        end
    end
    ha = rand(p.rng, best_nodes)
    a = t.a_labels[ha]

    if length(t.ha_children[ha]) <= p.solver.k_observation*t.n[ha]^p.solver.alpha_observation
        sp, o, r = generate_sor(p.problem, s, a, p.rng)

        hao = insert_obs_node!(t, p.problem, ha, o, sp, p.solver.enable_action_pw)
        v = BasicPOMCP.estimate_value(p.solved_estimator,
                                      p.problem,
                                      sp,
                                      PDPWObsNode(t, hao),
                                      steps-1)
        R = r + discount(p.problem)*v
    else
        hao = rand(p.rng, t.ha_children[ha])
        sp = rand(p.rng, t.B[hao])
        o = o_labels[hao]
        r = reward(s, a, sp)
        R = r + discount(p.problem)*simulate(p, sp, PDPWObsNode(t, hao), steps-1)
    end

    t.total_n[h] += 1
    t.n[ha] += 1
    t.v[ha] += (R-t.v[ha])/t.n[ha]

    return R
end


text:
function p52(data::Dict)

  # Setup basic dimensions of arrays
  # Parse & check FEdict data

  if :struc_el in keys(data)
    struc_el = data[:struc_el]
  else
    println("No fin_el type specified.")
    return
  end

  ndim = struc_el.ndim
  nst = struc_el.nst

  # Handle :r direction (implicit axisymmetric)
  if struc_el.direction == :r
    nre = struc_el.nxe::Int
    nze = struc_el.nye::Int
  end

  # Add radial stress
  if struc_el.axisymmetric
    nst = 4
  end

  fin_el = struc_el.fin_el
  @assert typeof(fin_el) <: FiniteElement

  if typeof(fin_el) == Line
    (nels, nn) = mesh_size(fin_el, struc_el.nxe)
  elseif typeof(fin_el) == Triangle || typeof(fin_el) == Quadrilateral
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye)
  elseif typeof(fin_el) == Hexahedron
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye, struc_el.nze)
  else
    println("$(typeof(fin_el)) is not a known finite element.")
    return
  end

  nodof = fin_el.nodof           # Degrees of freedom per node
  ndof = fin_el.nod * nodof      # Degrees of freedom per fin_el

  # Update penalty if specified in FEdict

  penalty = 1e20
  if :penalty in keys(data)
    penalty = data[:penalty]
  end

  # Allocate all arrays

  # Start with arrays to be initialized from FEdict

  if :properties in keys(data)
    prop = zeros(size(data[:properties], 1), size(data[:properties], 2))
    for i in 1:size(data[:properties], 1)
      prop[i, :] = data[:properties][i, :]
    end
  else
    println("No :properties key found in FEdict.")
    exit(1)
  end

  nf = ones(Int, nodof, nn)
  if :support in keys(data)
    for i in 1:size(data[:support], 1)
      nf[:, data[:support][i][1]] = data[:support][i][2]
    end
  else
    println("No :supports key found in FEdict.")
    exit(1)
  end

  x_coords = zeros(nn)
  if :x_coords in keys(data)
    x_coords = data[:x_coords]
  end

  y_coords = zeros(nn)
  if :y_coords in keys(data)
    y_coords = data[:y_coords]
  end

  z_coords = zeros(nn)
  if :z_coords in keys(data)
    z_coords = data[:z_coords]
  end

  r_coords = zeros(nn)
  if :r_coords in keys(data)
    r_coords = data[:r_coords]
  end

  etype = ones(Int, nels)
  if :etype in keys(data)
    etype = data[:etype]
  end

  @assert :lth in keys(data)
  @assert :iflag in keys(data)
  @assert :chi in keys(data)
  lth = data[:lth]::Int
  iflag = data[:iflag]::Int
  chi = data[:chi]*pi/180.0
  ca = cos(chi)
  sa = sin(chi)
  radius = 0.0

  # All other arrays

  points = zeros(struc_el.nip, ndim)
  g = zeros(Int, ndof)
  g_coord = zeros(ndim,nn)
  fun = zeros(fin_el.nod)
  coord = zeros(fin_el.nod, ndim)
  gamma = zeros(nels)
  jac = zeros(ndim, ndim)
  g_num = zeros(Int, fin_el.nod, nels)
  der = zeros(ndim, fin_el.nod)
  deriv = zeros(ndim, fin_el.nod)
  bee = zeros(nst,ndof)
  km = zeros(ndof, ndof)
  mm = zeros(ndof, ndof)
  gm = zeros(ndof, ndof)
  kg = zeros(ndof, ndof)
  eld = zeros(ndof)
  weights = zeros(struc_el.nip)
  g_g = zeros(Int, ndof, nels)
  num = zeros(Int, fin_el.nod)
  actions = zeros(ndof, nels)
  displacements = zeros(size(nf, 1), ndim)
  gc = ones(ndim)
  dee = zeros(nst,nst)
  sigma = zeros(nst)
  axial = zeros(nels)

  formnf!(nodof, nn, nf)
  neq = maximum(nf)
  println("There are $(neq) equations.\n")

  # Find global array sizes

  for iel in 1:nels
    geom_rect!(fin_el, iel, r_coords, z_coords, coord, num, struc_el.direction)
    num_to_g!(num, nf, g)
    g_num[:, iel] = num
    g_coord[:, num] = coord'
    g_g[:, iel] = g
  end

  sample!(fin_el, points, weights)

  gsm = spzeros(neq, neq)
  for iel in 1:nels
    deemat!(dee, prop[etype[iel], 1], prop[etype[iel], 2])
    num = g_num[:, iel]
    coord = g_coord[:, num]'              # Transpose
    g = g_g[:, iel]
    km = zeros(ndof, ndof)
    for i in 1:struc_el.nip
      shape_fun!(fun, points, i)
      shape_der!(der, points, i)
      jac = der*coord
      detm = det(jac)
      jac = inv(jac)
      deriv = jac*der
      (radius, bee) = bmat_nonaxi!(bee, radius, coord, deriv, fun, iflag, lth)
      km += (((bee')*dee)*bee)*detm*weights[i]*radius
    end
    fsparm!(gsm, g, km)
  end

  loads = OffsetArray(zeros(neq + 1), 0:neq)
  if :loaded_nodes in keys(data)
    for i in 1:size(data[:loaded_nodes], 1)
      loads[nf[:, data[:loaded_nodes][i][1]]] = data[:loaded_nodes][i][2]
    end
  end

  fixed_freedoms = 0
  if :fixed_freedoms in keys(data)
    fixed_freedoms = size(data[:fixed_freedoms], 1)
  end
  no = zeros(Int, fixed_freedoms)
  node = zeros(Int, fixed_freedoms)
  sense = zeros(Int, fixed_freedoms)
  value = zeros(Float64, fixed_freedoms)
  if :fixed_freedoms in keys(data) && fixed_freedoms > 0
    for i in 1:fixed_freedoms
      no[i] = nf[data[:fixed_freedoms][i][2], data[:fixed_freedoms][i][1]]
      value[i] = data[:fixed_freedoms][i][3]
      gsm[no[i], no[i]] += penalty
      loads[no[i]] = gsm[no[i], no[i]] * value[i]
    end
  end

  cfgsm = cholesky(gsm)
  loads[1:neq] = cfgsm \ loads[1:neq]

  displacements = zeros(size(nf))
  for i in 1:size(displacements, 1)
    for j in 1:size(displacements, 2)
      if nf[i, j] > 0
        displacements[i,j] = loads[nf[i, j]]
      end
    end
  end
  displacements = displacements'

  dis_df = DataFrame(
    r_disp = displacements[:, 1],
    z_disp = displacements[:, 2],
    t_disp = displacements[:, 3]
  )

  loads[0] = 0.0

  struc_el.nip = 1
  points = zeros(struc_el.nip, ndim)
  weights = zeros(struc_el.nip)

  sample!(fin_el, points, weights)

  gc1 = Vector{Float64}()
  gc2 = Vector{Float64}()
  s1 = Vector{Float64}()
  s2 = Vector{Float64}()
  s3 = Vector{Float64}()
  t1 = Vector{Float64}()
  t2 = Vector{Float64}()
  t3 = Vector{Float64}()

  for iel in 1:nels
    deemat!(dee, prop[etype[iel], 1], prop[etype[iel], 2])
    num = g_num[:, iel]
    coord = g_coord[:, num]'
    g = g_g[:, iel]
    eld = loads[g]
    for i in 1:struc_el.nip
      shape_fun!(fun, points, i)
      shape_der!(der, points, i)
      gc = fun'*coord
      jac = inv(der*coord)
      deriv = jac*der
      (radius, bee) = bmat_nonaxi!(bee,radius,coord,deriv,fun,iflag,lth)
      bee[1:4,:]=bee[1:4,:]*ca
      bee[5:6,:]=bee[5:6,:]*sa
      sigma = dee*(bee*eld)
      gc1 = append!(gc1, gc[1])
      gc2 = append!(gc2, gc[2])
      s1 = append!(s1, sigma[1])
      s2 = append!(s2, sigma[2])
      s3 = append!(s3, sigma[3])
      t1 = append!(t1, sigma[4])
      t2 = append!(t2, sigma[5])
      t3 = append!(t3, sigma[6])
    end
  end

  sigma_df = DataFrame(
    r_coord = gc1,
    z_coord = gc2,
    sig_r = s1,
    sig_z = s2,
    sig_t = s3,
    tau_rz = t1,
    tau_zt = t2,
    tau_tr = t3
  )

  fem = PtFEM.jFEM(struc_el, fin_el, ndim, nels, nst, ndof, nn, nodof,
    neq, penalty, etype, g, g_g, g_num, nf, no,
    node, num, sense, actions, bee, coord, gamma, dee,
    der, deriv, displacements, eld, fun, gc, g_coord, jac,
    km, mm, kg, cfgsm, loads, points, prop, sigma, value,
    weights, x_coords, y_coords, z_coords, axial)

    (fem, dis_df, sigma_df)
end

max_stars_repo_path: src/5 Elastic Solids/p52.jl

max_issues_repo_path: src/5 Elastic Solids/p52.jl

max_forks_repo_path: src/5 Elastic Solids/p52.jl


text:
mutable struct Rule
    left::Any
    right::Any
    expr::Expr # original expression
    mode::Symbol # can be :symbolic or :dynamic
end

import Base.==
==(a::Rule, b::Rule) = (a.expr == b.expr) && (a.mode == b.mode)

# operator symbols for simple term rewriting
const rewrite_syms = [:(=>)]
# operator symbols for regular pattern matching rules, "dynamic rules"
# that eval the right side at reduction time.
# might be used to implement big step semantics
const dynamic_syms = [:(|>)]

# symbols for bidirectional equality
const equational_syms = [:(==)]

# symbols for anti-rules
const inequality_syms = [:(!=), :(≠)]

function getmode(r::Expr)
    op = gethead(r)
    if op ∈ dynamic_syms # right hand execution, dynamic rules in egg
        :dynamic
    elseif op ∈ rewrite_syms # right side is quoted, symbolic replacement
        :symbolic
    elseif op ∈ equational_syms
        :equational
    elseif op ∈ inequality_syms
        :unequal
    else
        error(`rule "$r" is not in valid form.\n`)
    end
end

"""
Construct a `Rule` from a quoted expression.
You can also use the [`@rule`] macro to
create a `Rule`.
## Symbolic Rules

Rules defined as `left_hand => right_hand` are
called `symbolic` rules. Application of a `symbolic` Rule
is a replacement of the `left_hand` pattern with
the `right_hand` substitution, with the correct instantiation
of pattern variables. Function call symbols are not treated as pattern
variables, all other identifiers are treated as pattern variables.
Literals such as `5, :e, "hello"` are not treated as pattern
variables.

## Dynamic Rules

Rules defined as `left_hand |> right_hand` are
called `dynamic` rules. Dynamic rules behave like anonymous functions.
Instead of a symbolic substitution, the right hand of
a dynamic `|>` rule is evaluated during rewriting:
matched values are bound to pattern variables as in a
regular function call. This allows for dynamic computation
of

## Type Assertions

Type assertions are supported in the left hand of rules
to match and access literal values both when using classic
rewriting and EGraph based rewriting.
To use a type assertion pattern, add `::T` after
a pattern variable in the `left_hand` of a rule.

---

## Examples

Symbolic rule
```julia
Rule(:(a * b => b * a))
```

Equational rule
```julia
Rule(:(a * b == b * a))
```

Dynamic rule
```julia
Rule(:(a::Number * b::Number |> a*b))
```
"""
function Rule(e::Expr; mod::Module=@__MODULE__)
    e = rmlines(copy(e))
    mode = gethead(e)
    l, r = e.args[isexpr(e, :call) ? (2:3) : (1:2)]

    mode = getmode(e)

    l = interpolate_dollar(l, mod)
    l = df_walk(x -> eval_types_in_assertions(x, mod), l; skip_call=true)
    e.args[isexpr(e, :call) ? 2 : 1] = l
    return Rule(l, r, e, mode)
end

# Global Right Hand Side function cache for dynamic rules.
# Now we're talking.
# TODO use a LRUCache
const RHS_FUNCTION_CACHE = Dict{Tuple{Rule, Module}, Tuple{Vector{Symbol}, Function}}()
const RHS_FUNCTION_CACHE_LOCK = ReentrantLock()

function getrhsfun(r::Rule, m::Module)
    lock(RHS_FUNCTION_CACHE_LOCK) do
        p = (r,m)
        if !haskey(RHS_FUNCTION_CACHE, p)
            z = genrhsfun(r.left, r.right, m)
            RHS_FUNCTION_CACHE[p] = z
        end
        return RHS_FUNCTION_CACHE[p]
    end
end

macro rule(e)
    Rule(e; mod=__module__)
end

# string representation of the rule
function Base.show(io::IO, x::Rule)
    println(io, "Rule(:(", x.expr, "))")
end

"""
Generates a tuple containing the list of formal parameters (`Symbol`s)
and the [`RuntimeGeneratedFunction`](@ref) corresponding to the right hand
side of a `:dynamic` [`Rule`](@ref).
"""
function genrhsfun(left, right, mod::Module)
    # remove type assertions in left hand
    lhs = remove_assertions(left)
    # collect variable symbols in left hand
    lhs_vars = Set{Symbol}()
    df_walk( x -> (if x isa Symbol; push!(lhs_vars, x); end; x), left; skip_call=true )
    params = Expr(:tuple, :_lhs_expr, :_egraph, lhs_vars...)

    ex = :($params -> $right)
    (collect(lhs_vars), closure_generator(mod, ex))
end

# TODO is there anything better than eval to use here?
"""
When creating a theory, type assertions in the left hand contain symbols.
We want to replace the type symbols with the real type values, to fully support
the subtyping mechanism during pattern matching.
"""
function eval_types_in_assertions(x, mod::Module)
    if isexpr(x, :(::))
        !(x.args[1] isa Symbol) && error("Type assertion is not on metavariable")
        x.args[2] isa Type && (return x)
        Expr(:(::), x.args[1], getfield(mod, x.args[2]))
    else x
    end
end

iscond(e) = (isexpr(e, :call) && e.args[1] == :(==)) #||
    #(isexpr(e, :call) && e.args[1] == :(≠))

macro when(expr)
    @assert isexpr(expr, :call)
    op = expr.args[1]
    expr = rmlines(expr)
    if op == :⊢     # syntactical consequence. supported only in egraphs
        conditions = []
        env = expr.args[2]
        # extract conditions from env
        # support (a=b;c=d) [a=b c=d] and [a=b, c=d]
        if isexpr(env, :block) || isexpr(env, :vect) || isexpr(env, :hcat)
            for cond ∈ env.args
                if iscond(cond)
                    push!(conditions, cond)
                else
                    error("malformed condition $cond")
                end
            end
        # support a signle a=b or a≠b
        elseif iscond(env)
            push!(conditions, env)
        else
            error("malformed conditions $env")
        end
        # get rule on right
        rule = expr.args[3]
        op = gethead(rule)
        mode = getmode(rule)
        l, r = rule.args[isexpr(rule, :call) ? (2:3) : (1:2)]

        #
        # println(l)
        # println(r)
        # println(op)
        # println(mode)
        # println(conditions)

        if mode != :symbolic && mode != :dynamic
            error("only conditional dynamic or symbolic rules are supported")
        end

        cond = make_egraph_condition(conditions)
        # println(cond)

        if mode == :dynamic
            ret = :(($op)($l, if $cond
                $r
                else _lhs_expr end)) |> esc |> rmlines
            ret
        end

    end
end

"""
Generate a single expression containing an equality/inequality condition on the
current egraph for a conditional rule and the `@when` macro
"""
# TODO nun va bene, invece di find devi fare tutti gli addexpr
# e poi controllare le uguaglianze sulle EClass alla fine. usa gensym
function make_egraph_condition(conditions)
    egraphed_conditions = []
    for x ∈ conditions
        op, l, r = x.args[1:3]

        ll = df_walk( x -> (if x isa Symbol; dollar(x); end; x), l; skip_call=true )

        l |> dump
        ll |> dump

        op == :≠ && (op = :(!=))
        push!(egraphed_conditions,
            :(($op)(find(_egraph, $l), find(_egraph, $r))))
    end

    foldr((x,y) -> :($x && $y), egraphed_conditions)
end

export @when
export iscond


text:
using Base.depwarn

Base.@deprecate vecbind vcat
Base.@deprecate rbind vcat
Base.@deprecate cbind hcat
Base.@deprecate read_table readtable
Base.@deprecate print_table printtable
Base.@deprecate write_table writetable
import Base.merge
Base.@deprecate merge(df1::AbstractDataFrame, df2::AbstractDataFrame; on::Any = nothing, kind::Symbol = :inner) Base.join
Base.@deprecate colnames(adf::AbstractDataFrame) Base.names
Base.@deprecate colnames!(adf::AbstractDataFrame, vals) names!
Base.@deprecate coltypes(adf::AbstractDataFrame) eltypes
Base.@deprecate types(adf::AbstractDataFrame) eltypes(adf::AbstractDataFrame)
Base.@deprecate EachRow eachrow
Base.@deprecate EachCol eachcol
Base.@deprecate subset sub
Base.@deprecate drop_duplicates! unique!
Base.@deprecate duplicated nonunique
Base.@deprecate load_df loaddf
Base.@deprecate melt_df meltdf
Base.@deprecate stack_df stackdf
Base.@deprecate pivot_table pivottable

const DEFAULT_COLUMN_ELTYPE = Float64
function DataFrame(nrows::Integer, ncols::Integer)
    depwarn("DataFrame(::Integer, ::Integer) is deprecated", :DataFrame)
    columns = Array(Any, ncols)
    for i in 1:ncols
        columns[i] = DataArray(DEFAULT_COLUMN_ELTYPE, nrows)
    end
    cnames = gennames(ncols)
    return DataFrame(columns, Index(cnames))
end

function DataFrame(df::DataFrame)
    depwarn("DataFrame(::DataFrame) is deprecated, use convert(DataFrame, DataFrame) instead",
            :DataFrame)
    return df
end

function DataFrame(x::Union(Number, String))
    depwarn("DataFrame(::Union(Number, String)) is deprecated, use DataFrame(Vector{Any}) instead",
            :DataFrame)
    cols = {DataArray([x], falses(1))}
    colind = Index(gennames(1))
    return DataFrame(cols, colind)
end

# TODO: Replace this with convert call.
# Convert a standard Matrix to a DataFrame w/ pre-specified names
function DataFrame(x::Matrix, cn::Vector = gennames(size(x, 2)))
    depwarn("DataFrame(::Matrix, ::Vector)) is deprecated, use convert(DataFrame, Matrix) instead",
            :DataFrame)
    n = length(cn)
    cols = Array(Any, n)
    for i in 1:n
        cols[i] = DataArray(x[:, i])
    end
    return DataFrame(cols, Index(cn))
end

function DataFrame{T<:String}(columns::Vector{Any}, cnames::Vector{T})
    depwarn("DataFrame(::Vector{Any}, ::Vector{T<:String}) is deprecated, use DataFrame(::Vector{Any}, ::Vector{Symbol}) instead",
            :DataFrame)
    DataFrame(columns, map(symbol, cnames))
end

function DataFrame{D <: Associative, T <: String}(ds::Vector{D}, ks::Vector{T})
    depwarn("DataFrame(::Vector{D<:Associative}, ::Vector{T<:String}) is deprecated, use DataFrame(::Vector{D<:Associative}, ::Vector{Symbol}) instead",
            :DataFrame)
    DataFrame(ds, map(symbol, ks))
end

function DataFrame{K, V}(d::Associative{K, V})
    depwarn("DataFrame(::Vector{Associative}) with broadcasting is deprecated, please remove any reliance on broadcasting",
            :DataFrame)
    # Find the first position with maximum length in the Dict.
    lengths = map(length, values(d))
    max_length = maximum(lengths)
    maxpos = findfirst(lengths .== max_length)
    keymaxlen = keys(d)[maxpos]
    nrows = max_length
    # Start with a blank DataFrame
    df = DataFrame()
    for (k, v) in d
        if length(v) == nrows
            df[k] = v
        elseif rem(nrows, length(v)) == 0    # nrows is a multiple of length(v)
            df[k] = vcat(fill(v, div(nrows, length(v)))...)
        else
            vec = fill(v[1], nrows)
            j = 1
            for i = 1:nrows
                vec[i] = v[j]
                j += 1
                if j > length(v)
                    j = 1
                end
            end
            df[k] = vec
        end
    end
    df
end

# If we have a tuple, convert each value in the tuple to a
# DataVector and then pass the converted columns in, hoping for the best
function DataFrame(vals::Any...)
    depwarn("DataFrame(::Any...) is deprecated",
            :DataFrame)
    p = length(vals)
    columns = Array(Any, p)
    for j in 1:p
        if isa(vals[j], AbstractDataVector)
            columns[j] = vals[j]
        else
            columns[j] = convert(DataArray, vals[j])
        end
    end
    cnames = gennames(p)
    DataFrame(columns, Index(cnames))
end

# Iteration matches that of Associative types (experimental)
function Base.start(df::AbstractDataFrame)
    depwarn("Default AbstractDataFrame iterator is deprecated, use eachcol(::AbstractDataFrame) instead",
            :AbstractDataFrame)
    1
end
Base.done(df::AbstractDataFrame, i) = i > ncol(df)
Base.next(df::AbstractDataFrame, i) = ((names(df)[i], df[i]), i + 1)

##############################################################################
##
## Dict conversion
##
## Try to insure this invertible.
## Allow option to flatten a single row.
##
##############################################################################

function dict(adf::AbstractDataFrame, flatten::Bool = false)
    depwarn("dict(::AbstractDataFrame, ::Bool) is deprecated", :dict)
    res = Dict{Symbol, Any}()
    if flatten && size(adf, 1) == 1
        for colname in names(adf)
            res[colname] = adf[colname][1]
        end
    else
        for colname in names(adf)
            res[colname] = adf[colname]
        end
    end
    return res
end

function pool!(df::AbstractDataFrame, cname::String)
    depwarn("pool!(::AbstractDataFrame, ::String) is deprecated, use pool!(::AbstractDataFrame, ::Symbol) instead", :pool!)
    pool!(df, symbol(cname))
end

function pool!{T<:String}(df::AbstractDataFrame, cname::Vector{T})
    depwarn("pool!(::AbstractDataFrame, ::Vector{T<:String}) is deprecated, use pool!(::AbstractDataFrame, ::Vector{T<:Symbol}) instead", :pool!)
    pool!(df, map(symbol, cnames))
end

function Base.getindex(df::DataFrame, col_ind::String)
    depwarn("indexing DataFrames with strings is deprecated; use symbols instead", :getindex)
    getindex(df, symbol(col_ind))
end

function Base.getindex{T<:String}(df::DataFrame, col_inds::AbstractVector{T})
    depwarn("indexing DataFrames with strings is deprecated; use symbols instead", :getindex)
    getindex(df, map(symbol, col_inds))
end

function Base.getindex(df::DataFrame, row_ind, col_ind::String)
    depwarn("indexing DataFrames with strings is deprecated; use symbols instead", :getindex)
    getindex(df, row_ind, symbol(col_ind))
end

function Base.getindex{T<:String}(df::DataFrame, row_ind, col_inds::AbstractVector{T})
    depwarn("indexing DataFrames with strings is deprecated; use symbols instead", :getindex)
    getindex(df, row_ind, map(symbol, col_inds))
end

function Base.getindex(x::AbstractIndex, idx::String)
    depwarn("indexing DataFrames with strings is deprecated; use symbols instead", :getindex)
    getindex(x, symbol(idx))
end

function Base.getindex{T<:String}(x::AbstractIndex, idx::AbstractVector{T})
    depwarn("indexing DataFrames with strings is deprecated; use symbols instead", :getindex)
    getindex(x, map(symbol, idx))
end

function Base.setindex!(df::DataFrame, v, col_ind::String)
    depwarn("indexing DataFrames with strings is deprecated; use symbols instead", :setindex!)
    setindex!(df, v, symbol(col_ind))
end

function Base.setindex!{T<:String}(df::DataFrame, v, col_inds::AbstractVector{T})
    depwarn("indexing DataFrames with strings is deprecated; use symbols instead", :setindex!)
    setindex!(df, v, map(symbol, col_inds))
end

function Base.setindex!(df::DataFrame, v, row_ind, col_ind::String)
    depwarn("indexing DataFrames with strings is deprecated; use symbols instead", :setindex!)
    setindex!(df, v, row_ind, symbol(col_ind))
end

function Base.setindex!{T<:String}(df::DataFrame, v, row_ind, col_inds::AbstractVector{T})
    depwarn("indexing DataFrames with strings is deprecated; use symbols instead", :setindex!)
    setindex!(df, v, row_ind, map(symbol, col_inds))
end

# reorder! for factors by specifying a DataFrame
function DataArrays.reorder(fun::Function, x::PooledDataArray, df::AbstractDataFrame)
    depwarn("reordering DataFrames is deprecated", :reorder)
    dfc = copy(df)
    dfc["__key__"] = x
    gd = by(dfc, "__key__", df -> colwise(fun, without(df, "__key__")))
    idx = sortperm(gd[[2:ncol(gd)]])
    return PooledDataArray(x, dropna(gd[idx,1]))
end
function DataArrays.reorder(x::PooledDataArray, df::AbstractDataFrame)
    depwarn("reordering DataFrames is deprecated", :reorder)
    reorder(:mean, x, df)
end

function DataArrays.reorder(fun::Function, x::PooledDataArray, y::AbstractVector...)
    depwarn("reordering DataFrames is deprecated", :reorder)
    reorder(fun, x, DataFrame({y...}))
end

function Base.flipud(df::DataFrame)
    depwarn("flipud(DataFrame) is deprecated", :flipud)
    return df[reverse(1:nrow(df)), :]
end

function flipud!(df::DataFrame)
    depwarn("flipud!(DataFrame) is deprecated", :flipud!)
    df[1:nrow(df), :] = df[reverse(1:nrow(df)), :]
    return
end

function cleannames!(df::DataFrame)
    depwarn("cleannames!(::DataFrame) is deprecated (names are now required to be valid identifiers).", :cleannames!)
    oldnames = map(strip, names(df))
    newnames = map(n -> replace(n, r"\W", "_"), oldnames)
    names!(df, newnames)
    return
end

# rbind, cbind, vecbind

rbind(args...) = vcat(args...)

cbind(args...) = hcat(args...)

vecbind_type{T}(::Vector{T}) = Vector{T}
vecbind_type{T<:AbstractVector}(x::T) = Vector{eltype(x)}
vecbind_type{T<:AbstractDataVector}(x::T) = DataVector{eltype(x)}
vecbind_type{T}(::PooledDataVector{T}) = DataVector{T}
vecbind_type(v::StackedVector) = vecbind_promote_type(map(vecbind_type, v.components)...)
vecbind_type(v::RepeatedVector) = vecbind_type(v.parent)
vecbind_type(v::EachRepeatedVector) = vecbind_type(v.parent)

vecbind_promote_type{T1,T2}(x::Type{Vector{T1}}, y::Type{Vector{T2}}) = Array{promote_type(eltype(x), eltype(y)),1}
vecbind_promote_type{T1,T2}(x::Type{DataVector{T1}}, y::Type{DataVector{T2}}) = DataArray{promote_type(eltype(x), eltype(y)),1}
vecbind_promote_type{T1,T2}(x::Type{Vector{T1}}, y::Type{DataVector{T2}}) = DataArray{promote_type(eltype(x), eltype(y)),1}
vecbind_promote_type{T1,T2}(x::Type{DataVector{T1}}, y::Type{Vector{T2}}) = DataArray{promote_type(eltype(x), eltype(y)),1}
vecbind_promote_type(a, b, c, ds...) = vecbind_promote_type(a, vecbind_promote_type(b, c, ds...))
vecbind_promote_type(a, b, c) = vecbind_promote_type(a, vecbind_promote_type(b, c))

function vecbind_promote_type(a::AbstractVector)
    res = None
    if isdefined(a, 1)
        if length(a) == 1
            return a[1]
        else
            if isdefined(a, 2)
                res = vecbind_promote_type(a[1], a[2])
            else
                res = a[1]
            end
        end
    end
    for i in 3:length(a)
        if isdefined(a, i)
            res = vecbind_promote_type(res, a[i])
        end
    end
    return res
end

#constructor{T}(::Type{Vector{T}}, args...) = Array(T, args...)
#constructor{T}(::Type{DataVector{T}}, args...) = DataArray(T, args...)

function vecbind(xs::AbstractVector...)
    V = vecbind_promote_type(map(vecbind_type, {xs...}))
    len = sum(length, xs)
    res = constructor(V, len)
    k = 1
    for i in 1:length(xs)
        for j in 1:length(xs[i])
            res[k] = xs[i][j]
            k += 1
        end
    end
    res
end

function vecbind(xs::PooledDataVector...)
    vecbind(map(x -> convert(DataArray, x), xs)...)
end


text:
using Random
using Statistics
using Printf

mutable struct acqr
    a::Float64
    c::Float64
    q::Float64
    r::Float64
    m1::Float64
    P1::Float64
end
#
struct acqrEM
    y::Array{Float64,1}
    start::acqr
    fit::acqr
    loglik::Array{Float64,1}
    aic::Float64
    s::Bool
end
#
struct acqrKal
    y::Vector{Float64}
    model::acqr
    prediction
    predictionCov
    filter
    filterCov
    smoother
    smootherCov
    s::Bool
end

include("acqr_simula.jl")
include("acqr_kfilter.jl")
include("acqr_ksmoother.jl")
include("acqr_em_starting.jl")
include("acqr_em.jl")
include("acqr_kfilter_s.jl")
include("acqr_ksmoother_s.jl")
include("acqr_em_s.jl")
include("acqr_kalman.jl")


text:
"""
    Inclusivity(first::Bool, last::Bool) -> Inclusivity

Defines whether an `AbstractInterval` is open, half-open, or closed.
"""
struct Inclusivity
    first::Bool
    last::Bool
end

"""
    Inclusivity(i::Integer) -> Inclusivity

Defines whether an interval is open, half-open, or closed, using an integer code:

### Inclusivity Values

* `0`: Neither endpoint is included (the `AbstractInterval` is open)
* `1`: Only the lesser endpoint is included (the `AbstractInterval` is left-closed)
* `2`: Only the greater endpoint is included (the `AbstractInterval` is right-closed)
* `3`: Both endpoints are included (the `AbstractInterval` is closed)

Note that this constructor does not perform bounds-checking: instead it checks the values
of the two least-significant bits of the integer. This means that `Inclusivity(5)` is
equivalent to `Inclusivity(1)`.
"""
Inclusivity(i::Integer) = Inclusivity(i & 0b01 > 0, i & 0b10 > 0)

Base.copy(x::Inclusivity) = Inclusivity(x.first, x.last)

function Base.convert(::Type{I}, x::Inclusivity) where I <: Integer
    return I(x.last << 1 + x.first)
end

# https://github.com/JuliaLang/julia/pull/26601
if VERSION >= v"0.7.0-DEV.4743"
    Base.broadcastable(i::Inclusivity) = Ref(i)
end

Base.first(x::Inclusivity) = x.first
Base.last(x::Inclusivity) = x.last

isclosed(x::Inclusivity) = first(x) && last(x)
Base.isopen(x::Inclusivity) = !(first(x) || last(x))

Base.isless(a::Inclusivity, b::Inclusivity) = isless(convert(Int, a), convert(Int, b))

function Base.show(io::IO, x::T) where T <: Inclusivity
    if get(io, :compact, false)
        print(io, x)
    else
        print(io, "$T($(x.first), $(x.last))")
    end
end

function Base.print(io::IO, x::Inclusivity)
    open_char = first(x) ? '[' : '('
    close_char = last(x) ? ']' : ')'
    desc = isclosed(x) ? "Closed" : isopen(x) ? "Open" : first(x) ? "Left" : "Right"
    print(io, "Inclusivity ", open_char, desc, close_char)
end


text:
using BinaryBuilder

name = "Zstd"
version = v"1.5.2"

sources = [
    ArchiveSource("https://github.com/facebook/zstd/releases/download/v$version/zstd-$version.tar.gz",
                  "7c42d56fac126929a6a85dbc73ff1db2411d04f104fae9bdea51305663a83fd0"),
]

script = raw"""
cd ${WORKSPACE}/srcdir/zstd-*/
mkdir build-zstd && cd build-zstd
meson --cross-file="${MESON_TARGET_TOOLCHAIN}" ../build/meson/
ninja -j${nproc}
ninja install
"""

platforms = supported_platforms()

products = [
    LibraryProduct("libzstd", :libzstd),
    ExecutableProduct("zstd", :zstd),
    ExecutableProduct("zstdmt", :zstdmt),
]

dependencies = Dependency[]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat = "1.6")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


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

using BAT
using Test

using Distributions

@testset "named_tuple_transforms" begin
end


text:
using MultiScaleArrays
using OrdinaryDiffEq, DiffEqBase, Test, StochasticDiffEq, Statistics

struct Cell{B} <: AbstractMultiScaleArrayLeaf{B}
    values::Vector{B}
end
struct Population{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArray{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end
struct Tissue{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArray{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end
struct Embryo{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArrayHead{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end

cell1 = Cell([1.0; 2.0; 3.0])
cell2 = Cell([4.0; 5])
population = construct(Population, deepcopy([cell1, cell2])) # Make a Population from cells
cell3 = Cell([3.0; 2.0; 5.0])
cell4 = Cell([4.0; 6])
population2 = construct(Population, deepcopy([cell3, cell4]))
tissue1 = construct(Tissue, deepcopy([population, population2])) # Make a Tissue from Populations
tissue2 = construct(Tissue, deepcopy([population2, population]))
_embryo = construct(Embryo, deepcopy([tissue1, tissue2])) # Make an embryo from Tissues
embryo = deepcopy(_embryo)

cell_ode = function (dcell,cell,p,t)
    m = mean(cell)
    for i in eachindex(cell)
        dcell[i] = -m*cell[i]
    end
end

f = function (dembryo,embryo,p,t)
    for (cell, dcell) in LevelIter(3,embryo, dembryo)
        cell_ode(dcell,cell,p,t)
    end
end

tstop = [0.5]

condition = function (u, t, integrator)
    t ∈ tstop
end

affect! = function (integrator)
    add_node!(integrator, integrator.u[1, 1, 1], 1, 1)
end

growing_cb = DiscreteCallback(condition, affect!)

println("Do the ODE Part")

prob = ODEProblem(f, embryo, (0.0, 1.0))
test_embryo = deepcopy(embryo)

sol = solve(prob, Tsit5(), callback=growing_cb, tstops=tstop)
sol = solve(prob, Rosenbrock23(autodiff=false), tstops=tstop)
sol = solve(prob, Rosenbrock23(autodiff=false), callback=growing_cb, tstops=tstop)
sol = solve(prob, Rosenbrock23(), callback=growing_cb,tstops=tstop)

@test length(sol[end]) == 23

affect_del! = function (integrator)
    remove_node!(integrator, 1, 1, 1)
end

shrinking_cb = DiscreteCallback(condition, affect_del!)

sol = solve(prob, Tsit5(), callback=shrinking_cb, tstops=tstop)

sol = solve(prob, Rosenbrock23(autodiff=false),callback=shrinking_cb,tstops=tstop)

sol = solve(prob, Rosenbrock23(),callback=shrinking_cb,tstops=tstop)

@test length(sol[end]) == 17

println("Do the SDE Part")

g = function (du,u,p,t)
    for i in eachindex(u)
        du[i] = 0.1u[i]
    end
end
prob = SDEProblem(f, g, embryo, (0.0, 1.0))

@show SRIW1

sol = solve(prob, SRIW1(), callback=growing_cb, tstops=tstop)

@show SRA1

sol = solve(prob, SRA1(), callback=growing_cb, tstops=tstop)

@show RKMil

sol = solve(prob, RKMil(), callback=growing_cb, dt=1/10, tstops=tstop)

@show EM

sol = solve(prob, EM(), dt=1/20, callback=growing_cb, tstops=tstop)

@test length(sol[end]) == 23

@show SRIW1

sol = solve(prob, SRIW1(), callback=shrinking_cb, tstops=tstop)

@show SRA1

sol = solve(prob, SRA1(), callback=shrinking_cb, tstops=tstop)

@show RKMil

sol = solve(prob, RKMil(), dt=1/10, callback=shrinking_cb, tstops=tstop)

@show EM

sol = solve(prob, EM(), dt=1/10, callback=shrinking_cb, tstops=tstop)

@test length(sol[end]) == 17


text:
using Distributions

NORMAL_CONSTANT = sqrt(2*pi)

struct GaussianEstimator
    weight_sum::Float64
    mean::Float64
    variance_sum::Float64
    min_seen::Float64
    max_seen::Float64
end

function GaussianEstimator()
    return GaussianEstimator(0.0, 0.0, 0.0, typemax(Float64), typemin(Float64))
end

function add_weight(g::GaussianEstimator, value::Float64, weight::Float64)
    if g.weight_sum == 0
        return GaussianEstimator(weight, value, 0.0, value, value)
    else
        new_weight = g.weight_sum + weight
        new_mean = g.mean + weight * (value - g.mean) / new_weight
        new_variance = g.variance_sum + weight * (value - g.mean) * (value - new_mean)
        new_min = min(g.min_seen, value)
        new_max = max(g.max_seen, value)
        return GaussianEstimator(new_weight, new_mean, new_variance, new_min, new_max)
    end
end

function variance(g::GaussianEstimator)
    return g.weight_sum > 1.0 ? g.variance_sum / (g.weight_sum - 1.0) : 0.0
end
function stdev(g::GaussianEstimator)
    return sqrt(variance(g))
end

function pdf(g::GaussianEstimator, value::Float64)
    if g.weight_sum > 0.0
        std_dev = stdev(g)
        if std_dev > 0.0
            diff = value - g.mean
            return ((1.0 / (NORMAL_CONSTANT * std_dev)) * exp(-(diff^2 / (2.0 * std_dev^2))))
        elseif value == g.mean
            return 1.0
        end
    end
    return 0.0
end

"""
Estimates the amount of weight seen for a given class value for attribute
values less than, equal to and greater than the given value.

This is used to estimate the future distribution if a split point is added
at the value.
"""
function estimated_weight_split(g::GaussianEstimator, value::Float64)
    equalto_weight = pdf(g, value) * g.weight_sum
    std_dev = stdev(g)
    if std_dev > 0.0
        lessthan_weight = (cdf(Normal(g.mean, std_dev), value) * g.weight_sum) - equalto_weight
    else
        if value < g.mean
            lessthan_weight = g.weight_sum - equalto_weight
        else
            lessthan_weight = 0.0
        end
    end
    greaterthan_weight = max(g.weight_sum - equalto_weight - lessthan_weight, 0.0)

    return (lessthan_weight, equalto_weight, greaterthan_weight)
end


text:
#executable for launch tests
import Pkg

Pkg.activate(".")
Pkg.test()


text:
module EmergencyBrakingSystem

using Random
using Statistics
using AutomotiveDrivingModels
using AutoViz
using AutomotiveSensors
using AutomotivePOMDPs
using Parameters
using StaticArrays
using Reel

export
    EmergencySystem,
    generate_scenario,
    evaluateScenarioMetric,
    evaluateScenariosMetric

include("rendering.jl")
include("evaluation_scenario.jl")
include("emergency_system.jl")

end


text:
module Degeneracy

# """
# Degeneracy Module:
#     The local polytope has no preference for input, output, and party labels.
#     As a result, vertices and facet inequalities are computed which are degenerate
#     under permutation of inputs, outputs, or parties. This module houses methods
#     for identifying, enumerating, and eliminating such degeneracies.
# """

using ..Behavior

using LinearAlgebra
using QBase: permutation_matrices

# """
# symmetry_groups:
#     A symmetry group is generated by applying each the transforms onto a "canonical"
#     group element. This function constructs each of the unique groups that can
#     be generated from the set of elements.
#
# Inputs:
#     generators: Array of Matrices, generator of the group symmetry
#     elements: Array of row vectors, "canonical" elements on which a generator acts
#     type: optional String, "facet" or "vertex" denotes what the element is.
#
# Output:
#     symmetry_groups: Array, group elements are stored as arrays of row vectors.
# """
function symmetry_groups(generators, elements; type="facet")

    symmetry_groups = []

    num_generators = length(generators)
    num_transforms = 2^(num_generators)

    transforms = []
    for el in elements

        # check if element exists in symmetry group
        grouped_elements = collect(Iterators.flatten(symmetry_groups))

        if in(el,grouped_elements)
            continue
        end

        # form symmetry groups by applying all transforms to each element
        new_group = []
        for transform in generators

            if type == "facet"
                e = el*transform
            elseif type == "vertex"
                e = transform*el
            end

            if !in(e,new_group)
                push!(new_group, e)
            end
        end

        push!(symmetry_groups, new_group)
    end

    symmetry_groups
end

# """
# merge_symmetry_groups:
#     Combines two sets of groups such that groups with overlapping terms are merged.
#     The resulting sets are the group sets for the combined symmetry generators.
#
# Input:
#     group1/2_sets: Array, contains arrays of group elements. In practice, group1
#         and group2 should contain different partitions of the same elements.
#
# Output:
#     merge_group_sets: Array, contains arrays of group elements which are closed
#         under the generators of group1 and group2.
# """
function merge_symmetry_groups(group1_sets, group2_sets)

    merge_group_sets = group2_sets

    for group1_set in group1_sets

        # get the indices of the group2 sets which overlap with the group1 set
        overlap_ids = []

        merge_group_id = 1
        for merge_group_set in merge_group_sets

            group_overlap_exists = !all(
                (x) -> x === nothing,
                indexin(group1_set, merge_group_set)
            )
            if group_overlap_exists
                push!(overlap_ids, merge_group_id)
            end

            merge_group_id += 1
        end
        # reconstruct group2 s.t. group1[g1_id] is merged with group2[]
        merge_group = unique(cat(group1_set,  collect(Iterators.flatten(merge_group_sets[overlap_ids])), dims=1))
        merge_group_sets = merge_group_sets[ filter( (x) -> !in(x,overlap_ids), 1:end) ]
        push!(merge_group_sets, merge_group)
    end

    merge_group_sets
end

# """
# canonical_permutation_groups:
#     Finds the fundamental groups closed under input/output permutations for the
#     specified bell scenario.
#
# Input:
#     α/β_expt: Tuple, (# inputs, # outputs)
#     facets: Array, the complete set of facets (row vectors) for the bell scenario
#     rep: String, the representation used to write the facet vectors
#
# Output:
#     canonical_groups: Array, contains arrays of facets which represent a group
#         closed under bell permutation transformations.
# """
function canonical_facets(α_expt, β_expt, facets, rep="generalized")

    input_relabels = bipartite_input_relabels(α_expt, β_expt, rep)
    (α_out_relabels, β_out_relabels) = bipartite_output_relabels(α_expt, β_expt, rep)

    α_out_sym_groups = symmetry_groups(α_out_relabels, facets)
    β_out_sym_groups = symmetry_groups(β_out_relabels, facets)

    out_sym_groups = merge_symmetry_groups(α_out_sym_groups, β_out_sym_groups)
    in_sym_groups = symmetry_groups(input_relabels, facets)

    canonical_groups = merge_symmetry_groups(out_sym_groups, in_sym_groups)

    canonical_groups
end

# """
# canonical_vertices(α, β, vertices, rep):
#
#     Uses the permutation symmetries to find the vertex-transitive groups.
#
# Input:
#     α/β: Tuple (# inputs, # outputs)
#     vertices: Array, contains a complete set of polytope vertices
#     rep: String, specifies the representation of the vertices
#
# Output:
#     canonical_groups: Array, contains vertices grouped by vertex transitivity
#                       with respect to the provided permutation symmetry groups.
# """
function canonical_vertices(α, β, vertices, rep="generalized")

    input_relabels = bipartite_input_relabels(α, β, rep)
    (α_out_relabels, β_out_relabels) = bipartite_output_relabels(α, β, rep)

    α_out_sym_groups = symmetry_groups(α_out_relabels, vertices, type="vertex")
    β_out_sym_groups = symmetry_groups(β_out_relabels, vertices, type="vertex")

    out_sym_groups = merge_symmetry_groups(α_out_sym_groups, β_out_sym_groups)
    in_sym_groups = symmetry_groups(input_relabels, vertices, type="vertex")

    canonical_groups = merge_symmetry_groups(out_sym_groups, in_sym_groups)

    canonical_groups
end

# """
# bipartite_input_relabels:
#   Returns the set of input relabeling matrices for the specified bipartite
#   bell experiment. A relabeling matrix switches the index of an input such that
#   (x <-> x') and/or (y <-> y').
#
# Inputs:
#   (α/β)_expt: Tuple, experiment parameters (Int::num_inputs, Int::num_outputs)
#   rep: String, specify the matrix representation used
#
# Output:
#   maps: Array, contains the complete set of matrices which transform a behavior
#         to its relabeled counterpart (x <-> x') and/or (y <-> y')
# """
function bipartite_input_relabels(α_expt, β_expt, rep="generalized")
    (α_num_inputs, α_num_outputs) = α_expt
    (β_num_inputs, β_num_outputs) = β_expt

    gen_maps = generalized_input_relabels(α_expt, β_expt, true)

    maps = []

    if rep == "generalized"
        maps = gen_maps
    else
        proj = []
        proj_inv = []

        if rep == "non-signaling"
            proj = Behavior.gen_to_ns_proj(α_expt,β_expt)
            proj_inv = Behavior.ns_to_gen_proj(α_expt,β_expt)
        elseif rep == "fixed-direction"
            proj = Behavior.gen_to_fd_proj(α_expt,β_expt)
            proj_inv = Behavior.fd_to_gen_proj(α_expt,β_expt)
        elseif rep == "normalized"
            proj = Behavior.gen_to_norm_proj(α_expt,β_expt)
            proj_inv = Behavior.norm_to_gen_proj(α_expt,β_expt)
        end

        for map in gen_maps
            push!(maps, proj*map*proj_inv)
        end
    end

    maps
end

# """
# generalized_input_relabels:
#     Constructs relabeling maps for the generalized subspace. Maps are constructed
#     by enumerating all permutations of inputs.
#
# Inputs:
#   expt: Tuple, experiment paramterization (num_inputs, num_outputs)
#   include_constant: Bool, regularization constant is added to relabeling maps
#
# Output:
#   maps: Array, contains relabeling matrices for the generalized subspace
# """
function generalized_input_relabels(α_expt, β_expt, include_constant=false)
    (α_num_inputs, α_num_outputs) = α_expt
    (β_num_inputs, β_num_outputs) = β_expt

    α_maps = permutation_matrices(α_num_inputs)
    β_maps = permutation_matrices(β_num_inputs)

    α_id = diagm(0 => fill(1,α_num_outputs))
    β_id = diagm(0 => fill(1,β_num_outputs))

    maps = []
    for α_map in α_maps
        for β_map in β_maps

            map = kron(α_id, kron(β_id, kron(α_map,β_map)))
            if include_constant
                map = cat([[1][:,:], map]..., dims=(1,2))
            end

            push!(maps, map)
        end
    end

    maps
end

# """
# bipartite_output_relabels(α_expt, β_expt, rep)
#     Enumerates the generalized output relabeling maps and projects the set onto
#     the specified subspace representation.
#
# Inputs:
#     (α/β)_expt: Tuple, the relabeling map
#     rep: String, subspace representation
#
# Outputs:
#     relabels: Array, set of unique output relabeling maps
# """
function bipartite_output_relabels(α_expt, β_expt, rep="normalized")

    (α_num_inputs, α_num_outputs) = α_expt
    (β_num_inputs, β_num_outputs) = β_expt

    α_relabels = []
    β_relabels = []

    gen_relabels = generalized_output_relabels(α_expt,β_expt,true)

    if rep == "generalized"
        (α_relabels, β_relabels) = gen_relabels
    else
        proj = []
        proj_inv = []

        if rep == "non-signaling"
            proj = Behavior.gen_to_ns_proj(α_expt,β_expt)
            proj_inv = Behavior.ns_to_gen_proj(α_expt,β_expt)
        elseif rep == "fixed-direction"
            proj = Behavior.gen_to_fd_proj(α_expt,β_expt)
            proj_inv = Behavior.fd_to_gen_proj(α_expt,β_expt)
        elseif rep == "normalized"
            proj = Behavior.gen_to_norm_proj(α_expt,β_expt)
            proj_inv = Behavior.norm_to_gen_proj(α_expt,β_expt)
        end

        for α_relabel in gen_relabels[1]
            push!(α_relabels, proj*α_relabel*proj_inv)
        end

        for β_relabel in gen_relabels[2]
            push!(β_relabels, proj*β_relabel*proj_inv)
        end
    end

    relabels = (unique(α_relabels), unique(β_relabels))

    relabels
end

# """
# generalized_output_relabels:
#     enumerates the output relabeling maps in the generalized representation.
#
# Inputs:
#     α/β_expt: Tuple, experiment parameters
#     include_constant: Boolean, include regularization term
#
# Outputs:
#   α_maps: Array of output relabeling matrices for alice's outputs
#   β_maps: Array of output relabeling matrices for bob's outputs
# """
function generalized_output_relabels(α_expt,β_expt, include_constant)
    (α_num_inputs, α_num_outputs) = α_expt
    (β_num_inputs, β_num_outputs) = β_expt

    α_perm_maps = permutation_matrices(α_num_outputs)
    β_perm_maps = permutation_matrices(β_num_outputs)

    α_in_id = diagm( 0 => fill(1,α_num_inputs))
    β_in_id = diagm( 0 => fill(1,β_num_inputs))
    α_out_id = diagm( 0 => fill(1,α_num_outputs))
    β_out_id = diagm( 0 => fill(1,β_num_outputs))

    α_maps = []
    β_maps = []
    num_inputs = α_num_inputs*β_num_inputs

    # conditional output relabelings for alice
    α_map_index = []
    for x in 1:α_num_inputs

        x_m = zeros(Int64,(α_num_inputs,α_num_inputs))
        x_m[x,x] = 1

        x_rest = α_in_id - x_m
        αx_maps = []
        for α_perm in α_perm_maps

            map = kron(α_perm, kron(β_out_id,kron(x_m,β_in_id))) + kron(α_out_id, kron(β_out_id,kron(x_rest,β_in_id)))
            if include_constant
                map = cat([[1][:,:], map]..., dims=(1,2))
            end

            push!(αx_maps, map)
        end

        push!(α_map_index, αx_maps)
    end

    if α_num_outputs == 1
        push!(α_maps, α_map_index[1][1])
    else
        α_num_transforms = length(α_perm_maps)^α_num_inputs
        for i in 0:(α_num_transforms - 1)
            id_array = digits(i, base = length(α_perm_maps), pad = α_num_inputs) .+ 1
            push!(α_maps, prod(map(
                (id) -> begin
                    (α_map_index[id])[id_array[id]]
                end,
                1:α_num_inputs
            )))
        end
    end
    # conditional output relabelings for bob

    β_map_index = []
    for y in 1:β_num_inputs

        y_m = zeros(Int64,(β_num_inputs,β_num_inputs))
        y_m[y,y] = 1

        y_rest = β_in_id - y_m
        βy_maps = []
        for β_perm in β_perm_maps

            map = kron(α_out_id, kron(β_perm,kron(α_in_id,y_m))) + kron(α_out_id, kron(β_out_id,kron(α_in_id,y_rest)))
            if include_constant
                map = cat([[1][:,:], map]..., dims=(1,2))
            end

            push!(βy_maps, map)
        end

        push!(β_map_index, βy_maps)
    end

    β_num_transforms = length(β_perm_maps)^β_num_inputs
    if β_num_outputs == 1
        push!(β_maps, β_map_index[1][1])
    else
        for i in 0:(β_num_transforms - 1)
            id_array = digits(i, base = length(β_perm_maps), pad = β_num_inputs) .+ 1
            push!(β_maps, prod(map(
                (id) -> begin
                    (β_map_index[id])[id_array[id]]
                end,
                1:β_num_inputs
            )))
        end
    end

    (α_maps, β_maps)
end

end


text:
function get_from(d::Dict{DataType,SpeciesDynamics}, func)
    muddled = d |> values .|> func |> x -> filter(y -> !isnothing(y), x)
    isempty(muddled) && return Vector{Num}[]
    return reduce(vcat, muddled)::Vector{Num}
end

function vardivide(v::Num...)
    states = [filter(!ModelingToolkit.isparameter, v)...]
    params = [filter(ModelingToolkit.isparameter, v)...]
    return states, params
end


text:
#!/usr/bin/env julia --project
using ClimateMachine

using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.Diagnostics
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.Mesh.Filters
using Thermodynamics.TemperatureProfiles
using Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates

using ClimateMachine.BalanceLaws
import ClimateMachine.BalanceLaws:
    vars_state,
    prognostic_vars,
    flux,
    eq_tends,
    integral_load_auxiliary_state!,
    integral_set_auxiliary_state!,
    reverse_integral_load_auxiliary_state!,
    reverse_integral_set_auxiliary_state!

using ArgParse
using UnPack
using Distributions
using Random
using StaticArrays
using Test
using DocStringExtensions
using LinearAlgebra
using Printf

using CLIMAParameters
using CLIMAParameters.Planet: cp_d, MSLP, grav, LH_v0

using CLIMAParameters.Atmos.Microphysics

struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

# ------------------------ Begin Radiation Model ---------------------- #

"""
    DYCOMSRadiationModel <: RadiationModel

## References
 - [Stevens2005](@cite)
"""
struct DYCOMSRadiationModel{FT} <: RadiationModel
    "mass absorption coefficient `[m^2/kg]`"
    κ::FT
    "Troposphere cooling parameter `[m^(-4/3)]`"
    α_z::FT
    "Inversion height `[m]`"
    z_i::FT
    "Density"
    ρ_i::FT
    "Large scale divergence `[s^(-1)]`"
    D_subsidence::FT
    "Radiative flux parameter `[W/m^2]`"
    F_0::FT
    "Radiative flux parameter `[W/m^2]`"
    F_1::FT
    "is AtmosLES moisture model an equilibrium model"
    equilibrium_moisture_model::Bool
end

struct DYCOMSRadiation <: TendencyDef{Flux{FirstOrder}} end

eq_tends(pv::Energy, ::DYCOMSRadiationModel, ::Flux{FirstOrder}) =
    (DYCOMSRadiation(),)

function flux(::Energy, ::DYCOMSRadiation, atmos, args)
    @unpack state, aux = args
    m = radiation_model(atmos)
    FT = eltype(state)
    z = altitude(atmos, aux)
    Δz_i = max(z - m.z_i, -zero(FT))
    # Constants
    upward_flux_from_cloud = m.F_0 * exp(-aux.∫dnz.radiation.attenuation_coeff)
    upward_flux_from_sfc = m.F_1 * exp(-aux.∫dz.radiation.attenuation_coeff)
    param_set = parameter_set(atmos)
    free_troposphere_flux =
        m.ρ_i *
        FT(cp_d(param_set)) *
        m.D_subsidence *
        m.α_z *
        cbrt(Δz_i) *
        (Δz_i / 4 + m.z_i)
    F_rad =
        upward_flux_from_sfc + upward_flux_from_cloud + free_troposphere_flux
    ẑ = vertical_unit_vector(atmos, aux)
    return F_rad * ẑ
end

vars_state(m::DYCOMSRadiationModel, ::Auxiliary, FT) = @vars(Rad_flux::FT)

vars_state(m::DYCOMSRadiationModel, ::UpwardIntegrals, FT) =
    @vars(attenuation_coeff::FT)
function integral_load_auxiliary_state!(
    m::DYCOMSRadiationModel,
    integrand::Vars,
    state::Vars,
    aux::Vars,
)
    FT = eltype(state)

    if m.equilibrium_moisture_model
        integrand.radiation.attenuation_coeff =
            state.ρ * m.κ * aux.moisture.q_liq
    else
        integrand.radiation.attenuation_coeff = m.κ * state.moisture.ρq_liq
    end
end
function integral_set_auxiliary_state!(
    m::DYCOMSRadiationModel,
    aux::Vars,
    integral::Vars,
)
    integral = integral.radiation.attenuation_coeff
    aux.∫dz.radiation.attenuation_coeff = integral
end

vars_state(m::DYCOMSRadiationModel, ::DownwardIntegrals, FT) =
    @vars(attenuation_coeff::FT)
function reverse_integral_load_auxiliary_state!(
    m::DYCOMSRadiationModel,
    integrand::Vars,
    state::Vars,
    aux::Vars,
)
    FT = eltype(state)
    integrand.radiation.attenuation_coeff = aux.∫dz.radiation.attenuation_coeff
end
function reverse_integral_set_auxiliary_state!(
    m::DYCOMSRadiationModel,
    aux::Vars,
    integral::Vars,
)
    aux.∫dnz.radiation.attenuation_coeff = integral.radiation.attenuation_coeff
end

# -------------------------- End Radiation Model ------------------------ #

"""
  Initial Condition for DYCOMS_RF01 LES

## References
 - [Stevens2005](@cite)
"""
function init_dycoms!(problem, bl, state, aux, localgeo, t)
    FT = eltype(state)

    (x, y, z) = localgeo.coord
    param_set = parameter_set(bl)

    z = altitude(bl, aux)

    # These constants are those used by Stevens et al. (2005)
    qref = FT(9.0e-3)
    q_pt_sfc = PhasePartition(qref)
    Rm_sfc = gas_constant_air(param_set, q_pt_sfc)
    T_sfc = FT(290.4)
    _MSLP = FT(MSLP(param_set))
    _grav = FT(grav(param_set))

    # Specify moisture profiles
    q_liq = FT(0)
    q_ice = FT(0)
    zb = FT(600)         # initial cloud bottom
    zi = FT(840)         # initial cloud top

    if z <= zi
        θ_liq = FT(289.0)
        q_tot = qref
    else
        θ_liq = FT(297.0) + (z - zi)^(FT(1 / 3))
        q_tot = FT(1.5e-3)
    end

    ugeo = FT(7)
    vgeo = FT(-5.5)
    u, v, w = ugeo, vgeo, FT(0)

    # Perturb initial state to break symmetry and trigger turbulent convection
    r1 = FT(rand(Uniform(-0.001, 0.001)))
    if z <= 200.0
        θ_liq += r1 * θ_liq
    end

    # Pressure
    H = Rm_sfc * T_sfc / _grav
    p = _MSLP * exp(-z / H)

    # Density, Temperature

    ts = PhaseEquil_pθq(param_set, p, θ_liq, q_tot)
    ρ = air_density(ts)

    e_kin = FT(1 / 2) * FT((u^2 + v^2 + w^2))
    e_pot = gravitational_potential(bl.orientation, aux)
    E = ρ * total_energy(e_kin, e_pot, ts)

    state.ρ = ρ
    state.ρu = SVector(ρ * u, ρ * v, ρ * w)
    state.energy.ρe = E

    state.moisture.ρq_tot = ρ * q_tot

    if moisture_model(bl) isa NonEquilMoist
        q_init = PhasePartition(ts)
        state.moisture.ρq_liq = q_init.liq
        state.moisture.ρq_ice = q_init.ice
    end
    if precipitation_model(bl) isa RainModel
        state.precipitation.ρq_rai = FT(0)
    end

    return nothing
end

function config_dycoms(
    ::Type{FT},
    N,
    resolution,
    xmax,
    ymax,
    zmax,
    moisture_model = "equilibrium",
    precipitation_model = "noprecipitation",
) where {FT}
    # Reference state
    T_profile = DecayingTemperatureProfile{FT}(param_set)
    ref_state = HydrostaticState(T_profile)

    # Radiation model
    κ = FT(85)
    α_z = FT(1)
    z_i = FT(840)
    ρ_i = FT(1.13)

    D_subsidence = FT(3.75e-6)

    F_0 = FT(70)
    F_1 = FT(22)
    if moisture_model == "equilibrium"
        equilibrium_moisture_model = true
    else
        equilibrium_moisture_model = false
    end
    radiation = DYCOMSRadiationModel{FT}(
        κ,
        α_z,
        z_i,
        ρ_i,
        D_subsidence,
        F_0,
        F_1,
        equilibrium_moisture_model,
    )

    # Sources
    f_coriolis = FT(0.762e-4)
    u_geostrophic = FT(7.0)
    v_geostrophic = FT(-5.5)
    w_ref = FT(0)
    u_relaxation = SVector(u_geostrophic, v_geostrophic, w_ref)
    # Sponge
    c_sponge = 1
    # Rayleigh damping
    zsponge = FT(1000.0)
    rayleigh_sponge =
        RayleighSponge{FT}(zmax, zsponge, c_sponge, u_relaxation, 2)
    # Geostrophic forcing
    geostrophic_forcing =
        GeostrophicForcing{FT}(f_coriolis, u_geostrophic, v_geostrophic)

    # Boundary conditions
    # SGS Filter constants
    C_smag = FT(0.21) # 0.21 for stable testing, 0.18 in practice
    C_drag = FT(0.0011)
    LHF = FT(115)
    SHF = FT(15)
    moisture_flux = LHF / FT(LH_v0(param_set))

    source = (
        Gravity(),
        rayleigh_sponge,
        Subsidence{FT}(D_subsidence),
        geostrophic_forcing,
    )

    # moisture model and its sources
    if moisture_model == "equilibrium"
        moisture = EquilMoist(; maxiter = 4, tolerance = FT(1))
    elseif moisture_model == "nonequilibrium"
        source = (source..., CreateClouds())
        moisture = NonEquilMoist()
    else
        @warn @sprintf(
            """
%s: unrecognized moisture_model in source terms, using the defaults""",
            moisture_model,
        )
        moisture = EquilMoist(; maxiter = 4, tolerance = FT(1))
    end

    # precipitation model and its sources
    if precipitation_model == "noprecipitation"
        precipitation = NoPrecipitation()
    elseif precipitation_model == "rain"
        source = (source..., WarmRain_1M())
        precipitation = RainModel()
    else
        @warn @sprintf(
            """
%s: unrecognized precipitation_model in source terms, using the defaults""",
            precipitation_model,
        )
        precipitation = NoPrecipitation()
    end

    physics = AtmosPhysics{FT}(
        param_set;
        ref_state = ref_state,
        turbulence = Vreman{FT}(C_smag),
        moisture = moisture,
        precipitation = precipitation,
        radiation = radiation,
    )

    problem = AtmosProblem(
        boundaryconditions = (
            AtmosBC(
                physics;
                momentum = Impenetrable(DragLaw(
                    (state, aux, t, normPu) -> C_drag,
                )),
                energy = PrescribedEnergyFlux((state, aux, t) -> LHF + SHF),
                moisture = PrescribedMoistureFlux(
                    (state, aux, t) -> moisture_flux,
                ),
            ),
            AtmosBC(physics;),
        ),
        init_state_prognostic = init_dycoms!,
    )

    model = AtmosModel{FT}(
        AtmosLESConfigType,
        physics;
        problem = problem,
        source = source,
    )

    ode_solver = ClimateMachine.ExplicitSolverType(
        solver_method = LSRK144NiegemannDiehlBusch,
    )

    config = ClimateMachine.AtmosLESConfiguration(
        "DYCOMS",
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        param_set,
        init_dycoms!,
        model = model,
    )
    return config, ode_solver
end

function config_diagnostics(driver_config, timeend)
    ssecs = cld(timeend, 2) + 10
    interval = "$(ssecs)ssecs"
    dgngrp = setup_atmos_default_diagnostics(
        AtmosLESConfigType(),
        interval,
        driver_config.name,
    )
    return ClimateMachine.DiagnosticsConfiguration([dgngrp])
end

function main()
    # add a command line argument to specify the kind of
    # moisture and precipitation model you want
    # TODO: this will move to the future namelist functionality
    dycoms_args = ArgParseSettings(autofix_names = true)
    add_arg_group!(dycoms_args, "DYCOMS")
    @add_arg_table! dycoms_args begin
        "--moisture-model"
        help = "specify cloud condensate model"
        metavar = "equilibrium|nonequilibrium"
        arg_type = String
        default = "equilibrium"
        "--precipitation-model"
        help = "specify precipitation model"
        metavar = "noprecipitation|rain"
        arg_type = String
        default = "noprecipitation"
        "--check-asserts"
        help = "should asserts be checked at the end of the simulation"
        metavar = "yes|no"
        arg_type = String
        default = "no"
    end

    cl_args =
        ClimateMachine.init(parse_clargs = true, custom_clargs = dycoms_args)
    moisture_model = cl_args["moisture_model"]
    precipitation_model = cl_args["precipitation_model"]
    check_asserts = cl_args["check_asserts"]

    FT = Float64

    # DG polynomial order
    N = 4

    # Domain resolution and size
    Δh = FT(40)
    Δv = FT(20)
    resolution = (Δh, Δh, Δv)

    xmax = FT(1000)
    ymax = FT(1000)
    zmax = FT(1500)

    t0 = FT(0)
    timeend = FT(100) #FT(4 * 60 * 60)
    Cmax = FT(1.7)     # use this for single-rate explicit LSRK144

    driver_config, ode_solver_type = config_dycoms(
        FT,
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        moisture_model,
        precipitation_model,
    )
    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_solver_type = ode_solver_type,
        init_on_cpu = true,
        Courant_number = Cmax,
    )
    dgn_config = config_diagnostics(driver_config, timeend)

    if moisture_model == "equilibrium"
        filter_vars = ("moisture.ρq_tot",)
    elseif moisture_model == "nonequilibrium"
        filter_vars = ("moisture.ρq_tot", "moisture.ρq_liq", "moisture.ρq_ice")
    end
    if precipitation_model == "rain"
        filter_vars = (filter_vars..., "precipitation.ρq_rai")
    end

    cbtmarfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            filter_vars,
            solver_config.dg.grid,
            TMARFilter(),
        )
        nothing
    end

    result = ClimateMachine.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        user_callbacks = (cbtmarfilter,),
        check_euclidean_distance = true,
    )

    # some simple checks to ensure that rain and clouds exist in the CI runs
    if check_asserts == "yes"

        m = driver_config.bl
        Q = solver_config.Q
        ρ_ind = varsindex(vars_state(m, Prognostic(), FT), :ρ)

        if moisture_model == "nonequilibrium"

            ρq_liq_ind =
                varsindex(vars_state(m, Prognostic(), FT), :moisture, :ρq_liq)
            ρq_ice_ind =
                varsindex(vars_state(m, Prognostic(), FT), :moisture, :ρq_ice)

            min_q_liq = minimum(abs.(
                Array(Q[:, ρq_liq_ind, :]) ./ Array(Q[:, ρ_ind, :]),
            ))
            max_q_liq = maximum(abs.(
                Array(Q[:, ρq_liq_ind, :]) ./ Array(Q[:, ρ_ind, :]),
            ))

            min_q_ice = minimum(abs.(
                Array(Q[:, ρq_ice_ind, :]) ./ Array(Q[:, ρ_ind, :]),
            ))
            max_q_ice = maximum(abs.(
                Array(Q[:, ρq_ice_ind, :]) ./ Array(Q[:, ρ_ind, :]),
            ))

            @info(min_q_liq, max_q_liq)
            @info(min_q_ice, max_q_ice)

            # test that cloud condensate variables exist and are not NaN
            @test !isnan(max_q_liq)
            @test !isnan(max_q_ice)

            # test that there is reasonable amount of cloud water...
            @test abs(max_q_liq) > FT(5e-4)

            # ...and that there is no cloud ice
            @test isequal(min_q_ice, FT(0))
            @test isequal(max_q_ice, FT(0))

        end
        if precipitation_model == "rain"
            ρq_rai_ind = varsindex(
                vars_state(m, Prognostic(), FT),
                :precipitation,
                :ρq_rai,
            )

            min_q_rai = minimum(abs.(
                Array(Q[:, ρq_rai_ind, :]) ./ Array(Q[:, ρ_ind, :]),
            ))
            max_q_rai = maximum(abs.(
                Array(Q[:, ρq_rai_ind, :]) ./ Array(Q[:, ρ_ind, :]),
            ))

            @info(min_q_rai, max_q_rai)

            # test that rain variable exists and is not NaN
            @test !isnan(max_q_rai)

            # test that there is reasonable amount of rain water...
            @test abs(max_q_rai) > FT(1e-6)
        end
    end
end

main()


text:
function hello()
    println("Hello World!")
end

max_stars_repo_name: Devlious/LearningJulia

max_issues_repo_name: Devlious/LearningJulia

max_forks_repo_name: Devlious/LearningJulia


text:
using PastaQ
using Random
using ITensors
using Observers
import Optimisers

Random.seed!(1234)

# 1. Quantum state tomography on measurements generated
# as the output of a unitary quantum circuit

# Generate the sample data
N = 4
depth = 4
nshots = 10_000
circuit = randomcircuit(N, depth)
data, ψ = getsamples(circuit, nshots; local_basis=["X", "Y", "Z"])
writesamples(data, ψ, "data/qst_circuit.h5")

# Load target state and measurement data
data, Ψ = readsamples("data/qst_circuit.h5")

# Split data into a train and test dataset
train_data, test_data = split_dataset(data; train_ratio=0.9)

# Set parameters
N = length(Ψ)     # Number of qubits
χ = maxlinkdim(Ψ) # Bond dimension of variational MPS

# Initialize the variational MPS
ψ0 = randomstate(Ψ; χ=χ)

# Initialize stochastic gradient descent optimizer
opt = Optimisers.Descent(0.01)

# Initialize the observer for the fidelity
F(ψ::MPS; kwargs...) = fidelity(ψ, Ψ)
obs = Observer(["F" => F])

# Run quantum state tomography, where a variational MPS `|ψ(θ)⟩`
# is optimized to mimimize the cross entropy between the data and
# the tensor-network distribution `P(x) = |⟨x|ψ(θ)⟩|²`.
println("Running tomography to learn a pure state ψ:")
ψ = tomography(
  train_data,
  ψ0;
  test_data=test_data,
  optimizer=opt,
  batchsize=1000,
  epochs=5,
  (observer!)=obs,
  print_metrics="F",
)

@show maxlinkdim(ψ)
println()

# 2. Quantum state tomography on measurements generated
# as the output of a noisy quantum circuit

# Generate sample data
data, ρ = getsamples(
  circuit, nshots; local_basis=["X", "Y", "Z"], noise=("amplitude_damping", (γ=0.01,))
)
writesamples(data, ρ, "data/qst_circuit_noisy.h5")

# Load target state and measurement data
data, ϱ = readsamples("data/qst_circuit_noisy.h5")

# Split data into a train and test dataset
train_data, test_data = split_dataset(data; train_ratio=0.9)

# Set parameters
N = length(ϱ)     # Number of qubits
χ = maxlinkdim(ϱ) # Bond dimension of variational LPDO
ξ = 2             # Kraus dimension of variational LPDO

# Initialize the LPDO
ρ0 = randomstate(ϱ; mixed=true, χ=χ, ξ=ξ)

# Initialize stochastic gradient descent optimizer
opt = Optimisers.ADAM()

# Initialize the observer
F(ρ::LPDO; kwargs...) = fidelity(ρ, ϱ)

obs = Observer(["F" => F])

# Run quantum state tomography, where a variational LPDO `ρ(θ)`
# is optimized to mimimize the cross entropy between the data and
# the tensor-network distribution `P(x) = ⟨x|ρ(θ)|x⟩`.
println("Running tomography to learn a mixed state ρ:")
@disable_warn_order begin
  ρ = tomography(
    train_data,
    ρ0;
    test_data=test_data,
    optimizer=opt,
    batchsize=1000,
    epochs=5,
    (observer!)=obs,
    print_metrics="F",
  )
end
@show maxlinkdim(ρ.X)
println()


text:
#cd("/home/david/Documents/gitclones/AlphaZero.jl")
cd("/home/davidyang/Documents/gitclones/AlphaZero.jl")
using Pkg
Pkg.activate(".")
using AlphaZero
#include("/home/david/Documents/gitclones/AlphaZero.jl/games/amazon_full/main.jl")
include("/home/davidyang/Documents/gitclones/AlphaZero.jl/games/amazon_full/main.jl")

gs = Amazon_Full.GameSpec()
g = GI.init(gs)
GI.vectorize_state(gs,GI.current_state(g))
GI.play!(g,g|>GI.actions_mask|>findall|>rand)
GI.vectorize_state(gs,GI.current_state(g))

g = GI.init(gs)
GI.render(g)
while !g.finished
  action = g|>GI.actions_mask|>findall|>rand
  println(GI.action_string(gs,action))
  GI.play!(g,action)
  GI.render(g)
end

cd("/home/davidyang/Documents/gitclones/AlphaZero.jl")
#cd("/home/david/Documents/gitclones/AlphaZero.jl")
#@everywhere using Pkg
#@everywhere Pkg.activate(".")
#@everywhere using AlphaZero
using Pkg
Pkg.activate(".")
using AlphaZero
experiment = Examples.experiments["amazon_full"]
session = Session(experiment, dir="sessions/amazon_full4")
resume!(session)

duel1 = session.benchmark[1]
#duel1 equivalent to eval
progress = UI.Log.Progress(session.logger, duel1.sim.num_games)
env = session.env
params = env.params.self_play
net() = Network.copy(env.bestnn, on_gpu=duel1.sim.use_gpu, test_mode=true)
simulator = Simulator(net, record_trace) do net
  player = Benchmark.instantiate(duel1.player, env.gspec, net)
  baseline = Benchmark.instantiate(duel1.baseline, env.gspec, net)
  return TwoPlayers(player, baseline)
end
#=
samples, elapsed = @timed simulate(
    simulator, env.gspec, duel1.sim,
    game_simulated=(() -> next!(progress)))
=#
oracles = simulator.make_oracles()
spawn_oracles, done =
    AlphaZero.batchify_oracles(oracles; duel1.sim.num_workers, duel1.sim.batch_size, duel1.sim.fill_batches)
oracles = spawn_oracles()
player = simulator.make_player(oracles)

println("debugging play game")
trace = @time play_game(env.gspec, simulator.make_player(oracles), flip_probability=0.0)

colors_flipped = true
report = simulator.measure(trace, colors_flipped, player)

results, elapsed = @timed simulate_distributed(
    simulator, env.gspec, params.sim,
    game_simulated=()->Handlers.game_played(handler))

#seems stuck at the above function, so expand to see what's going on

simulate(simulator,
        gspec,
        SimParams(p; num_games=(w == workers[1] ? num_each + rem : num_each)),
        game_simulated=remote_game_simulated)

using Distributed
gspec = env.gspec
p = params.sim
chan = Distributed.RemoteChannel(()->Channel{Nothing}(1))
Util.@tspawn_main begin
  for i in 1:p.num_games
    take!(chan)
    game_simulated()
  end
end
remote_game_simulated() = put!(chan, nothing)
# Distributing the simulations across workers
num_each, rem = divrem(p.num_games, Distributed.nworkers())
@assert num_each >= 1
workers = Distributed.workers()
tasks = map(workers) do w
  Distributed.@spawnat w begin
    Util.@printing_errors begin
      simulate(
        simulator,
        gspec,
        SimParams(p; num_games=(w == workers[1] ? num_each + rem : num_each)),
        game_simulated=remote_game_simulated)
      end
  end
end
results = fetch.(tasks)
# If one of the worker raised an exception, we print it
for r in results
  if isa(r, Distributed.RemoteException)
    showerror(stderr, r, catch_backtrace())
  end
end

flip_probability = 0.0
player = simulator.make_player(oracles)
game = GI.init(env.gspec)
trace = Trace(GI.current_state(game))
#actions, π_target = think(player, game)
env = player.white.mcts
#MCTS.explore!(env,game,player.white.niters)
#η = MCTS.dirichlet_noise(game, env.noise_α)
#MCTS.run_simulation!(env,GI.clone(game),η=η)
#MCTS.state_info(env,GI.current_state(game))
#env.oracle(state)
#Network.evaluate(env.oracle,state)
nn = env.oracle
gspec = Network.game_spec(nn)
actions_mask = GI.actions_mask(GI.init(gspec, state))
x = GI.vectorize_state(gspec, state)
a = Float32.(actions_mask)
xnet, anet = Network.to_singletons.(Network.convert_input_tuple(nn, (x, a)))
#net_output = Network.forward_normalized(nn, xnet, anet)
p, v = Network.forward(nn, xnet)

while true
  if GI.game_terminated(game)
    return trace
  end
  if !iszero(flip_probability) && rand() < flip_probability
    GI.apply_random_symmetry!(game)
  end
  actions, π_target = think(player, game)
  τ = player_temperature(player, game, length(trace))
  π_sample = apply_temperature(π_target, τ)
  a = actions[Util.rand_categorical(π_sample)]
  GI.play!(game, a)
  push!(trace, π_target, GI.white_reward(game), GI.current_state(game))
end


text: # Stochastic Block model