text:
using MOCNeutronTransport
@testset "QuadraticSegment_2D" begin
    for F in [Float32, Float64, BigFloat]
        @testset "Constructors" begin
            # Constructor
            x⃗₁ = Point_2D{F}(0, 0)
            x⃗₂ = Point_2D{F}(2, 0)
            x⃗₃ = Point_2D{F}(1, 1)
            q = QuadraticSegment_2D(x⃗₁, x⃗₂, x⃗₃)
            @test q.points == SVector(x⃗₁, x⃗₂, x⃗₃)
        end

        @testset "Methods" begin
            # interpolation
            x⃗₁ = Point_2D{F}(0, 0)
            x⃗₂ = Point_2D{F}(2, 0)
            x⃗₃ = Point_2D{F}(1, 1)
            q = QuadraticSegment_2D(x⃗₁, x⃗₂, x⃗₃)
            for r = LinRange{F}(0, 1, 11)
                @test q(r) ≈ Point_2D{F}(2r, -(2r)^2 + 4r)
            end

            # gradient
            for r = LinRange{F}(0, 1, 11)
                @test ∇(q, r) ≈ Point_2D{F}(2, -(8r) + 4)
            end

            # arclength
            x⃗₁ = Point_2D{F}(0, 0)
            x⃗₂ = Point_2D{F}(2, 0)
            x⃗₃ = Point_2D{F}(1, 0)
            q = QuadraticSegment_2D(x⃗₁, x⃗₂, x⃗₃)
            # straight edge
            @test abs(arclength(q) - 2) < 1.0e-6
            # curved
            x⃗₃ = Point_2D{F}(1, 1)
            q = QuadraticSegment_2D(x⃗₁, x⃗₂, x⃗₃)
            @test abs(arclength(q) - 2.9578857151786138) < 1.0e-6

            # boundingbox
            x⃗₁ = Point_2D{F}(0, 0)
            x⃗₂ = Point_2D{F}(2, 0)
            x⃗₃ = Point_2D{F}(1, 1)
            q = QuadraticSegment_2D(x⃗₁, x⃗₂, x⃗₃)
            bb = boundingbox(q)
            @test bb.xmin ≈ 0
            @test bb.ymin ≈ 0
            @test bb.xmax ≈ 2
            @test bb.ymax ≈ 1
            x⃗₁ = Point_2D{F}(0, 0)
            x⃗₂ = Point_2D{F}(2, 2)
            x⃗₃ = Point_2D{F}(1, 1)
            q = QuadraticSegment_2D(x⃗₁, x⃗₂, x⃗₃)
            bb = boundingbox(q)
            @test bb.xmin ≈ 0
            @test bb.ymin ≈ 0
            @test bb.xmax ≈ 2
            @test bb.ymax ≈ 2
            x⃗₁ = Point_2D{F}(0, 0)
            x⃗₂ = Point_2D{F}(2, 0)
            x⃗₃ = Point_2D{F}(2.1, 1)
            q = QuadraticSegment_2D(x⃗₁, x⃗₂, x⃗₃)
            bb = boundingbox(q)
            @test bb.xmin ≈ 0
            @test bb.ymin ≈ 0
            @test bb.xmax ≈ 2.3272727272727276
            @test bb.ymax ≈ 1

            # isstraight
            x⃗₁ = Point_2D{F}(0, 0)
            x⃗₂ = Point_2D{F}(2, 0)
            x⃗₃ = Point_2D{F}(1, 0)
            q = QuadraticSegment_2D(x⃗₁, x⃗₂, x⃗₃)
            @test isstraight(q)
            x⃗₂ = Point_2D{F}(2, 0.0001)
            q = QuadraticSegment_2D(x⃗₁, x⃗₂, x⃗₃)
            @test !isstraight(q)

            # laplacian
            x⃗₁ = Point_2D{F}(0, 0)
            x⃗₂ = Point_2D{F}(2, 0)
            x⃗₃ = Point_2D{F}(1, 1)
            q = QuadraticSegment_2D(x⃗₁, x⃗₂, x⃗₃)
            p = ∇²(q, 0)
            @test p.x ≈ 0
            @test p.y ≈ -8
            p = ∇²(q, 1)
            @test p.x ≈ 0
            @test p.y ≈ -8

#            # closest_point
#            x⃗₁ = Point_2D{F}(0, 0)
#            x⃗₂ = Point_2D{F}(2, 0)
#            x⃗₃ = Point_2D{F}(1, 1)
#            q = QuadraticSegment_2D(x⃗₁, x⃗₂, x⃗₃)
#            p = Point_2D{F}(1, 1.1)
#            r, p_c = closest_point(p, q)
#            @test x⃗₃ ≈ p_c
#            p = Point_2D{F}(-0.1, 0)
#            r, p_c = closest_point(p, q)
#            @test x⃗₁ ≈ p_c
#
#            # intersect
#            x⃗₁ = Point_2D{F}(0, 0)
#            x⃗₂ = Point_2D{F}(2, 0)
#            x⃗₃ = Point_2D{F}(1, 1)
#            x⃗₄ = Point_2D{F}(1, 0)
#            x⃗₅ = Point_2D{F}(1, 2)
#
#            # 1 intersection
#            q = QuadraticSegment_2D(x⃗₁, x⃗₂, x⃗₃)
#            l = LineSegment_2D(x⃗₄, x⃗₅)
#            npoints, (point1, point2) = intersect(l, q)
#            @test npoints == 1
#            @test point1 ≈ Point_2D{F}(1, 1)
#
#            # 2 intersections
#            x⃗₄ = Point_2D{F}(0, 3//4)
#            x⃗₅ = Point_2D{F}(2, 3//4)
#            l = LineSegment_2D(x⃗₄, x⃗₅)
#            npoints, (point1, point2) = l ∩ q
#            @test npoints == 2
#            @test point1 ≈ Point_2D{F}(1//2, 3//4)
#            @test point2 ≈ Point_2D{F}(3//2, 3//4)
#
#            # 0 intersections
#            x⃗₄ = Point_2D{F}(0, 3)
#            x⃗₅ = Point_2D{F}(2, 3)
#            l = LineSegment_2D(x⃗₄, x⃗₅)
#            npoints, (point1, point2) = intersect(l, q)
#            @test npoints == 0
#            @test point1 ≈ Point_2D{F}(0)
#            @test point2 ≈ Point_2D{F}(0)
#
#            # is_left
#            x⃗₁ = Point_2D{F}(0, 0)
#            x⃗₂ = Point_2D{F}(2, 0)
#            x⃗₃ = Point_2D{F}(1, 1)
#            q = QuadraticSegment_2D(x⃗₁, x⃗₂, x⃗₃)
#            @test !is_left(Point_2D{F}(1, 0), q)
#            @test is_left(Point_2D{F}(1, 2), q)
        end
    end
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
#
# Geometric Sampler
#
# parameters:
# 	(1) p
#
# Note that this interprets p in the same way as the underlying Distributions.Geometric type, i.e.
# values close to 1.0 are unlikely to sample numbers above 0, while values close to 0.0 have samples with a large mean.
# This is the opposite interpretation to the Ruby implementation (and to Boltzmann Dists), and so is NOT the most
# convenient one for modification with the Decay Dist

type GeometricSampler <: DiscreteDistributionSampler
	paramranges::Vector{Tuple{Float64,Float64}}
	distribution::Geometric
	function GeometricSampler(params=Float64[])
		s = new([(0.0, 1.0)])
		setparams!(s, isempty(params) ? [0.5] : params)
		s
	end
end

function setparams!(s::GeometricSampler, params)
	checkparamranges(s, params)
	# parameter of Distribution.Geometric must be in the open interval (0,1), so silently adjust if necessary
	# note nextfloat(0.0), prevfloat(1.0) insufficient adjustment: gives rise to InexactValue error
	p = min(0.99999, max(0.00001, params[1]))
	s.distribution = Geometric(p)
end

getparams(s::GeometricSampler) = [s.distribution.p]

function estimateparams!(s::GeometricSampler, traces)
	samples = extractsamplesfromtraces(s, traces)
	minnumsamples = 1
	if length(samples) >= minnumsamples
		if all(samples.==0)
			# distribution gives error if all samples are 0
			s.distribution = Geometric(0.99999)
		else
			s.distribution = fit(primarydistributiontype(s.distribution), samples)
		end
	end
end

ext: jl


text:
"""
    select_top_taxa(se::SummarizedExperiment, assay_name::String; top_n::Int64 = 10, sel_method::Function = mean)

Sorts rows / features by one specific assay and returns the top elements.

# Arguments
- `se::SummarizedExperiment`: the experiment object of interest.
- `assay_name::String`: the name of the assay to base the evaluation on.
- `top_n::Int64`: the number of top elements that are returned.
- `sel_method::Function`: statistical measure to perform sorting (`mean` by default).
"""
function select_top_taxa(se::SummarizedExperiment, assay_name::String; top_n::Int64 = 10, sel_method::Function = mean)

    features = copy(se.rowdata)
    features[!, Symbol(assay_name)] = [sel_method(row) for row in eachrow(assay(se, assay_name))]

    sort(dropmissing(features, Symbol(assay_name)), Symbol(assay_name), rev = true)[1:top_n, :]

end

"""
    select_top_taxa(se::SummarizedExperiment, assay_names::Vector{String}; top_n::Int64 = 10, sel_method::Function = mean)

Sorts rows / features by multiple assays and returns the top elements.

# Arguments
- `se::SummarizedExperiment`: the experiment object of interest.
- `assay_names::Vector{String}`: a vector of multiple assay names to base the evaluation on.
- `top_n::Int64`: the number of top elements that are returned.
- `sel_method::Function`: statistical measure to perform sorting (`mean` by default).
"""
function select_top_taxa(se::SummarizedExperiment, assay_names::Vector{String}; top_n::Int64 = 10, sel_method::Function = mean)

    features = copy(se.rowdata)

    for ass in assay_names

        features[!, Symbol(ass)] = [sel_method(row) for row in eachrow(assay(se, ass))]

    end

    sort(dropmissing(features, Symbol.(assay_names)), Symbol.(assay_names), rev = true)[1:top_n, :]

end

"""
    dropmissing(se::SummarizedExperiment, assay_name::String)

Drop rows / features with missing elements for a specific assay.

# Arguments
- `se::SummarizedExperiment`: the experiment object of interest.
- `assay_name::String`: the name of the assay whose missing values should be dropped.
"""
function dropmissing(se::SummarizedExperiment, assay_name::String)

    keep_rows = map(x -> sum(ismissing.(x)) == 0, eachrow(assay(se, assay_name)))
    se[keep_rows, :]

end


text:
using Documenter, TetGen

function make_all()
    makedocs(

        sitename="TetGen.jl",
        modules = [TetGen],
        clean = true,
        doctest = false,
        authors = "Simon Danisch, Juergen Fuhrmann",
        repo="https://github.com/JuliaGeometry/TetGen.jl",
        pages=[
            "Home"=>"index.md"
        ]
    )

    if !isinteractive()
        deploydocs(repo = "github.com/JuliaGeometry/TetGen.jl.git")
    end
end

make_all()


text:
abstract type PreallocatedDataNewton{T <: Real, S} <: PreallocatedData{T, S} end

abstract type PreallocatedDataNewtonKrylov{T <: Real, S} <: PreallocatedDataNewton{T, S} end

include("K3Krylov.jl")
include("K3_5Krylov.jl")

abstract type PreallocatedDataNewtonStructured{T <: Real, S} <: PreallocatedDataNewton{T, S} end

include("K3_5Structured.jl")


text:
# # SingleLayerQG decaying 2D turbulence with and without finite Rossby radius of deformation
#
#md # This example can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/singlelayerqg_decaying_barotropic_equivalentbarotropic.ipynb).
#
# We use here the `SingleLayerQG` module to simulate decaying two-dimensional turbulence and
# investigate how does a finite Rossby radius of deformation affects its evolution.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add GeophysicalFlows, Printf, Random, Plots"
# ```

# ## Let's begin
# Let's load `GeophysicalFlows.jl` and some other needed packages.
#
using GeophysicalFlows, Printf, Random, Plots

using GeophysicalFlows: peakedisotropicspectrum
using Random: seed!

# ## Choosing a device: CPU or GPU

dev = CPU()     # Device (CPU/GPU)
nothing # hide

# ## Numerical, domain, and simulation parameters
#
# First, we pick some numerical and physical parameters for our model.

n, L  = 128, 2π             # grid resolution and domain length
deformation_radius = 0.35   # the deformation radius
nothing # hide

## Then we pick the time-stepper parameters
    dt = 1e-2  # timestep
nsteps = 4000  # total number of steps
 nsubs = 20    # number of steps between each plot
nothing # hide

# ## Problem setup
# We initialize two problems by providing a set of keyword arguments to the `Problem` constructor.
# The two problems are otherwise the same, except one has an infinite deformation radius, `prob_bqg`,
# and the other has finite deformation radius, `prob_eqbqg`.
prob_bqg = SingleLayerQG.Problem(dev; nx=n, Lx=L,
                                      dt=dt, stepper="FilteredRK4", aliased_fraction=0)
prob_eqbqg = SingleLayerQG.Problem(dev; nx=n, Lx=L, deformation_radius = deformation_radius,
                                      dt=dt, stepper="FilteredRK4", aliased_fraction=0)
nothing # hide

# ## Setting initial conditions

# For initial condition we construct a relative vorticity with energy most energy around total
# wavenumber ``k_0``.
seed!(1234)
k₀, E₀ = 6, 0.5
∇²ψ₀ = peakedisotropicspectrum(prob_bqg.grid, k₀, E₀, mask=prob_bqg.timestepper.filter)
nothing # hide

# `SingleLayerQG` allows us to set up the initial ``q`` for each problem via `set_q!()` function.
# To initialize both `prob_bqg` and `prob_eqbqg` with the same flow, we first use function
# `SingleLayerQG.streamfunctionfrompv!` to get the streamfunction that corresponds to the
# relative vorticity we computed above. This works in the purely barotropic problem, `prob_bqg`
# since in that case the QGPV is simply the relative vorticity.
∇²ψ₀h = rfft(∇²ψ₀)
ψ₀h = @. 0 * ∇²ψ₀h
SingleLayerQG.streamfunctionfrompv!(ψ₀h, ∇²ψ₀h, prob_bqg.params, prob_bqg.grid)
nothing # hide

# and then use the streamfunction to compute the corresponding ``q_0`` for each problem,
q₀_bqg   = irfft(-prob_bqg.grid.Krsq .* ψ₀h, prob_bqg.grid.nx)
q₀_eqbqg = irfft(-(prob_eqbqg.grid.Krsq .+ 1/prob_eqbqg.params.deformation_radius^2) .* ψ₀h, prob_bqg.grid.nx)
nothing # hide

# Now we can initialize our problems with the same flow.
SingleLayerQG.set_q!(prob_bqg, q₀_bqg)
SingleLayerQG.set_q!(prob_eqbqg, q₀_eqbqg)
nothing # hide

# Let's plot the initial vorticity field for each problem. A function that returns relative
# vorticity from each problem's state variable will prove useful. Note that when plotting, we
# decorate the variable to be plotted with `Array()` to make sure it is brought back on the
# CPU when the variable lives on the GPU.
relativevorticity(prob) = irfft(-prob.grid.Krsq .* prob.vars.ψh, prob.grid.nx)

x, y = prob_bqg.grid.x, prob_bqg.grid.y

p_bqg = heatmap(x, y, Array(relativevorticity(prob_bqg)'),
         aspectratio = 1,
                   c = :balance,
                clim = (-40, 40),
               xlims = (-L/2, L/2),
               ylims = (-L/2, L/2),
              xticks = -3:3,
              yticks = -3:3,
              xlabel = "x",
              ylabel = "y",
               title = "barotropic\n ∇²ψ, t=" * @sprintf("%.2f", prob_bqg.clock.t),
          framestyle = :box)

p_eqbqg = heatmap(x, y, Array(relativevorticity(prob_eqbqg)'),
         aspectratio = 1,
                   c = :balance,
                clim = (-40, 40),
               xlims = (-L/2, L/2),
               ylims = (-L/2, L/2),
              xticks = -3:3,
              yticks = -3:3,
              xlabel = "x",
              ylabel = "y",
               title = "equivalent barotropic; deformation radius: " * @sprintf("%.2f", prob_eqbqg.params.deformation_radius) * "\n ∇²ψ, t=" * @sprintf("%.2f", prob_eqbqg.clock.t),
          framestyle = :box)

l = @layout Plots.grid(1, 2)
p = plot(p_bqg, p_eqbqg, layout = l, size = (800, 380))

# ## Time-stepping the `Problem` forward

# Now we time-step both problems forward and animate the relative vorticity in each case.

startwalltime = time()

cfl(prob) = prob.clock.dt * maximum([maximum(prob.vars.u) / prob.grid.dx, maximum(prob.vars.v) / prob.grid.dy])

anim = @animate for j = 0:Int(nsteps/nsubs)
  if j % (1000 / nsubs) == 0
    log_bqg = @sprintf("barotropic; step: %04d, t: %d, cfl: %.2f, walltime: %.2f min",
        prob_bqg.clock.step, prob_bqg.clock.t, cfl(prob_bqg), (time()-startwalltime)/60)
    println(log_bqg)

    log_eqbqg = @sprintf("equivalent barotropic; step: %04d, t: %d, cfl: %.2f, walltime: %.2f min",
        prob_eqbqg.clock.step, prob_eqbqg.clock.t, cfl(prob_eqbqg), (time()-startwalltime)/60)
    println(log_eqbqg)
  end

  p[1][1][:z] = Array(relativevorticity(prob_bqg))
  p[1][:title] = "barotropic\n ∇²ψ, t=" * @sprintf("%.2f", prob_bqg.clock.t)
  p[2][1][:z] = Array(relativevorticity(prob_eqbqg))
  p[2][:title] = "equivalent barotropic; deformation radius: " * @sprintf("%.2f", prob_eqbqg.params.deformation_radius) * "\n ∇²ψ, t=" * @sprintf("%.2f", prob_eqbqg.clock.t)

  stepforward!(prob_bqg, nsubs)
  SingleLayerQG.updatevars!(prob_bqg)

  stepforward!(prob_eqbqg, nsubs)
  SingleLayerQG.updatevars!(prob_eqbqg)
end

mp4(anim, "singlelayerqg_barotropic_equivalentbarotropic.mp4", fps=18)


text:
InputFolder = './Images/HeLaMCF10A/';
OutputFolder = './Results/Images/HeLaMCF10A/';

@iA = '*.tif';

@fxg_mLapThrBinWatTiles [iA] > [L];
params.GRad = 7;
params.Thr = 0.02;
params.GaussianD = 5;
params.MinArea = 20;
/endf

/show iA > L;
/keep L > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
using AbstractTrees
using Base.Iterators: flatten, map as imap, filter as ifilter, product as iproduct
using Random: shuffle!

struct Tree{T}
    root::T
    subtrees
    Tree(el::T, subtrees=T[]) where T = new{T}(el, subtrees)
end

root(t::Tree) = t.root
subtrees(t::Tree) = t.subtrees

Base.show(io::IO, t::Tree) = print(io, "Tree(", t.root, ')')

AbstractTrees.children(t::Tree) = collect(subtrees(t))
Base.eltype(::Type{<:Tree{T}}) where {T} = T

unfold(f) = Base.Fix1(unfold, f)
function unfold(f, t::T) where {T}
    Tree(t, imap(unfold(f), f(t)))
end

# recursively shrinks and creates a tree
shrinkMap(root) = imap(t -> Tree(t, shrinkMap(t)), shrink(root))

function interleave(trees::Union{<:NTuple{N,Tree{T}},Vector{Tree{T}},<:Tuple}) where {T,N}
    # the root of interleaved trees is just all individual roots
    nRoot = map(root, trees)
    return Tree(nRoot, shrinkMap(nRoot)) # TODO: isn't this just `unfold`?
end

"""
    filter(p, t::Tree[, trim=false])

Filters `t` lazily such that all elements contained fulfill the predicate `p`, i.e. all elements for which `p` is `false` are removed.

The first-level subtrees produced by the returned tree will have unique roots amongst each other.

`trim` controls whether subtrees are removed completely when the root doesn't
fulfill the predicate or whether only that root should be skipped, still trying to
shrink its subtrees. This trades performance (less shrinks to check) for quality
(fewer/less diverse shrink values tried).
"""
function Base.filter(f, t::Tree{T}, trim=false) where {T}
    r = root(t)
    _filter(x) = filter(f, x, trim)
    flat = Flatten{Tree{T}}(imap(_filter, subtrees(t)))

    if f(r)
        Flatten{Tree{T}}(Ref(Ref(Tree(r, flat))))
    else
        if !trim
            flat
        else
            Flatten{Tree{T}}()
        end
    end
end

"""
    map(f, t::Tree)

Maps `f` lazily over all elements in `t`, producing a new tree.

The first-level subtrees produced by the returned tree will have unique roots amongst each other.
"""
Base.map(f, t::Tree) = imap(f, t)

# lazy mapping by default
_map(f) = Base.Fix1(imap, f)
function Base.Iterators.map(f, t::Tree{T}) where {T}
    r = root(t)
    lazySubtrees = iunique(imap(_map, subtrees(t)))
    return Tree(f(r), lazySubtrees)
end


text:
length(ARGS)==2 || error("Usage: julia zn11pvc.jl datafile gamma")
using HDF5,JLD,KUnet
KUnet.gpu(true)
@date @load ARGS[1]
g0 = eval(parse(ARGS[2]))

net = Layer[KPerceptron(size(ytrn,1), KUnet.kgauss, [g0])]
for epoch=1:20
    @show epoch
    @date train(net, xtrn, ytrn)
    @show size(net[1].s)
    @date KUnet.uniq!(net[1])
    @show size(net[1].s)
    @date zdev = predict(net,xdev)
    @date ztst = predict(net,xtst)
    @show (epoch, size(net[1].s), accuracy(ydev,zdev), accuracy(ytst,ztst))
end

# Use this to optimize gamma:
# isdefined(:xtrn) || (@date @load "zn11cpv.jld")
# niters=2000
# nbatch=100
# nc = size(ytrn,1)
# for g0=0.10:0.01:0.15
# #for g=2:.5:4; g0 = 2f0^-g
#     println("")
#     @show g0
#     net = Layer[KPerceptron(nc, KUnet.kgauss, [g0])]
#     @date train(net, xtrn, ytrn; iters=niters, batch=nbatch)
#     @show size(net[1].s)
#     @date y = predict(net,xdev)
#     @show accuracy(ydev,y)
# end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Quadratic decision function #
# ========================== #
struct QuadraticDecisionFunction{T} <: MOI.AbstractScalarFunction
    variable_part::MOI.ScalarQuadraticFunction{T}
    decision_part::MOI.ScalarQuadraticFunction{T}
    cross_terms::MOI.ScalarQuadraticFunction{T}
end

# Base overrides #
# ========================== #
function Base.copy(f::F) where F <: QuadraticDecisionFunction
    return F(copy(f.variable_part),
             copy(f.decision_part),
             copy(f.cross_terms))
end
function Base.convert(::Type{QuadraticDecisionFunction{T}}, α::T) where T
    return QuadraticDecisionFunction{T}(
        convert(MOI.ScalarQuadraticFunction{T}, α),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end
function Base.convert(::Type{QuadraticDecisionFunction{T}},
                      f::MOI.SingleVariable) where T
    return QuadraticDecisionFunction{T}(
        convert(MOI.ScalarQuadraticFunction{T}, f),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end
function Base.convert(::Type{QuadraticDecisionFunction{T}},
                      f::SingleDecision) where T
    return QuadraticDecisionFunction{T}(
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        MOI.ScalarQuadraticFunction{T}(
            [MOI.ScalarAffineTerm(one(T), f.decision)],
            MOI.ScalarQuadraticTerm{T}[],
            zero(T),
        ),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end
function Base.convert(::Type{QuadraticDecisionFunction{T}},
                      f::MOI.ScalarAffineFunction{T}) where T
    return QuadraticDecisionFunction{T}(
        convert(MOI.ScalarQuadraticFunction{T}, f),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end
function Base.convert(::Type{QuadraticDecisionFunction{T}},
                      f::MOI.ScalarQuadraticFunction{T}) where T
    return QuadraticDecisionFunction{T}(
        f,
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end
function Base.convert(::Type{QuadraticDecisionFunction{T}},
                      f::AffineDecisionFunction{T}) where T
    return QuadraticDecisionFunction{T}(
        convert(MOI.ScalarQuadraticFunction{T}, f.variable_part),
        convert(MOI.ScalarQuadraticFunction{T}, f.decision_part),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end
function Base.convert(::Type{MOI.SingleVariable},
                      f::QuadraticDecisionFunction{T}) where T
    return convert(MOI.SingleVariable, convert(AffineDecisionFunction{T}, f))
end
function Base.convert(::Type{SingleDecision},
                      f::QuadraticDecisionFunction{T}) where T
    return convert(SingleDecision, convert(AffineDecisionFunction{T}, f))
end

function Base.convert(::Type{AffineDecisionFunction{T}},
                      f::QuadraticDecisionFunction{T}) where T
    if !Base.isempty(f.cross_terms.quadratic_terms)
        throw(InexactError(:convert, AffineDecisionFunction{T}, f))
    end
    return AffineDecisionFunction{T}(
        convert(MOI.ScalarAffineFunction{T}, f.variable_part),
        convert(MOI.ScalarAffineFunction{T}, f.decision_part))
end

function Base.zero(F::Type{QuadraticDecisionFunction{T}}) where T
    return convert(F, zero(T))
end

function Base.iszero(f::QuadraticDecisionFunction)
    return iszero(MOI.constant(f)) && MOIU._is_constant(MOIU.canonical(f))
end

function Base.isone(f::QuadraticDecisionFunction)
    return isone(MOI.constant(f)) && MOIU._is_constant(MOIU.canonical(f))
end

function Base.one(F::Type{QuadraticDecisionFunction{T}}) where T
    return convert(F, one(T))
end

function MOI._dicts(f::QuadraticDecisionFunction)
    return (MOI.sum_dict(MOI.term_pair.(f.variable_part.affine_terms)),
            MOI.sum_dict(MOI.term_pair.(f.variable_part.quadratic_terms)),
            MOI.sum_dict(MOI.term_pair.(f.decision_part.affine_terms)),
            MOI.sum_dict(MOI.term_pair.(f.decision_part.quadratic_terms)),
            MOI.sum_dict(MOI.term_pair.(f.cross_terms.quadratic_terms)))
end

function Base.isapprox(f::QuadraticDecisionFunction, g::QuadraticDecisionFunction; kwargs...)
    return isapprox(MOI.constant(f), MOI.constant(g); kwargs...) && all(
        MOI.dict_compare.(
            MOI._dicts(f),
            MOI._dicts(g),
            (α, β) -> isapprox(α, β; kwargs...),
        ),
    )
end

# JuMP overrides #
# ========================== #
function QuadraticDecisionFunction(quad::DQE)
    JuMP._assert_isfinite(quad)
    # Variable part
    variable_part = JuMP.moi_function(quad.variables)
    # Decision part
    decision_affine_terms =  Vector{MOI.ScalarAffineTerm{Float64}}()
    for (coef, dvar) in linear_terms(quad.decisions.aff)
        # Any fixed decision value is set in the decision bridge
        push!(decision_affine_terms, MOI.ScalarAffineTerm(coef, index(dvar)))
    end
    decision_quad_terms = Vector{MOI.ScalarQuadraticTerm{Float64}}()
    for t in quad_terms(quad.decisions)
        # Any fixed decision value is set in the decision bridge
        push!(decision_quad_terms, JuMP._moi_quadratic_term(t))
    end
    # Cross terms
    cross_terms = Vector{MOI.ScalarQuadraticTerm{Float64}}()
    for (cross_term, coeff) in quad.cross_terms
        # Any fixed decision value is set in the decision bridge
        push!(cross_terms,
              MOI.ScalarQuadraticTerm(coeff, index(cross_term.decision), index(cross_term.variable)))
    end
    # Return QuadraticDecisionFunction with QuadraticPart
    return QuadraticDecisionFunction(
        variable_part,
        MOI.ScalarQuadraticFunction(decision_affine_terms,
                                    decision_quad_terms,
                                    0.0),
        MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[],
                                    cross_terms,
                                    0.0))
end
JuMP.moi_function(quad::DQE) = QuadraticDecisionFunction(quad)
function JuMP.moi_function_type(::Type{DecisionQuadExpr{T}}) where T
    return QuadraticDecisionFunction{T}
end

is_decision_type(::Type{<:QuadraticDecisionFunction}) = true

QuadraticDecisionFunction(quad::_DQE) = QuadraticDecisionFunction(convert(DQE, quad))

function _DecisionQuadExpr(m::Model, f::MOI.ScalarQuadraticFunction)
    quad = _DQE(_DecisionAffExpr(m, MOI.ScalarAffineFunction(f.affine_terms, 0.0)))
    for t in f.quadratic_terms
        v1 = t.variable_index_1
        v2 = t.variable_index_2
        coef = t.coefficient
        if v1 == v2
            coef /= 2
        end
        JuMP.add_to_expression!(quad, coef, DecisionRef(m, v1), DecisionRef(m, v2))
    end
    # There should not be any constants in the decision terms
    return quad
end

function DQE(model::Model, f::QuadraticDecisionFunction{T}) where T
    # Decision part
    cross_terms = OrderedDict{DecisionCrossTerm, Float64}()
    for term in f.cross_terms.quadratic_terms
        cross_term = DecisionCrossTerm(
            DecisionRef(model, term.variable_index_1),
            VariableRef(model, term.variable_index_2))
        JuMP._add_or_set!(cross_terms, cross_term, term.coefficient)
    end
    return DQE(QuadExpr(model,
                        convert(MOI.ScalarQuadraticFunction{T}, f.variable_part)),
               _DecisionQuadExpr(model,
                                convert(MOI.ScalarQuadraticFunction{T}, f.decision_part)),
               cross_terms)
end

function JuMP.jump_function_type(::Model,
                                 ::Type{QuadraticDecisionFunction{T}}) where T
    return DecisionQuadExpr{T}
end
function JuMP.jump_function(model::Model, f::QuadraticDecisionFunction{T}) where T
    return DecisionQuadExpr{T}(model, f)
end

# MOI Function interface #
# ========================== #
MOI.constant(f::QuadraticDecisionFunction) =
    MOI.constant(f.variable_part)
MOI.constant(f::QuadraticDecisionFunction, T::DataType) = MOI.constant(f)

function MOIU.eval_variables(varval::Function, f::QuadraticDecisionFunction)
    var_value = MOIU.eval_variables(varval, f.variable_part)
    dvar_value = MOIU.eval_variables(varval, f.decision_part)
    cross_value = MOIU.eval_variables(varval, f.cross_terms)
    return var_value + dvar_value + cross_value
end

function MOIU.map_indices(index_map::Function, f::QuadraticDecisionFunction{T}) where T
    return typeof(f)(MOIU.map_indices(index_map, f.variable_part),
                     MOIU.map_indices(index_map, f.decision_part),
                     MOIU.map_indices(index_map, f.cross_terms))
end

function MOIU.substitute_variables(variable_map::Function, f::QuadraticDecisionFunction{T}) where T
    g = QuadraticDecisionFunction(
        convert(MOI.ScalarQuadraticFunction{T}, MOI.constant(f.variable_part)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
    # Substitute variables
    for term in f.variable_part.affine_terms
        func::AffineDecisionFunction{T} = variable_map(term.variable_index)
        new_term = MOIU.operate(*, T, term.coefficient, func)::AffineDecisionFunction{T}
        MOIU.operate!(+, T, g, new_term)::QuadraticDecisionFunction{T}
    end
    for term in f.variable_part.quadratic_terms
        f1::AffineDecisionFunction{T} = variable_map(term.variable_index_1)
        f2::AffineDecisionFunction{T} = variable_map(term.variable_index_2)
        f12 = MOIU.operate(*, T, f1, f2)::QuadraticDecisionFunction{T}
        coeff = term.coefficient
        if term.variable_index_1 == term.variable_index_2
            coeff /= 2
        end
        new_term = MOIU.operate!(*, T, f12, coeff)::QuadraticDecisionFunction{T}
        MOIU.operate!(+, T, g, new_term)::QuadraticDecisionFunction{T}
    end
    # Substitute decisions
    for term in f.decision_part.affine_terms
        func::AffineDecisionFunction{T} = variable_map(term.variable_index)
        new_term = MOIU.operate(*, T, term.coefficient, func)::AffineDecisionFunction{T}
        MOIU.operate!(+, T, g, new_term)::QuadraticDecisionFunction{T}
    end
    for term in f.decision_part.quadratic_terms
        f1::AffineDecisionFunction{T} = variable_map(term.variable_index_1)
        f2::AffineDecisionFunction{T} = variable_map(term.variable_index_2)
        f12 = MOIU.operate(*, T, f1, f2)::QuadraticDecisionFunction{T}
        coeff = term.coefficient
        if term.variable_index_1 == term.variable_index_2
            coeff /= 2
        end
        new_term = MOIU.operate!(*, T, f12, coeff)::QuadraticDecisionFunction{T}
        MOIU.operate!(+, T, g, new_term)::QuadraticDecisionFunction{T}
    end
    # Substitute decisions in cross terms
    for term in f.cross_terms.quadratic_terms
        f1::AffineDecisionFunction{T} = variable_map(term.variable_index_1)
        f2::AffineDecisionFunction{T} = variable_map(term.variable_index_2)
        f12 = MOIU.operate(*, T, f1, f2)::QuadraticDecisionFunction{T}
        coeff = term.coefficient
        if term.variable_index_1 == term.variable_index_2
            coeff /= 2
        end
        new_term = MOIU.operate!(*, T, f12, coeff)::QuadraticDecisionFunction{T}
        MOIU.operate!(+, T, g, new_term)::QuadraticDecisionFunction{T}
    end
    return g
end

function MOIU.is_canonical(f::QuadraticDecisionFunction)
    return MOIU.is_canonical(f.variable_part) &&
        MOIU.is_canonical(f.decision_part) &&
        MOIU.is_canonical(f.cross_terms)
end

MOIU.canonical(f::QuadraticDecisionFunction) = MOIU.canonicalize!(copy(f))

function MOIU.canonicalize!(f::QuadraticDecisionFunction)
    MOIU.canonicalize!(f.variable_part)
    MOIU.canonicalize!(f.decision_part)
    MOIU.canonicalize!(f.cross_terms)
    return f
end

function MOIU._is_constant(f::QuadraticDecisionFunction)
    return MOIU._is_constant(f.variable_part) &&
        MOIU._is_constant(f.decision_part) &&
        MOIU._is_constant(f.cross_terms)
end

function MOIU.all_coefficients(p::Function, f::QuadraticDecisionFunction)
    return MOIU.all_coefficients(p, f.variable_part) &&
        MOIU.all_coefficients(p, f.decision_part) &&
        MOIU.all_coefficients(p, f.cross_terms)
end

function MOIU.isapprox_zero(f::QuadraticDecisionFunction, tol)
    return MOIU.isapprox_zero(f.variable_part, tol) &&
        MOIU.isapprox_zero(f.decision_part, tol) &&
        MOIU.isapprox_zero(f.cross_terms, tol)
end

function MOIU.filter_variables(keep::Function, f::QuadraticDecisionFunction)
    return typeof(f)(MOIU.filter_variables(keep, f.variable_part),
                     MOIU.filter_variables(keep, f.decision_part),
                     MOIU.filter_variables(keep, f.cross_terms))
end

function MOIU.operate_coefficients(f, func::QuadraticDecisionFunction)
    return typeof(func)(MOIU.operate_coefficients(f, func.variable_part),
                        MOIU.operate_coefficients(f, func.decision_part),
                        MOIU.operate_coefficients(f, func.cross_terms))
end

function add_term!(terms::Vector{MOI.ScalarQuadraticTerm{T}},
                   term::MOI.ScalarQuadraticTerm{T}) where T
    index_1 = term.variable_index_1
    index_2 = term.variable_index_2
    coefficient = term.coefficient
    i = something(findfirst(t -> t.variable_index_1 == index_1 &&
                            t.variable_index_2 == index_2,
                            terms), 0)
    if iszero(i)
        if !iszero(coefficient)
            # The variable was not already included in the terms
            push!(terms, MOI.ScalarQuadraticTerm(coefficient, index_1, index_2))
        end
    else
        # Variable already included, increment the coefficient
        new_coeff = terms[i].coefficient + coefficient
        if iszero(new_coeff)
            deleteat!(terms, i)
        else
            terms[i] = MOI.ScalarQuadraticTerm(new_coeff, index_1, index_2)
        end
    end
    return nothing
end

function remove_term!(terms::Vector{MOI.ScalarQuadraticTerm{T}},
                      term::MOI.ScalarQuadraticTerm{T}) where T
    index_1 = term.variable_index_1
    index_2 = term.variable_index_2
    coefficient = term.coefficient
    i = something(findfirst(t -> t.variable_index_1 == index_1 &&
                            t.variable_index_2 == index_2,
                            terms), 0)
    if iszero(i) && !iszero(coefficient)
        # The term was not already included in the terms
        push!(terms, MOI.ScalarQuadraticTerm(-coefficient, index_1, index_2))
    else
        # Term already included, increment the coefficient
        new_coeff = terms[i].coefficient - coefficient
        if iszero(new_coeff)
            deleteat!(terms, i)
        else
            terms[i] = MOI.ScalarQuadraticTerm(new_coeff, index_1, index_2)
        end
    end
    return nothing
end


text:
using Test
using Flux: flatten

@testset "helpers" begin
  @testset "flatten" begin
    x = randn(Float32, 10, 10, 3, 2)
    @test size(flatten(x)) == (300, 2)
  end
end


text:
include(joinpath(@__DIR__, "..", "src/dynamics", "flamingo", "visuals.jl"))
T = Float64
vis = Visualizer()
render(vis)
open(vis)

# get model
model_sim = deepcopy(flamingo_piecewise)

plot(x, model_sim.env.surf.(x))
dir = joinpath(module_dir(), "src/dynamics/flamingo")

path_base = joinpath(dir, "dynamics/base.jld2")
path_dyn = joinpath(dir, "dynamics/dynamics.jld2")
path_res = joinpath(dir, "piecewise/residual.jld2")
path_jac = joinpath(dir, "piecewise/sparse_jacobians.jld2")
path_linearized = joinpath(dir, "piecewise/linearized.jld2")

instantiate_base!(model_sim, path_base)

expr_dyn = generate_dynamics_expressions(model_sim, derivs = true)
save_expressions(expr_dyn, path_dyn, overwrite=true)
instantiate_dynamics!(model_sim, path_dyn, derivs = true)

expr_res, rz_sp, rθ_sp = generate_residual_expressions(model_sim, jacobians = :approx)
save_expressions(expr_res, path_res, overwrite=true)
@save path_jac rz_sp rθ_sp
instantiate_residual!(model_sim, path_res, jacobians = :approx)
model_sim.spa.rz_sp = copy(rz_sp)
model_sim.spa.rθ_sp = copy(rθ_sp)

model = get_model("flamingo", surf="flat")
nq = model.dim.q
nu = model.dim.u
nc = model.dim.c
nb = model.dim.b
nd = nq + nc + nb
nr = nq + nu + nc + nb + nd
nz = num_var(model)
nθ = num_data(model)

# get trajectory
# ref_traj = get_trajectory("flamingo", "gait1", load_type=:split_traj_alt, model=model)
# ref_traj = get_trajectory("flamingo", "gait_forward_36_3", load_type=:split_traj_alt, model=model)
ref_traj = get_trajectory("flamingo", "gait_forward_36_4", load_type=:split_traj_alt, model=model)

H = ref_traj.H
h = ref_traj.h
N_sample = 5
H_mpc = 15
h_sim = h / N_sample
H_sim = 2500#15000

# barrier parameter
κ_mpc = 1.0e-4

obj = TrackingVelocityObjective(H_mpc, model.dim,
    v = [Diagonal(1e-3 * [1e0,1,1e4,1,1,1,1,1e4,1e4]) for t = 1:H_mpc],
    q = [Diagonal(1e-1 * [3e2, 1e-6, 3e2, 1, 1, 1, 1, 0.1, 0.1]) for t = 1:H_mpc],
    u = [Diagonal(3e-1 * [0.1; 0.1; 0.3; 0.3; ones(nu-6); 2; 2]) for t = 1:H_mpc],
    γ = [Diagonal(1.0e-100 * ones(model.dim.c)) for t = 1:H_mpc],
    b = [Diagonal(1.0e-100 * ones(model.dim.b)) for t = 1:H_mpc])

p = linearized_mpc_policy(ref_traj, model, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
    n_opts = NewtonOptions(
        r_tol = 3e-4,
        max_iter = 5),
    mpc_opts = LinearizedMPCOptions(
        # live_plotting=true,
        altitude_update = true,
        altitude_impact_threshold = 0.02,
        altitude_verbose = true,
        )
    )

q1_ref = copy(ref_traj.q[2])
q0_ref = copy(ref_traj.q[1])
q1_sim = SVector{model.dim.q}(q1_ref)
q0_sim = SVector{model.dim.q}(copy(q1_sim - (q1_ref - q0_ref) / N_sample))
@assert norm((q1_sim - q0_sim) / h_sim - (q1_ref - q0_ref) / h) < 1.0e-8

sim = simulator(model_sim, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
    ip_opts = InteriorPointOptions(
        r_tol = 1.0e-8,
        κ_init = 1.0e-8,
        κ_tol = 2.0e-8),
    sim_opts = SimulatorOptions(warmstart = true)
    )

@time status = simulate!(sim)

# save trajectory
@save joinpath(module_dir(), "src/dynamics/flamingo/simulations/piecewise.jld2") sim
@load joinpath(module_dir(), "src/dynamics/flamingo/simulations/piecewise.jld2") sim

l = 9
lu = 1
plt = plot(layout=(3,1), legend=false)
plot!(plt[1,1], hcat(Vector.(vcat([fill(ref_traj.q[i], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[1,1], hcat(Vector.([q[l:l] for q in sim.traj.q])...)', color=:blue, linewidth=1.0)
plot!(plt[2,1], hcat(Vector.(vcat([fill(ref_traj.u[i][lu:lu], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[3,1], hcat(Vector.(vcat([fill(ref_traj.γ[i][1:nc], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[2,1], hcat(Vector.([u[lu:lu] for u in sim.traj.u]*N_sample)...)', color=:blue, linewidth=1.0)
# plot!(plt[3,1], hcat(Vector.([γ[1:nc] for γ in sim.traj.γ]*N_sample)...)', color=:blue, linewidth=1.0)
# plot!(plt[3,1], hcat(Vector.([b[1:nb] for b in sim.traj.b]*N_sample)...)', color=:red, linewidth=1.0)

plot_lines!(vis, model, sim.traj.q[1:N_sample:end], offset=-0.01)
plot_surface!(vis, model_sim.env, xlims=[-1, 9])
anim = visualize_meshrobot!(vis, model_sim, sim.traj, sample=10)
anim = visualize_force!(vis, model_sim, sim.traj, anim=anim, h=h_sim, sample=10)
open(vis)

filename = "flamingo_40_steps_sine"
MeshCat.convert_frames_to_video(
    "/home/simon/Downloads/$filename.tar",
    "/home/simon/Documents/$filename.mp4", overwrite=true)

convert_video_to_gif(
    "/home/simon/Documents/$filename.mp4",
    "/home/simon/Documents/$filename.gif", overwrite=true)

# Ghost
flamingo_ghost!(vis, sim, piecewise_smoothed)

# Animation
anim, shift_traj = flamingo_animation!(vis, sim, piecewise_smoothed)
anim = visualize_force!(vis, sim.model, shift_traj, anim=anim, h=h_sim, sample=10)


text:
pydict_setitem(x::Py, k, v) = errcheck(@autopy k v C.PyDict_SetItem(getptr(x), getptr(k_), getptr(v_)))

function pydict_fromiter(kvs)
    ans = pydict()
    for (k, v) in kvs
        pydict_setitem(ans, k, v)
    end
    return ans
end

function pystrdict_fromiter(kvs)
    ans = pydict()
    for (k, v) in kvs
        pydict_setitem(ans, string(k), v)
    end
    return ans
end

"""
    pydict(x)
    pydict(; x...)

Convert `x` to a Python `dict`. In the second form, the keys are strings.

If `x` is a Python object, this is equivalent to `dict(x)` in Python.
Otherwise `x` must iterate over key-value pairs.
"""
pydict(; kwargs...) = isempty(kwargs) ? pynew(errcheck(C.PyDict_New())) : pystrdict_fromiter(kwargs)
pydict(x) = ispy(x) ? pybuiltins.dict(x) : pydict_fromiter(x)
pydict(x::NamedTuple) = pydict(; x...)
export pydict


text:
### Response Parsing

immutable ResponseParser
    parser::Parser
    settings::ParserSettings

    function ResponseParser(r)
        parser = Parser()
        parser.data = r
        http_parser_init(parser,false)
        settings = ParserSettings(on_message_begin_cb, on_url_cb,
                          on_status_complete_cb, on_header_field_cb,
                          on_header_value_cb, on_headers_complete_cb,
                          on_body_cb, on_message_complete_cb)

        new(parser, settings)
    end
end

pd(p::Ptr{Parser}) = (unsafe_load(p).data)::ResponseStream

# All the `HttpParser` callbacks to be run in C land
# Each one adds data to the `RequestStream` until it is complete
#
function on_message_begin(parser)
    pd(parser).state = OnMessageBegin
    notify(pd(parser).state_change)
    return 0
end

function on_status_complete(parser)
    response_stream = pd(parser)
    response_stream.response.status = (unsafe_load(parser)).status_code
    response_stream.state = StatusComplete
    notify(response_stream.state_change)
    return 0
end

function on_header_field(parser, at, len)
    response_stream = pd(parser)
    header = unsafe_string(convert(Ptr{UInt8}, at))
    header_field = header[1:len]
    if response_stream.state == OnHeaderField
        field = string(get(response_stream.current_header, header_field))
        response_stream.current_header = Nullable(field)
    else
        response_stream.current_header = Nullable(header_field)
    end
    response_stream.state = OnHeaderField
    notify(response_stream.state_change)
    return 0
end

function parse_cookies!(response, cookie_strings)
    for cookie_str in split(cookie_strings, '\0')
        isempty(cookie_str) && continue
        parts = split(cookie_str, ';')
        isempty(parts) && continue
        nameval = split(parts[1], '=', limit=2)
        length(nameval)==2 || continue
        name, value = nameval
        c = Cookie(strip(name), strip(value))
        for part in parts[2:end]
            nameval = split(part, '=', limit=2)
            if length(nameval)==2
                name, value = nameval
                c.attrs[strip(name)] = strip(value)
            else
                c.attrs[strip(nameval[1])] = Compat.String("")
            end
        end
        response.cookies[c.name] = c
    end
    response
end

const is_set_cookie = r"set-cookie"i

function on_header_value(parser, at, len)
    response_stream = pd(parser)
    resp = response_stream.response
    s = unsafe_string(convert(Ptr{UInt8}, at), Int(len))
    current_header = get(response_stream.current_header)
    if response_stream.state == OnHeaderValue
        if is_set_cookie(current_header)
            write(response_stream.cookie_buffer, s)
        else
            resp.headers[current_header] = string(resp.headers[current_header], s)
        end
    else
        if is_set_cookie(current_header)
            write(response_stream.cookie_buffer, '\0', s)
            # maybe_cookie = parse_set_cookie(s)
            # if !isnull(maybe_cookie)
            #     cookie = get(maybe_cookie)
            #     response_stream.response.cookies[cookie.name] = cookie
            # end
        else
            response_stream.response.headers[current_header] = s
        end
    end
    # response_stream.current_header = Nullable()
    response_stream.state = OnHeaderValue
    notify(response_stream.state_change)
    return 0
end

function on_headers_complete(parser)
    response_stream = pd(parser)
    response = response_stream.response
    p = unsafe_load(parser)
    # get first two bits of p.type_and_flags
    ptype = p.type_and_flags & 0x03
    if ptype == 0
        response.method = http_method_str(convert(Int, p.method))
    elseif ptype == 1
        response.headers["status_code"] = string(convert(Int, p.status_code))
    end
    response.headers["http_major"] = string(convert(Int, p.http_major))
    response.headers["http_minor"] = string(convert(Int, p.http_minor))
    response.headers["Keep-Alive"] = string(http_should_keep_alive(parser))
    parse_cookies!(response, Compat.String(take!(response_stream.cookie_buffer)))

    response_stream.state = HeadersDone
    notify(response_stream.state_change)

    # From https://github.com/nodejs/http-parser/blob/master/http_parser.h
    # Comment starting line 72
    # On a HEAD method return 1 instead of 0 to indicate that the parser
    # should not expect a body.
    method = get( response.request ).method
    if method ∈ ("HEAD", "CONNECT")
        return 1  # Signal HTTP parser to not expect a body
    else
        return 0
    end
end

function on_body(parser, at, len)
    response_stream = pd(parser)
    append!(response_stream.buffer.data, unsafe_wrap(Array, convert(Ptr{UInt8}, at), (len,)))
    response_stream.buffer.size = length(response_stream.buffer.data)
    response_stream.state = OnBody
    notify(response_stream.state_change)
    return 0
end

function on_message_complete(parser)
    response_stream = pd(parser)
    if upgrade(unsafe_load(parser))
        response_stream.state = UpgradeConnection
    else
        response_stream.state = BodyDone
    end
    notify(response_stream.state_change)
    return 0
end

# Passes `request_data` into `parser`
function add_data(parser::ResponseParser, request_data)
    http_parser_execute(parser.parser, parser.settings, request_data)
end

# Datatype Tuples for the different `cfunction` signatures used by `HttpParser`
const HTTP_CB      = (Int, (Ptr{Parser},))
const HTTP_DATA_CB = (Int, (Ptr{Parser}, Ptr{Cchar}, Csize_t,))

function __init_parsing__()
    # Turn all the callbacks into C callable functions.
    global const on_message_begin_cb = cfunction(on_message_begin, HTTP_CB...)
    global const on_url_cb = C_NULL # callback valid only for Server, not for Request
    global const on_status_complete_cb = cfunction(on_status_complete, HTTP_CB...)
    global const on_header_field_cb = cfunction(on_header_field, HTTP_DATA_CB...)
    global const on_header_value_cb = cfunction(on_header_value, HTTP_DATA_CB...)
    global const on_headers_complete_cb = cfunction(on_headers_complete, HTTP_CB...)
    global const on_body_cb = cfunction(on_body, HTTP_DATA_CB...)
    global const on_message_complete_cb = cfunction(on_message_complete, HTTP_CB...)
end


text:
#############
# addresses #
#############
#
# Addresses are linked lists of non-Pair values, using Pair as cons.
# Symbols indicate statically known address fields. (TODO unify this)
#
# Therefore, Addresses can be conveniently constructed using Julia's => syntax
# for constructing pairs:
#
#   :a => 2 => :d => "asfd"
#
# which gives:
#
#  Pair(:a, Pair(2, Pair(:d, "asdf")))

###################
# address schemas #
###################

abstract type AddressSchema end

struct StaticAddressSchema <: AddressSchema
    # TODO can add type information
    leaf_nodes::Set{Symbol}
    internal_nodes::Set{Symbol}
end

leaf_node_keys(schema::StaticAddressSchema) = schema.leaf_nodes
internal_node_keys(schema::StaticAddressSchema) = schema.internal_nodes

struct VectorAddressSchema <: AddressSchema end
struct SingleDynamicKeyAddressSchema <: AddressSchema end
struct DynamicAddressSchema <: AddressSchema end
struct EmptyAddressSchema <: AddressSchema end

export AddressSchema
export StaticAddressSchema
export VectorAddressSchema
export SingleDynamicKeyAddressSchema
export DynamicAddressSchema
export EmptyAddressSchema

########################
# abstract address set #
########################

"""
    get_address_schema(T)

    Base.isempty(set)

Is the set empty?

Return a shallow, compile-time address schema.

    has_internal_node(set, addr)

    get_internal_node(set, addr)

Return the set of address set nodes in the given namespace (they are all nonempty?)

    get_internal_nodes(set)

Return an iterator over top-level internal nodes (pairs of keys and nodes)

    has_leaf_node(set, addr)

Is the address in the set

    get_leaf_nodes(set)

Return an iterator over top-level leaf nodes (keys)

    Base.in(set, addr)

Alias for has_leaf_node

    Base.getindex(set, addr)

Alias for get_internal_node

    Base.haskey(set, addr)

Alias for has_internal_node

"""
abstract type AddressSet end

has_internal_node(::AddressSet, addr) = false
get_internal_node(::AddressSet, addr) = throw(KeyError(addr))
has_leaf_node(::AddressSet, addr) = false
Base.in(addr, set::AddressSet) = has_leaf_node(set, addr)
Base.getindex(set::AddressSet, addr) = get_internal_node(set, addr)
Base.haskey(set::AddressSet, addr) = has_internal_node(set, addr)

export AddressSet

# has_internal_node()
#       - used in backprop_trace() OK
#       - used in regenerate() OK
# get_internal_node() [[ alias getindex, KeyError if it doesn't exist ]]
#       - used in backprop_trace() (check if there is one first, and return EmptySet o/w) OK
#       - used in regenerate() (check if there is one first, and return EmptySet o/w) OK
# has_leaf_node() [[ alias in ]]
#       - used in backprop_trace() OK
#       - used in regenerate() OK
# get_leaf_nodes() [[ return iterator over keys in this namespace ]]
# get_internal_nodes() [[ return iterator of pairs of keys and sub-address-set ]]

# for static address sets (with StaticAddressSchema schema)
# static_get_internal_node(set::StaticAddressSet, ::Val{A}) where {A}
#       - used in basic block backprop_trace, regenerate, etc.

# construction interface (for DynamicAddressSet)
# there is no construction interface for StaticAddressSet
# add_leaf_node!(set, addr) [[ alias push ]]
# set_internal_node!(set, addr, node::AddressSet) [[ alias setindex ]]

# there should by a copy constructor for StaticAddressSet from dynamic..
# later selection functions can generate static address sets directly..

#####################
# empty address set #
#####################

struct EmptyAddressSet <: AddressSet end
get_address_schema(::Type{EmptyAddressSet}) = EmptyAddressSchema()
Base.isempty(::EmptyAddressSet) = true

export EmptyAddressSet

######################
# static address set #
######################

# R is a tuple of symbols..
# T is a tuple of symbols
# U the tuple type of internal nodes
struct StaticAddressSet{R,T,U} <: AddressSet
    internal_nodes::NamedTuple{T,U}
end

function Base.isempty(set::StaticAddressSet{R,T,U}) where {R,T,U}
    length(R) == 0 && all(isempty(node) for node in set.internal_nodes)
end

function get_address_schema(::Type{StaticAddressSet{R,T,U}}) where {R,T,U}
    leaf_keys = Set{Symbol}()
    internal_keys = Set{Symbol}()
    for key in R
        push!(leaf_keys, key)
    end
    for (key, _) in zip(T, U.parameters)
        push!(internal_keys, key)
    end
    StaticAddressSchema(leaf_keys, internal_keys)
end

get_leaf_nodes(::StaticAddressSet{R,T,U}) where {R,T,U} = R

get_internal_nodes(set::StaticAddressSet) = pairs(set.internal_nodes)

function has_internal_node(set::StaticAddressSet, key::Symbol)
    haskey(set.internal_nodes, key)
end

function has_internal_node(set::StaticAddressSet, addr::Pair)
    (first, rest) = addr
    if haskey(set.internal_nodes, first)
        node = set.internal_nodes[first]
        has_internal_node(node, rest)
    else
        false
    end
end

function static_get_internal_node(set::StaticAddressSet, ::Val{A}) where {A}
    set.internal_nodes[A]
end

function get_internal_node(set::StaticAddressSet, key::Symbol)
    set.internal_nodes[key]
end

function get_internal_node(set::StaticAddressSet, addr::Pair)
    (first, rest) = addr
    node = set.internal_nodes[first]
    get_internal_node(node, rest)
end

function has_leaf_node(set::StaticAddressSet{R,T,U}, key::Symbol) where {R,T,U}
    key in R # TODO this is O(N)
end

function has_leaf_node(set::StaticAddressSet, addr::Pair)
    (first, rest) = addr
    if haskey(set.internal_nodes, first)
        node = set.internal_nodes[first]
        has_leaf_node(node, rest)
    else
        false
    end
end

function StaticAddressSet(other::AddressSet)
    leaf_keys = (get_leaf_nodes(other)...,)
    internal_keys_and_nodes = collect(get_internal_nodes(other))
    if length(internal_keys_and_nodes) > 0
        (internal_keys, internal_nodes) = collect(zip(internal_keys_and_nodes...))
    else
        (internal_keys, internal_nodes) = ((), ())
    end
    internal_types = map(typeof, internal_nodes)
    StaticAddressSet{leaf_keys,internal_keys,Tuple{internal_types...}}(
        NamedTuple{internal_keys}(internal_nodes))
end

export StaticAddressSet

#######################
# dynamic address set #
#######################

struct DynamicAddressSet <: AddressSet
    leaf_nodes::Set{Any} # set of keys
    internal_nodes::Dict{Any,AddressSet}
end

# invariant: all internal nodes are nonempty
function Base.isempty(set::DynamicAddressSet)
    isempty(set.leaf_nodes) && isempty(set.internal_nodes)
end

DynamicAddressSet() = DynamicAddressSet(Set{Any}(), Dict{Any,DynamicAddressSet}())

get_address_schema(::Type{DynamicAddressSet}) = DynamicAddressSchema()

has_leaf_node(set::DynamicAddressSet, addr) = (addr in set.leaf_nodes)

function has_leaf_node(set::DynamicAddressSet, addr::Pair)
    (first, rest) = addr
    if haskey(set.internal_nodes, first)
        internal_node = set.internal_nodes[first]
        has_leaf_node(internal_node, rest)
    else
        false
    end
end

function has_internal_node(set::DynamicAddressSet, addr)
    haskey(set.internal_nodes, addr)
end

function has_internal_node(set::DynamicAddressSet, addr::Pair)
    (first, rest) = addr
    if haskey(set.internal_nodes, first)
        node = set.internal_nodes[first]
        has_internal_node(node, rest)
    else
        false
    end
end

function get_internal_node(set::DynamicAddressSet, addr)
    set.internal_nodes[addr]
end

function get_internal_node(set::DynamicAddressSet, addr::Pair)
    (first, rest) = addr
    node = set.internal_nodes[first]
    get_internal_node(node, rest)
end

function push_leaf_node!(set::DynamicAddressSet, addr)
    if haskey(set.internal_nodes, addr)
        error("Tried to push_leaf_node! $addr but there is already a namespace rooted at $addr")
    end
    push!(set.leaf_nodes, addr)
end

function push_leaf_node!(set::DynamicAddressSet, addr::Pair)
    (first, rest) = addr
    if haskey(set.internal_nodes, first)
        node = set.internal_nodes[first]
    else
        node = DynamicAddressSet()
        set.internal_nodes[first] = node
    end
    push_leaf_node!(node, rest)
end

function set_internal_node!(set::DynamicAddressSet, addr, node)
    if isempty(node)
        return
    end
    if addr in set.leaf_nodes
        error("Tried to set_internal_node! $addr but that addres is already a leaf")
    else
        set.internal_nodes[addr] = node
    end
end

function set_internal_node!(set::DynamicAddressSet, addr::Pair, node)
    if isempty(node)
        return
    end
    (first, rest) = addr
    node = set.internal_nodes[first]
    set_internal_node!(node, rest, node)
end

get_leaf_nodes(addrs::DynamicAddressSet) = addrs.leaf_nodes

get_internal_nodes(addrs::DynamicAddressSet) = addrs.internal_nodes

export DynamicAddressSet
export push_leaf_node!, set_internal_node!

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
using MeshCat, MeshCatMechanisms, RigidBodyDynamics
using FileIO, MeshIO, GeometryBasics, CoordinateTransformations, Rotations, Meshing
using Plots

include(joinpath(pwd(),"models/kuka/kuka_utils.jl"))

urdf_path = joinpath(pwd(), "models/kuka/temp/kuka.urdf")

kuka = MeshCatMechanisms.parse_urdf(urdf_path, remove_fixed_tree_joints = true)
kuka_visuals = MeshCatMechanisms.URDFVisuals(urdf_path)

state = MechanismState(kuka)
state_cache = StateCache(kuka)
result = DynamicsResult(kuka)
result_cache = DynamicsResultCache(kuka)
nq = num_positions(kuka)
ee = findbody(kuka, "iiwa_link_7")
ee_point = Point3D(default_frame(ee), 0.0, 0.0, 0.0)
ee_jacobian_frame = ee_point.frame
ee_jacobian_path = path(kuka, root_body(kuka), ee)
world = root_frame(kuka)

# Kuka iiwa arm parsed from URDF using RigidBodyDynamics.jl
# + 3D particle
struct KukaParticle
	dim::Dimensions

	# particle
	mp

    state_cache1
    state_cache2
    state_cache3

    result_cache1
    result_cache2
    result_cache3

    world

	ee
    ee_point
    ee_jacobian_frame
    ee_jacobian_path
end

results_cache1 = DynamicsResultCache(kuka)
results_cache2 = DynamicsResultCache(kuka)
results_cache3 = DynamicsResultCache(kuka)

state_cache1 = StateCache(kuka)
state_cache2 = StateCache(kuka)
state_cache3 = StateCache(kuka)

function kuka_q(q)
	view(q, 1:7)
end

function particle_q(q)
	view(q, 8:10)
end

# Methods
function M_func(m::KukaParticle, q::AbstractVector{T}) where T
    # return Diagonal(m.mp * ones(3))
    state = m.state_cache3[T]
    result = m.result_cache3[T]
    set_configuration!(state, kuka_q(q))
    mass_matrix!(result.massmatrix, state)
    return SMatrix{10, 10}(cat(result.massmatrix, Diagonal(m.mp * ones(3)), dims = (1, 2)))
end

function C_func(m::KukaParticle, q, q̇)
    # return SVector{3}([0.0, 0.0, m.mp * 9.81])
    if eltype(q) <: ForwardDiff.Dual
		return SVector{10}([dynamics_bias_q(m, kuka_q(q), kuka_q(q̇))..., 0.0, 0.0, m.mp * 9.81])
	elseif eltype(q̇) <: ForwardDiff.Dual
		return SVector{10}([dynamics_bias_q̇(m, kuka_q(q), kuka_q(q̇))..., 0.0 ,0.0, m.mp * 9.81])
	else
		return SVector{10}([_dynamics_bias(m, kuka_q(q), kuka_q(q̇))..., 0.0, 0.0, m.mp * 9.81])
	end
end

function _dynamics_bias(m::KukaParticle, q, q̇)
	T = eltype(m.mp)
	state = m.state_cache1[T]
	result = m.result_cache1[T]
	set_configuration!(state, q)
	set_velocity!(state, q̇)
	dynamics_bias!(result, state)
    return result.dynamicsbias
end

function dynamics_bias_q(m::KukaParticle, q::AbstractVector{T}, q̇) where T
    state = m.state_cache2[T]
    result = m.result_cache2[T]
    set_configuration!(state, q)
    set_velocity!(state,q̇)
    dynamics_bias!(result, state)
    return result.dynamicsbias
end

function dynamics_bias_q̇(m::KukaParticle, q, q̇::AbstractVector{T}) where T
    state = m.state_cache3[T]
    result = m.result_cache3[T]
    set_configuration!(state, q)
    set_velocity!(state, q̇)
    dynamics_bias!(result, state)
    return result.dynamicsbias
end

function B_func(m::KukaParticle, q)
    SMatrix{10, 7}([Diagonal(ones(7)); zeros(3, 7)])
end

function kinematics_ee(m::KukaParticle, q::AbstractVector{T}) where T
    state = m.state_cache3[T]
    set_configuration!(state, kuka_q(q))
    return transform(state, m.ee_point, m.world).v
end

l_string = 0.25

function ϕ_func(m::KukaParticle, q::AbstractVector{T}) where T
	p_ee = kinematics_ee(m, q)
	p_p = q[8:10]
	diff = (p_ee - p_p)
	d_ee_p = sqrt(diff' * diff)
    SVector{1}([l_string - d_ee_p])
end

function P_func(m::KukaParticle, q::AbstractVector{T}) where T
    state = m.state_cache3[T]

    set_configuration!(state, kuka_q(q))

    pj1 = PointJacobian(transform(state, m.ee_point, m.world).frame,
		zeros(T, 3, 7))

    ee_in_world = transform(state, m.ee_point, m.world)

    point_jacobian!(pj1, state, m.ee_jacobian_path,ee_in_world) #TODO confirm this is correct

	p_ee = ee_in_world.v
	p_p = q[8:10]
	diff = (p_ee - p_p)

    return Array(-1.0 * 2.0 * diff' * [pj1.J[1:3,:] -1.0 * Diagonal(ones(3))])
end

# dimensions
nq = 7 + 3
nu = 7
nw = 0
nc = 1

# particle mass
mp = 0.01

model = KukaParticle(
	Dimensions(nq, nu, nw, nc),
	mp,
    state_cache1, state_cache2, state_cache3,
    results_cache1, results_cache2, results_cache3,
    world,
	ee,
	ee_point,
	ee_jacobian_frame,
	ee_jacobian_path)

function lagrangian_derivatives(model::KukaParticle, q, v)
	D1L = -1.0 * C_func(model, q, v)
    D2L = M_func(model, q) * v
	return D1L, D2L
end

function gravity_compensation(model::KukaParticle, h, q0, q1)
    q2 = q1
	# evalutate at midpoint
	qm1 = 0.5 * (q0 + q1)
    vm1 = (q1 - q0) / h[1]
    qm2 = 0.5 * (q1 + q2)
    vm2 = (q2 - q1) / h[1]

	D1L1, D2L1 = lagrangian_derivatives(model, qm1, vm1)
	D1L2, D2L2 = lagrangian_derivatives(model, qm2, vm2)

	# return 0.0
	return (0.5 * h[1] * D1L1 + D2L1 + 0.5 * h[1] * D1L2 - D2L2)[1:7]
end

function dynamics(model::KukaParticle, h, q0, q1, u1, λ1, q2)
	# evalutate at midpoint
	qm1 = 0.5 * (q0 + q1)
    vm1 = (q1 - q0) / h[1]
    qm2 = 0.5 * (q1 + q2)
    vm2 = (q2 - q1) / h[1]

	D1L1, D2L1 = lagrangian_derivatives(model, qm1, vm1)
	D1L2, D2L2 = lagrangian_derivatives(model, qm2, vm2)

    u_gravity = gravity_compensation(model, h, q0, q1)

	return (0.5 * h[1] * D1L1 + D2L1 + 0.5 * h[1] * D1L2 - D2L2
		+ B_func(model, qm2) * (u1 - u_gravity)
        + transpose(P_func(model, q2)) * λ1)
end

function residual(model::KukaParticle, z, θ, κ)
    nq = model.dim.q
    nu = model.dim.u
    nc = model.dim.c

	q0 = θ[1:nq]
    q1 = θ[nq .+ (1:nq)]
    u1 = θ[2nq .+ (1:nu)]
    h = θ[2nq + nu .+ (1:1)]

    q2 = z[1:nq]
    λ1 = z[nq .+ (1:nc)]
    s1 = z[nq + nc .+ (1:nc)]

    [
     dynamics(model, h, q0, q1, u1, λ1, q2);
     s1 .- ϕ_func(model, q2);
     λ1 .* s1 .- κ;
    ]
end

nz = nq + nc + nc
nθ = nq + nq + nu + 1

function r_func(r, z, θ, κ)
    r .= residual(model, z, θ, κ)
end

function rz_func(rz, z, θ)
    r(a) = residual(model, a, θ, 0.0)
    rz .= ForwardDiff.jacobian(r, z)
end

function rθ_func(rθ, z, θ)
    r(a) = residual(model, z, a, 0.0)
    rθ .= ForwardDiff.jacobian(r, θ)
end

idx_ineq = collect([nq .+ (1:(nc + nc))]...)
z_subset_init = 1.0 * ones(nz - nq)
rz_array = zeros(nz, nz)
rθ_array = zeros(nz, nθ)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
include("runtests-prerequisite.jl")

module Test_util_prettify_id_values_as_str

    using PostgresORM

    mutable struct MyMovie <: PostgresORM.IEntity

      title::Union{Missing,String}
      year::Union{Missing,Int32}
      description::Union{Missing,String}

      MyMovie(args::NamedTuple) = MyMovie(;args...)
      MyMovie(;title = missing,
               year = missing,
               description = missing) = (
                x = new(missing, missing);
                x.title = title;
                x.year = year;
                x.description = description;
                return x
                )
    end

    mutable struct MyNewStruct <: PostgresORM.IEntity

      actor_id::Union{Missing,Int32}
      first_name::Union{Missing,String}
      favorite_movie::Union{Missing, MyMovie}

      MyNewStruct(args::NamedTuple) = MyNewStruct(;args...)
      MyNewStruct(;
            actor_id = missing,
            first_name = missing,
            favorite_movie = missing,
     ) = (
                x = new(missing, missing, missing);
                x.actor_id = actor_id;
                x.first_name = first_name;
                x.favorite_movie = favorite_movie;
                return x
                )

    end
end # module

module MyNewStructORM
    data_type = Main.Test_util_prettify_id_values_as_str.MyNewStruct
    Main.PostgresORM.get_orm(x::Main.Test_util_prettify_id_values_as_str.MyNewStruct) = return(MyNewStructORM)
    get_table_name() = "public.my_new_struct"
    const columns_selection_and_mapping =
      Dict(
         :actor_id => "actor_id",
         :first_name => "first_name",
         :favorite_movie => ["favorite_movie_id"]
      )
    get_id_props() = return [:actor_id,:first_name,:favorite_movie]

    # A dictionnary of mapping between fields symbols and overriding types
    #   Left hanside is the field symbol ; right hand side is the type override
    const types_override = Dict()

    const track_changes = false
end

module MyMovieORM
    data_type = Main.Test_util_prettify_id_values_as_str.MyMovie
    Main.PostgresORM.get_orm(x::Main.Test_util_prettify_id_values_as_str.MyMovie) = return(MyMovieORM)
    get_table_name() = "public.my_movie"
    const columns_selection_and_mapping =
      Dict(
         :title => "title",
         :year => "year"
      )
    get_id_props() = return [:title,:year]

    # A dictionnary of mapping between fields symbols and overriding types
    #   Left hanside is the field symbol ; right hand side is the type override
    const types_override = Dict()

    const track_changes = false
end

@testset "Test `util_prettify_id_values_as_str(values::Vector{<:Any})`" begin

    mymovie =
        Test_util_prettify_id_values_as_str.MyMovie(title = "movie1",
                                                         year = 2007,
                                                         description = "description1")

    object = Test_util_prettify_id_values_as_str.MyNewStruct(actor_id = 2000,
                                                             first_name = "bob",
                                                             favorite_movie = mymovie)

    PostgresORM.Controller.util_prettify_id_values_as_str(object)

end


text:
#unit step function
function u(t::Float64)::Float64 #unit step function
    return ifelse( t < 0, zero(t), one(t) )
end

#unnormalized Gaussian function
function 𝒩ᵤ(x::Float64; μ::Float64, σ::Float64)::Float64
    return exp(-0.5 * ((x - μ) / σ)^2)
end

#normalized Gaussian function
function 𝒩(x::Float64; μ::Float64, σ::Float64)::Float64
    return 1 / (σ * sqrt(2 * π)) * exp(-0.5 * ((x - μ) / σ)^2)
end

#Dirac delta (Gaussian approximation) amplitude normalized
function δn(t::Float64,σ::Float64=1.0e-11)::Float64
   return 𝒩ᵤ(t, μ=0.0, σ=σ)
end

#Dirac delta (Gaussian approximation) unit area
function δ(t::Float64,σ::Float64=1.0e-10)::Float64
   return 𝒩(t, μ=0.0, σ=σ)
end


text:
"""
$(TYPEDEF)

The `Planar` joint type allows translation along two orthogonal vectors, referred to as ``x`` and ``y``,
as well as rotation about an axis ``z = x \\times y``.

The components of the 3-dimensional configuration vector ``q`` associated with a `Planar` joint are the
``x``- and ``y``-coordinates of the translation, and the angle of rotation ``\\theta`` about ``z``, in that order.

The components of the 3-dimension velocity vector ``v`` associated with a `Planar` joint are the
``x``- and ``y``-coordinates of the linear part of the joint twist, expressed in the frame after the joint,
followed by the ``z``-component of the angular part of this joint twist.

!!! warning

    For the `Planar` joint type, ``v \\neq \\dot{q}``! Although the angular parts of ``v`` and ``\\dot{q}``
    are the same, their linear parts differ. The linear part of ``v`` is the linear part of ``\\dot{q}``, rotated
    to the frame after the joint. This parameterization was chosen to allow the translational component of the
    joint transform to be independent of the rotation angle ``\\theta`` (i.e., the rotation is applied **after** the
    translation), while still retaining a constant motion subspace expressed in the frame after the joint.

"""
struct Planar{T} <: JointType{T}
    x_axis::SVector{3, T}
    y_axis::SVector{3, T}
    rot_axis::SVector{3, T}

    @doc """
    $(SIGNATURES)

    Construct a new `Planar` joint type with the ``xy``-plane in which translation is allowed defined
    by 3-vectors `x` and `y` expressed in the frame before the joint.
    """ ->
    function Planar{T}(x_axis::AbstractVector, y_axis::AbstractVector) where {T}
        x, y = map(axis -> normalize(SVector{3}(axis)), (x_axis, y_axis))
        @assert isapprox(x ⋅ y, 0; atol = 100 * eps(T))
        new{T}(x, y, x × y)
    end
end

Planar(x_axis::AbstractVector, y_axis::AbstractVector) = Planar{promote_eltype(x_axis, y_axis)}(x_axis, y_axis)

Base.show(io::IO, jt::Planar) = print(io, "Planar joint with x-axis $(jt.x_axis) and y-axis $(jt.y_axis)")

function Random.rand(::Type{Planar{T}}) where {T}
    x = normalize(randn(SVector{3, T}))
    y = normalize(randn(SVector{3, T}))
    y = normalize(y - (x ⋅ y) * x)
    Planar(x, y)
end

num_positions(::Type{<:Planar}) = 3
num_velocities(::Type{<:Planar}) = 3
has_fixed_subspaces(jt::Planar) = true
isfloating(::Type{<:Planar}) = false

@propagate_inbounds function rand_configuration!(q::AbstractVector, ::Planar)
    T = eltype(q)
    q[1] = rand() - T(0.5)
    q[2] = rand() - T(0.5)
    q[3] = randn()
    nothing
end

@propagate_inbounds function joint_transform(jt::Planar, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector)
    rot = RotMatrix(AngleAxis(q[3], jt.rot_axis[1], jt.rot_axis[2], jt.rot_axis[3], false))
    trans = jt.x_axis * q[1] + jt.y_axis * q[2]
    Transform3D(frame_after, frame_before, rot, trans)
end

@propagate_inbounds function joint_twist(jt::Planar, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector, v::AbstractVector)
    angular = jt.rot_axis * v[3]
    linear = jt.x_axis * v[1] + jt.y_axis * v[2]
    Twist(frame_after, frame_before, frame_after, angular, linear)
end

@propagate_inbounds function joint_spatial_acceleration(jt::Planar, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector, v::AbstractVector, vd::AbstractVector)
    S = promote_eltype(jt, q, v, vd)
    angular = jt.rot_axis * vd[3]
    linear = jt.x_axis * vd[1] + jt.y_axis * vd[2]
    SpatialAcceleration{S}(frame_after, frame_before, frame_after, angular, linear)
end

@inline function motion_subspace(jt::Planar, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector)
    S = promote_eltype(jt, q)
    angular = hcat(zero(SMatrix{3, 2, S}), jt.rot_axis)
    linear = hcat(jt.x_axis, jt.y_axis, zero(SVector{3, S}))
    GeometricJacobian(frame_after, frame_before, frame_after, angular, linear)
end

@inline function constraint_wrench_subspace(jt::Planar, joint_transform::Transform3D)
    S = promote_eltype(jt, joint_transform)
    angular = hcat(zero(SVector{3, S}), jt.x_axis, jt.y_axis)
    linear = hcat(jt.rot_axis, zero(SMatrix{3, 2, S}))
    WrenchMatrix(joint_transform.from, angular, linear)
end

@inline function bias_acceleration(jt::Planar, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector, v::AbstractVector)
    S = promote_eltype(jt, q, v)
    zero(SpatialAcceleration{S}, frame_after, frame_before, frame_after)
end

@propagate_inbounds function joint_torque!(τ::AbstractVector, jt::Planar, q::AbstractVector, joint_wrench::Wrench)
    τ[1] = dot(linear(joint_wrench), jt.x_axis)
    τ[2] = dot(linear(joint_wrench), jt.y_axis)
    τ[3] = dot(angular(joint_wrench), jt.rot_axis)
    nothing
end

@propagate_inbounds function configuration_derivative_to_velocity!(v::AbstractVector, jt::Planar, q::AbstractVector, q̇::AbstractVector)
    vlinear = RotMatrix(-q[3]) * SVector(q̇[1], q̇[2])
    v[1] = vlinear[1]
    v[2] = vlinear[2]
    v[3] = q̇[3]
    nothing
end

@propagate_inbounds function velocity_to_configuration_derivative!(q̇::AbstractVector, jt::Planar, q::AbstractVector, v::AbstractVector)
    q̇linear = RotMatrix(q[3]) * SVector(v[1], v[2])
    q̇[1] = q̇linear[1]
    q̇[2] = q̇linear[2]
    q̇[3] = v[3]
    nothing
end

@propagate_inbounds function configuration_derivative_to_velocity_adjoint!(out, jt::Planar, q::AbstractVector, f)
    outlinear = RotMatrix(q[3]) * SVector(f[1], f[2])
    out[1] = outlinear[1]
    out[2] = outlinear[2]
    out[3] = f[3]
    nothing
end

@propagate_inbounds function velocity_to_configuration_derivative_jacobian(::Planar, q::AbstractVector)
    # TODO: use SMatrix(RotZ(q[3]) once it's as fast
    rot = RotMatrix(q[3])
    @inbounds return @SMatrix([rot[1] rot[3] 0;
                               rot[2] rot[4] 0;
                               0         0   1])
end

@propagate_inbounds function configuration_derivative_to_velocity_jacobian(::Planar, q::AbstractVector)
    # TODO: use SMatrix(RotZ(-q[3]) once it's as fast
    rot = RotMatrix(-q[3])
    @inbounds return @SMatrix([rot[1] rot[3] 0;
                               rot[2] rot[4] 0;
                               0         0   1])
end


text:
@testset "types" begin
    @testset "Mag" begin
        @test Mag() isa Mag
        @test Mag(Mag()) isa Mag
        @test Mag(Arf()) isa Mag
        x = Mag()
        y = copy(x)
        Arblib.set!(x, 1)
        @test x == Mag(1)
        @test y == Mag(0)
    end

    @testset "Arf" begin
        @test Arf() isa Arf
        @test precision(Arf(prec = 80)) == 80
        @test Arf(UInt64(1)) isa Arf
        @test precision(Arf(UInt64(1), prec = 80)) == 80
        @test Arf(1) isa Arf
        @test precision(Arf(1, prec = 80)) == 80
        x = Arf()
        y = copy(x)
        Arblib.set!(x, 1)
        @test x == Arf(1)
        @test y == Arf(0)
    end

    @testset "Arb" begin
        @test Arb() isa Arb
        @test precision(Arb(prec = 80)) == 80
        x = Arb()
        y = copy(x)
        Arblib.set!(x, 1)
        @test x == Arb(1)
        @test y == Arb(0)
    end

    @testset "Acb" begin
        @test Acb() isa Acb
        @test precision(Acb(prec = 80)) == 80
        x = Acb()
        y = copy(x)
        Arblib.set!(x, 1)
        @test x == Acb(1)
        @test y == Acb(0)
    end

    @testset "ArbPoly" begin
        @test ArbPoly() isa ArbPoly
        @test precision(ArbPoly(prec = 80)) == 80
    end

    @testset "ArbSeries" begin
        @test ArbSeries() isa ArbSeries
        @test precision(ArbSeries(prec = 80)) == 80
        @test ArbSeries(degree = 3).degree == 3
    end

    @testset "AcbPoly" begin
        @test AcbPoly() isa AcbPoly
        @test precision(AcbPoly(prec = 80)) == 80
    end

    @testset "AcbSeries" begin
        @test AcbSeries() isa AcbSeries
        @test precision(AcbSeries(prec = 80)) == 80
        @test AcbSeries(degree = 3).degree == 3
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ d5a6bdfa-fbcd-436e-99cd-eadb8f5c63c8
using GMT

# ╔═╡ 89b3b5b0-d071-11eb-3994-7932442f1611
md"
# Ternary plots

To plot points on a ternary diagram at the positions listed in the file ternary.txt (that GMT knows where to find it), with default annotations and gridline spacings, using the specified labeling, do
"

# ╔═╡ 9ffa3892-5e30-4e46-a00a-a1682eebd888
ternary("@ternary.txt", labels=("Clay","Silt","Sand"), marker=:p, show=true)

# ╔═╡ be06e94b-ec5c-442b-9273-d55b6a6fede8
md"
-- Ok, good but a bit dull. What about coloring the points? And if I want to have the axes runing in clock-wise order? And what about adding a percentage symbol to the annotations?

Simple. First we create a colormap and to rotte the axes we use the option `clockwise=true`. Regarding the `%` sign, it requires using the `frame` option and that obliges to be explicit on the axes labels because we are no longer using handy defaults.
"

# ╔═╡ d9a81565-1dc7-482f-ac74-5102b5c79969
# Make use of the knowledge that z ranges berween 0 and 71 (gmtinfo module is a friend)
C = makecpt(T=(0,71));

# ╔═╡ 20979dbf-3d32-4764-975b-286469458b92
ternary("@ternary.txt", frame=(annot=:auto, grid=:a, ticks=:a, alabel="Clay", blabel="Silt", clabel="Sand", suffix=" %"), marker=:p, cmap=C, clockwise=true, show=true)

# ╔═╡ 6c7960c0-d0a4-11eb-3147-f92e4a7436e6
md"
-- Ah, much better, but now I would like to display the above data as an image.

Solution: use the `image=true` option. Note that we may skip the `cmap` option and an automatic `cmap` is compute for us (but one can use whatever cmap we want, just create a colormap with the wished colors)
"

# ╔═╡ ceff64b0-d0a4-11eb-1881-7d9a7b2be78e
ternary("@ternary.txt", frame=(annot=:auto, grid=:a, ticks=:a, alabel="Clay", blabel="Silt", clabel="Sand", suffix=" %"), marker=:p, image=true, clockwise=true, show=true)

# ╔═╡ 4700fb42-d0a5-11eb-3914-53b959143ec5
md"
-- And to overlay some contours?

Add the `contour` option. This option works either with automatically picked parameters or under the user full control (which contours to draw and which to annotate, etc). For simplicity we could use the automatic mode (just set `contour=true`) but the ternary plots may have several short contour lines that would not be annotated because they are too short for the default setting. So, and for demonstration sake, we will use the explicit `contour` form where we set also the distance between the labels.
"

# ╔═╡ c636e63e-d0a5-11eb-2d6a-ab44a74366a1
ternary("@ternary.txt", frame=(annot=:auto, grid=:a, ticks=:a, alabel="Clay", blabel="Silt", clabel="Sand", suffix=" %"), clockwise=true, image=true, contour=(annot=10,cont=5,labels=(distance=3,)), colorbar=true, show=true)

# ╔═╡ f9c30160-d0a5-11eb-18ea-81a47cb85b76
md"
And we can do a `contourf` style plot too, but in this case only the area inside the data cloud is imaged since the method used involves a Delaunay triangulation.
"

# ╔═╡ 161efc10-d0a6-11eb-3d7f-af1d658a3cf0
ternary("@ternary.txt", frame=(annot=:auto, grid=:a, ticks=:a, alabel="Clay", blabel="Silt", clabel="Sand", suffix=" %"), marker=:p, clockwise=true, contourf=(annot=10, cont=5), show=true)

# ╔═╡ Cell order:
# ╟─89b3b5b0-d071-11eb-3994-7932442f1611
# ╠═d5a6bdfa-fbcd-436e-99cd-eadb8f5c63c8
# ╠═9ffa3892-5e30-4e46-a00a-a1682eebd888
# ╟─be06e94b-ec5c-442b-9273-d55b6a6fede8
# ╠═d9a81565-1dc7-482f-ac74-5102b5c79969
# ╠═20979dbf-3d32-4764-975b-286469458b92
# ╟─6c7960c0-d0a4-11eb-3147-f92e4a7436e6
# ╠═ceff64b0-d0a4-11eb-1881-7d9a7b2be78e
# ╟─4700fb42-d0a5-11eb-3914-53b959143ec5
# ╠═c636e63e-d0a5-11eb-2d6a-ab44a74366a1
# ╟─f9c30160-d0a5-11eb-18ea-81a47cb85b76
# ╠═161efc10-d0a6-11eb-3d7f-af1d658a3cf0


text:
using LinearAlgebra

import Libdl
const depsfile = joinpath(@__DIR__, "deps.jl")

# If BLAS was compiled with MKL and the user wants MKL-based FFTs, we'll oblige.
# In that case, we have to do this little dance to get around having to use BinDeps
# for a library that's already linked to Julia.
settings = joinpath(first(DEPOT_PATH), "prefs", "FFTW")
mkpath(dirname(settings))
if haskey(ENV, "JULIA_FFTW_PROVIDER")
    provider = ENV["JULIA_FFTW_PROVIDER"]
    open(f -> println(f, provider), settings, "w")
elseif isfile(settings)
    provider = readchomp(settings)
else
    provider = "FFTW"
    open(f -> println(f, provider), settings, "w")
end
if provider == "MKL"
    if BLAS.vendor() === :mkl
        mklpath = Libdl.dlpath("libmkl_rt")
    else
        using Conda
        Conda.add("mkl_fft")
        mklpath = joinpath(Conda.lib_dir(Conda.ROOTENV), "libmkl_rt")
    end
    mklpath = escape_string(mklpath)
    isfile(depsfile) && rm(depsfile, force=true)
    open(depsfile, "w") do f
        println(f, """
            # This is an auto-generated file, do not edit
            import Libdl
            const libfftw3 = "$mklpath"
            const libfftw3f = libfftw3
            function check_deps()
                if Libdl.dlopen_e(libfftw3) == C_NULL
                    error("Unable to load MKL from '$mklpath'.\\n",
                          "Please rerun Pkg.build(\\"FFTW\\") and restart Julia.")
                end
            end
        """)
    end
elseif provider != "FFTW"
    error("Unrecognized JULIA_FFTW_PROVIDER \"$provider\".\n",
          "To fix this, set ENV[\"JULIA_FFTW_PROVIDER\"] to \"FFTW\" or \"MKL\"\n",
          "and rerun Pkg.build(\"FFTW\").")
else
    include("build_fftw.jl")
end


text:
using MMA
using Test

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


text:
export ft_writeee

function ft_writeee(ft_handle::Culong, wordoffset::Integer, value::UInt16)
  @assert wordoffset >= 0
  ft_status = ccall((:FT_WriteEE, d2xx),
                     Cuint,
                     (Culong, Cuint, Cushort),
                     ft_handle, wordoffset, value)
  checkstatus(ft_status)
  return nothing
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/HeatTransfer.jl/blob/master/LICENSE

"""
    HeatTransfer

Heat transfer analysis for JuliaFEM.

# Problem types

- `PlaneHeat` for 2d problems
- `Heat` for 3d problems

# Fields used in formulation

- `thermal conductivity`
- `heat source`
- `heat flux`
- `external temperature`
- `heat transfer coefficient`

# References

- https://en.wikipedia.org/wiki/Heat_equation
- https://en.wikipedia.org/wiki/Heat_capacity
- https://en.wikipedia.org/wiki/Heat_flux
- https://en.wikipedia.org/wiki/Thermal_conduction
- https://en.wikipedia.org/wiki/Thermal_conductivity
- https://en.wikipedia.org/wiki/Thermal_diffusivity
- https://en.wikipedia.org/wiki/Volumetric_heat_capacity

"""
module HeatTransfer

using FEMBase, LinearAlgebra, SparseArrays

struct PlaneHeat <: FieldProblem end
struct Heat <: FieldProblem end

FEMBase.get_unknown_field_name(::PlaneHeat) = "temperature"
FEMBase.get_unknown_field_name(::Heat) = "temperature"

function FEMBase.assemble_elements!(problem::Problem{P}, assembly::Assembly,
                                    elements::Vector{Element{B}}, time::Float64) where
                                    {B,P<:Union{PlaneHeat,Heat}}

    bi = BasisInfo(B)
    ndofs = length(bi)
    Ke = zeros(ndofs, ndofs)
    fe = zeros(ndofs)

    for element in elements
        fill!(Ke, 0.0)
        fill!(fe, 0.0)
        for ip in get_integration_points(element)
            J, detJ, N, dN = element_info!(bi, element, ip, time)
            s = ip.weight * detJ
            k = element("thermal conductivity", ip, time)
            Ke += s * k*dN'*dN
            if haskey(element, "heat source")
                f = element("heat source", ip, time)
                fe += s * N'*f
            end
        end
        if haskey(element, "temperature")
            T = [element("temperature", time)...]
            fe -= Ke*T
        end
        gdofs = get_gdofs(problem, element)
        add!(assembly.K, gdofs, gdofs, Ke)
        add!(assembly.f, gdofs, fe)
    end

end

function FEMBase.assemble_elements!(problem::Problem{PlaneHeat}, assembly::Assembly,
                                    elements::Vector{Element{B}}, time::Float64) where
                                    {B<:Union{Seg2,Seg3}}
    return assemble_boundary_elements!(problem, assembly, elements, time)
end

function FEMBase.assemble_elements!(problem::Problem{Heat}, assembly::Assembly,
                                    elements::Vector{Element{B}}, time::Float64) where
                                    {B<:Union{Tri3,Quad4,Tri6,Quad8,Quad9}}
    return assemble_boundary_elements!(problem, assembly, elements, time)
end

function assemble_boundary_elements!(problem::Problem, assembly::Assembly,
                                     elements::Vector{Element{B}}, time::Float64) where B

    bi = BasisInfo(B)
    ndofs = length(bi)
    Ke = zeros(ndofs, ndofs)
    fe = zeros(ndofs)

    for element in elements
        fill!(fe, 0.0)
        fill!(Ke, 0.0)
        for ip in get_integration_points(element, 2)
            J, detJ, N, dN = element_info!(bi, element, ip, time)
            s = ip.weight * detJ
            if haskey(element, "heat flux")
                g = element("heat flux", ip, time)
                fe += s * N'*g
            end
            if haskey(element, "heat transfer coefficient") && haskey(element, "external temperature")
                h = element("heat transfer coefficient", ip, time)
                Tu = element("external temperature", ip, time)
                Ke += s * h*N'*N
                fe += s * N'*h*Tu
            end
        end
        if haskey(element, "temperature")
            T = [element("temperature", time)...]
            fe -= Ke*T
        end
        gdofs = get_gdofs(problem, element)
        add!(assembly.K, gdofs, gdofs, Ke)
        add!(assembly.f, gdofs, fe)
    end

end

export Heat, PlaneHeat

end


text:
using Test, MusicAlbums

const DBFILE = joinpath(dirname(pathof(MusicAlbums)), "../test/albums2.sqlite")
const AUTHFILE = "file://" * joinpath(dirname(pathof(MusicAlbums)), "../resources/authkeys.json")

server = @async MusicAlbums.run(DBFILE, AUTHFILE)

Client.createUser("vanhp", "julia4evah!")
user = Client.loginUser("vanhp", "julia4evah!")

using HTTP; HTTP.CookieRequest.default_cookiejar[1]

alb1 = Client.createAlbum("Free Yourself Up", "Lake Street Dive", 2018, ["Baby Don't Leave Me Alone With My Thoughts", "Good Kisser"])
@test Client.pickAlbumToListen() == alb1
@test Client.pickAlbumToListen() == alb1

@test Client.getAlbum(alb1.id) == alb1

push!(alb1.songs, "Shame, Shame, Shame")
alb2 = Client.updateAlbum(alb1)
@test length(alb2.songs) == 3
@test length(Client.getAlbum(alb1.id).songs) == 3

Client.deleteAlbum(alb1.id)
# Client.pickAlbumToListen()

# alb2 = Client.createAlbum("Haunted Heart", "Charlie Haden Quartet West", 1991, ["Introduction", "Hello My Lovely"])
# @test Client.pickAlbumToListen() == alb2

# Client.createAlbum("Hazards Of Love", "The Decemberists", 2009, ["Prelude", "The Hazards Of Love 1"])
# Client.createAlbum("Tapestry", "Carole King", 1971, ["I Feel The Earth Move", "So Far Away"])

max_stars_repo_name: vanhp/MicroServiceJulia

max_issues_repo_name: vanhp/MicroServiceJulia

max_forks_repo_name: vanhp/MicroServiceJulia


text:
using LinearAlgebra
using CVChannel
using Test

#This script looks at the werner-holevo channel communication value for d=3
#It uses the general linear program to generate the constraints
println("\n Here we will see the Werner-Holevo channels are multiplicative")
println("for p less than 0.3, by searching over p in {0,0.01,...,0.35}.")

@testset "Non-multiplicativity of WH channels for p<0.35" begin
    λ_vals = [0:0.001:0.35;];
    ctr = 1;
    results = zeros(length(λ_vals),5)
    println("\nBeginning to get the data.")
    for λ_id in λ_vals
        if ctr == Int(round(length(λ_vals)-1)/2)
            println("Now halfway through.")
        end
        #Gets single copy cv
        A,B,g,a = generalWHLPConstraints(1,3,λ_id*ones(1))
        cv1, v1 =wernerHolevoCVPPT(1,3,A,B,g,a)
        #Gets parallely copy cv
        A,B,g,a = generalWHLPConstraints(2,3,λ_id*ones(2))
        cv2, v2 =wernerHolevoCVPPT(2,3,A,B,g,a)

        results[ctr,:] = [λ_id, cv1, cv2, cv1^2, cv2-cv1^2]
        ctr = ctr + 1
    end
    printed_30 = false
    print_array = zeros(1,5)
    for i in [1:Int(round(length(λ_vals)/10)):length(λ_vals);]
        if results[i,1]>0.3 && printed_30 == false
            print_array = [print_array; results[301,:]']
            printed_30 = true
        end
        print_array = [print_array; results[i,:]']
    end
    print_array = print_array[2:end,:]
    println("\n Here is a subset of the data. Column labels at bottom.")
    show(stdout, "text/plain", print_array)
    println("\n λ      cv(N)    cv(N^2)   diff")
    println("\nThus by looking at the diff column, we see what we were verifying.")
    info_vec = Vector{Union{Nothing,String}}(nothing, size(results)[1] + 1)
    info_vec[1] = "INFO:"
    info_vec[2] = "Generated by werner-holevo-non-multiplicativity.jl"
    header = ["λ" "cv(N)" "cv(N^2)" "cv(N)^2" "non-mult"]

    println("Please name the file you'd like to write the results to: \n")
    data_to_save = hcat(vcat(header,results),info_vec)
    file_name = readline()
    file_to_open = string(file_name,".csv")
    writedlm(file_to_open, data_to_save, ',')

    @test all(result -> result > 0 || isapprox(result, 0, atol=2e-6), results[:,5])
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
export divAlongDivisor

"""
    divAlongDivisor(::Array{Int64,1},::StackyFan,::Dict)

    Calculates the simplified divisorial index along a divisor of a given cone in a fan with divisorial rays. Takes the subcone consisting of all relevant residual rays along with the specified divisor, and counts the number of residual rays that are relevant with respect to that subcone.

# Examples
```jldoctest AlgDp
julia> F=makeStackyFan([1 0 0; 0 1 0; 1 0 2],[[0,1,2]],[1,2,1])

julia> div=Dict([1,0,0]=>0,[0,1,0]=>0,[1,0,2]=>1)

julia> divAlongDivisor([1,2,3],F,div,[1,0,2])
2
julia> div=Dict([1,0,0]=>0,[0,1,0]=>1,[1,0,2]=>1) #[1 0 0] [1 0 2] - [1 0 0] has non-zero projection

julia> divAlongDivisor([1,2,3],F,div,[1,0,2])
1
julia> div=Dict([1,0,0]=>1,[0,1,0]=>0,[1,0,2]=>1) #[0 1 0] has stacky value >1

julia> divAlongDivisor([1,2,3],F,div,[1,0,2])
1
julia> div=Dict([1,0,0]=>0,[0,1,0]=>1,[1,0,2]=>1) #[1 0 0] [0 1 0] - [1 0 0] is independent

julia> divAlongDivisor([1,2,3],F,div,[0,1,0])
0
```
"""
function divAlongDivisor(cone::Array{Int64,1}, F::StackyFan, div::Dict, divisor::Array{Int64,1})
    # This may not be a very clear name, doesn't mention "simplified" or "index"
    X=deepcopy(F)
    rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(X.fan.RAYS)))
    slicedRayMatrix=slicematrix(rayMatrix)

    # Construct the cone made up of the residual (non-divisorial) rays and the divisor
    residRays = Array{Int64, 1}[]
    residStack = Int64[]
    c = 0
    # Add the divisor (which should be divisorial)
    if div[divisor] == 1
        c += 1
        push!(residRays, divisor)
        push!(residStack, F.stacks[encode(divisor)])
    else
        error("The given ray is non-divisorial")
    end
    # Append the other rays
    for i in cone
        ray = slicedRayMatrix[i]
        stack = F.stacks[encode(ray)]
        if div[ray] == 0 && isRelevant(ray, cone, F)
            c += 1
            push!(residRays, ray)
            push!(residStack, stack)
        end
    end
    residConeZ = Array(0:c-1)
    residCone = Array(1:c)
    # Create the stacky fan structure containing only the cone
    residRayMatrix = Array{Int64}(transpose(hcat(residRays...)))
    subconeFan = makeStackyFan(residRayMatrix, [residConeZ], residStack)

    # Count the number of residual rays which are dependent in the above cone
    # Skip the first ray, which is the divisor
    divIndex = 0
    for i in 2:c
        ray = residRays[i]
        if isRelevant(ray, residCone, subconeFan)
            divIndex += 1
        end
    end

    return divIndex
end

"""
    positiveDivIndexCone()

    Tests whether there exists a cone with non-zero divisorial index along the given divisor.
"""
function positiveDivIndexCone(divisor::Array{Int64,1}, F::StackyFan, div::Dict)
    # Go through each maximal cone containing the ray, and testing the divisorial index
    rayMatrix = convert(Array{Int64,2}, Array(Polymake.common.primitive(X.fan.RAYS)))
    maxconeList = convertIncidenceMatrix(F.fan.MAXIMAL_CONES)
    divisorIndex = getIndex(divisor, rayMatrix)
    for maxCone in maxconeList
        if divisorIndex in maxCone
            coneDivIndex = divAlongDivisor(maxCone, F, div, divisor)
            # If the divisorial index along the divisor is non-zero, return true
            if coneDivIndex > 0
                return true
            end
        end
    end
    return false
end

max_stars_repo_name: nbenabla/destackification

max_issues_repo_name: nbenabla/destackification

max_forks_repo_name: nbenabla/destackification


text:
include("../../matplotlib_settings.jl")

set_config(:tab_compensated_summation, false)

using ChargedParticleDynamics.GuidingCenter4d: plot_integral_error, plot_poincare_loop, plot_poincare_trajectories

const coll_title = "Guiding Center 4D Fast 1st Poincaré Integral Invariant"

XMIN = 1.5
XMAX = 2.0
YMIN = -0.25
YMAX = +0.25

include("../guiding_center_4d_poincare_invariant_1st.jl")


text:
"""
`SavingsAccount2` is same as `SavingsAccount` but uses Lazy's `@forward` macro.
We will use Delegation pattern to gain existing `Account` functionality.
"""
struct SavingsAccount2
    acct::Account
    interest_rate::Float64

    SavingsAccount2(account_number, balance, date_opened, interest_rate) = new(
        Account(account_number, balance, date_opened),
        interest_rate
    )
end

# Forward accessors and functions
@forward SavingsAccount2.acct account_number, balance, date_opened
@forward SavingsAccount2.acct deposit!, withdraw!

# The @forward macro only works with first argument
transfer!(from::SavingsAccount2, to::SavingsAccount2, amount::Real) = transfer!(
    from.acct, to.acct, amount)

# new accessor
interest_rate(sa::SavingsAccount2) = sa.interest_rate

# new behavior
function accrue_daily_interest!(sa::SavingsAccount2)
    interest = balance(sa.acct) * interest_rate(sa) / 365
    deposit!(sa.acct, interest)
end


text:
struct MultiFieldCellField{DS<:DomainStyle} <: CellField
  single_fields::Vector{<:CellField}
  domain_style::DS

  function MultiFieldCellField(single_fields::Vector{<:CellField})
    @assert length(single_fields) > 0
    f1 = first(single_fields)

    # Find a suitable domain style
    if any( map(i->DomainStyle(i)==ReferenceDomain(),single_fields) )
      domain_style = ReferenceDomain()
    else
      domain_style = PhysicalDomain()
    end

    new{typeof(domain_style)}(single_fields,domain_style)
  end
end

function CellData.get_data(f::MultiFieldCellField)
  s = """
  Function get_data is not implemented for MultiFieldCellField at this moment.
  You need to extract the individual fields and then evaluate them separatelly.

  If ever implement this, evaluating a `MultiFieldCellField` directly would provide,
  at each evaluation point, a tuple with the value of the different fields.
  """
  @notimplemented s
end

function CellData.get_triangulation(f::MultiFieldCellField)
  s1 = first(f.single_fields)
  trian = get_triangulation(s1)
  @check all(map(i->trian===get_triangulation(i),f.single_fields))
  trian
end
CellData.DomainStyle(::Type{MultiFieldCellField{DS}}) where DS = DS()
num_fields(a::MultiFieldCellField) = length(a.single_fields)
Base.getindex(a::MultiFieldCellField,i::Integer) = a.single_fields[i]
Base.iterate(a::MultiFieldCellField)  = iterate(a.single_fields)
Base.iterate(a::MultiFieldCellField,state)  = iterate(a.single_fields,state)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule FLAC_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("FLAC")
JLLWrappers.@generate_main_file("FLAC", UUID("1d38b3a6-207b-531b-80e8-c83f48dafa73"))
end  # module FLAC_jll


text:
const O = AbstractArray

struct GymPOMDP{S, A} <: POMDP{S, A, O}
    env::GymEnv{S} # OpenAIGym environment
    pixel_observations::Bool # Whether or not to use an image as the observation
    γ::Float64 # Discount factor
    actions::Vector{A} # List of actions (Discrete only for now)
    special_render # Function for rendering image (takes in env state), nothing if want OpenAI render function
    sign_reward # Whether or not to normalize reward using sign(r)
    frame_stack # number of steps to take and frames to stack to make up observation
    frame_stack_dim
end

function GymPOMDP(env::GymEnv; pixel_observations = false,
                    γ = 0.99, actions = nothing, special_render = nothing, sign_reward = false, frame_stack = 1, frame_stack_dim = 3)
    s0 = reset!(env)
    if isnothing(actions)
        actions = actionset(env.pyenv.action_space)
    end
    GymPOMDP(env, pixel_observations, γ, actions, special_render, sign_reward, frame_stack, frame_stack_dim)
end

function GymPOMDP(environment::Symbol; version::Symbol = :v0, kwargs...)
    env = GymEnv(environment, version)
    GymPOMDP(env; kwargs...)
end

function POMDPs.initialstate(mdp::GymPOMDP, rng::AbstractRNG = Random.GLOBAL_RNG)
    ImplicitDistribution((rng) -> reset!(mdp.env))
end
POMDPs.initialobs(mdp::GymPOMDP, s) = Deterministic(stack_obs(mdp, [rand(observation(mdp, s)) for _=1:mdp.frame_stack]))

POMDPs.actions(mdp::GymPOMDP) = mdp.actions
POMDPs.actionindex(mdp::GymPOMDP, a) = a

function POMDPs.isterminal(mdp::GymPOMDP, s)
    try
        mdp.env.pyenv._elapsed_steps >= mdp.env.pyenv._max_episode_steps && return false
    catch
        try
            mdp.env.pyenv.num_steps >= mdp.env.pyenv.steps && return false
        catch
            return mdp.env.done
        end
    end
end
POMDPs.discount(mdp::GymPOMDP) = mdp.γ

function POMDPs.observation(mdp::GymPOMDP, s)
    if mdp.pixel_observations
        o = isnothing(mdp.special_render) ? permutedims(Float32.(channelview(render(mdp))), [2,3,1]) : mdp.special_render(s)
    else
        o = Float32.(s)
    end
    Deterministic(o)
end

torgb(o) = collect(colorview(RGB, permutedims(reinterpret.(N0f8, o), [3,1,2])))

function init_mujoco_render()
    py"""
    from mujoco_py import GlfwContext
    GlfwContext(offscreen=True)
    """
end

function render(mdp::GymPOMDP, s, a = nothing; kwargs...)
    @assert s == mdp.env.state
    render(mdp; kwargs...)
end

render(mdp::GymPOMDP; kwargs...) = torgb(render(mdp.env; mode=:rgb_array, kwargs...))

stack_obs(mdp, os) = length(os) == 1 ? os[1] : cat(os..., dims = mdp.frame_stack_dim)

function aggregate_info!(infos, info)
    for k in keys(infos[1])
        info[k] =  mean([i[k] for i in filter((x)->haskey(x, k), infos)])
    end
end

function POMDPs.gen(mdp::GymPOMDP, s, a, rng::AbstractRNG = Random.GLOBAL_RNG; info=Dict())
    @assert s == mdp.env.state
    a_py = (a isa Int) ? a - 1 : a # Python indexes from 0
    rtot, os, sp = 0, [], nothing
    infos = []
    for i=1:mdp.frame_stack
        r, sp, ret_info = step!(mdp.env, a_py)
        o = rand(rng, observation(mdp, sp))
        rtot += r
        push!(os, o)
        push!(infos, ret_info)
    end
    aggregate_info!(infos, info)

    mdp.sign_reward && (rtot = sign(rtot))
    o = stack_obs(mdp, os)
    return (sp=sp, o=o, r=rtot)
end

max_stars_repo_name: ancorso/POMDPGym

max_issues_repo_name: ancorso/POMDPGym

max_forks_repo_name: ancorso/POMDPGym