text:
using Test
using CartesianMesh
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCell

function displacement(alpha, x::V) where {V<:AbstractVector}
    u1 = alpha * x[2] * sin(pi * x[1])
    u2 = alpha * (x[1]^3 + cos(pi * x[2]))
    return [u1, u2]
end

function displacement(alpha, x::M) where {M<:AbstractMatrix}
    npts = size(x)[2]
    return hcat([displacement(alpha, x[:, i]) for i = 1:npts]...)
end

function stress_field(lambda, mu, alpha, x)
    s11 =
        (lambda + 2mu) * alpha * pi * x[2] * cos(pi * x[1]) -
        lambda * alpha * pi * sin(pi * x[2])
    s22 =
        -(lambda + 2mu) * alpha * pi * sin(pi * x[2]) +
        lambda * alpha * pi * x[2] * cos(pi * x[1])
    s12 = alpha * mu * (3x[1]^2 + sin(pi * x[1]))
    return [s11, s22, s12]
end

function body_force(lambda, mu, alpha, x)
    b1 = alpha * (lambda + 2mu) * pi^2 * x[2] * sin(pi * x[1])
    b2 =
        -alpha * (6mu * x[1] + (lambda + mu) * pi * cos(pi * x[1])) +
        alpha * (lambda + 2mu) * pi^2 * cos(pi * x[2])
    return [b1, b2]
end

function boundary_nodeids(femesh)
    bn = CutCell.bottom_boundary_node_ids(femesh)
    rn = CutCell.right_boundary_node_ids(femesh)
    tn = CutCell.top_boundary_node_ids(femesh)
    ln = CutCell.left_boundary_node_ids(femesh)
    return unique!(vcat(bn, rn, tn, ln))
end

function add_cell_error_squared!(err, interpolater, exactsolution, cellmap, quad, detjac)
    for (p, w) in quad
        numsol = interpolater(p)
        exsol = exactsolution(cellmap(p))
        err .+= (numsol - exsol) .^ 2 * detjac * w
    end
end

function mesh_L2_error(
    nodalsolutions,
    exactsolution,
    nodalconnectivity,
    cellmaps,
    basis,
    errorquad,
)

    ndofs = size(nodalsolutions)[1]
    detjac = CutCell.determinant_jacobian(cellmaps[1])
    err = zeros(ndofs)
    interpolater = InterpolatingPolynomial(ndofs, basis)
    for (cellid, cellmap) in enumerate(cellmaps)
        nodeids = nodalconnectivity[:, cellid]
        elementsolutions = nodalsolutions[:, nodeids]
        update!(interpolater, elementsolutions)
        add_cell_error_squared!(
            err,
            interpolater,
            exactsolution,
            cellmap,
            errorquad,
            detjac,
        )
    end
    return sqrt.(err)
end

function apply_displacement_boundary_condition!(matrix, rhs, displacement_function, femesh)

    boundarynodeids = boundary_nodeids(femesh)
    nodalcoordinates = CutCell.nodal_coordinates(femesh)
    boundarynodecoordinates = nodalcoordinates[:, boundarynodeids]
    boundarydisplacement = displacement_function(boundarynodecoordinates)
    CutCell.apply_dirichlet_bc!(matrix, rhs, boundarynodeids, boundarydisplacement)
end

function linear_system(basis, quad, stiffness, femesh, bodyforcefunc)

    cellmaps = CutCell.cell_maps(femesh)
    nodalconnectivity = CutCell.nodal_connectivity(femesh)

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    cellmatrix = CutCell.bilinear_form(basis, quad, stiffness, cellmaps[1])
    CutCell.assemble_bilinear_form!(sysmatrix, cellmatrix, nodalconnectivity, 2)
    CutCell.assemble_body_force_linear_form!(
        sysrhs,
        bodyforcefunc,
        basis,
        quad,
        cellmaps,
        nodalconnectivity,
    )

    K = CutCell.make_sparse(sysmatrix, femesh)
    R = CutCell.rhs(sysrhs, femesh)

    return K, R
end

function stress_linear_system(basis, quad, stiffness, nodalsolution, mesh)
    nodalconnectivity = CutCell.nodal_connectivity(mesh)

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    jac = CutCell.jacobian(mesh)
    detjac = CutCell.determinant_jacobian(mesh)
    cellmatrix = CutCell.stress_cell_mass_matrix(basis, quad, detjac)

    CutCell.assemble_bilinear_form!(sysmatrix, cellmatrix, nodalconnectivity, 3)
    CutCell.assemble_stress_linear_form!(
        sysrhs,
        basis,
        quad,
        stiffness,
        nodalsolution,
        nodalconnectivity,
        jac,
    )

    matrix = CutCell.make_sparse_stress_operator(sysmatrix, mesh)
    rhs = CutCell.stress_rhs(sysrhs, mesh)

    return matrix, rhs
end

function error_for_num_elements(numelements, polyorder, numqp, lambda, mu, alpha)
    x0 = [0.0, 0.0]
    widths = [1.0, 1.0]
    nelements = [numelements, numelements]
    mesh = UniformMesh(x0, widths, nelements)
    stiffness = plane_strain_voigt_hooke_matrix(lambda, mu)

    basis = TensorProductBasis(2, polyorder)
    quad = tensor_product_quadrature(2, numqp)
    errorquad = tensor_product_quadrature(2, numqp + 2)

    femesh = CutCell.Mesh(mesh, basis)

    matrix, rhs =
        linear_system(basis, quad, stiffness, femesh, x -> body_force(lambda, mu, alpha, x))
    apply_displacement_boundary_condition!(matrix, rhs, x -> displacement(alpha, x), femesh)

    sol = matrix \ rhs
    nodalsolutions = reshape(sol, 2, :)

    nodalconnectivity = CutCell.nodal_connectivity(femesh)
    cellmaps = CutCell.cell_maps(femesh)

    err = mesh_L2_error(
        nodalsolutions,
        x -> displacement(alpha, x),
        nodalconnectivity,
        cellmaps,
        basis,
        errorquad,
    )
    return err
end

function stress_error_for_num_elements(numelements, polyorder, numqp, lambda, mu, alpha)
    x0 = [0.0, 0.0]
    widths = [1.0, 1.0]
    nelements = [numelements, numelements]
    mesh = UniformMesh(x0, widths, nelements)
    stiffness = plane_strain_voigt_hooke_matrix(lambda, mu)

    basis = TensorProductBasis(2, polyorder)
    quad = tensor_product_quadrature(2, numqp)
    errorquad = tensor_product_quadrature(2, numqp + 2)

    femesh = CutCell.Mesh(mesh, basis)

    matrix, rhs =
        linear_system(basis, quad, stiffness, femesh, x -> body_force(lambda, mu, alpha, x))
    apply_displacement_boundary_condition!(matrix, rhs, x -> displacement(alpha, x), femesh)

    sol = matrix \ rhs
    nodalsolutions = reshape(sol, 2, :)

    matrix, rhs = stress_linear_system(basis, quad, stiffness, sol, femesh)
    stress = reshape(matrix \ rhs, 3, :)

    nodalconnectivity = CutCell.nodal_connectivity(femesh)
    cellmaps = CutCell.cell_maps(femesh)

    err = mesh_L2_error(
        stress,
        x -> stress_field(lambda, mu, alpha, x),
        nodalconnectivity,
        cellmaps,
        basis,
        errorquad,
    )
    return err
end

function required_quadrature_order(polyorder)
    ceil(Int, 0.5 * (2polyorder + 1))
end

function mean(v)
    return sum(v) / length(v)
end

function convergence_rate(v,dx)
    return mean(diff(log.(v)) ./ diff(log.(dx)))
end

function convergence(nelements, polyorder, lambda, mu, alpha)
    numqp = required_quadrature_order(polyorder)
    err = zeros(length(nelements), 2)
    for (idx, nelement) in enumerate(nelements)
        err[idx, :] = error_for_num_elements(nelement, polyorder, numqp, lambda, mu, alpha)
    end
    dx = 1.0 ./ nelements
    u1err = err[:, 1]
    u2err = err[:, 2]

    u1rate = convergence_rate(u1err,dx)
    u2rate = convergence_rate(u2err,dx)

    return round(u1rate, digits = 3), round(u2rate, digits = 3)
end

function stress_convergence(nelements, polyorder, lambda, mu, alpha)
    numqp = required_quadrature_order(polyorder)
    err = zeros(length(nelements), 3)
    for (idx, nelement) in enumerate(nelements)
        err[idx, :] =
            stress_error_for_num_elements(nelement, polyorder, numqp, lambda, mu, alpha)
    end
    dx = 1.0 ./ nelements
    return [convergence_rate(err[:,i],dx) for i = 1:3]
end

lambda = 1.0
mu = 2.0
alpha = 0.01

numelements = [2, 4, 8, 16, 32, 64]
u1rate, u2rate = convergence(numelements, 1, lambda, mu, alpha)
# println("Convergence of linear elements : ", u1rate, "    ", u2rate)
@test isapprox(u1rate, 2.0, atol = 0.05)
@test isapprox(u2rate, 2.0, atol = 0.05)
stressrate = stress_convergence(numelements,1,lambda,mu,alpha)
@test all(stressrate .>= 1.0)

u1rate, u2rate = convergence(numelements, 2, lambda, mu, alpha)
# println("Convergence of quadratic elements : ", u1rate, "    ", u2rate)
@test isapprox(u1rate, 3.0, atol = 0.05)
@test isapprox(u2rate, 3.0, atol = 0.05)
stressrate = stress_convergence(numelements,2,lambda,mu,alpha)
@test all(stressrate .>= 1.9)

# u1rate, u2rate = convergence(numelements, 3, lambda, mu, alpha)
# # println("Convergence of cubic elements : ", u1rate, "    ", u2rate)
# @test isapprox(u1rate, 4.0, atol = 0.05)
# @test isapprox(u2rate, 4.0, atol = 0.05)
# stressrate = stress_convergence(numelements,3,lambda,mu,alpha)
# @test all(stressrate .>= 3.0)
#
# u1rate, u2rate = convergence(numelements, 4, lambda, mu, alpha)
# # println("Convergence of quartic elements : ", u1rate, "    ", u2rate)
# @test isapprox(u1rate, 5.0, atol = 0.05)
# @test isapprox(u2rate, 5.0, atol = 0.05)
# stressrate = stress_convergence(numelements,4,lambda,mu,alpha)
# @test all(stressrate .>= 3.9)


text:
#===============================================================================

    The main object defined in this file is: Block.
    It is the smallest `composite unit` needed for blocking. It provides
    a window view into `SamplingUnit` that restricts the range of view
    to some range i:j that corresponds to a chosen block.

===============================================================================#

const TVIEW{T} = SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}

"""
    mutable struct Block{L,TGP,TGPl,TW,TWn,TX}
        PP::TVIEW{TGP}
        P_last::TVIEW{TGPl} # view into a single element
        WW::TVIEW{TW}
        Wnr::TWn
        XX::TVIEW{TX}
        ll::Float64
        ll_history::Vector{Float64}
    end

The smallest "containerless" unit that provides a view into `SamplingUnit`
restricted to a range `i:j` of a block. `L` is an important flag that indicates
whether it is a terminal block or not.

# Fields
---
- `PP`: a vector of views into relevant `GuidProp`
- `P_last`: a view into a single (or none) `GuidProp` corresponding to the
            terminal sub-interval.
- `WW`: a vector of views into containers for sampled Wiener process
- `Wnr`: a flag for sampling Wiener processes
- `XX`: a vector of views into containers for a sampled process
- `ll`: a placeholder for computed log-likelihood
- `ll_history`: a history of computed log-likelihoods (useful for MCMC)

    function Block(
        u::SamplingUnit,
        range::UnitRange{Int64},
        last_block=false,
        ll_hist_len=0
    )

Base constructor.
"""
mutable struct Block{L,TGP,TGPl,TW,TWn,TX}
    PP::TVIEW{TGP}
    P_excl::TVIEW{TGP}
    P_last::TVIEW{TGPl} # view to a single element
    Pb_excl::TVIEW{TGPl}
    WW::TVIEW{TW}
    Wnr::TWn
    XX::TVIEW{TX}
    ll::Float64
    ll_history::Vector{Float64}

    function Block(
            u::SamplingUnit{TGP,TGPl,TW,TWn,TX},
            range::UnitRange{Int64},
            last_block=false,
            ll_hist_len=0
        ) where {TGP,TGPl,TW,TWn,TX}
        PP = view(u.PP, range[1]:(range[end]-!last_block)) # omit the last law
        P_excl = view(u.PP, (range[end]+last_block):range[end])
        P_last = view(u.PPb, (range[end]+last_block):range[end])
        Pb_excl = view(u.PPb, range[1]:(range[end]-!last_block))

        XX = view(u.XX, range)
        WW = view(u.WW, range)

        new{last_block,TGP,TGPl,TW,TWn,TX}(
            PP, P_excl, P_last, Pb_excl, WW, u.Wnr, XX, -Inf,
            Vector{Float64}(undef, ll_hist_len)
        )
    end
end

"""
    set_ll!(b::Block, i::Int, v)

Set the internal log-likelihood history field `b.ll_history[i]` with a value `v`
"""
set_ll!(b::Block, i::Int, v) = (b.ll_history[i] = v)

"""
    save_ll!(b::Block, i::Int)

Commit the current log-likelihood field `ll` to history `b.ll_history` at index
`i`.
"""
save_ll!(b::Block, i::Int) = set_ll!(b, i, b.ll)

"""
    GP.recompute_guiding_term!(b::Block)

Recompute the guiding term.
"""
function GP.recompute_guiding_term!(b::Block) end

function GP.recompute_guiding_term!(b::Block{false})
    recompute_guiding_term!(b.PP, b.P_last[1])
end

function GP.recompute_guiding_term!(b::Block{true})
    recompute_guiding_term!(b.PP)
end

"""
    find_W_for_X!(b::Block)

Compute the Wiener process `b.WW` that is needed for obtaining path `b.XX` under
the law stored in `b`.
"""
function find_W_for_X!(b::Block) end

function find_W_for_X!(b::Block{false})
    for i in eachindex(b.PP)
        DD.invsolve!(b.XX[i], b.WW[i], b.PP[i])
    end
    DD.invsolve!(b.XX[end], b.WW[end], b.P_last[1])
end

function find_W_for_X!(b::Block{true})
    for i in eachindex(b.PP)
        DD.invsolve!(b.XX[i], b.WW[i], b.PP[i])
    end
end

"""
    GP.loglikhd(b::Block)

Compute the log-likelihood evaluated at a sampled path.
"""
function GP.loglikhd(b::Block) end

GP.loglikhd(b::Block{false}) = (
    loglikhd(b.PP, b.XX) + loglikhd(b.P_last[1], b.XX[end])
)

GP.loglikhd(b::Block{true}) = loglikhd(b.PP, b.XX)

"""
    loglikhd!(b::Block)

Compute the log-likelihood evaluated at a sampled path and store the result in
an internal field `ll`.
"""
loglikhd!(b::Block) = (b.ll = loglikhd(b))

"""
    recompute_path!(b::Block, WW=b.WW; skip=0)

Recompute the path `b.XX` for a given wiener process `WW`.
"""
function recompute_path!(b::Block, WW=b.WW; skip=0) end

function recompute_path!(b::Block{false}, WW=b.WW; skip=0)
    success = _recompute_path!(b, WW; skip=skip)
    success || return false
    y1 = b.XX[end-1].x[end]
    success, ll = GP.solve_and_ll!(
        b.XX[end], WW[end], b.P_last[1], y1; skip=skip
    )
    b.ll += ll
    success
end

function recompute_path!(b::Block{true}, WW=b.WW; skip=0)
    _recompute_path!(b, WW; skip=skip)
end

function _recompute_path!(b::Block, WW=b.WW; skip=0)
    y1 = b.XX[1].x[1]
    ll_tot = loglikhd_obs(b.PP[1], y1)
    for i in 1:length(b.PP)
        success, ll = GP.solve_and_ll!(b.XX[i], WW[i], b.PP[i], y1; skip=skip)
        success || (b.ll = ll; return false)
        ll_tot += ll
        y1 = b.XX[i].x[end]
    end
    b.ll = ll_tot
    true
end

ext: jl


text:
using HetaSimulator, Plots, CSV
# heta_update()

### creating platform

p = load_platform(".")  # aproximately the same as "heta build ."
models(p)    # models dictionary
scenarios(p) # scenarios dictionary (empty)

### creating Scenario with Julia code

model = models(p)[:nameless]

# default
scenario0 = Scenario(model, tspan = (0, 10), events_active = [:sw1=>false])
res0 = sim(scenario0)
plot(res0)

scenario0 = Scenario(model, saveat = [1, 4, 10], events_active = [:sw1=>false])
res0 = sim(scenario0)
plot(res0)

# update dose
scenario1 = Scenario(
    model,
    parameters = [:dose=>100.],
    events_active = [:sw1=>false, :sw2=>true],
    tspan = (0.,50.),
    observables = [:C1, :C2, :v_abs, :v_el, :v_distr]
)
res1 = sim(scenario1)
plot(res1)

# load scenarios into platform
add_scenarios!(p, [:scn0 => scenario0, :scn1 => scenario1])
p

# simulate all scenarios
res_all = sim(p)
plot(res_all)

### Creating scenarios from tables

scenarios_df = read_scenarios("scenarios.csv")
add_scenarios!(p, scenarios_df)
p

# simulate all scenarios
res_all = sim(p)
plot(res_all)

# simulate selected scenarios
res_all = sim(p, scenarios = [:dose_1, :dose_10, :dose_100])
plot(res_all)

### plot settings

plot(res_all, vars = [:C1, :C2])
plot(res_all, yscale=:log10, ylims=(1e-3,1e2))

### export result into DataFrame

res_df_all = DataFrame(res_all)
CSV.write("res_df_all.csv", res_df_all)

res_df = DataFrame(res_all, vars=[:C1,:C2])
plot(res_df[!,:t], res_df[!,:C1])

### chain syntax

Scenario(models(p)[:nameless], tspan=(0,100)) |> sim |> plot


text:
module StructuredGraphicalLinearAlgebra
export StructuredLinearFunctions, FreeStructuredLinearFunctions,
  Ob, Hom, dom, codom, compose, ⋅, ∘, id, oplus, ⊕, mzero, braid,
  mcopy, Δ, delete, ◊, plus, zero, plus, +, scalar, antipode, adjoint,
  ℝ, munit, →, diag

import Base: +
import LinearAlgebra: adjoint, diag

using ...Catlab, ...Theories
import ...Theories:
  Ob, Hom, dom, codom, compose, ⋅, ∘, id, oplus, ⊕, mzero, braid, munit
using ..GraphicalLinearAlgebra
import ..GraphicalLinearAlgebra:
  mcopy, Δ, delete, ◊, plus, zero, plus, +, scalar, antipode, adjoint

# Theories
###########

""" Theory of *structured linear maps*

Structured matrices have some properties that allow us to compute with them faster
than general dense matrices. Morphisms in this category represent structured matrices.
"""
@theory StructuredLinearFunctions{Ob,Hom} <: LinearFunctions{Ob,Hom} begin
  munit()::Ob
  @op (ℝ) := munit

  diag(v)::(A→A) ⊣ (A::Ob, v::(ℝ()→A))
  upperdiag(v)::(A⊕ℝ() → A⊕ℝ()) ⊣ (A::Ob, v::(ℝ() → A))
  lowerdiag(v)::(A⊕ℝ() → A⊕ℝ()) ⊣ (A::Ob, v::(ℝ() → A))
  # lowerdiag(v, u)::(A⊕ℝ() → A⊕ℝ()) ⊣ (A::Ob, v::(ℝ() → A⊕ℝ()), u::(ℝ() → A))
  bidiag(v, u)::(A⊕ℝ() → A⊕ℝ()) ⊣ (A::Ob, v::(ℝ() → A⊕ℝ()), u::(ℝ() → A))
  tridiag(v, u, w)::(A⊕ℝ() → A⊕ℝ()) ⊣ (A::Ob, v::(ℝ() → A⊕ℝ()), u::(ℝ() → A), w::(ℝ() → A))
  symtridiag(v, u)::(A⊕ℝ() → A⊕ℝ()) ⊣ (A::Ob, v::(ℝ() → A⊕ℝ()), u::(ℝ() → A))
  #A is n-k and b is k in a n×n matrix with v on the kth diagonal
  upperdiag(A, v)::(A⊕B→A⊕B) ⊣ (A::Ob, B::Ob, v::(ℝ()→B))
  # A is n-k and b is k in a n×n matrix with v on the kth lower diagonal
  lowerdiag(A, v)::(A⊕B→A⊕B) ⊣ (A::Ob, B::Ob, v::(ℝ()→B))
  # Axioms
  lowerdiag == ⊟(upperdiag)
  upperdiag(A,b) == ◊(A) ⊕ diag(b) ⊕ mzero(A)  ⊣ (A::Ob, b::(ℝ()→B))
  lowerdiag(A,b) == mzero(A)  ⊕ diag(b) ⊕ ◊(A) ⊣ (A::Ob, b::(ℝ()→B))
  # by default upperdiagonal means 1-upperdiagonal
  upperdiag(b) == upperdiag(ℝ(), b) ⊣ (b::(ℝ()→B))
  # by default lowerdiagonal means 1-lowerdiagonal
  lowerdiag(b) == lowerdiag(ℝ(), b) ⊣ (b::(ℝ()→B))

  bidiag(a,b) == diag(a) + upperdiag(b) ⊣ (A::Ob, a::(ℝ()→A⊗ℝ()), b::(ℝ()→A))
  tridiag(a,b,c) == diag(a) + lowerdiag(b) + upperdiag(c) ⊣ (A::Ob, a::(ℝ()→A⊗ℝ()), b::(ℝ()→A), c::(ℝ()→A))
  symtridiag(a,b) == tridiagonal(a,b,b) ⊣ (A::Ob, a::(ℝ()→A⊗ℝ()), b::(ℝ()→A))
end

@syntax FreeStructuredLinearFunctions{ObExpr,HomExpr} StructuredLinearFunctions begin
  oplus(A::Ob, B::Ob) = associate_unit(new(A,B), mzero)
  oplus(f::Hom, g::Hom) = associate(new(f,g))
  compose(f::Hom, g::Hom) = new(f,g; strict=true) # No normalization!
end

end


text:
function collatz_chain_length(x::Int)
	if x==1
		return 1
	end
	if x%2==0
		x=Int(x/2)
	else
		x=x*3+1
	end
	return collatz_chain_length(x)+1
end

function collatz_length(x::Int)
	if x==1
		return 1
	end
	count=1;
	while x!=1
		count+=1
		if x%2==0
			x=Int(x/2)
			if x==1
				break
			end
		else
			x=x*3+1
		end

	end
	return count
end

function compute()
    const N=1000000
    ans =0
    index=1
	for i=1:N
		if ans<collatz_chain_length(i)
			index=i
			ans=collatz_chain_length(i)
		end
	end

	return index
end

function compute2()
	 const N=1000000
    ans =0
    index=1
	for i=1:N
		if ans<collatz_length(i)
			index=i
			ans=collatz_length(i)
		end
	end

	return index
end

print(compute2())


text:
# Autogenerated wrapper script for ImarisWriter_jll for aarch64-linux-gnu-cxx03
export libbpImarisWriter96

using HDF5_jll
using Lz4_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("ImarisWriter")
JLLWrappers.@declare_library_product(libbpImarisWriter96, "libbpImarisWriter96.so")
function __init__()
    JLLWrappers.@generate_init_header(HDF5_jll, Lz4_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libbpImarisWriter96,
        "lib/libbpImarisWriter96.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# reusable functionality to implement code execution

export split_kwargs, assign_args!

## macro tools

# split keyword arguments expressions into groups. returns vectors of keyword argument
# values, one more than the number of groups (unmatched keywords in the last vector).
# intended for use in macros; the resulting groups can be used in expressions.
function split_kwargs(kwargs, kw_groups...)
    kwarg_groups = ntuple(_->[], length(kw_groups) + 1)
    for kwarg in kwargs
        # decode
        Meta.isexpr(kwarg, :(=)) || throw(ArgumentError("non-keyword argument like option '$kwarg'"))
        key, val = kwarg.args
        isa(key, Symbol) || throw(ArgumentError("non-symbolic keyword '$key'"))

        # find a matching group
        group = length(kwarg_groups)
        for (i, kws) in enumerate(kw_groups)
            if key in kws
                group = i
                break
            end
        end
        push!(kwarg_groups[group], kwarg)
    end

    return kwarg_groups
end

# assign arguments to variables, handle splatting
function assign_args!(code, args)
    # handle splatting
    splats = map(arg -> Meta.isexpr(arg, :(...)), args)
    args = map(args, splats) do arg, splat
        splat ? arg.args[1] : arg
    end

    # assign arguments to variables
    vars = Tuple(gensym() for arg in args)
    map(vars, args) do var,arg
        push!(code.args, :($var = $arg))
    end

    # convert the arguments, compile the function and call the kernel
    # while keeping the original arguments alive
    var_exprs = map(vars, args, splats) do var, arg, splat
         splat ? Expr(:(...), var) : var
    end

    return vars, var_exprs
end


text:
@discord_object struct Integration
    id::Snowflake
    name::String
    type::String
    enabled::Bool
    syncing::Bool
    role_id::Snowflake
    enable_emoticons::Bool
    expire_behavior::IntegrationExpireBehavior.IntegrationExpireBehaviorEnum
    expire_grace_period::Int
    user::User
    account::IntegrationAccount
end


text:
# N-dimensional face with normal in the positive D direction
struct Face{N} <: AbstractFace{N}
    cells::Tuple{AbstractTree{N},AbstractTree{N}}
    face_direction::Int
    state

    Face{N}() where N = new()
    Face{N}(cells, face_direction, state) where N = new(cells, face_direction, state)
end

function cells(face::Face{N}; with_fine_siblings = true) where N
    if !with_fine_siblings || !at_refinement(face)
        return face.cells
    else
        # TODO return iterator
        if level(face) == level(face.cells[2])
            return cat(face.cells[1], reshape(siblings(face.cells[2], face.face_direction, 1), 1<<(N-1)), dims=1)
        else
            return cat(reshape(siblings(face.cells[1], face.face_direction, 2), 1<<(N-1)), face.cells[2], dims=1)
        end
    end
end

# Initialize a face, here Val indicates the side (1 or 2) relative to cell of this face
function Face(cell::AbstractTree{N}, other_cell::AbstractTree{N}, face_direction, side, state::Function) where N
    side == 1 ? Face{N}((other_cell, cell), face_direction, state(centroid((other_cell, cell), face_direction))) : Face{N}((cell, other_cell), face_direction, state(centroid((cell, other_cell), face_direction)))
end
function Face(cell::AbstractTree{N}, other_cell::AbstractTree{N}, face_direction, side, state::Nothing) where N
    side == 1 ? Face{N}((other_cell, cell), face_direction, state) : Face{N}((cell, other_cell), face_direction, state)
end

@inline at_boundary(face::Face) = !initialized(face.cells[1]) || !initialized(face.cells[2])
@inline at_boundary(face::AbstractFace) = false
@inline at_refinement(face::Face) = initialized(face.cells[1]) && initialized(face.cells[2]) && level(face.cells[1]) != level(face.cells[2])
@inline at_refinement(face::AbstractFace) = false
@inline regular(face::Face) = level(face.cells[1]) == level(face.cells[2])
@inline regular(face::AbstractFace) = false

@inline active(face::Face) = (active(face.cells[1]) && active(face.cells[2]))# || (at_boundary(face) && (active(face.cells[1]) || active(face.cells[2])))
@inline active(face::AbstractFace) = false

@inline initialized(face::AbstractFace) = false
@inline initialized(face::Face) = true

@inline level(face::Face) = max(level(face.cells[1]), level(face.cells[2]))

function collect_faces(tree::AbstractTree{N}; filter::Function = face -> true) where N
    faces = Vector{Face{N}}()
    collect_faces!(faces, tree, filter)
    return faces
end

function collect_faces!(faces::Vector{Face{N}}, tree::AbstractTree{N}, filter::Function) where N
    @inbounds for dir=1:N, side=1:2
        face = tree.faces[dir, side]
        if (side==1 || !regular(face)) && filter(face)
            push!(faces, face)
        end
    end
    for child ∈ tree.children
        collect_faces!(faces, child, filter)
    end
end


text:
using SnoopCompile
invalidations = @snoopr using EAGO

@show length(invalidations)
trees = invalidation_trees(invalidations)
#methinvs = trees[end]
#root = methinvs.backedges[end]
#ascend(root)

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
nearestpsd_insts(use_nat::Bool) = [
    [(side, use_completable, false, use_nat, true)
    for side in vcat(30, 50:50:200, 300:100:800)] # includes compile run
    for use_completable in (false, true)
    ]

insts = OrderedDict()
insts["nat"] = (nothing, nearestpsd_insts(true))
insts["ext"] = (:SOCExpPSD, nearestpsd_insts(false))
return (NearestPSDJuMP, insts)


text:
module TeaFiles

using Dates

# Times in TeaFiles can either be a Julia DateTime or just a real number.
TimeLike = Union{Real, DateTime}

module Header
include("header/core.jl")
include("header/read.jl")
include("header/write.jl")
include("header/create.jl")
end  # Header

module Body
include("body/read.jl")
end  # Body

include("api.jl")

end  # TeaFiles


text:
using Images, ImageView, MLDatasets

train_x, train_y = MNIST.traindata();
test_x,  test_y  = MNIST.testdata();

# preview the images from train_x dataset
preview_img = zeros(size(train_x, 1), 0)

for i = 1:10
    preview_img = hcat(preview_img, train_x[:, :, i])
end

imshow(Gray.(preview_img))

# creating first neural network
using MXNet

train_length = 50000
validation_length = 10000

train_data_array  = mx.zeros((size(train_x, 1, 2)..., train_length...));
train_label_array = mx.zeros(train_length);

validation_data_array  = mx.zeros((size(train_x, 1, 2)..., validation_length...));
validation_label_array = mx.zeros(validation_length);

# The number of records we send to the training should be at least number of outcome
for idx = 1:train_length
    train_data_array[idx:idx] = reshape(train_x[:, :, idx], (size(train_x, 1, 2)..., 1...))
    train_label_array[idx:idx] = train_y[idx]
end

for idx = 1:validation_length
    validation_data_array[idx:idx] = reshape(train_x[:, :, train_length + idx], (size(train_x, 1, 2)..., 1...))
    validation_label_array[idx:idx] = train_y[train_length + idx]
end

arch = @mx.chain mx.Variable(:data) =>
  mx.FullyConnected(num_hidden=64) =>
  mx.FullyConnected(num_hidden=10) =>
  mx.SoftmaxOutput(mx.Variable(:label))

nnet = mx.FeedForward(arch, context = mx.cpu())

train_data_provider = mx.ArrayDataProvider(:data => train_data_array, :label => train_label_array, batch_size = 1000);
validation_data_provider = mx.ArrayDataProvider(:data => validation_data_array, :label => validation_label_array, batch_size = 1000);

mx.fit(nnet, mx.ADAM(), train_data_provider, eval_data = validation_data_provider, n_epoch = 50, callbacks = [mx.speedometer()]);

# more complicated network
arch = @mx.chain mx.Variable(:data) =>
  mx.FullyConnected(num_hidden=128) =>
  mx.Activation(act_type=:relu) =>
  mx.FullyConnected(num_hidden=64) =>
  mx.Activation(act_type=:relu) =>
  mx.FullyConnected(num_hidden=10) =>
  mx.SoftmaxOutput(mx.Variable(:label))

nnet = mx.FeedForward(arch, context = mx.cpu())

train_data_provider = mx.ArrayDataProvider(:data => train_data_array, :label => train_label_array, batch_size = 1000);
validation_data_provider = mx.ArrayDataProvider(:data => validation_data_array, :label => validation_label_array, batch_size = 1000);

mx.fit(nnet, mx.ADAM(), train_data_provider, eval_data = validation_data_provider, n_epoch = 30, callbacks = [mx.speedometer()]);

# predicting
data_array  = mx.zeros((size(test_x, 1, 2)..., 10));
mx.copy!(data_array, test_x[:, :, 1:10]);

data_provider = mx.ArrayDataProvider(:data => data_array);
results = round(mx.predict(nnet, data_provider; verbosity = 0), 2)

preview_img = zeros(size(test_x, 1), 0)

for i = 1:10
    preview_img = hcat(preview_img, test_x[:, :, i])
end

map(x -> findmax(results[:, x])[2] - 1, 1:10)
imshow(preview_img)


text:
# Calculator in Julia

function calculator()
    println("------- Welcome to the Calculator app written in Julia -------")
    while true
        # Print the Calculator app Menu|
        println("==== Menu ====")
        println("1. Addition ➕")
        println("2. Subtraction ➖")
        println("3  Multiplication ✖")
        println("4. Division ➗")
        println("5. Exit the Calculator app ↗")
        operation = parse(Float64, readline())
        # Convert the string value input to a float

        # Depending on the user's choice, do the corresponding operation
        if operation == 1
            print("Enter the first number: ")
            num1 = parse(Float64, readline())
            print("Enter the second number: ")
            num2 = parse(Float64, readline())
            println("Output: $(num1 + num2)")
        elseif operation == 2
            print("Enter the first number: ")
            num1 = parse(Float64, readline())
            print("Enter the second number: ")
            num2 = parse(Float64, readline())
            println("Output: $(num1 - num2)")
        elseif operation == 3
            print("Enter the first number: ")
            num1 = parse(Float64, readline())
            print("Enter the second number: ")
            num2 = parse(Float64, readline())
            println("Output: $(num1 * num2)")
        elseif operation == 4
            print("Enter the first number: ")
            num1 = parse(Float64, readline())
            print("Enter the second number: ")
            num2 = parse(Float64, readline())
            println("Output: $(num1 / num2)")
        elseif operation == 5
            println("Exiting calculator")
            break # Stop the while loop if the user decides to exit
        else
            # If the user enters another number, then ask to try again
            println("Invalid input, please try again.")
        end

    end
end


text:
################################################################################
#
#  The NumFieldMor type
#
################################################################################

mutable struct LocalFieldMor{S, T, U, V, W} <: Map{S, T, HeckeMap, LocalFieldMor}
  header::MapHeader{S, T}
  image_data::U
  inverse_data::V
  absolute_basis::Vector{W}
  absolute_basis_matrix_image::Generic.MatSpaceElem{padic}
  rref::Tuple{Generic.MatSpaceElem{padic}, Generic.MatSpaceElem{padic}}
  pivots_of_rref::Vector{Int}

  function LocalFieldMor{S, T, U, V}() where {S, T, U, V}
    z = new{S, T, U, V, elem_type(S)}()
    return z
  end

  function LocalFieldMor(K::S, L::T) where {S <: Union{FlintPadicField, FlintQadicField, LocalField}, T <: Union{FlintPadicField, FlintQadicField, LocalField}}
    z = new{typeof(K), typeof(L), map_data_type(K, L), map_data_type(L, K), elem_type(K)}()
    z.header = MapHeader(K, L)
    return z
  end

  function LocalFieldMor{S, T, U, V}(h::MapHeader{S, T}, i::U, p::V) where {S, T, U, V}
    z = new{S, T, U, V, elem_type(S)}(h, i, p)
    return z
  end
end

################################################################################
#
#   Identity
#
################################################################################

function id_hom(K::T) where T <: Union{LocalField, FlintQadicField}
  z = LocalFieldMor{typeof(K), typeof(K), map_data_type(K, K), map_data_type(K, K)}(MapHeader(K, K), map_data(K, K, true), map_data(K, K, true))
end

################################################################################
#
#  Morphism type
#
################################################################################

morphism_type(K::LocalField) = morphism_type(typeof(K), typeof(K))

morphism_type(K::FlintQadicField) = morphism_type(typeof(K), typeof(K))

morphism_type(K::Type{T}) where T <: Union{LocalField, FlintQadicField} = morphism_type(T, T)

morphism_type(K::S, L::T) where {S <: Union{LocalField, FlintQadicField}, T <: Union{LocalField, FlintQadicField}} = morphism_type(S, T)

morphism_type(::Type{S}, ::Type{T}) where {S <: Union{LocalField, FlintQadicField}, T <: Union{LocalField, FlintQadicField}} = LocalFieldMor{S, T, map_data_type(S, T), map_data_type(T, S), elem_type(S)}
################################################################################
#
#  The hom function
#
################################################################################

function hom(K::S, L::T, x...; inverse = nothing,
                               check::Bool = true,
                               compute_inverse = false) where {S <: Union{LocalField, FlintQadicField},
                                                               T <: Union{LocalField, FlintQadicField}}
  header = MapHeader(K, L)

  # Check if data defines a morphism
  image_data = map_data(K, L, x..., check = check)

  if inverse !== nothing
    # Check if data defines a morphism
    # This goes through _validata_data, since we don't want to split the
    # argument if for example the argument is a Vector
    inverse_data = map_data(L, K, inverse, check = check)

    z = LocalFieldMor{S, T, typeof(image_data),
                       typeof(inverse_data)}(header, image_data, inverse_data)

  else
    z = LocalFieldMor{S, T, typeof(image_data), map_data_type(L, K)}()
  end

  z.header = header
  z.image_data = image_data

  if compute_inverse
    _assure_has_inverse_data(z)
  end

  return z
end

################################################################################
#
#  Image and preimage function
#
################################################################################

function image(f::LocalFieldMor, x)
  @assert domain(f) === parent(x)
  return image(f.image_data, codomain(f), x)
end

################################################################################
#
#  Types to store the data for the maps
#
################################################################################

mutable struct MapDataFromPadicField{S}
  codomain::S
end

function map_data_type(T::Type{FlintPadicField}, L::Type{S}) where S <: Union{LocalField, FlintQadicField}
  MapDataFromPadicField{S}
end

map_data_type(T::FlintPadicField, L::Union{LocalField, FlintQadicField}) = map_data_type(typeof(T), typeof(L))

# Test if data u, v specfiying a map K -> L define the same morphism
function _isequal(K, L, u::MapDataFromPadicField{T}, v::MapDataFromPadicField{T}) where {T, S}
  return true
end

function image(f::MapDataFromPadicField, L, y)
  return L(y)
end

function map_data(K::FlintPadicField, L, ::Bool)
  return MapDataFromPadicField{typeof(L)}(L)
end

function map_data(K::FlintPadicField, L; check::Bool = true)
  return MapDataFromPadicField{typeof(L)}(L)
end

# From LocalField into something
mutable struct MapDataFromLocalField{T, S}
  prim_image::T
  base_field_map_data::S
  isid::Bool

  function MapDataFromLocalField{T, S}(x::T, y::S) where {T, S}
    z = new{T, S}(x, y, false)
    return z
  end

  function MapDataFromLocalField{T, S}(x::Bool) where {T, S}
    @assert x
    z = new{T, S}()
    z.isid = true
    return z
  end
end

function map_data_type(T::Type{<:LocalField}, L::Type{<:Union{LocalField, FlintQadicField}})
  MapDataFromLocalField{elem_type(L), map_data_type(base_field_type(T), L)}
end

map_data_type(T::LocalField, L::Union{LocalField, FlintQadicField}) = map_data_type(typeof(T), typeof(L))

# Test if data u, v specfiying a map K -> L define the same morphism
function _isequal(K, L, u::MapDataFromLocalField{T, S}, v::MapDataFromLocalField{T, S}) where {T, S}
  if u.isid && v.isid
    return true
  end

  return image(u, L, gen(K)) == image(v, L, gen(K)) &&
         _isequal(base_field(K), L, u.base_field_map_data, v.base_field_map_data)
end

function image(f::MapDataFromLocalField, L, y)
  f.isid && return L(y)
  # TODO: Cache the polynomial ring
  Ly, = PolynomialRing(L, "y", cached = false)
  z = map_coefficients(t -> image(f.base_field_map_data, L, t), y.data, parent = Ly)
  return evaluate(z, f.prim_image)
end

function map_data(K::LocalField, L, ::Bool)
  z = MapDataFromLocalField{elem_type(L), map_data_type(base_field(K), L)}(true)
  z.base_field_map_data = map_data(base_field(K), L, true)
  return z
end

function map_data(K::LocalField, L, x...; check = true)
  z = map_data(base_field(K), L, Base.front(x)...; check = check)

  local yy::elem_type(L)

  if Base.last(x) isa LocalFieldMor
    domain(Base.last(x)) !== K && error("Data does not define a morphism")
    _y = image_primitive_element(Base.last(x))
    if parent(_y) === L
      yy = _y
    else
      yy = L(_y)::elem_type(L)
    end
  else
    y = Base.last(x)
    if parent(y) === L
      yy = y
    else
      yy = L(y)::elem_type(L)
    end
  end

  if check
    y = evaluate(map_coefficients(w -> image(z, L, w), defining_polynomial(K), cached = false), yy)
    !iszero(y) && error("Data does not define a morphism")
  end

  @assert typeof(yy) == elem_type(L)
  @assert typeof(z) == map_data_type(base_field(K), L)

  return MapDataFromLocalField{typeof(yy), typeof(z)}(yy, z)
end

# From QadicField into something
mutable struct MapDataFromQadicField{T}
  prim_image::T
  isid::Bool

  function MapDataFromQadicField{T}(x::T) where {T}
    z = new{T}(x, false)
    return z
  end

  function MapDataFromQadicField{T}(x::Bool) where {T}
    @assert x
    z = new{T}()
    z.isid = true
    return z
  end
end

function map_data_type(T::Type{FlintQadicField}, L::Type{<:Union{LocalField, FlintQadicField}})
  MapDataFromQadicField{elem_type(L)}
end

map_data_type(::FlintQadicField, L::Union{LocalField, FlintQadicField}) = map_data_type(FlintQadicField, typeof(L))

# Test if data u, v specfiying a map K -> L define the same morphism
function _isequal(K, L, u::MapDataFromQadicField{T}, v::MapDataFromQadicField{T}) where {T}
  if u.isid && v.isid
    return true
  end

  return image(u, L, gen(K)) == image(v, L, gen(K))
end

function image(f::MapDataFromQadicField, L, y)
  f.isid && return L(y)
  # TODO: Cache the polynomial ring
  Qpx = PolynomialRing(base_field(parent(y)), cached = false)[1]
  return evaluate(Qpx(y), f.prim_image)
end

function map_data(K::FlintQadicField, L, f::LocalFieldMor; check::Bool = true)
  if check
    K !== domain(f) && error("Data does not define a morphism")
  end

  z = L(image_primitive_element(f))

  return map_data(K, L, z; check = false)
end

function map_data(K::FlintQadicField, L, ::Bool)
  return MapDataFromQadicField{elem_type(L)}(true)
end

map_data(K::FlintQadicField, L; check = true) = map_data(K, L, L(gen(K)), check = check)

function map_data(K::FlintQadicField, L, x; check = true)

  if check
    y = evaluate(defining_polynomial(K), x)
    !iszero(y) && error("Data does not define a morphism")
  end

  return MapDataFromQadicField{typeof(x)}(x)
end

function image_primitive_element(f::LocalFieldMor)
  if f.image_data.isid
    return codomain(f)(gen(domain(f)))
  end
  return f.image_data.prim_image
end

################################################################################
#
#  Preimage computation
#
# ##############################################################################

function _assert_has_preimage_data(f::LocalFieldMor)
  if isdefined(f, :absolute_basis_matrix_image)
    return nothing
  end
  K = domain(f)
  L = codomain(f)
  b = absolute_basis(K)
  d = absolute_degree(K)
  n = absolute_degree(L)
  M = zero_matrix(absolute_base_field(K), n, d)
  for i in 1:d
    c = f(b[i])
    cc = absolute_coordinates(c)
    for j in 1:length(cc)
      M[j, i] = cc[j]
    end
  end

  r, R, U =  _rref_with_trans(M)
  pivots = _get_pivots_ut(R)

  f.absolute_basis_matrix_image = M
  f.absolute_basis = b
  f.pivots_of_rref = pivots
  f.rref = R, U

  return nothing
end

function haspreimage(f::LocalFieldMor, g::Union{qadic, LocalFieldElem})
  if isdefined(f, :inverse_data)
    return true, image(f.inverse_data, domain(f), g)
  end
  @assert parent(g) === codomain(f)
  K = domain(f)
  d = absolute_degree(K)
  cc = absolute_coordinates(g)
  _assert_has_preimage_data(f)
  fl, s = can_solve_given_rref(f.rref[1], f.rref[2], f.pivots_of_rref, cc)
  if !fl
    return false, zero(K)
  else
    b = f.absolute_basis
    # This is suboptimal
    prim_preimg = zero(K)
    for i = 1:d
      prim_preimg += s[i, 1] * b[i]
    end
    return true, prim_preimg
  end
end

function preimage(f::LocalFieldMor, g::Union{qadic, LocalFieldElem})
  fl, y = haspreimage(f, g)
  @assert fl
  return y
end

################################################################################
#
#  Computation of the inverse (data)
#
################################################################################

function _assure_has_inverse_data(f::LocalFieldMor)
  if isdefined(f, :inverse_data)
    return nothing
  else
    pr = _compute_inverse_data(f, domain(f), codomain(f))
    f.inverse_data = pr
    return nothing
  end
end

function inv(f::LocalFieldMor{S, T}) where {S, T}
  _assure_has_inverse_data(f)
  pr = f.inverse_data
  hd = MapHeader(codomain(f), domain(f))
  g = NumFieldMor{T, S, map_data_type(T, S), map_data_type(S, T)}(hd, pr, f.image_data)
  return g
end

function _compute_inverse_data(f#= image data =#, K, L::LocalField)
  return _compute_inverse_data(f#= image data =#, K, L, L)
end

function _compute_inverse_data(f#= image data =#, K, L::QadicField)
  return _compute_inverse_data(f#= image data =#, K, L, L)
end

function _compute_inverse_data(f#= image data =#, K, LL, L::LocalField)
  g = LL(gen(L))
  fl, prim_preimg = haspreimage(f, LL(g))
  @assert fl
  return MapDataFromLocalField{typeof(prim_preimg)}(prim_preimg)
end

function _compute_inverse_data(f#= image data =#, K, LL, L::FlintQadicField)
  g = LL(gen(L))
  fl, prim_preimg = haspreimage(f, LL(g))
  @assert fl
  return MapDataFromQadicField{typeof(prim_preimg)}(prim_preimg)
end

################################################################################
#
#  Composition
#
################################################################################

# f : K -> L, g : L -> M
function _compose(f::MapDataFromQadicField, g#= map data =#, K, L, M)
  return map_data_type(K, M)(image(g, M, image(f, L, gen(K))))
end

function _compose(f::MapDataFromLocalField, g#= map data =#, K, L, M)
  return map_data_type(K, M)(image(g, M, image(f, L, gen(K))),
                             _compose(f.base_field_map_data, g, base_field(K), L, M))
end

function Base.:(*)(f::LocalFieldMor, g::LocalFieldMor)
  @req codomain(f) === domain(g) "Composition: Maps are not compatible"
  z = LocalFieldMor(domain(f), codomain(g))
  z.image_data = _compose(f.image_data, g.image_data, domain(f), codomain(f), codomain(g))
  if isdefined(f, :inverse_data) && isdefined(g, :inverse_data)
    z.inverse_data = _compose(g.inverse_data, f.inverse_data, codomain(g), domain(g), domain(f))
  end
  return z
end

################################################################################
#
#  Equality
#
################################################################################

function Base.:(==)(u::T, v::T) where T <: LocalFieldMor
  if (domain(u) != domain(v)) || (codomain(u) != codomain(v))
    return false
  end
  return _isequal(domain(u), codomain(u), u.image_data, v.image_data)
end

################################################################################
#
#  Powering
#
################################################################################

function ^(f::LocalFieldMor, b::Int)
  K = domain(f)
  @assert K == codomain(f)
  d = absolute_degree(K)
  b = mod(b, d)
  if b == 0
    return id_hom(K)
  elseif b == 1
    return f
  else
    bit = ~((~UInt(0)) >> 1)
    while (UInt(bit) & b) == 0
      bit >>= 1
    end
    z = f
    bit >>= 1
    while bit != 0
      z = z * z
      if (UInt(bit) & b) != 0
        z = z * f
      end
      bit >>= 1
    end
    return z
  end
end

^(a::LocalFieldMor, n::Union{Integer, fmpz})  = _generic_power(a, n)

################################################################################
#
#  Hashing
#
################################################################################

function Base.hash(f::MapDataFromQadicField, K, L, h::UInt)
  if f.isid
    return xor(hash(L, h), hash(K, h))
  else
    return hash(f.prim_image, h)
  end
end

function Base.hash(f::MapDataFromLocalField, K, L, h::UInt)
  if f.isid
    h = xor(hash(L, h), hash(K, h))
  else
    h = hash(f.prim_image, h)
  end
  h = hash(f.base_field_map_data, base_field(K), L, h)
  return h
end

Base.hash(f::LocalFieldMor, h::UInt) = hash(f.image_data, domain(f), codomain(f), h)

################################################################################
#
#  Restriction
#
################################################################################

function restrict(f::LocalFieldMor, K::Union{FlintQadicField, LocalField})
  k = domain(f)
  return hom(K, k, k(gen(k)))*f
end

################################################################################
#
#   To make the automorphism group work
#
################################################################################

function GrpGenToNfMorSet(G::GrpGen, K::T) where T <: Union{LocalField, FlintQadicField}
  return GrpGenToNfMorSet(automorphisms(K), G, NfMorSet(K))
end

function GrpGenToNfMorSet(G::GrpGen, aut::Vector{S}, K::T) where {S <: LocalFieldMor, T <: Union{LocalField, FlintQadicField}}
  return GrpGenToNfMorSet(aut, G, NfMorSet(K))
end


text:
module Assemblers

using Gridap

using Gridap.Helpers

using SparseArrays
using SparseMatricesCSR

export Assembler
export SparseMatrixAssembler
export assemble
export assemble!
export sparse_from_coo

"""
Abstract assembly operator
Parametrized by the type of returned matrix and vector
"""
abstract type Assembler{M<:AbstractMatrix,V<:AbstractVector} end

"""
Assembly of a vector allocating output
"""
function assemble(
  ::Assembler{M,V},
  ::Vararg{Tuple{<:CellVector,<:CellNumber}})::V where {M,V}
  @abstractmethod
end

"""
Assembly of a matrix allocating output
"""
function assemble(
  ::Assembler{M,V},
  ::Vararg{Tuple{<:CellMatrix,<:CellNumber,<:CellNumber}})::M where {M,V}
  @abstractmethod
end

"""
In-place assembly of a vector (allows optimizations)
"""
function assemble!(
  ::V,
  ::Assembler{M,V},
  ::Vararg{Tuple{<:CellVector,<:CellNumber}})::V where {M,V}
  @abstractmethod
end

"""
In-place assembly of a matrix (allows a LOT of optimizations)
"""
function assemble!(
  ::M,
  ::Assembler{M,V},
  ::Vararg{Tuple{<:CellMatrix,<:CellNumber,<:CellNumber}})::M where {M,V}
  @abstractmethod
end

function assemble(a::Assembler,cv::CellVector)
  l = length(cv)
  ide = IdentityCellNumber(Int,l)
  assemble(a,(cv,ide))
end

function assemble(a::Assembler,cv::CellMatrix)
  l = length(cv)
  ide = IdentityCellNumber(Int,l)
  assemble(a,(cv,ide,ide))
end

function assemble!(r,a::Assembler,cv::CellVector)
  l = length(cv)
  ide = IdentityCellNumber(Int,l)
  assemble!(r,a,(cv,ide))
end

function assemble!(r,a::Assembler,cv::CellMatrix)
  l = length(cv)
  ide = IdentityCellNumber(Int,l)
  assemble!(r,a,(cv,ide,ide))
end

"""
Assembler that produces SparseMatrices from the SparseArrays and SparseMatricesCSR packages
"""
struct SparseMatrixAssembler{E,M} <: Assembler{M,Vector{E}}
  testfesp::FESpace
  trialfesp::FESpace

  function SparseMatrixAssembler(
    ::Type{M},
    testfesp::FESpace{D,Z,T},
    trialfesp::FESpace{D,Z,T}) where {M<:AbstractSparseMatrix,D,Z,T}
    E = eltype(T)
    new{E,M}(testfesp,trialfesp)
  end
end

function SparseMatrixAssembler(test::FESpace{D,Z,T}, trial::FESpace{D,Z,T}) where {D,Z,T}
  SparseMatrixAssembler(SparseMatrixCSC, trial,test)
end

function assemble(
  this::SparseMatrixAssembler{E},
  vals::Vararg{Tuple{<:CellVector,<:CellNumber}}) where {E}

  n = num_free_dofs(this.testfesp)
  vec = zeros(E,n)
  assemble!(vec,this,vals...)
  vec
end

function assemble!(
  vec::Vector{E},
  this::SparseMatrixAssembler{E},
  allvals::Vararg{Tuple{<:CellVector,<:CellNumber}}) where {E}

  vec .= zero(E)
  rows = celldofids(this.testfesp)
  for (vals, cellids) in allvals
    _vals = apply_constraints(this.testfesp, vals, cellids)
    _rows = reindex(rows,cellids)
    _assemble_vector!(vec, _vals, _rows)
  end
end

function _assemble_vector!(vec,vals,rows)
  for (rows_c,vals_c) in zip(rows,vals)
    for (i,gid) in enumerate(rows_c)
      if gid > 0
        vec[gid] += vals_c[i]
      end
    end
  end
end

function assemble(
  this::SparseMatrixAssembler{E,M},
  allvals::Vararg{Tuple{<:CellMatrix,<:CellNumber,<:CellNumber}}) where {E,M}

  I = Int
  aux_row = I[]; aux_col = I[]; aux_val = E[]

  _rows_m = celldofids(this.testfesp)
  _cols_m = celldofids(this.trialfesp)
  for (vals, cellids_row, cellids_col) in allvals
    _vals = apply_constraints_rows(this.testfesp, vals, cellids_row)
    rows_m = reindex(_rows_m, cellids_row)
    vals_m = apply_constraints_cols(this.trialfesp, _vals, cellids_col)
    cols_m = reindex(_cols_m, cellids_col)
    _assemble_sparse_matrix_values!(M,
      aux_row,aux_col,aux_val,vals_m,rows_m,cols_m)
  end
  num_rows = num_free_dofs(this.testfesp)
  num_cols = num_free_dofs(this.trialfesp)
  finalize_coo!(M,aux_row,aux_col,aux_val,num_rows,num_cols)
  sparse_from_coo(M,aux_row,aux_col,aux_val)
end

function _assemble_sparse_matrix_values!(::Type{M},aux_row,aux_col,aux_val,vals,rows,cols) where {M}
  for (rows_c, cols_c, vals_c) in zip(rows,cols,vals)
     for (j,gidcol) in enumerate(cols_c)
       if gidcol > 0
        for (i,gidrow) in enumerate(rows_c)
          if gidrow > 0
            push_coo!(M,aux_row, aux_col, aux_val, gidrow, gidcol, vals_c[i,j])
          end
        end
      end
    end
  end
end

function assemble!(
  mat::AbstractSparseMatrix{E,Int},
  this::SparseMatrixAssembler{E,M},
  vals::Vararg{Tuple{<:CellMatrix,<:CellNumber,<:CellNumber}}) where {E,M}
  # This routine can be optimized a lot taking into a count the sparsity graph of mat
  # For the moment we create an intermediate matrix and then transfer the nz values
  m = assemble(this,vals...)
  nonzeros(mat) .= nonzeros(m)
end

function sparse_from_coo(::Type{<:SparseMatrixCSC}, args...)
  sparse(args...)
end

function sparse_from_coo(::Type{<:SparseMatrixCSR}, args...)
  sparsecsr(args...)
end

function sparse_from_coo(::Type{<:SymSparseMatrixCSR}, args...)
  symsparsecsr(args...)
end

end # module


text:
macro load_turing_fmi()
    expr = quote
        @info "Loading foreign model interface to \u001b[3m\u001b[34;1mTuring.jl\u001b[0m\n\n          \u001b[34;1mhttps://turing.ml/dev/\n\n "

        import Jaynes: has_top, get_top, has_sub, get_sub, get_score, collect!
        using Turing
    end

    expr = MacroTools.prewalk(unblock ∘ rmlines, expr)
    esc(expr)
end


text:
"""
    DisjointSetsWithBirth{T}

Almost identical to `DataStructures.IntDisjointSets`, but keeps track of vertex birth times
and birth vertices.
Has no `num_groups` method.
"""
struct DisjointSetsWithBirth{
    I<:AbstractArray,A<:AbstractArray,B<:AbstractArray,C<:AbstractArray
}
    vertices::I
    parents::A
    ranks::B
    births::C

    function DisjointSetsWithBirth(vertices, births)
        parents = collect(vertices)
        ranks = similar(parents, Int)
        ranks .= 0
        births = collect(zip(births, vertices))
        return new{typeof(vertices),typeof(parents),typeof(ranks),typeof(births)}(
            vertices, parents, ranks, births
        )
    end
end

Base.length(s::DisjointSetsWithBirth) = length(s.parents)

function DataStructures.find_root!(s::DisjointSetsWithBirth, x)
    parents = s.parents
    p = parents[x]
    @inbounds if parents[p] != p
        parents[x] = p = find_root!(s, p)
    end
    return p
end

"""
    find_leaves!(s::DisjointSetsWithBirth, x)

Find all leaves below `x`, i.e. vertices that have `x` as root.
"""
function find_leaves!(s::DisjointSetsWithBirth, x)
    leaves = eltype(s.parents)[]
    for i in s.vertices
        find_root!(s, i) == x && push!(leaves, i)
    end
    return leaves
end

function Base.union!(s::DisjointSetsWithBirth, x, y)
    parents = s.parents
    xroot = find_root!(s, x)
    yroot = find_root!(s, y)
    return xroot != yroot ? root_union!(s, xroot, yroot) : xroot
end

function DataStructures.root_union!(s::DisjointSetsWithBirth, x, y)
    parents = s.parents
    ranks = s.ranks
    births = s.births
    @inbounds xrank = ranks[x]
    @inbounds yrank = ranks[y]

    if xrank < yrank
        x, y = y, x
    elseif xrank == yrank
        ranks[x] += 1
    end
    @inbounds parents[y] = x
    @inbounds births[x] = min(births[x], births[y])
    return x
end

birth(dset::DisjointSetsWithBirth, i) = dset.births[i]

function interval(dset::DisjointSetsWithBirth, filtration, vertex, edge, cutoff, reps)
    birth_time, birth_vertex = birth(dset, vertex)
    death_time = isnothing(edge) ? Inf : birth(edge)
    if death_time - birth_time > cutoff
        birth_simplex::simplex_type(filtration, 0) = simplex(
            filtration, Val(0), (birth_vertex,)
        )
        if reps
            rep = (;
                representative=sort!([
                    simplex(filtration, Val(0), (v,)) for v in find_leaves!(dset, vertex)
                ]),
            )
        else
            rep = NamedTuple()
        end
        meta = (; birth_simplex=birth_simplex, death_simplex=edge, rep...)
        return PersistenceInterval(birth_time, death_time, meta)
    else
        return nothing
    end
end

"""
    zeroth_intervals(filtration, cutoff, verbose, field_type, reps)

Compute 0-dimensional persistent homology using Kruskal's Algorithm.

Only keep intervals with desired birth/death `cutoff`. Compute homology with coefficients in
`field_type`. If `reps` is `true`, compute representative cocycles. Show a progress bar if
`verbose` is set.
"""
function zeroth_intervals(filtration, cutoff, verbose, ::Type{F}, reps) where {F}
    V = simplex_type(filtration, 0)
    CE = chain_element_type(V, F)
    dset = DisjointSetsWithBirth(vertices(filtration), births(filtration))

    intervals = PersistenceInterval[]

    to_skip = simplex_type(filtration, 1)[]
    to_reduce = simplex_type(filtration, 1)[]
    simplices = sort!(edges(filtration))
    if verbose
        progbar = Progress(
            length(simplices) + nv(filtration); desc="Computing 0d intervals... "
        )
    end
    for edge in simplices
        u, v = vertices(edge)
        i = find_root!(dset, u)
        j = find_root!(dset, v)
        if i ≠ j
            # According to the elder rule, the vertex with the higer birth will die first.
            last_vertex = birth(dset, i) > birth(dset, j) ? i : j
            int = interval(dset, filtration, last_vertex, edge, cutoff, reps)
            !isnothing(int) && push!(intervals, int)

            union!(dset, i, j)
            push!(to_skip, edge)
        else
            push!(to_reduce, edge)
        end
        verbose && next!(progbar; showvalues=((:n_intervals, length(intervals)),))
    end
    for v in vertices(filtration)
        if find_root!(dset, v) == v && !isnothing(simplex(filtration, Val(0), (v,)))
            int = interval(dset, filtration, v, nothing, cutoff, reps)
            push!(intervals, int)
        end
        verbose && next!(progbar; showvalues=((:n_intervals, length(intervals)),))
    end
    reverse!(to_reduce)

    thresh = Float64(threshold(filtration))
    diagram = PersistenceDiagram(
        sort!(intervals; by=persistence);
        threshold=thresh,
        dim=0,
        field=F,
        filtration=filtration,
    )
    return (postprocess_diagram(filtration, diagram), to_reduce, to_skip)
end


text:
using Mosek, Base.Test

function test_lo1()
    printstream(msg::String) = nothing # print(msg)

    bkc = [MSK_BK_FX, MSK_BK_LO, MSK_BK_UP]
    blc = [30.0, 15.0, -Inf]
    buc = [30.0, +Inf, 25.0]
    bkx = [ MSK_BK_LO, MSK_BK_RA, MSK_BK_LO, MSK_BK_LO ]
    blx = [   0.0,  0.0,    0.0,    0.0]
    bux = [+Inf, 10.0, +Inf, +Inf]
    numvar = length(bkx)
    numcon = length(bkc)
    c = [ 3.0, 1.0, 5.0, 1.0 ]
    A = sparse([1, 2, 1, 2, 3, 1, 2, 2, 3],
               [1, 1, 2, 2, 2, 3, 3, 4, 4],
               [3.0, 2.0, 1.0, 1.0, 2.0, 2.0, 3.0, 1.0, 3.0 ],
               numcon,numvar)

    task = maketask()
    putstreamfunc(task,MSK_STREAM_LOG,printstream)
    putobjname(task,"lo1")
    appendcons(task,numcon)
    for i=1:numcon
        putconname(task,i,@sprintf("c%02d",i))
    end
    appendvars(task,numvar)
    for j=1:numvar
        putvarname(task,j,@sprintf("x%02d",j))
    end
    putclist(task,[1,2,3,4], c)

    putacolslice(task,
                 1, numvar+1,
                 A.colptr[1:numvar], A.colptr[2:numvar+1],
                 A.rowval,
                 A.nzval)
    putvarboundslice(task, 1, numvar+1, bkx,blx,bux)
    putconboundslice(task,1,numcon+1,bkc,blc,buc)
    putobjsense(task,MSK_OBJECTIVE_SENSE_MAXIMIZE)
    optimize(task)
    solutionsummary(task,MSK_STREAM_MSG)

    solsta = getsolsta(task,MSK_SOL_BAS)
    prosta = getprosta(task,MSK_SOL_BAS)

    @test getsolsta(task,MSK_SOL_ITR) in [ MSK_SOL_STA_OPTIMAL, MSK_SOL_STA_NEAR_OPTIMAL ]
    @test solsta in [ MSK_SOL_STA_OPTIMAL, MSK_SOL_STA_NEAR_OPTIMAL ]
    @test prosta in [MSK_PRO_STA_PRIM_AND_DUAL_FEAS,MSK_PRO_STA_NEAR_PRIM_AND_DUAL_FEAS]

    if solsta in     [ MSK_SOL_STA_OPTIMAL,
                       MSK_SOL_STA_NEAR_OPTIMAL ]
        xx = getxx(task,MSK_SOL_BAS)
        # check feasibility and optimality of solution
    end
end

function test_qo1()
    printstream(msg::String) = nothing # print(msg)
    bkc   = [ MSK_BK_LO ]
    blc   = [ 1.0 ]
    buc   = [ Inf ]
    bkx   = [ MSK_BK_LO, MSK_BK_LO, MSK_BK_LO ]
    blx   = [ 0.0,  0.0, 0.0 ]
    bux   = [ Inf,  Inf, Inf ]
    numvar = length(bkx)
    numcon = length(bkc)
    c     = [ 0.0, -1.0, 0.0 ]
    A     = sparse( [ 1, 1, 1 ],
                    [ 1, 2, 3 ],
                    [ 1.0, 1.0, 1.0 ],
                    numcon, numvar )
    task = maketask()
    putstreamfunc(task,MSK_STREAM_LOG,printstream)
    appendcons(task,numcon)
    appendvars(task,numvar)
    putclist(task,[1:numvar;],c)
    putvarboundslice(task,1,numvar+1,bkx,blx,bux)
    putacolslice(task,1,numvar+1,
                 A.colptr[1:numvar],A.colptr[2:numvar+1],
                 A.rowval,A.nzval)
    qsubi = [ 1,   2,    3,   3   ]
    qsubj = [ 1,   2,    1,   3   ]
    qval  = [ 2.0, 0.2, -1.0, 2.0 ]
    putqobj(task,qsubi,qsubj,qval)
    putobjsense(task,MSK_OBJECTIVE_SENSE_MINIMIZE)
    optimize(task)
    solutionsummary(task,MSK_STREAM_MSG)

    prosta = getprosta(task,MSK_SOL_ITR)
    solsta = getsolsta(task,MSK_SOL_ITR)

    @test solsta in (MSK_SOL_STA_OPTIMAL, MSK_SOL_STA_NEAR_OPTIMAL)
    @test prosta in (MSK_PRO_STA_PRIM_AND_DUAL_FEAS,MSK_PRO_STA_NEAR_PRIM_AND_DUAL_FEAS)

    if solsta in     [ MSK_SOL_STA_OPTIMAL,
                       MSK_SOL_STA_NEAR_OPTIMAL ]
        xx = getxx(task,MSK_SOL_ITR)
        # check feasibility and optimality of solution
    end
end

function test_qcqo1()
    printstream(msg::String) = nothing # print(msg)
    task = maketask()
    putstreamfunc(task,MSK_STREAM_LOG,m -> nothing)
    bkc   = [ MSK_BK_LO ]
    blc   = [ 1.0 ]
    buc   = [ Inf ]
    bkx   = [ MSK_BK_LO, MSK_BK_LO, MSK_BK_LO ]
    blx   = [ 0.0,  0.0, 0.0 ]
    bux   = [ Inf,  Inf, Inf ]
    c     = [ 0.0, -1.0, 0.0 ]
    asub  = [ 1 ,2, 3 ]
    aval  = [ 1.0, 1.0, 1.0 ]
    numvar = length(bkx)
    numcon = length(bkc)
    appendcons(task,numcon)
    appendvars(task,numvar)
    putcfix(task,0.0)
    putclist(task,[1:numvar;],c)
    putvarboundslice(task,1,numvar+1,bkx,blx,bux)
    putarow(task,1,asub,aval)
    putconbound(task,1,bkc[1],blc[1],buc[1])
    qsubi = [ 1,   2,    3,   3   ]
    qsubj = [ 1,   2,    1,   3   ]
    qval  = [ 2.0, 0.2, -1.0, 2.0 ]
    putqobj(task,qsubi,qsubj,qval)
    qsubi = [  1,    2,    3,   3   ]
    qsubj = [  1,    2,    3,   1   ]
    qval  = [ -2.0, -2.0, -0.2, 0.2 ]
    putqconk(task,1, qsubi,qsubj, qval)
    putobjsense(task,MSK_OBJECTIVE_SENSE_MINIMIZE)
    optimize(task)
    solutionsummary(task,MSK_STREAM_MSG)

    prosta = getprosta(task,MSK_SOL_ITR)
    solsta = getsolsta(task,MSK_SOL_ITR)

    @test solsta in (MSK_SOL_STA_OPTIMAL, MSK_SOL_STA_NEAR_OPTIMAL)
    @test prosta in (MSK_PRO_STA_PRIM_AND_DUAL_FEAS,MSK_PRO_STA_NEAR_PRIM_AND_DUAL_FEAS)

    if solsta in     [ MSK_SOL_STA_OPTIMAL,
                       MSK_SOL_STA_NEAR_OPTIMAL ]
        xx = getxx(task,MSK_SOL_ITR)
        # check feasibility and optimality of solution
    end
end

function test_milo1()
    printstream(msg::String) = print(msg)  #print(msg)
    task = maketask()
    putstreamfunc(task,MSK_STREAM_LOG,printstream)
    bkc = [ MSK_BK_UP, MSK_BK_LO  ]
    blc = [      -Inf,      -4.0  ]
    buc = [     250.0,       Inf  ]
    bkx = [ MSK_BK_LO, MSK_BK_LO  ]
    blx = [       0.0,       0.0  ]
    bux = [       Inf,       Inf  ]
    c   = [       1.0,      0.64 ]
    A    = sparse( [ 1, 1, 2, 2],
                   [ 1, 2, 1, 2],
                   [ 50.0, 3.0, 31.0, -2.0])
    numvar = length(bkx)
    numcon = length(bkc)
    appendcons(task,numcon)
    appendvars(task,numvar)
    putclist(task,[1:numvar;],c)
    putvarboundslice(task,1,numvar+1,bkx,blx,bux)
    putacolslice(task,1,numvar+1, A.colptr[1:numvar],A.colptr[2:numvar+1],A.rowval,A.nzval)
    putconboundslice(task,1,numcon+1,bkc,blc,buc)
    putobjsense(task,MSK_OBJECTIVE_SENSE_MAXIMIZE)
    putvartypelist(task,[ 1, 2 ], [ MSK_VAR_TYPE_INT, MSK_VAR_TYPE_INT ])
    optimize(task)
    solutionsummary(task,MSK_STREAM_MSG)
    prosta = getprosta(task,MSK_SOL_ITG)
    solsta = getsolsta(task,MSK_SOL_ITG)

    @test solsta in (MSK_SOL_STA_INTEGER_OPTIMAL, MSK_SOL_STA_NEAR_INTEGER_OPTIMAL)
    @test prosta in (MSK_PRO_STA_PRIM_FEAS,MSK_PRO_STA_NEAR_PRIM_FEAS)

    if solsta in     [ MSK_SOL_STA_INTEGER_OPTIMAL,
                       MSK_SOL_STA_NEAR_INTEGER_OPTIMAL ]
        xx = getxx(task,MSK_SOL_ITG)
    end
end

function test_cqo()
        printstream(msg::String) = nothing # print(msg)
    callback(where,dinf,iinf,liinf) = 0
    bkc = [ MSK_BK_FX ]
    blc = [ 1.0 ]
    buc = [ 1.0 ]
    c   = [               0.0,              0.0,              0.0,
                          1.0,              1.0,              1.0 ]
    bkx = [ MSK_BK_LO,MSK_BK_LO,MSK_BK_LO,
            MSK_BK_FR,MSK_BK_FR,MSK_BK_FR ]
    blx = [               0.0,              0.0,              0.0,
                          -Inf,             -Inf,             -Inf ]
    bux = [               Inf,              Inf,              Inf,
                          Inf,              Inf,              Inf ]
    asub  = [ 1 ,2, 3 ]
    aval  = [ 1.0, 1.0, 1.0 ]
    numvar = length(bkx)
    numcon = length(bkc)
    task = maketask()
    putstreamfunc(task,MSK_STREAM_LOG,printstream)
    putcallbackfunc(task,callback)
    appendcons(task,numcon)
    appendvars(task,numvar)
    putclist(task,[1:6;],c)
    putvarboundslice(task,1,numvar+1,bkx,blx,bux)
    putarow(task,1,asub,aval)
    putconbound(task,1,bkc[1],blc[1],buc[1])
    appendcone(task,MSK_CT_QUAD, 0.0, [ 4, 1, 2 ])
    appendcone(task,MSK_CT_RQUAD, 0.0, [ 5, 6, 3 ])
    putobjsense(task,MSK_OBJECTIVE_SENSE_MINIMIZE)
    optimize(task)
    solutionsummary(task,MSK_STREAM_MSG)
    prosta = getprosta(task,MSK_SOL_ITR)
    solsta = getsolsta(task,MSK_SOL_ITR)

    @test solsta in (MSK_SOL_STA_OPTIMAL, MSK_SOL_STA_NEAR_OPTIMAL)
    @test prosta in (MSK_PRO_STA_PRIM_AND_DUAL_FEAS,MSK_PRO_STA_NEAR_PRIM_AND_DUAL_FEAS)

    if solsta in     [ MSK_SOL_STA_OPTIMAL,
                       MSK_SOL_STA_NEAR_OPTIMAL ]
        xx = getxx(task,MSK_SOL_ITR)
        # check feasibility and optimality of solution
    end
end

function test_sdo1()
    task = maketask()
    putstreamfunc(task,MSK_STREAM_LOG,m -> nothing)
    bkc = [MSK_BK_FX,
           MSK_BK_FX]
    blc = [1.0, 0.5]
    buc = [1.0, 0.5]
    A = sparse( [1,2,2],[1,2,3],[1.0, 1.0, 1.0])
    conesub = [1, 2, 3]
    barci = [1, 2, 2, 3, 3]
    barcj = [1, 1, 2, 2, 3]
    barcval = [2.0, 1.0, 2.0, 1.0, 2.0]
    barai   = Any[ [1, 2, 3],
                   [1, 2, 3, 2, 3, 3] ]
    baraj   = Any[ [1, 2, 3],
                   [1, 1, 1, 2, 2, 3] ]
    baraval = Any[ [1.0, 1.0, 1.0],
                   [1.0, 1.0, 1.0, 1.0, 1.0, 1.0] ]
    numvar = 3
    numcon = length(bkc)
    barvardim = [3]
    appendvars(task,numvar)
    appendcons(task,numcon)
    appendbarvars(task,barvardim)
    putcj(task, 1, 1.0)
    putvarboundslice(task,1,numvar+1,
                     [ MSK_BK_FR for i in 1:numvar ],
                     [ -Inf      for i in 1:numvar ],
                     [ +Inf      for i in 1:numvar ])
    putconboundslice(task,1,numcon+1, bkc,blc,buc)
    putacolslice(task,1,numvar+1,
                 A.colptr[1:numvar], A.colptr[2:numvar+1],
                 A.rowval,A.nzval)
    appendcone(task,MSK_CT_QUAD, 0.0, conesub)
    symc  = appendsparsesymmat(task,barvardim[1],
                               barci,
                               barcj,
                               barcval)

    syma0 = appendsparsesymmat(task,
                               barvardim[1],
                               barai[1],
                               baraj[1],+
                               baraval[1])
    syma1 = appendsparsesymmat(task,barvardim[1],
                               barai[2],
                               baraj[2],
                               baraval[2])
    putbarcj(task,1, [symc], [1.0])
    putbaraij(task,1, 1, [syma0], [1.0])
    putbaraij(task,2, 1, [syma1], [1.0])
    putobjsense(task,MSK_OBJECTIVE_SENSE_MINIMIZE)
    optimize(task)
    solutionsummary(task,MSK_STREAM_MSG)
    prosta = getprosta(task,MSK_SOL_ITR)
    solsta = getsolsta(task,MSK_SOL_ITR)

    @test solsta in (MSK_SOL_STA_OPTIMAL, MSK_SOL_STA_NEAR_OPTIMAL)
    @test prosta in (MSK_PRO_STA_PRIM_AND_DUAL_FEAS,MSK_PRO_STA_NEAR_PRIM_AND_DUAL_FEAS)

    if solsta in     [ MSK_SOL_STA_OPTIMAL,
                       MSK_SOL_STA_NEAR_OPTIMAL ]
        xx = getxx(task,MSK_SOL_ITR)
        # check feasibility and optimality of solution
    end
end

function test_nlo1()
    t = maketask()
    putstreamfunc(t,MSK_STREAM_LOG,m -> nothing)
    appendvars(t,3)
    appendcons(t,1)
    putvarbound(t,1,MSK_BK_LO,0.0, +Inf)  # x0
    putvarbound(t,2,MSK_BK_LO,0.0, +Inf)  # x1
    putvarbound(t,3,MSK_BK_LO,0.0, +Inf)  # x2
    putarow(t,1,[1,2,3],[1.0,1.0,1.0])   # x0 + x1 + x2
    putconbound(t,1,MSK_BK_FX,1.0, 1.0)  # = 1.0
    putclist(t,[1], [-1.0])

    function indexof(v,a)
        for i in 1:length(a)
            if a[i] == v return i
            end
        end
        return 0
    end
    evalobj(x :: Vector{Float64}) = log(x[2]+x[3])
    evalconi(x:: Vector{Float64},i:: Int32) = 0.0
    function grdobj(x :: Vector{Float64},sub:: Vector{Int32}, val:: Vector{Float64})
        for j in 1:length(sub)
            if sub[j] == 2 || sub[j] == 3
                val[j] = -1.0 / (x[2]+x[3])
            end
        end
    end
    function grdconi(x  :: Vector{Float64},
                     i  :: Int32,
                     sub:: Vector{Int32},
                     val:: Vector{Float64})
        none
    end
    function grdlag(x ::   Vector{Float64},
                    yo::   Float64,
                    yc::   Vector{Float64},
                    subi:: Vector{Int32},
                    val::  Vector{Float64})
        val[2] = - yo * 1.0 / (x[2]+x[3])
        val[3] = - yo * 1.0 / (x[2]+x[3])
    end
    function heslag(x ::      Vector{Float64},
                    yo::      Float64,
                    yc::      Vector{Float64},
                    subi::    Vector{Int32},
                    hessubi:: Vector{Int32},
                    hessubj:: Vector{Int32},
                    hesval::  Vector{Float64})

        hessubi[1] = 2; hessubj[1] = 2; hesval[1] = (x[2]+x[3])^(-2)
        hessubi[2] = 3; hessubj[2] = 2; hesval[2] = (x[2]+x[3])^(-2)
        hessubi[3] = 3; hessubj[3] = 3; hesval[3] = (x[2]+x[3])^(-2)
    end

    putnlcallbacks(t,
                   [2,3], # subscripts of non-zeros in the gradient of the objective
                   Int[], # subscripts of non-zeros in the gradient of the constraints
                   [1,1], # rowptr for subscripts of non-zeros in the gradient of the constraints
                   [2,3,3], # hessubi
                   [2,2,3], # hessubj
                   evalobj,
                   evalconi,
                   grdlag,
                   grdobj,
                   grdconi,
                   heslag)

    optimize(t)

    solsta = getsolsta(t,MSK_SOL_ITR)
    prosta = getprosta(t,MSK_SOL_ITR)

    @test solsta in (MSK_SOL_STA_OPTIMAL, MSK_SOL_STA_NEAR_OPTIMAL)
    @test prosta in (MSK_PRO_STA_PRIM_AND_DUAL_FEAS,MSK_PRO_STA_NEAR_PRIM_AND_DUAL_FEAS)
end

@testset "[apitest]" begin
    @testset "lo1" begin
        test_lo1()
    end

    @testset "qo1" begin
        test_qo1()
    end

    @testset "qcqo1" begin
        test_qcqo1()
    end

    @testset "milo1" begin
        test_milo1()
    end

    @testset "cqo" begin
        test_cqo()
    end

    @testset "sdo1" begin
        test_sdo1()
    end

    @testset "nlo1" begin
        test_nlo1()
    end
end


text:
using Base.Test
# blob.jl
import WaveletScattering: ScatteredBlob
# node.jl
import WaveletScattering: Node
# path.jl
import WaveletScattering: Path, PathKey

# show
io = IOBuffer()
node = Node([1.0, 2.0], (PathKey(:time)=>1:1:2,))
blob = ScatteredBlob(DataStructures.SortedDict((Path() => node,)))
show(io, blob)
@test String(take!(io)) == "ScatteredBlob(1 node)"
blob = ScatteredBlob(DataStructures.SortedDict(
    (Path((:γ, :time) => 0) => node, Path((:γ, :time) => 1) => node)))
show(io, blob)
@test String(take!(io)) == "ScatteredBlob(2 nodes)"


text:
using PlantBiophysics
using Documenter
using PlantBiophysics
using DataFrames
using CSV
using Plots

DocMeta.setdocmeta!(PlantBiophysics, :DocTestSetup, :(using PlantBiophysics,Plots,DataFrames,CSV); recursive=true)

makedocs(;
    modules=[PlantBiophysics],
    authors="remi.vezy <VEZY@users.noreply.github.com> and contributors",
    repo="https://github.com/VEZY/PlantBiophysics.jl/blob/{commit}{path}#L{line}",
    sitename="PlantBiophysics.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://VEZY.github.io/PlantBiophysics.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Getting started" => "get_started.md",
        "Models" => Any[
            "Photosynthesis" => "photosynthesis.md",
            "Stomatal conductance" => "gs.md",
            "Energy balance" => "energy_balance.md",
            "Light interception" => "light.md"
            ],
        "Details" => Any[
            "Package design" => "concepts.md",
            "Implement a model" => "implement_a_model.md",
            "Uncertainty propagation" => "uncertainty_propagation.md",
        ],
        "Functions" => "functions.md"
    ],
)

deploydocs(;
    repo="github.com/VEZY/PlantBiophysics.jl.git",
)


text:
function disparate_impact(ft::Fairness.FairTensor; priv_grps::Array{String, 1})
	num, den = 0, 0
	priv_count, unpriv_count = 0, 0
	for grp in ft.labels
		if grp in priv_grps
			den += sum(ft[Fairness._ftIdx(ft, grp), 2, :])
			priv_count += sum(ft[Fairness._ftIdx(ft, grp), :, :])
		else
			num += sum(ft[Fairness._ftIdx(ft, grp), 1, :])
			unpriv_count += sum(ft[Fairness._ftIdx(ft, grp), :, :])
		end
	end
	return (num/unpriv_count)/(den/priv_count)
end

struct DisparateImpact <: Fairness.Measure
	grp::Symbol
	priv_grps::Array{String, 1}
end

function (di::DisparateImpact)(ŷ, X, y)
	if typeof(ŷ[1]) <: MLJBase.UnivariateFinite
		ŷ = StatsBase.mode.(ŷ)
	end
	grps = X[:, di.grp]
	ft = fair_tensor(categorical(ŷ), categorical(y), categorical(grps))
	return disparate_impact(ft, priv_grps=di.priv_grps)
end

MLJBase.name(::Type{<:DisparateImpact}) = "DisparateImpact"
# MLJBase.target_scitype(::Type{<:DisparateImpact}) = AbstractArray{Multiclass{2},1}
MLJBase.supports_weights(::Type{<:DisparateImpact}) = false # for now
# MLJBase.prediction_type(::Type{<:DisparateImpact}) = :deterministic # Not specifying it to have check_measures false
MLJBase.orientation(::Type{<:DisparateImpact}) = :other # other options are :score, :loss
MLJBase.reports_each_observation(::Type{<:DisparateImpact}) = false
MLJBase.aggregation(::Type{<:DisparateImpact}) = Fairness.Mean()
MLJBase.is_feature_dependent(::Type{<:DisparateImpact}) = true


text:
using LispSyntax
using Stage

# ----------------------------------------------------------------------------------------------------------------------
# Setup
# ----------------------------------------------------------------------------------------------------------------------
macro incr(x)
  quote
    $(esc(x)) = $(esc(x)) + 1
    $(esc(x))
  end
end

# ----------------------------------------------------------------------------------------------------------------------
# Reader
# ----------------------------------------------------------------------------------------------------------------------
@expect LispSyntax.read("1.1f") == 1.1f0
@expect LispSyntax.read("1.2f") == 1.2f0
@expect LispSyntax.read("2f")   == 2f0

@expect LispSyntax.read("3.0d") == 3.0

@expect LispSyntax.read("4")    == 4

@expect LispSyntax.read("\\u2312") == '\u2312'

@expect LispSyntax.read("\\040") == ' '

@expect LispSyntax.read("\\c") == 'c'

@expect LispSyntax.read("\"test\"") == "test"

@expect LispSyntax.read("true") == true
@expect LispSyntax.read("false") == false

@expect LispSyntax.read("test") == :test

@expect LispSyntax.read("()") == sx()
@expect LispSyntax.read("(1.1f)") == sx(1.1f0)
@expect LispSyntax.read("(1.1f 2.2f)") == sx(1.1f0, 2.2f0)
@expect LispSyntax.read("(+ 1.1f 2)") == sx(:+, 1.1f0, 2)
@expect LispSyntax.read("(this (+ 1.1f 2))") == sx(:this, sx(:+, 1.1f0, 2))
@expect LispSyntax.read("(this (+ 1.1f 2) )") == sx(:this, sx(:+, 1.1f0, 2))

@expect LispSyntax.read("#{1 2 3 4}") == Set([1, 2, 3, 4])

@expect LispSyntax.read("""#{
                        1 2
                        3 4
                        }""") == Set([1, 2, 3, 4])

@expect LispSyntax.read("{a 2 b 3}") == Dict(:a => 2, :b => 3)

@expect LispSyntax.read("""{
                        a 2
                        b 3
                        }""") == Dict(:a => 2, :b => 3)

@expect LispSyntax.read("[1 2 3 4]")  == sx(1, 2, 3, 4)

@expect LispSyntax.read("""[
                        1 2
                        3 4
                        ]""")  == sx(1, 2, 3, 4)

@expect LispSyntax.read("[]")         == sx()
@expect LispSyntax.read("[1]")        == sx(1)

@expect LispSyntax.read("'test")      == sx(:quote, :test)

@expect LispSyntax.read("`test")      == sx(:quasi, :test)

@expect LispSyntax.read("~test")      == sx(:splice, :test)
@expect LispSyntax.read("~@(1 2 3)")  == sx(:splice_seq, sx(1, 2, 3))

@expect LispSyntax.read("`~test")     == sx(:quasi, sx(:splice, :test))

@expect desx(sx(:splice_seq, sx(1, 2, 3))) == Any[ :splice_seq, [1, 2, 3] ]
@expect desx(sx(:splice_seq, sx(1, 2, sx(3)))) == Any[ :splice_seq, Any[ 1, 2, [3] ] ]

@expect LispSyntax.read("""(defn multiline
                           [x]
                           (+ x 1))""") == sx(:defn, :multiline, sx(:x), sx(:+, :x, 1))

@expect LispSyntax.read("""
(defn f1 [n]
   (if (< n 2)
       1
       (+ (f1 (- n 1))
          (f1 (- n 2)))))
""") == sx(:defn, :f1, sx(:n), sx(:if, sx(:<, :n, 2), 1, sx(:+, sx(:f1, sx(:-, :n, 1)), sx(:f1, sx(:-, :n, 2)))))

assign_reader_dispatch(:sx, x -> sx(x.vector...))
assign_reader_dispatch(:hash, x -> [ x.vector[i] => x.vector[i+1] for i = 1:2:length(x.vector) ])
@expect LispSyntax.read("#sx[a b c]") == sx(:a, :b, :c)
@expect LispSyntax.read("#sx [ 1 2 3 ]") == sx(1, 2, 3)

# ----------------------------------------------------------------------------------------------------------------------
# Code generation
# ----------------------------------------------------------------------------------------------------------------------
@expect codegen(desx(LispSyntax.read("(if true a)"))) == :(true && $(esc(:a)))
@expect codegen(desx(LispSyntax.read("(if true a b)"))) == :(true ? $(esc(:a)) : $(esc(:b)))

@expect codegen(desx(LispSyntax.read("(call)"))) == :($(esc(:call))())
@expect codegen(desx(LispSyntax.read("(call a)"))) == :($(esc(:call))($(esc(:a))))
@expect codegen(desx(LispSyntax.read("(call a b)"))) == :($(esc(:call))($(esc(:a)), $(esc(:b))))
@expect codegen(desx(LispSyntax.read("(call a b c)"))) == :($(esc(:call))($(esc(:a)), $(esc(:b)), $(esc(:c))))

@expect codegen(desx(LispSyntax.read("(lambda (x) (call x))"))) == Expr(:function, :((x,)), Expr(:block, :($(esc(:call))(x))))
@expect codegen(desx(LispSyntax.read("(def x 3)"))) == :($(esc(:x)) = 3)
@expect codegen(desx(LispSyntax.read("(def x (+ 3 1))"))) == :($(esc(:x)) = $(esc(:+))(3, 1))

@expect codegen(desx(LispSyntax.read("test"))) == :($(esc(:test)))
@expect codegen(desx(LispSyntax.read("'test"))) == QuoteNode(:test)
@expect codegen(desx(LispSyntax.read("'(1 2)"))) == :(construct_sexpr(1, 2))
@expect codegen(desx(LispSyntax.read("'(1 x)"))) == :(construct_sexpr(1, :x))
@expect codegen(desx(LispSyntax.read("'(1 (1 2))"))) == :(construct_sexpr(1, construct_sexpr(1, 2)))
@expect codegen(desx(LispSyntax.read("'(1 (test x))"))) == :(construct_sexpr(1, construct_sexpr(:test, :x)))
@expect codegen(desx(LispSyntax.read("(call 1 '2)"))) == :($(esc(:call))(1, 2))

# ----------------------------------------------------------------------------------------------------------------------
# Scope and variables
# ----------------------------------------------------------------------------------------------------------------------
global x = 10
@expect @lisp("x") == 10
@expect lisp"x" == 10

lisp"(def w (+ 3 1))"
@expect w == 4

# ----------------------------------------------------------------------------------------------------------------------
# Quoting and splicing
# ----------------------------------------------------------------------------------------------------------------------
@expect @lisp("`~x") == 10
@expect lisp"`~x" == 10
@expect lisp"'test" == :test
@expect lisp"'(1 2)" == Any[1, 2]
@expect lisp"'(1 x)" == Any[1, :x]
@expect lisp"'(1 (1 2))" == Any[1, Any[1, 2]]
@expect lisp"'(1 (test x))" == Any[1, Any[:test, :x]]
@expect @lisp("`(test ~x)") == Any[ :test, 10 ]
@expect lisp"`(test ~x)" == Any[ :test, 10 ]
@expect @lisp("`(~x ~x)") == Any[ 10, 10 ]
global y = Any[ 1, 2 ]
@expect @lisp("`(~x ~@y)") == Any[ 10, 1, 2 ]
@expect @lisp("`(~x ~y)") == Any[ 10, Any[1, 2] ]

@expect @lisp("`(10 ~(+ 10 x))") == Any[10, 20]

@expect lisp"(quote (+ 1 2))" == Any[:+, 1, 2]

# ----------------------------------------------------------------------------------------------------------------------
# Functions
# ----------------------------------------------------------------------------------------------------------------------
@lisp("(defn xxx [a b] (+ a b))")
@expect @lisp("(xxx 1 2)") == 3

global z = 10
@lisp("(defn yyy [a] (+ a z))")
@expect @lisp("(yyy 1)") == 11
@expect @lisp("(yyy z)") == 20

# recursion
lisp"(defn fib [a] (if (< a 2) a (+ (fib (- a 1)) (fib (- a 2)))))"
@expect lisp"(fib 2)" == 1
@expect lisp"(fib 4)" == 3
@expect lisp"(fib 30)" == 832040
@expect lisp"(fib 40)" == 102334155

# Note this version is very slow due to the anonymous function
lisp"(def fib2 (lambda [a] (if (< a 2) a (+ (fib2 (- a 1)) (fib2 (- a 2))))))"
@expect lisp"(fib2 2)" == 1
@expect lisp"(fib2 4)" == 3
@expect lisp"(fib2 30)" == 832040

lisp"(defn dostuff [a] (@incr a) (@incr a) (@incr a))"
@expect lisp"(dostuff 3)" == 6
@expect lisp"(dostuff 6)" == 9

lisp"(def dostuff2 (lambda [a] (@incr a) (@incr a) (@incr a)))"
@expect lisp"(dostuff2 3)" == 6
@expect lisp"(dostuff2 6)" == 9

lisp"(def dostuff3 (fn [a] (@incr a) (@incr a) (@incr a)))"
@expect lisp"(dostuff3 3)" == 6
@expect lisp"(dostuff3 6)" == 9
@expect lisp"((lambda [x] (+ x 1)) 5)" == 6
@expect lisp"#{1 2 z}" == Set([1, 2, 10])
@expect lisp"{1 2 2 z}" == Dict(1 => 2, 2 => 10)
@expect lisp"#sx[+ 1 2]" == 3
@expect lisp"#hash['+ 1 '- z]" == Dict(:+ => 1, :- => 10)

# ----------------------------------------------------------------------------------------------------------------------
# Macros
# ----------------------------------------------------------------------------------------------------------------------
lisp"(defn fact [a] (if (< a 1) 1 (* a (fact (- a 1)))))"
lisp"(defmacro fapply [f a] `(~f ~a))"
@expect @fapply(fib2, 2) == 1
@expect @fapply(fact, 3 + 1) == 24
@expect lisp"(@fapply fib2 2)" == 1
@expect lisp"(@fapply fact (+ 3 1))" == 24

fcount = 0
lisp"(defmacro fapply_trace [f a] (global fcount) (@incr fcount) `(~f ~a))"
@expect @fapply_trace(fib2, 2) == 1
@expect fcount == 1
@expect @fapply_trace(fact, 3 + 1) == 24
@expect fcount == 2

# ----------------------------------------------------------------------------------------------------------------------
# Loops
# ----------------------------------------------------------------------------------------------------------------------
number = 0
output = 0

lisp"(while (< number 2) (@incr number) (@incr output))"
@expect number == 2
@expect output == 2
r = output
lisp"(for [i (range 1 10)] (@incr r))"
@expect r == 12

r = 0
lisp"(for [i (range 1 10) j (range 1 10)] (@incr r))"
@expect r == 100

# ----------------------------------------------------------------------------------------------------------------------
# Let and do
# ----------------------------------------------------------------------------------------------------------------------
@expect lisp"(let [x 10] x)" == 10
@expect lisp"(let [x 10 y 20] (+ x y))" == 30
@expect lisp"(let [x 10 y 20 z 20] (+ x y z))" == 50
@expect lisp"(let [x 10 y 20 z 20] (+ x y z number))" == 52
@expect lisp"(let [x 10 y 20 z 20 number 10] (+ x y z number))" == 60
@expect lisp"(let [x 10 y 20 z 20] (- (+ x y z number) output))" == 50

lisp"(do (@incr r) (@incr number))"
@expect number == 3
@expect r == 101

# ----------------------------------------------------------------------------------------------------------------------
# Import
# ----------------------------------------------------------------------------------------------------------------------
lisp"(import ParserCombinator)"
@expect lisp"(@E_str \"S\")" == E"S"

# ----------------------------------------------------------------------------------------------------------------------
# Bug reports
# ----------------------------------------------------------------------------------------------------------------------
@expect lisp"""(def game_map (Dict
          (=> 'living_room
              '((you are in the living room
                 of a wizards house - there is a wizard
                 snoring loudly on the couch -)
                (west door garden)
                (upstairs stairway attic)))))""" == Dict(:living_room =>
                                                         Any[ Any[ :you, :are, :in, :the, :living, :room, :of, :a, :wizards, :house, :-,
                                                                   :there, :is, :a, :wizard, :snoring, :loudly, :on, :the, :couch, :- ],
                                                              Any[ :west, :door, :garden ],
                                                              Any[ :upstairs, :stairway, :attic ] ])


text:
const car_width = 0.05
const car_height = car_width / 2.0
const clearance = 0.2 * car_height
const flag_height = 0.05

struct MountainCarState <: State
    position::Float64
    speed::Float64
end

function vec(state::MountainCarState)::Array{Float64}
    [state.position, state.speed]
end

struct MountainCarDiscState <: State
    disc_position::Int64
    disc_speed::Int64
end

function vec(disc_state::MountainCarDiscState)::Array{Int64}
    [disc_state.disc_position, disc_state.disc_speed]
end

function unvec(state::Array{Float64}; cont::Bool)
    if cont
        return MountainCarState(state[1], state[2])
    else
        return MountainCarDiscState(state[1], state[2])
    end
end

struct MountainCarAction <: Action
    id::Int64
end

struct MountainCarParameters <: Parameters
    theta1::Float64
    theta2::Float64
end

function vec(params::MountainCarParameters)::Array{Float64}
    [params.theta1, params.theta2]
end

function unvec_params(params::Array{Float64})::MountainCarParameters
    MountainCarParameters(params[1], params[2])
end

struct MountainCarTransition <: Transition
    initial_state::MountainCarDiscState
    action::MountainCarAction
    cost::Float64
    final_state::MountainCarDiscState
end

struct MountainCarContTransition <: Transition
    initial_state::MountainCarState
    action::MountainCarAction
    cost::Float64
    final_state::MountainCarState
end

struct MountainCar
    min_position::Float64
    max_position::Float64
    max_speed::Float64
    goal_position::Float64
    goal_speed::Float64
    force::Float64
    position_discretization::Int64
    speed_discretization::Int64
    position_grid_cell_size::Float64
    speed_grid_cell_size::Float64
    start_state::MountainCarState
    rock_c::Float64
    rock_position::Float64
    noise::Normal
end

function MountainCar(rock_c::Float64; position_sigma::Float64 = 0.0)
    start_state = MountainCarState(-π / 6, 0)
    position_discretization = 150 # 300 # 500
    speed_discretization = 150 # 500 # 250 # new param
    min_position = -1.2
    max_position = 0.5 # 0.6
    max_speed = 0.07  # 0.2
    goal_position = 0.5
    goal_speed = 0
    force = 0.001
    rock_position = 0.25
    position_grid_cell_size = (max_position - min_position) / (position_discretization - 1)
    speed_grid_cell_size = (max_speed * 2) / (speed_discretization - 1)
    noise = Normal(0.0, position_sigma)
    MountainCar(
        min_position,
        max_position,
        max_speed,
        goal_position,
        goal_speed,
        force,
        position_discretization,
        speed_discretization,
        position_grid_cell_size,
        speed_grid_cell_size,
        start_state,
        rock_c,
        rock_position,
        noise,
    )
end

function init(mountaincar::MountainCar; rng = nothing, cont = false)
    start_state = mountaincar.start_state
    if !isnothing(rng)
        random_position = rand(rng, Uniform(mountaincar.min_position, mountaincar.max_position))
        random_speed = rand(rng, Uniform(-mountaincar.max_speed, mountaincar.max_speed))
        start_state = MountainCarState(random_position, random_speed)
    end
    if cont
        return start_state
    end
    cont_state_to_disc(mountaincar, start_state)
end

function init_random(mountaincar::MountainCar, rng::MersenneTwister)
    start_state = mountaincar.start_state
    noise = rand(rng, Uniform(-0.1, 0.1))
    MountainCarState(start_state.position + noise, start_state.speed)
end

function getActions(mountaincar::MountainCar)
    [MountainCarAction(0), MountainCarAction(1)]
end

function checkGoal(mountaincar::MountainCar, state::MountainCarState)
    state.position >= mountaincar.goal_position
end

function checkGoal(mountaincar::MountainCar, state::MountainCarDiscState)
    checkGoal(mountaincar, disc_state_to_cont(mountaincar, state))
end

function step(
    mountaincar::MountainCar,
    disc_state::MountainCarDiscState,
    action::MountainCarAction,
    params::Array{Float64};
    debug = false,
    rng = nothing,
)
    step(
        mountaincar,
        disc_state,
        action,
        MountainCarParameters(params[1], params[2]),
        debug = debug,
        rng = rng,
    )
end

function step(
    mountaincar::MountainCar,
    state::MountainCarState,
    action::MountainCarAction,
    params::Array{Float64};
    debug = false,
    rng = nothing,
)
    step(
        mountaincar,
        state,
        action,
        MountainCarParameters(params[1], params[2]),
        debug = debug,
        rng = rng,
    )
end

function step(
    mountaincar::MountainCar,
    disc_state::MountainCarDiscState,
    action::MountainCarAction,
    params::MountainCarParameters;
    debug = false,
    rng = nothing,
)
    state = disc_state_to_cont(mountaincar, disc_state)
    next_slipped_state, cost = step(
        mountaincar,
        state,
        action,
        params;
        debug = debug,
        rng = rng,
    )
    next_slipped_disc_state = cont_state_to_disc(mountaincar, next_slipped_state)
    next_slipped_disc_state, cost
end

function step_discrete(
    mountaincar::MountainCar,
    state::MountainCarState,
    action::MountainCarAction,
    params::MountainCarParameters;
    debug = false,
    rng = nothing,
)
    cost = getCost(mountaincar, state)
    new_position = position_dynamics(mountaincar, state; rng = rng)
    new_speed = speed_dynamics(mountaincar, state, action, params)
    slip = 0

    disc_state = cont_state_to_disc(mountaincar, state)
    disc_rock_position =
        cont_state_to_disc(
            mountaincar,
            MountainCarState(mountaincar.rock_position, 0),
        ).disc_position
    new_disc_position =
        cont_state_to_disc(mountaincar, MountainCarState(new_position, 0)).disc_position
    if mountaincar.rock_c > 0
        # Check for rock in discrete grid to ensure that no discretization errors happen
        if sign(disc_state.disc_position - disc_rock_position) !=
           sign(new_disc_position - disc_rock_position)
            if debug
                println("Hit rock")
            end
            if new_speed > 0
                slip = max(-mountaincar.rock_c, -new_speed)
            else
                slip = min(mountaincar.rock_c, -new_speed)
            end
        end
    end

    slipped_speed = clamp(new_speed + slip, -mountaincar.max_speed, mountaincar.max_speed)
    @assert abs(slipped_speed) <= abs(new_speed)
    MountainCarState(new_position, slipped_speed), cost
end

function step(
    mountaincar::MountainCar,
    state::MountainCarState,
    action::MountainCarAction,
    params::MountainCarParameters;
    debug = false,
    rng = nothing,
)
    cost = getCost(mountaincar, state)
    new_position = position_dynamics(mountaincar, state; rng = rng)
    new_speed = speed_dynamics(mountaincar, state, action, params)
    slip = 0

    if mountaincar.rock_c > 0
        if sign(state.position - mountaincar.rock_position) !=
           sign(new_position - mountaincar.rock_position)
            if debug
                println("Hit rock")
            end
            if new_speed > 0
                slip = max(-mountaincar.rock_c, -new_speed)
            else
                slip = min(mountaincar.rock_c, -new_speed)
            end
        end
    end

    slipped_speed = clamp(new_speed + slip, -mountaincar.max_speed, mountaincar.max_speed)
    @assert abs(slipped_speed) <= abs(new_speed)
    MountainCarState(new_position, slipped_speed), cost
end

function position_dynamics(
    mountaincar::MountainCar,
    state::MountainCarState;
    rng = nothing,
)
    new_position = state.position + state.speed
    if !isnothing(rng)
        new_position += rand(rng, mountaincar.noise)
    end
    clamp(new_position, mountaincar.min_position, mountaincar.max_position)
end

function speed_dynamics(
    mountaincar::MountainCar,
    state::MountainCarState,
    action::MountainCarAction,
    params::MountainCarParameters,
)
    new_speed =
        state.speed +
        (2 * action.id - 1) * mountaincar.force +
        cos(params.theta2 * state.position) * params.theta1
    # clamp(new_speed, -mountaincar.max_speed, mountaincar.max_speed)
    new_speed
end

function getCost(mountaincar::MountainCar, state::MountainCarState)
    if state.position < mountaincar.goal_position
        return abs(state.position - mountaincar.goal_position)^4
        # return 1.0
    else
        return 0.0
    end
end

function getCost(mountaincar::MountainCar, disc_state::MountainCarDiscState)
    getCost(mountaincar, disc_state_to_cont(mountaincar, disc_state))
end

function cont_state_to_disc(mountaincar::MountainCar, state::MountainCarState)
    position = state.position - mountaincar.min_position
    speed = state.speed + mountaincar.max_speed
    disc_position = clamp(
        cont_to_disc(position, mountaincar.position_grid_cell_size),
        0,
        mountaincar.position_discretization - 1,
    )
    disc_speed = clamp(
        cont_to_disc(speed, mountaincar.speed_grid_cell_size),
        0,
        mountaincar.speed_discretization - 1,
    )
    MountainCarDiscState(disc_position, disc_speed)
end

function absorbing_state(mountaincar::MountainCar)
    MountainCarDiscState(Inf, Inf)
end

function absorbing_state_idx(mountaincar::MountainCar)
    mountaincar.position_discretization * mountaincar.speed_discretization + 1
end

function disc_state_to_idx(mountaincar::MountainCar, disc_state::MountainCarDiscState)
    disc_position = disc_state.disc_position
    disc_speed = disc_state.disc_speed
    # Adding 1 for julia indexing
    idx = disc_position * mountaincar.speed_discretization + disc_speed + 1
    @assert idx >= 1 &&
            idx <= mountaincar.position_discretization * mountaincar.speed_discretization
    idx
end

function cont_state_to_idx(mountaincar::MountainCar, cont_state::MountainCarState)
    disc_state = cont_state_to_disc(mountaincar, cont_state)
    disc_state_to_idx(mountaincar, disc_state)
end

function idx_to_disc_state(mountaincar::MountainCar, idx::Int64)
    disc_position = (idx - 1) ÷ mountaincar.speed_discretization
    disc_speed = (idx - 1) % mountaincar.speed_discretization
    @assert disc_position >= 0 && disc_position < mountaincar.position_discretization
    @assert disc_speed >= 0 && disc_speed < mountaincar.speed_discretization
    MountainCarDiscState(disc_position, disc_speed)
end

function idx_to_cont_state(mountaincar::MountainCar, idx::Int64)
    disc_state = idx_to_disc_state(mountaincar, idx)
    disc_state_to_cont(mountaincar, disc_state)
end

function cont_to_disc(x::Float64, cell_size::Float64)
    disc_x = floor(Int64, (x + 0.5 * cell_size) / cell_size)
    if x >= 0
        return disc_x
    else
        return disc_x - 1
    end
end

function disc_state_to_cont(mountaincar::MountainCar, state::MountainCarDiscState)
    disc_position = state.disc_position
    disc_speed = state.disc_speed
    position = disc_to_cont(disc_position, mountaincar.position_grid_cell_size)
    speed = disc_to_cont(disc_speed, mountaincar.speed_grid_cell_size)
    position = clamp(
        position + mountaincar.min_position,
        mountaincar.min_position,
        mountaincar.max_position,
    )
    speed =
        clamp(speed - mountaincar.max_speed, -mountaincar.max_speed, mountaincar.max_speed)
    MountainCarState(position, speed)
end

function disc_to_cont(x::Int64, cell_size::Float64)
    x * cell_size
end