text:
# # Surface-grid operations

#md # ```@meta
#md # CurrentModule = ImmersedLayers
#md # ```

#=
Here, we will discuss the various surface-grid operators available in the package.
We will start by generating the cache, just as we did in [Immersed layer caches](@ref)
=#

using ImmersedLayers
using Plots
using LinearAlgebra

#=
### Set up the grid, shape, and cache
We do this just as we did in [Immersed layer caches](@ref)
=#
Δx = 0.01
Lx = 4.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx)
RadC = 1.0
Δs = 1.4*cellsize(g)
body = Circle(RadC,Δs)
cache = SurfaceScalarCache(body,g,scaling=GridScaling)

#=
## Basic regularization and interpolation
Let's perform an example in which we regularize the x coordinate of
the surface points onto the grid. We can get the surface coordinates
by using the [`points`](@ref) function. This can be applied either to
`body` directly or to the `cache`. As for any `VectorData` type, the components
of `pts` are `pts.u` and `pts.v`.
=#
pts = points(cache);

#=
Now, we set up some blank grid data on which to regularize onto
=#
gx = zeros_grid(cache);

#=
Now regularize. We will time it to show that it is fast and memory-efficient:
=#
regularize!(gx,pts.u,cache) #hide
@time regularize!(gx,pts.u,cache);

#=
Let's plot this to look at it. This also gives us a chance to highlight
the plot recipe for grid data associated with the cache, which is achieved by
simply supplying the cache to the `plot` function in `Plots.jl`. By default,
this plots the immersed points, as well, but this can be suppressed by
adding the keyword `layers=false`.
=#
plot(gx,cache)

#=
This shows how the regularization spreads the data over a couple of cells
around the surface. In the parlance of potential theory, this is a **single layer**.
=#
#=
If we wish to interpolate, then we do so with the [`interpolate!`]
function. For example, suppose we have a uniform field that we wish to interpolate:
=#
oc = 2.5*ones_grid(cache)

#=
Now set up some surface data to receive the interpolated data, and
interpolate:
=#
f = zeros_surface(cache);
interpolate!(f,oc,cache)

#=
It is clear that the interpolation preserves the value of the
field. This is also true for linearly varying fields, since the DDF
is built to ensure this. Let's try this, using the $x$ coordinate of the grid. Here, we use the `coordinates`
function of `CartesianGrids.jl`, which gets the coordinates of the grid,
and set the values of grid data to the $x$ coordinate. We interpolate and plot,
comparing to the actual $x$ coordinate of the points on the body:
=#
#x, y = coordinates(oc,g)
#xg = similar(oc)
#xg .= x
xg = x_grid(cache)
interpolate!(f,xg,cache)
plot(f,ylim=(-2,2),label="Interpolated from grid",ylabel="x",xlabel="Index")
plot!(pts.u,label="Actual body coordinate")

#=
## A double layer
Now we will generate a double layer. Mathematically, this takes the form

$$D_s f = \nabla\cdot \left( \delta(\chi) \mathbf{n} f \right)$$

for some scalar data $f$ on the surface. (See [Background](@ref) for an example.)
Notice that it maps scalar data on the surface ($f$) to scalar data in space. So to
calculate this using our discrete tools, we set up some grid data to receive the
result. Then, we use the function [`surface_divergence!`](@ref) to compute the
double layer. Here, we will demonstrate this on the $y$ coordinate of the
surface points:
=#
dl = zeros_grid(cache)
surface_divergence!(dl,pts.v,cache)
plot(dl,cache)

#=
If the surface data are vectors, $\mathbf{f}$, then this operation is a little
different:

$$D_s \mathbf{f} = \nabla\cdot \delta(\chi) \left(  \mathbf{n} \mathbf{f} + \mathbf{f} \mathbf{n} \right)$$

This maps $\mathbf{f}$ to a vector field. We use this in conjunction with a cache
generated with [`SurfaceVectorCache`](@ref).
=#

#=
The transpose of the double layer, $D_s$, is the operation

$$G_s u = \mathbf{n}\cdot \delta^{T}(\chi) \nabla u$$

for some scalar field data $u$. This operation computes the gradient of the field
data, interpolates this gradient onto the surface, and obtains the normal
component of that surface data. As such, it represents an important tool
for computing the normal derivative of field data on a surface. In the
package, we use [`surface_grad!`](@ref) for this operation.

The vector field version of this is

$$G_s \mathbf{u} = \mathbf{n}\cdot \delta^{T}(\chi) (\nabla \mathbf{u} + \nabla^{T} \mathbf{u})$$

which maps vector field data $\mathbf{u}$ to vector-valued surface data.
=#

#=
## A curl layer
We also sometimes need to take the curl of the regularized surface data,

$$C_s f = \nabla\times \left( \delta(\chi) \mathbf{n} f \right)$$

For this, we use the [`surface_curl!`](@ref) operator. Let's demonstrate this
on a uniform field on the surface.
=#
gc = zeros_gridcurl(cache)
f = ones_surface(cache)
surface_curl!(gc,f,cache)
plot(gc,cache)

#=
The continuous version of this operation is actually zero. It's not quite
zero in discrete form. However, its norm is much smaller than that of
the double layer.
=#
norm(gc,cache)/norm(dl,cache)

#=
Finally, a pair of operations that are used in support of the
previous ones, or occasionally on their own, are

$$R_n f = \delta(\chi)\mathbf{n}\circ f$$

for scalar surface data $f$, which maps to a vector field, effectively a
field of doublet strengths; and its transpose

$$R_n^T \mathbf{u} = \mathbf{n}\cdot \delta^{T}(\chi)\mathbf{u}$$

which maps vector field data $\mathbf{u}$ to a scalar surface field, the normal
component of the vector field on the surface. These are
provided by [`regularize_normal!`](@ref) and [`normal_interpolate!`](@ref),
respectively.
=#

#=
## Masks
Masks are grid data that take the value 1 in one region (e.g., the interior of a surface)
and 0 in the other (e.g., the exterior). The functions [`mask`](@ref)
and [`complementary_mask`](@ref) achieve this
=#
m = mask(cache)
cm = complementary_mask(cache)
plot(
    surface(m,cache,layers=false),
    surface(cm,cache,layers=false)
    )

#=
One can apply a mask to some grid data by multiplying it, using, e.g.,
the `product!` function in `CartesianGrids.jl`. Let's demonstrate that
with the grid data of $x$ coordinates:
=#
xmask = zeros_grid(cache)
xcmask = zeros_grid(cache)
product!(xmask,xg,m)
product!(xcmask,xg,cm)
plot(
  plot(xmask,cache),
  plot(xcmask,cache)
  )

#=
The mask and complementary mask effectively partition the field into two parts.
We can also apply masks in place, using [`mask!`](@ref) and [`complementary_mask!`](@ref):
=#
xmask .= xg
mask!(xmask,cache)
plot(xmask,cache)

#md # ## Surface-grid operator functions
#md # ```@docs
#md # regularize!
#md # interpolate!
#md # regularize_normal!
#md # normal_interpolate!
#md # regularize_normal_cross!
#md # normal_cross_interpolate!
#md # surface_divergence!
#md # surface_grad!
#md # surface_curl!
#md # surface_divergence_cross!
#md # surface_grad_cross!
#md # surface_curl_cross!
#md # mask!
#md # mask
#md # complementary_mask!
#md # complementary_mask
#md # ```


text:
module VanFoFy

include("all_modules.jl")

using .Input: CellData
using .Ellipticals: Weierstrass, Theta, theta
using .Symbolic: VarLinForm, variables, add!, add_conjugated!, mul!

function construct_problem(cell::CellData)
    wei = Weierstrass(cell)
    #Q = Q_special(wei)
    #lshear = construct_lshear(cell, wei)
    #plane  = construct_plane(cell, wei, Q)
    #lext   = construct_lext(cell, plane)
end

export Theta, theta
export VarLinForm, variables, add!, add_conjugated!, mul!

end # module

ext: jl


text:
# This file was generated, do not modify it. # hide
mao = dataset("gap", "mao")
describe(mao, :nmissing)


text:
# Master REPL functions
using Spark

# Load a list of workers from a file, contact them
function load(master::Master, configfile)
    f = open(configfile)
    config = JSON.parse(readall(f))
    for worker in config
        # Read the configuration file and try to connect to all the workers
        hostname = worker[1]
        port = worker[2]
        socket = None
        active = false
        try
            # Try to connect to all the clients
            socket = connect(hostname, port)
            active = true
        catch e
            println("Couldn't connect to $hostname:$port...")
            showerror(STDERR, e)
        end
        wf = WorkerRef(hostname, port, active)
        master.workers = cat(1, master.workers, [wf])
        master.sockets[wf] = socket
    end
    identify(master)
end

# Set up the local RPC server for worker->master RDD requests
function initserver(master::Master)
    server = listen(IPv4(0), master.port)
    println("Starting server")
    @async while true
        sock = accept(server)
        @async while true
            try
                line = readline(sock)
                response = handle(master, line)
                println(sock, response)
            catch e
                showerror(STDERR, e)
                break
            end
        end
    end
end

# General handle, in case we want more Worker -> Master RPC like ping
function handle(master::Master, line::ASCIIString)
    # General format of a message:
    # {:call => "funcname", :args => {anything}}
    # this is dispatched to any function call - fine since we're assuming
    # a private non-adversarial network.
    msg = JSON.parse(strip(line))
    r = false
    if "call" in keys(msg)
        r = eval(Expr(:call, symbol(msg["call"]), master, msg["args"]))
    end
    return r
end

#recovers given partition
function recover_partition(master::Master, rdd_id::Int64, partition_id::Int64)

    #Check if rdd comes directly from disk. If so, recover partition and return
    if length(master.rdds[rdd_id].dependencies) == 0
        recover_part(master, master.rdds[rdd_id], partition_id)
        return
    end

    lost_partitions::Array{(Int64, Int64)} = Array((Int64, Int64), 0)
    #finds lost predecessors
    for rdd in master.rdds[rdd_id].dependencies
        for partition in rdd[2]
            if !partition[2].active
                push!(lost_partitions, (rdd[1], partition[1]))
            end
        end
    end

    if length(lost_partitions) == 0
        return
    else
        for partition in lost_partitions
            recover_partition(master, partition[1], partition[2])
        end
        recover_part(master, master.rdds[rdd_id], partition_id)
    end
end

#recovers all lost partitions of the given rdd
function recover(master::Master, rdd::RDD)
    rdd_id = rdd.ID
    lost_partitions::Array{Int64} = Array(Int64, 0)
    for worker in master.rdds[rdd_id].partitions
        if !worker[2].active
            push!(lost_partitions, worker[1])
        end
    end

    for partition_id in lost_partitions
        recover_partition(master, rdd_id, partition_id)
    end
    return partition_by(master, master.rdds[rdd_id], HashPartitioner())
end


text:
using FastTransforms
import Base: *

#=
N = 2.^(7:10)

for n in N
    A = sphones(Float64, n, n)
    B = deepcopy(A)

    @time PC = c_plan_sph2fourier(n)
    @time c_sph2fourier(PC, B)

    @time PJ = SlowSphericalHarmonicPlan(A)
    @time A = PJ*A;

    @show vecnorm(A-B)/vecnorm(A)

    println()
    println()
end
=#

#=
α, β, γ = 0.0, -0.5, -0.5
#α, β, γ = -0.5, -0.5, -0.5
#α, β, γ = 0.0, 0.0, 0.0

for n in N
    A = trirand(Float64, n, n)
    B = deepcopy(A)

    @time PC = c_plan_tri2cheb(n, α, β, γ)
    @time c_tri2cheb(PC, B)

    @time PJ = SlowTriangularHarmonicPlan(A, α, β, γ);
    @time A = PJ*A;

    @show vecnorm(A-B)/vecnorm(B)

    println()
    println()
end

CF = zero(F); CF[1,2] = 1.0; CF[2,2] = 1.0; CF[3,2] = 1.0;
c_tri2cheb(cp, CF)
CF
c_cheb2tri(cp, CF)
CF

cprot = c_plan_rottriangle(size(F, 1), α, β, γ)

CF = zero(F); CF[1,2] = 1.0; CF[2,2] = 1.0; CF[3,2] = 1.0;
c_execute_tri_hi2lo(cprot, CF)
CF
c_execute_tri_lo2hi(cprot, CF)
CF
=#

using ApproxFun, FastTransforms

import ApproxFun: jacobip

jacobinorm(n,a,b) = if n ≠ 0
    sqrt((2n+a+b+1))*exp((lgamma(n+a+b+1)+lgamma(n+1)-log(2)*(a+b+1)-lgamma(n+a+1)-lgamma(n+b+1))/2)
else
    sqrt(exp(lgamma(a+b+2)-log(2)*(a+b+1)-lgamma(a+1)-lgamma(b+1)))
end

njacobip(n,a,b,x) = jacobinorm(n,a,b) * jacobip(n,a,b,x)

α = 0.0#0.25
β = 0.0#-0.25
γ = 0.0#-0.125

P = (ℓ,m,x,y) -> (2*(1-x))^m*njacobip(ℓ-m,2m+β+γ+1,α,2x-1)*njacobip(m,γ,β,2y/(1-x)-1)

p_T = chebyshevpoints(10)

f = (x,y) -> P(0,0,x,y) + P(1,1,x,y) + P(2,2,x,y) + P(3,3,x,y) + P(4,4,x,y) +
             P(1,0,x,y) + P(2,1,x,y) + P(3,2,x,y) +
             P(2,0,x,y) + P(3,1,x,y) +
             P(3,0,x,y)

f̃ = (s,t) -> f((s+1)/2, (t+1)/2*(1-(s+1)/2))

F = f̃.(p_T, p_T')

for j = 1:size(F,2)
    F[:,j] = chebyshevtransform(F[:,j])
end

for k = 1:size(F,1)
    F[k,:] = chebyshevtransform(F[k,:])
end

P = Tri2Cheb(

cp = c_plan_tri2cheb(size(F, 1), α, β, γ)

CF̌ = copy(F)

c_cheb2tri(cp, CF̌)

CF̌


text:
using Test
using StochasticDiffEq
using SparseArrays
using LinearAlgebra

#https://github.com/JuliaDiffEq/SparseDiffTools.jl/blob/master/test/test_integration.jl
function f(dx,x,p,t)
  for i in 2:length(x)-1
    dx[i] = x[i-1] - 2x[i] + x[i+1]
  end
  dx[1] = -2x[1] + x[2]
  dx[end] = x[end-1] - 2x[end]
  nothing
end

function g(dx,x,p,t)
  dx .= 0
  nothing
end

function second_derivative_stencil(N)
  A = zeros(N,N)
  for i in 1:N, j in 1:N
    (j-i==-1 || j-i==1) && (A[i,j]=1)
    j-i==0 && (A[i,j]=-2)
  end
  A
end

function generate_sparsity_pattern(N::Integer)
  dl = repeat([1.0],N-1)
  du = repeat([1.0],N-1)
  d = repeat([-2.0],N)
  return Tridiagonal(dl,d,du)
end

jac_sp = sparse(generate_sparsity_pattern(10))
#jac = second_derivative_stencil(10)
colorvec = repeat(1:3,10)[1:10]
u0=[1.,2.,3,4,5,5,4,3,2,1]
tspan=(0.,10.)
sdefun_sp= SDEFunction(f,g,colorvec=colorvec,jac_prototype=jac_sp)
prob_sp = SDEProblem(sdefun_sp,g,u0,tspan)
prob_std = SDEProblem(f,g,u0,tspan)

sol_sp=solve(prob_sp,SKenCarp(autodiff=false))
@test sol_sp.retcode==:Success#test sparse finitediff
sol=solve(prob_std,SKenCarp(autodiff=false))
@test sol_sp.u[end]≈sol.u[end] atol=1e-4
@test length(sol_sp.t)==length(sol.t)

sol_sp=solve(prob_sp,SKenCarp())
sol=solve(prob_std,SKenCarp())
@test sol_sp.u[end]≈sol.u[end] atol=1e-4
@test length(sol_sp.t)==length(sol.t)

#sol_sp=solve(prob_sp,SKenCarp(autodiff=false),abstol=1e-10,reltol=1e-10)
#sol=solve(prob_std,SKenCarp(autodiff=false),abstol=1e-10,reltol=1e-10)
#@test sol_sp.u[end]≈sol.u[end]
#@test length(sol_sp.t)==length(sol.t)


text:
using Random

"""
    function generate_players(;
        dataset_size::Int64,
        ranking_schemes::Vector{String},
        verbose::Bool
    )::Vector{Player}

Generates an array of players.

Arguments
=========
- `dataset_size`:
    the number of players to create.
- `ranking_schemes`:
    the ranking algorithms that will be applied.
- `verbose`:
    conditional on whether to print results.
"""

function generate_players(;
    dataset_size::Int64,
    ranking_schemes::Vector{String},
    verbose::Bool
)::Vector{Player}
    players = Vector{Player}([])
    burn_in_rankings = shuffle!(Vector{Int64}(collect(1:dataset_size)))
    for i in 1:dataset_size
        ground_truth_ranking = i
        burn_in_ranking = burn_in_rankings[i]
        p = Player(
            id = i,
            ground_truth_ranking = ground_truth_ranking,
            burn_in_ranking = burn_in_ranking,
            ranking_schemes = ranking_schemes,
        )
        push!(players, p)
    end
    if verbose
        players_count = length(players)
        println("########")
        println("Vector of $players_count created.")
        println("########")
        println()
    end
    return players
end

max_stars_repo_name: Danial-Hussain/Tennis-Ranking

max_issues_repo_name: Danial-Hussain/Tennis-Ranking

max_forks_repo_name: Danial-Hussain/Tennis-Ranking


text:
# In this tutorial, we will learn
#  - How to solve a simple PDE with a DG method
#  - How to compute jumps and averages of quantities on the mesh skeleton
#  - How to implement the method of manufactured solutions
#  - How to integrate error norms
#  - How to generate Cartesian meshes in arbitrary dimensions
#
#
# ## Problem statement
#
# The goal of this tutorial is to solve a PDE using a DG formulation. For simplicity, we take the Poisson equation on the unit cube $\Omega \doteq (0,1)^3$ as the model problem, namely

# ```math
# \left\lbrace
# \begin{aligned}
# -\Delta u = f  \ &\text{in} \ \Omega,\\
# u = g \ &\text{on}\ \partial\Omega,\\
# \end{aligned}
# \right.
# ```
# where $f$ is the source term and $g$ is the prescribed Dirichlet boundary function. In this tutorial, we follow the method of manufactured solutions since we want to illustrate how to compute discretization errors. We take $u(x) = 3 x_1 + x_2 + 2 x_3$ as the exact solution of the problem, for which $f=0$ and $g(x) = u(x)$. The selected manufactured solution $u$ is a first order multi-variate polynomial, which can be represented exactly by the FE interpolation that we are going to define below. In this scenario, the discretization error has to be close to the machine precision. We will use this result to validate the proposed implementation.
#
# ## Numerical Scheme
#
# We consider a DG formulation to approximate
# the problem. In particular, we consider  the symmetric
# interior penalty method (see, e.g. [1], for specific details). For this formulation, the approximation space is made of discontinuous piece-wise polynomials, namely
#
# ```math
# V \doteq \{ v\in L^2(\Omega):\ v|_{T}\in Q_p(T) \text{ for all } T\in\mathcal{T}  \},
# ```
# where $\mathcal{T}$ is the set of all cells $T$ of the FE mesh, and $Q_p(T)$ is a polynomial space of degree $p$ defined on a generic cell $T$. For simplicity, we consider Cartesian meshes in this tutorial. In this case, the space $Q_p(T)$ is made of multi-variate polynomials up to degree $p$ in each spatial coordinate.
#
# In order to write the weak form of the problem, we need to introduce some notation. The sets of interior and boundary facets associated with the FE mesh $\mathcal{T}$ are denoted here as $\mathcal{F}_\Gamma$ and $\mathcal{F}_{\partial\Omega}$ respectively. In addition, for a given function $v\in V$ restricted to the interior facets $\mathcal{F}_\Gamma$, we introduce the well known jump and mean value operators,
# ```math
# \begin{aligned}
# \lbrack\!\lbrack v\ n \rbrack\!\rbrack &\doteq v^+\ n^+ + v^- n^-,\\
# \{\! \!\{\nabla v\}\! \!\} &\doteq \dfrac{ \nabla v^+ + \nabla v^-}{2},
# \end{aligned}
# ```
# with $v^+$, and $v^-$ being the restrictions of $v\in V$ to the cells $T^+$, $T^-$ that share a generic interior facet in $\mathcal{F}_\Gamma$, and $n^+$, and $n^-$ are the facet outward unit normals from either the perspective of $T^+$ and $T^-$ respectively.
#
# With this notation, the weak form associated with the interior penalty formulation of our problem reads: find $u\in V$ such that $a(u,v) = b(v)$ for all $v\in V$.  The bilinear and linear forms  $a(\cdot,\cdot)$ and $b(\cdot)$ have contributions associated with the bulk of $\Omega$,  the boundary facets $\mathcal{F}_{\partial\Omega}$, and the interior facets  $\mathcal{F}_\Gamma$, namely
#  ```math
# \begin{aligned}
# a(u,v) &= a_{\Omega}(u,v) + a_{\partial\Omega}(u,v) + a_{\Gamma}(u,v),\\
# b(v) &= b_{\Omega}(v) + b_{\partial\Omega}(v).
# \end{aligned}
# ```
# These contributions are defined as
# ```math
# \begin{aligned}
# a_{\Omega}(u,v) &\doteq \sum_{T\in\mathcal{T}} \int_{T} \nabla v \cdot \nabla u \ {\rm d}T,\\
# b_{\Omega}(v) &\doteq \int_{\Omega} v\ f \ {\rm d}\Omega,
# \end{aligned}
# ```
# for the volume,
# ```math
# \begin{aligned}
# a_{\partial\Omega}(u,v) &\doteq \sum_{F\in\mathcal{F}_{\partial\Omega}} \dfrac{\gamma}{|F|} \int_{F} v\ u \ {\rm d}F \\ & -  \sum_{F\in\mathcal{F}_{\partial\Omega}} \int_{F} v\ (\nabla u \cdot n)  \ {\rm d}F \\ & -  \sum_{F\in\mathcal{F}_{\partial\Omega}} \int_{F} (\nabla v \cdot n)\ u  \ {\rm d}F, \\
# b_{\partial\Omega} &\doteq \sum_{F\in\mathcal{F}_{\partial\Omega}} \dfrac{\gamma}{|F|} \int_{F} v\ g \ {\rm d}F \\ & -  \sum_{F\in\mathcal{F}_{\partial\Omega}} \int_{F} (\nabla v \cdot n)\ g  \ {\rm d}F,
# \end{aligned}
# ```
# for the boundary facets and,
# ```math
# \begin{aligned}
# a_{\Gamma}(u,v) &\doteq \sum_{F\in\mathcal{F}_{\Gamma}} \dfrac{\gamma}{|F|} \int_{F} \lbrack\!\lbrack v\ n \rbrack\!\rbrack\cdot \lbrack\!\lbrack u\ n \rbrack\!\rbrack \ {\rm d}F \\ &-  \sum_{F\in\mathcal{F}_{\Gamma}} \int_{F} \lbrack\!\lbrack v\ n \rbrack\!\rbrack\cdot \{\! \!\{\nabla u\}\! \!\} \ {\rm d}F \\ & -  \sum_{F\in\mathcal{F}_{\Gamma}} \int_{F} \{\! \!\{\nabla v\}\! \!\}\cdot \lbrack\!\lbrack u\ n \rbrack\!\rbrack \ {\rm d}F,
# \end{aligned}
# ```
#  for the interior facets. In previous expressions, $|F|$ denotes the diameter of the face $F$ (in our Cartesian grid, this is equivalent to the characteristic mesh size $h$), and $\gamma$ is a stabilization parameter that should be chosen large enough such that the bilinear form $a(\cdot,\cdot)$ is stable and continuous. Here, we take $\gamma = p\ (p+1)$ as done in the numerical experiments in reference [2].
#
# ## Manufactured solution
#
# We start by loading the Gridap library and defining the manufactured solution $u$ and the associated source term $f$ and Dirichlet function $g$.

using Gridap
u(x) = 3*x[1] + x[2] + 2*x[3]
f(x) = 0
g(x) = u(x)

# We also need to define the gradient of $u$ since we will compute the $H^1$ error norm later. In that case, the gradient is simply defined as

∇u(x) = VectorValue(3,1,2)

# In addition, we need to tell the Gridap library that the gradient of the function `u` is available in the function `∇u` (at this moment `u` and `∇u` are two standard Julia functions without any connection between them). This is done by adding an extra method to the function `gradient` (aka `∇`) defined in Gridap:

import Gridap: ∇
∇(::typeof(u)) = ∇u

#  Now, it is possible to recover function `∇u` from function `u` as `∇(u)`. You can check that the following expression evaluates to `true`.

∇(u) === ∇u

# ## Cartesian mesh generation
#  In order to discretize the geometry of the unit cube, we use the Cartesian mesh generator available in Gridap.

L = 1.0
domain = (0.0, L, 0.0, L, 0.0, L)
n = 4
partition = (n,n,n)
model = CartesianDiscreteModel(domain,partition)

# The type `CartesianDiscreteModel` is a concrete type that inherits from `DiscreteModel`, which is specifically designed for building Cartesian meshes. The `CartesianDiscreteModel` constructor takes a tuple containing limits of the box we want to discretize  plus a tuple with the number of cells to be generated in each direction (here $4\times4\times4$ cells). You can write the model in vtk format to visualize it (see next figure).

 writevtk(model,"model")

# ![](../assets/dg_discretization/model.png)
#
#  Note that the `CaresianDiscreteModel` is implemented for arbitrary dimensions. For instance, the following lines build a `CartesianDiscreteModel`  for the unit square $(0,1)^2$ with 4 cells per direction

domain2D = (0.0, L, 0.0, L)
partition2D = (n,n)
model2D = CartesianDiscreteModel(domain2D,partition2D)

# You could also generate a mesh for the unit tesseract $(0,1)^4$ (i.e., the unit cube in 4D). Look how the 2D and 3D models are built and just follow the sequence.
#
# ## FE spaces
#
# On top of the discrete model, we create the discontinuous space $V$ as follows

order = 3
V = TestFESpace(
  reffe=:Lagrangian, valuetype=Float64, order=order,
  conformity=:L2, model=model)

# We have select a Lagrangian, scalar-valued interpolation of order $3$ within the cells of the discrete model. Since the cells are hexahedra, the resulting Lagrangian shape functions are tri-cubic polynomials. In contrast to previous tutorials, where we have constructed $H^1$-conforming (i.e., continuous) FE spaces, here we construct a $L^2$-conforming (i.e., discontinuous) FE space. That is, we do not impose any type of continuity of the shape function on the cell boundaries, which leads to the discontinuous FE space $V$ of the DG formulation. Note also that we do not pass any information about the Dirichlet boundary to the `TestFESpace` constructor since the Dirichlet boundary conditions are not imposed strongly in this example.
#
# From the `V` object we have constructed in previous code snippet, we build the trial FE space as usual.

U = TrialFESpace(V)

# Note that we do not pass any Dirichlet function to the `TrialFESpace` constructor since we do not impose Dirichlet boundary conditions strongly here.
#
# ## Numerical integration
#
# Once the FE spaces are ready, the next step is to set up  the numerical integration. In this example, we need to integrate in three different domains: the volume covered by the cells $\mathcal{T}$  (i.e., the computational domain $\Omega$), the surface covered by the boundary facets $\mathcal{F}_{\partial\Omega}$ (i.e., the boundary $\partial\Omega$), and the surface covered by the interior facets $\mathcal{F}_{\Gamma}$ (i.e. the so-called mesh skeleton). In order to integrate in $\Omega$ and on its boundary $\partial\Omega$, we use `Triangulation` and `BoundaryTriangulation` objects as already discussed in previous tutorials.

trian = Triangulation(model)
btrian = BoundaryTriangulation(model)

# Here, we do not pass any boundary identifier to the `BoundaryTriangulation` constructor. In this case, an integration mesh for the entire boundary $\partial\Omega$ is constructed by default (which is just what we need in this example).
#
# In order to generate an integration mesh for the interior facets $\mathcal{F}_{\Gamma}$, we use a new type of `Triangulation` referred to as `SkeletonTriangulation`. It can be constructed from a `DiscreteModel` object as follows:

strian = SkeletonTriangulation(model)

# As any other type of `Triangulation`, an `SkeletonTriangulation` can be written into a vtk file for its visualization (see next figure, where the interior facets $\mathcal{F}_\Gamma$ are clearly observed).

writevtk(strian,"strian")

# ![](../assets/dg_discretization/skeleton_trian.png)
#
# Once we have constructed the triangulations needed in this example, we define the corresponding quadrature rules.

degree = 2*order
quad = CellQuadrature(trian,degree)
bquad = CellQuadrature(btrian,degree)
squad = CellQuadrature(strian,degree)

# We still need a way to represent the unit outward normal vector to the boundary $\partial\Omega$, and the unit normal vector on the interior faces $\mathcal{F}_\Gamma$. This is done with the `get_normal_vector` getter.

nb = get_normal_vector(btrian)
ns = get_normal_vector(strian)

# The `get_normal_vector` getter takes either a boundary or a skeleton triangulation and returns an object representing the normal vector to the corresponding surface. For boundary triangulations, the returned normal vector is the unit outwards one, whereas for skeleton triangulations the orientation of the returned normal is arbitrary. In the current implementation (Gridap v0.5.0), the unit normal is outwards to the cell with smaller id among the two cells that share an interior facet in $\mathcal{F}_\Gamma$.
#
# ## Weak form
#
# With these ingredients we can define the different terms in the weak form. First, we start with the terms $a_\Omega(\cdot,\cdot)$ , and $b_\Omega(\cdot)$ associated with integrals in the volume $\Omega$. This is done as in the tutorial for the Poisson equation.

a_Ω(u,v) = ∇(v)⊙∇(u)
b_Ω(v) = v*f
t_Ω = AffineFETerm(a_Ω,b_Ω,trian,quad)

# The terms $a_{\partial\Omega}(\cdot,\cdot)$ and $b_{\partial\Omega}(\cdot)$ associated with integrals on the boundary $\partial\Omega$ are defined using an analogous approach. First, we define two functions representing the integrands of the forms $a_{\partial\Omega}(\cdot,\cdot)$ and $b_{\partial\Omega}(\cdot)$. Then, we build an `AffineFETerm` from these functions and the boundary triangulation and its corresponding quadrature rule:

h = L / n
γ = order*(order+1)
a_∂Ω(u,v) = (γ/h)*v*u - v*(∇(u)⋅nb) - (∇(v)⋅nb)*u
b_∂Ω(v) = (γ/h)*v*g - (∇(v)⋅nb)*g
t_∂Ω = AffineFETerm(a_∂Ω,b_∂Ω,btrian,bquad)

# Note that in the definition of the functions `a_∂Ω` and `b_∂Ω`, we have used the object `nb` representing the outward unit normal to the boundary $\partial\Omega$. The code definition of  `a_∂Ω` and `b_∂Ω` is indeed very close to the mathematical definition of the forms  $a_{\partial\Omega}(\cdot,\cdot)$ and $b_{\partial\Omega}(\cdot)$.
#
# Finally, we need to define the term $a_\Gamma(\cdot,\cdot)$ integrated on the interior facets $\mathcal{F}_\Gamma$. In this case, we use a `LinearFETerm` since the terms integrated on the interior facets only contribute to the system matrix and not to the right-hand-side vector.

a_Γ(u,v) = (γ/h)*jump(v*ns)⊙jump(u*ns) - jump(v*ns)⊙mean(∇(u)) - mean(∇(v))⊙jump(u*ns)
t_Γ = LinearFETerm(a_Γ,strian,squad)

# Note that the arguments `v`, `u` of function  `a_Γ` represent a test and trial function *restricted* to the interior facets $\mathcal{F}_\Gamma$. As mentioned before in the presentation of the DG formulation, the restriction of a function $v\in V$ to the interior faces leads to two different values $v^+$ and $v^-$ . In order to compute jumps and averages of the quantities $v^+$ and $v^-$, we use the functions `jump` and `mean`, which represent the jump and mean value operators $\lbrack\!\lbrack \cdot \rbrack\!\rbrack$ and $\{\! \!\{\cdot\}\! \!\}$ respectively. Note also that we have used the object `ns` representing the unit normal vector on the interior facets. As a result, the notation used to define function `a_Γ` is very close to the mathematical definition of the terms in the bilinear form $a_\Gamma(\cdot,\cdot)$.
#
# Once the different terms of the weak form have been defined, we build and solve the FE problem.

op = AffineFEOperator(U,V,t_Ω,t_∂Ω,t_Γ)
uh = solve(op)

# ## Discretization error
#
# We end this tutorial by quantifying the discretization error associated with  the computed numerical solution `uh`. In DG methods a simple error indicator is the jump of the computed (discontinuous) approximation on the interior faces. This quantity can be easily computed in Gridap as follows. First, we need to restrict the computed solution `uh` to the skeleton triangulation.

uh_Γ = restrict(uh,strian)

# The resulting object `uh_Γ` is an object which represents the two values $u^+_h$, $u^-_h$ of the solution $u_h$ restricted to the interior facets $\mathcal{F}_\Gamma$. We compute and visualize the jump of these values as follows (see next figure):

writevtk(strian,"jumps",cellfields=["jump_u"=>jump(uh_Γ)])

# Note that the jump of the numerical solution is very small, close to the machine precision (as expected in this example with manufactured solution).
# ![](../assets/dg_discretization/jump_u.png)
#
#  A more rigorous way of quantifying the error is to measure it with a norm. Here, we use the $L^2$ and $H^1$ norms, namely
#  ```math
# \begin{aligned}
#  \| w \|_{L^2}^2 & \doteq \int_{\Omega} w^2 \ \text{d}\Omega, \\
#  \| w \|_{H^1}^2 & \doteq \int_{\Omega} w^2 + \nabla w \cdot \nabla w \ \text{d}\Omega.
# \end{aligned}
# ```
#
# The discretization error can be computed in this example as the difference of the manufactured and numerical solutions.

e = u - uh

# We compute the error norms as follows. First, we implement the integrands of the norms we want to compute.

l2(u) = u*u
h1(u) = a_Ω(u,u) + l2(u)

# Then, we compute the corresponding integrals with the `integrate` function.

el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))

# The `integrate` function returns a lazy object representing the contribution to the integral of each cell in the underlying triangulation. To end up with the desired error norms, one has to sum these contributions and take the square root. You can check that the computed error norms are close to machine precision (as one would expect).

tol = 1.e-10
@assert el2 < tol
@assert eh1 < tol

# ## References
#
# [1] D. N. Arnold, F. Brezzi, B. Cockburn, and L. Donatella Marini. Unified analysis of discontinuous Galerkin methods for elliptic problems. *SIAM Journal on Numerical Analysis*, 39 (5):1749–1779, 2001. doi:[10.1137/S0036142901384162](http://dx.doi.org/10.1137/S0036142901384162).
#
# [2] B. Cockburn, G. Kanschat, and D. Schötzau. An equal-order DG method for the incompressible Navier-Stokes equations. *Journal of Scientific Computing*, 40(1-3):188–210, 2009. doi:[10.1007/s10915-008-9261-1](http://dx.doi.org/10.1007/s10915-008-9261-1).
#

max_stars_repo_name: andre-massing/Tutorials

max_issues_repo_name: andre-massing/Tutorials

max_forks_repo_name: andre-massing/Tutorials


text:
module SurvivorSelection
    using ..Representation
    # Rank selection for (μ + λ)
    function rank_selection(offspring::Representation.Population,
                                parents::Representation.Population, μ_selected::Int)
        sorted_population = sort([offspring.members; parents.members], by=i->i.fitness)
        return Representation.Population(sorted_population[1:μ_selected])
    end
    # Rank selection for (μ, λ)
    function rank_selection(offspring::Representation.Population, μ_selected::Int)
        sorted_population = sort(offspring.members, by=i->i.fitness)
        return Representation.Population(sorted_population[1:μ_selected])
    end

end # module SurvivorSelection


text:
mutable struct GeoArray{T<:Union{Real, Union{Missing, Real}}} <: AbstractArray{T, 3}
    A::AbstractArray{T, 3}
    f::AffineMap
    crs::WellKnownText{GeoFormatTypes.CRS, <:String}
end
GeoArray(A::AbstractArray{T, 3} where T<:Union{Real, Union{Missing, Real}}) = GeoArray(A, geotransform_to_affine(SVector(0.,1.,0.,0.,0.,1.)), "")
GeoArray(A::AbstractArray{T, 3} where T<:Union{Real, Union{Missing, Real}}, f::AffineMap) = GeoArray(A, f, GFT.WellKnownText(GFT.CRS(), ""))
GeoArray(A::AbstractArray{T, 3} where T<:Union{Real, Union{Missing, Real}}, f::AffineMap, crs::String) = GeoArray(A, f, GFT.WellKnownText(GFT.CRS(), crs))
GeoArray(A::AbstractArray{T, 2} where T<:Union{Real, Union{Missing, Real}}, args...) = GeoArray(reshape(A, size(A)..., 1), args...)
function GeoArray(A::AbstractArray{T, 3} where T<:Union{Real, Union{Missing, Real}}, x::AbstractRange, y::AbstractRange, args...)
    size(A)[1:2] != (length(x), length(y)) && error("Size of `GeoArray` $(size(A)) does not match size of (x,y): $((length(x),length(y))). Note that this function takes *center coordinates*.")
    f = unitrange_to_affine(x, y)
    GeoArray(A, f, args...)
end

Base.size(ga::GeoArray) = size(ga.A)
Base.IndexStyle(::Type{T}) where {T<:GeoArray} = IndexLinear()
Base.getindex(ga::GeoArray, i::Int) = getindex(ga.A, i)
Base.getindex(ga::GeoArray, I::Vararg{Int, 2}) = getindex(ga.A, I..., :)
Base.getindex(ga::GeoArray, I::Vararg{Int, 3}) = getindex(ga.A, I...)
Base.iterate(ga::GeoArray) = iterate(ga.A)
Base.length(ga::GeoArray) = length(ga.A)
Base.parent(ga::GeoArray) = ga.A
# Base.map(f, ga::GeoArray) = GeoArray(map(f, ga.A), ga.f, ga.crs)
# Base.convert(::Type{Array{T, 3}}, A::GeoArray{T}) where {T} = convert(Array{T,3}, ga.A)
Base.eltype(::Type{GeoArray{T}}) where {T} = T

Base.show(io::IO, ::MIME"text/plain", ga::GeoArray) = show(io, ga)
function Base.show(io::IO, ga::GeoArray)
    crs = GeoFormatTypes.val(ga.crs)
    wkt = length(crs) == 0 ? "undefined CRS" : "CRS $crs"
    print(io, "$(join(size(ga), "x")) $(typeof(ga.A)) with $(ga.f) and $(wkt)")
end

# Generate upper left coordinates for specic index
function coords(ga::GeoArray, p::SVector{2, Int})
    ga.f(p.-1)
end
coords(ga::GeoArray, p::Vector{Int}) = coords(ga, SVector{2}(p))

# Generate center coordinates for specific index
function centercoords(ga::GeoArray, p::SVector{2, Int})
    ga.f(p.-0.5)
end
centercoords(ga::GeoArray, p::Vector{Int}) = centercoords(ga, SVector{2}(p))

# Convert coordinates back to indices
function indices(ga::GeoArray, p::SVector{2, Float64})
    map(x->round(Int, x), inv(ga.f)(p)::SVector{2, Float64}).+1
end
indices(ga::GeoArray, p::Vector{Float64}) = indices(ga, SVector{2}(p))
indices(ga::GeoArray, p::Tuple{Float64, Float64}) = indices(ga, SVector{2}(p))

# Overload indexing directly into GeoArray
# TODO This could be used for interpolation instead
function Base.getindex(ga::GeoArray, I::SVector{2, Float64})
    (i, j) = indices(ga, I)
    return ga[i, j, :]
end
Base.getindex(ga::GeoArray, I::Vararg{Float64, 2}) = Base.getindex(ga, SVector{2}(I))

# Generate coordinates for complete GeoArray
function coords(ga::GeoArray)
    (ui, uj) = size(ga)[1:2]
    ci = [coords(ga, SVector{2}(i,j)) for i in 1:ui+1, j in 1:uj+1]
end
function centercoords(ga::GeoArray)
    (ui, uj) = size(ga)[1:2]
    ci = [centercoords(ga, SVector{2}(i,j)) for i in 1:ui, j in 1:uj]
end
function centercoordsnotmissing(ga::GeoArray)
    (ui, uj) = size(ga)[1:2]
    ci = [centercoords(ga, SVector{2}(i,j)) for i in 1:ui, j in 1:uj if ~ismissing(ga.A[i, j])]
end
function centercoordsmissing(ga::GeoArray)
    (ui, uj) = size(ga)[1:2]
    ci = [centercoords(ga, SVector{2}(i,j)) for i in 1:ui, j in 1:uj if ismissing(ga.A[i, j])]
end
function indexmissing(ga::GeoArray)
    (ui, uj) = size(ga)[1:2]
    ci = [[i,j] for i in 1:ui, j in 1:uj if ismissing(ga.A[i, j])]
end

# Generate coordinates for one dimension of a GeoArray
function coords(ga::GeoArray, dim::Symbol)
    if is_rotated(ga)
        error("This method cannot be used for a rotated GeoArray")
    end
    if dim==:x
        ui = size(ga,1)
        ci = [coords(ga, SVector{2}(i,1))[1] for i in 1:ui+1]
    elseif dim==:y
        uj = size(ga,2)
        ci = [coords(ga, SVector{2}(1,j))[2] for j in 1:uj+1]
    else
        error("Use :x or :y as second argument")
    end
    return ci
end
function centercoords(ga::GeoArray, dim::Symbol)
    if is_rotated(ga)
        error("This method cannot be used for a rotated GeoArray")
    end
    if dim==:x
        ui = size(ga,1)
        ci = [centercoords(ga, SVector{2}(i,1))[1] for i in 1:ui]
    elseif dim==:y
        uj = size(ga,2)
        ci = [centercoords(ga, SVector{2}(1,j))[2] for j in 1:uj]
    else
        error("Use :x or :y as second argument")
    end
    return ci
end

"""Set AffineMap of `GeoArray` by specifying the *center coordinates* for each x, y dimension by a `UnitRange`."""
function centercoords!(ga, x::AbstractUnitRange, y::AbstractUnitRange)
    size(ga)[1:2] != (length(x), length(y)) && error("Size of `GeoArray` $(size(ga)) does not match size of (x,y): $((length(x),length(y))). Note that this function takes *center coordinates*.")
    ga.f = unitrange_to_affine(x, y)
    ga
end


text:
module LinearSVM

include("DualCD.jl")
include("Pegasos.jl")
include("ADMMDriver.jl")
#include("ADMMWorker.jl")
include("Utilities.jl")

# DualCD.jl:
export dualcd
export pegasos

# ADMMDriver.jl:
export distributeData
export admm

# Utilities.jl:
export loadDataset
export testResult

end # module


text:
# PROJECT_AFFINE    Project a point into an affine set.
#
#   project_affine(v,A,b) is the projection of v onto
#   the affine set { x | Ax = b }.
#
#   You can also call the function as
#
#     [x C d] = project_affine(v,A,b);
#
#   and then call it again with different argument v2 via
#
#     x2 = project_affine(v2,A,b,C,d);
#
#   If calling the function repeatedly with the same A and b,
#   all evaluations after the initial one will be much faster
#   when passing in the cached values C and d.

# TODO: Determine dimensions of all objects
function project_affine(
	v::Array,
	A::Array,
	b::Array,
	C::Array = Array(Float64, 0),
	d::Array = Array(Float64, 0)
)
	# x = v - pinv(A) * (A * v - b)
    if isempty(C) || isempty(d)
        pA = pinv(A)
        C = eye(length(v)) - pA * A
        d = pA * b
    end

	x = C * v + d

	return x, C, d
end


text:
# This file calculates scattering from the gravity Helmholtz equation via the adaptive spectral method.
# Δu + (E+x_2)u = 0,
# u|Γ = 0,
# u^i = Φ(⋅,⋅,E)|Γ
# u = u^i + u^s.
# The normal derivative ∂u/∂n of the entire wave is calculated on the sound-soft boundaries.
# The scattered wave is calculated by convolving ∂u/∂n with the fundamental solution.
# Then, the total wave is obtained by summing the incident and the scattered waves.

using ApproxFun, SingularIntegralEquations
include("Scatteraux.jl")

E = 20.
ω = 2π
ui = (x,y) ->  lhelmfs(complex(x,y),-5.0im,E)

# The gravity Helmholtz Green's function.
g3 = (x,y) ->  lhelmfs(x,y,E)
r = (x,y) ->  lhelm_riemann(x,y,E)

dom = ∪(Segment.([-10.0-3.0im,5.0,-2+5im],[-5.0+0.0im,10.0-3im,2+5im]))
sp = Space(dom)
cwsp = CauchyWeight(sp⊗sp,0)
uiΓ,⨍ = Fun(t->ui(real(t),imag(t)),sp),DefiniteLineIntegral(dom)

@time G = GreensFun(g3,cwsp;method=:unsplit)

@time ∂u∂n = ⨍[G]\uiΓ
println("The length of ∂u∂n is: ",ncoefficients(∂u∂n))
us = (x,y) ->  -linesum(g3,∂u∂n,complex(x,y))


text: multsum(n, m, lim) = sum(0:n:lim-1) + sum(0:m:lim-1) - sum(0:lcm(n,m):lim-1)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
Read/write a region of image data from multiple bands.

This method allows reading a region of one or more `GDALRasterBands` from this
dataset into a buffer, or writing data from a buffer into a region of the
`GDALRasterBands`. It automatically takes care of data type translation if the
data type (`eBufType`) of the buffer is different than that of the
`GDALRasterBand`. The method also takes care of image decimation / replication
if the buffer size (`nBufXSize x nBufYSize`) is different than the size of the
region being accessed (`nXSize x nYSize`).

The `nPixelSpace`, `nLineSpace` and `nBandSpace` parameters allow reading into
or writing from various organization of buffers.

For highest performance full resolution data access, read and write on "block
boundaries" as returned by `GetBlockSize()`, or use the `ReadBlock()` and
`WriteBlock()` methods.

### Parameters
* `eRWFlag`     Either `GF_Read` to read a region of data, or `GF_Write` to
write a region of data.
* `nXOff`       The pixel offset to the top left corner of the region of the
band to be accessed. This would be zero to start from the left side.
* `nYOff`       The line offset to the top left corner of the region of the
band to be accessed. This would be zero to start from the top.
* `nXSize`      The width of the region of the band to be accessed in pixels.
* `nYSize`      The height of the region of the band to be accessed in lines.
* `pData`       The buffer into which the data should be read, or from which it
should be written. It must contain ≥`nBufXSize * nBufYSize * nBandCount` words
of type `eBufType`. It is organized in left to right, top to bottom pixel
order. Spacing is controlled by the `nPixelSpace`, and `nLineSpace` parameters
* `nBufXSize`   The width of the buffer image into which the desired region is
to be read, or from which it is to be written.
* `nBufYSize`   The height of the buffer image into which the desired region is
to be read, or from which it is to be written.
* `eBufType`    The type of the pixel values in `pBuffer`. The pixel values
will be auto-translated to/from the `GDALRasterBand` data type as needed.
* `nBandCount`  The number of bands being read or written.
* `panBandMap`  The list of nBandCount band numbers being read/written. Note
band numbers are 1 based. May be `NULL` to select the first `nBandCount` bands
* `nPixelSpace` The byte offset from the start of one pixel value in `pBuffer`
to the start of the next pixel value within a scanline. If defaulted (0) the
size of the datatype `eBufType` is used.
* `nLineSpace`  The byte offset from the start of one scanline in pBuffer to
the start of the next. If defaulted (0) the size of the datatype
`eBufType * nBufXSize` is used.
* `nBandSpace`  the byte offset from the start of one bands data to the start
of the next. If defaulted (0) the value will be `nLineSpace * nBufYSize`
implying band sequential organization of the data buffer.
* `psExtraArg`  (new in GDAL 2.0) pointer to a GDALRasterIOExtraArg structure
with additional arguments to specify resampling and progress callback, or
`NULL` for default behaviour. The `GDAL_RASTERIO_RESAMPLING` configuration
option can also be defined to override the default resampling to one of
`BILINEAR`, `CUBIC`, `CUBICSPLINE`, `LANCZOS`, `AVERAGE` or `MODE`.
"""
function rasterio!{T <: Real}(
            dataset::Dataset,
            buffer::Array{T, 3},
            bands::Vector{Cint},
            _width::Integer,
            _height::Integer,
            access::GDAL.GDALRWFlag = GDAL.GF_Read,
            xoffset::Integer = 0,
            yoffset::Integer = 0,
            nPixelSpace::Integer = 0,
            nLineSpace::Integer = 0,
            nBandSpace::Integer = 0)
    xsize, ysize, zsize = size(buffer)
    _nband = length(bands)
    @assert _nband == zsize
    GDAL.datasetrasterio(dataset, access, xoffset, yoffset, _width,
                         _height, pointer(buffer), xsize, ysize, _gdaltype[T],
                         _nband, pointer(bands), nPixelSpace, nLineSpace,
                         nBandSpace)
    buffer
end

function rasterio!{T <: Real}(
            dataset::Dataset,
            buffer::Array{T, 3},
            bands::Vector{Cint},
            access::GDAL.GDALRWFlag = GDAL.GF_Read,
            xoffset::Integer = 0,
            yoffset::Integer = 0,
            nPixelSpace::Integer = 0,
            nLineSpace::Integer = 0,
            nBandSpace::Integer = 0)
    rasterio!(dataset, buffer, bands, width(dataset), height(dataset), access,
              xoffset, yoffset, nPixelSpace, nLineSpace, nBandSpace)
end

function rasterio!{T <: Real, U <: Integer}(
            dataset::Dataset,
            buffer::Array{T, 3},
            bands::Vector{Cint},
            rows::UnitRange{U},
            cols::UnitRange{U},
            access::GDAL.GDALRWFlag = GDAL.GF_Read,
            nPixelSpace::Integer = 0,
            nLineSpace::Integer = 0,
            nBandSpace::Integer = 0)
    _width = cols[end] - cols[1] + 1
    _width < 0 && error("invalid window width")
    _height = rows[end] - rows[1] + 1
    _height < 0 && error("invalid window height")
    rasterio!(dataset, buffer, bands, _width, _height, access, cols[1], rows[1],
              nPixelSpace, nLineSpace, nBandSpace)
end

# """
#     GDALDatasetRasterIOEx(GDALDatasetH hDS,
#                           GDALRWFlag eRWFlag,
#                           int nXOff,
#                           int nYOff,
#                           int nXSize,
#                           int nYSize,
#                           void * pData,
#                           int nBufXSize,
#                           int nBufYSize,
#                           GDALDataType eBufType,
#                           int nBandCount,
#                           int * panBandMap,
#                           GSpacing nPixelSpace,
#                           GSpacing nLineSpace,
#                           GSpacing nBandSpace,
#                           GDALRasterIOExtraArg * psExtraArg) -> CPLErr
# Read/write a region of image data from multiple bands.
# """
# function datasetrasterioex{T <: GDALDatasetH}(hDS::Ptr{T},eRWFlag::GDALRWFlag,nDSXOff::Integer,nDSYOff::Integer,nDSXSize::Integer,nDSYSize::Integer,pBuffer,nBXSize::Integer,nBYSize::Integer,eBDataType::GDALDataType,nBandCount::Integer,panBandCount,nPixelSpace::GSpacing,nLineSpace::GSpacing,nBandSpace::GSpacing,psExtraArg)
#     ccall((:GDALDatasetRasterIOEx,libgdal),CPLErr,(Ptr{GDALDatasetH},GDALRWFlag,Cint,Cint,Cint,Cint,Ptr{Void},Cint,Cint,GDALDataType,Cint,Ptr{Cint},GSpacing,GSpacing,GSpacing,Ptr{GDALRasterIOExtraArg}),hDS,eRWFlag,nDSXOff,nDSYOff,nDSXSize,nDSYSize,pBuffer,nBXSize,nBYSize,eBDataType,nBandCount,panBandCount,nPixelSpace,nLineSpace,nBandSpace,psExtraArg)
# end

"""
Read/write a region of image data for this band.

This method allows reading a region of a `GDALRasterBand` into a buffer, or
writing data from a buffer into a region of a `GDALRasterBand`. It
automatically takes care of data type translation if the data type (`eBufType`)
of the buffer is different than that of the `GDALRasterBand`. The method also
takes care of image decimation / replication if the buffer size
`(nBufXSize x nBufYSize)` is different than the size of the region being
accessed `(nXSize x nYSize)`.

The `nPixelSpace` and `nLineSpace` parameters allow reading into or writing
from unusually organized buffers. This is primarily used for buffers containing
more than one bands raster data in interleaved format.

Some formats may efficiently implement decimation into a buffer by reading from
lower resolution overview images.

For highest performance full resolution data access, read and write on "block
boundaries" returned by `GetBlockSize()`, or use the `ReadBlock()` and
`WriteBlock()` methods.

### Parameters
* `eRWFlag`     Either GF_Read to read a region of data, or GF_Write to write a
region of data.
* `nXOff`       The pixel offset to the top left corner of the region of the
band to be accessed. This would be zero to start from the left side.
* `nYOff`       The line offset to the top left corner of the region of the
band to be accessed. This would be zero to start from the top.
* `nXSize`      The width of the region of the band to be accessed in pixels.
* `nYSize`      The height of the region of the band to be accessed in lines.
* `pData`       The buffer into which the data should be read, or from which it
should be written. This buffer must contain at least `(nBufXSize * nBufYSize)`
words of type `eBufType`. It is organized in left to right, top to bottom pixel
order. Spacing is controlled by the `nPixelSpace`, and `nLineSpace` parameters.
* `nBXSize`     The width of the buffer image into which the desired region is
to be read, or from which it is to be written.
* `nBYSize`     The height of the buffer image into which the desired region is
to be read, or from which it is to be written.
* `eBufType`    The type of the pixel values in the `buffer`. The pixel values
will be auto-translated to/from the `GDALRasterBand` data type as needed.
* `nPixelSpace` The byte offset from the start of one pixel value in `buffer`
to the start of the next pixel value within a scanline. If defaulted (0) the
size of the datatype `eBufType` is used.
* `nLineSpace`  The byte offset from the start of one scanline in `buffer` to
the start of the next. If defaulted (0) the size of the datatype
`(eBufType * nBufXSize)` is used.
* `psExtraArg`  (new in GDAL 2.0) pointer to a GDALRasterIOExtraArg structure
with additional arguments to specify resampling and progress callback, or
`NULL` for default behaviour. The `GDAL_RASTERIO_RESAMPLING` configuration
option can also be defined to override the default resampling to one of
`BILINEAR`, `CUBIC`, `CUBICSPLINE`, `LANCZOS`, `AVERAGE` or `MODE`.

### Returns
`CE_Failure` if the access fails, otherwise `CE_None`.
"""
function rasterio!{T <: Real}(rasterband::RasterBand,
                              buffer::Array{T,2},
                              _width::Integer,
                              _height::Integer,
                              xoffset::Integer = 0,
                              yoffset::Integer = 0,
                              access::GDAL.GDALRWFlag = GDAL.GF_Read,
                              nPixelSpace::Integer = 0,
                              nLineSpace::Integer = 0)
    xsize, ysize = size(buffer)
    GDAL.rasterio(rasterband, access, xoffset, yoffset, _width, _height,
                  pointer(buffer), xsize, ysize, _gdaltype[T],
                  nPixelSpace, nLineSpace)
    buffer
end

function rasterio!{T <: Real}(rasterband::RasterBand,
                              buffer::Array{T,2},
                              xoffset::Integer = 0,
                              yoffset::Integer = 0,
                              access::GDAL.GDALRWFlag = GDAL.GF_Read,
                              nPixelSpace::Integer = 0,
                              nLineSpace::Integer = 0)
    rasterio!(rasterband, buffer, width(rasterband), height(rasterband),
              xoffset, yoffset, access, nPixelSpace, nLineSpace)
end

function rasterio!{T <: Real, U <: Integer}(rasterband::RasterBand,
                                            buffer::Array{T,2},
                                            rows::UnitRange{U},
                                            cols::UnitRange{U},
                                            access::GDAL.GDALRWFlag = GDAL.GF_Read,
                                            nPixelSpace::Integer = 0,
                                            nLineSpace::Integer = 0)
    _width = length(cols)
    _width < 1 && error("invalid window width")
    _height = length(rows)
    _height < 1 && error("invalid window height")
    rasterio!(rasterband, buffer, _width, _height, cols[1]-1, rows[1]-1, access,
              nPixelSpace, nLineSpace)
end

# """
#     GDALRasterIOEx(GDALRasterBandH hBand,
#                    GDALRWFlag eRWFlag,
#                    int nXOff,
#                    int nYOff,
#                    int nXSize,
#                    int nYSize,
#                    void * pData,
#                    int nBufXSize,
#                    int nBufYSize,
#                    GDALDataType eBufType,
#                    GSpacing nPixelSpace,
#                    GSpacing nLineSpace,
#                    GDALRasterIOExtraArg * psExtraArg) -> CPLErr
# Read/write a region of image data for this band.
# """
# function rasterioex{T <: GDALRasterBandH}(hRBand::Ptr{T},eRWFlag::GDALRWFlag,nDSXOff::Integer,nDSYOff::Integer,nDSXSize::Integer,nDSYSize::Integer,pBuffer,nBXSize::Integer,nBYSize::Integer,eBDataType::GDALDataType,nPixelSpace::GSpacing,nLineSpace::GSpacing,psExtraArg)
#     ccall((:GDALRasterIOEx,libgdal),CPLErr,(Ptr{GDALRasterBandH},GDALRWFlag,Cint,Cint,Cint,Cint,Ptr{Void},Cint,Cint,GDALDataType,GSpacing,GSpacing,Ptr{GDALRasterIOExtraArg}),hRBand,eRWFlag,nDSXOff,nDSYOff,nDSXSize,nDSYSize,pBuffer,nBXSize,nBYSize,eBDataType,nPixelSpace,nLineSpace,psExtraArg)
# end

function fetch!{T <: Real}(band::RasterBand, buffer::Array{T,2})
    (band == C_NULL) && error("Can't read invalid rasterband")
    rasterio!(band, buffer, GDAL.GF_Read)
end

function fetch!{T <: Real}(band::RasterBand,
                           buffer::Array{T,2},
                           _width::Integer,
                           _height::Integer,
                           xoffset::Integer,
                           yoffset::Integer)
    (band == C_NULL) && error("Can't read invalid rasterband")
    rasterio!(band, buffer, _width, _height, xoffset, yoffset)
end

function fetch!{T <: Real, U <: Integer}(band::RasterBand,
                                         buffer::Array{T,2},
                                         rows::UnitRange{U},
                                         cols::UnitRange{U})
    (band == C_NULL) && error("Can't read invalid rasterband")
    rasterio!(band, buffer, rows, cols)
end

function fetch(band::RasterBand)
    (band == C_NULL) && error("Can't read invalid rasterband")
    rasterio!(band, Array(getdatatype(band), width(band), height(band)))
end

function fetch(band::RasterBand,
               _width::Integer,
               _height::Integer,
               xoffset::Integer,
               yoffset::Integer)
    (band == C_NULL) && error("Can't read invalid rasterband")
    buffer = Array(getdatatype(band), width(band), height(band))
    rasterio!(band, buffer, _width, _height, xoffset, yoffset)
end

function fetch{U <: Integer}(band::RasterBand,
                             rows::UnitRange{U},
                             cols::UnitRange{U})
    (band == C_NULL) && error("Can't read invalid rasterband")
    buffer = Array(getdatatype(band), length(cols), length(rows))
    rasterio!(band, buffer, rows, cols)
end

function update!{T <: Real}(band::RasterBand, buffer::Array{T,2})
    (band == C_NULL) && error("Can't write invalid rasterband")
    rasterio!(band, buffer, GDAL.GF_Write)
end

function update!{T <: Real}(band::RasterBand,
                            buffer::Array{T,2},
                            _width::Integer,
                            _height::Integer,
                            xoffset::Integer,
                            yoffset::Integer)
    (band == C_NULL) && error("Can't write invalid rasterband")
    rasterio!(band, buffer, _width, _height,
              xoffset, yoffset, GDAL.GF_Write)
end

function update!{T <: Real, U <: Integer}(band::RasterBand,
                                          buffer::Array{T,2},
                                          rows::UnitRange{U},
                                          cols::UnitRange{U})
    (band == C_NULL) && error("Can't write invalid rasterband")
    rasterio!(band, buffer, rows, cols, GDAL.GF_Write)
end

function fetch!{T <: Real}(dataset::Dataset, buffer::Array{T,2}, i::Integer)
    (dataset == C_NULL) && error("Can't read closed dataset")
    fetch!(fetchband(dataset, i), buffer)
end

function fetch!{T <: Real}(dataset::Dataset,
                           buffer::Array{T,3},
                           indices::Vector{Cint})
    (dataset == C_NULL) && error("Can't read closed dataset")
    rasterio!(dataset, buffer, indices, GDAL.GF_Read)
end

function fetch!{T <: Real}(dataset::Dataset, buffer::Array{T,3})
    (dataset == C_NULL) && error("Can't read closed dataset")
    _nband = nraster(dataset)
    @assert size(buffer, 3) == _nband
    rasterio!(dataset, buffer, collect(Cint, 1:_nband), GDAL.GF_Read)
end

function fetch!{T <: Real}(dataset::Dataset,
                           buffer::Array{T,2},
                           i::Integer,
                           _width::Integer,
                           _height::Integer,
                           xoffset::Integer,
                           yoffset::Integer)
    (dataset == C_NULL) && error("Can't read closed dataset")
    fetch!(fetchband(dataset, i), buffer, _width, _height, xoffset, yoffset)
end

function fetch!{T <: Real}(dataset::Dataset,
                           buffer::Array{T,3},
                           indices::Vector{Cint},
                           _width::Integer,
                           _height::Integer,
                           xoffset::Integer,
                           yoffset::Integer)
    (dataset == C_NULL) && error("Can't read closed dataset")
    rasterio!(dataset, buffer, indices, _width, _height, xoffset, yoffset)
end

function fetch!{T <: Real, U <: Integer}(dataset::Dataset,
                                         buffer::Array{T,2},
                                         i::Integer,
                                         rows::UnitRange{U},
                                         cols::UnitRange{U})
    (dataset == C_NULL) && error("Can't read closed dataset")
    fetch!(fetchband(dataset, i), buffer, rows, cols)
end

function fetch!{T <: Real, U <: Integer}(dataset::Dataset,
                                         buffer::Array{T,3},
                                         indices::Vector{Cint},
                                         rows::UnitRange{U},
                                         cols::UnitRange{U})
    (dataset == C_NULL) && error("Can't read closed dataset")
    rasterio!(dataset, buffer, indices, rows, cols)
end

function fetch(dataset::Dataset, i::Integer)
    (dataset == C_NULL) && error("Can't read closed dataset")
    fetch(fetchband(dataset, i))
end

function fetch(dataset::Dataset, indices::Vector{Cint})
    (dataset == C_NULL) && error("Can't read closed dataset")
    buffer = Array(getdatatype(fetchband(dataset, indices[1])),
                   width(dataset), height(dataset), length(indices))
    rasterio!(dataset, buffer, indices)
end

function fetch(dataset::Dataset)
    (dataset == C_NULL) && error("Can't read closed dataset")
    buffer = Array(getdatatype(fetchband(dataset, 1)),
                   width(dataset), height(dataset), nraster(dataset))
    fetch!(dataset, buffer)
end

function fetch(dataset::Dataset,
               i::Integer,
               _width::Integer,
               _height::Integer,
               xoffset::Integer,
               yoffset::Integer)
    (dataset == C_NULL) && error("Can't read closed dataset")
    fetch(fetchband(dataset, i), _width, _height, xoffset, yoffset)
end

function fetch{T <: Integer}(dataset::Dataset,
                             indices::Vector{T},
                             _width::Integer,
                             _height::Integer,
                             xoffset::Integer,
                             yoffset::Integer)
    (dataset == C_NULL) && error("Can't read closed dataset")
    buffer = Array(getdatatype(fetchband(dataset, indices[1])),
                   width(dataset), height(dataset), length(indices))
    rasterio!(dataset, buffer, indices, _width, _height, xoffset, yoffset)
end

function fetch{U <: Integer}(dataset::Dataset,
                             i::Integer,
                             rows::UnitRange{U},
                             cols::UnitRange{U})
    (dataset == C_NULL) && error("Can't read closed dataset")
    fetch(fetchband(dataset, i), rows, cols)
end

function fetch{U <: Integer}(dataset::Dataset,
                             indices::Vector{Cint},
                             rows::UnitRange{U},
                             cols::UnitRange{U})
    (dataset == C_NULL) && error("Can't read closed dataset")
    buffer = Array(getdatatype(fetchband(dataset, indices[1])),
                   width(dataset), height(dataset), length(indices))
    rasterio!(dataset, buffer, indices, rows, cols)
end

function update!{T <: Real}(dataset::Dataset, buffer::Array{T,2}, i::Integer)
    (dataset == C_NULL) && error("Can't write closed dataset")
    update!(fetchband(dataset, i), buffer)
end

function update!{T <: Real}(dataset::Dataset,
                            buffer::Array{T,3},
                            indices::Vector{Cint})
    (dataset == C_NULL) && error("Can't write closed dataset")
    rasterio!(dataset, buffer, indices, GDAL.GF_Write)
end

function update!{T <: Real}(dataset::Dataset,
                            buffer::Array{T,2},
                            i::Integer,
                            _width::Integer,
                            _height::Integer,
                            xoffset::Integer,
                            yoffset::Integer)
    (dataset == C_NULL) && error("Can't write closed dataset")
    update!(fetchband(dataset, i), buffer, _width, _height,
            xoffset, yoffset)
end

function update!{T <: Real}(dataset::Dataset,
                            buffer::Array{T,3},
                            indices::Vector{Cint},
                            _width::Integer,
                            _height::Integer,
                            xoffset::Integer,
                            yoffset::Integer)
    (dataset == C_NULL) && error("Can't write closed dataset")
    rasterio!(dataset, buffer, indices, _width, _height,
              xoffset, yoffset, GDAL.GF_Write)
end

function update!{T <: Real, U <: Integer}(dataset::Dataset,
                                          buffer::Array{T,2},
                                          i::Integer,
                                          rows::UnitRange{U},
                                          cols::UnitRange{U})
    (dataset == C_NULL) && error("Can't write closed dataset")
    update!(fetchband(dataset, i), buffer, rows, cols)
end

function update!{T <: Real, U <: Integer}(dataset::Dataset,
                                          buffer::Array{T,3},
                                          indices::Vector{Cint},
                                          rows::UnitRange{U},
                                          cols::UnitRange{U})
    (dataset == C_NULL) && error("Can't write closed dataset")
    rasterio!(dataset, buffer, indices, rows, cols, GDAL.GF_Write)
end


text:
module TheoCrystal

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Theo Crystal" => Ahorn.EntityPlacement(
        Maple.TheoCrystal
    )
)

sprite = "characters/theoCrystal/idle00.png"

function Ahorn.selection(entity::Maple.TheoCrystal)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y - 10)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.TheoCrystal, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, -10)

end


text:
test_file_dir() = joinpath(@__DIR__, "test_files");

test_fig_path(fName :: String) = joinpath(test_file_dir(), fName * ".pdf");

function save_test_figure(p, fName :: String)
    isdir(test_file_dir())  ||  mkpath(test_file_dir());
    savefig(p, test_fig_path(fName));
    println("Saved test figure ", fName);
end

delete_test_figure(fName :: String) =
    isfile(test_fig_path(fName))  &&  rm(test_fig_path(fName));

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


text:
module Viznet

using Compose

include("ordered_dict.jl")
include("mathjax.jl")
include("lattices/lattices.jl")
#include("compose.jl")
include("brush.jl")
include("styles.jl")
include("intersection.jl")
include("deprecations.jl")
end # module


text:
module BoundingBoxes

import LinearAlgebra: norm

const ZERO_FLOAT32 = zero(Float32)

export BoundingBox

mutable struct BoundingBox
    minCorner   ::NTuple{3,Float32}
    maxCorner   ::NTuple{3,Float32}
end

function BoundingBox()
    minCorner = (Inf32, Inf32, Inf32)
    maxCorner = (ZERO_FLOAT32, ZERO_FLOAT32, ZERO_FLOAT32)
    BoundingBox( minCorner, maxCorner )
end

function BoundingBox( minCorner::NTuple{3,Int}, maxCorner::NTuple{3,Int} )
	BoundingBox( 	map(x->convert(Float32,x), minCorner),
					map(x->convert(Float32,x), maxCorner) )
end

function BoundingBox( minCorner::Vector{Float32}, maxCorner::Vector{Float32} )
    BoundingBox((minCorner...,), (maxCorner...,))
end

"""
get bounding box from a node list
"""
function BoundingBox(nodeList::Union{Vector, Set})
    minCorner = fill(Inf32,3)
    maxCorner = fill(ZERO_FLOAT32, 3)
    for node in nodeList
        minCorner = map(min, minCorner, node[1:3])
        maxCorner = map(max, maxCorner, node[1:3])
    end
    BoundingBox(minCorner, maxCorner)
end

function get_unit_range( self::BoundingBox )
    map((x,y)->floor(Int,x):ceil(Int,y), self.minCorner, self.maxCorner)
end

function Base.size(self::BoundingBox)
    map(length, get_unit_range(self))
end

function Base.isequal(self::BoundingBox, other::BoundingBox)
    self.minCorner==other.minCorner && self.maxCorner==other.maxCorner
end

function Base.:(==)(self::BoundingBox, other::BoundingBox)
    isequal(self, other)
end

"""
    Base.union(self::BoundingBox, other::BoundingBox)
get the bounding box containing these two bounding boxes
"""
function Base.union(self::BoundingBox, other::BoundingBox)
    minCorner = map(min, self.minCorner, other.minCorner)
    maxCorner = map(max, self.maxCorner, other.maxCorner)
    BoundingBox(minCorner, maxCorner)
end

function isinside(self::BoundingBox, point::Union{Tuple, Vector})
    all( map((x,y)->x>y, self.maxCorner, point[1:3]) ) &&
    all( map((x,y)->x<y, self.minCorner, point[1:3]) )
end

"""
    to_voxel_space(self::BoundingBox, voxelSize::Union{Tuple, Vector})

"""
function to_voxel_space(self::BoundingBox, voxelSize::Union{Tuple, Vector})
    minCorner = map(fld, self.minCorner, voxelSize)
    maxCorner = map(cld, self.maxCorner, voxelSize)
    BoundingBox(minCorner, maxCorner)
end

"""
    distance_from(self::BoundingBox, point::Union{Tuple, Vector})

compute the distance from bounding box using a smart way
https://stackoverflow.com/questions/5254838/calculating-distance-between-a-point-and-a-rectangular-box-nearest-point

code in JS
function distance(rect, p) {
  var dx = Math.max(rect.min.x - p.x, 0, p.x - rect.max.x);
  var dy = Math.max(rect.min.y - p.y, 0, p.y - rect.max.y);
  return Math.sqrt(dx*dx + dy*dy);
}
"""
function distance_from(self::BoundingBox, point::Union{Tuple, Vector})
    @assert length(point) == 3 || length(point)==4
    d = map((cmin, p, cmax)-> max(cmin-p, Float32(0), p-cmax),
                            self.minCorner, point[1:3], self.maxCorner )
    norm([d...])
end

end # module


text:
module Nilpotent

export He, BHPTile, a, b, bhp_subst

using StructEquality
using Base.Iterators

QQ = Rational{Int}

using ...CoreDefs

@def_structequal struct He <: DGroupElem
    a::QQ
    b::QQ
    c::QQ
end

@enum BHPTile begin
    A=1
    B=2
end

function Base.:*(x :: He, y :: He)
    return He( x.a + y.a, x.b + y.b, x.c + y.c + (x.a*y.b-y.a*x.b)//2)
end

function CoreDefs.dilate(λ, x :: He)
    return He(λ*x.a, λ * x.b, λ^2 * x.c )
end

function Base.inv(x :: He)
    return He(-x.a, -x.b, -x.c)
end

function CoreDefs.id(::Type{He})
    return He(0,0,0)
end

function a(a,b,c)
    return He(a,b,c) => A
end
function a()
    return a(0,0,0)
end
function b(a,b,c)
    return He(a,b,c) => B
end
function b()
    return b(0,0,0)
end

function bhp_subst()
    heisenberg_subst = Dict([
        (A, [
            a(-2,-2,-8),
            a(0,-2,-8),
            a(2,-2,-8),
            a(-2,0,-8),
            a(0,0,-8),
            a(2,0,-8),
            a(-2,2,-8),
            a(0,2,-8),
            a(2,2,-8),
            a(-2,-2,-6),
            b(0,-2,-6),
            b(2,-2,-6),
            a(-2,0,-6),
            a(0,0,-6),
            a(2,0,-6),
            b(-2,2,-6),
            a(0,2,-6),
            b(2,2,-6),
            b(-2,-2,-4),
            a(0,-2,-4),
            a(2,-2,-4),
            b(-2,0,-4),
            a(0,0,-4),
            b(2,0,-4),
            a(-2,2,-4),
            a(0,2,-4),
            a(2,2,-4),
            a(-2,-2,-2),
            a(0,-2,-2),
            a(2,-2,-2),
            a(-2,0,-2),
            a(0,0,-2),
            a(2,0,-2),
            a(-2,2,-2),
            a(0,2,-2),
            a(2,2,-2),
            b(-2,-2,0),
            b(0,-2,0),
            b(2,-2,0),
            b(-2,0,0),
            a(0,0,0),
            b(2,0,0),
            b(-2,2,0),
            b(0,2,0),
            b(2,2,0),
            a(-2,-2,2),
            a(0,-2,2),
            a(2,-2,2),
            a(-2,0,2),
            a(0,0,2),
            a(2,0,2),
            a(-2,2,2),
            a(0,2,2),
            a(2,2,2),
            a(-2,-2,4),
            a(0,-2,4),
            b(2,-2,4),
            a(-2,0,4),
            a(0,0,4),
            a(2,0,4),
            a(-2,2,4),
            a(0,2,4),
            b(2,2,4),
            b(-2,-2,6),
            b(0,-2,6),
            a(2,-2,6),
            b(-2,0,6),
            a(0,0,6),
            b(2,0,6),
            b(-2,2,6),
            a(0,2,6),
            a(2,2,6),
            a(-2,-2,8),
            a(0,-2,8),
            a(2,-2,8),
            a(-2,0,8),
            a(0,0,8),
            a(2,0,8),
            a(-2,2,8),
            a(0,2,8),
            a(2,2,8),
        ]),
        (B, [
            a(-2,-2,-8),
            a(0,-2,-8),
            a(2,-2,-8),
            a(-2,0,-8),
            a(0,0,-8),
            a(2,0,-8),
            a(-2,2,-8),
            a(0,2,-8),
            a(2,2,-8),
            a(-2,-2,-6),
            b(0,-2,-6),
            b(2,-2,-6),
            b(-2,0,-6),
            a(0,0,-6),
            a(2,0,-6),
            b(-2,2,-6),
            a(0,2,-6),
            b(2,2,-6),
            b(-2,-2,-4),
            a(0,-2,-4),
            a(2,-2,-4),
            a(-2,0,-4),
            a(0,0,-4),
            b(2,0,-4),
            a(-2,2,-4),
            a(0,2,-4),
            a(2,2,-4),
            a(-2,-2,-2),
            a(0,-2,-2),
            a(2,-2,-2),
            a(-2,0,-2),
            a(0,0,-2),
            a(2,0,-2),
            a(-2,2,-2),
            a(0,2,-2),
            a(2,2,-2),
            b(-2,-2,0),
            b(0,-2,0),
            b(2,-2,0),
            b(-2,0,0),
            a(0,0,0),
            b(2,0,0),
            b(-2,2,0),
            b(0,2,0),
            b(2,2,0),
            a(-2,-2,2),
            a(0,-2,2),
            a(2,-2,2),
            a(-2,0,2),
            a(0,0,2),
            a(2,0,2),
            a(-2,2,2),
            a(0,2,2),
            a(2,2,2),
            a(-2,-2,4),
            b(0,-2,4),
            a(2,-2,4),
            a(-2,0,4),
            a(0,0,4),
            a(2,0,4),
            b(-2,2,4),
            a(0,2,4),
            b(2,2,4),
            b(-2,-2,6),
            a(0,-2,6),
            b(2,-2,6),
            b(-2,0,6),
            a(0,0,6),
            b(2,0,6),
            a(-2,2,6),
            a(0,2,6),
            a(2,2,6),
            a(-2,-2,8),
            a(0,-2,8),
            a(2,-2,8),
            a(-2,0,8),
            a(0,0,8),
            a(2,0,8),
            a(-2,2,8),
            a(0,2,8),
            a(2,2,8),
        ]),
    ])

    return SubSystem(heisenberg_subst, 3//1)
end

function CoreDefs.collar_in(tiling, t ::He)
    collar_shape = t*[He(x,y,z) for x=[-2,0,2] for y=[-2,0,2] for z=[-2,0,2]]
    collar = Dict([])
    tiling_dict = Dict(tiling)
    for s in collar_shape
        if !haskey(tiling_dict, s)
            throw(UnrecognizedCollar)
        end
        collar[s] = tiling_dict[s]
    end
    return collar
end

function CoreDefs.is_interior(tiling :: Dict, t :: He)
    collar_shape = t*[He(x,y,z) for x=[-2,0,2] for y=[-2,0,2] for z=[-2,0,2]]
    return all(g -> haskey(tiling, g), collar_shape)
end

end


text:
const CHAIN_INDEX_DOC = """
`chain_index::Vector{Int}`: An optional vector of integers specifying which chain each step
belongs to. For instance, `chain_index[step]` should return `2` if `log_likelihood[:, step]`
belongs to the second chain.
"""

const DATA_ARG = """
`data`: An array of data points used to estimate the parameters of the model.
"""

const LIKELIHOOD_FUNCTION_ARG = """
`ll_fun::Function`: A function taking a single data point and returning the log-likelihood
of that point. This function must take the form `f(θ[1], ..., θ[n], data)`, where `θ` is the
parameter vector. See also the `splat` keyword argument.
"""

const LOG_LIK_ARR = """
`log_likelihood::Array`: A matrix or 3d array of log-likelihood values indexed as
`[data, step, chain]`. The chain argument can be left off if `chain_index` is provided
or if all posterior samples were drawn from a single chain.
"""

const R_EFF_DOC = """
`r_eff::AbstractVector`: An (optional) vector of relative effective sample sizes used in ESS
calculations. If left empty, calculated automatically using the FFTESS method from
InferenceDiagnostics.jl. See `relative_eff` to calculate these values.
"""

const ARGS = """`args...`: Positional arguments to be passed to"""

const KWARGS = """`kwargs...`: Keyword arguments to be passed to"""

###############
## FUNCTIONS ##
###############

"""
Throw a warning if any `pareto_k` values exceed 0.7.
"""
function _throw_pareto_k_warning(ξ)
    if any(ξ .≥ 1)
        @warn "Some Pareto k values are extremely high (>1). PSIS will not produce " *
              "consistent estimates."
    elseif any(ξ .≥ 0.7)
        @warn "Some Pareto k values are high (>.7), indicating PSIS has failed to " *
              "approximate the true distribution."
    elseif any(ξ .≥ 0.5)
        @info "Some Pareto k values are slightly high (>.5); some pointwise estimates " *
              "may be slow to converge or have high variance."
    end
end

"""
Assume that all objects belong to a single chain if chain index is missing. Inform user.
"""
function _assume_one_chain(matrix)
    @info "Chain information was not provided; " *
          "all samples are assumed to be drawn from a single chain."
    return ones(length(matrix))
end

"""
Convert a matrix+chain_index representation to a 3d array representation.
"""
function _convert_to_array(matrix::AbstractMatrix, chain_index::AbstractVector)
    indices = unique(chain_index)
    biggest_idx = maximum(indices)
    dims = size(matrix)
    if dims[2] ≠ length(chain_index)
        throw(ArgumentError("Some entries do not have a chain index."))
    elseif !issetequal(indices, 1:biggest_idx)
        throw(
            ArgumentError(
                "Indices must be numbered from 1 through the total number of chains."
            )
        )
    else
        # Check how many elements are in each chain, assign to "counts"
        counts = count.(eachslice(chain_index .== indices'; dims=2))
        # check if all inputs are the same length
        if !all(==(counts[1]), counts)
            throw(ArgumentError("All chains must be of equal length."))
        end
    end
    new_ratios = similar(matrix, dims[1], dims[2] ÷ biggest_idx, biggest_idx)
    for i in 1:biggest_idx
        new_ratios[:, :, i] .= matrix[:, chain_index .== i]
    end
    return new_ratios
end

max_stars_repo_name: ParadaCarleton/JuLOOa

max_issues_repo_name: ParadaCarleton/JuLOOa

max_forks_repo_name: ParadaCarleton/JuLOOa


text:
using Distributed
using CalibrateEDMF
using CalibrateEDMF.ReferenceModels
using CalibrateEDMF.ReferenceStats
using CalibrateEDMF.TurbulenceConvectionUtils
using TurbulenceConvection
using CalibrateEDMF.HelperFuncs

export run_SCM_parallel, eval_single_ref_model, versioned_model_eval_parallel

function run_SCM_parallel(
    ME::ModelEvaluator;
    error_check::Bool = false,
    namelist_args = nothing,
    failure_handler = "high_loss",
) where {FT <: Real}
    return run_SCM_parallel(
        ME.param_cons,
        ME.param_names,
        ME.ref_models,
        ME.ref_stats,
        error_check = error_check,
        namelist_args = namelist_args,
        failure_handler = failure_handler,
    )
end

function run_SCM_parallel(
    u::Vector{FT},
    u_names::Vector{String},
    RM::Vector{ReferenceModel},
    RS::ReferenceStatistics;
    error_check::Bool = false,
    namelist_args = nothing,
    failure_handler = "high_loss",
) where {FT <: Real}

    mkpath(joinpath(pwd(), "tmp"))
    result_arr = pmap(x -> eval_single_ref_model(x..., RS, u, u_names, namelist_args), enumerate(RM))
    # Unpack
    sim_dirs, g_scm, g_scm_pca, sim_errors = (getindex.(result_arr, i) for i in 1:4)
    g_scm = vcat(g_scm...)
    g_scm_pca = vcat(g_scm_pca...)
    model_error = any(sim_errors)

    # penalize nan-values in output
    if model_error && failure_handler == "high_loss"
        g_scm .= 1e5
        g_scm_pca .= 1e5
    end

    @info "Length of g_scm (full, pca): ($(length(g_scm)), $(length(g_scm_pca)))"
    if error_check
        return sim_dirs, g_scm, g_scm_pca, model_error
    else
        return sim_dirs, g_scm, g_scm_pca
    end
end

function eval_single_ref_model(
    m_index::IT,
    m::ReferenceModel,
    RS::ReferenceStatistics,
    u::Vector{FT},
    u_names::Vector{String},
    namelist_args = nothing,
) where {FT <: Real, IT <: Int}

    # create temporary directory to store SCM data in
    tmpdir = mktempdir(joinpath(pwd(), "tmp"))
    # run TurbulenceConvection.jl. Get output directory for simulation data
    sim_dir, model_error = run_SCM_handler(m, tmpdir, u, u_names, namelist_args)
    filename = get_stats_path(sim_dir)
    if model_error
        g_scm = fill(NaN, length(get_z_obs(m)) * length(m.y_names))
    else
        g_scm = get_profile(m, filename, z_scm = get_z_obs(m))
        g_scm = normalize_profile(g_scm, length(m.y_names), RS.norm_vec[m_index])
    end
    # perform PCA reduction
    g_scm_pca = RS.pca_vec[m_index]' * g_scm
    return sim_dir, g_scm, g_scm_pca, model_error
end

function versioned_model_eval_parallel(
    version::Union{String, Int},
    outdir_path::String,
    mode::String,
    config::Dict{Any, Any},
)
    @assert mode in ["train", "validation"]
    # Omits validation if unsolicited
    if mode == "validation" && isnothing(get(config, "validation", nothing))
        return
    elseif mode == "validation"
        input_path = scm_val_init_path(outdir_path, version)
        output_path = scm_val_output_path(outdir_path, version)
    else
        input_path = scm_init_path(outdir_path, version)
        output_path = scm_output_path(outdir_path, version)
    end
    # Load inputs
    scm_args = load(input_path)
    namelist_args = get_entry(config["scm"], "namelist_args", nothing)
    model_evaluator = scm_args["model_evaluator"]
    batch_indices = scm_args["batch_indices"]
    # Eval
    failure_handler = get_entry(config["process"], "failure_handler", "high_loss")

    sim_dirs, g_scm, g_scm_pca =
        run_SCM_parallel(model_evaluator, namelist_args = namelist_args, failure_handler = failure_handler)

    # Store output and delete input
    jldsave(output_path; sim_dirs, g_scm, g_scm_pca, model_evaluator, version, batch_indices)
    rm(input_path)
end


text:
sourcecode = readlines("a18.txt")

function shuntyard(line, prio)
  ops = ['(']
  out = []
  function shunt(op)
    op == ' ' && return
    while op != '(' && first(ops) != '(' && (op == ')' || prio(ops[1]) <= prio(op))
      pushfirst!(out, Expr(:call, Symbol(popfirst!(ops)), popfirst!(out), popfirst!(out)))
    end
    if (op == ')')
      popfirst!(ops)
    else
      pushfirst!(ops, op)
    end
  end
  numberat = findfirst(r"\d+", line)
  while !isnothing(numberat)
    if first(numberat) != 1
      shunt.(SubString(line, 1, first(numberat)-1) |> collect)
    end
    pushfirst!(out, parse(Int, SubString(line, first(numberat), last(numberat))))
    line = SubString(line, last(numberat)+1)
    numberat = findfirst(r"\d+", line)
  end
  shunt.(line |> collect)
  shunt(')')
  out[1]
end

function equalprio(op)
  1
end

#part 1
println("$(sum(eval.(shuntyard.(sourcecode, equalprio))))")

function plusfirst(op)
  if op == '+'
    1
  elseif op == '*'
    2
  end
end

#part 2
println("$(sum(eval.(shuntyard.(sourcecode, plusfirst))))")


text:
using Images, ImageView, ImageFiltering
br = load("bridge.jpg")
imshow(br)
#= You can find misty.jpg and bridge.jpg in the resources folder,
make sure to put it in the same directory you run the code =#

misty = load("misty.jpg")
imshow(misty)
channel_view = channelview(misty); #converts misty to an array that can be modified
channel_view = permuteddimsview(channel_view, (2, 3, 1)); #changes channel locations for convenience: height, width, rgb

x_coords = 1:599 #you can specify any range of x coords here, since width is 1200 this is right half of image now
x_coords_2 = 600:1200

y_coords = 1:799
y_coords_2 = 800:1600

channel_view[y_coords, x_coords, 1] = min.(channel_view[y_coords, x_coords, 1] .* 1.8, 1);
# multiplies red pixels (top left) by 1.8 but restricts their maximum value to 1, using the min function.

channel_view[y_coords, x_coords_2 , 2] = min.(channel_view[y_coords, x_coords_2, 2] .* 1.2, 1);
# multiplies green pixels (top right) by 1.2 but restricts their maximum value to 1, using the min function

channel_view[y_coords_2, x_coords_2, 3] = min.(channel_view[y_coords_2, x_coords_2, 3] .* 1.4, 1);
# multiplies blue pixels (bottom right) in  by 1.4 but restricts their maximum value to 1, using the min function

# bottom left of image is left untouched

imshow(misty)

# individual pixel values can be modified like this to incrase saturation.
# now let's do something even cooler!
# we can create a color filter using the same principle

img = load("bridge.jpg")
gray_img = RGB.(Gray.(img)) # creates a gray img

channel_view = channelview(img);
channel_view_gray = channelview(gray_img);

img_array = permuteddimsview(channel_view, (2, 3, 1)); # creates array of img with height, width, rgb
gray_img_array = permuteddimsview(channel_view_gray, (2, 3, 1)); # creates array of gray img with height, width, grayvalue

mask = fill(true, size(img)); # creates a grid of true values with the size of img

retained_colors_x = 500:508; # the pixels that are in this range are going to be
retained_colors_y = 400:470; # the colors that will stay in the final image (as a range, not specific values)

# img[retained_colors_x, retained_colors_y] # uncomment this to see the selected color range

for channel_id = 1:3

    # views the image array only in specific channel (for loop)
    current_channel = view(img_array, :, :, channel_id)
    # takes the retained color part only in specified channel
    retained_color_area = current_channel[retained_colors_y, retained_colors_x, :]

    #gets minimum and maximum channel values from the retained color area
    channel_min = minimum(retained_color_area)
    channel_max = maximum(retained_color_area)

    # only the pixels from the full image that are in the retained color area return true to the image mask
    channel_mask = channel_min .< current_channel .< channel_max
    mask = mask .& channel_mask
end
# only the pixels that correspond to true in the general mask retain their pixel values
# the rest is replaced with the gray image
img_masked = img_array .* mask .+ gray_img_array .* .~(mask);
# converts the array to an image that we can view (also puts the channels back into place)
final_image = colorview(RGB, permutedims(img_masked, (3, 1, 2)))

imshow(final_image)

max_stars_repo_name: gunalpms/imageprojl

max_issues_repo_name: gunalpms/imageprojl

max_forks_repo_name: gunalpms/imageprojl


text:
function car_false_positive(false_positive_rate::Float64, roadway::Roadway, scene::Scene, obstacles::Vector{ConvexPolygon}, rng::AbstractRNG)
    FP = rand(rng) < false_positive_rate
    if FP
        return random_car(roadway, scene, rng)
    else
        return nothing
    end
end

function ped_false_positive(false_positive_rate::Float64, roadway::Roadway, scene::Scene, obstacles::Vector{ConvexPolygon}, rng::AbstractRNG)
    FP = rand(rng) < false_positive_rate
    if FP
        return random_pedestrian(roadway, scene, rng)
    else
        return nothing
    end
end

function random_pedestrian(roadway::Roadway, scene::Scene, rng::AbstractRNG, ped_max_speed::Float64=2.0, ped_def::VehicleDef= VehicleDef(AgentClass.PEDESTRIAN, 1., 1.) )
    crosswalk = rand(rng, get_crosswalks(roadway))

    # position along the crosswalk
    t0 = rand(rng, Uniform(-crosswalk.width/2, crosswalk.width/2))
    s0 = rand(rng, Uniform(0., crosswalk.curve[end].s))
    ϕ0 = rand(rng, [0.,float(π)])

    #velocity
    v0 = rand(rng, Uniform(0., ped_max_speed))
    posF = Frenet(crosswalk, s0, t0, ϕ0)

    ped_state = VehicleState(posF, roadway, v0);

    return Vehicle(ped_state, ped_def, 101)
end

function random_car(roadway::Roadway, scene::Scene, rng::AbstractRNG, car_max_speed::Float64=2.0, car_def::VehicleDef=VehicleDef())
    lane = rand(rng, get_car_lanes(roadway))

    s0 = rand(rng, Uniform(0., lane.curve[end].s))
    v0 = rand(rng, Uniform(0., car_max_speed))

    car_state = VehicleState(Frenet(lane, s0), roadway, v0)

    return Vehicle(car_state, car_def, 2)
end

function get_car_lanes(roadway::Roadway)
    lanes = Lane[]
    for i=1:length(roadway.segments)
        for lane in roadway.segments[i].lanes
            if !is_crosswalk(lane) && !(lane.tag ∈ [LaneTag(15,1), LaneTag(16, 1), LaneTag(6,1), LaneTag(13, 1), LaneTag(14,1)]) #XXX hack
                push!(lanes, lane)
            end
        end
    end
    return lanes
end

function get_crosswalks(roadway::Roadway)
    lanes = Lane[]
    for i=1:length(roadway.segments)
        for lane in roadway.segments[i].lanes
            if is_crosswalk(lane)
                push!(lanes, lane)
            end
        end
    end
    return lanes
end

is_crosswalk(lane::Lane) = lane.width > 3.0  #XXX hack


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""FlowSchema defines the schema of a group of flows. Note that a flow is made up of a set of inbound API requests with similar attributes and is identified by a pair of strings: the name of the FlowSchema and a \"flow distinguisher\".

    IoK8sApiFlowcontrolV1alpha1FlowSchema(;
        apiVersion=nothing,
        kind=nothing,
        metadata=nothing,
        spec=nothing,
        status=nothing,
    )

    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - metadata::IoK8sApimachineryPkgApisMetaV1ObjectMeta : `metadata` is the standard object's metadata. More info: https://git.k8s.io/community/contributors/devel/api-conventions.md#metadata
    - spec::IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec : `spec` is the specification of the desired behavior of a FlowSchema. More info: https://git.k8s.io/community/contributors/devel/api-conventions.md#spec-and-status
    - status::IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus : `status` is the current status of a FlowSchema. More info: https://git.k8s.io/community/contributors/devel/api-conventions.md#spec-and-status
"""
mutable struct IoK8sApiFlowcontrolV1alpha1FlowSchema <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec } # spec name: spec
    status::Any # spec type: Union{ Nothing, IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus } # spec name: status

    function IoK8sApiFlowcontrolV1alpha1FlowSchema(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing, status=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1FlowSchema, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiFlowcontrolV1alpha1FlowSchema, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiFlowcontrolV1alpha1FlowSchema, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiFlowcontrolV1alpha1FlowSchema, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        validate_property(IoK8sApiFlowcontrolV1alpha1FlowSchema, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1FlowSchema

const _property_map_IoK8sApiFlowcontrolV1alpha1FlowSchema = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"), Symbol("status")=>Symbol("status"))
const _property_types_IoK8sApiFlowcontrolV1alpha1FlowSchema = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec", Symbol("status")=>"IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1FlowSchema }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1FlowSchema))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1FlowSchema }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1FlowSchema[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1FlowSchema }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1FlowSchema[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1FlowSchema)
    true
end

function validate_property(::Type{ IoK8sApiFlowcontrolV1alpha1FlowSchema }, name::Symbol, val)
end


text:
using GLVisualize
include(GLVisualize.dir("examples", "ExampleRunner.jl"))
using ExampleRunner
import ExampleRunner: flatten_paths

files = [
    "introduction/rotate_robj.jl",
    "introduction/screens.jl",
    "plots/3dplots.jl",
    "plots/lines_scatter.jl",
    "plots/hybrid.jl",
    "camera/camera.jl",
    "gui/color_chooser.jl",
    "gui/image_processing.jl",
    "gui/buttons.jl",
    "gui/fractal_lines.jl",
    "gui/mandalas.jl",
    "plots/drawing.jl",
    "interactive/graph_editing.jl",
    "interactive/mario_game.jl",
    "text/text_particle.jl",
]

map!(x-> GLVisualize.dir("examples", x), files, files)
files = union(files, flatten_paths(GLVisualize.dir("examples")))
#push!(files, GLVisualize.dir("test", "summary.jl"))
files = unique(normpath.(files))
# Create an examplerunner, that displays all examples in the example folder, plus
# a runtest specific summary.
config = ExampleRunner.RunnerConfig(
    record = false,
    files = files,
    resolution = (210, 210)
)

ExampleRunner.run(config)


text:
using Pkg
Pkg.activate("./")
Pkg.instantiate()
using FaceRecognition

#%%
#   generate_exampleS("./data/images/", 1000, 50, 1001:1010, "example1")

function generate_example(
      image_dir::AbstractString,
      n::Int,
      d::Int,
      test_range::UnitRange{Int};
      filename::AbstractString,
      save_eigenfaces::Bool=false
)
      training_images = load_images(image_dir, n)
      test_images = load_images(image_dir, test_range)

      model = train_model(training_images, d; normalize=true)

      reconstructed_images = reconstruct_images(test_images, model)
      example = vcat(
            reduce(hcat, test_images),
            reduce(hcat, reconstructed_images),
            get_difference(test_images, reconstructed_images).*3
      )
      save(string("./examples/", filename, ".png"),  example)
      if save_eigenfaces
            eigenfaces = reduce(vcat, (reduce(hcat, get_eigenfaces(model, d))[:, i*128*10+1:(i+1)*128*10] for i in 0:(d ÷ 10)-1))
            save(string("./examples/eigenfaces.png"), eigenfaces)
      end
      return example
end

#generate_example("./data/images/", 1000, 50, 1001:1010, "example2")

#%%


text:
using QPSReader, QuadraticModels, QuadraticModelsCPLEX
using Test

@testset "QuadraticModelsCPLEX.jl" begin
  qps1 = readqps("QAFIRO.SIF") #lower bounds
  stats1 = cplex(QuadraticModel(qps1))
  @test isapprox(stats1.objective, -1.59078179, atol=1e-2)
  @test stats1.status == :acceptable

  qps2 = readqps("HS21.SIF") # low/upp bounds
  stats2 = cplex(QuadraticModel(qps2))
  @test isapprox(stats2.objective, -9.99599999e1, atol=1e-2)
  @test stats2.status == :acceptable

  qps3 = readqps("HS52.SIF") # free bounds
  stats3 = cplex(QuadraticModel(qps3))
  @test isapprox(stats3.objective, 5.32664756, atol=1e-2)
  @test stats3.status == :acceptable
end

@testset "presolve" begin
  qps1 = readqps("QAFIRO.SIF") #lower bounds
  presolve_to_file(QuadraticModel(qps1), path = "QAFIRO_PS.mps")
  qpsp1 = readqps("QAFIRO_PS.mps")
  stats1 = cplex(QuadraticModel(qpsp1))
  @test isapprox(stats1.objective, -1.59078179, atol=1e-2)
  @test stats1.status == :acceptable

  qps2 = readqps("HS21.SIF") #lower bounds
  @test_throws CPXPresolveException presolve_to_file(QuadraticModel(qps2), path = "HS21_PS.mps")

  qps3 = readqps("HS52.SIF") #lower bounds
  presolve_to_file(QuadraticModel(qps3), path = "HS52_PS.mps")
  qpsp3 = readqps("HS52_PS.mps")
  stats3 = cplex(QuadraticModel(qpsp3))
  @test isapprox(stats3.objective, 5.32664756, atol=1e-2)
  @test stats3.status == :acceptable
end


text:
export OfflinePolicy, RLTransition

struct RLTransition
    state
    action
    reward
    terminal
    next_state
end

Base.@kwdef struct OfflinePolicy{L,T} <: AbstractPolicy
    learner::L
    dataset::T
    continuous::Bool
    batch_size::Int
end

(π::OfflinePolicy)(env) = π(env, ActionStyle(env), action_space(env))

function (π::OfflinePolicy)(env, ::MinimalActionSet, ::Base.OneTo)
    if π.continuous
        π.learner(env)
    else
        findmax(π.learner(env))[2]
    end
end
(π::OfflinePolicy)(env, ::FullActionSet, ::Base.OneTo) = findmax(π.learner(env), legal_action_space_mask(env))[2]

function (π::OfflinePolicy)(env, ::MinimalActionSet, A)
    if π.continuous
        π.learner(env)
    else
        A[findmax(π.learner(env))[2]]
    end
end
(π::OfflinePolicy)(env, ::FullActionSet, A) = A[findmax(π.learner(env), legal_action_space_mask(env))[2]]

function RLBase.update!(
    p::OfflinePolicy,
    traj::AbstractTrajectory,
    ::AbstractEnv,
    ::PreExperimentStage,
)
    l = p.learner
    if in(:pretrain_step, fieldnames(typeof(l)))
        println("Pretrain...")
        for _ in 1:l.pretrain_step
            inds, batch = sample(l.rng, p.dataset, p.batch_size)
            update!(l, batch)
        end
    end
end

function RLBase.update!(
    p::OfflinePolicy,
    traj::AbstractTrajectory,
    ::AbstractEnv,
    ::PreActStage,
)
    l = p.learner
    l.update_step += 1

    if in(:target_update_freq, fieldnames(typeof(l))) && l.update_step % l.target_update_freq == 0
        copyto!(l.target_approximator, l.approximator)
    end

    l.update_step % l.update_freq == 0 || return

    inds, batch = sample(l.rng, p.dataset, p.batch_size)

    update!(l, batch)
end

function StatsBase.sample(rng::AbstractRNG, dataset::Vector{T}, batch_size::Int) where {T}
    valid_range = 1:length(dataset)
    inds = rand(rng, valid_range, batch_size)
    batch_data = dataset[inds]
    s_length = size(batch_data[1].state)[1]

    s = Array{Float32}(undef, s_length, batch_size)
    s′ = Array{Float32}(undef, s_length, batch_size)
    a = []
    r = []
    t = []
    for (i, data) in enumerate(batch_data)
        s[:, i] = data.state
        push!(a, data.action)
        s′[:, i] = data.next_state
        push!(r, data.reward)
        push!(t, data.terminal)
    end
    batch = NamedTuple{SARTS}((s, a, r, t, s′))
    inds, batch
end

"""
    calculate_CQL_loss(q_value, action; method)
See paper: [Conservative Q-Learning for Offline Reinforcement Learning](https://arxiv.org/abs/2006.04779)
"""
function calculate_CQL_loss(q_value::Matrix{T}, action::Vector{R}; method = "CQL(H)") where {T, R}
    if method == "CQL(H)"
        cql_loss = mean(log.(sum(exp.(q_value), dims=1)) .- q_value[action])
    else
        @error Wrong method parameter
    end
    return cql_loss
end