text:
# # Introduction to ClimaCore.jl
#
# ### What is ClimaCore?
#
# A suite of tools for constructing spatial discretizations
#
# - primarily aimed at climate and weather models
# - initial aim:
#   - spectral element discretization in the horizontal
#   - staggered finite difference in the vertical
# - currently under development

using ClimaCore,
    ClimaCorePlots, LinearAlgebra, IntervalSets, UnPack, Plots, OrdinaryDiffEq
#----------------------------------------------------------------------------

# ## 1. Constructing a discretization

# ### 1.1 Domains
#
# A _domain_ a region of space (think of a mathematical domain).

column_domain = ClimaCore.Domains.IntervalDomain(
    ClimaCore.Geometry.ZPoint(0.0) .. ClimaCore.Geometry.ZPoint(10.0),
    boundary_tags = (:bottom, :top),
)
#----------------------------------------------------------------------------

rectangle_domain = ClimaCore.Domains.RectangleDomain(
    ClimaCore.Geometry.XPoint(-2π) .. ClimaCore.Geometry.XPoint(2π),
    ClimaCore.Geometry.YPoint(-2π) .. ClimaCore.Geometry.YPoint(2π),
    x1periodic = true,
    x2periodic = true,
)
#----------------------------------------------------------------------------

# ### 1.2 Meshes
#
# A _mesh_ is a division of a domain into elements

column_mesh = ClimaCore.Meshes.IntervalMesh(column_domain, nelems = 32)
#----------------------------------------------------------------------------

rectangle_mesh = ClimaCore.Meshes.RectilinearMesh(rectangle_domain, 16, 16)
#----------------------------------------------------------------------------

# ### 1.3 Topologies
#
# A _topology_ determines the ordering and connections between elements of a mesh
# - At the moment, this is only required for 2D meshes

rectangle_topology = ClimaCore.Topologies.Topology2D(rectangle_mesh)
#----------------------------------------------------------------------------

# ### 1.4 Spaces
#
# A _space_ represents a discretized function space over some domain. Currently two discretizations are supported.
#

# #### 1.4.1 Staggered finite difference discretization
#
# This discretizes an interval domain by approximating the function by a value at either the center of each element (`CenterFiniteDifferenceSpace`), or the faces between elements (`FaceFiniteDifferenceSpace`).
#
# You can construct either the center or face space from the mesh, then construct the opposite space from the original one (this is to avoid allocating additional memory).

column_center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(column_mesh)
## construct the face space from the center one
column_face_space =
    ClimaCore.Spaces.FaceFiniteDifferenceSpace(column_center_space)
#----------------------------------------------------------------------------

# #### 1.4.2 Spectral element discretization
#
# A spectral element space approximates the function with polynomials in each element. The polynomials are represented using a *nodal discretization*, which stores the values of the polynomials at particular points in each element (termed *nodes* or *degrees of freedom*).
#
# These nodes are chosen by a particular *quadrature rule*, which allows us to integrate functions over the domain. The only supported choice for now is a Gauss-Legendre-Lobatto rule.

## Gauss-Legendre-Lobatto quadrature with 4 nodes in each direction, so 16 in each element
quad = ClimaCore.Spaces.Quadratures.GLL{4}()
rectangle_space =
    ClimaCore.Spaces.SpectralElementSpace2D(rectangle_topology, quad)
#----------------------------------------------------------------------------

# ### 1.5 Fields
#
# Finally, we can construct a *field*: a function in a space. A field is simply a space and the values at each node in the space.
#
# The easiest field to construct is the _coordinate field_

coord = ClimaCore.Fields.coordinate_field(rectangle_space)
#----------------------------------------------------------------------------

# This is a *struct-value field*: it contains coordinates in a struct at each point. We can extract just the `x` coordinate, to get a *scalar field*:

x = coord.x
#----------------------------------------------------------------------------

# Although you can't index directly into a field, it can be used in some other ways similar to a Julia `Array`. For example, broadcasting can be used to define new fields in terms of other ones:

sinx = sin.(x)
#----------------------------------------------------------------------------

# Fields can be easily vizualized with Plots.jl:

import Plots
Plots.plot(sinx)
#----------------------------------------------------------------------------

# If you're using the terminal, `UnicodePlots` is also supported.

# This works similarly for finite difference discretizations

column_center_coords = ClimaCore.Fields.coordinate_field(column_center_space)
column_face_coords = ClimaCore.Fields.coordinate_field(column_face_space)
#----------------------------------------------------------------------------

plot(sin.(column_center_coords.z), ylim = (0.0, 10.0))
plot!(cos.(column_face_coords.z), ylim = (0.0, 10.0))
#----------------------------------------------------------------------------

# Reduction operations are defined anologously:
#
# - `sum` will give the integral of the function
# ```math
# \int_D f(x) dx
# ```
# - `norm` will give the L² function norm
# ```math
# \sqrt{\int_D |f(x)|^2 dx}
# ```
#

sum(sinx) ## integral
#----------------------------------------------------------------------------

norm(sinx) ## L² norm
#----------------------------------------------------------------------------

# ### 1.6 Vectors and vector fields
#
# A *vector field* is a field with vector-valued quantity, i.e. at every point in space, you have a vector.
#
# However one of the key requirements of ClimaCore is to support vectors specified in curvilinear or non-Cartesian coordinates. We will discuss this in a bit further, but for now, you can define a 2-dimensional vector field using `Geometry.UVVector`:

v = ClimaCore.Geometry.UVVector.(coord.y, .-coord.x)
#----------------------------------------------------------------------------

# ## 2. Operators
#
# _Operators_ can compute spatial derivative operations.
#
#  - for performance reasons, we need to be able to "fuse" multiple operators and function applications
#  - Julia provides a tool for this: **broadcasting**, with a very flexible API
#
# Can think of operators are "pseudo-functions": can't be called directly, but act similar to functions in the context of broadcasting.

# ### 2.1 Spectral element operators
#
# The `Gradient` operator takes the gradient of a scalar field, and returns a vector field.

grad = ClimaCore.Operators.Gradient()
∇sinx = grad.(sinx)
#----------------------------------------------------------------------------

plot(∇sinx.components.data.:1, clim = (-1, 1))
#----------------------------------------------------------------------------

# This returns the gradient in [_covariant_](https://en.wikipedia.org/wiki/Covariance_and_contravariance_of_vectors) coordinates
# ```math
# (\nabla f)_i = \frac{\partial f}{\partial \xi^i}
# ```
# where $(\xi^1,\xi^2)$ are the coordinates in the *reference element*: a square $[-1,1]^2$.
#
# This can be converted to a local orthogonal basis by multiplying by the partial derivative matrix
# ```math
# \frac{\partial \xi}{\partial x}
# ```
# This can be done calling `ClimaCore.Geometry.LocalVector:

∇sinx_cart = ClimaCore.Geometry.LocalVector.(∇sinx)
#----------------------------------------------------------------------------

plot(∇sinx_cart.components.data.:1, clim = (-1, 1))
#----------------------------------------------------------------------------

plot(∇sinx_cart.components.data.:2, clim = (-1, 1))
#----------------------------------------------------------------------------

∇sinx_ref = ClimaCore.Geometry.UVVector.(cos.(x), 0.0)
norm(∇sinx_cart .- ∇sinx_ref)
#----------------------------------------------------------------------------

# Similarly, the `Divergence` operator takes the divergence of vector field, and returns a scalar field.
#
# If we take the divergence of a gradient, we can get a Laplacian:

div = ClimaCore.Operators.Divergence()
∇²sinx = div.(grad.(sinx))
plot(∇²sinx)
#----------------------------------------------------------------------------

# *Note*: In curvilinear coordinates, the divergence is defined in terms of the _contravariant_ components $u^i$:
# ```math
# \nabla \cdot u = \frac{1}{J} \sum_i \frac{\partial}{\partial \xi^i} (J u^i)
# ```
# The `Divergence` operator handles this conversion internally.

# #### 2.1.1 Direct stiffness summation
#
# Spectral element operators only operate _within_ a single element, and so the result may be discontinuous. To address this, the usual fix is _direct stiffness summation_ (DSS), which averages the values at the element boundaries.
#
# This corresponds to the $L^2$ projection onto the subset of continuous functions in our function space.

∇²sinx_dss = ClimaCore.Spaces.weighted_dss!(copy(∇²sinx))
plot(∇²sinx_dss)
#----------------------------------------------------------------------------

plot(∇²sinx_dss .- ∇²sinx)
#----------------------------------------------------------------------------

# ### 2.2 Finite difference operators
#
# Finite difference operators are similar with some subtle differences:
# - they can change staggering (center to face, or vice versa)
# - they can span multiple elements
#   - no DSS is required
#   - boundary handling may be required
#
# We use the following convention:
#  - centers are indexed by integers `1, 2, ..., n`
#  - faces are indexed by half integers `½, 1+½, ..., n+½`

# **Face to center gradient**
#
# An finite-difference operator defines a _stencil_. For example, the gradient operator
#
# ```math
# \nabla\theta[i] = \frac{\theta [i+\tfrac{1}{2}] - \theta[i-\tfrac{1}{2}]}{\Delta z}
# ```
# (actually, a little more complicated as it gives a vector in a covariant basis)
#
#
# ```
#         ...
#       /
# θ[2+½]
#       \
#         ∇θ[2]
#       /
# θ[1+½]
#       \
#         ∇θ[1]
#       /
# θ[½]
# ```
#
# Every center value is well-defined, so boundary handling is optional.
#

cosz = cos.(column_face_coords.z)
gradf2c = ClimaCore.Operators.GradientF2C()
∇cosz = gradf2c.(cosz)
#----------------------------------------------------------------------------

plot(map(x -> x.w, ClimaCore.Geometry.WVector.(∇cosz)), ylim = (0, 10))
#----------------------------------------------------------------------------

# **Center to face gradient**
#
# Uses the same stencil, but doesn't work directly:

sinz = sin.(column_center_coords.z)
gradc2f = ClimaCore.Operators.GradientC2F()
## ∇sinz = gradc2f.(sinz) ## this would throw an error
#----------------------------------------------------------------------------

# This throws an error because face values at the boundary are _not_ well-defined:
#
# ```
# ...
#       \
#         ∇θ[2+½]
#       /
# θ[2]
#       \
#         ∇θ[1+½]
#       /
# θ[1]
#       \
#         ????
# ```
#
# To handle boundaries we need to *modify the stencil*. Two options:
# - provide the _value_ $\theta^*$ of $\theta$ at the boundary:
# ```math
# \nabla\theta[\tfrac{1}{2}] = \frac{\theta[1] - \theta^*}{\Delta z /2}
# ```
#
# - provide the *gradient* $\nabla\theta^*$ of $\theta$ at the boundary:
# ```math
# \nabla\theta[\tfrac{1}{2}] = \nabla\theta^*
# ```
#
# These modified stencils are provided as keyword arguments to the operator (based on the boundary label names):

sinz = sin.(column_center_coords.z)
gradc2f = ClimaCore.Operators.GradientC2F(
    bottom = ClimaCore.Operators.SetValue(sin(0.0)),
    top = ClimaCore.Operators.SetGradient(
        ClimaCore.Geometry.WVector(cos(10.0)),
    ),
)
∇sinz = gradc2f.(sinz)
#----------------------------------------------------------------------------

plot(map(x -> x.w, ClimaCore.Geometry.WVector.(∇sinz)), ylim = (0, 10))
#----------------------------------------------------------------------------

# As before, multiple operators (or functions) can be fused together with broadcasting.
#
# One extra advantage of this is that boundaries of the inner operators only need to be specified if they would affect the final result.
#
# Consider the center-to-center Laplacian:
#
# ```
# ...
#       \       /
#         ∇θ[2+½]
#       /       \
# θ[2]            ∇⋅∇θ[2]
#       \       /
#         ∇θ[1+½]
#       /       \
# θ[1]            ∇⋅∇θ[1]
#               /
#          ∇θ*
# ```

sinz = sin.(column_center_coords.z)
## we don't need to specify boundaries, as the stencil won't reach that far
gradc2f = ClimaCore.Operators.GradientC2F()
divf2c = ClimaCore.Operators.DivergenceF2C(
    bottom = ClimaCore.Operators.SetValue(ClimaCore.Geometry.WVector(cos(0.0))),
    top = ClimaCore.Operators.SetValue(ClimaCore.Geometry.WVector(cos(10.0))),
)
∇∇sinz = divf2c.(gradc2f.(sinz))
#----------------------------------------------------------------------------

plot(∇∇sinz, ylim = (0, 10))
#----------------------------------------------------------------------------

# # 3. Solving PDEs
#
# ClimaCore can be used for spatial discretizations of PDEs. For temporal discretization, we can use the OrdinaryDiffEq package, which we aim to be compatibile with.

using OrdinaryDiffEq
#----------------------------------------------------------------------------

# ### 3.1 Heat equation using finite differences
#
# We will use a cell-center discretization of the heat equation:
# ```math
# \frac{\partial y}{\partial t} = \alpha \nabla \cdot \nabla y
# ```
#
# At the bottom we will use a Dirichlet condition ``y(0) = 1``` at the bottom: since we don't actually have a value located at the bottom, we will use a `SetValue` boundary modifier on the inner gradient.
#
# At the top we will use a Neumann condition ``\frac{\partial y}{\partial z}(10) = 0``. We can do this two equivalent ways:
#  - a `SetGradient` on the gradient operator
#  - a `SetValue` on the divergence operator
#
# either will work.

y0 = zeros(column_center_space)

## define the tendency function
function heat_fd_tendency!(dydt, y, α, t)
    gradc2f = ClimaCore.Operators.GradientC2F(
        bottom = ClimaCore.Operators.SetValue(1.0),
        top = ClimaCore.Operators.SetGradient(ClimaCore.Geometry.WVector(0.0)),
    )
    divf2c = ClimaCore.Operators.DivergenceF2C()
    ## the @. macro "dots" the whole expression
    ## i.e.  dydt .= α .* divf2c.(gradc2f.(y))
    @. dydt = α * divf2c(gradc2f(y))
end

heat_fd_prob = ODEProblem(heat_fd_tendency!, y0, (0.0, 5.0), 0.1)
heat_fd_sol = solve(heat_fd_prob, SSPRK33(), dt = 0.1, saveat = 0.25)
#----------------------------------------------------------------------------

anim = Plots.@animate for u in heat_fd_sol.u
    plot(u, xlim = (0, 1), ylim = (0, 10))
end
mp4(anim)
#----------------------------------------------------------------------------

# ### 3.2 Heat equation using continuous Galerkin (CG) spectral element
#
#

function heat_cg_tendency!(dydt, y, α, t)
    grad = ClimaCore.Operators.Gradient()
    wdiv = ClimaCore.Operators.WeakDivergence()
    ## apply element operators
    @. dydt = α * wdiv(grad(y))

    ## direct stiffness summation (DSS): project to continuous function space
    ClimaCore.Spaces.weighted_dss!(dydt)
    return dydt
end

y0 = exp.(.-(coord.y .^ 2 .+ coord.x .^ 2) ./ 2)

heat_cg_prob = ODEProblem(heat_cg_tendency!, y0, (0.0, 5.0), 0.1)
heat_cg_sol = solve(heat_cg_prob, SSPRK33(), dt = 0.1, saveat = 0.5)
#----------------------------------------------------------------------------

anim = Plots.@animate for u in heat_cg_sol.u
    Plots.plot(u, c = :thermal)
end
mp4(anim)
#----------------------------------------------------------------------------

# ### 3.3 Shallow water equations
#
#
# The shallow water equations in vector invariant form can be written as
# ```math
# \begin{align*}
#     \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho u) &= 0\\
#     \frac{\partial u_i}{\partial t} + \nabla (\Phi + \tfrac{1}{2}\|u\|^2)_i  &= J (u \times (\nabla \times u))_i
# \end{align*}
# ```
# where ``J`` is the Jacobian determinant, and ``\Phi = g \rho``.
#
# Note that the velocity ``u`` is specified in _covariant_ coordinates ``u_i``.
#
# For vizualization purposes, we can model a passive tracer $\theta$ as
# ```math
# \frac{\partial \rho \theta}{\partial t} + \nabla \cdot (\rho \theta u) = 0
# ```

using ClimaCore.Geometry

parameters = (
    ϵ = 0.1,  ## perturbation size for initial condition
    l = 0.5, ## Gaussian width
    k = 0.5, ## Sinusoidal wavenumber
    ρ₀ = 1.0, ## reference density
    c = 2,
    g = 10,
    D₄ = 1e-4, ## hyperdiffusion coefficient
)

function init_state(local_geometry, p)
    coord = local_geometry.coordinates
    @unpack x, y = coord
    ## set initial state
    ρ = p.ρ₀

    ## set initial velocity
    U₁ = cosh(y)^(-2)

    ## Ψ′ = exp(-(x2 + p.l / 10)^2 / 2p.l^2) * cos(p.k * x) * cos(p.k * y)
    ## Vortical velocity fields (u₁′, u₂′) = (-∂²Ψ′, ∂¹Ψ′)
    ϕ = exp(-(y + p.l / 10)^2 / 2p.l^2)
    u₁′ = ϕ * (y + p.l / 10) / p.l^2 * cos(p.k * x) * cos(p.k * y)
    u₁′ += p.k * ϕ * cos(p.k * x) * sin(p.k * y)
    u₂′ = -p.k * ϕ * sin(p.k * x) * cos(p.k * y)

    u = Geometry.Covariant12Vector(
        Geometry.UVVector(U₁ + p.ϵ * u₁′, p.ϵ * u₂′),
        local_geometry,
    )

    ## set initial tracer
    θ = sin(p.k * y)
    return (ρ = ρ, u = u, ρθ = ρ * θ)
end

y0 =
    init_state.(
        ClimaCore.Fields.local_geometry_field(rectangle_space),
        Ref(parameters),
    )

## plot initial tracer
Plots.plot(y0.ρθ)
#----------------------------------------------------------------------------

function shallow_water_tendency!(dydt, y, _, t)

    @unpack D₄, g = parameters

    sdiv = ClimaCore.Operators.Divergence()
    wdiv = ClimaCore.Operators.WeakDivergence()
    grad = ClimaCore.Operators.Gradient()
    wgrad = ClimaCore.Operators.WeakGradient()
    curl = ClimaCore.Operators.Curl()
    wcurl = ClimaCore.Operators.WeakCurl()

    ## compute hyperviscosity first
    @. dydt.u =
        wgrad(sdiv(y.u)) -
        Geometry.Covariant12Vector(wcurl(Geometry.Covariant3Vector(curl(y.u))))
    @. dydt.ρθ = wdiv(grad(y.ρθ))

    ClimaCore.Spaces.weighted_dss!(dydt)

    @. dydt.u =
        -D₄ *
        (wgrad(sdiv(dydt.u)) - Geometry.Covariant12Vector(wcurl(curl(dydt.u))))
    @. dydt.ρθ = -D₄ * wdiv(grad(dydt.ρθ))

    ## comute rest of tendency
    @. begin
        dydt.ρ = -wdiv(y.ρ * y.u)
        dydt.u += -grad(g * y.ρ + norm(y.u)^2 / 2) + y.u × curl(y.u)
        dydt.ρθ += -wdiv(y.ρθ * y.u)
    end
    ClimaCore.Spaces.weighted_dss!(dydt)
    return dydt
end

#----------------------------------------------------------------------------

shallow_water_prob = ODEProblem(shallow_water_tendency!, y0, (0.0, 20.0))
@time shallow_water_sol =
    solve(shallow_water_prob, SSPRK33(), dt = 0.05, saveat = 1.0)
anim = Plots.@animate for u in shallow_water_sol.u
    Plots.plot(u.ρθ, clim = (-1, 1))
end
mp4(anim)
#----------------------------------------------------------------------------

#----------------------------------------------------------------------------


text:
include("./exportData.jl")
include("./fenkar_ode.jl")
using .exportData, .fenkar_ode, DifferentialEquations, DiffEqParamEstim, PyPlot, BlackBoxOptim, Optim, DelimitedFiles

#solve and take data
prob = ODEProblem(fenkar!, u0, tspan, p)
sol = solve(prob, Tsit5(); saveat=tspan[begin]:5.0:tspan[end])

# getExpData
ind = 78
t, data, p[15] = getExpData(ind, tInds=1:500)
p[14] = 0.0 #t[findfirst(data[1,:] .> 0.5)]-p[15]/2

# remake problem
prob = remake(prob, tspan=(0.0, t[end]), p=p)

# use the save_idxs to select the mid-point of the u-variable
sol = solve(prob, Tsit5(); saveat=t, save_idxs=1:1)

fig = plt.figure(constrained_layout=true)
plt.plot(t, data[:], "-", label="Exp. Sample")
plt.plot(sol.t, sol[1,:], "--", label="Initial Parameters")
plt.xlabel("\$ t \$")

# define the cost function
cost_function = build_loss_objective(prob, Tsit5(), L2Loss(t,data), maxiters=10000, verbose_opt=true, verbose_steps=50; save_idxs=1:1);

# define limits for parameters
#          tsi,  tv1m,   tv2m,   tvp,   twm,    twp,  td,    to,    tr,   xk,  uc,    uv, ucsi,        t0,     TI,   IA
#   p = [ 29.0,  19.6, 1250.0,   3.33, 41.0,  870.0, 0.25, 12.5,  33.3, 10.0, 0.13, 0.04,  0.85,      1.0,   BCL, 0.05]
lower = [  1.0,  10.0,  400.0,   2.33, 10.0,  100.0, 0.10, 10.0,  10.0,  9.0, 0.10, 0.01,  0.10,      0.0,   40.0, 0.00]
upper = [ 50.0,  50.0, 2000.0,   4.33, 50.0, 1000.0, 0.50, 50.0,  50.0, 11.0, 0.90, 0.90,  0.90,    p[15], 1000.0, 1.00]

@assert all(p .>= lower) && all(p .<= upper)

bound = [(lower[n], upper[n]) for n in 1:length(p)]

# fit
result = bboptimize(cost_function; SearchRange = bound, MaxSteps = 1e5)#, NThreads=Threads.nthreads()-1)
@info "Initial p = $(p) \n"
@info "Minimizer = $(result.archive_output.best_candidate) \n"
prob = remake(prob, p=result.archive_output.best_candidate)
sol = solve(prob, Tsit5(); saveat=t, save_idxs=1:1)
plt.plot(sol.t, sol[1,:], "--", label="BBO Parameters")

result = optimize(cost_function, result.archive_output.best_candidate, NelderMead())
@info "Initial p = $(p) \n"
@info "Minimizer = $(result.minimizer) \n"
prob = remake(prob, p=result.minimizer)
sol = solve(prob, Tsit5(); saveat=t, save_idxs=1:1)
plt.plot(sol.t, sol[1,:], "--", label="Nelder-Mead Parameters")

plt.savefig("./ode/bbo/$(ind)_fittings.pdf")
plt.close()

open("./ode/bbo/$(ind)_params.txt", "w") do io
	writedlm(io, ["tsi" "tv1m" "tv2m" "tvp" "twm" "twp" "td" "to" "tr" "xk" "uc" "uv" "ucsi" "t0" "TI" "IA"])
	writedlm(io, "\n")
	writedlm(io, transpose(p))
	writedlm(io, "\n")
	#writedlm(io, transpose(result.archive_output.best_candidate))
	writedlm(io, transpose(result.minimizer))
end

ext: jl


text:
#predevelopment and consenting

function OMCost!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)
    key = (options.postprocess ? "postprocessing" : "evaluating")
    println(key*" operation cost")
    operation!(options,sizes,data,output)
    println(key*" maintenance cost")
    maintenance!(options,sizes,data,output)
    nothing

end

function operation!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)

    if options.operation == 0
        if options.postprocess
            nothing
        else
            Crent = (data.l .* data.Pe) .* output.AEP * data.duration .* transpose(data.turbineQuantity)
            CinsO = zeros(sizes.winds,sizes.turbines)
            CtransO = zeros(sizes.winds,sizes.turbines)
            for i in 1:sizes.winds,j in 1:sizes.turbines
                CinsO[i,j] += data.ComIns[i] * data.IC[j]
                CtransO[i,j] += data.Ctransunit[i] * data.IC[j]
            end
            output.OPEX += Crent + CinsO + CtransO
            output.operation += Crent + CinsO + CtransO
        end
    else
        error("operation mode not found.")
    end
    nothing

end

function maintenance!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)

    if options.maintenance == 0
        if options.postprocess
            nothing
        else
            Ccm = zeros(sizes.winds,sizes.turbines)
            Cmdir = zeros(sizes.winds,sizes.turbines)
            Cmind = zeros(sizes.winds,sizes.turbines)

            Ctrans = 2*(data.d .* data.tc)
            Clab = data.N1dOM .* data.LrOM
            Ccm .+= Ctrans .+ transpose(Clab) .+ data.consumCost

            for i in 1:sizes.winds,j in 1:sizes.turbines
                Cmdir[i,j] = (1-data.Pd[j])*data.lambda[j]*Ccm[i,j] + data.Pd[j]*data.Csm[j]
            end

            Cmind .+= (data.Cindport + data.Cindves + data.Cindlab).*transpose(data.turbineQuantity)

            output.OPEX += Cmdir + Cmind
            output.maintenance += Cmdir + Cmind
        end
    else
        error("maintenance mode not found.")
    end
    nothing

end

max_stars_repo_name: Thiago-NovaesB/zefiro

max_issues_repo_name: Thiago-NovaesB/zefiro

max_forks_repo_name: Thiago-NovaesB/zefiro


text:
"""
    SeisDelay(d ; <keyword arguments>)

Apply a time delay to a 2D array input.

# Arguments
* `d`: Input 2D data array in tx domain.

# Keyword arguments
* `delay=0.1`: time delay in seconds.
* `dt=0.001`: time sampling

# Output
* `d2`: Delayed data in time domain

# Example
```julia
julia> d = SeisLinearEvents(d); SeisPlot(d);
julia> d2 = SeisDelay(d;dt=0.004); SeisPlot(d);
"""

function SeisDelay(d; delay=0.1, dt=0.001)

    nt = size(d,1)
    nf = nt
    D = fft(d,1);
    dw = 2*pi/nt/dt
    nw = round(Int,nt/2) + 1
    for iw = 1 : nw
	w = (iw-1)*dw
	D[iw,:] = D[iw,:]*exp(-1im*w*delay)
    end
    # symmetries
    for iw=nw+1:nf
	D[iw,:] = conj(D[nf-iw+2,:])
    end
    d2 = ifft(D,1)
    d2 = real(d2[1:nt,:])
    return d2

end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
# These functions are not included yet because I am not yet sure which plotting library to put my bets on.
# Makie&StatsMakie could be supercool, but I ran into some issues
import Plots
using SparseArrays
Plots.plot(m::unfold.UnfoldModel)  = plot_results(m.results)

function plot_results(results::DataFrame)
    Plots.plot(results.time,results.estimate,
            group=(results.term,results.group),
            legend=:outerbottom)

end

Plots.heatmap(X::SparseMatrixCSC) = heatmap(Matrix(X))


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

using Markdown
using InteractiveUtils

# ╔═╡ cae0afe4-1cd9-11eb-1043-ebcfa0d7db13
# Nota: Aquí el bloque begin-end es necesario por la forma en que opera nuestro
#       actual IDE, Pluto Notebook. En un script independiente de Julia no es
#       necesario de utilizarse, aunque es válido para encapsulamiento.
begin
	import Pkg
	Pkg.activate(".")
	using DifferentialEquations, Plots
end

# ╔═╡ 33d5d3f0-1cd8-11eb-2a7c-a1e41bd39d42
md"# ¿Qué es y qué no es este taller?
## Qué no es...
* No será una introducción completa al lenguaje de programación, Julia (Muy poco tiempo). Aunque **no** se asume conocimiento previo de todas maneras.
* No es un taller que pueda reemplazar completamente uno de análisis numérico
## Qué es...
* Un vistazo a la realidad de lo que puede ser crear modelos y utilizar métodos numéricos para cómputo científico en el lenguaje Julia.
* Una introducción a tópicos aplicados en ciencias e ingeniería.
# ¿Cómo usar este proyecto?
Si se desea utilizar como proyecto autocontenido, el cuaderno actual provee ya la interfaz necesaria para funcionar al abrir. No obstante, si se desean crear otros archivos y manipularlos se proveen el `Project.toml` y `Manifest.toml` para utilizar una vez dentro del directorio, el comando
```Julia
Pkg> activate .
Pkg> instantiate
```
y de forma opcional:
```Julia
Pkg> update; precompile
```
Para inicializar la paquetería (descargar y precompilar) necesaria. Esto ya se realiza dentro de este cuaderno a continuación para facilidad de uso:
"

# ╔═╡ 9de1d230-1da5-11eb-1111-ab0b443361e1
md"# Breve introducción al lenguaje
## Cualidades y filosofías
Comenzamos haciendo una comparación con otros lenguajes populares. Julia toma las siguientes características altamente buscadas de lenguajes existentes:
* Interactivo (como Python)
* Rápido, compitiendo con velocidades de C y Fortran
* De código abierto (a diferencia de C\# y el lenguaje Wolfram en Mathematica)
* Especializable a limpieza de datos y modelaje estadístico (similar a R)
* Optimizado para cómputo científico (como Octave)
* Realiza cómputo y modelaje simbólico (como Mathematica)
* Poderoso en metaprogramación (como Ruby)

Algunas de las cosas que julia busca priorizar son:
* Legibilidad del código (hasta llegar a ser tan expresivo como lo que escribiríamos en un cuaderno)
* Reducir la brecha entre el concepto y el código: Escribimos como pensamos.
* La fácil optimización hacia código eficiente
* Enfoque a \"verbos\" (métodos) en lugar de \"sustantivos\" (clases). Esto principalmente por lo que llamamos Multiple Dispatch

## Multiple Dispatch
En Julia, cualquier tipo de dato puede realizar cualquier acción (método) una vez que se le sea definido. Esto significa que podemos enfocarnos en crear comportamientos que necesitemos sin pensar en que deben estar asociados a una clase para existir.

Esto se lleva un paso más allá con **multiple dispatch**, que nos permite definir polimorficamente dicho método. Es decir que una acción, digamos `correr`, puede actuar de diferentes maneras y llamarse igual, decidiéndose su acción en el contexto de sus argumentos. Ejemplo, `correr(server::Server)` y `correr(bruno::Perro)` son la misma función, pero despachada con distinto comportamiento al momento de evaluarse.
"

# ╔═╡ c0523206-1da5-11eb-1b19-cb7103b63d28
md"# Modelos
## Preliminares de ecuaciones diferenciales
Gracias a la gran eficiencia de métodos iterativos y de álgebra lineal, la forma preferencial de pensar en ecuaciones diferenciales es en su forma de *sistema de ecuaciones de primer órden*. Esto siempre es posible de obtener mediante [reducción de órden](https://en.wikipedia.org/wiki/Ordinary_differential_equation#Reduction_of_order).

Una vez realizado, la forma que obtenemos, para ecuaciones autónomas (no dependientes de manera explícita del tiempo), es:

$$
x'(t) = F(x(t))
$$

donde $x: \mathbb{R} \rightarrow \mathbb{R}^n$ es una trayectoria en el espacio para la cual deseamos resolver y obedece la dinámica del campo vectorial:

$$F: \mathbb{R}^n \rightarrow \mathbb{R}^n$$

Esto define un sistema de $n$ ecuaciones a resolver. A continuación veremos algunos ejemplos.

## 1. Péndulo caótico
El péndulo doble es un sistema muy famosamente estudiado por exhibir un comportamiento **caótico** (cuya definición matemática es rigurosa). Sus ecuaciones del movimiento son las siguientes:

$$\frac{d}{dt}
\begin{pmatrix}
\alpha \\ l_\alpha \\ \beta \\ l_\beta
\end{pmatrix}=
\begin{pmatrix}
2\frac{l_\alpha - (1+\cos\beta)l_\beta}{3-\cos 2\beta} \\
-2\sin\alpha - \sin(\alpha + \beta) \\
2\frac{-(1+\cos\beta)l_\alpha + (3+2\cos\beta)l_\beta}{3-\cos2\beta}\\
-\sin(\alpha+\beta) - 2\sin(\beta)\frac{(l_\alpha-l_\beta)l_\beta}{3-\cos2\beta} + 2\sin(2\beta)\frac{l_\alpha^2-2(1+\cos\beta)l_\alpha l_\beta + (3+2\cos\beta)l_\beta^2}{(3-\cos2\beta)^2}
\end{pmatrix}$$

Resolveremos este sistema utilizando los métodos del paquetes `DifferentialEquations.jl`.

El movimiento generado por esta ecuaciones se ve similar a:
"

# ╔═╡ 70ac4432-1dad-11eb-2ffc-bdb991592480
html"

"

# ╔═╡ 02926624-1da9-11eb-1b45-b5268512ebb9
# El prefijo const es opcional y no significa VALOR constante, si no tipo constante.
begin
	const m₁, m₂, L₁, L₂, g = 1, 2, 1, 2, 9.81
	initial = [0, π/3, 0, 3pi/5]
	tspan = (0., 50.)
end;

# ╔═╡ 80b59058-1da9-11eb-3e67-8170a6930d2b
md"Se define una función auxiliar para transformar de coordenadas polares a cartesianas:"

# ╔═╡ 7449c578-1da9-11eb-0373-d529c7d1aa51
function polar2cart(sol; dt=0.02, l1=L₁, l2=L₂, vars=(2,4))
    u = sol.t[1]:dt:sol.t[end]

    p1 = l1*map(x->x[vars[1]], sol.(u))
    p2 = l2*map(y->y[vars[2]], sol.(u))

    x1 = l1*sin.(p1)
    y1 = l1*-cos.(p1)
    (u, (x1 + l2*sin.(p2),
     y1 - l2*cos.(p2)))
end

# ╔═╡ 8bbe3dec-1da9-11eb-3376-47c99c6bd393
function double_pendulum(xdot,x,p,t)
    xdot[1] = x[2]
    xdot[2] = - ((g*(2*m₁+m₂)*sin(x[1]) + m₂*(g*sin(x[1]-2*x[3]) +
				2*(L₂*x[4]^2+L₁*x[2]^2*cos(x[1]-x[3]))*sin(x[1]-x[3])))/
		        (2*L₁*(m₁+m₂-m₂*cos(x[1]-x[3])^2)))
    xdot[3] = x[4]
    xdot[4] = (((m₁+m₂)*(L₁*x[2]^2+g*cos(x[1])) +
			   L₂*m₂*x[4]^2*cos(x[1]-x[3]))*sin(x[1]-x[3]))/
			   (L₂*(m₁+m₂-m₂*cos(x[1]-x[3])^2))
end

# ╔═╡ d2a2721e-1da9-11eb-1eea-d3a584fafa8e
double_pendulum_problem = ODEProblem(double_pendulum, initial, tspan);

# ╔═╡ ec5d56a6-1da9-11eb-1fc4-57840baed98d
sol = solve(double_pendulum_problem, Vern7(), abs_tol=1e-10, dt=0.05)

# ╔═╡ 3317fbdc-1daa-11eb-1d17-6791df36b6e4
begin
	ts, ps = polar2cart(sol, l1=L₁, l2=L₂, dt=0.01)
	plot(ps...)
end

# ╔═╡ 61376792-1dae-11eb-02fa-93e6782be478
md"## 2. Sistema de Hénon-Heiles
Es un sistema dinámico que modela el movimiento de una estrella orbitando alrededor de su centro galáctico mientras yace restringido en un plano. Éste es un ejemplo de un sistema Hamiltoniano, y tiene la forma:

$$
\begin{align}
\frac{d^2x}{dt^2}&=-\frac{\partial V}{\partial x}\\
\frac{d^2y}{dt^2}&=-\frac{\partial V}{\partial y}
\end{align}
$$

donde

$$V(x,y)={\frac {1}{2}}(x^{2}+y^{2})+\lambda \left(x^{2}y-{\frac {y^{3}}{3}}\right).$$

es conocido como el **potencial de Hénon–Heiles**. De éste puede derivarse su Hamiltoniano:

$$H={\frac {1}{2}}(p_{x}^{2}+p_{y}^{2})+{\frac {1}{2}}(x^{2}+y^{2})+\lambda \left(x^{2}y-{\frac {y^{3}}{3}}\right).$$

Esta cantidad representa un invariante del sistema dinámico: Una cantidad conservada. En este caso, es la energía total del sistema.
"

# ╔═╡ cc7905be-1daf-11eb-33d3-9d354f8d9419
begin
	# Parámetros
	initial₂ = [0.,0.1,0.5,0]
	tspan₂ = (0,100.)
	# V: Potencial, T: Energía cinética total, E: Energía total
	V(x,y) = 1//2 * (x^2 + y^2 + 2x^2*y - 2//3 * y^3)
	E(x,y,dx,dy) = V(x,y) + 1//2 * (dx^2 + dy^2);
end;

# ╔═╡ ab323b54-1db0-11eb-0d16-2b909c5a4a16
md"Definimos el modelo en una función:"

# ╔═╡ 4c4d367a-1db0-11eb-0363-27634674bbfe
function Hénon_Heiles(du,u,p,t)
    x  = u[1]
    y  = u[2]
    dx = u[3]
    dy = u[4]
    du[1] = dx
    du[2] = dy
    du[3] = -x - 2x*y
    du[4] = y^2 - y -x^2
end

# ╔═╡ a9f31cc2-1db0-11eb-3df5-2f65f8d96d67
md"Resolvemos el problema:"

# ╔═╡ 279aff64-1db1-11eb-2ce5-7955f39dd250
begin
	prob₂ = ODEProblem(Hénon_Heiles, initial₂, tspan₂)
	sol₂ = solve(prob₂, Vern9(), abs_tol=1e-16, rel_tol=1e-16);
end

# ╔═╡ e47793e6-1db0-11eb-2ac5-511dbb2dbaa8
plot(sol₂, vars=(1,2), title = "Órbita del sistema de Hénon-Heiles", xaxis = "x", yaxis = "y", leg=false)

# ╔═╡ 3a2ad3a6-1db2-11eb-136e-93feebbe9377
md"Parece estar correctamente resuelto pero... examinando la evolución de la energía total/Hamiltoniano podemos encontrar lo siguiente:"

# ╔═╡ 0fb21260-1db2-11eb-211b-d7a3cc0c3aa9
begin
	energy = map(x->E(x...), sol₂.u)
	plot(sol₂.t, energy .- energy[1], title = "Cambio de la energía en el tiempo", xaxis = "Número de iteraciones", yaxis = "Cambio en energía")
end

# ╔═╡ 168fa38a-1db3-11eb-02f8-b547454f2511
md"¡La energía total cambia! Eso quiere decir que la ley de conservación que esperamos que se cumpla físicamente no parece cumplirse en nuestra simulación. Esto delata un error en la resolución de la ecuación, específicamente por el método utilizado.

Para evitar eso, podemos utilizar algo conocido como un **integrador simplético** que considere la estructura de sistema Hamiltoniano que tiene nuestras ecuaciones.

Éste lo implementamos a continuación:"

# ╔═╡ 64f336f8-1db4-11eb-3d08-45ce58c5a127
function HH_acceleration!(dv,v,u,p,t)
    x,y  = u
    dx,dy = dv
    dv[1] = -x - 2x*y
    dv[2] = y^2 - y -x^2
end;

# ╔═╡ 7375b80e-1db4-11eb-395d-4f018f0271dc
md"Debemos ahora definir la condición inicial por separado, pues nuestro sistema, al ser Hamiltoniano, tiene una segmentación natural en esos pares de variables."

# ╔═╡ 67f80a2c-1db4-11eb-0f13-4bacb5d52f03
begin
	initial_positions = [0.0,0.1]
	initial_velocities = [0.5,0.0]
end;

# ╔═╡ b4ae1924-1db4-11eb-2f1e-b398cd26c703
md"Resolvemos el problema ahora como una ecuación de segundo orden pero con estructura simpléctica detectada:"

# ╔═╡ 72822478-1db4-11eb-3b89-13dc7ce09df1
begin
	prob₃ = SecondOrderODEProblem(HH_acceleration!,
								 initial_velocities,
								 initial_positions,
								 tspan₂)
	sol₃ = solve(prob₃, KahanLi8(), dt=1/10)
end;

# ╔═╡ 0ffa749e-1db5-11eb-072f-67d20fd86af5
md"Ahora podemos observar cómo la energía se mantiene muy cercana a cero, oscilando solamente por errores de precisión numérica pero sin existir una tendencia de crecimiento anómala."

# ╔═╡ e20eb43c-1db4-11eb-3fcb-7fbd0ff32e0c
begin
	energy₂ = map(x->E(x[3], x[4], x[1], x[2]), sol₃.u)
	plot(sol₃.t, energy₂ .- energy₂[1], title = "Cambio de la energía en el tiempo", xaxis = "Número de iteraciones", yaxis = "Cambio en energía")
end

# ╔═╡ 4812b46e-1daf-11eb-0ff3-4f5706cd25b3
md"
# Bibliografía
* [Documentación de SciML](https://sciml.ai/)
* [Tutoriales de SciML](https://github.com/SciML/SciMLTutorials.jl)
* Hénon, Michel (1983), \"Numerical exploration of Hamiltonian Systems\", in Iooss, G. (ed.), Chaotic Behaviour of Deterministic Systems, Elsevier Science Ltd, pp. 53–170, ISBN 044486542X
* Aguirre, Jacobo; Vallejo, Juan C.; Sanjuán, Miguel A. F. (2001-11-27). \"Wada basins and chaotic invariant sets in the Hénon-Heiles system\". Physical Review E. American Physical Society (APS). 64 (6): 066208. doi:10.1103/physreve.64.066208. ISSN 1063-651X.
* Steven Johnson. 18.335J Introduction to Numerical Methods . Spring 2019. Massachusetts Institute of Technology: MIT OpenCourseWare, https://ocw.mit.edu. License: Creative Commons BY-NC-SA.
"

# ╔═╡ Cell order:
# ╟─33d5d3f0-1cd8-11eb-2a7c-a1e41bd39d42
# ╠═cae0afe4-1cd9-11eb-1043-ebcfa0d7db13
# ╟─9de1d230-1da5-11eb-1111-ab0b443361e1
# ╟─c0523206-1da5-11eb-1b19-cb7103b63d28
# ╟─70ac4432-1dad-11eb-2ffc-bdb991592480
# ╠═02926624-1da9-11eb-1b45-b5268512ebb9
# ╟─80b59058-1da9-11eb-3e67-8170a6930d2b
# ╠═7449c578-1da9-11eb-0373-d529c7d1aa51
# ╠═8bbe3dec-1da9-11eb-3376-47c99c6bd393
# ╠═d2a2721e-1da9-11eb-1eea-d3a584fafa8e
# ╠═ec5d56a6-1da9-11eb-1fc4-57840baed98d
# ╟─3317fbdc-1daa-11eb-1d17-6791df36b6e4
# ╟─61376792-1dae-11eb-02fa-93e6782be478
# ╠═cc7905be-1daf-11eb-33d3-9d354f8d9419
# ╟─ab323b54-1db0-11eb-0d16-2b909c5a4a16
# ╠═4c4d367a-1db0-11eb-0363-27634674bbfe
# ╟─a9f31cc2-1db0-11eb-3df5-2f65f8d96d67
# ╠═279aff64-1db1-11eb-2ce5-7955f39dd250
# ╟─e47793e6-1db0-11eb-2ac5-511dbb2dbaa8
# ╟─3a2ad3a6-1db2-11eb-136e-93feebbe9377
# ╟─0fb21260-1db2-11eb-211b-d7a3cc0c3aa9
# ╟─168fa38a-1db3-11eb-02f8-b547454f2511
# ╠═64f336f8-1db4-11eb-3d08-45ce58c5a127
# ╟─7375b80e-1db4-11eb-395d-4f018f0271dc
# ╠═67f80a2c-1db4-11eb-0f13-4bacb5d52f03
# ╟─b4ae1924-1db4-11eb-2f1e-b398cd26c703
# ╠═72822478-1db4-11eb-3b89-13dc7ce09df1
# ╟─0ffa749e-1db5-11eb-072f-67d20fd86af5
# ╟─e20eb43c-1db4-11eb-3fcb-7fbd0ff32e0c
# ╟─4812b46e-1daf-11eb-0ff3-4f5706cd25b3

max_stars_repo_name: anathompson/taller-julia-macaotech

max_issues_repo_name: anathompson/taller-julia-macaotech

max_forks_repo_name: anathompson/taller-julia-macaotech


text:
module FaceLabelsModule

using Gridap
using Gridap.Helpers
using Gridap.CellValuesGallery

export FaceLabels
export labels_on_dim
export labels_on_tag
export ntags
export tag_from_name
export name_from_tag
export first_tag_on_face
export add_tag_from_tags!

"""
Classification of nfaces into geometrical and physical labels
"""
struct FaceLabels
  dim_to_nface_to_label::Vector{IndexCellValue}
  tag_to_labels::Vector{Vector{Int}}
  tag_to_name::Vector{String}
end

function FaceLabels(
  dim_to_nface_to_label::Vector{Vector{Int}},
  physlabel_to_labels::Vector{Vector{Int}},
  tag_to_name::Vector{String})

  cv = [ CellValueFromArray(v) for v in dim_to_nface_to_label ]
  FaceLabels(cv, physlabel_to_labels, tag_to_name)
end

"""
Returns an AbstractVector{Int} that represent the label for
each nface of dimension dim
"""
labels_on_dim(l::FaceLabels,dim::Integer) = l.dim_to_nface_to_label[dim+1]

"""
Returns a Vector{Int} with the labels associated with a given phystag
"""
labels_on_tag(l::FaceLabels,tag::Integer) = l.tag_to_labels[tag]

ntags(l::FaceLabels) = length(l.tag_to_labels)

function tag_from_name(l::FaceLabels,name::String)
  for tag in 1:ntags(l)
    if l.tag_to_name[tag] == name
      return tag
    end
  end
  @unreachable
  0
end

name_from_tag(l::FaceLabels,tag::Integer) = l.tag_to_name[tag]

function first_tag_on_face(labels::FaceLabels,dim::Integer)
  face_to_label = labels_on_dim(labels,dim)
  nfaces = length(face_to_label)
  face_to_tag = zeros(Int,nfaces)
  _first_tag_on_face!(face_to_tag,face_to_label,labels.tag_to_labels)
  CellValueFromArray(face_to_tag)
end

function _first_tag_on_face!(face_to_tag,face_to_label,tag_to_labels)
  nfaces = length(face_to_tag)
  for face in 1:nfaces
    label = face_to_label[face]
    for (tag,labels) in enumerate(tag_to_labels)
      if label in labels
        face_to_tag[face] = tag
      end
    end
  end
end

function add_tag_from_tags!(facelabels::FaceLabels, name::String, tags::Vector{Int})
  @assert ! (name in facelabels.tag_to_name)
  labels = Int[]
  for tag in tags
     for label in labels_on_tag(facelabels,tag)
       push!(labels,label)
     end
  end
  push!(facelabels.tag_to_labels,sort(collect(Set(labels))))
  push!(facelabels.tag_to_name,name)
end

function add_tag_from_tags!(
  facelabels::FaceLabels, name::String, names::Vector{String})
  tags = [tag_from_name(facelabels,name) for name in names ]
  add_tag_from_tags!(facelabels,name,tags)
end

# Pretty printing

import Base: show

function show(io::IO,self::FaceLabels)
  print(io,"FaceLabels object")
end

function show(io::IO,::MIME"text/plain",labels::FaceLabels)
  show(io,labels)
  print(io,":")
  D = length(labels.dim_to_nface_to_label)-1
  for d = 0:D
    print(io,"\n $d-faces: $(length(labels_on_dim(labels,d)))")
  end
  print(io,"\n tags: $(ntags(labels))")
end

end #module


text:
function bfs(maze::Maze)::Array{Tuple{Int, Int}, 1}
	que::Queue{Tuple{Int, Int}} = Queue{Tuple{Int, Int}}()
	if maze.entr[1] == 1
		start = (maze.entr[1] + 1, maze.entr[2])
		enqueue!(que, start)
		set_mazecolor(maze, start, 0xffafafaa)
	elseif maze.entr[2] == 1
		start = (maze.entr[1], maze.entr[2] + 1)
		enqueue!(que, start)
		set_mazecolor(maze, start, 0xffafafab)
	elseif maze.entr[1] == maze.w
		start = (maze.entr[1] - 1, maze.entr[2])
		enqueue!(que, start)
		set_mazecolor(maze, start, 0xffafafac)
	else
		start = (maze.entr[1], maze.entr[2] - 1)
		enqueue!(que, start)
		set_mazecolor(maze, start, 0xffafafad)
	end
	set_mazecolor(maze, maze.entr, 0xff00ff00)

	maxque = 1

	while length(que) > 0
		cur_pos = dequeue!(que)
		if cur_pos == maze.exit
			break
		end
		dirs = get_dirs(maze, cur_pos)
		if dirs & 1 != 0
			chk_pos = (cur_pos[1] + 1, cur_pos[2])
			set_mazecolor(maze, chk_pos, 0xffafafaa)
			enqueue!(que, chk_pos)
		end
		if dirs & 2 != 0
			chk_pos = (cur_pos[1], cur_pos[2] + 1)
			set_mazecolor(maze, chk_pos, 0xffafafab)
			enqueue!(que, chk_pos)
		end
		if dirs & 4 != 0
			chk_pos = (cur_pos[1] - 1, cur_pos[2])
			set_mazecolor(maze, chk_pos, 0xffafafac)
			enqueue!(que, chk_pos)
		end
		if dirs & 8 != 0
			chk_pos = (cur_pos[1], cur_pos[2] - 1)
			set_mazecolor(maze, chk_pos, 0xffafafad)
			enqueue!(que, chk_pos)
		end
		maxque = maxque < length(que) ? length(que) : maxque
	end

	# println(maxque)

	if get_mazecolor(maze, maze.exit) == 0xffffffff
		return []
	end

	ret = [maze.exit]

	cur_x, cur_y = maze.exit

	while (cur_x, cur_y) != maze.entr
		cur_col = get_mazecolor(maze, (cur_x, cur_y))
		set_mazecolor(maze, (cur_x, cur_y), 0xffff9400)
		if cur_col == 0xffafafaa
			cur_x -= 1
		elseif cur_col == 0xffafafab
			cur_y -= 1
		elseif cur_col == 0xffafafac
			cur_x += 1
		elseif cur_col == 0xffafafad
			cur_y += 1
		end
		push!(ret, (cur_x, cur_y))
	end

	set_mazecolor(maze, maze.exit, 0xff0000ff)

	reverse!(ret)

	return ret
end


text:
function RLCore.Experiment(
    ::Val{:JuliaRL},
    ::Val{:VPG},
    ::Val{:CartPole},
    ::Nothing;
    save_dir = nothing,
    seed = 123,
)
    if isnothing(save_dir)
        t = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
        save_dir = joinpath(pwd(), "checkpoints", "JuliaRL_VPG_CartPole_$(t)")
    end

    lg = TBLogger(joinpath(save_dir, "tb_log"), min_level = Logging.Info)
    rng = StableRNG(seed)
    env = CartPoleEnv(; T = Float32, rng = rng)
    ns, na = length(get_state(env)), length(get_actions(env))

    agent = Agent(
        policy = VPGPolicy(
            approximator = NeuralNetworkApproximator(
                model = Chain(
                    Dense(ns, 128, relu; initW = glorot_uniform(rng)),
                    Dense(128, 128, relu; initW = glorot_uniform(rng)),
                    Dense(128, na; initW = glorot_uniform(rng)),
                ),
                optimizer = ADAM(),
            ) |> cpu,
            baseline = NeuralNetworkApproximator(
                model = Chain(
                    Dense(ns, 128, relu; initW = glorot_uniform(rng)),
                    Dense(128, 128, relu; initW = glorot_uniform(rng)),
                    Dense(128, 1; initW = glorot_uniform(rng)),
                ),
                optimizer = ADAM(),
            ) |> cpu,
            action_space = get_actions(env),
            dist = Categorical,
            γ = 0.99f0,
            rng = rng,
        ),
        trajectory = ElasticSARTTrajectory(
            state = Vector{Float32} => (ns,),
        ),
    )
    # VPG is updated after each episode
    stop_condition = StopAfterEpisode(500)

    total_reward_per_episode = TotalRewardPerEpisode()
    time_per_step = TimePerStep()
    hook = ComposedHook(
        total_reward_per_episode,
        time_per_step,
        DoEveryNEpisode() do t, agent, env
            with_logger(lg) do
                @info(
                    "training",
                    loss = agent.policy.loss,
                    baseline_loss = agent.policy.baseline_loss,
                    reward = total_reward_per_episode.rewards[end],
                )
            end
        end,
    )

    description = "# Play CartPole with VPG"

    Experiment(agent, env, stop_condition, hook, description)
end


text:
# This example predicts the temporal evolution of a field
# consisting of coupled discrete systems, U,
# which has to be represented as vectors of vectors.
# Where the field comes from does not matter, but to make the example
# runnable we load one of the test systems of TimeseriesPrediction.
#
# In this example we use many coupled henon maps, which can be thought of
# as a field, where each map is a point in space.
#
# Importantly, the results are compared with the "real" evolution of the
# system.
#
# In the plots, the x axis is space and y axis is time.

# ### Simulate a test system
using PyPlot
using TimeseriesPrediction

testdir = dirname(dirname(pathof(TimeseriesPrediction)))*"/test"
@assert isdir(testdir)
include(testdir*"/system_defs.jl")

M = 100 # how many maps to couple, i.e. how much "space" to cover
N = 1000
p = 20
U,V = coupled_henon1D(M, N+p)

# ### Temporal prediction of field U
pool = U[1:N]
test  = U[N:N+p]

γ = 2; τ = 1; B = 1; k = 1;
em = cubic_shell_embedding(pool,γ,τ,B,k,ConstantBoundary(10.))

pred = temporalprediction(pool,em, p; progress = false)

err = [abs.(test[i]-pred[i]) for i=1:p+1]

println("Maximum error: ", maximum(maximum(e) for e in err))

# ### Plot prediction

# Deduce field maximum values
vmax = max(maximum(maximum(s) for s in test),
           maximum(maximum(s) for s in pred))
vmin = min(minimum(minimum(s) for s in test),
           minimum(minimum(s) for s in pred))

# plot plot plot
fig = figure(figsize=(8,8))
ax1 = subplot2grid((3,1), (0,0))
ax2 = subplot2grid((3,1), (1,0))
ax3 = subplot2grid((3,1), (2,0))

ptest = reverse(Matrix(cat(test..., dims = 2)'), dims = 1)
ppred = reverse(Matrix(cat(pred..., dims = 2)'), dims = 1)
perr = reverse(Matrix(cat(err..., dims = 2)'), dims = 1)

im1 = ax1[:imshow](ppred, cmap="viridis", vmin = vmin, vmax = vmax)
im2 = ax2[:imshow](ptest, cmap="viridis", vmin = vmin, vmax = vmax)
im3 = ax3[:imshow](perr, cmap="inferno", vmin = 0, vmax = vmax-vmin)

for (im, ax) in zip([im1,im2,im3], [ax1,ax2,ax3])
    colorbar(im, ax = ax, fraction=0.02, pad=0.01)# format="%.1f")
end
ax1[:set_title]("Prediction")
ax2[:set_title]("Real evolution")
ax3[:set_title]("Absolute error")

ax2[:set_ylabel]("time")
ax3[:set_xlabel]("space")
tight_layout(w_pad=0.1, h_pad=0.00001)


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %%
using Plots
gr()
X = 1100 .+ 200randn(10^6)
histogram(X; norm=true, alpha=0.3, linealpha=0.3, linecolor=:match)

# %%
stephist(X; norm=true)

# %%

max_stars_repo_path: 0023/linecolor and stephist.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0023/linecolor and stephist.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0023/linecolor and stephist.jl

max_forks_repo_name: genkuroki/public


text:
# NOTE: all command entries are restricted to have only one method
# thus we always have at most one doc string.

function read_content(md::Markdown.MD)
    if md.content[1] isa Markdown.MD
        return read_content(md.content[1])
    else
        return md.content
    end
end

read_content(x) = x

"""
    rm_format(md)

Remove Markdown DOM and flatten to strings.
"""
function rm_format end
rm_format(x::Markdown.Paragraph) = join(map(rm_format, x.content))
rm_format(x::Markdown.Code) = x.code
rm_format(x::String) = x
rm_format(x::Markdown.MD) = rm_format(x.content[1])

"""
    read_doc(markdown)

Read CLI documentation from markdown format doc strings.
"""
function read_doc(doc::Markdown.MD)
    has_docstring(doc) || return "", read_args(nothing), read_flags(nothing), read_options(nothing)
    intro = read_intro(doc)
    args = read_args(read_section(doc, "Arguments"))
    flags = read_flags(read_section(doc, "Flags"))
    options = read_options(read_section(doc, "Options"))
    return intro, args, flags, options
end

function read_intro(md::Markdown.MD)
    intro = []
    for line in read_content(md)
        if line isa Markdown.Header{1} && line.text[1] in ["Arguments", "Options", "Flags"]
            break
        else
            push!(intro, line)
        end
    end
    return join(map(rm_format, intro), "\n")
end

function read_section(md::Markdown.MD, title)
    ct = read_content(md)
    nlines = length(ct)
    for k in 1:nlines
        line = ct[k]
        if line isa Markdown.Header{1} && line.text[1] == title
            if k + 1 <= nlines
                return ct[k+1]
            end
        end
    end
    return
end

read_args(::Nothing) = Dict{String,String}()
read_flags(::Nothing) = Dict{String,Tuple{String,Bool}}()
read_options(::Nothing) = Dict{String,Tuple{String,String,Bool}}()

function read_args(md::Markdown.List)
    args = Dict{String,String}()
    for each in md.items
        name, doc = read_item(each[1])
        args[name] = doc
    end
    return args
end

function read_options(md::Markdown.List)
    options = Dict{String,Tuple{String,String,Bool}}()
    # (name, doc, short)
    for each in md.items
        name, doc = read_item(each[1])
        m = match(r"^(-.*) +<(.+)>$", name)

        if m === nothing
            err_m = match(r"^-.*$", strip(name))
            if err_m === nothing
                throw(Meta.ParseError("invalid option: $name"))
            else
                throw(Meta.ParseError("invalid option: $name, expect option argument doc, got only flag"))
            end
        end

        flags = split(m[1], ",")
        if length(flags) == 1
            startswith(name, "--") || throw(Meta.ParseError("invalid option/flag: $name"))
            options[lstrip(m[1], '-')] = (m[2], doc, false)
        elseif length(flags) == 2
            names = strip.(flags)
            name = startswith(names[1], "--") ? lstrip(names[1], '-') :
                startswith(names[2], "--") ? lstrip(names[2], '-') :
                throw(Meta.ParseError("invalid option/flag: $name"))

            options[name] = (m[2], doc, true)
        else
            throw(Meta.ParseError("invalid option/flag syntax: $name"))
        end
    end
    return options
end

function read_flags(md::Markdown.List)
    options = Dict{String,Tuple{String,Bool}}()
    # (doc, short)
    for each in md.items
        name, doc = read_item(each[1])
        flags = split(name, ",")
        if length(flags) == 1
            startswith(name, "--") || throw(Meta.ParseError("invalid flag: $name, might be --$name"))
            options[name[3:end]] = (doc, false)
        elseif length(flags) == 2
            names = strip.(flags)
            name = startswith(names[1], "--") ? lstrip(names[1], '-') :
                startswith(names[2], "--") ? lstrip(names[2], '-') :
                throw(Meta.ParseError("invalid flag: $name"))

            options[name] = (doc, true)
        else
            throw(Meta.ParseError("invalid flag syntax: $name"))
        end
    end
    return options
end

function read_item(raw::Markdown.Paragraph)
    length(raw.content) == 2 || throw(Meta.ParseError("invalid command entry argument doc syntax"))
    raw.content[1] isa Markdown.Code ||
        throw(Meta.ParseError("command argument name should be marked by inline code"))
    name = raw.content[1].code
    doc = read_docstring(raw.content[2])
    return name, doc
end

function read_docstring(doc::String)
    m = match(r"^: *(.*)", strip(doc))
    m === nothing && throw(Meta.ParseError("invalid docstring format: $doc"))
    return String(m[1])
end

function has_docstring(doc::Markdown.MD)
    paragraph = first(read_content(doc))
    flag = paragraph isa Markdown.Paragraph && paragraph.content == Any["No documentation found."]
    return !flag
end


text:
using LinearAlgebra, Krylov, LinearOperators, SparseArrays, Stopping

#Krylov @kdot
macro kdot(n, x, y)
  return esc(:(Krylov.krylov_dot($n, $x, 1, $y, 1)))
end

using SolverTools, Logging
"""
Randomized coordinate descent

Sect. 3.7 in Gower, R. M., & Richtárik, P. (2015).
Randomized iterative methods for linear systems.
SIAM Journal on Matrix Analysis and Applications, 36(4), 1660-1690.
"""
function RandomizedCD(
  A::AbstractMatrix,
  b::AbstractVector{T};
  is_zero_start::Bool = true,
  x0::AbstractVector{T} = zeros(T, size(A, 2)),
  atol::AbstractFloat = 1e-7,
  rtol::AbstractFloat = 1e-15,
  max_iter::Int = size(A, 2)^2,
  max_time::Float64 = 60.0,
  max_cntrs = Stopping._init_max_counters_linear_operators(quick = 20000),
  verbose::Int = 100,
  kwargs...,
) where {T <: AbstractFloat}
  m, n = size(A)
  x = copy(x0)
  res = is_zero_start ? b : b - A * x
  nrm0 = norm(res)

  time_init = time()
  elapsed_time = time_init

  cntrs = LACounters()
  max_f = false

  OK = nrm0 <= atol
  k = 0

  #@info log_header([:iter, :un, :time], [Int, T, T])
  #@info log_row(Any[0, res[1], elapsed_time])
  while !OK && (k <= max_iter) && (elapsed_time - time_init <= max_time) && !max_f

    #rand a number between 1 and n
    #224.662 ns (4 allocations: 79 bytes) - independent of the n
    i = mod(k, n) + 1#Int(floor(rand() * n) + 1)
    Ai = A[:, i]

    #ei = zeros(n); ei[i] = 1.0 #unit vector in R^n
    #xk  = Ai == 0 ? x0 : x0 - dot(Ai,res)/norm(Ai,2)^2 * ei
    Aires = @kdot(m, Ai, res)
    nAi = @kdot(m, Ai, Ai)
    x[i] -= Aires / nAi

    #res = b - A*x
    res += Ai * Aires / nAi
    cntrs.nprod += 1
    nrm = norm(res, Inf)
    OK = nrm <= atol + nrm0 * rtol

    sum, max_f = 0, false
    for f in [:nprod, :ntprod, :nctprod]
      ff = getfield(cntrs, f)
      max_f = max_f || (ff > max_cntrs[f])
      sum += ff
    end
    max_f = max_f || sum > max_cntrs[:neval_sum]
    k += 1
    elapsed_time = time()
    if mod(k, verbose) == 0 #print every 20 iterations
      # @info log_row(Any[k, res[1], elapsed_time])
    end
  end

  return x, OK, k
end


text:
module GLPKInterfaceMIP

import GLPK
importall MathProgBase.SolverInterface
importall ..GLPKInterfaceBase

export
    GLPKSolverMIP,
    GLPKCallbackData,
    optimize!,
    loadproblem!,
    writeproblem,
    getvarLB,
    setvarLB!,
    getvarUB,
    setvarUB!,
    getconstrLB,
    setconstrLB!,
    getconstrUB,
    setconstrUB!,
    getobj,
    setobj!,
    addvar!,
    addconstr!,
    setsense!,
    getsense,
    numvar,
    numconstr,
    setvartype!,
    status,
    getobjval,
    getsolution,
    getconstrsolution,
    getreducedcosts,
    getconstrduals,
    getrawsolver,
    setlazycallback!,
    setcutcallback!,
    setheuristiccallback!,
    cbaddlazy!,
    cbgetlpsolution,
    cbgetmipsolution,
    cbgetstate,
    cbgetobj,
    cbgetbestbound,
    cbgetexplorednodes,
    cbaddcut!,
    cbaddsolution

type GLPKMathProgModelMIP <: GLPKMathProgModel
    inner::GLPK.Prob
    param::GLPK.IntoptParam
    smplxparam::GLPK.SimplexParam
    lazycb::Union{Function,Void}
    cutcb::Union{Function,Void}
    heuristiccb::Union{Function,Void}
    infocb::Union{Function,Void}
    objbound::Float64
    cbdata::MathProgCallbackData
    binaries::Vector{Int}
    userlimit::Bool
    function GLPKMathProgModelMIP()
        lpm = new(GLPK.Prob(), GLPK.IntoptParam(), GLPK.SimplexParam(),
                  nothing, nothing, nothing, nothing, -Inf)
        lpm.cbdata = GLPKCallbackData(lpm)
        lpm.binaries = Int[]
        lpm.userlimit = false
        return lpm
    end
end

type GLPKCallbackData <: MathProgCallbackData
    model::GLPKMathProgModelMIP
    tree::Ptr{Void}
    state::Symbol
    reason::Cint
    sol::Vector{Float64}
    vartype::Vector{Symbol}
    GLPKCallbackData(model::GLPKMathProgModelMIP) = new(model, C_NULL, :Other, -1, Float64[], Char[])
end

type GLPKSolverMIP <: AbstractMathProgSolver
    presolve::Bool
    opts
    GLPKSolverMIP(;presolve::Bool=false, opts...) = new(presolve, opts)
end

callback_abort(stat, tree) = (stat == :Exit && GLPK.ios_terminate(tree))

function _internal_callback(tree::Ptr{Void}, info::Ptr{Void})
    cb_data = unsafe_pointer_to_objref(info)::GLPKCallbackData
    lpm = cb_data.model
    cb_data.tree = tree

    reason = GLPK.ios_reason(tree)
    cb_data.reason = reason

    if reason == GLPK.ISELECT
        #println("reason=SELECT")
        cb_data.state = :Intermediate
    elseif reason == GLPK.IPREPRO
        #println("reason=PREPRO")
        cb_data.state = :MIPNode
    elseif reason == GLPK.IROWGEN
        #println("reason=ROWGEN")
        cb_data.state = :MIPNode
        # if the current solution is actually integer feasible, then
        # return MIPSol status.
        _initsolution!(cb_data)
        cbgetlpsolution(cb_data, cb_data.sol)
        # tol_int = 1e-5 by default
        # TODO: query from GLPK
        all_integer = true
        for i in 1:length(cb_data.sol)
            v::Symbol = cb_data.vartype[i]
            (v == :Int || v == :Bin) || continue
            if abs(cb_data.sol[i]-round(cb_data.sol[i])) > 1e-5
                all_integer = false
                break
            end
        end
        if all_integer
            cb_data.state = :MIPSol
        end
        fill!(cb_data.sol, NaN)

        if lpm.lazycb != nothing
            stat = lpm.lazycb(cb_data)
            callback_abort(stat,tree)
        end
    elseif reason == GLPK.IHEUR
        #println("reason=HEUR")
        cb_data.state = :MIPNode
        if lpm.heuristiccb != nothing
            stat = lpm.heuristiccb(cb_data)
            callback_abort(stat,tree)
        end
    elseif reason == GLPK.ICUTGEN
        #println("reason=CUTGEN")
        cb_data.state = :MIPNode
        if lpm.cutcb != nothing
            stat = lpm.cutcb(cb_data)
            callback_abort(stat,tree)
        end
    elseif reason == GLPK.IBRANCH
        #println("reason=BRANCH")
        cb_data.state = :MIPNode
    elseif reason == GLPK.IBINGO
        #println("reason=BINGO")
        cb_data.state = :MIPSol
    else
        error("internal library error")
    end

    # Doesn't seem like there's a natural "reason" to put this with,
    # so let's just call it everywhere for now
    if lpm.infocb != nothing
        stat = lpm.infocb(cb_data)
        callback_abort(stat,tree)
    end

    bn = GLPK.ios_best_node(tree)
    bn != 0 && (lpm.objbound = GLPK.ios_node_bound(tree, bn))

    return
end

function LinearQuadraticModel(s::GLPKSolverMIP)
    lpm = GLPKMathProgModelMIP()
    lpm.param.msg_lev = GLPK.MSG_ERR
    lpm.smplxparam.msg_lev = GLPK.MSG_ERR
    if s.presolve
        lpm.param.presolve = GLPK.ON
    end

    lpm.param.cb_func = cfunction(_internal_callback, Void, (Ptr{Void}, Ptr{Void}))
    lpm.param.cb_info = pointer_from_objref(lpm.cbdata)

    for (k,v) in s.opts
        if k in [:cb_func, :cb_info]
            warn("ignored option: $(string(k)); use the MathProgBase callback interface instead")
            continue
        end
        i = findfirst(x->x==k, fieldnames(typeof(lpm.param)))
        s = findfirst(x->x==k, fieldnames(typeof(lpm.smplxparam)))
        if !(i > 0 || s > 0)
            warn("Ignored option: $(string(k))")
            continue
        end
        if i > 0
            t = typeof(lpm.param).types[i]
            setfield!(lpm.param, i, convert(t, v))
        end
        if s > 0
            t = typeof(lpm.smplxparam).types[s]
            setfield!(lpm.smplxparam, s, convert(t, v))
        end
    end

    return lpm
end

setlazycallback!(m::GLPKMathProgModel, f::Union{Function,Void}) = (m.lazycb = f)
setcutcallback!(m::GLPKMathProgModel, f::Union{Function,Void}) = (m.cutcb = f)
setheuristiccallback!(m::GLPKMathProgModel, f::Union{Function,Void}) = (m.heuristiccb = f)
setinfocallback!(m::GLPKMathProgModel, f::Union{Function,Void}) = (m.infocb = f)

_check_tree(d::GLPKCallbackData, funcname::AbstractString) =
    (d.tree != C_NULL && d.reason != -1) || error("$funcname can only be called from within a callback")

cbgetstate(d::GLPKCallbackData) = d.state

function cbgetlpsolution(d::GLPKCallbackData, output::Vector)
    _check_tree(d, "cbgetlpsolution")
    lp = GLPK.ios_get_prob(d.tree)
    n = GLPK.get_num_cols(lp)
    length(output) >= n || error("output vector is too short")

    for c = 1:n
        output[c] = GLPK.get_col_prim(lp, c)
    end
    return output
end

function cbgetlpsolution(d::GLPKCallbackData)
    _check_tree(d, "cbgetlpsolution")
    lp = GLPK.ios_get_prob(d.tree)
    n = GLPK.get_num_cols(lp)
    output = Vector{Float64}(n)

    for c = 1:n
        output[c] = GLPK.get_col_prim(lp, c)
    end
    return output
end

function cbgetmipsolution(d::GLPKCallbackData, output::Vector)
    # assuming we're in the lazy callback where
    # the LP solution is actually integral.
    # If we add an informational callback for GLPK.IBINGO,
    # then this will need to be modified.
    return cbgetlpsolution(d, output)
end
cbgetmipsolution(d::GLPKCallbackData) = cbgetlpsolution(d)

function cbgetbestbound(d::GLPKCallbackData)
    _check_tree(d, "cbbestbound")
    lpm = d.model
    return lpm.objbound
end

function cbgetobj(d::GLPKCallbackData)
    _check_tree(d, "cbgetobj")
    lp = GLPK.ios_get_prob(d.tree)
    return GLPK.mip_obj_val(lp)
end

function cbgetexplorednodes(d::GLPKCallbackData)
    _check_tree(d, "cbgetexplorednodes")
    a, _, t = GLPK.ios_tree_size(d.tree)
    return t - a
end

function cbaddlazy!(d::GLPKCallbackData, colidx::Vector, colcoef::Vector, sense::Char, rhs::Real)
    #println("Adding lazy")
    (d.tree != C_NULL && d.reason == GLPK.IROWGEN) ||
        error("cbaddlazy! can only be called from within a lazycallback")
    length(colidx) == length(colcoef) || error("colidx and colcoef have different legths")
    if sense == '='
        bt = GLPK.FX
        rowlb = rhs
        rowub = rhs
    elseif sense == '<'
        bt = GLPK.UP
        rowlb = -Inf
        rowub = rhs
    elseif sense == '>'
        bt = GLPK.LO
        rowlb = rhs
        rowub = Inf
    else
        error("sense must be '=', '<' or '>'")
    end
    # allocating a new vector is not efficient
    solution = cbgetmipsolution(d)
    # if the cut does not exclude the current solution, ignore it
    val = dot(colcoef,solution[colidx])
    if (rowlb - 1e-8 <= val <= rowub + 1e-8)
        # would be better to use GLPK's internal tolerances
        #Base.warn_once("Ignoring lazy constraint which is already satisfied")
        return
    end

    lp = GLPK.ios_get_prob(d.tree)
    GLPK.add_rows(lp, 1)
    m = GLPK.get_num_rows(lp)
    GLPK.set_mat_row(lp, m, colidx, colcoef)
    GLPK.set_row_bnds(lp, m, bt, rowlb, rowub)
    return
end

function cbaddcut!(d::GLPKCallbackData, colidx::Vector, colcoef::Vector, sense::Char, rhs::Real)
    #println("Adding cut")
    (d.tree != C_NULL && d.reason == GLPK.ICUTGEN) ||
        error("cbaddcut! can only be called from within a cutcallback")
    if sense == '<'
        bt = GLPK.UP
    elseif sense == '>'
        bt = GLPK.LO
    elseif sense == '='
        error("unsupported sense in cut plane '='")
    else
        error("sense must be '<' or '>'")
    end
    GLPK.ios_add_row(d.tree, "", 101, colidx, colcoef, bt, rhs)
    return
end

function _initsolution!(d::GLPKCallbackData)
    lp = GLPK.ios_get_prob(d.tree)
    n = GLPK.get_num_cols(lp)
    length(d.sol) == n && return
    resize!(d.sol, n)
    fill!(d.sol, NaN)
    return
end

function _fillsolution!(d::GLPKCallbackData)
    lp = GLPK.ios_get_prob(d.tree)
    n = GLPK.get_num_cols(lp)
    sol = d.sol
    for c = 1:n
        isnan(sol[c]) || continue
        sol[c] = GLPK.mip_col_val(lp, c)
    end
end

function cbaddsolution!(d::GLPKCallbackData)
    #println("Adding sol")
    (d.tree != C_NULL && d.reason == GLPK.IHEUR) ||
        error("cbaddsolution! can only be called from within a heuristiccallback")
    _initsolution!(d)
    _fillsolution!(d)
    # test feasibility of solution, would be better if GLPK supported this
    l = getvarLB(d.model)
    u = getvarUB(d.model)
    for i in 1:length(l)
        if d.sol[i] < l[i] - 1e-6 || d.sol[i] > u[i] + 1e-6
            warn("Ignoring infeasible solution from heuristic callback")
            return
        end
    end
    A = getconstrmatrix(d.model)
    lb = getconstrLB(d.model)
    ub = getconstrUB(d.model)
    y = A*d.sol
    for i in 1:length(lb)
        if y[i] < lb[i] - 1e-6 || y[i] > ub[i] + 1e-6
            warn("Ignoring infeasible solution from heuristic callback")
            return
        end
    end
    GLPK.ios_heur_sol(d.tree, d.sol)
    fill!(d.sol, NaN)
end

function cbsetsolutionvalue!(d::GLPKCallbackData, idx::Integer, val::Real)
    _check_tree(d, "cbsetsolutionvalue!")
    _initsolution!(d)
    d.sol[idx] = val
end

function setsense!(lpm::GLPKMathProgModelMIP, sense)
    lp = lpm.inner
    if sense == :Min
        GLPK.set_obj_dir(lp, GLPK.MIN)
        lpm.objbound = -Inf
    elseif sense == :Max
        GLPK.set_obj_dir(lp, GLPK.MAX)
        lpm.objbound = Inf
    else
        error("unrecognized objective sense: $sense")
    end
end

function setvartype!(lpm::GLPKMathProgModelMIP, vartype::Vector{Symbol})
    lp = lpm.inner
    lpm.binaries = Int[]
    ncol = numvar(lpm)
    @assert length(vartype) == ncol
    for i in 1:ncol
        if vartype[i] == :Int
            coltype = GLPK.IV
        elseif vartype[i] == :Cont
            coltype = GLPK.CV
        elseif vartype[i] == :Bin
            push!(lpm.binaries, i)
            coltype = GLPK.IV
        else
            error("invalid variable type: $(vartype[i])")
        end
        GLPK.set_col_kind(lp, i, coltype)
    end
end

const vartype_map = Dict(
    GLPK.CV => :Cont,
    GLPK.IV => :Int,
    GLPK.BV => :Bin
)

function getvartype(lpm::GLPKMathProgModelMIP)
    lp = lpm.inner
    ncol = numvar(lpm)
    coltype = Array(Symbol, ncol)
    for i in 1:ncol
        ct = GLPK.get_col_kind(lp, i)
        coltype[i] = vartype_map[ct]
        if i in lpm.binaries
            coltype[i] = :Bin
        elseif coltype[i] == :Bin # GLPK said it was binary, but we didn't tell it
            coltype[i] = :Int
        end
    end
    return coltype
end

function optimize!(lpm::GLPKMathProgModelMIP)
    vartype = getvartype(lpm)
    lb = getvarLB(lpm)
    ub = getvarUB(lpm)
    old_lb = copy(lb)
    old_ub = copy(ub)
    for c in 1:length(vartype)
        vartype[c] in [:Int,:Bin] && (lb[c] = ceil(lb[c]); ub[c] = floor(ub[c]))
        vartype[c] == :Bin && (lb[c] = max(lb[c],0.0); ub[c] = min(ub[c],1.0))
    end
    lpm.cbdata.vartype = vartype
    try
        setvarLB!(lpm, lb)
        setvarUB!(lpm, ub)
        if lpm.param.presolve == GLPK.OFF
            ret_ps = GLPK.simplex(lpm.inner, lpm.smplxparam)
            ret_ps != 0 && return ret_ps
        end
        ret = GLPK.intopt(lpm.inner, lpm.param)
        if ret == GLPK.EMIPGAP || ret == GLPK.ETMLIM || ret == GLPK.ESTOP
            lpm.userlimit = true
        end
    finally
        setvarLB!(lpm, old_lb)
        setvarUB!(lpm, old_ub)
    end
end

function status(lpm::GLPKMathProgModelMIP)
    if lpm.userlimit
        return :UserLimit
    end
    s = GLPK.mip_status(lpm.inner)
    if s == GLPK.UNDEF
        if lpm.param.presolve == GLPK.OFF && GLPK.get_status(lpm.inner) == GLPK.NOFEAS
            return :Infeasible
        else
            return :Error
        end
    end
    if s == GLPK.OPT
        return :Optimal
    elseif s == GLPK.INFEAS
        return :Infeasible
    elseif s == GLPK.UNBND
        return :Unbounded
    elseif s == GLPK.FEAS
        return :Feasible
    elseif s == GLPK.NOFEAS
        return :Infeasible
    elseif s == GLPK.UNDEF
        return :Undefined
    else
        error("internal library error")
    end
end

getobjval(lpm::GLPKMathProgModelMIP) = GLPK.mip_obj_val(lpm.inner)

getobjbound(lpm::GLPKMathProgModelMIP) = lpm.objbound

function getsolution(lpm::GLPKMathProgModelMIP)
    lp = lpm.inner
    n = GLPK.get_num_cols(lp)

    x = Array(Float64, n)
    for c = 1:n
        x[c] = GLPK.mip_col_val(lp, c)
    end
    return x
end

function getconstrsolution(lpm::GLPKMathProgModelMIP)
    lp = lpm.inner
    m = GLPK.get_num_rows(lp)

    x = Array(Float64, m)
    for r = 1:m
        x[r] = GLPK.mip_row_val(lp, r)
    end
    return x
end

end


text:
struct MeanRewardBeliefMDP{P<:POMDP, R, S, B, A} <: MDP{B, A}
    pomdp::P
    resample::R
    max_frac_replaced::Float64
    _pm::Vector{S}
    _rm::Vector{Float64}
    _wm::Vector{Float64}
end

function MeanRewardBeliefMDP(pomdp::P, resampler, max_frac_replaced) where {P<:POMDP}
    S = state_type(pomdp)
    b0 = resample(resampler, initial_state_distribution(pomdp), MersenneTwister(14))
    MeanRewardBeliefMDP{P, typeof(resampler), S, typeof(b0), action_type(pomdp)}(
        pomdp,
        resampler,
        max_frac_replaced,
        S[],
        Float64[],
        Float64[]
       )
end

function generate_sr(bmdp::MeanRewardBeliefMDP, b, a, rng::AbstractRNG)
    s = rand(rng, filter(s->!isterminal(bmdp.pomdp, s), particles(b)))
    sp, o = generate_so(bmdp.pomdp, s, a, rng)

    if n_particles(b) > 2*bmdp.resample.n
        b = resample(bmdp.resample, b, rng)
    end

    ps = particles(b)
    pm = bmdp._pm
    rm = bmdp._rm
    wm = bmdp._wm
    resize!(pm, 0)
    resize!(rm, 0)
    resize!(wm, 0)

    all_terminal = true
    for i in 1:n_particles(b)
        s = ps[i]
        if !isterminal(bmdp.pomdp, s)
            all_terminal = false
            sp, r = generate_sr(bmdp.pomdp, s, a, rng)
            push!(pm, sp)
            push!(rm, r)
            od = observation(bmdp.pomdp, s, a, sp)
            push!(wm, pdf(od, o))
        end
    end
    ws = sum(wm)
    if all_terminal || ws < eps(1.0/length(wm))
        # warn("All states in particle collection were terminal.")
        return resample(bmdp.resample, reset_distribution(bmdp.pomdp, b, a, o), rng), 0.0
    end

    pc = resample(bmdp.resample, WeightedParticleBelief{state_type(bmdp.pomdp)}(pm, wm, ws, nothing), rng)
    ps = particles(pc)

    mpw = max_possible_weight(bmdp.pomdp, a, o)
    frac_replaced = bmdp.max_frac_replaced*max(0.0, 1.0 - maximum(wm)/mpw)
    n_replaced = floor(Int, frac_replaced*length(ps))
    is = randperm(rng, length(ps))[1:n_replaced]
    for i in is
        ps[i] = new_particle(bmdp.pomdp, b, a, o, rng)
    end

    return pc, dot(wm, rm)/ws
end

function initial_state(bmdp::MeanRewardBeliefMDP, rng::AbstractRNG)
    return resample(bmdp.resample, initial_state_distribution(bmdp.pomdp), rng)
end

actions(bmdp::MeanRewardBeliefMDP{P,U,B,A}, b::B) where {P,U,B,A} = actions(bmdp.pomdp, b)
actions(bmdp::MeanRewardBeliefMDP) = actions(bmdp.pomdp)

isterminal(bmdp::MeanRewardBeliefMDP, b) = all(isterminal(bmdp.pomdp, s) for s in iterator(b))

discount(bmdp::MeanRewardBeliefMDP) = discount(bmdp.pomdp)


text:
import Base: ==

function ==(sys1::AbstractSystem, sys2::AbstractSystem)
    if typeof(sys1) != typeof(sys2)
        return false
    end
    for field in fieldnames(typeof(sys1))
        if getfield(sys1, field) != getfield(sys2, field)
            return false
        end
    end
    return true
end

function Base.isapprox(sys1::AbstractSystem, sys2::AbstractSystem; kwargs...)
    if typeof(sys1) != typeof(sys2)
        return false
    end
    for field in fieldnames(typeof(sys1))
        if !isapprox(getfield(sys1, field), getfield(sys2, field); kwargs...)
            return false
        end
    end
    return true
end


text:
using Images, Colors, GeometryTypes
using Reactive, FileIO, GLVisualize, Compat
using GLAbstraction, GeometryTypes, GLWindow

if !isdefined(:runtests)
	window = glscreen()
end
const record_interactive = true

# loadasset loads data from the GLVisualize asset folder and is defined as
# FileIO.load(assetpath(name))
doge = loadasset("doge.png")
# Convert to RGBA{Float32}. Float for filtering and 32 because it fits the GPU better
img = map(RGBA{Float32}, doge)
# create a slider that goes from 1-20 in 0.1 steps
slider_s, slider = widget(Signal(1f0), range=1f0:0.1f0:20f0, window)

# performant conversion to RGBAU8, implemted with a functor
# in 0.5 anonymous functions offer the same speed, so this wouldn't be needed
immutable ClampRGBAU8 end
@compat (::ClampRGBAU8)(x) = RGBA{U8}(clamp(comp1(x), 0,1), clamp(comp2(x), 0,1), clamp(comp3(x), 0,1), clamp(alpha(x), 0,1))

"""
Applies a gaussian filter to `img` and converts it to RGBA{U8}
"""
function myfilter(img, sigma)
	img = Images.imfilter_gaussian(img, [sigma, sigma])
	map(ClampRGBAU8(), img).data
end

startvalue = myfilter(img, value(slider_s))
# Use Reactive.async_map, to filter the image without blocking the main process
task, imgsig = async_map(myfilter, startvalue, Signal(img), slider_s)
# visualize the image signal
image_renderable = visualize(
    imgsig,
    model=translationmatrix(Vec3f0(50,100,0)),
    is_fully_opaque=false
)
_view(image_renderable)

vec2i(a,b,x...) = Vec{2,Int}(round(Int, a), round(Int, b))
vec2i(vec::Vec) = vec2i(vec...)
"""
creates a rectangle around `robj`
"""
function screen(robj)
	bb = value(boundingbox(robj))
	m  = vec2i(minimum(bb))
	area = SimpleRectangle{Float32}(0,0, ((vec2i(maximum(bb))-m)+30)...)

	_view(visualize((area, [Point2f0(0)]),
        color=RGBA{Float32}(0,0,0,0), stroke_color=RGBA{Float32}(0,0,0,0.7),
        stroke_width=3f0),
        camera=:fixed_pixel
    )
	robj.children[][:model] = translationmatrix(Vec3f0(15,15,0)-minimum(bb))
	_view(robj, camera=:fixed_pixel)
end
screen(slider)

if !isdefined(:runtests)
	renderloop(window)
end


text:
export HamiltonianSplittingBoris

"""
    HamiltonianSplittingBoris( maxwell_solver,
                               kernel_smoother_0, kernel_smoother_1,
                               particle_group,
                               e_dofs_1, e_dofs_2, b_dofs)

Boris pusher in GEMPIC framework (spline finite elements)

- `mid` describing the magnetic field at time ``t_{n+1/2}`` (used for push)
- `j_dofs` for kernel representation of current density.
- `maxwell_solver`    : Maxwell solver
- `kernel_smoother_0` : Kernel smoother
- `kernel_smoother_1` : Kernel smoother
- `particle_group`    : Particle group
- `efield_dofs`       : array for the coefficients of the efields
- `bfield_dofs`       : array for the coefficients of the bfield
- `x_min`             : Lower bound of x domain
- `Lx`                : Length of the domain in x direction.

"""
struct HamiltonianSplittingBoris <: AbstractSplitting
    maxwell_solver::AbstractMaxwellSolver
    kernel_smoother_0::ParticleMeshCoupling1D #(order p+1)
    kernel_smoother_1::ParticleMeshCoupling1D #(order p)
    particle_group::ParticleGroup
    spline_degree::Int64
    Lx::Float64
    x_min::Float64 # Lower bound for x domain
    delta_x::Float64 # Grid spacing

    cell_integrals_0::SVector
    cell_integrals_1::SVector

    e_dofs::Array{Array{Float64,1},1}
    e_dofs_mid::Array{Array{Float64,1},1}
    b_dofs::Array{Float64,1}
    b_dofs_mid::Array{Float64,1}
    j_dofs::Array{Array{Float64,1}}

    function HamiltonianSplittingBoris(
        maxwell_solver::Maxwell1DFEM,
        kernel_smoother_0::ParticleMeshCoupling1D,
        kernel_smoother_1::ParticleMeshCoupling1D,
        particle_group::ParticleGroup,
        e_dofs::Array{Vector{Float64},1},
        b_dofs::Vector{Float64},
    )
        e_dofs_mid = [
            zeros(Float64, kernel_smoother_1.n_dofs),
            zeros(Float64, kernel_smoother_1.n_dofs),
        ]

        j_dofs = [
            zeros(Float64, kernel_smoother_0.n_dofs),
            zeros(Float64, kernel_smoother_0.n_dofs),
        ]

        b_dofs_mid = zeros(Float64, kernel_smoother_1.n_dofs)

        x_min = maxwell_solver.xmin
        Lx = maxwell_solver.Lx
        spline_degree = 3
        delta_x = Lx / kernel_smoother_1.n_dofs

        cell_integrals_1 = SVector{3}([0.5, 2.0, 0.5] ./ 3.0)
        cell_integrals_0 = SVector{4}([1.0, 11.0, 11.0, 1.0] ./ 24.0)

        return new(
            maxwell_solver,
            kernel_smoother_0,
            kernel_smoother_1,
            particle_group,
            spline_degree,
            Lx,
            x_min,
            delta_x,
            cell_integrals_0,
            cell_integrals_1,
            e_dofs,
            e_dofs_mid,
            b_dofs,
            b_dofs_mid,
            j_dofs,
        )
    end
end

export staggering!
"""
    staggering(splitting, dt)

Propagate ``E_0`` to ``E_{1/2}`` and ``x_0`` to ``x_{1/2}`` to initialize
the staggering
- `splitting` : time splitting object
- `dt`   : time step
"""
function staggering!(splitting, dt)
    push_x_accumulate_j!(splitting, dt * 0.5)

    # (4) Compute E_{n+1}
    splitting.e_dofs_mid[1] .= splitting.e_dofs[1]
    splitting.e_dofs_mid[2] .= splitting.e_dofs[2]
    splitting.j_dofs[1] .= 0.5dt .* splitting.j_dofs[1]
    splitting.j_dofs[2] .= 0.5dt .* splitting.j_dofs[2]
    # Ex part:

    compute_e_from_j!(
        splitting.e_dofs_mid[1], splitting.maxwell_solver, splitting.j_dofs[1], 1
    )

    #todo "Combine the two steps for efficiency"

    compute_e_from_b!(
        splitting.e_dofs_mid[2], splitting.maxwell_solver, 0.5dt, splitting.b_dofs
    )

    return compute_e_from_j!(
        splitting.e_dofs_mid[2], splitting.maxwell_solver, splitting.j_dofs[2], 2
    )
end

"""
    strang_splitting(splitting, dt, number_steps)

Second order Boris pusher using staggered grid
- `splitting` : time splitting object
- `dt`   : time step
- `number_steps` : number of time steps
"""
function strang_splitting!(
    splitting::HamiltonianSplittingBoris, dt::Float64, number_steps::Int64
)
    for i_step in 1:number_steps

        # (1) Compute B_{n+1/2} from B_n
        splitting.b_dofs_mid .= splitting.b_dofs

        compute_b_from_e!(
            splitting.b_dofs, splitting.maxwell_solver, dt, splitting.e_dofs_mid[2]
        )

        splitting.b_dofs_mid .= (splitting.b_dofs_mid .+ splitting.b_dofs) * 0.5

        # (2) Propagate v: v_{n-1/2} -> v_{n+1/2}
        # (2a) Half time step with E-part
        push_v_epart!(splitting, 0.5dt)
        # (2b) Full time step with B-part
        push_v_bpart!(splitting, dt)
        # (2c) Half time step with E-part
        push_v_epart!(splitting, 0.5dt)

        # (3) Propagate x: x_n -> x_{n+1}. Includes also accumulation of j_x, j_y
        push_x_accumulate_j!(splitting, dt)

        # (4) Compute E_{n+1}
        splitting.e_dofs[1] .= splitting.e_dofs_mid[1]
        splitting.e_dofs[2] .= splitting.e_dofs_mid[2]
        splitting.j_dofs[1] .= dt * splitting.j_dofs[1]
        splitting.j_dofs[2] .= dt * splitting.j_dofs[2]

        # Ex part:
        compute_e_from_j!(
            splitting.e_dofs_mid[1], splitting.maxwell_solver, splitting.j_dofs[1], 1
        )

        # Ey part:
        compute_e_from_b!(
            splitting.e_dofs_mid[2], splitting.maxwell_solver, dt, splitting.b_dofs
        )

        compute_e_from_j!(
            splitting.e_dofs_mid[2], splitting.maxwell_solver, splitting.j_dofs[2], 2
        )
    end
end

"""
    push_v_epart(splitting, dt)

Pusher for ``E \\nabla_v `` part

```math
V_{new} = V_{old} + dt ⋅ E
```
"""
function push_v_epart!(splitting, dt)
    qm = splitting.particle_group.q_over_m

    for i_part in 1:(splitting.particle_group.n_particles)
        xi = get_x(splitting.particle_group, i_part)

        efield_1 = evaluate(splitting.kernel_smoother_1, xi[1], splitting.e_dofs_mid[1])

        efield_2 = evaluate(splitting.kernel_smoother_0, xi[1], splitting.e_dofs_mid[2])

        v_new = get_v(splitting.particle_group, i_part)
        v_new .= v_new .+ dt * qm .* [efield_1, efield_2]

        set_v!(splitting.particle_group, i_part, v_new)
    end
end

"""
    push_v_bpart!(splitting, dt)

  Pusher for vxB part
"""
function push_v_bpart!(splitting::HamiltonianSplittingBoris, dt::Float64)
    qmdt = splitting.particle_group.q_over_m * 0.5 * dt
    v_new = zeros(3)

    for i_part in 1:(splitting.particle_group.n_particles)
        vi = get_v(splitting.particle_group, i_part)
        xi = get_x(splitting.particle_group, i_part)

        bfield = evaluate(splitting.kernel_smoother_1, xi[1], splitting.b_dofs_mid)

        bfield = qmdt * bfield

        M11 = 1.0 / (1.0 + bfield^2)
        M12 = M11 * bfield * 2.0
        M11 = M11 * (1 - bfield^2)

        v_new[1] = M11 * vi[1] + M12 * vi[2]
        v_new[2] = -M12 * vi[1] + M11 * vi[2]
        v_new[3] = 0.0

        set_v!(splitting.particle_group, i_part, v_new)
    end
end

"""
    push_x_accumulate_j!(splitting, dt)

Pusher for x and accumulate current densities

For each particle compute the index of the first DoF on the grid it
contributes to and its position (normalized to cell size one).

Then update particle position:  ``X_{new} = X_{old} + dt ⋅ V``

!!! note
    `j_dofs` does not hold the values for j itself but
    for the integrated j.

"""
function push_x_accumulate_j!(splitting, dt)
    n_cells = splitting.kernel_smoother_0.n_dofs

    fill!(splitting.j_dofs[1], 0.0)
    fill!(splitting.j_dofs[2], 0.0)

    for i_part in 1:(splitting.particle_group.n_particles)

        # Read out particle position and velocity
        x_old = get_x(splitting.particle_group, i_part)
        vi = get_v(splitting.particle_group, i_part)

        # Then update particle position:  X_new = X_old + dt * V
        x_new = x_old .+ dt * vi

        # Get charge for accumulation of j
        wi = get_charge(splitting.particle_group, i_part)
        qoverm = splitting.particle_group.q_over_m

        #todo "Check here also second posibility with sum of two accumulations"
        # Accumulate jx
        add_charge!(
            splitting.j_dofs[1],
            splitting.kernel_smoother_1,
            (x_old[1] + x_new[1]) * 0.5,
            wi[1] * vi[1],
        )
        # Accumulate jy
        add_charge!(
            splitting.j_dofs[2],
            splitting.kernel_smoother_0,
            (x_old[1] + x_new[1]) * 0.5,
            wi[1] * vi[2],
        )

        x_new[1] = mod(x_new[1], splitting.Lx)
        set_x!(splitting.particle_group, i_part, x_new)
    end
end


text:
b₀_of(problem) =
    σ₀_of(problem.physic) * ρ_of(problem.atoms) * geometric_integral(problem.atoms)

σ₀_of(problem::Scalar) = 4π / k₀^2
σ₀_of(problem::Vectorial) = 6π / k₀^2

"""
    ρ_of(atoms::Shape{T})
Computes Volumetric Density
"""
ρ_of(atoms::Atom{T}) where {T} = atoms.N / Volume_of(atoms)

"""
    Volume_of(atoms::Shape{Circle})
"""
Volume_of(atoms::Atom{Circle}) = π * atoms.sizes^2 # sizes == circle radius

"""
    Volume_of(atoms::Shape{Cube})
"""
Volume_of(atoms::Atom{Cube}) = atoms.sizes^3 # sizes == length of Cube's Sides

"""
    Volume_of(atoms::Shape{Cylinder})
"""
Volume_of(atoms::Atom{Cylinder}) = atoms.sizes[:h] * π * atoms.sizes[:R]^2

"""
    Volume_of(atoms::Shape{Sphere})
"""
Volume_of(atoms::Atom{Sphere}) = (4π / 3) * atoms.sizes^3 # sizes == radius of Sphere

#=
These "geometric_integral" values holds only for homogeneous distributions
- They represent the light path acroos the cloud -
For example, to cross a sphere, we need the diameter of the sphere (diameter = 2*radius)
=#
geometric_integral(atoms::Atom{Sphere}) = 2 * atoms.sizes # diameter = 2*radius
geometric_integral(atoms::Atom{Cube}) = atoms.sizes # diameter = cube side = L

#=
(usually) we want to know the physics in the laser direction, that is, at z-direction,
thus, we want to know the light path across the height.
=#
function geometric_integral(atoms::Atom{Cylinder})
    @info "using the cylinder height as light path"
    return atoms.sizes[:h]
end

"""
    estimate_E₀(laser)
"""
estimate_E₀(laser) = √(laser.s * (1 + 4(laser.Δ / Γ)^2) / 2)

"""
    sph2cart([θ, ϕ, r])

    θ = "elevation" or "Polar" = projection on Z-axis (in radians)
    ϕ = azimuth = projection on XY-plane (in radians)
    r = radius

	ref: https://doc.sagemath.org/html/en/thematic_tutorials/vector_calculus/vector_calc_change.html
"""
function sph2cart(spherical_coordinate)
    θ = spherical_coordinate[1]
    ϕ = spherical_coordinate[2]
    r = spherical_coordinate[3]
    x = r * cos(ϕ) * sin(θ)
    y = r * sin(ϕ) * sin(θ)
    z = r * cos(θ)
    return [x, y, z]
end
function cart2sph(cartesian_coordinate)
    x = cartesian_coordinate[1]
    y = cartesian_coordinate[2]
    z = cartesian_coordinate[3]
    ϕ = atan(y, x)
    θ = atan(sqrt(x^2 + y^2), z)
    r = sqrt(x^2 + y^2 + z^2)
    return [θ, ϕ, r]
end


text:
using Test, SimpleHypergraphs
using Random
import LightGraphs

h1 = Hypergraph{Float64}(5,4)
h1[1:3,1] .= 1.5
h1[3,4] = 2.5
h1[2,3] = 3.5
h1[4,3:4] .= 4.5
h1[5,4] = 5.5
h1[5,2] = 6.5

@testset "SimpleHypergraphs Hypergraph      " begin

    h = hg_load("data/test1.hgf", Int)
    @test size(h) == (4, 4)
    @test nhv(h) == 4
    @test nhe(h) == 4
    m = Matrix(h)
    @test m == h
    @test h == [1       nothing 4       nothing
                2       3       nothing nothing
                nothing nothing 5       nothing
                nothing nothing 6       nothing]
    mktemp("data") do path, _
        println(path)
        hg_save(path, h)

        loaded_hg = replace(read(path, String), r"\n*$" => "")

        @test loaded_hg ==
            reduce(replace,
                ["\r\n"=>"\n",
                r"^\"\"\"(?s).*\"\"\"\n"=>"", #remove initial comments
                r"\n*$"=>""], #remove final \n*
                init=read("data/test1.hgf", String)) #no comments

        @test loaded_hg ==
            reduce(replace,
                ["\r\n"=>"\n",
                r"^\"\"\"(?s).*\"\"\"\n"=>"", #remove initial comments
                r"\n*$"=>""], #remove final \n*
                init=read("data/test_singlelinecomment.hgf", String)) #single line comment

        @test loaded_hg ==
            reduce(replace,
                ["\r\n"=>"\n",
                r"^\"\"\"(?s).*\"\"\"\n"=>"", #remove initial comments
                r"\n*$"=>""], #remove final \n*
                init=read("data/test_multiplelinescomment.hgf", String)) #multiple lines comment
    end

    @test_throws ArgumentError hg_load("data/test_malformedcomment.hgf", Int)

    h2 = Hypergraph{Float64}(0,0)
    for i in 1:4 add_vertex!(h2) end
    add_hyperedge!(h2;vertices=Dict(1:3 .=> 1.5))
    add_hyperedge!(h2)
    add_vertex!(h2;hyperedges=Dict(2=>6.5))
    add_hyperedge!(h2;vertices=Dict(2 => 3.5, 4 => 4.5))
    add_hyperedge!(h2;vertices=Dict(3:5 .=> (2.5,4.5,5.5)))
    @test h1 == h2
    m = Matrix(h1)
    @test  m == Matrix(h2)
    @test h1 == Hypergraph(m)
    @test h1 == Hypergraph{Nothing}(m)
    @test h1 == Hypergraph{Nothing, Nothing}(m)
    @test getindex(h1,3,1) == 1.5

    h3 = Hypergraph{Float64,String,Nothing}(1,1)
    @test add_vertex!(h3;v_meta="test") == 2
    @test set_vertex_meta!(h3,"t",1) == ["t","test"]
    @test get_vertex_meta(h3,2) == "test"
    @test get_hyperedge_meta(h3,1) == nothing
    @test_throws BoundsError get_hyperedge_meta(h3,2)

    h4 = Hypergraph{Float64,Nothing,String}(1,1)
    @test add_hyperedge!(h4;he_meta="test") == 2
    @test set_hyperedge_meta!(h4,"t",1) == ["t","test"]
    @test get_hyperedge_meta(h4,2) == "test"
    @test get_vertex_meta(h4,1) == nothing
    @test_throws BoundsError get_vertex_meta(h4,2)

    h1_0 = deepcopy(h1)
    @test add_vertex!(h1_0) == 6
    h1_0[6,:] = h1_0[5,:]
    @test remove_vertex!(h1_0,5) == h1

end;

@testset "SimpleHypergraphs BipartiteView   " begin
    h2 = deepcopy(h1)

    @test LightGraphs.nv(LightGraphs.zero(BipartiteView{Int})) == 0

    b = BipartiteView(h2)
    @test LightGraphs.edgetype(b) == LightGraphs.SimpleGraphs.SimpleEdge{Int}
    @test LightGraphs.has_vertex(b, 0) == false
    @test LightGraphs.has_vertex(b, 1) == true
    @test LightGraphs.nv(LightGraphs.zero(b)) == 0

    @test LightGraphs.is_directed(b) == false
    @test LightGraphs.is_directed(typeof(b)) == false
    @test LightGraphs.eltype(b) == Int

    @test sum(LightGraphs.adjacency_matrix(LightGraphs.SimpleGraph(b))) == 18

    @test sort(collect(LightGraphs.outneighbors(b,5))) == [7,9]
    @test sort(collect(LightGraphs.outneighbors(b,1))) == [6]
    @test sort(collect(LightGraphs.inneighbors(b,9))) == [3,4,5]

    @test Set(LightGraphs.vertices(b)) == Set(1:LightGraphs.nv(b))

    @test shortest_path(b,1,5) == [1,3,5]
    @test LightGraphs.is_weakly_connected(b) == true

    @test add_vertex!(h2) == 6
    @test add_hyperedge!(h2) == 5
    h2[5,5] = 1
    h2[6,5] = 1

    @test shortest_path(b,1,6) == [1,3,5,6]

    bipartite_graph = LightGraphs.SimpleGraph(b)

    @test LightGraphs.SimpleGraphs.fadj(bipartite_graph)==LightGraphs.SimpleGraphs.fadj(b)
    @test LightGraphs.nv(b) == 11
    @test LightGraphs.ne(b) == 11

    @test sort!(LightGraphs.SimpleGraphs.fadj(b,1)) == [7]
    @test sort!(LightGraphs.SimpleGraphs.fadj(b,2)) == [7,9]
end;

@testset "SimpleHypergraphs TwoSectionView  " begin

    ht = Hypergraph{Float64}(3,3)
    ht[1:2,1:2] .= 2.
    ht[:, 3] .= 2.

    add_vertex!(h1)
    add_hyperedge!(h1)
    h1[5,5] = 1
    h1[6,5] = 1

    @test LightGraphs.nv(LightGraphs.zero(TwoSectionView{Int})) == 0

    t = TwoSectionView(h1)
    @test LightGraphs.edgetype(t) == LightGraphs.SimpleGraphs.SimpleEdge{Int}
    @test LightGraphs.has_vertex(t, 0) == false
    @test LightGraphs.has_vertex(t, 1) == true
    @test LightGraphs.nv(LightGraphs.zero(t)) == 0

    @test LightGraphs.is_directed(t) == false
    @test LightGraphs.is_directed(typeof(t)) == false
    @test LightGraphs.eltype(t) == Int

    @test LightGraphs.nv(t) == 6
    @test LightGraphs.ne(t) == 8

    @test sort(LightGraphs.all_neighbors(t, 1)) == [2,3]
    @test sort(LightGraphs.outneighbors(t, 5)) == [3,4,6]
    @test sort(LightGraphs.inneighbors(t, 4)) == [2,3,5]
    @inferred LightGraphs.all_neighbors(t, 1)

    @test LightGraphs.has_edge(t, 1, 2) == true
    @test LightGraphs.has_edge(t, 1, 5) == false

    @test sum(LightGraphs.adjacency_matrix(LightGraphs.SimpleGraph(t))) == 16
    @test shortest_path(t,1,5) == [1,3,5]
    @test LightGraphs.is_weakly_connected(t) == true

    @test SimpleHypergraphs.add_vertex!(h1) == 7
    h1[7,5] = 1

    @test shortest_path(t,1,6) == [1,3,5,6]

    @test LightGraphs.nv(t) == 7
    @test LightGraphs.ne(t) == 10
    @test sort(LightGraphs.outneighbors(t, 5)) == [3,4,6,7]

    @test sum(LightGraphs.adjacency_matrix(LightGraphs.SimpleGraph(t))) == 20

    Random.seed!(0);
    g = LightGraphs.erdos_renyi(8, 0.3)
    h_from_g = Hypergraph(g)
    @test LightGraphs.adjacency_matrix(g) == LightGraphs.adjacency_matrix(TwoSectionView(h_from_g))
    @test minimum([sum((h_from_g .== true)[:,n]) for n in 1:6] .== 2)
    @test LightGraphs.modularity(g,[1,1,2,2,3,3,4,4]) ≈ modularity(h_from_g, Set.([[1,2],[3,4],[5,6],[7,8]]))
    @test LightGraphs.SimpleGraphs.fadj(g) == LightGraphs.SimpleGraphs.fadj(TwoSectionView(h_from_g))
end;

@testset "SimpleHypergraphs Modularity      " begin
    Random.seed!(1234);
    hg = Hypergraph{Bool}(10, 12)
    for i in eachindex(hg)
        if rand() < 0.2
            hg[i] = true
        end
    end

    cfmr = CFModularityRandom(3,10000)

    @test findcommunities(hg,cfmr) ==
         (bp = Set.([[4, 5, 9], [1, 3, 6, 7], [2, 8, 10]]), bm = 0.21505688117829677)
    @test modularity(hg,  Set.([1:10])) == 0.0
    Random.seed!(1234);
    @test randompartition(hg, 2) == Set.([[1, 5, 6, 7, 9], [2, 3, 4, 8, 10]])

    hh = Hypergraph{Bool}(7,4)
    hh[1,1] = true
    hh[2,1:2] .= true
    hh[3,1:3] .= true
    hh[4,4] = true
    hh[5:6,3] .= true
    @test nhv(hh) == 7
    @test nhe(hh) == 4

    @test modularity(hh,Set.([[1,2,3],[4],[5],[6],[7]])) ≈ 223/972
    @test modularity(hh,Set.([[1,2,3],[4,5,6,7]])) ≈ 14/72
    @test modularity(hh,Set.([[1,2,3,5,6],[4,7]])) ≈ 16/81
    @test modularity(hh,Set.([[1,2,3,5,6],[4],[7]])) ≈ 16/81
    @test modularity(hh, Set.([1:nhv(hh)])) == 0.0
    ha = SimpleHypergraphs.HypergraphAggs(hh)
    @test ha.hes == [3, 2, 3, 1]
    @test ha.max_hes == 3
    @test ha.deg_vs == [1, 2, 3, 1, 1, 1, 0]
    @test ha.volV == 9
    @test modularity(hh, Set.([[1,2,3],[4],[5],[6],[7]]), ha) ≈ 223/972
    cfmr = CFModularityRandom(2,10000)
    @test cfmr.n==2
    @test cfmr.reps == 10000
    @test findcommunities(hh,cfmr).bm ≈ 16/81
    Random.seed!(1234);
    cnm = CFModularityCNMLike(100)
    @test cnm.reps == 100
    @test findcommunities(hh, CFModularityRandom(4,10000)).bm ≈ findcommunities(hh, cnm).bm
    Random.seed!(0);
    @test findcommunities(hh, cnm).bm ≈ 223/972

end;
                                     #
@testset "SimpleHypergraphs randomized tests" begin
    Random.seed!(0)
    N = 100
    res = Vector{Bool}(undef, N)
    for i in 1:N
        m1 = CFModularityCNMLike(100)
        m2 = CFModularityRandom(2,100)
        r = rand([repeat([nothing],6)..., true], 12, 8)
        hh = Hypergraph(r)
        bm1 = findcommunities(hh, m1).bm
        bm2 = findcommunities(hh, m2).bm
        res[i] = (bm1 > bm2)
    end
    @test sum(res) >= N*0.80
end


text:
using Revise
using
    BeamPropagation, StaticArrays, DelimitedFiles, BenchmarkTools, Plots, StructArrays, StatsBase

# Define constants
const h = @with_unit 1 "h"
const ħ = @with_unit 1 "ħ"
const μ = @with_unit 1 "μB"
const M = @with_unit 190 "u"
const λ = @with_unit 577 "nm"
const Γ = @with_unit 2π * 8.3e6 "MHz"
const Iₛ = @with_unit 5.0 "mW/cm^2"
const w = @with_unit 0.5 "cm"
const P = @with_unit 250.0 "mW"

# Define parameters
n           = 1000000
dt          = 5e-7
max_steps   = 1e5

# Generate initial positions and velocities
r_transform(r) = @. 2.5e-3 * (2r - 1)
rx = r_transform.(rand(Float64, n))
ry = r_transform.(rand(Float64, n))
rz = zeros(n)

v_transform(v) = @. 4.0 * (2v - 1)
vx = v_transform.(rand(Float64, n))
vy = v_transform.(rand(Float64, n))
vz = randn(n) .* (30 / 2.35) .+ 70

rs = SVector.(rx, ry, rz)
vs = SVector.(vx, vy, vz)
as = SVector.(zeros(n), zeros(n), zeros(n))
gs = rand(1:12, n)

# Define a function `magnetic_force` that maps particle → force (based on a combination or `v` and `v`)
Barr = readdlm("examples/Bfieldmap.txt")
xlut   = Barr[:,1];         ylut   = Barr[:,2];         zlut   = Barr[:,3];
Bxlut  = Barr[:,4];         Bylut  = Barr[:,5];         Bzlut  = Barr[:,6];
dBxlut = Barr[:,7] ./ 100;  dBylut = Barr[:,8] ./ 100;  dBzlut = Barr[:,9] ./ 100;

dFlut = Array{SVector{3, Float64}}(undef, 51, 51, 201)
Bnormlut = zeros(Float64, (51, 51, 201))
for (i, (x, y, z)) in enumerate(zip(xlut, ylut, zlut))

    ix = round(Int, 1e4 * x + 25) + 1
    iy = round(Int, 1e4 * y + 25) + 1
    iz = mod(round(Int, z / 20e-5), 201) + 1

    dFlut[ix, iy, iz] = SVector(dBxlut[i], dBylut[i], dBzlut[i]) * (μ / M)
    Bnormlut[ix, iy, iz] = sqrt(Bxlut[i]^2 + Bylut[i]^2 + Bzlut[i]^2)
end

@inline function is_dead(r)
    x, y, z = r
    return (x^2 + y^2 > 5.76e-6) | (z > 1.2)
end

num_gs = 12; num_es = 4;

TItab = [0 0.0003 0 0.0833 0 0.1656 0 0.0833 0 0.0008 0 0;
         0.0835 0 0.0830 0 0.0003 0 0.0007 0 0 0 0.1657 0;
         0 0.1662 0 0.0007 0 0.0004 0 0 0 0.0827 0 0.0833;
         0 0 0.0007 0 0.1662 0 0.0829 0 0.0831 0 0.0004 0]

BRtab = [0 0.001 0 0.167 0 0.331 0 0.167 0.333 0.002 0 0;
         0.167 0 0.166 0 0.001 0.002 0.001 0 0 0.332 0.331 0;
         0.333 0.332 0 0.001 0 0.001 0 0 0 0.165 0 0.167;
         0 0.001 0.001 0.332 0.332 0 0.166 0 0.166 0 0.001 0]

const g_to_e_TotalProb = sum(TItab, dims=1)
const end_gs = collect(1:n)
const μ_signs = SVector(-1, -1, -1, -1, -1, 1, 1, -1, 1, 1, 1, 1)
const wvs = Array{Weights{Float64,Float64,Array{Float64,1}}}(undef, num_gs)

TIBRtab = zeros(Float64, num_gs, num_gs);
for g in 1:num_gs, g′ in 1:num_gs
    for e in 1:num_es
        TIBRtab[g, g′] += TItab[e, g] * BRtab[e, g′] # Probability of gstate -> estate -> gstate′
    end
end

for g in 1:num_gs
    wvs[g] = Weights(TIBRtab[g, :])
end

δs = [-13.75e9, 2.5e9];
p = @params (dFlut, Bnormlut, δs)

@inline Rsc(s, δ) =  (Γ / 2) * (s / (1 + s + 4(δ/Γ)^2))
@inline detuning(Δ0, UB, vz) = 2π * (Δ0 + UB / h + vz / λ)
@inline Prob(g, s, δ, dt) = Rsc(s, δ) * dt * g_to_e_TotalProb[g]

@inline function scatterphoton(g, s, δ, dt)::Tuple{Bool, Int64}
    if rand() > Prob(g, s, δ, dt)
        return false, g
    else
        g′ = sample(end_gs, wvs[g])
        return true, g′
    end
end

@inline function f(r, v, a, g, dt, p)
    x, y, z = r

    # Compute "lookup" indices corresponding to the current particle position
    # These are used to index into any arrays utilized -- here `Bnormlut` and `dFlut`
    ix = round(Int, 1e4 * x + 25) + 1
    iy = round(Int, 1e4 * y + 25) + 1
    iz = mod(round(Int, z / 20e-5), 201) + 1

    P0 = 2P / (Iₛ * π * w^2)
    s = P0 * exp(-(x^2 + y^2) / (2w^2))
    UB = μ * μ_signs[g] * p.Bnormlut[ix, iy, iz]
    d1 = detuning(p.δs[1], UB, v[3])
    scat1, g′  = scatterphoton(g, s, d1, dt)
    d2 = detuning(p.δs[2], UB, v[3])
    scat2, g′′ = scatterphoton(g′, s, d2, dt)
    a′ = μ_signs[g] * p.dFlut[ix, iy, iz]

    return (a′, g′′)
end

η = 35; ξ = 100;

@btime propagate!(rs_, vs_, as_, gs_, $f, $is_dead, $η, $ξ, $dt, $max_steps, $p) setup=(rs_=deepcopy($rs); vs_=deepcopy($vs); as_=deepcopy($as); gs_=deepcopy($gs)) evals=1;

rs_=deepcopy(rs); vs_=deepcopy(vs); as_=deepcopy(as); gs_=deepcopy(gs);
rs_traj, vs_traj = propagate!(rs_, vs_, as_, gs_, f, is_dead, η, ξ, dt, max_steps, p);

@inline function f_nothing(r, v, a, g, dt, p)
    return (a, g)
end

rs_=deepcopy(rs); vs_=deepcopy(vs); as_=deepcopy(as); gs_=deepcopy(gs);
rs_traj, vs_traj = propagate!(rs_, vs_, as_, gs_, f_nothing, is_dead, η, ξ, dt, max_steps, p);

# plot()
# anim = @animate for i in 1:100
#     r = rs_traj[i]
#     scatter!(
#         [r_[3] for r_ in r],
#         [r_[1] for r_ in r],
#         [r_[2] for r_ in r],
#         xlims=(0, 1.2),
#         ylims=(-0.005, 0.005),
#         zlims=(-0.005, 0.005),
#         alpha=0.30,
#         legend=false,
#         markersize=2
#         )
# end
# gif(anim, "some_file_name.gif", fps=10)
#
# rz_final = [r[end][3] for r in rs_traj]
# vz_final = [v[end][3] for v in vs_traj]
# idxs = rz_final .> 1.14
# histogram(vz_final[idxs])