text:
function comb_sort(seq::Vector{Int64})::Vector{Int64}
    gap = length(seq)
    swap = true
    while gap > 1 || swap
        gap = max(1, floor(Int, gap / 1.25))
        swap = false
        for i = 0:length(seq)-gap-1
            if seq[i+1] > seq[i+gap+1]
                seq[i+1], seq[i+gap+1] = (seq[i+gap+1], seq[i+1])
                swap = true
            end
        end
    end
    return seq
end

if abspath(PROGRAM_FILE) == @__FILE__
    unsorted = [14, 11, 19, 5, 16, 10, 19, 12, 5, 12]
    expected = [5, 5, 10, 11, 12, 12, 14, 16, 19, 19]
    @assert(comb_sort(unsorted) == expected)
    println("OK")
end


text:
# Autogenerated wrapper script for libmseed_jll for i686-w64-mingw32
export libmseed

JLLWrappers.@generate_wrapper_header("libmseed")
JLLWrappers.@declare_library_product(libmseed, "libmseed.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libmseed,
        "bin\\libmseed.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
module ExampleProblems

export SlidingPuzzle, CityNavigation, TSP, GNP

include("npuzzle.jl")
include("cnp.jl")
include("tsp.jl")
include("grid_nav.jl")

end


text:
# Load Julia packages (libraries) needed  for the snippets in chapter 0

using StatisticalRethinking, CmdStan, LinearAlgebra
#gr(size=(600,600));

# CmdStan uses a tmp directory to store the output of cmdstan

ProjDir = rel_path("..", "scripts", "04")
cd(ProjDir)

# ### snippet 4.7

howell1 = CSV.read(rel_path("..", "data", "Howell1.csv"), delim=';')
df = convert(DataFrame, howell1);

# Use only adults

df2 = filter(row -> row[:age] >= 18, df);
first(df2, 5)

# Use data from m4.1s

# Check if the m4.1s.jls file is present. If not, run the model.

!isfile(joinpath(ProjDir, "m4.1s.jls")) && include(joinpath(ProjDir, "m4.1s.jl"))

chn = deserialize(joinpath(ProjDir, "m4.1s.jls"))

# Describe the draws

MCMCChains.describe(chn)

# Plot the density of posterior draws

density(chn, lab="All heights", xlab="height [cm]", ylab="density")

# Compute cor

mu_sigma = hcat(chn.value[:, 2, 1], chn.value[:,1, 1])
LinearAlgebra.diag(cov(mu_sigma))

# Compute cov

cor(mu_sigma)

# End of `clip_07.0s.jl`


text:
@testset "122.best-time-to-buy-and-sell-stock-ii.jl" begin
    @test max_profit([7, 1, 5, 3, 6, 4]) == 7
    @test max_profit([1, 2, 3, 4, 5]) == 4
    @test max_profit([7, 6, 4, 3, 1]) == 0
end


text:
using TopOpt, Zygote, ChainRulesCore
# using Makie
Nonconvex.@load Ipopt

E = 1.0 # Young’s modulus
v = 0.3 # Poisson’s ratio
f = 1.0 # downward force

problem = PointLoadCantilever(Val{:Linear}, (160, 40), (1.0, 1.0), E, v, f)

V = 0.5
xmin = 0.0001
rmin = 4.0
p = 3.0
convcriteria = KKTCriteria()
solver = FEASolver(Direct, problem; xmin)
nvars = length(solver.vars)
x0 = fill(V, nvars)

penalty = TopOpt.PowerPenalty(p)
solver = FEASolver(Direct, problem; xmin, penalty)
filter = DensityFilter(solver; rmin)
comp = Compliance(problem, solver)
volfrac = Volume(problem, solver)

m = 20
act = relu
nn = Chain(Dense(2, m, act), Dense(m, m, act), Dense(m, m, act), Dense(m, 1, sigmoid))
nn_model = NeuralNetwork(nn, problem)
tf = TrainFunction(nn_model)
p0 = nn_model.init_params
tf(p0)

obj = p -> comp(filter(tf(p)))
constr = p -> volfrac(filter(tf(p))) - V

alg = MMA()
options = MMAOptions()

##
alg = IpoptAlg()
options = IpoptOptions(; max_iter=200)
model1 = Model()
nparams = length(p0)
addvar!(model1, fill(-100.0, nparams), fill(100.0, nparams))
set_objective!(model1, p -> (constr(p) + 0.1)^2)
res1 = optimize(model1, alg, p0; options)
##

struct LoggerFunction{F} <: Function
    f::F
end
function (lf::LoggerFunction)(x)
    @info "Solution extrema: $(extrema(x))"
    val = lf.f(x)
    @info "Objective value: $(val)"
    return val
end
function ChainRulesCore.rrule(lf::LoggerFunction, x)
    @info "Solution extrema: $(extrema(x))"
    val, pb = Zygote.pullback(lf.f, x)
    @info "Objective value: $(val)"
    return val, Δ -> begin
        temp = pb(Δ)
        return (NoTangent(), temp[1])
    end
end

# alg = IpoptAlg()
# nparams = length(p0)
# res2 = res1
# for _ in 1:10
#     global res2
#     # Incremental optimization
#     model2 = Model()
#     addvar!(model2, res2.minimizer .- 1.0, res2.minimizer .+ 1.0)
#     noise = randn(nparams)
#     set_objective!(model2, LoggerFunction(p -> obj(p) + 0.1 * norm(p - res2.minimizer) + dot(noise, p)))
#     # set_objective!(model2, LoggerFunction(obj))
#     add_ineq_constraint!(model2, p -> constr(p))
#     # res2 = optimize(model2, alg, res1.minimizer; options)
#     options = IpoptOptions(max_iter = 10)
#     res2 = optimize(model2, alg, res2.minimizer; options)

#     # Feasibility restoration
#     # model2 = Model()
#     # addvar!(model2, res2.minimizer .- 2.0, res2.minimizer .+ 2.0)
#     # set_objective!(model2, p -> norm(p - res2.minimizer))
#     # add_ineq_constraint!(model2, p -> 100*constr(p))
#     # options = IpoptOptions(max_iter = 5)
#     # res2 = optimize(model2, alg, res2.minimizer; options)
# end

# model2 = Model()
# addvar!(model2, fill(-10.0, nparams), fill(10.0, nparams))
# noise = randn(nparams)
# # set_objective!(model2, obj)
# set_objective!(model2, p -> obj(p) + 0.1 * norm(p))
# add_ineq_constraint!(model2, constr)
# # add_ineq_constraint!(model2, p -> 10*constr(p))
# options = IpoptOptions(max_iter = 300)
# res2 = optimize(model2, alg, res1.minimizer; options)

μ = 1.0
res2 = res1
# Increase the number of iterations to get a good design
for _ in 1:3
    global μ *= 2
    global res2
    model2 = Model()
    addvar!(model2, fill(-100.0, nparams), fill(100.0, nparams))
    noise = randn(nparams)
    # set_objective!(model2, obj)
    set_objective!(model2, p -> μ * obj(p) - log(max(0, -constr(p))))
    # Increase the following number e.g. to 100 for a better design
    options = IpoptOptions(; max_iter=20)
    res2 = optimize(model2, alg, res2.minimizer; options)
    @show extrema(tf(res2.minimizer))
    @show obj(res2.minimizer)
    @show constr(res2.minimizer)
end

Zygote.gradient(obj, p0)
Zygote.gradient(constr, p0)

tf(res2.minimizer)

#=
using Makie, CairoMakie

topology = filter(tf(res2.minimizer))
using TopOpt.TopOptProblems.Visualization: visualize
visualize(problem; topology)
=#


text:
using DiscreteEventsLite

scheduler = Scheduler()
register!(scheduler, s->s, every, 1.0, 2)
@time run!(scheduler, 10^5)


text:
module ADatasets
using DataFrames, MLDataPattern, FileIO, Distances, StatsBase, Printf, Statistics, Random, JLD2, CSV

surveydatasets = ["breast-cancer-wisconsin", "cardiotocography", "ecoli", "magic-telescope", "waveform-1", "statlog-segment", "wall-following-robot", "yeast", "sonar"]
easydatasets = ["abalone", "blood-transfusion", "breast-cancer-wisconsin", "breast-tissue", "cardiotocography", "ecoli", "glass", "haberman", "ionosphere", "iris", "isolet", "letter-recognition", "libras", "magic-telescope", "miniboone", "multiple-features", "musk-2", "page-blocks", "parkinsons", "pendigits", "pima-indians", "sonar", "spect-heart", "statlog-satimage", "statlog-segment", "statlog-shuttle", "statlog-vehicle", "synthetic-control-chart", "wall-following-robot", "waveform-1", "waveform-2", "wine", "yeast"]

include("roc.jl")
include("perfmeasures.jl")
include("densitylevel.jl")
include("makesets.jl")
include("evaluation.jl")

"""
		function append2file(fname::String,dname::String,d::DataFrame)

		add dataframe `d` to the file `fname` under the name `dname`.
		If dataframe exists, it is appended to the end
"""
function append2file(fname::String,dname::String,d::DataFrame)
  !isdir(dirname(fname)) && mkpath(dirname(fname))
  if isfile(fname)
    a = load(fname)
    if haskey(a,dname)
    	a[dname] = vcat(a[dname],d)
    else
    	a[dname] = d;
    end
  	save(fname,a)
  else
  	save(fname,dname,d)
  end
end

function append2file(fname::String,d::DataFrame)
  !isdir(dirname(fname)) && mkpath(dirname(fname))
  d = isfile(fname) ? vcat(load(fname), d) : d
  save(fname,d)
end

end # module


text:
using Test
include("Wasserstein.jl")

#	test Wasserstein
println("\n____________ Test Wasserstein.jl")
p1= [.3, 0.7]
s1= [-.7, 1.1]
p2= [0.2, 0.3, 0.2, 0.1, 0.2]
s2= [-.3, .9, 1.7, 3.2, 1.9]
d= distFunction(s1, s2)

@show A= Wasserstein(p1, p2, distFunction(s1, s2), 2.)

@test A.distance ≈ 0.9679876032264049
@test A.π ≈ [0.2 0.1 0.0 0.0 0.0; 0.0 0.2 0.2 0.1 0.2]

println()
@show A= Sinkhorn(p1, p2, distFunction(s1, s2), 2., 4.3)
@test A.distance ≈ 0.967988 atol= 1e-4
@test A.π ≈ [0.2 0.1 0.0 0.0 0.0; 0.0 0.2 0.2 0.1 0.2] atol = 1e-6

println()
@show A= Sinkhorn(p1, p2, distFunction(s1, s2), 2., 1.)
@test A.distance ≈ 0.9923876770989377 atol = .01
@test A.π ≈ [0.2 0.1 0.0 0.0 0.0; 0.0 0.2 0.2 0.1 0.2] atol = 0.02


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 940b1996-fe9d-11ea-2fa4-8b72bee62b76
md"""
# Grid creation and visualization in PDELib.jl

This notebook shows how to perform grid creation and visualization with the assistance of the packages [ExtendableGrids.jl](https://github.com/j-fu/ExtendableGrids.jl) and [SimplexGridFactory.jl](https://github.com/j-fu/SimplexGridFactory.jl) which are part of the  [PDELib.jl](https://github.com/WIAS-BERLIN/PDElib.jl) meta package.

Visualization in this notebook is done using the [GridVisualize.jl](https://github.com/j-fu/GridVisualize.jl) package.
"""

# ╔═╡ d432ad64-f91f-11ea-2e48-4bc7472ac64c
begin
	using SimplexGridFactory
	using ExtendableGrids
	using Triangulate
	using TetGen
	using GridVisualize
        using PlutoVista
	using PlutoUI
        isdefined(Main,:PlutoRunner) ? 	default_plotter!(PlutoVista) : default_plotter!(nothing)
end

# ╔═╡ 47103e8f-a4ff-46ed-a632-572a2e194a50
md"""
## 1D grids

1D grids are created just from arrays of montonically increasing  coordinates
using the [simplexgrid](https://j-fu.github.io/ExtendableGrids.jl/stable/simplexgrid/#ExtendableGrids.simplexgrid-Tuple{AbstractVector{T}%20where%20T}) method.
"""

# ╔═╡ 93a0c45e-d6a3-415a-a82c-e4f7e2a09d22
X1=range(0,1;length=11)

# ╔═╡ 4622a1fc-fda7-4211-9cc0-4eb1a1584aa6
g1=simplexgrid(X1); @show g1

# ╔═╡ 3212b930-194d-422e-9d06-65885b25cc6d
md"""
We can plot a grid with a method from `GridVisualize.jl`
"""

# ╔═╡ 5520b8c0-0874-4790-a956-224e6c43d9cf
gridplot(g1; resolution=(500,150),legend=:rt)

# ╔═╡ 13aef2a1-5745-4fe5-9659-6b7c0e7267fc
md"""
We see some  additional information:

- `cellregion`: each grid cell (interval, triangle, tetrahedron) as an integer region marker attached
- `bfaceregion`: boundary faces (points, lines, triangles) have an interger boundary region marker attached

We can also have a look into the grid structure:
"""

# ╔═╡ f19d20c8-4d69-442d-99c8-10874fa0a6d3
g1.components

# ╔═╡ 28e2e3b0-c168-481b-b467-29e6a5407431
md"""
Components can be accessed via `[ ]`. In fact the keys in the dictionary of components are [types](https://j-fu.github.io/ExtendableGrids.jl/stable/tdict/).
"""

# ╔═╡ f04017d7-1c55-4118-8467-1e134259e35d
g1[Coordinates]

# ╔═╡ 98016b49-8786-46b9-aca6-01d15c253b3f
g1[CellNodes]

# ╔═╡ 2d392a98-bf32-4145-ae93-a8e218367277
md"""
### Modifying region markers

The `simplexgrid` method provides a default distribution of markers, but we would like to be able to change them. This can be done by putting masks on cells or faces (points in 1D):
"""

# ╔═╡ 88247350-aa6c-4876-82c3-3534036d5702
g2=deepcopy(g1)

# ╔═╡ 42f8d91d-14c7-488f-bdb9-b3d22705186f
cellmask!(g2, [0.0], [0.5], 2);

# ╔═╡ db7e0233-4e08-41b8-9ebe-5570e6d32264
bfacemask!(g2, [0.5],[0.5], 3);

# ╔═╡ 74641f73-2efe-4df8-bebb-ed97e77d869e
gridplot(g2; resolution=(500, 150),legend=:rt)

# ╔═╡ 14fb7977-93cf-4f74-a8ec-b6ee25dbdf86
md"""
### Creating locally refined grids

For this purpose, we just need to create arrays with the corresponding coordinate values. This can be done programmatically.

Two support metods are provided for this purpose.
"""

# ╔═╡ 2d5cb9e1-2d14-415e-b792-c3124901011d
hmin=0.01 ; hmax=0.1

# ╔═╡ b1f903b3-29d7-4909-b7e2-8ef3528c9965
md"""
The `geomspace` method creates an array such that the smallest interval size is `hmin` and the largest interval size is not larger but close to `hmax`, and the interval sizes constitute a geometric sequence.
"""

# ╔═╡ 19125aed-5c46-4968-bcfc-0b469628b68e
X2L=geomspace(0,0.5,hmax,hmin)

# ╔═╡ f9debaf8-efe8-44d1-9671-4aa5bffa9bb8
DX2=X2L[2:end].-X2L[1:end-1]

# ╔═╡ 4f7f9eaf-61c5-4543-ae39-e58ca14fb89a
DX2[1:end-1]./DX2[2:end]

# ╔═╡ f873cb89-da6e-4e11-b3f6-2bf0b766ce5f
X2R=geomspace(0.5,1,hmin,hmax)

# ╔═╡ a2e6c397-4a7e-47eb-ad84-df6e9d3fdd43
md"""
We can glue these arrays together and create a grid from them:
"""

# ╔═╡ 805c2a1e-a6c6-47fc-b719-31c2a671a8d0
X2=glue(X2L,X2R)

# ╔═╡ dd7d8e17-f825-4047-9386-d9e2bfd0a48d
gridplot(simplexgrid(X2); resolution=(500,150),legend=:rt)

# ╔═╡ 3d7db57f-2864-4984-a979-609e1d838a9f
md"""
### Plotting functions

We assume that functions can be represented by their node values an plotted via their piecewise linear interpolants. E.g. they could come from some simulation.
"""

# ╔═╡ 4cf840e6-86c5-4af1-a780-6fc78b60716b
g1d2=simplexgrid(range(-10,10,length=201))

# ╔═╡ 151cc4b8-c5ed-4f5e-8d5f-2f708c9c7fae
fsin=map(sin,g1d2)

# ╔═╡ 605581c4-3c6f-4c31-bd90-7645d3f70315
fcos=map(cos,g1d2)

# ╔═╡ bd40e614-00ef-426d-aa98-eca2ef48320e
fsinh=map(x->sinh(0.2*x), g1d2)

# ╔═╡ 71af99ab-612a-4821-8e9c-efc8766e2e3e
let
	vis=GridVisualizer(;resolution=(600,300),legend=:lt)

	scalarplot!(vis, g1d2, fsinh, label="sinh", markershape=:dtriangle, color=:red,markevery=5,clear=false)

	scalarplot!(vis, g1d2, fcos, label="cos", markershape=:xcross, color=:green, linestyle=:dash, clear=false,markevery=20)

	scalarplot!(vis, g1d2, fsin, label="sin", markershape=:none, color=:blue, linestyle=:dot, clear=false, markevery=20)

	reveal(vis)
end

# ╔═╡ 6dfa1d73-8baa-4589-b2e6-547834c9e444
md"""
## 2D grids
### Tensor product grids

For 2D tensor product grids, we can again use the `simplexgrid` method
and apply the mask methods for modifying cell and boundary region markers.
"""

# ╔═╡ f9599246-8238-432c-a315-300d74abfa2c
begin
	g2d1=simplexgrid(X1,X2)
	cellmask!(g2d1, [0.0,0.0], [0.5, 0.5], 2)
	cellmask!(g2d1, [0.5,0.5], [1.0, 1.0], 3)
	bfacemask!(g2d1, [0.0, 0.0], [0.0, 0.5],5)
end

# ╔═╡ ba5af291-fb16-4f21-8b74-664284bf7bd9
gridplot(g2d1,resolution=(600,400),linewidth=0.5,legend=:lt)

# ╔═╡ dfbaacea-4cb4-4147-9f1b-4424d7a7e89b
md"""
To interact with the plot, you can use the mouse wheel or double toch to zoom,
"shift-mouse-left" to pan, and "alt-mouse-left" or "ctrl-mouse-left" to reset.
"""

# ╔═╡ 7d1698dd-3bb7-4b38-9c6d-a88652749eee
md"""
We can also have a look into the components of a 2D grid:
"""

# ╔═╡ 81249f7c-abdf-43cc-b57f-2915b09da009
g2d1.components

# ╔═╡ 0765641a-8ed9-4579-bd9b-90bb02a55792
md"""
### Unstructured grids

For the triangulation of unstructured grids, we use the mesh generator Triangle via the [Triangulate.jl](https://github.com/JuliaGeometry/Triangulate.jl)  and [SimplexGridFactory.jl](https://github.com/j-fu/SimplexGridFactory.jl) packages.

The later package exports the `SimplexGridBuilder` which shall help to simplify the creation of the input for `Triangulate`.
"""

# ╔═╡ 884a11a2-15cf-40fc-a1ca-66ea23c6094e
builder2=let
	b=SimplexGridBuilder(Generator=Triangulate)
	p1=point!(b,0,0)
	p2=point!(b,1,0)
	p3=point!(b,1,1)

	# Specify outer boundary
	facetregion!(b,1)
	facet!(b,p1,p2)
	facetregion!(b,2)
	facet!(b,p2,p3)
	facetregion!(b,3)
	facet!(b,p3,p1)

	cellregion!(b,1)
	regionpoint!(b,0.75,0.25)

	options!(b,maxvolume=0.01)
	b
end

# ╔═╡ 8fbd238c-723e-4bce-af69-9cabfc03f8d9
md"""
We can plot the current state of the builder:
"""

# ╔═╡ 342788d4-e5d0-4239-be87-7658bb67c999
grid2d2=simplexgrid(builder2;maxvolume=0.001)

# ╔═╡ 8f7d958e-5dc5-4324-af21-ad829d7d77eb
gridplot(grid2d2, resolution=(400,300),linewidth=0.5)

# ╔═╡ fd27b44a-f923-11ea-2afb-d79f7e62e214
md"""
### More complicated grids

More complicated grids include:

- local refinement
- interior boundaries
- different region markers
- holes

The particular way to describe these things is due to Jonathan Shewchuk and his mesh generator [Triangle](https://www.cs.cmu.edu/~quake/triangle.html) via its Julia wrapper package  [Triangulate.jl](https://github.com/JuliaGeometry/Triangulate.jl).
"""

# ╔═╡ 4a289b23-46b9-495d-b19c-42b3da71b242
md"""
#### Local refinement
"""

# ╔═╡ b12838f0-fe9c-11ea-2939-155ed907322d
refinement_center=[0.8,0.2]

# ╔═╡ d5d8a1d6-fe9d-11ea-0fd8-df6e81492cb5
md"""
For local refimenent, we define  a function, which is able to
tell if a triangle is to be refined ("unsuitable") or can be kept as it is.

The function measures the distance between the refinement center and the  triangle barycenter. We require that the area increases with the distance from the refinement center.
"""

# ╔═╡ aae2e82a-fe9c-11ea-0427-593f8d2c7746
function unsuitable(x1,y1,x2,y2,x3,y3,area)
        bary_x=(x1+x2+x3)/3.0
        bary_y=(y1+y2+y3)/3.0
        dx=bary_x-refinement_center[1]
        dy=bary_y-refinement_center[2]
        qdist=dx^2+dy^2
        area>0.1*max(1.0e-2,qdist)
end;

# ╔═╡ 1ae86964-fe9e-11ea-303b-65bb128384a5
md"""
#### Interior boundaries

Interior boundaries are described in a similar as exterior ones - just by facets connecting points.
"""

# ╔═╡ 3b8cd906-bc4e-44af-bcf4-8836d597ed4c
md"""
#### Subregions

Subregions are defined as regions surrounded by interior boundaries.
By placing a "region point" into such a region and specifying a "region number",
we can set the cell region marker for all triangles created in the subregion.
"""

# ╔═╡ 9d240ef9-6639-4bde-a463-ea78480a970d
md"""
#### Holes
Holes are defined in a similar way as subregions, but a "hole point" is places into the place which shall become the hole.
"""

# ╔═╡ 511b26c6-f920-11ea-1228-51c3750f495c
builder3=let
	b=SimplexGridBuilder(Generator=Triangulate;tol=1.0e-10)

	#  Specify points
	p1=point!(b,0,0)
	p2=point!(b,1,0)
	p3=point!(b,1,1)
	p4=point!(b,0,0.7)

	# Specify outer boundary
	facetregion!(b,1)
	facet!(b,p1,p2)
	facetregion!(b,2)
	facet!(b,p2,p3)
	facetregion!(b,3)
	facet!(b,p3,p4)
	facetregion!(b,4)
	facet!(b,p1,p4)

	# Activate unsuitable callback
	options!(b,unsuitable=unsuitable)

	# Specify interior boundary
	facetregion!(b,5)
	facet!(b,p1,p3)

	# Coarse elements in upper left region #1
	cellregion!(b,1)
	maxvolume!(b,0.1)
	regionpoint!(b,0.1,0.5)

	# Fine elements in lower right region #2
	cellregion!(b,2)
	maxvolume!(b,0.01)
	regionpoint!(b,0.9,0.5)

	# Hole
	hp1=point!(b,0.4,0.1)
	hp2=point!(b,0.6,0.1)
	hp3=point!(b,0.5,0.3)
	holepoint!(b,0.5,0.2)
	facetregion!(b,6)
	facet!(b,hp1,hp2)
	facet!(b,hp2,hp3)
	facet!(b,hp3,hp1)

	b
end;

# ╔═╡ d2129483-285b-49a2-a11d-886956146b85
md"""
__Create a simplex grid from the builder__
"""

# ╔═╡ ac93589b-6315-4677-9542-c0a2333f1755
grid2d3=simplexgrid(builder3)

# ╔═╡ 59a6c8b5-25aa-47aa-9489-a803672013df
gridplot(grid2d3,legend=:lt, resolution=(400,400))

# ╔═╡ 4c99c40f-cf93-4cba-bef1-0c4ffcbf6833
md"""
### Plotting of functions

Functions defined on the nodes of a triangular grid can be seen as piecewise linear functions from the P1 finite element space defined by the triangulation.
"""

# ╔═╡ bad736d7-875c-4bc0-9ec4-494a90a508f7
fsin2=map((x,y)-> sin(x)*y, grid2d2)

# ╔═╡ a375c23f-6b8c-4b2c-a8b5-d38e6b5a8f6d
fsin3=map((x,y)-> sin(y)*x, grid2d3)

# ╔═╡ c3ef9067-3cdb-4bfd-9406-6ded64539978
scalarplot(grid2d2, fsin2, label="grid2d2")

# ╔═╡ 4dfb2e0f-3e3a-4053-8a76-765546e96992
	scalarplot(grid2d3, fsin3, label="grid2d3",colormap=:spring,isolines=10)

# ╔═╡ 2682df92-5955-4b17-ae4f-8e99c5b17980
md"""
## 3D Grids

### Tensor product grids

Please note that "masking" is not yet implemented.
Furthermore, PyPlot visualization is slow, with GLMakie it is way faster.
"""

# ╔═╡ 265fe6c7-d1cc-48a6-8295-f8f55acf677c
X3=range(0.,10.1,length=11)

# ╔═╡ b357395f-2a6e-476f-b008-02802c85a541
grid3d1=simplexgrid(X3,X3,X3)

# ╔═╡ af449be7-aab6-4de5-a059-3f8508502676
func3=map((x,y,z)-> sin(x/2)*cos(y/2)*z/10,grid3d1)

# ╔═╡ 38e2b4a8-2480-40e7-bde3-6d1775201aae
p3dg=GridVisualizer(dim=3,resolution=(200,200))

# ╔═╡ f97d085c-e7bf-4561-8183-673912bdeab6
gridplot!(p3dg,grid3d1,zplanes=[zplane],yplanes=[yplane], xplanes=[xplane], resolution=(200,200),show=true)

# ╔═╡ ef1fde48-fe90-4714-ac86-614ae3451aa7
p3ds=GridVisualizer(dim=3,resolution=(400,400))

# ╔═╡ d73d18e7-bcf9-4cc1-9154-b70dc1ff5524

	scalarplot!(p3ds,grid3d1, func3, zplanes=[zplane], yplanes=[yplane],xplanes=[xplane],levels=[flevel],colormap=:spring,resolution=(200,200),show=true,levelalpha=0.5,outlinealpha=0.1)

# ╔═╡ a3844fda-5725-4d95-894b-051a5f6c2faa
md"""
f=$(@bind flevel Slider(range(extrema(func3)...,length=20),default=mean(func3),show_value=true))

x=$(@bind xplane Slider(X3[1]:0.1:X3[end],default=X3[end],show_value=true))

y=$(@bind yplane Slider(X3[1]:0.1:X3[end],default=X3[end],show_value=true))

z=$(@bind zplane Slider(X3[1]:0.1:X3[end],default=X3[end],show_value=true))

"""

# ╔═╡ 04041481-0f03-41e1-a7de-1b3fd033c952
mean(x)=sum(x)/length(x)

# ╔═╡ 6cad87eb-1c59-4000-b688-a6f6d41f9413
md"""
### Unstructured grids

The SimplexGridBuilder API supports creation of three-dimensional grids in way very similar to the 2D case. Just define points with three coordinates and planar (!) facets with at least three points to describe the geometry.

The backend for mesh generation in this case is the [TetGen](http://tetgen.org) mesh generator by Hang Si from WIAS Berlin and its Julia wrapper [TetGen.jl](https://github.com/JuliaGeometry/TetGen.jl).
"""

# ╔═╡ fefc7587-8e25-4080-b934-90c0e1afc56a
builder3d=let

    b=SimplexGridBuilder(Generator=TetGen)

    p1=point!(b,0,0,0)
    p2=point!(b,1,0,0)
    p3=point!(b,1,1,0)
    p4=point!(b,0,1,0)
    p5=point!(b,0,0,1)
    p6=point!(b,1,0,1)
    p7=point!(b,1,1,1)
    p8=point!(b,0,1,1)

    facetregion!(b,1)
    facet!(b,p1 ,p2 ,p3 ,p4)
    facetregion!(b,2)
    facet!(b,p5 ,p6 ,p7 ,p8)
    facetregion!(b,3)
    facet!(b,p1 ,p2 ,p6 ,p5)
    facetregion!(b,4)
    facet!(b,p2 ,p3 ,p7 ,p6)
    facetregion!(b,5)
    facet!(b,p3 ,p4 ,p8 ,p7)
    facetregion!(b,6)
    facet!(b,p4 ,p1 ,p5 ,p8)

	hp1=point!(b,0.4,0.4,0.4)
    hp2=point!(b,0.6,0.4,0.4)
    hp3=point!(b,0.6,0.6,0.4)
    hp4=point!(b,0.4,0.6,0.4)
    hp5=point!(b,0.4,0.4,0.6)
    hp6=point!(b,0.6,0.4,0.6)
    hp7=point!(b,0.6,0.6,0.6)
    hp8=point!(b,0.4,0.6,0.6)

    facetregion!(b,7)
    facet!(b,hp1 ,hp2 ,hp3 ,hp4)
    facet!(b,hp5 ,hp6 ,hp7 ,hp8)
    facet!(b,hp1 ,hp2 ,hp6 ,hp5)
    facet!(b,hp2 ,hp3 ,hp7 ,hp6)
    facet!(b,hp3 ,hp4 ,hp8 ,hp7)
    facet!(b,hp4 ,hp1 ,hp5 ,hp8)
	holepoint!(b, 0.5,0.5,0.5)

	b

end;

# ╔═╡ 065735f7-c799-4284-bd59-fe6383bb987c
grid3d2=simplexgrid(builder3d,maxvolume=0.0001)

# ╔═╡ 329992a0-e352-468b-af8b-0b190315fc61
gridplot(grid3d2,zplane=0.1,azim=20,elev=20,linewidth=0.5,outlinealpha=0.3)

# ╔═╡ a7965a6e-2e83-47eb-aee2-d366246a8637
html"""<hr>"""

# ╔═╡ 7ad541b1-f40f-4cdd-b7b5-b792a8e63d71
TableOfContents(depth=4)

# ╔═╡ 071b8834-d3d1-4d08-979f-56b05bc1e0d3
md"""
    begin
       using Pkg
       Pkg.activate(mktempdir())
       Pkg.add(["PlutoUI","Revise","Triangulate","TetGen"])
       using Revise
	   Pkg.develop(["ExtendableGrids","SimplexGridFactory",
           "GridVisualize","PlutoVista"])
    end
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PlutoVista = "646e1f28-b900-46d7-9d87-d554eb38a413"
SimplexGridFactory = "57bfcd06-606e-45d6-baf4-4ba06da0efd5"
TetGen = "c5d3f3f7-f850-59f6-8a2e-ffc6dc1317ea"
Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"

[compat]
ExtendableGrids = "~0.8.7"
GridVisualize = "~0.3.9"
PlutoUI = "~0.7.16"
PlutoVista = "~0.8.6"
SimplexGridFactory = "~0.5.9"
TetGen = "~1.3.0"
Triangulate = "~2.1.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[AbstractTrees]]
git-tree-sha1 = "03e0550477d86222521d254b741d470ba17ea0b5"
uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
version = "0.3.4"

[[Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "a851fec56cb73cfdf43762999ec72eff5b86882a"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.15.0"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3f3a2501fa7236e9b911e0f7a588c657e822bb6d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.2.3+0"

[[ElasticArrays]]
deps = ["Adapt"]
git-tree-sha1 = "a0fcc1bb3c9ceaf07e1d0529c9806ce94be6adf9"
uuid = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4"
version = "1.2.9"

[[ExtendableGrids]]
deps = ["AbstractTrees", "Dates", "DocStringExtensions", "ElasticArrays", "InteractiveUtils", "LinearAlgebra", "Printf", "Random", "SparseArrays", "Test"]
git-tree-sha1 = "1e8e50f054057f23e908fbd6935766dca6293cc2"
uuid = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
version = "0.8.7"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "58bcdf5ebc057b085e58d95c138725628dd7453c"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.4.1"

[[GridVisualize]]
deps = ["ColorSchemes", "Colors", "DocStringExtensions", "ElasticArrays", "ExtendableGrids", "GeometryBasics", "LinearAlgebra", "Observables", "OrderedCollections", "PkgVersion", "Printf", "Requires", "StaticArrays"]
git-tree-sha1 = "925ba2f11df005d894b113292d32fca9afe3f8c8"
uuid = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
version = "0.3.9"

[[Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[HypertextLiteral]]
git-tree-sha1 = "5efcf53d798efede8fee5b2c8b09284be359bf24"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.2"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Observables]]
git-tree-sha1 = "fe29afdef3d0c4a8286128d4e45cc50621b1e43d"
uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
version = "0.4.0"

[[OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "f19e978f81eca5fd7620650d7dbea58f825802ee"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PkgVersion]]
deps = ["Pkg"]
git-tree-sha1 = "a7a7e1a88853564e551e4eba8650f8c38df79b37"
uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688"
version = "0.1.1"

[[PlutoUI]]
deps = ["Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "4c8a7d080daca18545c56f1cac28710c362478f3"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.16"

[[PlutoVista]]
deps = ["ColorSchemes", "Colors", "DocStringExtensions", "GridVisualize", "UUIDs"]
git-tree-sha1 = "34fc7e41e6eefa58fef0786ab62a20262df88764"
uuid = "646e1f28-b900-46d7-9d87-d554eb38a413"
version = "0.8.6"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SimplexGridFactory]]
deps = ["DocStringExtensions", "ElasticArrays", "ExtendableGrids", "GridVisualize", "LinearAlgebra", "Printf", "Test"]
git-tree-sha1 = "af52ec74a4b6cfcc5b6d60d259099fa0596de2c1"
uuid = "57bfcd06-606e-45d6-baf4-4ba06da0efd5"
version = "0.5.9"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3c76dde64d03699e074ac02eb2e8ba8254d428da"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.13"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StructArrays]]
deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
git-tree-sha1 = "2ce41e0d042c60ecd131e9fb7154a3bfadbf50d3"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.3"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TetGen]]
deps = ["DocStringExtensions", "GeometryBasics", "LinearAlgebra", "Printf", "StaticArrays", "TetGen_jll"]
git-tree-sha1 = "2f1d87ccacd2a7faf9e0bade918946ec4d90bfdf"
uuid = "c5d3f3f7-f850-59f6-8a2e-ffc6dc1317ea"
version = "1.3.0"

[[TetGen_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "bc67a0d0b799fe248b1f199a5c893ccf316f0e60"
uuid = "b47fdcd6-d2c1-58e9-bbba-c1cee8d8c179"
version = "1.5.2+0"

[[Triangle_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "bfdd9ef1004eb9d407af935a6f36a4e0af711369"
uuid = "5639c1d2-226c-5e70-8d55-b3095415a16a"
version = "1.6.1+0"

[[Triangulate]]
deps = ["DocStringExtensions", "Libdl", "Printf", "Test", "Triangle_jll"]
git-tree-sha1 = "2b4f716b192c0c615d96d541ee029e85666388cb"
uuid = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
version = "2.1.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─940b1996-fe9d-11ea-2fa4-8b72bee62b76
# ╠═d432ad64-f91f-11ea-2e48-4bc7472ac64c
# ╟─47103e8f-a4ff-46ed-a632-572a2e194a50
# ╠═93a0c45e-d6a3-415a-a82c-e4f7e2a09d22
# ╠═4622a1fc-fda7-4211-9cc0-4eb1a1584aa6
# ╟─3212b930-194d-422e-9d06-65885b25cc6d
# ╠═5520b8c0-0874-4790-a956-224e6c43d9cf
# ╟─13aef2a1-5745-4fe5-9659-6b7c0e7267fc
# ╠═f19d20c8-4d69-442d-99c8-10874fa0a6d3
# ╟─28e2e3b0-c168-481b-b467-29e6a5407431
# ╠═f04017d7-1c55-4118-8467-1e134259e35d
# ╠═98016b49-8786-46b9-aca6-01d15c253b3f
# ╟─2d392a98-bf32-4145-ae93-a8e218367277
# ╠═88247350-aa6c-4876-82c3-3534036d5702
# ╠═42f8d91d-14c7-488f-bdb9-b3d22705186f
# ╠═db7e0233-4e08-41b8-9ebe-5570e6d32264
# ╠═74641f73-2efe-4df8-bebb-ed97e77d869e
# ╟─14fb7977-93cf-4f74-a8ec-b6ee25dbdf86
# ╠═2d5cb9e1-2d14-415e-b792-c3124901011d
# ╟─b1f903b3-29d7-4909-b7e2-8ef3528c9965
# ╠═19125aed-5c46-4968-bcfc-0b469628b68e
# ╠═f9debaf8-efe8-44d1-9671-4aa5bffa9bb8
# ╠═4f7f9eaf-61c5-4543-ae39-e58ca14fb89a
# ╠═f873cb89-da6e-4e11-b3f6-2bf0b766ce5f
# ╟─a2e6c397-4a7e-47eb-ad84-df6e9d3fdd43
# ╠═805c2a1e-a6c6-47fc-b719-31c2a671a8d0
# ╠═dd7d8e17-f825-4047-9386-d9e2bfd0a48d
# ╟─3d7db57f-2864-4984-a979-609e1d838a9f
# ╠═4cf840e6-86c5-4af1-a780-6fc78b60716b
# ╠═151cc4b8-c5ed-4f5e-8d5f-2f708c9c7fae
# ╠═605581c4-3c6f-4c31-bd90-7645d3f70315
# ╠═bd40e614-00ef-426d-aa98-eca2ef48320e
# ╠═71af99ab-612a-4821-8e9c-efc8766e2e3e
# ╟─6dfa1d73-8baa-4589-b2e6-547834c9e444
# ╠═f9599246-8238-432c-a315-300d74abfa2c
# ╠═ba5af291-fb16-4f21-8b74-664284bf7bd9
# ╟─dfbaacea-4cb4-4147-9f1b-4424d7a7e89b
# ╟─7d1698dd-3bb7-4b38-9c6d-a88652749eee
# ╠═81249f7c-abdf-43cc-b57f-2915b09da009
# ╟─0765641a-8ed9-4579-bd9b-90bb02a55792
# ╠═884a11a2-15cf-40fc-a1ca-66ea23c6094e
# ╟─8fbd238c-723e-4bce-af69-9cabfc03f8d9
# ╠═342788d4-e5d0-4239-be87-7658bb67c999
# ╠═8f7d958e-5dc5-4324-af21-ad829d7d77eb
# ╟─fd27b44a-f923-11ea-2afb-d79f7e62e214
# ╟─4a289b23-46b9-495d-b19c-42b3da71b242
# ╠═b12838f0-fe9c-11ea-2939-155ed907322d
# ╟─d5d8a1d6-fe9d-11ea-0fd8-df6e81492cb5
# ╠═aae2e82a-fe9c-11ea-0427-593f8d2c7746
# ╟─1ae86964-fe9e-11ea-303b-65bb128384a5
# ╟─3b8cd906-bc4e-44af-bcf4-8836d597ed4c
# ╟─9d240ef9-6639-4bde-a463-ea78480a970d
# ╠═511b26c6-f920-11ea-1228-51c3750f495c
# ╟─d2129483-285b-49a2-a11d-886956146b85
# ╠═ac93589b-6315-4677-9542-c0a2333f1755
# ╠═59a6c8b5-25aa-47aa-9489-a803672013df
# ╟─4c99c40f-cf93-4cba-bef1-0c4ffcbf6833
# ╠═bad736d7-875c-4bc0-9ec4-494a90a508f7
# ╠═a375c23f-6b8c-4b2c-a8b5-d38e6b5a8f6d
# ╠═c3ef9067-3cdb-4bfd-9406-6ded64539978
# ╠═4dfb2e0f-3e3a-4053-8a76-765546e96992
# ╟─2682df92-5955-4b17-ae4f-8e99c5b17980
# ╠═265fe6c7-d1cc-48a6-8295-f8f55acf677c
# ╠═b357395f-2a6e-476f-b008-02802c85a541
# ╠═af449be7-aab6-4de5-a059-3f8508502676
# ╠═38e2b4a8-2480-40e7-bde3-6d1775201aae
# ╠═f97d085c-e7bf-4561-8183-673912bdeab6
# ╠═ef1fde48-fe90-4714-ac86-614ae3451aa7
# ╠═d73d18e7-bcf9-4cc1-9154-b70dc1ff5524
# ╠═a3844fda-5725-4d95-894b-051a5f6c2faa
# ╠═04041481-0f03-41e1-a7de-1b3fd033c952
# ╟─6cad87eb-1c59-4000-b688-a6f6d41f9413
# ╠═fefc7587-8e25-4080-b934-90c0e1afc56a
# ╠═065735f7-c799-4284-bd59-fe6383bb987c
# ╠═329992a0-e352-468b-af8b-0b190315fc61
# ╟─a7965a6e-2e83-47eb-aee2-d366246a8637
# ╠═7ad541b1-f40f-4cdd-b7b5-b792a8e63d71
# ╠═071b8834-d3d1-4d08-979f-56b05bc1e0d3
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
"""
** Task **

- `n_trials`: number of trials
- `trial`: current trial
- `width`: screen width
- `height`: screen height
- `scheduler`: event scheduler
- `screen`: visual objects on screen
- `canvas`: GTK canvas
- `window`: GTK window
- `visible`: GUI visible
- `speed`: real time speed
- `study_words`: words presented in study phase
- `test_words`: words presented in study phase

Function Signature

````julia
Task(;n_trials=10, trial=1, lb=2.0, ub=10.0, width=600.0, height=600.0, scheduler=nothing,
    screen=Vector{VisualObject}(), window=nothing, canvas=nothing, visible=false, speed=1.0)
````
"""
@concrete mutable struct Task <: AbstractTask
    study_trial::Int
    test_trial::Int
    n_blocks::Int
    block::Int
    width::Float64
    hight::Float64
    scheduler
    screen::Vector{VisualObject}
    canvas
    window
    visible::Bool
    realtime::Bool
    speed::Float64
    study_words
    test_words
    data
    test_phase
    test_word::NamedTuple
end

function Task(;
    study_trial = 1,
    test_trial = 1,
    n_blocks = 1,
    block = 1,
    width = 600.0,
    height = 600.0,
    scheduler = nothing,
    screen = Vector{VisualObject}(),
    window = nothing,
    canvas = nothing,
    visible = false,
    realtime = false,
    speed = 1.0,
    study_words = String[],
    test_words = String[],
    data = nothing,
    test_phase = false,
    test_word = (a=1,)
    )
    visible ? ((canvas,window) = setup_window(width)) : nothing
    visible ? Gtk.showall(window) : nothing
    return Task(study_trial, test_trial, n_blocks, block, width, height, scheduler, screen, canvas, window, visible,
        realtime, speed, study_words, test_words, data, test_phase, test_word)
end

function draw_object!(task, word)
    c = task.canvas
    w = task.width
	x = w/2
	y = w/2
    @guarded draw(c) do widget
        ctx = getgc(c)
        select_font_face(ctx, "Arial", Cairo.FONT_SLANT_NORMAL,
             Cairo.FONT_WEIGHT_BOLD);
        set_font_size(ctx, 36)
        set_source_rgb(ctx, 0, 0, 0)
        extents = text_extents(ctx, word)
        x′ = x - (extents[3]/2 + extents[1])
        y′ = y - (extents[4]/2 + extents[2])
        move_to(ctx, x′, y′)
        show_text(ctx, word)
    end
    Gtk.showall(c)
    return nothing
end

function start!(task::Task, actr)
    sort!(task.study_words)
    sort!(task.test_words)
    study_trial!(task, actr)
end

function present_stimulus(task, actr, word)
    vo = VisualObject(;text=word)
    add_to_visicon!(actr, vo; stuff=true)
    push!(task.screen, vo)
    task.visible ? draw_object!(task, word) : nothing
end

function study_trial!(task, actr)
    @unpack study_words,study_trial = task
    isi = .5
    description = "present stimulus"
    word = study_words[study_trial].word
    register!(task, present_stimulus, after, isi, task, actr, word;
        description)
    task.visible ? register!(task, clear!, after, 1.5, task) : nothing
    register!(task, update_task!, after, 1.5, task, actr)
end

function start_test(task, actr)
    isi = .5
    description = "start test"
    word = "start test"
    register!(task, present_stimulus, after, isi, task, actr, word;
        description)
    register!(task, empty!, after, 1.0, task.screen)
    task.visible ? register!(task, clear!, after, 1.0, task) : nothing
    register!(task, update_task!, after, 1.0, task, actr)
end

function test_trial!(task, actr)
    @unpack test_words,test_trial = task
    isi = .5
    description = "present stimulus"
    task.test_word = test_words[test_trial]
    word = task.test_word.word
    register!(task, present_stimulus, after, isi, task, actr, word;
        description)
end

function update_block!(task)
    task.study_trial = 0
    task.block += 1
    shuffle!(task.study_words)
end

function update_task!(task, actr)
    # run through training trials
    if task.study_trial < length(task.study_words)
        task.study_trial += 1
        study_trial!(task, actr)
        return nothing
    # repeat training block
    elseif task.block < task.n_blocks
        update_block!(task)
        update_task!(task, actr)
        return nothing
    end
    # signal switch to test trials
    if !task.test_phase
        start_test(task, actr)
        task.test_phase = true
        return nothing
    end
    # run through test trails
    if task.test_trial < length(task.test_words)
        task.test_trial += 1
        test_trial!(task, actr)
        return nothing
    else
        stop!(task.scheduler)
    end
end

function press_key!(task::Task, actr, key)
    @unpack data,test_word = task
    push!(data, [test_word.word test_word.type key])
    empty!(task.screen)
    task.visible ? clear!(task) : nothing
    update_task!(task, actr)
end

function populate_lists()
    study_words = [
        (word="chair",type="target"),
        (word="duck",type="target"),
        (word="book",type="target"),
        (word="pan",type="target"),
        (word="phone",type="target"),
        (word="shirt",type="target"),
        (word="tree",type="target")
    ]
    test_words = [
        (word="chair",type="target"),
        (word="duck",type="target"),
        (word="book",type="target"),
        (word="pan",type="target"),
        (word="phone",type="target"),
        (word="shirt",type="target"),
        (word="tree",type="target"),

        (word="rug",type="foil"),
        (word="car",type="foil"),
        (word="pool",type="foil"),
        (word="lump",type="foil"),
        (word="hair",type="foil"),
        (word="face",type="foil"),
        (word="harp",type="foil"),
        (word="ape",type="foil")
    ]
    stimuli = (study_words=study_words, test_words=test_words)
    return stimuli
end


text:
#-----------------------------------------------------------------------
# MetadataTools
# https://github.com/IainNZ/MetadataTools.jl
#-----------------------------------------------------------------------
# Copyright (c) 2015: Iain Dunning
# Licensed under the MIT License
#-----------------------------------------------------------------------
# graph_stats.jl
# Calculate some rankings for package connectedness
#----------------------------------------------------------------------

using MetadataTools

# Load in all the package data from METADATA
pkgmeta = get_all_pkg()

# Build a PkgGraph with edges A->B if A directly requires B
pg = make_dep_graph(pkgmeta)
total_pkgs = size(pg)

# Count size of every connected subgraph obtained by traversing the
# dependency graph starting from each package
pkg_req_counts = Any[]
for pkg in keys(pkgmeta)
    # Don't warn about cycles, we just want total counts anyway
    g = get_pkg_dep_graph(pkg, pg)
    req_count = size(g)
    push!(pkg_req_counts, (req_count,pkg))
end
sort!(pkg_req_counts, rev=true)  # Sort descending
print_with_color(:blue, "Top packages by number of packages depended\n")
print_with_color(:blue, "on, directly or indirectly:\n")
for (count, pkg) in pkg_req_counts[1:10]
    print(rpad(pkg,30," "))
    @printf("%4d / %4d\n", count-1, total_pkgs)
end
println()

# Do the same, but only counting immediate dependencies
pkg_req_counts = Any[]
for pkg in keys(pkgmeta)
    # Don't warn about cycles, we just want total counts anyway
    g = get_pkg_dep_graph(pkg, pg, depth=2)
    req_count = size(g)
    push!(pkg_req_counts, (req_count,pkg))
end
sort!(pkg_req_counts, rev=true)  # Sort descending
print_with_color(:blue, "Top packages by number of packages depended\n")
print_with_color(:blue, "on, directly:\n")
for (count, pkg) in pkg_req_counts[1:10]
    print(rpad(pkg,30," "))
    @printf("%4d / %4d\n", count-1, total_pkgs)
end
println()

# Build a graph with edges A->B if A is directly required by B
pgrev = make_dep_graph(pkgmeta, reverse=true)

# Count size of every subgraphs like above
pkg_req_counts = Any[]
for pkg in keys(pkgmeta)
    g = get_pkg_dep_graph(pkg, pgrev)
    req_count = size(g)
    push!(pkg_req_counts, (req_count,pkg))
end
sort!(pkg_req_counts, rev=true)  # Sort descending
print_with_color(:blue, "Top packages by number of packages that depend\n")
print_with_color(:blue, "on them, directly or indirectly:\n")
for (count, pkg) in pkg_req_counts[1:10]
    print(rpad(pkg,30," "))
    @printf("%4d / %4d\n", count-1, total_pkgs)
end
println()

# Same, but only direct dependencies
pkg_req_counts = Any[]
for pkg in keys(pkgmeta)
    g = get_pkg_dep_graph(pkg, pgrev, depth=2)
    req_count = size(g)
    push!(pkg_req_counts, (req_count,pkg))
end
sort!(pkg_req_counts, rev=true)  # Sort descending
print_with_color(:blue, "Top packages by number of packages that depend\n")
print_with_color(:blue, "on them, directly:\n")
for (count, pkg) in pkg_req_counts[1:10]
    print(rpad(pkg,30," "))
    @printf("%4d / %4d\n", count-1, total_pkgs)
end


text:
"""
Creates console and file loggers.

**Note:** Log messages may not be written to the file until flush() or close() is called on
the returned logger.

# Arguments
- `console_level = Logging.Error`: level for console messages
- `file_level = Logging.Info`: level for file messages
- `filename::Union{Nothing, AbstractString} = "power-systems.log"`: log file; pass nothing
  to disable file logging

# Example
```julia
logger = configure_logging(console_level = Logging.Info)
@info "log message"
close(logger)
```
"""
function configure_logging(;
    console_level = Logging.Error,
    file_level = Logging.Info,
    filename::Union{Nothing, AbstractString} = "power-systems.log",
)
    return IS.configure_logging(
        console = true,
        console_stream = stderr,
        console_level = console_level,
        file = filename !== nothing,
        filename = filename,
        file_level = file_level,
        file_mode = "w+",
        tracker = nothing,
        set_global = true,
    )
end


text:
my_f(x,y) = 2x + y

greetings() = "hello"


text:
module StructJuMP

import MPI
import JuMP # To reexport, should be using (not import)
import MathProgBase
import MathProgBase.MathProgSolverInterface
import ReverseDiffSparse
import StructJuMPSolverInterface

export StructuredModel, getStructure, getparent, getchildren, getProcIdxSet,
       num_scenarios, @second_stage, getprobability, getMyRank

# ---------------
# StructureData
# ---------------
type MPIWrapper
    comm::MPI.Comm
    init::Function

    function MPIWrapper()
        instance = new(MPI.Comm(-1))
        finalizer(instance, freeMPIWrapper)

        instance.init = function(ucomm::MPI.Comm)
            if isdefined(:MPI) && MPI.Initialized() && ucomm.val == -1
                instance.comm = MPI.COMM_WORLD
            elseif isdefined(:MPI) && !MPI.Initialized()
                MPI.Init()
                instance.comm = MPI.COMM_WORLD
            elseif isdefined(:MPI) && MPI.Initialized() && ucomm.val != -1
                instance.comm = ucomm
            elseif isdefined(:MPI) && MPI.Finalized()
                error("MPI is already finalized!")
            else
                #doing nothing
            end
        end

        return instance
    end

end
function freeMPIWrapper(instance::MPIWrapper)
    if isdefined(:MPI) && MPI.Initialized() && !MPI.Finalized()
        MPI.Finalize()
    end
end

const mpiWrapper = MPIWrapper();

type StructureData
    probability::Vector{Float64}
    children::Dict{Int,JuMP.Model}
    parent
    num_scen::Int
    othermap::Dict{JuMP.Variable,JuMP.Variable}
    mpiWrapper::MPIWrapper
end

default_probability(m::JuMP.Model) = 1 / num_scenarios(m)
default_probability(::Void) = 1.0

# ---------------
# StructuredModel
# ---------------

# Constructor with the number of scenarios
function StructuredModel(;solver=JuMP.UnsetSolver(), parent=nothing, same_children_as=nothing, id=0, comm=MPI.Comm(-1), num_scenarios::Int=0, prob::Float64=default_probability(parent))
    m = JuMP.Model(solver=solver)
    if parent === nothing
        id = 0
        mpiWrapper.init(comm)
        JuMP.setsolvehook(m,StructJuMPSolverInterface.sj_solve)
    else
        @assert id != 0
        stoch = getStructure(parent)
        stoch.children[id] = m
        push!(stoch.probability, prob)
    end

    if same_children_as !== nothing
      if !isa(same_children_as, JuMP.Model) || !haskey(same_children_as.ext, :Stochastic)
        error("The JuMP model given for the argument `same_children_as' is not valid. Please create it using the `StructuredModel' function.")
      end
      probability = same_children_as.ext[:Stochastic].probability
      children = same_children_as.ext[:Stochastic].children
    else
      probability = Float64[]
      children = Dict{Int, JuMP.Model}()
    end
    m.ext[:Stochastic] = StructureData(probability, children, parent, num_scenarios, Dict{JuMP.Variable,JuMP.Variable}(), mpiWrapper)
    m
end

# -------------
# Get functions
# -------------

getStructure(m::JuMP.Model)  = m.ext[:Stochastic]::StructureData
getparent(m::JuMP.Model)      = getStructure(m).parent
getchildren(m::JuMP.Model)    = getStructure(m).children::Dict{Int,JuMP.Model}
getprobability(m::JuMP.Model) = getStructure(m).probability::Vector{Float64}
num_scenarios(m::JuMP.Model)  = getStructure(m).num_scen::Int

function getMyRank()
    myrank = 0;
    mysize = 1;
    if isdefined(:MPI) && MPI.Initialized() && !MPI.Finalized()
        comm = MPI.COMM_WORLD
        mysize = MPI.Comm_size(comm)
        myrank = MPI.Comm_rank(comm)
    end
    return myrank,mysize
end

function getProcIdxSet(numScens::Integer)
    mysize = 1;
    myrank = 0;
    if isdefined(:MPI) == true && MPI.Initialized() == true
        comm = MPI.COMM_WORLD
        mysize = MPI.Comm_size(comm)
        myrank = MPI.Comm_rank(comm)
    end
    # Why don't we just take a round-and-robin?
    proc_idx_set = Int[];
    for s = myrank:mysize:(numScens-1)
        push!(proc_idx_set, s+1);
    end
    return proc_idx_set;
end

function getProcIdxSet(m::JuMP.Model)
    numScens = num_scenarios(m)
    return getProcIdxSet(numScens)
end

macro second_stage(m,ind,code)
    return quote
        proc_idx_set = getProcIdxSet($(esc(m)))
        for $(esc(ind)) in proc_idx_set
            $(esc(code))
        end
    end
end

include("nlp.jl")

end


text:
immutable FNN <: Model; net; FNN(a...;o...)=new(Net(a...;o...)); end
params(m::FNN)=params(m.net)

function train(m::FNN, data, loss; gclip=0, gcheck=false, maxnorm=nothing, losscnt=nothing, o...)
    for item in data
        (x,ygold) = item2xy(item)
        reset!(m.net; o...)
        ypred = forw(m.net, x...; trn=true, o...)
        losscnt != nothing && (losscnt[1] += loss(ypred, ygold); losscnt[2] += 1)
        back(m.net, ygold, loss; o...)
        gcheck && break
        g = (gclip > 0 || maxnorm!=nothing ? gnorm(m) : 0)
        update!(m; gclip=(g > gclip > 0 ? gclip/g : 0), o...)
        if maxnorm != nothing
            w=wnorm(m)
            w > maxnorm[1] && (maxnorm[1]=w)
            g > maxnorm[2] && (maxnorm[2]=g)
        end
    end
end

function test(m::FNN, data, loss; gcheck=false, o...)
    sumloss = numloss = 0
    for item in data
        (x,ygold) = item2xy(item)
        ypred = forw(m.net, x...; trn=false, o...)
        sumloss += loss(ypred, ygold); numloss += 1
        gcheck && break
    end
    sumloss / numloss
end

function predict(m::FNN, data; o...)
    y = Any[]
    for x in data
        ypred = forw(m.net, x...; trn=false, o...)
        ycopy = isa(ypred, Array) ? copy(ypred) : convert(Array, ypred)
        push!(y, ycopy)
    end
    return y
end

# function gradcheck(m::FNN, data, loss; delta=1e-4, rtol=eps(Float64)^(1/5), atol=eps(Float64)^(1/5), gcheck=10, o...)
#     x = ygold = loss0 = nothing
#     for item in data
#         (x,ygold) = item2xy(item)
#         ypred = forw(m.net, x...; trn=true, o...)
#         loss0 = loss(ypred, ygold)
#         back(m.net, ygold, loss; o...)
#         break
#     end
#     pp = params(m)
#     for n=1:length(pp)
#         p = pp[n]
#         psave = copy(p.out)
#         pdiff = convert(Array, p.dif)
#         wlen = length(p.out)
#         irange = (wlen <= gcheck ? (1:wlen) : rand(gradcheck_rng, 1:wlen, gcheck))
#         for i in irange
#             wi0 = p.out[i]
#             wi1 = (wi0 >= 0 ? wi0 + delta : wi0 - delta)
#             p.out[i] = wi1
#             ypred = forw(m.net, x...; trn=false, o...)
#             loss1 = loss(ypred, ygold)
#             p.out[i] = wi0
#             dwi = (loss1 - loss0) / (wi1 - wi0)
#             if !isapprox(pdiff[i], dwi; rtol=rtol, atol=atol)
#                 println(tuple(:gc, n, i, pdiff[i], dwi))
#             end
#         end
#         @assert isequal(p.out, psave)
#     end
# end


text:
const SGKmer{K} = BioSequences.DNAMer{K}

Base.write(io::IO, k::BioSequences.DNAMer) = write(io, reinterpret(UInt64, k))
Base.read(io::IO, ::Type{T}) where {T <: BioSequences.DNAMer} = convert(T, Base.read(io, UInt64))

# Conversion
# ----------

kmer_index( kmer::BioSequences.DNAMer{K} ) where {K} = Int(reinterpret(UInt64, kmer)) + 1

kmer_index( seq::BioSequence ) = Int(kmer_index_trailing( seq )) + 1

# calculate kmer index directly
@inline function kmer_index_trailing( ::Type{T}, seq ) where T <: Unsigned
   x     = T(0)
   for nt in seq
      ntint = convert(UInt8, trailing_zeros(nt))
      if ntint > 0x03 || isambiguous(nt)
         return zero(T)
      else
         x = x << 2 | ntint
      end
   end
   x
end

kmer_index_trailing( seq ) = kmer_index_trailing(UInt64, seq )

kmer(seq::SGSequence) = sgkmer(seq)
sgkmer(seq::String) = SGKmer(seq)
sgkmer(seq::SGSequence) = SGKmer(seq)


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

mpi_init()

function kappa(x, y)
    return 1+x^2+y^2
end

function rhs(x, y)
    -2*x*(1 - x)*(x^2 + y^2 + 1) + 2*x*(-x*y*(1 - y) + y*(1 - x)*(1 - y)) -
             2*y*(1 - y)*(x^2 + y^2 + 1) + 2*y*(-x*y*(1 - x) + x*(1 - x)*(1 - y))
end

function ufunc(x, y)
    x * (1-x) * y * (1-y)
end

mc = MPIConfig(50)
global_to_local, local_to_global = dofmap(mc)
X, Y = get_xy(mc)
f_local = rhs.(X, Y)
using Random; Random.seed!(233)
θ = Variable(fc_init([2,20,20,20,1]))
θ_shared = mpi_bcast(θ)
κ_local = (fc([X'[:] Y'[:]], [20,20,20,1], θ_shared) + 5.0)|>squeeze
κ_local = reshape(κ_local, (mc.n, mc.n))
u_local = poisson_solver(κ_local, f_local, mc)
u = mpi_gather(u_local)[global_to_local]

loss = sum(u^2)
g = gradients(loss, θ)

sess = Session(); init(sess)

function calculate_loss(x)
    run(sess, loss, θ=>x)
end

function calculate_gradients(x)
    run(sess, g, θ=>x)
end

f(x) = (calculate_loss(x), calculate_gradients(x))

test_gradients(f, run(sess, θ), mpi=true)

if mpi_size()>1
    mpi_finalize()
end


text:
# TODO

include("SiemensRaw.jl")

struct SiemensFile <: MRIFile
  filename::String
end

function MRIBase.RawAcquisitionData(f::SiemensFile)
  error("Not yet implemented!")
end

export readRawSiemens
function readRawSiemens(filename)

  # See https://github.com/ismrmrd/siemens_to_ismrmrd/blob/master/main.cpp

  open(filename) do fd
    buffer = Vector{UInt8}(undef, sizeof(MrParcRaidFileHeader))
    readbytes!(fd, buffer)
    ParcRaidHead =  reinterpret(MrParcRaidFileHeader, buffer)[1]

    @info ParcRaidHead.hdSize_
    @info ParcRaidHead.count_

    VBFILE = false

    if ParcRaidHead.hdSize_ > 32
        VBFILE = true

        # Rewind, we have no raid file header.
        # siemens_dat.seekg(0, std::ios::beg);
        seekstart(fd)

        ParcRaidHead.hdSize_ = ParcRaidHead.count_
        ParcRaidHead.count_ = 1
    elseif ParcRaidHead.hdSize_ != 0
        # This is a VB line data file
        error("Only VD line files with MrParcRaidFileHeader.hdSize_ == 0 (MR_PARC_RAID_ALLDATA) supported.")
    end

  end

end


text:
# ---
# title: 172. Factorial Trailing Zeroes
# id: problem172
# author: zhwang
# date: 2022-03-25
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/factorial-trailing-zeroes/description/>
# hidden: true
# ---
#
# Given an integer `n`, return _the number of trailing zeroes in`n!`_.
#
# **Follow up:** Could you write a solution that works in logarithmic time
# complexity?
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 3
#     Output: 0
#     Explanation:  3! = 6, no trailing zero.
#
#
# **Example 2:**
#
#
#
#     Input: n = 5
#     Output: 1
#     Explanation:  5! = 120, one trailing zero.
#
#
# **Example 3:**
#
#
#
#     Input: n = 0
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 104`
#
#
## @lc code=start
using LeetCode

## simulate Mathematica's function NestWhileList
function nest_while_list(f::Function, val::T, chk::Function)::Vector{T} where T
    res = [val]
    while chk(val)
        val = f(val)
        push!(res, val)
    end
    res
end

trailing_zeroes(n::Int) = sum(nest_while_list(i->i÷5, n÷5, >(1)))
## equivalent to the following
## function trailing_zeroes(n::Int)::Int
##     res = 0
##     while n >=5
##         n ÷= 5
##         res += n
##     end
##     res
## end

## @lc code=end

## @lc test=start
@testset "172.factorial-trailing-zeroes.jl" begin
    @test all(trailing_zeroes(i) == 0 for i in 0:4)
    @test trailing_zeroes(5) == 1
    @test trailing_zeroes(1999) == 496
    @test trailing_zeroes(2000) == trailing_zeroes(2001) == 499
end
## @lc test=end

## @lc info=start
# link: [solution to 172]()
## @lc info=end


text:
using SequenceReduction

soccer_events=CSV.read("C://Users//tbunc//OneDrive//Data//soccer_events//events.csv"; copycols=true)

# describe(soccer_events)

# Ok so for my first cut I will combine side (home or away) with the basic event.

# so keep just a few of the cols.

soccer_events2 = soccer_events[1:Int(floor(nrow(soccer_events)/100)) ,[:id_odsp, :sort_order, :sequence_event]]

@time basic_lev(soccer_events2, sequence_id=:id_odsp, sequence_order=:sort_order);
@time basic_lev_multi(soccer_events2, sequence_id=:id_odsp, sequence_order=:sort_order);
# Ok so potentially could do some form of outlier removal for master sequence
soccer_events

first_trace = hcat(soccer_lev[2][1][3]...)
soccer_lev[2][1][3]
reduce(hcat,soccer_lev[2][1][3])
# Lets turn this into an array

using RecursiveArrayTools
using Clustering

VA = VectorOfArray(soccer_lev[2][1][3])
arr = convert(Array,VA)
trace_vals_1 = convert(Array{Float64,2},arr[1,:,:])

# R  = kmeans(trace_vals_1, 2);
# counts(R)
# assignments(R)

dbscan(trace_vals_1,.05)
trace_vals_1[1,:]


text:
# code to try to figure out problem with gammaz inequality in bad diamond I
# Claudia March 2015

# test optBL with Case F Bad Diamond I
# Claudia January 2015

## include("../examples/case_f_example.jl");
## parameters!(net)
## q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
## q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
## q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
## q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
## q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

## d = DataCF([q1,q2,q3,q4,q5]);
## extractQuartet!(net,d)

## df = writeExpCF(d.quartet)
## writetable("CaseF_output.csv",df)

## include("../examples/case_f_example2.jl");
## parameters!(net)
## q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
## q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
## q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
## q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
## q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

## d = DataCF([q1,q2,q3,q4,q5]);
## extractQuartet!(net,d)

## df = writeExpCF(d.quartet)
## writetable("CaseF_output2.csv",df)

include("../src/types.jl")
include("../src/functions.jl")

df = readtable("CaseF_output.csv")
df2 = readtable("CaseF_output2.csv") #longer branches
d = readTableCF(df)
d2 = readTableCF(df2)

# starting ht (gamma,t4,t5,t9)
ht = [0.1,1.,1.,1.]
ht = [0.3,10.,10.,1.] # crashed, with check inside logPseudoLik, does not crash!

tree = string("(((6:0.1,(4)11#H1:::",string(1-ht[1]),")1:",string(ht[3]),",(11#H1:::",string(ht[1]),",7))5:",string(ht[4])",8:0.1,10:0.1);") # Case F: bad diamond I
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
net = readTopologyUpdate("prueba_tree.txt");
printEdges(net)

net.ht
realht1 = [0.1,0.127,0.0285]
realht2 = [1.1,0.49,0.2]

@time optBL!(net,d,true)

# ----------------- NLopt example website -----------------------
# to do here: put website example and change starting point to close to the boundary
# also, do a small example of optimization with inequality x_n+x_n-1<=1 to send to the author of nlopt

using NLopt

count = 0 # keep track of # function evaluations

function myfunc(x::Vector, grad::Vector)
    println("myfunc with x $(x)")
    if length(grad) > 0
        grad[1] = 0
        grad[2] = 0.5/sqrt(x[2])
    end
    global count
    count::Int += 1
    println("f_$count($x)")
    sqrt(x[2])
end

function myconstraint(x::Vector, grad::Vector, a, b)
    println("myconstraint with x $(x)")
    if length(grad) > 0
        grad[1] = 3a * (a*x[1] + b)^2
        grad[2] = -1
    end
    val = (a*x[1] + b)^3 - x[2] #should be negative
    println("value of inequality: $(val)")
    #val<0 || error("val in myconstraint not negative")
    return val
end

#opt = Opt(:LD_MMA, 2)
opt = Opt(:LN_COBYLA, 2)
lower_bounds!(opt, [-Inf, 0.])
xtol_rel!(opt,1e-4)

min_objective!(opt, myfunc)
inequality_constraint!(opt, (x,g) -> myconstraint(x,g,2,0), 1e-8)
inequality_constraint!(opt, (x,g) -> myconstraint(x,g,-1,1), 1e-8)

x=[1.234,5.678]
x=[1.,2.]
(minf,minx,ret) = optimize(opt,x)
println("got $minf at $minx after $count iterations (returned $ret)")


text:
abstract type AbstractTerm{T} <: AbstractOperator{T} end

term_type(op) = typeof(op)
term_type(t::AbstractTerm) = throw(ArgumentError("`term_type(::$(typeof(t)))` not defined"))

# TODO Workaround to be removed later
function Rimu.BasisSetRep(term::AbstractTerm, add)
    sm, basis = Rimu.Hamiltonians.build_sparse_matrix_from_LO(term, add; nnzs=0)
    BasisSetRep(sm, basis, term)
end

"""
    ParticleCountTerm(fun)

The particle count term:

```math
\\sum_{σ,p} f(σ,p) \\hat{n}_{p,σ},
```

where ``f`` is the `fun`, ``σ`` the spin (component) index, and ``p`` the mode.
"""
struct ParticleCountTerm{F,T<:Real} <: AbstractTerm{T}
    fun::F
end
function ParticleCountTerm(fun::F) where {F}
    T = float(typeof(fun(1, 1)))
    return ParticleCountTerm{F,T}(fun)
end

term_type(::ParticleCountTerm) = ParticleCountTerm

CompositeAction(::ParticleCountTerm) = NoCompositeAction()
LOStructure(::ParticleCountTerm) = IsHermitian()

num_offdiagonals(::ParticleCountTerm, args...) = 0

function diagonal_element(op::ParticleCountTerm, address, map, comp=1)
    return sum(map) do index
        index.occnum * op.fun(comp, index.mode)
    end
end

"""
     OnsiteInteraction(fun)

The onsite interaciton term:

```math
\\sum_{σ,τ,p} f(σ,τ) \\hat{a}^†_{p,σ} \\hat{a}^†_{p,τ} \\hat{a}_{p,τ} \\hat{a}_{p,σ},
```

where ``f`` is the `fun`, ``σ`` and ``τ`` the spin (component) indices, and ``p`` the mode.
"""
struct OnsiteInteractionTerm{F,T<:Real} <: AbstractTerm{T}
    fun::F
end
function OnsiteInteractionTerm(fun::F) where {F}
    T = float(typeof(fun(1, 1)))
    return OnsiteInteractionTerm{F,T}(fun)
end

term_type(::OnsiteInteractionTerm) = OnsiteInteractionTerm

CompositeAction(::OnsiteInteractionTerm) = OneWayCompositeAction()
LOStructure(::OnsiteInteractionTerm) = IsDiagonal()

num_offdiagonals(op::OnsiteInteractionTerm, _, _) = 0
num_offdiagonals(op::OnsiteInteractionTerm, _, _, _, _) = 0

function diagonal_element(op::OnsiteInteractionTerm, ::BoseFS, map, comp=1)
    return op.fun(comp, comp) * sum(n -> n.occnum * (n.occnum - 1), map) / 2
end
function diagonal_element(op::OnsiteInteractionTerm, ::FermiFS, _, _)
    return 0.0
end
function diagonal_element(op::OnsiteInteractionTerm, _, _, map_a, map_b, (c_a, c_b))
    N1, N2 = length(map_a), length(map_b)
    i = j = 1
    value = 0.0
    while i ≤ N1 && j ≤ N2
        a = map_a[i]
        b = map_b[j]
        if a.mode == b.mode
            i += 1
            j += 1
            value += a.occnum * b.occnum
        elseif a.mode < b.mode
            i += 1
        else
            j += 1
        end
    end
    return op.fun(c_a, c_b) * value
end

"""
    NeighbourOneBodyTerm(fun)

The neighbour hopping term in 1D:

```math
\\sum_{σ,⟨j, i⟩} f(σ) \\hat{a}^†_j \\hat{a}_i,
```

where ``f`` is the `fun` and ``sigma`` the spin component index.
"""
struct NeighbourOneBodyTerm{F,T<:Real} <: AbstractTerm{T}
    fun::F
end
function NeighbourOneBodyTerm(fun::F) where {F}
    T = float(typeof(fun(1)))
    return NeighbourOneBodyTerm{F,T}(fun)
end

term_type(::NeighbourOneBodyTerm) = NeighbourOneBodyTerm

CompositeAction(::NeighbourOneBodyTerm) = NoCompositeAction()
LOStructure(::NeighbourOneBodyTerm) = IsDiagonal()

function num_offdiagonals(op::NeighbourOneBodyTerm, _, map)
    return 2 * length(map)
end

function diagonal_element(op::NeighbourOneBodyTerm, args...)
    return 0.0
end

function get_offdiagonal(op::NeighbourOneBodyTerm, add, map, i, comp=1)
    M = num_modes(add)
    p, k = fldmod1(i, 2)
    src = map[p]
    dst = find_mode(add, mod1(src.mode + ifelse(k == 1, -1, 1), M))
    new_add, value = excitation(add, (dst,), (src,))
    return new_add, -op.fun(comp) * value
end

"""
    FullOneBodyTerm(fun)

The full one-body term:

```math
\\sum_{σ,p,q} f(σ,p,q) \\hat{a}^†_{σ,q} \\hat{a}_{σ,p}
```

where ``f`` is the `fun`, ``σ`` the spin (component) index, and ``p``, ``q`` the source and
target modes.
"""
struct FullOneBodyTerm{F,T,ADJ} <: AbstractTerm{T}
    fun::F
end
function FullOneBodyTerm(fun::F) where {F}
    T = float(typeof(fun(1, 1, 1)))
    return FullOneBodyTerm{F,T,false}(fun)
end

term_type(::FullOneBodyTerm) = FullOneBodyTerm

CompositeAction(::FullOneBodyTerm) = NoCompositeAction()

# TODO: break this out or find a way to signal that the function is Hermitian?
function LOStructure(op::FullOneBodyTerm{F}) where {F}
    if F <: ConstFunction{<:Real} || F <: ParameterColumnFunction{<:Any,<:Real}
        return IsHermitian()
    else
        return AdjointKnown()
    end
end

isadjoint(::FullOneBodyTerm{<:Any,<:Any,Adjoint}) where {Adjoint} = Adjoint

function Base.adjoint(op::FullOneBodyTerm{F,T,Adjoint}) where {F,T,Adjoint}
    FullOneBodyTerm{F,T,!Adjoint}(op.fun)
end

num_offdiagonals(::FullOneBodyTerm, add, map) = length(map) * (num_modes(add) - 1)

function diagonal_element(op::FullOneBodyTerm, add, map, comp=1)
    result = sum(map) do index
        op.fun(comp, index.mode, index.mode) * float(index.occnum)
    end
    if isadjoint(op)
        return conj(result)
    else
        return result
    end
end

function get_offdiagonal(op::FullOneBodyTerm, add, map, i, comp=1)
    M = num_modes(add)
    p, q = fldmod1(i, M - 1)
    p_index = map[p] # p-th occupied mode in add
    q += q ≥ p_index.mode # leave out diagonal matrix element
    q_index = find_mode(add, q) # q-th mode in add (not counting p)
    new_add, value = excitation(add, (q_index,), (p_index,)) # a_q^† a_p |add⟩
    if !isadjoint(op)
        return new_add, op.fun(comp, p_index.mode, q) * value
    else
        return new_add, conj(op.fun(comp, q, p_index.mode)) * value
    end
end

"""
    MomentumTwoBodyTerm(f; fold=true)

The momentum transfer term:

```math
\\sum_{σ,τ,p,q,k} f(σ,τ,p,q,k) a^†_{p + k,σ} a^†{q - k,τ} a_{q,τ} a_{p,σ}
```

In the sum above, ``σ`` and ``τ`` are spin indices, while ``p``, ``q``, and ``k`` are
integer mode indices.

If `fold` is set, transfers that would go out of the Brillouin zone are folded back in.
"""
struct MomentumTwoBodyTerm{F,T,Fold,Adjoint} <: AbstractTerm{T}
    fun::F
end
function MomentumTwoBodyTerm(fun::F; fold=true) where {F}
    T = float(typeof(fun(1, 1, 1, 1, 1)))
    return MomentumTwoBodyTerm{F,T,fold,false}(fun)
end
function MomentumTwoBodyTerm(val::Number=1; kwargs...)
    return MomentumTwoBodyTerm(ConstFunction(float(val)); kwargs...)
end

term_type(::MomentumTwoBodyTerm) = MomentumTwoBodyTerm
_fold(::MomentumTwoBodyTerm{<:Any,<:Any,Fold}) where {Fold} = Fold
isadjoint(::MomentumTwoBodyTerm{<:Any,<:Any,<:Any,Adjoint}) where {Adjoint} = Adjoint

function Base.show(io::IO, op::MomentumTwoBodyTerm)
    print(IOContext(io, :compact => true), "MomentumTwoBodyTerm(")
    if op.fun isa ConstFunction
        print(io, "$(op.fun.value);")
    else
        print(io, "$(op.fun);")
    end
    print(io, "fold=$(_fold(op)), isadjoint=$(isadjoint(op))")
end

CompositeAction(::MomentumTwoBodyTerm) = OneWayCompositeAction()

function LOStructure(::MomentumTwoBodyTerm{F}) where {F}
    if F <: ConstFunction{<:Real} || F <: InteractionMatrixFunction{<:Any,<:Real}
        return IsHermitian()
    else
        return AdjointKnown()
    end
end
function Base.adjoint(op::MomentumTwoBodyTerm{F,T,Fold,Adjoint}) where {F,T,Fold,Adjoint}
    return MomentumTwoBodyTerm{F,T,Fold,!Adjoint}(op.fun)
end

# Signle bosonic component
function num_offdiagonals(::MomentumTwoBodyTerm, bs::BoseFS, map)
    singlies = length(map)
    doublies = count(i -> i.occnum ≥ 2, map)
    M = num_modes(bs)
    return singlies * (singlies - 1) * (M - 2) + doublies * (M - 1)
end
function get_offdiagonal(op::MomentumTwoBodyTerm, bs::BoseFS, map, i, comp=1)
    new_add, val, p, q, k = momentum_transfer_excitation(bs, i, map; fold=_fold(op))
    val = if !isadjoint(op)
        op.fun(comp, comp, p, q, k) * val
    else
        error("adjoint is not supported for bosons")
        conj(op.fun(comp, comp, q + k, p - k, -k)) * val
    end

    return new_add, val/2
end
function diagonal_element(op::MomentumTwoBodyTerm, bs::BoseFS, map, comp=1)
    onproduct_zero = zero(eltype(op))
    onproduct_nonzero = zero(eltype(op))
    for i in 1:length(map)
        occ_i = map[i].occnum
        p = map[i].mode
        onproduct_zero += occ_i * (occ_i - 1) * op.fun(comp, comp, p, p, 0)
        for j in 1:i-1
            occ_j = map[j].occnum
            q = map[j].mode
            k = p - q
            onproduct_nonzero += op.fun(comp, comp, p, q, k) * 4 * occ_i * occ_j
        end
    end
    if isadjoint(op)
        return conj(onproduct_zero + onproduct_nonzero) / 2
    else
        return (onproduct_zero + onproduct_nonzero) / 2
    end
end

# Single fermionic component has no contributions.
num_offdiagonals(::MomentumTwoBodyTerm, ::FermiFS, _) = 0
diagonal_element(::MomentumTwoBodyTerm{<:Any,T}, ::FermiFS, map, comp=1) where {T} = zero(T)

# Cross-component part
function num_offdiagonals(::MomentumTwoBodyTerm, add_a, add_b, map_a, map_b)
    M = num_modes(add_a)
    return length(map_a) * length(map_b) * (M - 1)
end
function get_offdiagonal(op::MomentumTwoBodyTerm, add_a, add_b, map_a, map_b, i, (c_a, c_b))
    new_add_a, new_add_b, val, p, q, k = momentum_transfer_excitation(
        add_a, add_b, i, map_a, map_b, fold=_fold(op)
    )
    val = if !isadjoint(op)
        op.fun(c_a, c_b, p, q, k) * val
    else
        conj(op.fun(c_b, c_a, p - k, q + k, -k)) * val
    end
    return new_add_a, new_add_b, val
end
function diagonal_element(op::MomentumTwoBodyTerm, _, _, map_a, map_b, (c_a, c_b))
    onproduct_zero = 0.0
    onproduct_nonzero = 0.0
    for i in 1:length(map_a)
        occ_i = map_a[i].occnum
        p = map_a[i].mode
        for j in 1:length(map_b)
            occ_j = map_b[j].occnum
            q = map_b[j].mode
            onproduct_nonzero += (
                (op.fun(c_a, c_b, p, q, 0) + op.fun(c_b, c_a, q, p, 0))/2 * occ_i * occ_j
            )
        end
    end
    return onproduct_nonzero
end

"""
    FullTwoBodyTerm(fun)

```math
\\sum_{σ,τ,p,q,r,s} f(σ,τ,p,q,r,s)
    \\hat{a}^†_{s,σ} \\hat{a}^†{r,τ} \\hat{a}_{q,τ} \\hat{a}_{p,σ}
```

where ``f`` is the `fun`, ``σ`` and ``tau`` the spin component indices, and ``p,q,r,s`` the
modes.

If ``σ = τ`` and the component is fermionic, this term produces no elements.
"""
struct FullTwoBodyTerm{F,T,Adjoint} <: AbstractTerm{T}
    fun::F
end
function FullTwoBodyTerm(fun::F) where {F}
    T = float(typeof(fun(1,1,1,1,1,1)))
    return FullTwoBodyTerm{F,T,false}(fun)
end

term_type(::FullTwoBodyTerm) = FullTwoBodyTerm
isadjoint(::FullTwoBodyTerm{<:Any,<:Any,Adjoint}) where {Adjoint} = Adjoint

function Base.show(io::IO, op::FullTwoBodyTerm)
    print(IOContext(io, :compact => true), "FullTwoBodyTerm(")
    if op.fun isa ConstFunction
        print(io, "$(op.fun.value); ")
    else
        print(io, "$(op.fun); ")
    end
    print(io, "fold=$(_fold(op)), isadjoint=$(isadjoint(op)))")
end

CompositeAction(::FullTwoBodyTerm) = OneWayCompositeAction()

function LOStructure(::FullTwoBodyTerm{F}) where {F}
    if F <: ConstFunction{<:Real} || F <: InteractionMatrixFunction{<:Any,<:Real}
        return IsHermitian()
    else
        return AdjointKnown()
    end
end
function Base.adjoint(op::FullTwoBodyTerm{T,F,Adjoint}) where {T,F,Adjoint}
    return FullTwoBodyTerm{T,F,!Adjoint}(op.fun)
end

# TODO: Signle bosonic component not implemented

# Single fermionic component has no contributions.
num_offdiagonals(::FullTwoBodyTerm, ::FermiFS, map, comp=1) = 0
diagonal_element(::FullTwoBodyTerm{<:Any,T}, ::FermiFS, _, _) where {T} = zero(T)

function full_two_body_excitation(add_a, add_b, i, map_a, map_b)
    N1 = length(map_a)
    N2 = length(map_b)
    M = num_modes(add_a)
    s, r, q, p = Tuple(CartesianIndices((M, M, N2, N1))[i])
    s_index = find_mode(add_a, s)
    r_index = find_mode(add_b, r)
    q_index = map_b[q]
    p_index = map_a[p]

    new_add_a, val_a = excitation(add_a, (s_index,), (p_index,))
    new_add_b, val_b = excitation(add_b, (r_index,), (q_index,))

    return new_add_a, new_add_b, val_a * val_b, p_index.mode, q_index.mode, r, s
end

# Cross-component part
function num_offdiagonals(::FullTwoBodyTerm, add_a, add_b, map_a, map_b)
    M = num_modes(add_a)
    return length(map_a) * length(map_b) * (M - 1) * (M - 1)
end
function get_offdiagonal(op::FullTwoBodyTerm, add_a, add_b, map_a, map_b, i, (c_a, c_b))
    new_add_a, new_add_b, val, p, q, r, s = full_two_body_excitation(
        add_a, add_b, i, map_a, map_b,
    )
    val = if !isadjoint(op)
        op.fun(c_a, c_b, p, q, r, s) * val
    else
        error("TODO")
        #conj(op.fun(c_b, c_a, r, s, p, q)) * val ?
    end
    return new_add_a, new_add_b, val
end
function diagonal_element(
    op::FullTwoBodyTerm, _, _, map_a, map_b, (c_a, c_b)
)
    # This is equivalent to ∑_{p,q} fun(σ,τ,p,q,q,p) n_{p,σ} n_{q,τ}
    onproduct = zero(eltype(op))
    for i in map_a
        for j in map_b
            p = i.mode
            q = j.mode
            onproduct += op.fun(c_a, c_b, p, q, q, p) * i.occnum * j.occnum
        end
    end
    return onproduct
end

using Rimu.Hamiltonians: transcorrelated_three_body_excitation

"""
    MomentumThreeBodyTerm(fun)

The momentum-preserving three body term:

```math
\\sum_{σ,τ,p,q,s,k,k'} f(σ,τ,k,k')
    \\hat{a}^†_{p-k,σ} \\hat{a}^†_{q+k',σ} \\hat{a}^†_{s+k-k',τ}
    \\hat{a}_{s,τ} \\hat{a}_{q,σ} \\hat{a}_{p,σ},
```

Where...

TODO: make `f` accept `p`, `q`, and `s`?
"""
struct MomentumThreeBodyTerm{T,F,Adjoint} <: AbstractTerm{T}
    fun::F
end

function MomentumThreeBodyTerm(fun; fold=false, isadjoint=false)
    if fold
        throw(ArgumentError("folding three body terms not implemented"))
    end
    T = float(typeof(fun(1, 1, 1, 1)))
    return MomentumThreeBodyTerm{T,typeof(fun),isadjoint}(fun)
end

term_type(::MomentumThreeBodyTerm) = MomentumThreeBodyTerm

LOStructure(::MomentumThreeBodyTerm) = IsHermitian() # TODO: or is it?
CompositeAction(::MomentumThreeBodyTerm) = TwoWayCompositeAction()

# Single bosonic component not implemented

# Single fermionic component has no contributions
num_offdiagonals(::MomentumThreeBodyTerm, ::FermiFS, _) = 0
diagonal_element(::MomentumThreeBodyTerm{T}, ::FermiFS, map, c=1) where {T} = zero(T)

# Cross-component part
function num_offdiagonals(::MomentumThreeBodyTerm, add_a, add_b, map_a, map_b)
    M = num_modes(add_a)
    N1 = length(map_a)
    N2 = length(map_b)

    return (N1 - 1) * N1 * N2 * M^2
end

function get_offdiagonal(
    op::MomentumThreeBodyTerm, add_a, add_b, map_a, map_b, i, (c_a, c_b)
)
    # TODO should be p,q,s,k,l - but rimu only returns k and l (momentum changes)
    new_add_a, new_add_b, val, k, l = transcorrelated_three_body_excitation(
        add_a, add_b, i, map_a, map_b
    )
    val *= op.fun(c_a, c_b, k, l)
    # TODO conj if we need it.
    return new_add_a, new_add_b, val
end

function diagonal_element(op::MomentumThreeBodyTerm, _, _, map_a, map_b, (c_a, c_b))
    value = zero(eltype(op))
    for p in 1:length(map_a)
        for q in 1:p-1
            k = map_a[p].mode - map_a[q].mode
            qkk = op.fun(c_a, c_b, -k, k)
            value += 2 * qkk * length(map_b)
        end
    end
    return value
end


text:
module MaxHelpingHandCustomizableGlassExitBlock

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/CustomizableGlassExitBlock" CustomizableGlassExitBlock(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight, behindFgTiles::Bool=false)

const placements = Ahorn.PlacementDict(
    "Customizable Glass Exit Block (max480's Helping Hand)" => Ahorn.EntityPlacement(
        CustomizableGlassExitBlock,
        "rectangle"
    )
)

Ahorn.minimumSize(entity::CustomizableGlassExitBlock) = 8, 8
Ahorn.resizable(entity::CustomizableGlassExitBlock) = true, true

Ahorn.selection(entity::CustomizableGlassExitBlock) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CustomizableGlassExitBlock, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (1.0, 1.0, 1.0, 0.5), (1.0, 1.0, 1.0, 0.5))
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
using HarmonicBalance
using Symbolics
using Test

@variables Ω,γ,λ,F, x,θ,η,α, ω0, ω,t,T, ψ
@variables x(t)

natural_equation =  d(d(x,t),t) + γ*d(x,t) + Ω^2*(1-λ*cos(2*ω*t+ψ))*x + α * x^3 +η *d(x,t) * x^2
forces =  F*cos(ω*t+θ)
dEOM = HarmonicBalance.DifferentialEquation(natural_equation + forces, x)
HarmonicBalance.add_harmonic!(dEOM, x, ω)
harmonic_eq = HarmonicBalance.get_harmonic_equations(dEOM, slow_time=T, fast_time=t);
p = HarmonicBalance.Problem(harmonic_eq);

fixed_parameters = (Ω => 1.0,γ => 1E-2, λ => 5E-2, F => 1E-3,  α => 1.,  η=>0.3, θ => 0, ψ => 0)
sweep = ω => LinRange(0.9, 1.1, 100)
soln = HarmonicBalance.get_steady_states(p, sweep, fixed_parameters)

# save the result, try and load in the next step
#current_path = @__DIR__
#HarmonicBalance.save(current_path * "/parametron_result.jld2", soln)

###
# COMPARE TO KNOWN RESULTS
###
@variables u1, v1
ref1 = (Ω^2)*u1 + F*cos(θ) + γ*Differential(T)(u1) + (3//4)*α*(u1^3) + γ*ω*v1 + (2//1)*ω*Differential(T)(v1) + (1//4)*η*ω*(v1^3) + (3//4)*η*(u1^2)*Differential(T)(u1) + (1//4)*η*(v1^2)*Differential(T)(u1) + (3//4)*α*(v1^2)*u1 + (1//4)*η*ω*(u1^2)*v1 + (1//2)*η*u1*v1*Differential(T)(v1) + (1//2)*λ*(Ω^2)*v1*sin(ψ) - (ω^2)*u1 - (1//2)*λ*(Ω^2)*u1*cos(ψ)
ref2 = γ*Differential(T)(v1) + (Ω^2)*v1 + (3//4)*α*(v1^3) + (3//4)*α*(u1^2)*v1 + (1//4)*η*(u1^2)*Differential(T)(v1) + (3//4)*η*(v1^2)*Differential(T)(v1) + (1//2)*λ*(Ω^2)*v1*cos(ψ) + (1//2)*η*u1*v1*Differential(T)(u1) + (1//2)*λ*(Ω^2)*u1*sin(ψ) - F*sin(θ) - (ω^2)*v1 - (2//1)*ω*Differential(T)(u1) - γ*ω*u1 - (1//4)*η*ω*(u1^3) - (1//4)*η*ω*(v1^2)*u1

averaged = HarmonicBalance._remove_brackets(harmonic_eq)
@test isequal(simplify(expand(averaged[1]-ref1)), 0)
@test isequal(simplify(expand(averaged[2]-ref2)), 0)


text:
#delta function
function kd(x,y)
  if (x-y)==0
    return 1
  else
    return 0
  end
end


text:
"""

"""
abstract type AbstractAreaOfInfluence end

struct AOIProject2d <: AbstractAreaOfInfluence end


text:
# Sequencing process particle centered.

push!(LOAD_PATH, "./models")
push!(LOAD_PATH, ".")
push!(LOAD_PATH, "./SparseInverseProblems/src")

using SuperResModels
using SparseInverseProblemsMod
using Distributions
using QuadGK
using Roots

# To see the progress in pmap
# To use PmapProgressMeter you need to clone it manually
# do in Julia: Pkg.clone("https://github.com/slundberg/PmapProgressMeter.jl")
# if this is no longer existent, uncoment the alternative pmap in this code.
@everywhere using ProgressMeter
@everywhere using PmapProgressMeter
@everywhere using Interpolations
@everywhere using Vessel_network

using PyCall
@pyimport numpy as np
@pyimport matplotlib.pyplot as plt

@everywhere begin
    include("2d_parameters.jl")
    parameters = SuperResModels.Conv2dParameters(x_max, x_max, filter, filter_dx, filter_dy, sigma, sigma, sigma, sigma)
    model_static = SuperResModels.Conv2d(parameters)
    dynamic_parameters = SuperResModels.DynamicConv2dParameters(K, tau, v_max, 20)
    model_dynamic = SuperResModels.DynamicConv2d(model_static, dynamic_parameters)
    n_x = model_static.n_x
    n_y = model_static.n_y
end

# Get time of script start
now_str = string(now())
now_str = replace(now_str, ":", "-")
now_str = replace(now_str, ".", "-")

# A single static particle.
thetas = reshape([x_max/2; x_max/2], 2, 1)
weights = [1.0]

# Phi represents the forward operator.
single_particle_norm = norm(phi(model_static, thetas, weights), lp_norm)
@everywhere single_particle_norm = $single_particle_norm

# The examples we seek to emulate correspond to two sets of particles moving along a two stripes on each side.
# So particles_m stands for those located on the west side, moving towards the south.
# The vector particles_m basically are the locations at time 0 of the particles, it also defines quantity.
# Particles_p are the symmetric equivalent, on the east and moving towards the north.

Nparticles = 500
activation_prob = 0.0002
expected_time = 34

### Generate sequence ###
println("Generating sequence... ")
video = SharedArray{Float64}(n_x * n_y, n_im)

function activate_particles(Nparticles,activation_prob, expected_time, Vessel)
    # particles are tuples of position,velocity, weight, survival time.
    # The position is a one dimensional variable, that once passed to the vessel
    # function and displaced accordingly, would return the respective position.
    B = Binomial(Nparticles,activation_prob)
    P = Binomial(1,0.5)
    N = Poisson(expected_time)

    N_new_particles = rand(B)
    new_particles = []
    for i in 1:N_new_particles
        # decide if it is flowing downwards (=0) or upwards (=1)
        direction = rand(P)
        position = Uniform_sampler(Vessel)
        weight = 1
	    velocity = (direction-0.5)*2*particle_velocity
        time = rand(N)
        push!(new_particles,(position, velocity, weight, time))
    end
    return new_particles
end

function time_step(Vessel,particles)
    for j in length(particles):-1:1
	particle = particles[j]
	# update position
	new_position = displacement(Vessel, particle[1][1],
						   particle[1][2], tau*particle[2])
	# update time
	new_time = particle[4]-1
	# discard if time is over or displaced outside the network
	if (new_time <= 0) | isnan(new_position[1])
	    deleteat!(particles,j)
	else
	    particles[j]=(new_position,particle[2], particle[3], new_time)
	end
    end
end
particles = []

for i in 1:n_im
    # Move particles forward in time
    time_step(Vessel, particles)
    # Activate new particles and push them inside the list of particles
    new_particles = activate_particles(Nparticles,activation_prob,expected_time, Vessel)
    for new_particle in new_particles
	push!(particles, new_particle)
    end
    # Represent the static particles to input into the phi function
    # (the forward operator). Thetas are 2d vectors representing position
    thetas = zeros(2,length(particles))
    weights = zeros(length(particles))
    for j in 1:length(particles)
    # The case in the upward or downward vessel
    coords = evaluation(Vessel, particles[j][1][1] , particles[j][1][2])
    if particles[j][2]>0
	    thetas[1,j] = coords[1]
	    thetas[2,j] = coords[2]
    else
        thetas[1,j] = coords[1]+dx
        thetas[2,j] = coords[2]-dx
    end
	weights[j] = particles[j][3]
    end
    ## Image the particles
    # include noise
    video[:,i] = sigma_noise*randn(size(video[:,i]))
    if (length(weights) > 0)
        video[:,i] += phi(model_static, thetas, weights)
    end
end

### Obtaining time frames in which the quantity of particles remained constant
println("Getting short sequences without jump...")
# get the total mass at each time step
frame_norms = [norm(video[:, i], lp_norm) for i in 1:n_im]
@everywhere frame_norms = $frame_norms
# Find locations in which there was a significative mass difference between
# two consecutive time steps.
jumps = find(abs.(frame_norms[2:end] - frame_norms[1:end-1]) .>
                                            jump_threshold*single_particle_norm)
jumps = [0; jumps; n_im]
# Obtain non-overlapping intervals of at least 2K+1 consecutive time samples
# in which the mass didn't changed.
short_seqs = []
for i in 1:length(jumps)-1
    append!(short_seqs, [(jumps[i] + 5*j + 1):(jumps[i] + 5*j + 5)
                                for j in 0:(div(jumps[i+1] - jumps[i], 5) -1)])
end
println("The number of jumps is: ", length(jumps))

### Function that given the a subquence of the total video, will estimate the
### locations and weights of the involved particles.
@everywhere function posvel_from_seq(video, seq)
    assert(length(seq) == 5)
    target = video[:,seq][:]
    # estimated number of particles in the sequence.
    est_num_particles = div(frame_norms[seq[1]], single_particle_norm*0.95)
    if (est_num_particles == 0)
        return Matrix{Float64}(5,0)
    end
    # Function required to use the SparseInverseProblems.ADCG method.
    function callback(old_thetas, thetas, weights, output, old_obj_val)
        #evalute current OV
        new_obj_val,t = SparseInverseProblemsMod.loss(
                             SparseInverseProblemsMod.LSLoss(), output - target)
        #println("gap = $(old_obj_val - new_obj_val)")
        if old_obj_val - new_obj_val < 1E-4
            return true
        end
        return false
    end
    # It uses the ACDG algorithm to estimate the location and weights of the
    # particles. With the estimated number of particles we can bound the total
    # variation on the solutions.
    (thetas_est,weights_est) = SparseInverseProblemsMod.ADCG(model_dynamic,
                                    SparseInverseProblemsMod.LSLoss(), target,
                                        frame_norms[seq[1]], callback=callback,
                                                                  max_iters=100)
    if length(thetas_est) > 0
        println("est_num = ", est_num_particles)
        println("thetas = ", thetas_est)
        println("weights = ", weights_est)
        return [thetas_est; weights_est']
    else
        print("Failed to obtain a suitable solution")
        return Matrix{Float64}(5,0)
    end
end

println("Inverting...")

all_thetas = pmap(seq -> begin sleep(1); posvel_from_seq(video, seq) end,
                                       Progress(length(short_seqs)), short_seqs)

println("Reprojecting...")

errors = zeros(length(short_seqs))
### Measurements error, we simulate the measurements that would be obtained
### with our reconstructed values

for i in 1:length(short_seqs)
    seq=  short_seqs[i]
    target = video[:, seq][:]
    if length(all_thetas[i]) > 0
        reprojection = phi(model_dynamic, all_thetas[i][1:4,:]
                                            ,all_thetas[i][5,:])
        println("error = ", norm(target-reprojection))
        errors[i] = norm(target-reprojection)
    else
        errors[i] = norm(target)
    end
end

### Save the simulated data ###
data_folder = "data/2Dsimulations/"*now_str

using PyCall
@pyimport numpy as np
mkdir(data_folder)
cp("2d_parameters.jl", string(data_folder, "/2d_parameters.jl"))
cd(data_folder)
short_seq_array = hcat(short_seqs...)
np.save("video", video)
np.save("frame_norms", frame_norms)
np.save("jumps", jumps)
np.save("short_seq_array", short_seq_array)
for i in 1:length(short_seqs)
    np.save(string("thetas-", i), all_thetas[i])
end
np.save("errors", errors)
cd("..")

#
# if false
# for i = 1:500
# plt.figure()
# plt.pcolormesh(np.linspace(0, 0.01, n_x), np.linspace(0, 0.01, n_x), reshape(video[:,i],n_x,n_y))
# plt.colorbar()
# plt.show()
# end
# end

# showFrames = 1:500
# plt.figure()
# # Compute the L2 norm at each time sample
# L2norms = sqrt.(sum(video[:,showFrames].^2,1))[:]
# plt.plot(1:length(L2norms), L2norms)
# # Paint the close to constant cases
# minSnapshot = 5
# tolerance = 0.05
# j = 1
# i = 1
# while i <= length(L2norms)-minSnapshot
#     while (j+i <= length(L2norms)-minSnapshot) & (abs.(L2norms[i]-L2norms[i+j])<=tolerance)
# 	j = j +1
#     end
#     if j-1 >= minSnapshot
# 	plt.plot([i, i+j-1], [mean(L2norms[i:i+j-1]), mean(L2norms[i:i+j-1])],color="r")
#     end
#     i = i+j
#     j = 1
# end
# plt.show()


text:
@testset "1163.last-substring-in-lexicographical-order.jl" begin
    @test last_substring("arpzrkqizlvqouudiggnezejhhuaanmaickndesynbkscubboihueuatpbhclzifkrcjhncllfdyqeawmdcefyhcnzciypnrezcu") ==
          "zrkqizlvqouudiggnezejhhuaanmaickndesynbkscubboihueuatpbhclzifkrcjhncllfdyqeawmdcefyhcnzciypnrezcu"
    @test last_substring("abab") == "bab"
    @test last_substring("leetcode") == "tcode"
    @test last_substring("cccc") == "cccc"
    @test last_substring("ccacccacccc") == "cccc"
    @test last_substring("caccacccaa") == "cccaa"
    s = "a" ^ 4000000
    @test last_substring(s) == s
end