text:
@testset "Valuation" begin
    @testset "Example 1" begin
        # setup
        @var x
        f = [(x - 10)^5]
        endgame_tracker, starts =
            total_degree(System(f); tracker_options = (automatic_differentiation = 3,))
        S = collect(starts)
        tracker = endgame_tracker.tracker
        val = HC.Valuation(1)
        HC.init!(val)
        track!(tracker, S[1], 1, 1e-13)

        t = tracker.state.t
        HC.update!(val, tracker.predictor, real(t))

        @test val.val_x[1] ≈ 0 atol = 10 * (1e-13)^(1 / 5)
        @test val.val_tẋ[1] ≈ 1 / 5 atol = 10 * (1e-13)^(1 / 5)
    end

    @testset "Example 2" begin
        @var x y
        f = [2.3 * x^2 + 1.2 * y^2 + 3x - 2y + 3, 2.3 * x^2 + 1.2 * y^2 + 5x + 2y - 5]
        endgame_tracker, starts = total_degree(
            System(f);
            gamma = 1.3im + 0.4,
            tracker_options = TrackerOptions(
                parameters = :conservative,
                automatic_differentiation = 3,
            ),
        )
        S = collect(starts)
        tracker = endgame_tracker.tracker
        val = HC.Valuation(2)
        tf = 1e-10
        track!(tracker, S[3], 1, tf)

        t = tracker.state.t
        HC.update!(val, tracker.predictor, real(t))

        @test val.val_x[1] ≈ -1 atol = 10 * tf^(1 / 2)
        @test val.val_x[2] ≈ -1 atol = 10 * tf^(1 / 2)
        @test val.val_tẋ[1] ≈ -1 atol = sqrt(tf)
        @test val.val_tẋ[2] ≈ -1 atol = sqrt(tf)
    end

    @testset "Example 3" begin
        a = [0.257, -0.139, -1.73, -0.199, 1.79, -1.32]
        @var x y
        f1 = (a[1] * x^2 + a[2] * y) * (a[3] * x + a[4] * y) + 1
        f2 = (a[1] * x^2 + a[2] * y) * (a[5] * x + a[6] * y) + 1
        endgame_tracker, starts = total_degree(
            System([f1, f2]);
            gamma = 1.3im + 0.4,
            tracker_options = TrackerOptions(
                parameters = :conservative,
                automatic_differentiation = 3,
            ),
        )
        S = collect(starts)
        tracker = endgame_tracker.tracker
        val = HC.Valuation(2)
        tf = 1e-10
        HC.track!(tracker, S[3], 1, tf)

        t = tracker.state.t
        HC.update!(val, tracker.predictor, real(t))

        @test val.val_x[1] ≈ -1 / 6 atol = tf^(1 / 6)
        @test val.val_x[2] ≈ -2 / 6 atol = tf^(1 / 6)
        @test val.val_tẋ[1] ≈ -1 / 6 atol = tf^(1 / 6)
        @test val.val_tẋ[2] ≈ -2 / 6 atol = tf^(1 / 6)
        @test !isempty(sprint(show, val))
    end
end


text:
using ABM
using Revise
using Distributions
using BenchmarkTools

simstate = SimState()
myschedule = Schedule(simstate)
width = 150.0
height = 150.0
neighborhood_distance = 10.0
jump = 0.7

global field = Field2D(width,height,neighborhood_distance/1.5,true)

addfield!(simstate,field)

function fstep(state::SimState, agent::Agent)
    mypos = agent.state.pos #getObjectLocation(field, agent)
    #@time
   neighborhood = getNeighborsWithinDistance(field, mypos, neighborhood_distance)
  # println("n: ",length(neighborhood))
    #for neighboring in neighborhood
    #    println(agent.state.name," ",agent.id," ",neighboring.id)
    #end
    agent.state.orientation = mypos
    newpos =  Real2D(rand(Uniform(0, width)),rand(Uniform(0, height)))
    setObjectLocation!(field, agent, newpos)
end

mutable struct BoidData
    name::String
    pos::Real2D
    orientation::Real2D
end

for i in 1:10000
    pos = Real2D(rand(Uniform(0, width)), rand(Uniform(0, height)))
    d = BoidData("Boid", pos, Real2D(rand(Uniform(0, width)),rand(Uniform(0, height))))
    boid = Agent(fstep,d)
    setObjectLocation!(field, boid, pos)
    #if i == 1
        scheduleRepeating!(myschedule,boid)
    #end
end

#@time simulate!(myschedule,10)

using Profile
#using ProfileView
using Juno
Profile.init(delay=0.0005, n = 10^7)

Profile.clear()
# @time Profile.@profile simulate!(myschedule,10);
# Juno.profiler()

# @time step!(myschedule)
# println("size ",length(field.fO))
# @time step!(myschedule)
# println("size ",length(field.fO))
# @time step!(myschedule)
# println("size ",length(field.fO))
#

@time  simulate!(myschedule,10);

"""@time while myschedule.steps < 3
    println("[",myschedule.steps,"] time: ",myschedule.time)
    @time step!(myschedule)
    #Swap the fields status to the new one A = B
    @time for field in myschedule.simstate.fields
        swapState!(field)
    end
end"""

"""
function consistency(neighborhood::Vector{Union{Agent,Patch}})
    if length(neighborhood) == 0 return Real2D(0,0) end
    x = 0
    y = 0
    count = 0
    for neighboring in neighborhood
        count++
        x += neighboring.state.orientation.x
        y += neighboring.state.orientation.y
    end
    if count > 0
        x /= count
        y /= count
    end
    return Real2D(x,y)
end

function cohesion(boid::Agent, neighborhood::Vector{Union{Agent,Patch}})
    if length(neighborhood) == 0 return Real2D(0,0) end
    x = 0
    y = 0
    count = 0
    for neighboring in neighborhood
        count++
        tpos = toroidal(boid.state.pos, field.width, field.height)
        ntpos = toroidal(neighboring.state.pos, field.width, field.height)
        x += dx;
                y += dy;
    end
    if count > 0
        x /= count
        y /= count
    end
    return Real2D(-x/10.0,-y/10.0)
end"""

ext: jl


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ b9e63e09-ec0b-4a21-8fea-5f670ea616d8
using ColorTypes

# ╔═╡ 741aafd4-bca2-11eb-0c10-9bc782988f10
md"# Pluto Webcam"

# ╔═╡ 2d57aaa3-e03e-4d2e-9a09-d5f0d8db535c
md"""
- Based on [MIT Introduction to Computational Thinking](https://computationalthinking.mit.edu/Spring21/images/)
- which itself is based on [printi-static](https://github.com/fonsp/printi-static)
"""

# ╔═╡ afb13dfd-6ddb-42bf-bfb9-4e5165b5d11b
function process_raw_camera_data(raw_camera_data)
	if ismissing(raw_camera_data)
		return
	end

	# the raw image data is a long byte array, we need to transform it into something
	# more "Julian" - something with more _structure_.

	# The encoding of the raw byte stream is:
	# every 4 bytes is a single pixel
	# every pixel has 4 values: Red, Green, Blue, Alpha
	# (we ignore alpha for this notebook)

	# So to get the red values for each pixel, we take every 4th value, starting at
	# the 1st:
	reds_flat = UInt8.(raw_camera_data["data"][1:4:end])
	greens_flat = UInt8.(raw_camera_data["data"][2:4:end])
	blues_flat = UInt8.(raw_camera_data["data"][3:4:end])

	# but these are still 1-dimensional arrays, nicknamed 'flat' arrays
	# We will 'reshape' this into 2D arrays:

	width = raw_camera_data["width"]
	height = raw_camera_data["height"]

	# shuffle and flip to get it in the right shape
	reds = reshape(reds_flat, (width, height))' / 255.0
	greens = reshape(greens_flat, (width, height))' / 255.0
	blues = reshape(blues_flat, (width, height))' / 255.0

	# we have our 2D array for each color
	# Let's create a single 2D array, where each value contains the R, G and B value of
	# that pixel

	RGB.(reds, greens, blues)
end

# ╔═╡ cdabe012-4e0e-4ef0-917b-a512537f553b
function camera_input(;max_size=150, default_url="https://i.imgur.com/SUmi94P.png")
"""

			<button id="stop" title="Stop video">✖</button>
			<button id="pop-out" title="Pop out/pop in">⏏</button>

		<video playsinline autoplay></video>

		<button id="shutter" title="Click to take a picture">📷</button>

		Enable webcam

""" |> HTML
end

# ╔═╡ a1f47136-d662-4420-bd08-5a0af9edf474
md"## Usage"

# ╔═╡ b206c86f-546c-4a33-93d2-1324e162ab83
@bind image_bytes camera_input()

# ╔═╡ 5d5cb174-d808-4a51-aedf-d6c64ba92653
my_face = process_raw_camera_data(image_bytes)

# ╔═╡ Cell order:
# ╟─741aafd4-bca2-11eb-0c10-9bc782988f10
# ╟─2d57aaa3-e03e-4d2e-9a09-d5f0d8db535c
# ╠═b9e63e09-ec0b-4a21-8fea-5f670ea616d8
# ╠═afb13dfd-6ddb-42bf-bfb9-4e5165b5d11b
# ╠═cdabe012-4e0e-4ef0-917b-a512537f553b
# ╟─a1f47136-d662-4420-bd08-5a0af9edf474
# ╠═b206c86f-546c-4a33-93d2-1324e162ab83
# ╠═5d5cb174-d808-4a51-aedf-d6c64ba92653


text:
module ApproximateDPs

import FSimZoo: SingleIntegrator
using FSimZoo
using FSimBase

using DifferentialEquations
using ComponentArrays, UnPack
using JLD2, FileIO, DataFrames
using DynamicPolynomials: @polyvar, PolyVar, AbstractPolynomialLike
using Combinatorics: multiexponents
using Transducers
using LinearAlgebra
using NLopt
using NumericalIntegration: integrate

## algorithms
# adp
export CTValueIterationADP
# irl
export CTLinearValueIterationIRL
## utils
export AbstractApproximator, LinearApproximator
export PolynomialBasis
# ## environments
# # integrated environments

include("utils/utils.jl")
include("algorithms/algorithms.jl")
# include("environments/environments.jl")

end


text:
@doc raw"""
	Spinor(α,β)

Spinor(α,β) with Cayley-Klein parameters α and β.

Based on "Introduction to the Shinnar-Le Roux algorithm", Patrick Le Roux (1995).

A spinor is a way to represent 3D rotations, the underlying representation is a 2 X 2 complex unitary matrix (``\alpha,\beta\in\mathbb{C}``):

```math
R=\left[\begin{array}{cc}
\alpha & -\beta^{*}\\
\beta & \alpha^{*}
\end{array}\right],
```
with ``|\alpha|^2+|\beta|^2 = 1``.

This later operates on the ``2\times2`` representation of ``(x,y,z)`` as follows ``V^{+} = R V R^{*}``.
"""
struct Spinor
	α::Complex
	β::Complex
end
Base.show(io::IO,s::Spinor) = print(io, "Spinor(α = ",round(s.α,digits=3),", β = ",round(s.β,digits=3),")")

"""
Spinor multiplication identity.
	(α1,β1)×(α2,β2) = (α1 α2 - β2⋆ β1 , β2 α1 + α2⋆ β1)
"""
*(s1::Spinor,s2::Spinor) = begin
	Spinor(s1.α*s2.α - conj(s2.β)*s1.β,
		   s1.α*s2.β + conj(s2.α)*s1.β)
end

#Rotation matrices
"Spinor clockwise rotation matrix with angle φ with respect to z-axis."
Rz(φ) = Spinor(exp(-im*φ/2), 0)
"Spinor clockwise rotation matrix with angle Θ with respect to y-axis."
Ry(θ) = Spinor(cos(θ/2), sin(θ/2))
"""
Spinor clockwise rotation matrix with angle Θ with respect to x-axis.

Rx(θ) = Rz(-π/2) Ry(θ) Rz(π/2)
"""
Rx(θ) = Spinor(cos(θ/2), -im*sin(θ/2))
"""
Spinor rotation matrix.

Rg(φ1,θ,φ2) = Rz(φ2) Ry(θ) Rz(φ1)
"""
Rg(φ1,θ,φ2) = Spinor(cos(θ/2)*exp(-im*(φ1+φ2)/2), sin(θ/2)*exp(-im*(φ1-φ2)/2))
"""
Spinor rotation matrix with angle θ with axis in the xy plane u=(cosφ,sinφ).

Rφ(φ,θ) = Rg(-φ,θ,φ) = Rz(φ) Ry(θ) Rz(-φ)
"""
Rφ(φ,θ) = Spinor(cos(θ/2),exp(im*φ)*sin(θ/2))

@doc raw"""
Pauly, J., Le Roux, P., Nishimura, D., & Macovski, A. (1991).
Parameter relations for the Shinnar-Le Roux selective excitation pulse design algorithm (NMR imaging).
IEEE Transactions on Medical Imaging, 10(1), 53–65. doi:10.1109/42.75611

Spinor rotation matrix.

Rotation of φ with respect to the axis of rotation n=(nx,ny,nz).

```math
\varphi=-\gamma\Delta t\sqrt{\left|B_{1}\right|^{2}+\left(\boldsymbol{G}\cdot\boldsymbol{x}\right)^{2}}=-\gamma\Delta t\left\Vert \boldsymbol{B}\right\Vert
```
```math
\boldsymbol{n}=\boldsymbol{B}/\left\Vert \boldsymbol{B}\right\Vert
```
"""
Q(φ, nxy, nz) = begin
	Spinor( cos(φ/2)-im*nz*sin(φ/2), -im*nxy*sin(φ/2) )
end
"""
It calculates |\\alpha|^2+|\\beta|^2 of the Cayley-Klein parameters
"""
abs(s::MRIsim.Spinor) = abs(s.α)^2 + abs(s.β)^2

"""RF Object"""
mutable struct RF
	A::Complex # Amplitud B1x + i B1y [T]
	T::Float64 # Duration [s]
end
"Duration `T` [s] of the RF array Array{RF,1}."
dur(x::Array{RF,1}) = sum(x[i].T for i=1:size(x,1))
"Duration `T` [s] of the RF array Array{RF,2}."
dur(x::Array{RF,2}) = maximum(sum([x[i,j].T for i=1:size(x,1),j=1:size(x,2)],dims=2))
"Generate an RF sequence with amplitudes sampled from a function."
RF_fun(f::Function,T::Real,N::Int64=300) = begin
	RFs = [RF(f(t),T/N) for t = range(0,stop=T,length=N)]
end
one(T::Spinor) = Spinor(1,0)
getproperty(x::Vector{RF}, f::Symbol) = getproperty.(x,f)
getproperty(x::Matrix{RF}, f::Symbol) = begin
	if f == :A
		getproperty.(x,:A)
	elseif f == :T
		x[1,:].T
	end
end


text:
using CurveProximityQueries
using Plots
import Random: seed!

pgfplots()
seed!(19)

a = rand(Bernstein{2,12})
cp = a.f.control_points .- 0.5
a = Bernstein(cp..., 2*cp[1]-cp[2], cp[1])

b = Bernstein([[0.95, -0.2], [0.85, -0.2], [0.85, -0.1], [0.9125, -0.1], [0.975, -0.1], [1.0375, -0.1], [1.1, -0.1], [1.1, -0.1875], [1.1, -0.275], [1.1, -0.3625], [1.1, -0.45], [0.975, -0.45], [0.85, -0.45], [0.725, -0.45], [0.6, -0.45], [0.6, -0.3375], [0.6, -0.225], [0.6, -0.1125], [0.6, 0.0], [0.675, 0.0], [0.75, 0.0], [0.825, 0.0], [0.9, 0.0], [0.9, 0.1], [0.8, 0.1], [0.75, 0.1], [0.7, 0.1], [0.7, 0.05], [0.75, 0.0], [0.8, -0.05], [0.9, -0.05]])

pts = closest_points(a, b)

plot(a,
     linecolor=:black,
     linewidth=1.0,
     fill=(0, Gray(0.7)),
    )
plot!(b,
      linewidth=2.0,
      linecolor=:black,
     )
plot!(pts,
      linestyle=:dash,
      markersize=8.0,
      linewidth=2.0,
     )
annotate!(0.07, -0.2, text("\$\\mathcal{B}\$", 25))
annotate!(0.8, -0.2, text("\$\\Psi_{\\mathcal{I}}\$", 25))
annotate!(pts[1][1]+0.05, pts[1][2]+0.05, text("\$b^*\$", 22))
annotate!(pts[2][1]-0.10, pts[2][2]+0.05, text("\$\\psi(t^*)\$", 22))
savefig("out/fig2a.pdf")


text:
@testset "test_connections" begin
## STATE SPACE ##
# CONTINUOUS
C_111 = ss([1], [2], [3], [4])
C_211 = ss(eye_(2), [1; 2], [1 0], [0])
C_212 = ss(eye_(2), [1; 2], eye_(2), [0; 0])
C_221 = ss(eye_(2), [1 0; 0 2], [1 0], [0 0])
C_222 = ss(eye_(2), [1 0; 0 2], eye_(2), zeros(2,2))
C_022 = ss(4*eye_(2))

# DISCRETE
D_111 = ss([1], [2], [3], [4], 0.005)
D_211 = ss(eye_(2), [1; 2], [1 0], [0], 0.005)
D_212 = ss(eye_(2), [1;2], eye_(2), [0; 0], 0.005)
D_221 = ss(eye_(2), [1 0; 0 2], [1 0], [0 0], 0.005)
D_222 = ss(eye_(2), [1 0; 0 2], eye_(2), zeros(2,2), 0.005)
D_022 = ss(4*eye_(2), 0.005)

@test [C_111 C_221] == ss(eye_(3), [2 0 0; 0 1 0; 0 0 2], [3 1 0], [4 0 0])
@test [C_111; C_212] == ss(eye_(3), [2; 1; 2], [3 0 0; 0 1 0; 0 0 1], [4; 0; 0])
@test append(C_111, C_211) == ss(eye_(3), [2 0; 0 1; 0 2], [3 0 0; 0 1 0], [4 0; 0 0])
@test [C_022 C_222] == ss(eye_(2), [0 0 1 0; 0 0 0 2], [1 0; 0 1], [4 0 0 0; 0 4 0 0])
@test [C_022; C_222] == ss(eye_(2), [1 0; 0 2], [0 0; 0 0; 1 0; 0 1], [4 0; 0 4; 0 0; 0 0])

@test [D_111 D_221] == ss(eye_(3), [2 0 0; 0 1 0; 0 0 2], [3 1 0], [4 0 0], 0.005)
@test [D_111; D_212] == ss(eye_(3), [2; 1; 2], [3 0 0; 0 1 0; 0 0 1], [4; 0; 0], 0.005)
@test append(D_111, D_211) == ss(eye_(3), [2 0; 0 1; 0 2], [3 0 0; 0 1 0], [4 0; 0 0], 0.005)
@test [D_022 D_222] == ss(eye_(2), [0 0 1 0; 0 0 0 2], [1 0; 0 1], [4 0 0 0; 0 4 0 0], 0.005)
@test [D_022; D_222] == ss(eye_(2), [1 0; 0 2], [0 0; 0 0; 1 0; 0 1], [4 0; 0 4; 0 0; 0 0], 0.005)

@test series(C_111, C_212) == C_212*C_111
@test parallel(C_111, C_211) == C_111 + C_211

# Errors
@test_throws ErrorException [C_111 D_111]                 # Sampling time mismatch
@test_throws ErrorException [C_111; D_111]                # Sampling time mismatch
@test_throws ErrorException append(C_111, D_111)          # Sampling time mismatch
@test_throws ErrorException [C_111 C_212]                 # Dimension mismatch
@test_throws ErrorException [C_111; C_221]                # Dimension mismatch

## TRANSFER FUNCTION ##
# CONTINUOUS
Ctf_111 = tf([1, 2], [1, 5])
Ctf_211 = tf([1, 2, 3], [1, 8, 15])
Ctf_212 = tf(vecarray(2, 1, [1, 2, 3], [1, 2]), vecarray(2, 1, [1, 8, 15], [1, 8, 15]))
Ctf_221 = tf(vecarray(1, 2, [1, 2, 3], [1, 2]), vecarray(1, 2, [1, 8, 15], [1, 8, 15]))
Ctf_222 = [Ctf_221; Ctf_221]
Ctf_022 = tf(4*eye_(2))

# DISCRETE
Dtf_111 = tf([1, 2], [1, 5], 0.005)
Dtf_211 = tf([1, 2, 3], [1, 8, 15], 0.005)
Dtf_212 = tf(vecarray(2, 1, [1, 2, 3], [1, 2]), vecarray(2, 1, [1, 8, 15], [1, 8, 15]), 0.005)
Dtf_221 = tf(vecarray(1, 2, [1, 2, 3], [1, 2]), vecarray(1, 2, [1, 8, 15], [1, 8, 15]), 0.005)
Dtf_222 = [Dtf_221; Dtf_221];
Dtf_022 = tf(4*eye_(2), 0.005)

s = tf("s")
@test [Ctf_111 Ctf_221] == tf(vecarray(1, 3, [1,2], [1,2,3], [0,1,2]),
    vecarray(1, 3, [1,5], [1,8,15], [1,8,15]))
@test [Ctf_111; Ctf_212] == tf(vecarray(3, 1, [1,2], [1,2,3], [0,1,2]),
    vecarray(3, 1, [1,5], [1,8,15], [1,8,15]))
@test append(Ctf_111, Ctf_211) == tf(vecarray(2, 2, [1,2], [0], [0], [1,2,3]),
    vecarray(2, 2, [1,5], [1], [1], [1,8,15]));
@test [Ctf_022 Ctf_222] == tf(vecarray(2, 4, [4], [0], [1,2,3], [0,1,2], [0], [4], [1,2,3], [0,1,2]),
    vecarray(2, 4, [1], [1], [1,8,15], [1,8,15], [1], [1], [1,8,15], [1,8,15]))
@test [Ctf_022; Ctf_222] == tf(vecarray(4, 2, [4], [0], [0], [4], [1,2,3], [0,1,2], [1,2,3], [0,1,2]),
    vecarray(4, 2, [1], [1], [1], [1], [1,8,15], [1,8,15], [1,8,15], [1,8,15]))
@test [Ctf_022 Ctf_022] == [[tf(4) 0;0 4] 4*eye_(2)]

@test [Dtf_111 Dtf_221] == tf(vecarray(1, 3, [1,2], [1,2,3], [0,1,2]),
    vecarray(1, 3, [1,5], [1,8,15], [1,8,15]), 0.005)
@test [Dtf_111; Dtf_212] == tf(vecarray(3, 1, [1,2], [1,2,3], [0,1,2]),
    vecarray(3, 1, [1,5], [1,8,15], [1,8,15]), 0.005)
@test append(Dtf_111, Dtf_211) == tf(vecarray(2, 2, [1,2], [0], [0], [1,2,3]),
    vecarray(2, 2, [1,5], [1], [1], [1,8,15]), 0.005);
@test [Dtf_022 Dtf_222] == tf(vecarray(2, 4, [4], [0], [1,2,3], [0,1,2], [0], [4], [1,2,3], [0,1,2]),
    vecarray(2, 4, [1], [1], [1,8,15], [1,8,15], [1], [1], [1,8,15], [1,8,15]), 0.005)
@test [Dtf_022; Dtf_222] == tf(vecarray(4, 2, [4], [0], [0], [4], [1,2,3], [0,1,2], [1,2,3], [0,1,2]),
    vecarray(4, 2, [1], [1], [1], [1], [1,8,15], [1,8,15], [1,8,15], [1,8,15]), 0.005)

@test series(Ctf_111, Ctf_212) == tf(vecarray(2, 1, [1,4,7,6], [0,1,4,4]),
    vecarray(2, 1, [1,13,55,75], [1,13,55,75]));
@test parallel(Ctf_111, Ctf_211) == tf([2,17,44,45], [1,13,55,75])

# Combination tf and ss
@test [C_111 Ctf_221] == [C_111 ss(Ctf_221)]
@test [C_111; Ctf_212] == [C_111; ss(Ctf_212)]
@test append(C_111, Ctf_211) == append(C_111, ss(Ctf_211))
@test [D_111 Dtf_221] == [D_111 ss(Dtf_221)]
@test [D_111; Dtf_212] == [D_111; ss(Dtf_212)]
@test append(D_111, Dtf_211) == append(D_111, ss(Dtf_211))

# hcat and vcat for StateSpace and Matrix
A = [-1.1 -1.2; -1.3 -1.4]
B = [1 2; 3 4]
C = [5 6; 7 8]
D = [1 0; 0 1]
P = ss(A, B, C, D)
@test [P fill(2.5, 2, 1)] == ss(A, [B fill(0, 2, 1)], C, [D fill(2.5, 2, 1)])
@test [fill(2.5, 2, 1) P] == ss(A, [fill(0, 2, 1) B], C, [fill(2.5, 2, 1) D])
@test [P; fill(2.5, 1, 2)] == ss(A, B, [C; fill(0, 1, 2)], [D; fill(2.5, 1, 2)])
@test [fill(2.5, 1, 2); P] == ss(A, B, [fill(0, 1, 2); C], [fill(2.5, 1, 2); D])

# hcat and vcat for StateSpace and Number
P = ss(-1.0, 2.0, 3.0, 4.0)
@test [P 2.5] == ss(-1.0, [2.0 0.0], 3.0, [4.0 2.5])
@test [2.5 P] == ss(-1.0, [0.0 2.0], 3.0, [2.5 4.0])
@test [P; 2.5] == ss(-1.0, 2.0, [3.0; 0.0], [4.0; 2.5])
@test [2.5; P] == ss(-1.0, 2.0, [0.0; 3.0], [2.5; 4.0])

@test [2.5 P 3.5] == ss(-1.0, [0.0 2.0 0.0], 3.0, [2.5 4.0 3.5])
@test [2.5; P; 3.5] == ss(-1.0, 2.0, [0.0; 3.0; 0.0], [2.5; 4.0; 3.5])

# Combination tfRational and sisoZpk
Czpk_111 = zpk([-2],[-5],1)
Czpk_211 = zpk([-1+sqrt(2)im,-1-sqrt(2)im], [-5,-3],1)
Czpk_212 = zpk(vecarray(2, 1, [-1+sqrt(2)im,-1-sqrt(2)im], [-2]), vecarray(2, 1, [-5,-3], [-5,-3]), fill(1, 2, 1))
Czpk_221 = zpk(vecarray(1, 2, [-1+sqrt(2)im,-1-sqrt(2)im], [-2]), vecarray(1, 2, [-5,-3], [-5,-3]), fill(1, 1, 2))
Czpk_222 = [Czpk_221; Czpk_221]
Czpk_022 = [zpk([],[],4) 0; 0 zpk([],[],4)]

#Make sure that we get a vector
arr = Array{typeof(zpk(tf(1))),1}(undef, 2)
arr[1] = zpk(tf(1)); arr[2] = zpk(2);
@test [tf(1), zpk(2)] == arr
arr2 = Array{typeof(tf(1, 0.1)),1}(undef, 2)
arr2[1] = tf(1, 0.1); arr2[2] = tf(2, 0.1);
@test [tf(1, 0.1), tf(2, 0.1)] == arr2
arr3 = Array{typeof(ss(0.)),1}(undef, 3)
arr3[1] = ss(0.); arr3[2] = ss(1.); arr3[3] = ss(2.)
@test [0., zpk(1), ss(2.)] == arr3

arr4 = Array{typeof(ss(0)),1}(undef, 3)
arr4[1] = ss(0); arr4[2] = ss(1); arr4[3] = ss(2)
@test [0., zpk(1), ss(2)] == arr4

@test Czpk_111 ≈ Ctf_111
@test Czpk_211 ≈ Ctf_211
@test Czpk_212 ≈ Ctf_212
@test Czpk_221 ≈ Ctf_221
@test Czpk_222 ≈ Ctf_222
@test Czpk_022 ≈ Ctf_022

@test Czpk_222 ≈ [Ctf_221; Czpk_221]

#This might fail depending on if minreal is used or not
@test (Czpk_211+1) ≈ (Ctf_211+1)

@test_broken [D_111 1.0] # Concatenation of discrete system with constant
@test_broken [D_222 fill(1.5, 2, 2)] # Concatenation of discrete system with matrix

end


text:
println("loading packages... ")

using ArgParse
using OrthographicProjection

println("packages OK")

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
		"source"
			help = "A text file with Potree directories list or a single Potree directory"
			required = true
		"--output", "-o"
            help = "Output image"
			required = true
		"--bbin"
            help = "Bounding box as 'x_min y_min z_min x_max y_max z_max' or Potree JSON volume model"
			required = true
		"--po"
            help = "Projection plane: XY+, XY-, XZ+, XZ-, YZ+, YZ-"
			arg_type = String
            default = "XY+"
        "--gsd"
            help = "Resolution"
            arg_type = Float64
            default = 0.3
		"--quote"
			help = "Distance of plane from origin"
			arg_type = Float64
		"--thickness"
			help = "Section thickness"
			arg_type = Float64
		"--pc"
			help = "If true a pc of extracted model is saved in a LAS file"
			action = :store_true
		"--ucs"
			help = "Path to UCS JSON file. If not provided is the Identity matrix"
			arg_type = String
		"--bgcolor"
			help = "Background color"
			arg_type = String
		"--epsg"
			help = "EPSG code"
			arg_type = Int
    end

    return parse_args(s)
end

function main()
    args = parse_commandline()

	txtpotreedirs = args["source"]
	outputimage = args["output"]
	bbin = args["bbin"]
	PO = args["po"]
	GSD = args["gsd"]
	bgcolor = args["bgcolor"]
	altitude = args["quote"]
	thickness = args["thickness"]
	pc = args["pc"]
	ucs = args["ucs"]
	epsg = args["epsg"]

	OrthographicProjection.flushprintln(" ")
	OrthographicProjection.flushprintln("== params ==")
	OrthographicProjection.flushprintln("Sources  =>  $txtpotreedirs")
	OrthographicProjection.flushprintln("Output image  =>  $outputimage")

	b = tryparse.(Float64,split(bbin, " "))
	if length(b) == 6
		#bbin = (hcat([b[1],b[2],b[3]]),hcat([b[4],b[5],b[6]]))
		bbin = Common.AABB(b[4],b[1],b[5],b[2],b[6],b[3])
	end

	OrthographicProjection.flushprintln("Bounding Box  =>")
	OrthographicProjection.flushprintln(bbin)
	OrthographicProjection.flushprintln("Point of View  =>  $PO")
	OrthographicProjection.flushprintln("GSD  =>  $GSD")

	if !isnothing(bgcolor)
		col = tryparse.(Float64,split(bgcolor, " "))
		background = [col[1],col[2],col[3]]
	else
		background = [1.0,1.0,1.0]
	end

	OrthographicProjection.flushprintln("Background color  =>  $background")

	if !isnothing(altitude)
		OrthographicProjection.flushprintln("Altitude  =>  $altitude")
		OrthographicProjection.flushprintln("Thickness  =>  $thickness")
	end

	OrthographicProjection.flushprintln("Extract point cloud  =>  $pc")

	if isnothing(ucs)
		ucs = Matrix{Float64}(Common.I,4,4)
	else
		OrthographicProjection.flushprintln("User coordinates system  =>  $ucs")
	end

	if !isnothing(epsg)
		OrthographicProjection.flushprintln("EPSG  =>  $epsg")
	end

	OrthographicProjection.orthophoto(txtpotreedirs, outputimage, bbin, GSD, PO, altitude, thickness, ucs, background, pc, epsg)
end

@time main()


text:
using GR

include("noise.jl")

const n = 100000
const radius = 200
const width, height = 500, 500
const length = 50
const scale = 0.005

function main()

    T = rand(n) * 2pi
    R = sqrt.(rand(n))
    X = R .* cos.(T)
    Y = R .* sin.(T)
    intensity = (1.001 .- sqrt.(X.^2 .+ Y.^2)) .^ 0.75
    X = X .* radius .+ div(width, 2)
    Y = Y .* radius .+ div(height, 2)

    w, h, I = readimage("julia-logo.png")
    P = zeros(UInt32, width, height)

    settextfontprec(128, 0)
    settextalign(2, 3)
    setcharheight(0.036)

    for time in 0:0.004:1
        P .= 0xffffffff
        cos_t = 1.5 * cos(2pi * time)
        sin_t = 1.5 * sin(2pi * time)
        for i in 1:n
            x, y = X[i], Y[i]
            dx = noise4(scale * x, scale * y, cos_t, sin_t, 2)
            dx *= intensity[i] * length
            dy = noise4(100 + scale * x, 200 + scale * y, cos_t, sin_t, 2)
            dy *= intensity[i] * length
            P[round(Int, x + dx), round(Int, y + dy)] = I[round(Int, x), round(Int, y)]
        end

        clearws()
        setviewport(0, 1, 0, 1)
        drawimage(0, 1, 0, 1, width, height, P)
        textext(0.5, -0.2 + 1.5*time, "Julia Noise\\n\\nMade with GR.jl\\n\\nOriginal idea by Necessary Disorder")
        updatews()
    end
end

main()


text:
module LiftedMaps

export LiftedMap, isondiagonal, isblockdiagonal

using LinearMaps
using LinearAlgebra
using BlockArrays

struct LiftedMap{T,TI,TJ,TM,TN} <: LinearMap{T}
    A::LinearMap{T}
    I::TI
    J::TJ
    M::TM
    N::TN
end

LiftedMap(A::AbstractMatrix, I, J, M, N) = LiftedMap(LinearMap(A), I, J, M, N)
LinearMaps.MulStyle(A::LiftedMap) = LinearMaps.FiveArg()

Base.size(A::LiftedMap) = (length(A.M), length(A.N))
Base.axes(A::LiftedMap) = (A.M, A.N)
Base.axes(A::LiftedMap, i::Int) = axes(A)[i]

LinearAlgebra.adjoint(A::LiftedMap) = LiftedMap(adjoint(A.A), A.J, A.I, A.N, A.M)
LinearAlgebra.transpose(A::LiftedMap) = LiftedMap(transpose(A.A), A.J, A.I, A.N, A.M)

function LinearAlgebra.mul!(y::AbstractVector, L::LiftedMap,
    x::AbstractVector, α::Number, β::Number)

    bvy = PseudoBlockVector(y, blocksizes(L.M)...)
    bvx = PseudoBlockVector(x, blocksizes(L.N)...)

    yI = view(bvy, L.I)
    xJ = view(bvx, L.J)
    AIJ = L.A

    y .*= β
    LinearAlgebra.mul!(yI, AIJ, xJ, α, 1)
    return y
end

# function LinearAlgebra.mul!(y::AbstractVector,
#     Lt::LinearMaps.TransposeMap{<:Any,<:LiftedMap},
#     x::AbstractVector, α::Number, β::Number)

#     L = Lt.lmap

#     bvy = PseudoBlockVector(y, blocksizes(L.N)...)
#     bvx = PseudoBlockVector(x, blocksizes(L.M)...)

#     yJ = view(bvy, L.J)
#     xI = view(bvx, L.I)
#     AIJ = L.A

#     y .*= β
#     LinearAlgebra.mul!(yJ, transpose(AIJ), xI, α, 1)
#     return y
# end

function LinearAlgebra.mul!(y::AbstractVector, L::LiftedMap, x::AbstractVector)

    bvy = PseudoBlockVector(y, blocksizes(L.M)...)
    bvx = PseudoBlockVector(x, blocksizes(L.N)...)

    yI = view(bvy, L.I)
    xJ = view(bvx, L.J)
    AIJ = L.A

    fill!(y,0)
    LinearAlgebra.mul!(yI, AIJ, xJ)
    return y
end

# function LinearAlgebra.mul!(y::AbstractVector, Lt::LinearMaps.TransposeMap{<:Any,<:LiftedMap}, x::AbstractVector)

#     L = Lt.lmap

#     bvy = PseudoBlockVector(y, blocksizes(L.N)...)
#     bvx = PseudoBlockVector(x, blocksizes(L.M)...)

#     yJ = view(bvy, L.J)
#     xI = view(bvx, L.I)
#     AIJ = L.A

#     fill!(y,0)
#     LinearAlgebra.mul!(yJ, transpose(AIJ), xI)
#     return y
# end

function Base.:(*)(A::LiftedMap, x::AbstractVector)
    axes(A,2) == axes(x,1) ||
        throw(DimensionMismatch("second dimension of left factor, $(size(A, 2)), " *
                "does not match first dimension of right factor, $(size(x, 1))"))

    T = promote_type(eltype(A), eltype(x))
    y = PseudoBlockVector{T}(undef, BlockArrays.blocksizes(A,1))
    fill!(y, zero(T))
    LinearAlgebra.mul!(y, A, x)
end

function Base.Matrix{T}(A::LiftedMap) where {T}
    # M = zeros(T, size(A))
    # M = PseudoBlockMatrix{T}(undef, BlockArrays.blocksizes(A)...)
    M = similar(Array{T}, axes(A))
    fill!(M,0)
    m = Matrix(A.A)
    M[A.I,A.J] .= m
    return Matrix{T}(M)
end

function Base.Matrix{T}(A::LinearMaps.ScaledMap{<:Any, <:Any, <:LiftedMap}) where {T}
    M = Matrix{T}(A.lmap)
    M .*= A.λ
    return M
end

function isondiagonal(A::LiftedMaps.LiftedMap)
    return A.I == A.J
end

function isondiagonal(A::LinearMaps.ScaledMap)
    return isondiagonal(A.lmap)
end

function isondiagonal(A::LinearMaps.LinearMap)
    return false
end

function isblockdiagonal(A::LinearMaps.LinearCombination)
    for map in A.maps
        isondiagonal(map) || return false
    end
    return true
end

function isblockdiagonal(A::LinearMaps.LinearMap)
    return false
end

end


text:
@doc """Result of a capture
  Result.val has form:
  Processor id -> (lensname -> (varname -> Vector of values captured at lens)""" ->
immutable Result
  values::Dict{Int,Dict{Symbol,Dict{Symbol,Vector{Any}}}}
end

Result() = Result(Dict{Int,Dict{Symbol,Vector{Any}}}())
convert(::Type{Vector{Result}}, x::Vector{Any}) =
  (rs = similar(x,Result); for i = 1:length(x) rs[i] = x[i] end)

# Convenience functions for extracting data from a Result
function get(r::Result; proc_id::Int=1, lensname=nothing, capturename=nothing)
  entries = r.values[proc_id]
  if lensname == nothing
    length(entries) != 1 && error("No lensname specified and more than one lens captured")
    lensname = first(entries)[1]
  end
  if capturename == nothing
    length(entries[lensname]) != 1 && error("No capture name specified and more than one captured <found></found>")
    capturename = first(entries[lensname])[1]
  end
  entries[lensname][capturename]
end

get{T}(@compat r::Tuple{T,Result}) = get(r[2])


text:
####################
## old code below ##
####################

function subsampleProducer(nObs, minObs; rolling=false)
    ## produce indices for changing subsamples over time
    ##
    ## Inputs:
    ## 	nObs			scalar value denoting overall number of
    ## 					observations
    ## 	minObs 		scalar value denoting minimum subsample size
    ## 	rolling		optional - boolean value, indicating whether
    ## 					subsample is each day or rolled over. Default is
    ## 					increasing subsamples.
    ##
    ## Outputs:
    ## 	(begInd, endInd) tuple of start and end indices of subsample,
    ## 	with last subsample ending at second last index, since
    ## 	estimated values of given subsample can only be used for the
    ## 	strategy at the next day

    for ii=minObs:(nObs-1)
        if rolling
            produce((ii-minObs+1, ii))
        else
            produce((1, ii))
        end
    end
end


text:
@def F_YF begin #(V, U, Ax, R, SA)		# fRont tire force

  if length(Ax) != 1
    Fyf = zeros(Float64, (length(Ax),1))
    for ii in eachindex(Ax)
      Fyf[ii] = ((PD2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff))*sin(PC1*atan((((PK1*sin(2*atan(PK2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff))))/(((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)) + ((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)))/(((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff))^2 + EP^2)^(0.5))*0.001)+EP))*((atan((V[ii] + la*R[ii])/(U[ii]+EP)) - SA[ii]) + PH2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff) + PH1)) - ((PE2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff) + PE1)*(1 - PE3)*(((atan((V[ii] + la*R[ii])/(U[ii]+EP)) - SA[ii]) + PH2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff) + PH1))/((((atan((V[ii] + la*R[ii])/(U[ii]+EP)) - SA[ii]) + PH2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff) + PH1)^2 + EP^2)^(0.5)))*((((PK1*sin(2*atan(PK2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff))))/(((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)) + ((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)))/(((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff))^2 + EP^2)^(0.5))*0.001)+EP))*((atan((V[ii] + la*R[ii])/(U[ii]+EP)) - SA[ii]) + PH2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff) + PH1)) - atan((((PK1*sin(2*atan(PK2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff))))/(((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)) + ((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)))/(((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff))^2 + EP^2)^(0.5))*0.001)+EP))*((atan((V[ii] + la*R[ii])/(U[ii]+EP)) - SA[ii]) + PH2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff) + PH1)))))) + (PV2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PV1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)));
    end
  else
    Fyf = Float64
    Fyf = ((PD2*(FzF0 - (Ax - V*R)*dFzx_coeff)^2 + PD1*(FzF0 - (Ax - V*R)*dFzx_coeff))*sin(PC1*atan((((PK1*sin(2*atan(PK2*(FzF0 - (Ax - V*R)*dFzx_coeff))))/(((PD2*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff)) + ((PD2*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff)))/(((PD2*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff))^2 + EP^2)^(0.5))*0.001)+EP))*((atan((V + la*R)/(U+EP)) - SA) + PH2*(FzF0 - (Ax - V*R)*dFzx_coeff) + PH1)) - ((PE2*(FzF0 - (Ax - V*R)*dFzx_coeff) + PE1)*(1 - PE3)*(((atan((V + la*R)/(U+EP)) - SA) + PH2*(FzF0 - (Ax - V*R)*dFzx_coeff) + PH1))/((((atan((V + la*R)/(U+EP)) - SA) + PH2*(FzF0 - (Ax - V*R)*dFzx_coeff) + PH1)^2 + EP^2)^(0.5)))*((((PK1*sin(2*atan(PK2*(FzF0 - (Ax - V*R)*dFzx_coeff))))/(((PD2*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff)) + ((PD2*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff)))/(((PD2*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff))^2 + EP^2)^(0.5))*0.001)+EP))*((atan((V + la*R)/(U+EP)) - SA) + PH2*(FzF0 - (Ax - V*R)*dFzx_coeff) + PH1)) - atan((((PK1*sin(2*atan(PK2*(FzF0 - (Ax - V*R)*dFzx_coeff))))/(((PD2*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff)) + ((PD2*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff)))/(((PD2*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*dFzx_coeff))^2 + EP^2)^(0.5))*0.001)+EP))*((atan((V + la*R)/(U+EP)) - SA) + PH2*(FzF0 - (Ax - V*R)*dFzx_coeff) + PH1)))))) + (PV2*(FzF0 - (Ax - V*R)*dFzx_coeff)^2 + PV1*(FzF0 - (Ax - V*R)*dFzx_coeff)));
  end
  Fyf
end


text:
@testset "c\"enum ...\"" begin
	@eval c``

	@eval c"""
	enum OpaqueEnum;
	"""
	@eval c"""
	enum OpaqueEnum;
	"""
	@test_throws ErrorException sizeof(c"enum OpaqueEnum")
	@test sizeof(c"enum OpaqueEnum *") == sizeof(Ptr{Cvoid})

	@eval c"""
	enum SimpleEnum;
	"""
	@eval c"""
	enum SimpleEnum {
		SE1,
		SE2,
		SE3,
	};
	"""
	@test sizeof(c"enum SimpleEnum") == 4
	@test c"SE1" == 0
	@test c"SE2" == 1
	@test c"SE3" == 2

	@eval c"""
	enum LessSimpleEnum {
		LSE1,
		LSE2,
		LSE3 = -1,
	};
	"""
	@test sizeof(c"enum LessSimpleEnum") == 4
	@test c"LSE1" == 0
	@test c"LSE2" == 1
	@test c"LSE3" == -1

	@eval c"""
	enum ComplexEnum {
		CE1 = 1,
		CE2 = CE1,
		CE3 = 0,
		CE4,
	};
	"""
	@test sizeof(c"enum ComplexEnum") == 4
	@test c"CE1" == 1
	@test c"CE2" == 1
	@test c"CE3" == 0
	@test c"CE4" == 1

	@eval c"""
	enum TrickyEnum {
		TE1 = -1,
		TE2 = 0xffffffff,
	};
	"""
	@test sizeof(c"enum TrickyEnum") == 8
	@test c"TE1" == -1
	@test c"TE2" == 0xffffffff

	@eval c"""
	enum PackedSimpleEnum {
		PSE1,
		PSE2,
		PSE3,
	} __attribute__((packed));
	"""
	@test sizeof(c"enum PackedSimpleEnum") == 1
	@test c"PSE1" == 0
	@test c"PSE2" == 1
	@test c"PSE3" == 2

	@eval c"""
	enum PackedLessSimpleEnum {
		PLSE1,
		PLSE2,
		PLSE3 = -1,
	} __attribute__((packed));
	"""
	@test sizeof(c"enum PackedLessSimpleEnum") == 1
	@test c"PLSE1" == 0
	@test c"PLSE2" == 1
	@test c"PLSE3" == -1

	@eval c"""
	enum PackedComplexEnum {
		PCE1 = 1,
		PCE2 = PCE1,
		PCE3 = 0,
		PCE4,
	} __attribute__((packed));
	"""
	@test sizeof(c"enum PackedComplexEnum") == 1
	@test c"PCE1" == 1
	@test c"PCE2" == 1
	@test c"PCE3" == 0
	@test c"PCE4" == 1

	@eval c"""
	enum PackedTrickyEnum {
		PTE1 = -1,
		PTE2 = 0xffffffff,
	} __attribute__((packed));
	"""
	@test sizeof(c"enum PackedTrickyEnum") == 8
	@test c"PTE1" == -1
	@test c"PTE2" == 0xffffffff
end


text:
#!/usr/bin/env julia --project
using ClimateMachine
ClimateMachine.init(parse_clargs = true)

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

using Distributions
using StaticArrays
using Test
using DocStringExtensions
using LinearAlgebra

using CLIMAParameters
using CLIMAParameters.Planet: R_d, cp_d, cv_d, MSLP, grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

# -------------------- Surface Driven Bubble ----------------- #
# Rising thermals driven by a prescribed surface heat flux.
# 1) Boundary Conditions:
#       Laterally periodic with no flow penetration through top
#       and bottom wall boundaries.
#       Momentum: Impenetrable(FreeSlip())
#       Energy:   Spatially varying non-zero heat flux up to time t₁
# 2) Domain: 1250m × 1250m × 1000m
# Configuration defaults are in `src/Driver/Configurations.jl`

"""
  Surface Driven Thermal Bubble
"""
function init_surfacebubble!(problem, bl, state, aux, localgeo, t)
    (x, y, z) = localgeo.coord
    param_set = parameter_set(bl)
    FT = eltype(state)
    R_gas::FT = R_d(param_set)
    c_p::FT = cp_d(param_set)
    c_v::FT = cv_d(param_set)
    p0::FT = MSLP(param_set)
    _grav::FT = grav(param_set)
    γ::FT = c_p / c_v

    xc::FT = 1250
    yc::FT = 1250
    zc::FT = 1250
    θ_ref::FT = 300
    Δθ::FT = 0

    #Perturbed state:
    θ = θ_ref + Δθ # potential temperature
    π_exner = FT(1) - _grav / (c_p * θ) * z # exner pressure
    ρ = p0 / (R_gas * θ) * (π_exner)^(c_v / R_gas) # density

    q_tot = FT(0)
    ts = PhaseEquil_ρθq(param_set, ρ, θ, q_tot)
    q_pt = PhasePartition(ts)

    ρu = SVector(FT(0), FT(0), FT(0))
    # energy definitions
    e_kin = FT(0)
    e_pot = gravitational_potential(bl.orientation, aux)
    ρe_tot = ρ * total_energy(e_kin, e_pot, ts)
    state.ρ = ρ
    state.ρu = ρu
    state.energy.ρe = ρe_tot
    state.moisture.ρq_tot = ρ * q_pt.tot
end

function config_surfacebubble(FT, N, resolution, xmax, ymax, zmax)
    # Boundary conditions
    # Heat Flux Peak Magnitude
    F₀ = FT(100)
    # Time [s] at which `heater` turns off
    t₁ = FT(500)
    # Plume wavelength scaling
    x₀ = xmax
    function energyflux(state, aux, t)
        x = aux.coord[1]
        y = aux.coord[2]
        MSEF = F₀ * (cospi(2 * x / x₀))^2 * (cospi(2 * y / x₀))^2
        t < t₁ ? MSEF : zero(MSEF)
    end

    C_smag = FT(0.23)

    physics = AtmosPhysics{FT}(
        param_set;
        turbulence = SmagorinskyLilly{FT}(C_smag),
        moisture = EquilMoist(),
    )

    problem = AtmosProblem(
        boundaryconditions = (
            AtmosBC(physics; energy = PrescribedEnergyFlux(energyflux)),
            AtmosBC(physics;),
        ),
        init_state_prognostic = init_surfacebubble!,
    )
    model = AtmosModel{FT}(
        AtmosLESConfigType,
        physics;
        problem = problem,
        source = (Gravity(),),
    )
    config = ClimateMachine.AtmosLESConfiguration(
        "SurfaceDrivenBubble",
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        param_set,
        init_surfacebubble!,
        model = model,
    )

    return config
end

function config_diagnostics(
    driver_config,
    xmax::FT,
    ymax::FT,
    zmax::FT,
    resolution,
    interval,
) where {FT}
    dgngrp = StdDiagnostics.AtmosLESDefault(interval, driver_config.name)
    boundaries = [
        FT(0) FT(0) FT(0)
        xmax ymax zmax
    ]
    interpol = ClimateMachine.InterpolationConfiguration(
        driver_config,
        boundaries,
        resolution,
    )
    pdgngrp = setup_atmos_default_perturbations(
        AtmosLESConfigType(),
        interval,
        driver_config.name,
        interpol = interpol,
    )
    return ClimateMachine.DiagnosticsConfiguration([dgngrp, pdgngrp])
end

function main()
    FT = Float64
    # DG polynomial order
    N = 4
    # Domain resolution and size
    Δh = FT(50)
    Δv = FT(50)
    resolution = (Δh, Δh, Δv)
    xmax = FT(2000)
    ymax = FT(2000)
    zmax = FT(2000)
    t0 = FT(0)
    timeend = FT(2000)

    driver_config = config_surfacebubble(FT, N, resolution, xmax, ymax, zmax)

    ode_solver_type = ClimateMachine.ExplicitSolverType(
        solver_method = LSRK144NiegemannDiehlBusch,
    )

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        init_on_cpu = true,
        ode_solver_type = ode_solver_type,
    )

    dgn_ssecs = (timeend / 2) + 10
    dgn_interval = "$(dgn_ssecs)ssecs"
    dgn_config = config_diagnostics(
        driver_config,
        xmax,
        ymax,
        zmax,
        resolution,
        dgn_interval,
    )

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

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

    @test isapprox(result, FT(1); atol = 1.5e-3)
end

main()


text:
"""
    set path to SNOPT MathOptInterface interface.
        suggested: https://github.com/thowell/SNOPT7.jl
            follow installation instructions from README
        or: modify official repository (https://github.com/snopt/SNOPT7.jl)
            to have increased workspace size.

    SNOPT license is required.
"""
function include_snopt(path = joinpath(dirname(pwd()),"SNOPT7.jl/src/SNOPT7.jl"))
    include(path)
end


text:
@sumProductRule(:node_type     => Nonlinear{Unscented},
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Nothing, Message{Gaussian}),
                :name          => SPNonlinearUTOutNG)

@sumProductRule(:node_type     => Nonlinear{Unscented},
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Message{Gaussian}, Nothing),
                :name          => SPNonlinearUTIn1GG)

mutable struct SPNonlinearUTOutNGX <: SumProductRule{Nonlinear{Unscented}} end
outboundType(::Type{SPNonlinearUTOutNGX}) = Message{GaussianMeanVariance}
function isApplicable(::Type{SPNonlinearUTOutNGX}, input_types::Vector{<:Type})
    total_inputs = length(input_types)
    (total_inputs > 2) || return false
    (input_types[1] == Nothing) || return false

    for input_type in input_types[2:end]
        matches(input_type, Message{Gaussian}) || return false
    end

    return true
end

mutable struct SPNonlinearUTInGX <: SumProductRule{Nonlinear{Unscented}} end
outboundType(::Type{SPNonlinearUTInGX}) = Message{Gaussian}
function isApplicable(::Type{SPNonlinearUTInGX}, input_types::Vector{<:Type})
    total_inputs = length(input_types)
    (total_inputs > 2) || return false
    (input_types[1] != Nothing) || return false

    nothing_inputs = 0
    gaussian_inputs = 0
    for input_type in input_types
        if input_type == Nothing
            nothing_inputs += 1
        elseif matches(input_type, Message{Gaussian})
            gaussian_inputs += 1
        end
    end

    return (nothing_inputs == 1) && (gaussian_inputs == total_inputs-1)
end

mutable struct MNonlinearUTInGX <: MarginalRule{Nonlinear{Unscented}} end
function isApplicable(::Type{MNonlinearUTInGX}, input_types::Vector{<:Type})
    total_inputs = length(input_types)
    (total_inputs > 2) || return false
    (input_types[1] == Nothing) || return false # Indicates marginalization over outbound variable

    for input_type in input_types[2:end]
        matches(input_type, Message{Gaussian}) || return false
    end

    return true
end


text:
#module Visualize
#using Bridge, PyPlot, StaticVector
#end
using JLD2, FileIO
const d = 2

sim = [:pendulum, :pendulum][simid]
simname = String(sim)
println(ENV["BRIDGE_OUTDIR"])

mkpath(joinpath("output", simname))
try # save cp of this file as documentation
    cp(@__FILE__(), joinpath("output",simname,"$simname.jl"); force=true)
catch
end

iterations = 100000
saveit = div(iterations, 20)
alpha = 0.5 # 1 - exp(-alpha Exp) is AR(1) coefficient of Brownian motion valued random walk
independent = false # independent proposals
adaptit = 1000 # adapt every `it`th step
adaptmax = iterations

if simid == 1
    partial = true
    adaptive = true
    initnu = :tilde
    alpha = 0.5
    burnin = 2000
elseif simid == 2
    partial = true
    adaptive = true
    initnu = :tilde
    alpha = 0.5
    burnin = 2000
else
    error("provide `simid in 1:2`")
end
Σ_ =  [1.0, 0.001][simid]

include("pendulum.jl")

XX = Vector{typeof(X)}(m)
XXmean = Vector{typeof(X)}(m)
XXᵒ = Vector{typeof(X)}(m)
WW = Vector{typeof(W)}(m)
WWᵒ = Vector{typeof(W)}(m)

# Create linear noise approximations

Pt = Vector(m)
Pᵒ = Vector(m)
H♢, v = Bridge.gpupdate(πH*one(SM), zero(SV), L, Σ, V.yy[end])

for i in m:-1:1
    tt_ = range(V.tt[i], stop=V.tt[i+1], length=M+1)
    XX[i] = Bridge.samplepath(tt_, zero(SV))
    WW[i] = Bridge.samplepath(tt_, zero(RV))

    a_ = Bridge.a(XX[i].tt[end], v, P)

     if initnu == :brown
        Pt[i] = Bridge.linearappr(SamplePath(XX[i].tt, zeros(SV, M+1)), Bridge.NoDrift(P))
    elseif initnu == :backward
        Pt[i] = Bridge.LinearNoiseAppr(XX[i].tt, P, v, a_, :backward)
    elseif initnu == :tilde
        Pt[i] = Bridge.linearappr(SamplePath(XX[i].tt, fill(zero(SV), M+1)), Tilde(P))
    end

    Pᵒ[i] = Bridge.GuidedBridge(XX[i].tt, P, Pt[i], v, H♢)
    H♢, v = Bridge.gpupdate(Pᵒ[i], L, Σ, V.yy[i])
end

π0 = Bridge.Gaussian(v, Hermitian(H♢))

# initialize

y = π0.μ
for i in 1:m
    sample!(WW[i], Wiener{RV}())
    y = Bridge.bridge!(XX[i], y, WW[i], Pᵒ[i])
end
XXmean = [zero(XX[i]) for i in 1:m]

function smooth(π0, XX, WW, P, Pᵒ, iterations, alpha; verbose = true,adaptive = true, adaptmax = iterations, adaptit = 5000, saveit = 500, independent = false,
    smoothmean = false, hwindow=20)
    Paths = []

    m = length(XX)

    # create workspace
    XXᵒ = deepcopy(XX)
    WWᵒ = deepcopy(WW)
    Wnr = Wiener{valtype(WW[1])}()

    # initialize
    mcstate = [mcstart(XX[i].yy) for i in 1:m]
    acc = 0
    y0 = π0.μ
    X0 = typeof(y0)[]
    Xt = typeof(y0)[]
    newblock = false
    for it in 1:iterations
        doaccept = false
        if it % saveit == 0
            push!(Paths, collect(Iterators.flatten(XX[i].yy[1:end-1] for i in 1:m)))
        end

        if isdefined(:burnin) && it == burnin
            mcstate = [mcstart(XX[i].yy) for i in 1:m]
        end

        if adaptive && it < adaptmax && it % adaptit == 0 # adaptive smoothing
            H♢, v = Bridge.gpupdate(πH*one(SM), zero(SV), L, Σ, V.yy[end])
            for i in m:-1:1
                xx = mcstate[i][1]
                if isa(Pᵒ[i].Pt, Bridge.LinearAppr)
                    if smoothmean
                        Y = SamplePath(XX[i].tt, [mean(xx[max(1, j-hwindow):min(end, j+hwindow)]) for j in 1:length(xx)])
                    else
                        Y = SamplePath(XX[i].tt, xx)
                    end
                    Bridge.linearappr!(Pᵒ[i].Pt, Y, P)
                else
                    error("Adaption not implemented for proposals Pᵒ.Pt")
                end
                Pᵒ[i] = Bridge.GuidedBridge(XX[i].tt, P, Pᵒ[i].Pt, v, H♢)
                H♢, v = Bridge.gpupdate(Pᵒ[i], L, Σ, V.yy[i])
            end
            π0 = Bridge.Gaussian(v, Hermitian(H♢))
            #y0ᵒ = rand(π0)
            newblock = true
            if it == adaptit # first adaptive proposal needs some help sometimes
                doaccept = true
            end

        end

        push!(X0, y0)
        push!(Xt, XX[div(end,2)].yy[1])
        rho_ = exp(-alpha*randexp())
        if newblock
            #y0ᵒ = rand(π0)
            y0ᵒ = y0
        elseif !independent

            y0ᵒ = π0.μ + sqrt(rho_)*(rand(π0) - π0.μ) + sqrt(1-rho_)*(y0 - π0.μ)
        else
            y0ᵒ = rand(π0)
        end
        y = y0ᵒ
        for i in 1:m
            sample!(WWᵒ[i], Wnr)
            if !independent
                WWᵒ[i].yy[:] = sqrt(rho_)*WWᵒ[i].yy + sqrt(1-rho_)*WW[i].yy
            end
            y = Bridge.bridge!(XXᵒ[i], y, WWᵒ[i], Pᵒ[i])
        end

        ll = 0.0
        for i in 1:m
            ll += llikelihood(LeftRule(), XXᵒ[i],  Pᵒ[i]) - llikelihood(LeftRule(), XX[i],  Pᵒ[i])
        end
        verbose && print(@sprintf("%7i", it), " ll ",  @sprintf("%10.3f", ll))

        #if true
        if doaccept || rand() < exp(ll)
            acc += 1
            verbose && print("\t ✓")
            y0 = y0ᵒ
            for i in 1:m
                XX[i], XXᵒ[i] = XXᵒ[i], XX[i]
                WW[i], WWᵒ[i] = WWᵒ[i], WW[i]
            end
            newblock = false
        else
            verbose && print("\t .")
        end
        verbose && println([" ", "N"][1 + newblock], "\t", round(acc/it, digits=2), "\t", round.(y0 - x0, 2))
        for i in 1:m
            mcstate[i] = Bridge.mcnext!(mcstate[i],XX[i].yy)
        end
    end
    Paths, X0, Xt, mcstate, acc
end

Paths, X0, Xt, mcstates, acc = smooth(π0, XX, WW, P, Pᵒ, iterations, alpha;
adaptive = adaptive,
adaptit = adaptit,
adaptmax = adaptmax,
saveit = saveit,
verbose = true, independent = independent)

open(joinpath("output", simname,"info.txt"), "w") do f
    println(f, "acc")
    println(f, round(acc/iterations, digits=3))
end

writedlm(joinpath("output", simname, "x0n$simid.csv"), [1:iterations Bridge.mat(X0)'])
writedlm(joinpath("output", simname, "xtn$simid.csv"), [1:iterations Bridge.mat(Xt)'])

# Plot result

Xmeanm, Xstdm = Bridge.mcmarginalstats(mcstates)

save(joinpath(ENV["BRIDGE_OUTDIR"],  "$simname$(simid)paths.jld2"), "Path", Paths)
save(joinpath(ENV["BRIDGE_OUTDIR"],  "$simname$(simid)states.jld2"), "mcstates", mcstates)

for i in 1:2
    writedlm(joinpath(ENV["BRIDGE_OUTDIR"],  "$simname$(simid)pathsx$i.csv"),  hcat(map(x->getindex.(x, i), Paths)...))
    writedlm(joinpath(ENV["BRIDGE_OUTDIR"],  "$simname$(simid)meanx$i.csv"),  getindex.(Xmeanm, i))
    writedlm(joinpath(ENV["BRIDGE_OUTDIR"],  "$simname$(simid)stdx$i.csv"),  getindex.(Xstdm, i))
    writedlm(joinpath(ENV["BRIDGE_OUTDIR"],  "$simname$(simid)x0$i.csv"),  getindex.(X0, i))
    writedlm(joinpath(ENV["BRIDGE_OUTDIR"],  "$simname$(simid)xt$i.csv"),  getindex.(Xt, i))
    writedlm(joinpath(ENV["BRIDGE_OUTDIR"],  "$simname$(simid)Xtrue$i.csv"),  [Xtrue.tt getindex.(Xtrue.yy, i)])

end
writedlm(joinpath(ENV["BRIDGE_OUTDIR"],  "$simname$(simid)V1.csv"),  [V.tt V.yy])


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

module Rotations3D

using StaticArrays
using LinearAlgebra: norm, rank, svd, Diagonal, dot
using ACE.SphericalHarmonics: index_y
using Combinatorics: permutations

export ClebschGordan, Rot3DCoeffs, ri_basis, rpi_basis

"""
`ClebschGordan: ` storing precomputed Clebsch-Gordan coefficients; see
`?clebschgordan` for the convention that is use.
"""
struct ClebschGordan{T}
	vals::Dict{Tuple{Int, Int, Int, Int, Int, Int}, T}
end

"""
`ERot3DCoeffs: ` storing recursively precomputed matrix-valued coefficients for a
rotation-equivariant basis. L = 0 equiv invariant, L = 1 equiv equivariant
"""
struct Rot3DCoeffs{T,L}
   vals::Vector{Dict}
   cg::ClebschGordan{T}
end

# -----------------------------------
# iterating over an m collection
# -----------------------------------

_mvec(::CartesianIndex{0}) = SVector(Int(0))

_mvec(mpre::CartesianIndex) = SVector(Tuple(mpre)..., - sum(Tuple(mpre)))

struct MRange{N, T2}
   ll::SVector{N, Int}
   cartrg::T2
   L::Int
end

Base.length(mr::MRange) = sum(_->1, _mrange(mr.ll,mr.L))

_mrange(ll, L) = MRange(ll, Iterators.Stateful(
               filter((x) -> abs(sum(x))<= L, Tuple.(CartesianIndices(ntuple(i -> -ll[i]:ll[i], length(ll)))))
                     ),L)

function Base.iterate(mr::MRange, args...)
	while true
   		if isempty(mr.cartrg)
      		return nothing
   		end
   		mpre = popfirst!(mr.cartrg)
   		return SVector(mpre), nothing
	end
end

# ----------------------------------------------------------------------
#     ClebschGordan code
# ----------------------------------------------------------------------

cg_conditions(j1,m1, j2,m2, J,M) =
	cg_l_condition(j1, j2, J)   &&
	cg_m_condition(m1, m2, M)   &&
	(abs(m1) <= j1) && (abs(m2) <= j2) && (abs(M) <= J)

cg_l_condition(j1, j2, J) = (abs(j1-j2) <= J <= j1 + j2)

cg_m_condition(m1, m2, M) = (M == m1 + m2)

"""
`clebschgordan(j1, m1, j2, m2, J, M, T=Float64)` :

A reference implementation of Clebsch-Gordon coefficients based on

https://hal.inria.fr/hal-01851097/document
Equation (4-6)

This heavily uses BigInt and BigFloat and should therefore not be employed
for performance critical tasks, but only precomputation.

The ordering of parameters corresponds to the following convention:
```
clebschgordan(j1, m1, j2, m2, J, M) = C_{j1m1j2m2}^{JM}
```
where
```
   D_{m1k1}^{l1} D_{m2k2}^{l2}}
	=
	∑_j  C_{l1m1l2m2}^{j(m1+m2)} C_{l1k1l2k2}^{j2(k1+k2)} D_{(m1+m2)(k1+k2)}^{j}
```
"""
function clebschgordan(j1, m1, j2, m2, J, M, T=Float64)
	if !cg_conditions(j1, m1, j2, m2, J, M)
		return zero(T)
	end

   N = (2*J+1) *
       factorial(big(j1+m1)) * factorial(big(j1-m1)) *
       factorial(big(j2+m2)) * factorial(big(j2-m2)) *
       factorial(big(J+M)) * factorial(big(J-M)) /
       factorial(big( j1+j2-J)) /
       factorial(big( j1-j2+J)) /
       factorial(big(-j1+j2+J)) /
       factorial(big(j1+j2+J+1))

   G = big(0)
   # 0 ≦ k ≦ j1+j2-J
   # 0 ≤ j1-m1-k ≤ j1-j2+J   <=>   j2-J-m1 ≤ k ≤ j1-m1
   # 0 ≤ j2+m2-k ≤ -j1+j2+J  <=>   j1-J+m2 ≤ k ≤ j2+m2
   lb = (0, j2-J-m1, j1-J+m2)
   ub = (j1+j2-J, j1-m1, j2+m2)
   for k in maximum(lb):minimum(ub)
      bk = big(k)
      G += (-1)^k *
           binomial(big( j1+j2-J), big(k)) *
           binomial(big( j1-j2+J), big(j1-m1-k)) *
           binomial(big(-j1+j2+J), big(j2+m2-k))
   end

   return T(sqrt(N) * G)
end

ClebschGordan(T=Float64) =
	ClebschGordan{T}(Dict{Tuple{Int,Int,Int,Int,Int,Int}, T}())

_cg_key(j1, m1, j2, m2, J, M) = (j1, m1, j2, m2, J, M)
	# Int.((index_y(j1,m1), index_y(j2,m2), index_y(J,M)))

function (cg::ClebschGordan{T})(j1, m1, j2, m2, J, M) where {T}
	if !cg_conditions(j1,m1, j2,m2, J,M)
		return zero(T)
	end
	key = _cg_key(j1, m1, j2, m2, J, M)
	if haskey(cg.vals, key)
		return cg.vals[key]
	end
	val = clebschgordan(j1, m1, j2, m2, J, M, T)
	cg.vals[key] = val
	return val
end

# ----------------------------------------------------------------------
#     IRot3DCoeffs code: generalized cg coefficients
#
#  Note: in this section kk is a tuple of m-values, it is not
#        related to the k index in the 1-p basis (or radial basis)
# ----------------------------------------------------------------------

function get0(L::Int,T=Float64)
	if L== 0
		return T(0)
	elseif L == 1
		return @SArray zeros(Complex{T},3,3)
		#SMatrix{3, 3, Complex{T}, 9}(0, 0, 0, 0, 0, 0, 0, 0, 0)
	else
		#For L== 2 @SArray zeros(Complex{T},3,3,3)
		ErrorException("Only types for L in {0,1} implemented")
	end
end
dicttype(A::Rot3DCoeffs,N::Integer) = dicttype(A::Rot3DCoeffs,Val(N))

dicttype(A::Rot3DCoeffs{T,L},::Val{N}) where {T,L,N} =
   Dict{Tuple{SVector{N,Int}, SVector{N,Int}, SVector{N,Int}}, typeof(get0(L,T))}

Rot3DCoeffs(L,T=Float64) = Rot3DCoeffs{T,L}(Dict[], ClebschGordan(T))

function get_vals(A::Rot3DCoeffs, valN::Val{N}) where {N}
	if length(A.vals) < N
		for n = length(A.vals)+1:N
			push!(A.vals, dicttype(A,n)())
		end
	end
   	return A.vals[N]::dicttype(A,valN)
end

_key(ll::StaticVector{N}, mm::StaticVector{N}, kk::StaticVector{N}) where {N} =
      (SVector{N, Int}(ll), SVector{N, Int}(mm), SVector{N, Int}(kk))

  function (A::Rot3DCoeffs{T,L})(ll::StaticVector{N},
                              mm::StaticVector{N},
                              kk::StaticVector{N}) where {T,L, N}
     if       abs(sum(mm)) > L ||
              abs(sum(kk)) > L ||
              !all(abs.(mm) .<= ll) ||
              !all(abs.(kk) .<= ll)
        return get0(L,T)
     end
     vals = get_vals(A, Val(N))  # this should infer the type!
     key = _key(ll, mm, kk)
     if haskey(vals, key)
        val  = vals[key]
     else
        val = _compute_val(A, key...)
        vals[key] = val
     end
     return val
  end

  function (A::Rot3DCoeffs{T,0})(ll::StaticVector{1},
                              mm::StaticVector{1},
                              kk::StaticVector{1}) where {T}
     if ll[1] == mm[1] == kk[1] == 0
        return T(8 * pi^2)
     else
        return get0(0,T)
     end
  end

  function (A::Rot3DCoeffs{T,1})(ll::StaticVector{1},
                              mm::StaticVector{1},
                              kk::StaticVector{1}) where {T}
     if ll[1] == 1 && abs(mm[1]) <= 1 && abs(kk[1]) <= 1
        return  rmatrices[(mm[1],kk[1])]
     else
        return get0(1,T)
     end
  end

#function get0val(m,k,L)

#end
rmatrices = Dict(
    (-1,-1) => SMatrix{3, 3, ComplexF64, 9}(1/6, 1im/6, 0, -1im/6, 1/6, 0, 0, 0, 0),
    (-1,0) => SMatrix{3, 3, ComplexF64, 9}(0, 0, 0, 0, 0, 0, 1/(3*sqrt(2)), 1im/(3*sqrt(2)), 0),
    (-1,1) => SMatrix{3, 3, ComplexF64, 9}(-1/6, -1im/6, 0, -1im/6, 1/6, 0, 0, 0, 0),
    (0,-1) => SMatrix{3, 3, ComplexF64, 9}(0, 0, 1/(3*sqrt(2)), 0, 0, -1im/(3*sqrt(2)), 0, 0, 0),
    (0,0) => SMatrix{3, 3, ComplexF64, 9}(0, 0, 0, 0, 0, 0, 0, 0, 1/3),
    (0,1) => SMatrix{3, 3, ComplexF64, 9}(0, 0, -1/(3*sqrt(2)), 0, 0, -1im/(3*sqrt(2)), 0, 0, 0),
    (1,-1) => SMatrix{3, 3, ComplexF64, 9}(-1/6, 1im/6, 0, 1im/6, 1/6, 0, 0, 0, 0),
    (1,0) => SMatrix{3, 3, ComplexF64, 9}(0, 0, 0, 0, 0, 0, -1/(3*sqrt(2)), 1im/(3*sqrt(2)), 0),
    (1,1) => SMatrix{3, 3, ComplexF64, 9}(1/6, -1im/6, 0, 1im/6, 1/6, 0, 0, 0, 0)
    )

function (A::Rot3DCoeffs{T,L})(ll::StaticVector{1},
							 mm::StaticVector{1},
							 kk::StaticVector{1}) where {T,L}
	ErrorException("Not implemented for L = " + string(L))
end

function _compute_val(A::Rot3DCoeffs{T,L}, ll::StaticVector{N},
                                        mm::StaticVector{N},
                                        kk::StaticVector{N}) where {T,L,N}

   val = get0(L,T)
   llp = ll[1:N-2]
   mmp = mm[1:N-2]
   kkp = kk[1:N-2]
   for j = abs(ll[N-1]-ll[N]):(ll[N-1]+ll[N])
      if abs(kk[N-1]+kk[N]) > j || abs(mm[N-1]+mm[N]) > j
         continue
      end
		cgk = try
			A.cg(ll[N-1], kk[N-1], ll[N], kk[N], j, kk[N-1]+kk[N])
		catch
			@show (ll[N-1], kk[N-1], ll[N], kk[N], j, kk[N-1]+kk[N])
			get0(L,T)
		end
		cgm = A.cg(ll[N-1], mm[N-1], ll[N], mm[N], j, mm[N-1]+mm[N])
		if cgk * cgm  != 0
			val += cgk * cgm * A( SVector(llp..., j),
								       SVector(mmp..., mm[N-1]+mm[N]),
								       SVector(kkp..., kk[N-1]+kk[N]) )
		end
   end
   return val
end

function re_basis(A::Rot3DCoeffs{T,0}, ll::SVector{N, Int}; ordered=false) where {T,N}
	"""computes coefficients for a rotational 0-equinvariant (invariant) tensor basis"""
	CC = compute_Al(A, ll, Val(ordered))
	svdC = svd(CC)
	rk = rank(Diagonal(svdC.S))
	return svdC.U[:, 1:rk]'
end

function compute_Al(A::Rot3DCoeffs{T,0}, ll::SVector, ::Val{false}) where {T}
	len = length(_mrange(ll,0))
   CC = zeros(T, len, len)
   for (im, mm) in enumerate(_mrange(ll,0)), (ik, kk) in enumerate(_mrange(ll,0))
      CC[ik, im] = A(ll, mm, kk)
   end
   return CC
end

function re_basis(A::Rot3DCoeffs{T,L}, ll::SVector{N, Int}) where {T,L,N}
	"""return coefficients for a rotational 1-equinvariant (covariant) tensor basis"""
	#Mre = collect(Iterators.product(1:3, _mrange(ll,1)))
	GG = _gramian(A, ll)
   	S = svd(GG)
   	rk = rank(GG; rtol =  1e-7)
	Ure = S.U[:, 1:rk]'

	return Diagonal(1. ./ sqrt.(S.S[1:rk])) * Ure
end

imag_tol = 1e-15

function _gramian(A::Rot3DCoeffs{T,1}, ll::SVector) where {T}
	len = length(_mrange(ll,1))*3
	GG = zeros(ComplexF64, len, len) # = Gramian of rotational eqiuvariant span
	# kk = mprime
    for (i, (k,mm)) in enumerate(Iterators.product(1:3, _mrange(ll,1))), (ip, (kp,mmp)) in enumerate(Iterators.product(1:3, _mrange(ll,1)))
		for mu in _mrange(ll,1)
		  	GG[i,ip] += dot(A(ll, mm, mu)[:,k],A(ll, mmp, mu)[:,kp])
		end
	end
	@assert all(abs.(imag(GG)) .<= imag_tol)
    return real.(GG)
end

rpi_basis(A::Rot3DCoeffs, zz, nn, ll) =
			rpi_basis(A, SVector(zz...), SVector(nn...), SVector(ll...))

function rpi_basis(A::Rot3DCoeffs{T,L},
						 zz::SVector{N},
						 nn::SVector{N, Int},
						 ll::SVector{N, Int}) where {T,L,N}
	Ure = re_basis(A, ll)
	Mre = collect( _mrange(ll, L))   # rows...
	G = _gramian(zz, nn, ll, Ure, Mre)
	#print("G type :", typeof(G),"\n")
   	S = svd(G)
   	rk = rank(G; rtol =  1e-7)
	Urpi = S.U[:, 1:rk]'
	#print("Urpi type :", typeof(Urpi),"\n")
	return Diagonal(sqrt.(S.S[1:rk])) * Urpi * Ure, Mre
end

function rpi_basis(A::Rot3DCoeffs{T,1},
						 zz::SVector{N},
						 nn::SVector{N, Int},
						 ll::SVector{N, Int}) where {T,N}
	Mre = collect( _mrange(ll, 1))   # rows...
	G = _gramian(A, zz, nn, ll, Mre)
	#print("G type :", typeof(G),"\n")
   	S = svd(G)
   	rk = rank(G; rtol =  1e-7)
	Urepi = Diagonal(1/sqrt.(S.S[1:rk])) * S.U[:, 1:rk]'
	Utilde = zeros(SArray{Tuple{3},Complex{T},1,3}, rk, length(Mre))
	for α in 1:rk
		for (imu, mu) in enumerate(Mre)
			for (i, (k,mm)) in enumerate(Iterators.product(1:3, Mre))
				Utilde[α,imu] += Urepi[α,i] * A(ll,mu,mm)[:,k]
			end
		end
	end
	#print("Urpi type :", typeof(Urpi),"\n")
	return Utilde, Mre
end

function _gramian(A::Rot3DCoeffs{T,1}, zz, nn, ll::SVector{N}, Mre) where{T,N}
  	len = length(Mre)*3
	GG = zeros(Complex{T}, len, len)

	for (i1, (k1,mm1)) in enumerate(Iterators.product(1:3, _mrange(ll, 1))), (i2, (k2,mm2)) in enumerate(Iterators.product(1:3, _mrange(ll, 1)))
		for σ_ in permutations(1:N)
			σ = SVector{N,Int}(σ_)
	      	if ( zz[σ] != zz) || (nn[σ] != nn) || (ll[σ] != ll); continue; end
			for mu1 in  _mrange(ll, 1), mu2 in _mrange(ll, 1)
	        	if mu1[σ] == mu2
					GG[i1, i2] += dot(A(ll, mu2, mm2)[:,k2],A(ll[σ], mu1[σ], mm1[σ])[:,k1])
	         	end
	      	end
	   	end
	end
   	return GG
	#@assert all(abs.(imag(GG)) .<= imag_tol)
	#return real.(GG)
end

function _gramian(zz, nn, ll, Uri, Mri)
   N = length(nn)
   nri = size(Uri, 1)
   @assert size(Uri, 1) == nri
   G = zeros(Float64,nri, nri)

   for σ in permutations(1:N)
      if (zz[σ] != zz) || (nn[σ] != nn) || (ll[σ] != ll); continue; end
      for (iU1, mm1) in enumerate(Mri), (iU2, mm2) in enumerate(Mri)
         if mm1[σ] == mm2
            for i1 = 1:nri, i2 = 1:nri
               G[i1, i2] += conj(Uri[i1, iU1]) * Uri[i2, iU2]
            end
         end
      end
   end
   return G
   #@assert all(abs.(imag(G)) .<= imag_tol)
   #return real.(G)
end

end


text:
cases = ["B-splines","Cheb","Lin"]

# construct basis
holder = (
    Basis(SplineParams(25, -1, 1, 1),SplineParams(20, -1, 2, 3)),
    Basis(ChebParams(25, -1, 1),ChebParams(20, -1, 2)),
    Basis(LinParams(25, -1, 1),LinParams(20, -1, 2))
)

@testset for (i, case) in enumerate(cases)
    basis = holder[i]

    # get nodes
    X, x12 = nodes(basis)

    # function to interpolate
    f(x1, x2) = cos.(x2) ./ exp.(x1)
    f(X::Matrix) = f.(X[:, 1], X[:, 2])

    # function at nodes
    y = f(X)
    yprime2 = - sin.(X[:, 2]) ./ exp.(X[:, 1])

    # benchmark coefficients
    c_direct, bs_direct = funfitxy(basis, X, y)

    # this order isn't included bs_direct
    @test_throws ErrorException funeval(c_direct, convert(Expanded, bs_direct), [0 1])

    @testset "test funfitxy for tensor and direct agree on coefs" begin
        c_tensor, bs_tensor = funfitxy(basis, x12, y)
        @test maximum(abs, c_tensor -  c_direct) <=  1e-12
    end

    @testset "test funfitf" begin
        c = funfitf(basis, f)
        @test maximum(abs, c -  c_direct) <=  1e-12
    end

    @testset "test funeval methods" begin
        # single point
        sp = @inferred funeval(c_direct, basis, X[5:5, :])[1]
        @test maximum(abs, sp -  y[5]) <= 1e-12

        # multiple points using tensor directly
        mp = @inferred funeval(c_direct, basis, x12)
        @test maximum(abs, mp -  y) <=  1e-12

        # multiple points using direct
        mp = @inferred funeval(c_direct, basis, X)
        @test maximum(abs, mp -  y) <=  1e-12

        # multiple points giving basis in direct form
        mpd = @inferred funeval(c_direct, bs_direct)
        @test maximum(abs, mpd -  y) <=  1e-12

        # multiple points giving basis in expanded form
        Phiexp = Base.convert(Expanded, bs_direct)
        mpe = @inferred funeval(c_direct, Phiexp)
        @test maximum(abs, mpe -  y) <=  1e-12

        # order != 0. Note for Spline err is 7e-5. For Cheb 9e-14 and Lin 2e-1
        if case != "Lin"
            der1 = @inferred funeval(c_direct, basis, X, [0 1])
            @test maximum(abs, der1 - yprime2) <= 2e-1
            der2 = @inferred funeval(c_direct, basis, x12, [0 1])
            @test maximum(abs, der1 - yprime2) <= 2e-1
            @test maximum(abs, der1 - der2) <= 1e-12
        end

    end

    @testset "other funeval methods" begin
        b1 = basis[1]
        y1 = sin.(x12[1])
        c = funfitxy(b1, x12[1], y1)[1]
        for (i, x) in enumerate(x12[1])
            @test ≈(y1[i], @inferred(funeval(c, b1, x)), atol=1e-14)
        end

        for (i, x) in enumerate(x12[1])
            @test ≈([y1[i], y1[i]], @inferred(funeval([c c], b1, x)), atol=1e-14)
        end

        @test ≈(y1, @inferred(funeval(c, b1, x12[1])), atol=1e-14)
        @test ≈([y1 y1], @inferred(funeval([c c], b1, x12[1])), atol=1e-14)

    end

    @testset "test interpoland methods" begin
        # (Basis, BasisMatrix, Array)
        intp1 = Interpoland(basis, BasisMatrix(basis, Tensor()), y)
        @test maximum(abs, intp1(X) - y) <= 1e-12

        # (Basis, Array,..)
        intp2 = Interpoland(basis, y)
        @test maximum(abs, intp2(X) - y) <= 1e-12

        # (BasisParams, Function)
        intp3 = Interpoland(basis, f)
        @test maximum(abs, intp3(X) - y) <= 1e-12
    end

    @testset "Printing" begin
        iob = IOBuffer()
        show(iob, Interpoland(basis, y))
    end

    # errors
    @test_throws ErrorException funeval(c_direct, basis, x12, 1)
    @test_throws ErrorException funeval(c_direct, basis, X, 1)

end # testset


text:
# -- PRIVATE METHODS HERE ------------------------------------------------------- #
function check_if_user_dictionary_contains_correct_keys(user_data_dictionary::Dict{String,Any})::PSResult{Bool}

    # the user dictionary should contain a user_data root, and a alpha_vantage_api_key child -
    # TODO: fill me in ...
    return PSResult{Bool}(true)
end

# -- PUBLIC METHODS HERE ------------------------------------------------------- #
"""
    build_api_user_model(path)

Returns either a PSResult{PSError} if something went wrong, or a PSResult{PSUserModel} object holding the user email and AlphaVantage API key.
The PSError and PSUserModel can be accessed using the `value` field on the Result return wrapper
"""
function build_api_user_model(path_to_configuration_file::String)::(Union{PSResult{T}, Nothing} where T<:Any)

    # some users checks -
    # did the user pass in a legit path?
    check_result = is_path_valid(path_to_configuration_file)
    if (typeof(check_result.value) == Bool && check_result.value == false)
        error_message = "error: $(path_to_configuration_file) in not a valid path"
        return PSResult{PSError}(PSError(error_message))
    end

    # ok, path seems legit - load the default user information from the config.json file -
    user_json_dictionary = JSON.parsefile(path_to_configuration_file)

    # does the user dictionary contain the correct keys?
    if (check_if_user_dictionary_contains_correct_keys(user_json_dictionary) == false)
        error_message = "error: missing keys in user configuration dictionary"
        return PSResult{PSError}(PSError(error_message))
    end

    # -- DO NOT EDIT BELOW THIS LINE ------------------------------------------#

    # grab the user data -
    alpha_vantage_api_key = user_json_dictionary["user_data"]["alpha_vantage_api_key"]
    alpha_vantage_api_email = user_json_dictionary["user_data"]["alpha_vantage_api_email"]

    # build APIUserModel -
    api_user_model = PSUserModel(alpha_vantage_api_email, alpha_vantage_api_key)

    # return the user_data_dictionary -
    return PSResult{PSUserModel}(api_user_model)
    # -------------------------------------------------------------------------#
end


text:
module Example
    function entryPoint()
        return "Hello"
    end
end

max_stars_repo_name: ELIFE-ASU/dglmoore

max_issues_repo_name: ELIFE-ASU/dglmoore

max_forks_repo_name: ELIFE-ASU/dglmoore


text:
#!/usr/bin/env julia

using RandomTree

main()


text:
using CVGSurfaceFile
using Test

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

max_stars_repo_name: powerzf/CVGSurfaceFile

max_issues_repo_name: powerzf/CVGSurfaceFile

max_forks_repo_name: powerzf/CVGSurfaceFile


text:
using Documenter, FITSIO

makedocs(
    modules = [FITSIO],
    sitename = "FITSIO.jl",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages    = [
        "Introduction" => "index.md",
        "API Reference" => "api.md"
    ]
)

deploydocs(repo = "github.com/JuliaAstro/FITSIO.jl.git", push_preview=true)


text:
using DataFrames
using VectorAutoregression
using Plots
pyplot()
default(show=true, reuse=false)

qdata = readtable("importable_data/qdata_processed.csv")
mdata = readtable("importable_data/mdata_processed.csv")

# # Aggregate monthly data to quarterly by taking final month values
# finaldata = mdata[collect(range(3,3,length(GDP))),:]
# finaldata[:GDP] = GDP

mfvar_out = VectorAutoregression.var_mf_ols(qdata, mdata, 3, 1)
irfs = VectorAutoregression.oirf(mfvar_out, 20)
varnames = VectorAutoregression.varnames(mfvar_out)

ipshock1 = irfs[1][10, :, 1]
ipshock2 = irfs[1][10, :, 4]
ipshock3 = irfs[1][10, :, 7]
plot(ipshock1, legend=false, linewidth=2)
plot!(ipshock2, legend=false, linewidth=2)
plot!(ipshock3, legend=false, linewidth = 2)
savefig("ipshock.pdf")

inflshock1 = irfs[1][10, :, 2]
inflshock2 = irfs[1][10, :, 5]
inflshock3 = irfs[1][10, :, 8]
plot(inflshock1, legend=false, linewidth=2)
plot!(inflshock2, legend=false, linewidth=2)
plot!(inflshock3, legend=false, linewidth = 2)

unempshock1 = irfs[1][10, :, 3]
unempshock2 = irfs[1][10, :, 6]
unempshock3 = irfs[1][10, :, 9]
plot(unempshock1, legend=false, linewidth=2)
plot!(unempshock2, legend=false, linewidth=2)
plot!(unempshock3, legend=false, linewidth = 2)

max_stars_repo_name: lnsongxf/MF-VAR

max_issues_repo_name: lnsongxf/MF-VAR

max_forks_repo_name: lnsongxf/MF-VAR


text:
module SpinSymmetry

using SparseArrays: sparse

export Flip, Shift, Swap, SpatialReflection, GenericSymmetry
export zbasis, FullZBasis, ZBlockBasis, basissize
export SymmetrizedBasis, symmetrized_basis, symmetrize_state, symmetrize_operator, transformationmatrix

include("abstract.jl")
include("basis.jl")
include("symmetries.jl")

end


text:
using BenchmarkTools
import ParticleSwarmOptimizer
pso = ParticleSwarmOptimizer

pso.get_localbest([0.5, 0.5, 0.5, 0.4, 0.5], [2, 4, 5], <)

@benchmark pso.get_localbest([0.5, 0.5, 0.5, 0.4, 0.5], [2, 4, 5], <)
@benchmark pso.get_localbest([0.5, 0.5, 0.5, 0.4, 0.5], [2], <)


text:
# WARNING: This file was generated from fastq-parser.rl using ragel. Do not edit!
const fastqparser_start  = 18
const fastqparser_first_final  = 18
const fastqparser_error  = 0
const fastqparser_en_main  = 18
const _fastqparser_nfa_targs = Int8[ 0, 0 ,  ]
const _fastqparser_nfa_offsets = Int8[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ,  ]
const _fastqparser_nfa_push_actions = Int8[ 0, 0 ,  ]
const _fastqparser_nfa_pop_trans = Int8[ 0, 0 ,  ]
Ragel.@generate_read!_function(
"fastqparser",
FASTQParser,
SeqRecord,
begin
begin
if ( p == pe  )
	@goto _test_eof

end
if ( cs  == 18 )
	@goto st_case_18
elseif ( cs  == 0 )
	@goto st_case_0
elseif ( cs  == 1 )
	@goto st_case_1
elseif ( cs  == 2 )
	@goto st_case_2
elseif ( cs  == 3 )
	@goto st_case_3
elseif ( cs  == 4 )
	@goto st_case_4
elseif ( cs  == 5 )
	@goto st_case_5
elseif ( cs  == 6 )
	@goto st_case_6
elseif ( cs  == 7 )
	@goto st_case_7
elseif ( cs  == 8 )
	@goto st_case_8
elseif ( cs  == 9 )
	@goto st_case_9
elseif ( cs  == 10 )
	@goto st_case_10
elseif ( cs  == 11 )
	@goto st_case_11
elseif ( cs  == 12 )
	@goto st_case_12
elseif ( cs  == 19 )
	@goto st_case_19
elseif ( cs  == 13 )
	@goto st_case_13
elseif ( cs  == 14 )
	@goto st_case_14
elseif ( cs  == 15 )
	@goto st_case_15
elseif ( cs  == 16 )
	@goto st_case_16
elseif ( cs  == 17 )
	@goto st_case_17
end
@goto st_out
@label ctr39
begin
state.linenum += 1
end
@goto st18
@label st18
p+= 1;
if 	( p == pe  )
	@goto _test_eof18

end
@label st_case_18
if ( (data[1+(p )]) == 10 )
	begin
	@goto ctr39

end
elseif 	( (data[1+(p )]) == 32 )
	begin
	@goto st18

end
elseif 	( (data[1+(p )]) == 64 )
	begin
	ck = 0;
	if ( (length(input.qualbuf) == length(input.seqbuf)) )
		ck += 1;

	end
	if ( 1 <= ck  )
		begin
		@goto st1

	end

end
@goto st0

end
end
if ( 9 <= (data[1+(p )])&& (data[1+(p )])<= 13  )
	begin
	@goto st18

end

end
begin
@goto st0

end
@label st_case_0
@label st0
cs = 0;
@goto _out
@label ctr41
begin
if length(input.seqbuf) != length(input.qualbuf)
error("Error parsing FASTQ: sequence and quality scores must be of equal length")
end

# if a name and s description are repeated (i.e. old-fashioned fastq)
# make sure they match.
if (!isempty(input.name2buf) && input.name2buf != output.name) ||
(!isempty(input.desc2buf) && input.desc2buf != output.metadata.description)
error("Error parsing FASTQ: sequence and quality scores have non-matching identifiers")
end

# sequence
resize!(output.seq, input.seqbuf.position - 1)
encode_copy!(output.seq, 1, input.seqbuf.buffer, 1, input.seqbuf.position - 1)

# quality
encoding, input.quality_encodings =
infer_quality_encoding(input.qualbuf.buffer, 1,
input.qualbuf.position - 1,
input.quality_encodings)

decode_quality_string!(encoding, input.qualbuf.buffer,
output.metadata.quality, 1,
input.qualbuf.position - 1)

# reset temporaries for the next run
empty!(input.qualbuf)
empty!(input.seqbuf)
empty!(input.name2buf)
empty!(input.desc2buf)

Ragel.@yield 1
end
@goto st1
@label st1
p+= 1;
if ( p == pe  )
@goto _test_eof1

end
@label st_case_1
if ( (data[1+(p )])== 32  )
begin
@goto st0

end

end
if ( 9 <= (data[1+(p )])&& (data[1+(p )])<= 13  )
begin
@goto st0

end

end
begin
@goto ctr0

end
@label ctr0
begin
Ragel.@anchor!
end
@goto st2
@label st2
p+= 1;
if ( p == pe  )
@goto _test_eof2

end
@label st_case_2
if ( (data[1+(p )]) == 10 )
begin
@goto ctr4

end
elseif ( (data[1+(p )]) == 12 )
begin
@goto st0

end
elseif ( (data[1+(p )]) == 13 )
begin
@goto ctr5

end
elseif ( (data[1+(p )]) == 32 )
begin
@goto ctr3

end
end
if ( 9 <= (data[1+(p )])&& (data[1+(p )])<= 11  )
begin
@goto ctr3

end

end
begin
@goto st2

end
@label ctr3
begin
Ragel.@copy_from_anchor!(output.name)
end
@goto st3
@label st3
p+= 1;
if ( p == pe  )
@goto _test_eof3

end
@label st_case_3
if ( (data[1+(p )]) == 10 )
begin
@goto st0

end
elseif ( (data[1+(p )]) == 32 )
begin
@goto st3

end
end
if ( (data[1+(p )])> 11  )
begin
if ( (data[1+(p )])<= 13  )
begin
@goto st0

end

end

end

elseif ( (data[1+(p )])>= 9  )
begin
@goto st3

end

end
begin
@goto ctr6

end
@label ctr6
begin
Ragel.@anchor!
end
@goto st4
@label st4
p+= 1;
if ( p == pe  )
@goto _test_eof4

end
@label st_case_4
if ( (data[1+(p )]) == 10 )
begin
@goto ctr9

end
elseif ( (data[1+(p )]) == 13 )
begin
@goto ctr10

end
end
begin
@goto st4

end
@label ctr4
begin
Ragel.@copy_from_anchor!(output.name)
end
begin
state.linenum += 1
end
@goto st5
@label ctr9
begin
Ragel.@copy_from_anchor!(output.metadata.description)
end
begin
state.linenum += 1
end
@goto st5
@label ctr37
begin
state.linenum += 1
end
@goto st5
@label st5
p+= 1;
if ( p == pe  )
@goto _test_eof5

end
@label st_case_5
if ( (data[1+(p )]) == 10 )
begin
@goto ctr11

end
elseif ( (data[1+(p )]) == 13 )
begin
@goto st7

end
end
if ( 65 <= (data[1+(p )])&& (data[1+(p )])<= 122  )
begin
@goto ctr13

end

end
begin
@goto st0

end
@label ctr11
begin
state.linenum += 1
end
@goto st6
@label ctr34
begin
Ragel.@append_from_anchor!(input.seqbuf)
end
begin
state.linenum += 1
end
@goto st6
@label st6
p+= 1;
if ( p == pe  )
@goto _test_eof6

end
@label st_case_6
if ( (data[1+(p )]) == 10 )
begin
@goto ctr11

end
elseif ( (data[1+(p )]) == 13 )
begin
@goto st7

end
elseif ( (data[1+(p )]) == 43 )
begin
@goto st8

end
end
if ( 65 <= (data[1+(p )])&& (data[1+(p )])<= 122  )
begin
@goto ctr13

end

end
begin
@goto st0

end
@label ctr35
begin
Ragel.@append_from_anchor!(input.seqbuf)
end
@goto st7
@label st7
p+= 1;
if ( p == pe  )
@goto _test_eof7

end
@label st_case_7
if ( (data[1+(p )])== 10  )
begin
@goto ctr11

end

end
begin
@goto st0

end
@label st8
p+= 1;
if ( p == pe  )
@goto _test_eof8

end
@label st_case_8
if ( (data[1+(p )]) == 10 )
begin
@goto ctr16

end
elseif ( (data[1+(p )]) == 13 )
begin
@goto st15

end
elseif ( (data[1+(p )]) == 32 )
begin
@goto st0

end
end
if ( 9 <= (data[1+(p )])&& (data[1+(p )])<= 12  )
begin
@goto st0

end

end
begin
@goto ctr15

end
@label ctr15
begin
Ragel.@anchor!
end
@goto st9
@label st9
p+= 1;
if ( p == pe  )
@goto _test_eof9

end
@label st_case_9
if ( (data[1+(p )]) == 10 )
begin
@goto ctr20

end
elseif ( (data[1+(p )]) == 12 )
begin
@goto st0

end
elseif ( (data[1+(p )]) == 13 )
begin
@goto ctr21

end
elseif ( (data[1+(p )]) == 32 )
begin
@goto ctr19

end
end
if ( 9 <= (data[1+(p )])&& (data[1+(p )])<= 11  )
begin
@goto ctr19

end

end
begin
@goto st9

end
@label ctr19
begin
Ragel.@copy_from_anchor!(input.name2buf)
end
@goto st10
@label st10
p+= 1;
if ( p == pe  )
@goto _test_eof10

end
@label st_case_10
if ( (data[1+(p )]) == 10 )
begin
@goto st0

end
elseif ( (data[1+(p )]) == 32 )
begin
@goto st10

end
end
if ( (data[1+(p )])> 11  )
begin
if ( (data[1+(p )])<= 13  )
begin
@goto st0

end

end

end

elseif ( (data[1+(p )])>= 9  )
begin
@goto st10

end

end
begin
@goto ctr22

end
@label ctr22
begin
Ragel.@anchor!
end
@goto st11
@label st11
p+= 1;
if ( p == pe  )
@goto _test_eof11

end
@label st_case_11
if ( (data[1+(p )]) == 10 )
begin
@goto ctr25

end
elseif ( (data[1+(p )]) == 13 )
begin
@goto ctr26

end
end
begin
@goto st11

end
@label ctr16
begin
state.linenum += 1
end
@goto st12
@label ctr20
begin
Ragel.@copy_from_anchor!(input.name2buf)
end
begin
state.linenum += 1
end
@goto st12
@label ctr25
begin
Ragel.@copy_from_anchor!(input.desc2buf)
end
begin
state.linenum += 1
end
@goto st12
@label st12
p+= 1;
if ( p == pe  )
@goto _test_eof12

end
@label st_case_12
if ( (data[1+(p )]) == 10 )
begin
@goto ctr27

end
elseif ( (data[1+(p )]) == 13 )
begin
@goto st13

end
end
if ( 33 <= (data[1+(p )])&& (data[1+(p )])<= 126  )
begin
ck = 0;
if ( (length(input.qualbuf) + input.qualcount < length(input.seqbuf)) )
ck += 1;

end
if ( 1 <= ck  )
begin
@goto ctr29

end

end
@goto st0

end

end
begin
@goto st0

end
@label ctr27
begin
state.linenum += 1
end
@goto st19
@label ctr31
begin
Ragel.@append_from_anchor!(input.qualbuf)
input.qualcount = 0
end
begin
state.linenum += 1
end
@goto st19
@label st19
p+= 1;
if ( p == pe  )
@goto _test_eof19

end
@label st_case_19
if ( (data[1+(p )]) == 10 )
begin
@goto ctr27

end
elseif ( (data[1+(p )]) == 13 )
begin
@goto st13

end
elseif ( (data[1+(p )]) == 64 )
begin
ck = 0;
if ( (length(input.qualbuf) + input.qualcount < length(input.seqbuf)) )
ck += 1;

end
if ( (length(input.qualbuf) == length(input.seqbuf)) )
ck += 2;

end
if ( ck < 2  )
begin
if ( 1 <= ck  )
begin
@goto ctr29

end

end

end

elseif ( ck > 2  )
begin
@goto ctr29

end

else
begin
@goto ctr41

end

end
@goto st0

end
end
if ( 33 <= (data[1+(p )])&& (data[1+(p )])<= 126  )
begin
ck = 0;
if ( (length(input.qualbuf) + input.qualcount < length(input.seqbuf)) )
ck += 1;

end
if ( 1 <= ck  )
begin
@goto ctr29

end

end
@goto st0

end

end
begin
@goto st0

end
@label ctr32
begin
Ragel.@append_from_anchor!(input.qualbuf)
input.qualcount = 0
end
@goto st13
@label st13
p+= 1;
if ( p == pe  )
@goto _test_eof13

end
@label st_case_13
if ( (data[1+(p )])== 10  )
begin
@goto ctr27

end

end
begin
@goto st0

end
@label ctr29
begin
Ragel.@anchor!
end
begin
input.qualcount += 1
end
@goto st14
@label ctr33
begin
input.qualcount += 1
end
@goto st14
@label st14
p+= 1;
if ( p == pe  )
@goto _test_eof14

end
@label st_case_14
if ( (data[1+(p )]) == 10 )
begin
@goto ctr31

end
elseif ( (data[1+(p )]) == 13 )
begin
@goto ctr32

end
end
if ( 33 <= (data[1+(p )])&& (data[1+(p )])<= 126  )
begin
ck = 0;
if ( (length(input.qualbuf) + input.qualcount < length(input.seqbuf)) )
ck += 1;

end
if ( 1 <= ck  )
begin
@goto ctr33

end

end
@goto st0

end

end
begin
@goto st0

end
@label ctr21
begin
Ragel.@copy_from_anchor!(input.name2buf)
end
@goto st15
@label ctr26
begin
Ragel.@copy_from_anchor!(input.desc2buf)
end
@goto st15
@label st15
p+= 1;
if ( p == pe  )
@goto _test_eof15

end
@label st_case_15
if ( (data[1+(p )])== 10  )
begin
@goto ctr16

end

end
begin
@goto st0

end
@label ctr13
begin
Ragel.@anchor!
end
@goto st16
@label st16
p+= 1;
if ( p == pe  )
@goto _test_eof16

end
@label st_case_16
if ( (data[1+(p )]) == 10 )
begin
@goto ctr34

end
elseif ( (data[1+(p )]) == 13 )
begin
@goto ctr35

end
end
if ( 65 <= (data[1+(p )])&& (data[1+(p )])<= 122  )
begin
@goto st16

end

end
begin
@goto st0

end
@label ctr5
begin
Ragel.@copy_from_anchor!(output.name)
end
@goto st17
@label ctr10
begin
Ragel.@copy_from_anchor!(output.metadata.description)
end
@goto st17
@label st17
p+= 1;
if ( p == pe  )
@goto _test_eof17

end
@label st_case_17
if ( (data[1+(p )])== 10  )
begin
@goto ctr37

end

end
begin
@goto st0

end
@label st_out
@label _test_eof18
cs = 18;
@goto _test_eof
@label _test_eof1
cs = 1;
@goto _test_eof
@label _test_eof2
cs = 2;
@goto _test_eof
@label _test_eof3
cs = 3;
@goto _test_eof
@label _test_eof4
cs = 4;
@goto _test_eof
@label _test_eof5
cs = 5;
@goto _test_eof
@label _test_eof6
cs = 6;
@goto _test_eof
@label _test_eof7
cs = 7;
@goto _test_eof
@label _test_eof8
cs = 8;
@goto _test_eof
@label _test_eof9
cs = 9;
@goto _test_eof
@label _test_eof10
cs = 10;
@goto _test_eof
@label _test_eof11
cs = 11;
@goto _test_eof
@label _test_eof12
cs = 12;
@goto _test_eof
@label _test_eof19
cs = 19;
@goto _test_eof
@label _test_eof13
cs = 13;
@goto _test_eof
@label _test_eof14
cs = 14;
@goto _test_eof
@label _test_eof15
cs = 15;
@goto _test_eof
@label _test_eof16
cs = 16;
@goto _test_eof
@label _test_eof17
cs = 17;
@goto _test_eof
@label _test_eof
begin

end
if ( p == eof  )
begin
if ( cs  == 19 )
begin
if length(input.seqbuf) != length(input.qualbuf)
error("Error parsing FASTQ: sequence and quality scores must be of equal length")
end

# if a name and s description are repeated (i.e. old-fashioned fastq)
# make sure they match.
if (!isempty(input.name2buf) && input.name2buf != output.name) ||
(!isempty(input.desc2buf) && input.desc2buf != output.metadata.description)
error("Error parsing FASTQ: sequence and quality scores have non-matching identifiers")
end

# sequence
resize!(output.seq, input.seqbuf.position - 1)
encode_copy!(output.seq, 1, input.seqbuf.buffer, 1, input.seqbuf.position - 1)

# quality
encoding, input.quality_encodings =
infer_quality_encoding(input.qualbuf.buffer, 1,
input.qualbuf.position - 1,
input.quality_encodings)

decode_quality_string!(encoding, input.qualbuf.buffer,
output.metadata.quality, 1,
input.qualbuf.position - 1)

# reset temporaries for the next run
empty!(input.qualbuf)
empty!(input.seqbuf)
empty!(input.name2buf)
empty!(input.desc2buf)

Ragel.@yield 0
end

break;
end

end

end
@label _out
begin

end

end
end)


text:
struct Table <: Drawable
    cell_fn::Function
    num_rows::Integer
    num_cols::Integer
    line_color::Color
    line_thickness::Integer
    title::Union{Nothing,Label}
    title_right::Union{Nothing,Label}
    header::Union{Nothing,Label}
    header_right::Union{Nothing,Label}
    margins::Margins
    font_size::Integer
    valign::Union{Nothing,Symbol}
    halign::Union{Nothing,Symbol}
end

function Table(
    num_rows::Integer,
    num_cols::Integer,
    cell_fn::Function;
    line_thickness=1,
    line_color="#333",
    title=nothing,
    title_right=nothing,
    header=nothing,
    header_right=nothing,
    margin_top=nothing,
    margin_bottom=nothing,
    margin_left=nothing,
    margin_right=nothing,
    font_size=22,
    valign=nothing,
    halign=nothing,
    )
    Table(
        cell_fn,
        num_rows,
        num_cols,
        Color(line_color),
        line_thickness,
        Label(title),
        Label(title_right),
        Label(header),
        Label(header_right),
        Margins(
            top=margin_top,
            bottom=margin_bottom,
            left=margin_left,
            right=margin_right,
        ),
        font_size,
        valign,
        halign,
    )
end

getnumrows(t::Table) = t.num_rows
getnumcols(t::Table) = t.num_cols
getlinecolor(t::Table) = string(t.line_color.name)

function getwidth(t::Table)
    sum([
        sum(map(col->getcolwidth(t,col), 1:t.num_cols)),
        t.margins.left,
        t.margins.right,
    ])
end

function getheight(t::Table)
    sum([
        sum(map(row->getrowheight(t,row), 1:t.num_rows)),
        t.margins.top,
        t.margins.bottom,
    ])
end

function getcellheights(t::Table, row)
    getheight.(
        map(col->t.cell_fn(row, col), 1:t.num_cols)
    )
end

getwidth(::Nothing) = 0
getheight(::Nothing) = 0

function getcellwidths(t::Table, col)
    getwidth.(
        map(row->t.cell_fn(row, col), 1:t.num_rows)
    )
end

function getcolwidth(t::Table, col::Int)
    getcellwidths(t, col) |> maximum
end

function getrowheight(t::Table, row::Int)
    getcellheights(t, row) |> maximum
end

function getrowheights(t::Table)
    map(row->getrowheight(t, row), 1:t.num_rows)
end
function getcolwidths(t::Table)
    map(col->getcolwidth(t, col), 1:t.num_cols)
end

function Base.show(io::IO, m::MIME"image/svg+xml", t::Table)
    Drawing(getwidth(t), getheight(t), "./output/table.svg")
    origin()
    sethue("#fff")
    box(O, getwidth(t), getheight(t), :fill)
    draw(t)
    finish()
    preview()
    print(io, read("./output/table.svg", String))
end

struct TableHead
    t::Table
    width::Integer
end
getwidth(x::TableHead) = x.width
getheight(x::TableHead) = x.t.margins.top
function draw(x::TableHead, position)
    draw(
        x.t.title,
        position,
        x.width;
        halign=:left
    )
end

struct TableFooter
    t::Table
    width::Integer
end
getwidth(x::TableFooter) = x.width
getheight(x::TableFooter) = x.t.margins.bottom
function draw(x::TableFooter, position)
    draw(
        x.t.header,
        position,
        x.width;
        halign=:left
    )
end

struct TableBody
    t::Table
    width::Integer
    height::Integer
    valign::Union{Nothing,Symbol}
    halign::Union{Nothing,Symbol}
end
getwidth(x::TableBody) = x.width
getheight(x::TableBody) = x.height
function draw(x::TableBody, position)
    sethue("#fff")
    #box(O, getwidth(x), getheight(x), :fill)
    row_heights = getrowheights(x.t)
    col_widths = getcolwidths(x.t)
    lt = Luxor.Table(row_heights, col_widths)
    num_rows = length(row_heights)
    num_cols = length(col_widths)
    valign = x.valign
    halign = x.halign
    for r in 1:num_rows, c in 1:num_cols
        row_height = row_heights[r]
        col_width = col_widths[c]
        cell_size = (col_width, row_height)
        draw(tablecell(x.t.cell_fn(r,c); valign=valign, halign=halign), lt[r,c], cell_size)
    end
end

function TableBody(
    t,
    width,
    height;
    valign=nothing,
    halign=nothing,
    )
    TableBody(
        t,
        width,
        height,
        valign,
        halign,
    )
end

function draw(t::Table)
    width = getwidth(t)
    thead = TableHead(t, width)
    tbody = TableBody(t, width, sum(getrowheights(t)); valign=t.valign, halign=t.halign)
    tfoot = TableFooter(t, width)
    partitions = []
    for partition in [thead,tbody,tfoot]
        if !isnothing(partition)
            push!(partitions, partition)
        end
    end
    row_heights = getheight.(partitions)
    col_widths = [width]
    lt = Luxor.Table(row_heights, col_widths)
    for partition_idx in 1:length(lt)
        partition = partitions[partition_idx]
        position = lt[partition_idx]
        draw(partition, position)
    end
end


text:
"""
A mutable type that bundles multi-component records at multiple receivers.
```julia
data=Data(tgrid, ageom, [:P, :Vx])

```
Here, we initialized records, of `:P` and `:Vx` fields, in time domain for receivers and supersources in `ageom`.

## Indexing
* `data[i]` : the records due to ith supersource of all fields
* `data[i][:P]` : extract field `:P`
* `data[i][:nr]` : number of receivers (same as in `ageom`)
* `data.grid` : returns `tgrid`
As mutable objects in Julia are like containers that might hold different values over time, `data` can be modified.
"""
Data=Array{NamedD{Recs},1}

function Array{NamedD{Recs},1}(grid::StepRangeLen, ageom::AGeom, fields::Vector{Symbol})
	return [NamedD(grid,Recs(ageom[i].nr),fields) for i in 1:length(ageom)]
end


text:
module Promisables

  export Promise, Then;
  export Pending, Rejected, Fulfilled;
  export Resolve, Fulfill, Reject;
  export @pawait;

  abstract type Status end;
  type Fulfilled <: Status end;
  type Pending <: Status end;
  type Rejected <: Status end;

  macro pawait(body)
    return quote
      resultingPromise = $(esc(body));
      c = Channel{Any}(1)
      fin = (v) -> put!(c, 1);
      Then(fin, fin, resultingPromise);
      take!(c);
    end
  end

  mutable struct Promise
    status::Status
    channel::Channel{Any}
    value::Any
    Promise() = begin
      c = Channel{Any}(32)
      new(Pending(), c, nothing)
    end
    Promise(c::Channel{Any}) = begin
      new(Pending(), c, nothing)
    end
  end;

  function Fulfill(p::Promise, value::Any)
    if (typeof(value) == Promise)
      success = (x) -> Fulfill(p, x);
      error = (err) -> Reject(p, err);
      return Then(success, error, value);
    else
      p.value = value;
      p.status = Fulfilled();
      put!(p.channel, value);
    end
    return p;
  end

  function Reject(err::Exception)
    p = Promise();
    Reject(p, err);
    return p;
  end

  function Reject(p::Promise, err::Exception)
    put!(p.channel, err);
    p.status = Rejected();
    return p;
  end

  T1 = T where T<:Function
  T2 = T where T<:Function

  function Then(f::T1, p::Promise)::Promise
    Then(f, identity, p);
  end

  function Then(f::T1, err::T2, p::Promise)::Promise
    if (typeof(p.status) == Fulfilled())
      return f(p.value);
    end
    if (typeof(p.status) == Rejected())
      return err(p.value);
    end
    newChan = Channel{Any}(32)
    np = Promise(newChan);
    @schedule begin
      value = take!(p.channel);
      if (typeof(value) <: Exception)
        value = err(value);
        if (typeof(value) == Promise)
          success = (x) -> Fulfill(np, x);
          error = (nerr) -> Reject(np, nerr);
          Then(success, error, value);
          return;
        end
        Reject(np, value);
        return;
      end
      Fulfill(np, f(value));
    end
    return np;
  end

  function Resolve(value)::Promise
    np = Promise();
    Fulfill(np, value);
    return np;
  end

end


text:
module MomentClosure

using ModelingToolkit
using ModelingToolkit: get_noiseeqs, states, get_eqs

using Catalyst
import Catalyst: species, reactionparams, reactions, speciesmap, paramsmap, numspecies,
				 numreactions, netstoichmat, get_iv, get_ps
using Catalyst: substoichmat, prodstoichmat

using Symbolics: value, var_from_nested_derivative, map_subscripts, hessian, gradient, setmetadata, scalarize
#import Symbolics: degree

using SciMLBase, SciMLBase.EnsembleAnalysis
using DiffEqJump
using Random
using Distributions: Geometric

using SymbolicUtils.Rewriters: Chain, PassThrough, Prewalk, Fixpoint
using SymbolicUtils: Symbolic, Term, Real, expand, simplify, operation,
					 arguments, @rule, @acrule, isnotflat, flatten_term,
					 istree, FnType

using DataStructures: OrderedDict
using TupleTools: sort
using Combinatorics
using Cumulants
using Latexify

using DocStringExtensions

export generate_central_moment_eqs, generate_raw_moment_eqs, bernoulli_moment_eqs,
       propensities, ReactionSystemMod, species, reactionparams, speciesmap, paramsmap,
	   numspecies, numreactions, netstoichmat,
       moment_closure, deterministic_IC, JumpProblem,
	   get_raw_moments, get_central_moments, get_cumulants, get_moments_FSP,
	   linear_mapping_approximation

# needed for ReactionSystemMod definition
export @parameters, @variables

include("reaction_systems.jl")
include("moment_equations.jl")
include("symbolic.jl")
include("moment_convert.jl")
include("stochastic_stoichiometry.jl")
include("central_moment_equations.jl")
include("central_moment_equations_SDE.jl")
include("raw_moment_equations.jl")
include("raw_moment_equations_SDE.jl")
include("closure_methods/closure.jl")
include("bernoulli.jl")
include("utils.jl")
include("closure_methods/zero_closure.jl")
include("closure_methods/normal_closure.jl")
include("closure_methods/log_normal_closure.jl")
include("closure_methods/poisson_closure.jl")
include("closure_methods/gamma_closure.jl")
include("closure_methods/derivative_matching.jl")
include("closure_methods/conditional_gaussian.jl")
include("closure_methods/conditional_derivative_matching.jl")
include("closure_methods/linear_mapping_approximation.jl")
include("stochastic_simulation.jl")

end