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

using Markdown
using InteractiveUtils

# ╔═╡ c9762a40-3119-11eb-021f-3526881bee53
begin
   using MLJ, MLJFlux, Flux
   import Random.seed!; seed!(123)
   using Plots
end

# ╔═╡ ffc72f20-3184-11eb-0f84-f57d92c16b65
md"### Aqui construímos nosso próprio modelo usando MLJ+Flux.

É necessário a criação de uma estrutura que herda de um builder. Adicionalmente, devemos sobrecarregar a função build, pois esta será chamada internamente."

# ╔═╡ 563a4730-311b-11eb-0b9f-45f44f995aa6
begin

	mutable struct MNISTBuilder <: MLJFlux.Builder
		n_hidden::Int
	end

	function MLJFlux.build(builder::MNISTBuilder, (n1, n2), m, n_channels)
		return Chain(
			flatten,
			Dense(n1 * n2, builder.n_hidden, relu),
			Dense(builder.n_hidden, m),
		)
	end

end

# ╔═╡ c7f7cf90-3185-11eb-0f94-f5e182753280
md"### Instanciando Classificador"

# ╔═╡ 5c511db0-311b-11eb-3eec-2fb8f22a5a53
begin
	@load ImageClassifier
	clf = ImageClassifier(;
		builder=MNISTBuilder(128),
		optimiser=ADAM(0.001),
		loss=Flux.crossentropy,
		epochs=6,
		batch_size=128,
	)
end

# ╔═╡ aa9142ae-3185-11eb-22a3-f39504976138
md"### O problema

url: http://yann.lecun.com/exdb/mnist/

The MNIST database of handwritten digits, available from this page, has a training set of 60,000 examples, and a test set of 10,000 examples. It is a subset of a larger set available from NIST. The digits have been size-normalized and centered in a fixed-size image.
It is a good database for people who want to try learning techniques and pattern recognition methods on real-world data while spending minimal efforts on preprocessing and formatting."

# ╔═╡ 5b631430-311b-11eb-1388-6f434e9cc354
begin
	X = Flux.Data.MNIST.images()
	y = Flux.Data.MNIST.labels()
	y = coerce(y, Multiclass)
end

# ╔═╡ 130f6240-311d-11eb-0efd-7b4822376ef8
X[100]

# ╔═╡ e40bc4c0-3185-11eb-0794-63633f2c5269
md"### Particionando os dados

Realizamos aqui o que conhecemos como hold-out. Particionamos em treino e teste, que é a forma mais simples de se fazer uma validação. Sabemos que existem outras no curso. No caso, deixei 70% para treino e 30% para teste. Como podem notar, a função gera índices qu correspondem às linhas dos exemplos. Ele não tira cópia dos dados."

# ╔═╡ fa3a71d0-3124-11eb-3b83-014b67169868
train, test = partition(eachindex(y), 0.01, shuffle=true, rng=1234)

# ╔═╡ d6fd7bc0-3185-11eb-13c0-7381eb7a588a
md"### Treinando a Rede

Um conceito interessante do MLJ é o machine. Esta função faz uma mapeamento do preditor ao dado. Quando aplicamos machine, estamos gerando um nó num grafo de computação e um estado associado. Existem algmas vantagens neste paradigma como: associar um mesmo preditor a vários dados distintos que podem ser processados em diferentes recursos; evitar re execução de um mesmo preditor e um mesmo dado."

# ╔═╡ 5d2b5110-311b-11eb-3dc3-0b7fdfd060fd
mach = machine(clf, X, y)

# ╔═╡ 177cb6c2-3186-11eb-2d73-39f8a0ac039b
fit!(mach,rows=train)

# ╔═╡ 3d71d0e0-3186-11eb-3a9c-d17db573ba54
md"### Predição

Aqui usamos o predict_mode para pegar as previsões discretas sem as probabilidades. às vezes é preciso usar apenas o predict para pegar ambas as classes e probabilidades, pois certas métricas necessitam"

# ╔═╡ c5fd5630-3124-11eb-2844-7fc187c8a99b
ŷ = predict_mode(mach, rows=test)

# ╔═╡ d0212970-3124-11eb-2041-91a20d4349a2
cm = confusion_matrix(ŷ,y[test])

# ╔═╡ d1e7b8a0-3124-11eb-254e-83e4e57c4982
accuracy(cm)

# ╔═╡ dfab9600-3124-11eb-2459-6f59178a5d84
balanced_accuracy(ŷ,y[test])

# ╔═╡ 9d2bb6ae-31ed-11eb-2878-074b2adae35e
md"### Encontrando o melhor modelo"

# ╔═╡ 2e03a3a2-31ee-11eb-03eb-35d112cba7a6
clf.builder.n_hidden

# ╔═╡ a894b420-31ed-11eb-2de1-6b5a80d19ea9
begin

	# defining hyperparams for tunning
	r1 = range(clf, :lambda, lower=0.001, upper=1.0)
	r2 = range(clf, :epochs, lower=10, upper=300)
	r3 = range(clf, :(builder.n_hidden), lower=64, upper=256)

	# attaching tune
	self_tuning_ann_model = TunedModel(model        =  clf,
									   resampling   = CV(nfolds = 3),
									   tuning       = RandomSearch(),
									   range        = [r1,r2,r3],
									   measure      = accuracy,
									   operation    = predict_mode,
									   n=8)

end

# ╔═╡ bd1883d0-31ee-11eb-0929-9f9e106252e9
begin
    # putting into the machine
	self_tuning_ann = machine(self_tuning_ann_model, X, y)

	# fitting with tunning
	fit!(self_tuning_ann,rows=train)
end

# ╔═╡ 2629c250-329c-11eb-332d-7debeb62dabe
fitted_params(self_tuning_ann)

# ╔═╡ 7f796040-329c-11eb-0cba-2f3c6648a92e
report(self_tuning_ann)

# ╔═╡ 7b4d6320-3186-11eb-3228-616d0e73b93c
md"### Aplicando o melhor modelo

Depois que testamos combinações com variações de parâmetros, vamos aplicar o modelo aprendido no conjunto de teste."

# ╔═╡ 5e5eef12-311b-11eb-23a9-6b715d078655
ŷb = predict_mode(self_tuning_ann,rows=test)

# ╔═╡ 8238a550-32b2-11eb-1d18-ed7f4b004b27
cmb = confusion_matrix(ŷb,y[test])

# ╔═╡ b35e47c0-32b2-11eb-1fce-cf420134aed0
accuracy(cmb)

# ╔═╡ b9377900-32b2-11eb-3478-e11f07450ec5
balanced_accuracy(ŷb,y[test])

# ╔═╡ Cell order:
# ╠═c9762a40-3119-11eb-021f-3526881bee53
# ╟─ffc72f20-3184-11eb-0f84-f57d92c16b65
# ╠═563a4730-311b-11eb-0b9f-45f44f995aa6
# ╟─c7f7cf90-3185-11eb-0f94-f5e182753280
# ╠═5c511db0-311b-11eb-3eec-2fb8f22a5a53
# ╟─aa9142ae-3185-11eb-22a3-f39504976138
# ╠═5b631430-311b-11eb-1388-6f434e9cc354
# ╠═130f6240-311d-11eb-0efd-7b4822376ef8
# ╟─e40bc4c0-3185-11eb-0794-63633f2c5269
# ╠═fa3a71d0-3124-11eb-3b83-014b67169868
# ╟─d6fd7bc0-3185-11eb-13c0-7381eb7a588a
# ╠═5d2b5110-311b-11eb-3dc3-0b7fdfd060fd
# ╠═177cb6c2-3186-11eb-2d73-39f8a0ac039b
# ╟─3d71d0e0-3186-11eb-3a9c-d17db573ba54
# ╠═c5fd5630-3124-11eb-2844-7fc187c8a99b
# ╠═d0212970-3124-11eb-2041-91a20d4349a2
# ╠═d1e7b8a0-3124-11eb-254e-83e4e57c4982
# ╠═dfab9600-3124-11eb-2459-6f59178a5d84
# ╟─9d2bb6ae-31ed-11eb-2878-074b2adae35e
# ╠═2e03a3a2-31ee-11eb-03eb-35d112cba7a6
# ╠═a894b420-31ed-11eb-2de1-6b5a80d19ea9
# ╠═bd1883d0-31ee-11eb-0929-9f9e106252e9
# ╠═2629c250-329c-11eb-332d-7debeb62dabe
# ╠═7f796040-329c-11eb-0cba-2f3c6648a92e
# ╟─7b4d6320-3186-11eb-3228-616d0e73b93c
# ╠═5e5eef12-311b-11eb-23a9-6b715d078655
# ╠═8238a550-32b2-11eb-1d18-ed7f4b004b27
# ╠═b35e47c0-32b2-11eb-1fce-cf420134aed0
# ╠═b9377900-32b2-11eb-3478-e11f07450ec5

max_stars_repo_name: lalvim/machinelearning

max_issues_repo_name: lalvim/machinelearning

max_forks_repo_name: lalvim/machinelearning


text:
@testset "Mincut" begin

    g = lg.CompleteDiGraph(5)
    cap1 = [
       0.0 2.0 2.0 0.0 0.0
       0.0 0.0 0.0 0.0 3.0
       0.0 1.0 0.0 3.0 0.0
       0.0 0.0 0.0 0.0 1.0
       0.0 0.0 0.0 0.0 0.0
    ];
    (part1, part2, value) = LightGraphsFlows.mincut(g,1,5,cap1,LightGraphsFlows.PushRelabelAlgorithm())
    @test value ≈ 4.0
    @test part1 == [1]
    @test sort(part2) == collect(2:5)
    cap2 = [
       0.0 3.0 2.0 0.0 0.0
       0.0 0.0 0.0 0.0 3.0
       0.0 1.0 0.0 3.0 0.0
       0.0 0.0 0.0 0.0 1.5
       0.0 0.0 0.0 0.0 0.0
    ];
    (part1, part2, value) = LightGraphsFlows.mincut(g,1,5,cap2,LightGraphsFlows.PushRelabelAlgorithm())
    @test value ≈ 4.5
    @test sort(part1) == collect(1:4)
    @test part2 == [5]

    g2 = lg.DiGraph(5)
    lg.add_edge!(g2,1,2)
    lg.add_edge!(g2,1,3)
    lg.add_edge!(g2,3,4)
    lg.add_edge!(g2,3,2)
    lg.add_edge!(g2,2,5)

    (part1, part2, value) = LightGraphsFlows.mincut(g2,1,5,cap1,LightGraphsFlows.PushRelabelAlgorithm())
    @test value ≈ 3.0
    @test sort(part1) == [1,3,4]
    @test sort(part2) == [2,5]

end

ext: jl


text:
crc3, parameters, environment, alphas, stepsymbol, cl_stall = MT.build_crc3(; nwings)
crc3_model = MT.AircraftDynamics(crc3, parameters, environment, stepsymbol, cl_stall)
ts = range(0, stop=0.1, length=20)
us = ones(length(ts)) .* 3000.0 * 2  * pi / 60
alpha = 5.0 * pi/180
gamma = 3.0 * pi/180
vinf = 5.0
x0 = MT.steady_climb_state(alpha, gamma, vinf)
nx = length(x0)

xs = MT.integrate_dynamics(crc3_model, x0, us, ts)
MT.plot_states(ts, xs, nx)

max_stars_repo_name: byuflowlab/MultifidelityTrajectories

max_issues_repo_name: byuflowlab/MultifidelityTrajectories

max_forks_repo_name: byuflowlab/MultifidelityTrajectories


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module Ontology
export Annotation, AnnotationID, ObAnnotation, HomAnnotation

using AutoHashEquals
using Reexport

using ..Doctrine

# Data types
############

""" Semantic annotation of computer program.

This type is agnostic to the programming language of the computer program. All
language-specific information is stored in the `language` dictionary.
"""
abstract type Annotation end

""" Unique identifer of annotation.
"""
@auto_hash_equals struct AnnotationID
  language::String
  package::String
  id::String
end

struct ObAnnotation <: Annotation
  name::AnnotationID
  language::Dict{Symbol,Any}
  definition::Monocl.Ob
  slots::Vector{Monocl.Hom}
end

struct HomAnnotation <: Annotation
  name::AnnotationID
  language::Dict{Symbol,Any}
  definition::Monocl.Hom
end

# Modules
#########

include("OntologyJSON.jl")
include("OntologyDBs.jl")
include("rdf/OntologyRDF.jl")

@reexport using .OntologyJSON
@reexport using .OntologyDBs
@reexport using .OntologyRDF

end

max_stars_repo_name: cclauss/semanticflowgraph

max_issues_repo_name: cclauss/semanticflowgraph

max_forks_repo_name: cclauss/semanticflowgraph


text:
######################################################################
# Custom Futures Unit Tests
# -----
# Licensed under Apache License 2.0
module TestFutures
using Test
using Pia.Futures
import ExtraFun: TimeoutError
import Pia: CancellationError

@testset "Futures" begin
    @testset "instantiate" begin
        @test Futures.ispending(Future())
        @test Futures.typevar(Future()) === Any
        @test Futures.issuccess(Future(42))
        @test Futures.typevar(Future(42)) === Int
        @test Futures.typevar(Future{Number}(42)) === Number
    end

    @testset "states" begin
        @testset "success {Any}" begin
            let future = Future()
                @test  Futures.ispending(future)
                @test !Futures.issuccess(future)
                @test !Futures.isfailure(future)
                @test !Futures.iscancelled(future)

                # Note: Julia automatically returns RHS on assignment regardless of setindex! return value
                @test ((future[] = 42);       future[] == 42)
                @test ((future[] = "foobar"); future[] == 42)

                @test !Futures.ispending(future)
                @test  Futures.issuccess(future)
                @test !Futures.isfailure(future)
                @test !Futures.iscancelled(future)
            end
        end

        @testset "success {Number}" begin
            let future = Future{Number}()
                @test  Futures.ispending(future)
                @test !Futures.issuccess(future)
                @test !Futures.isfailure(future)
                @test !Futures.iscancelled(future)

                # Note: Julia automatically returns RHS on assignment regardless of setindex! return value
                @test_throws MethodError future[] = "foobar"
                @test ((future[] = 42);    future[] == 42)
                @test ((future[] = 69.69); future[] == 42)

                @test !Futures.ispending(future)
                @test  Futures.issuccess(future)
                @test !Futures.isfailure(future)
                @test !Futures.iscancelled(future)
            end
        end
    end

    @testset "cancel & timeout" begin
        let future = Future()
            @async begin; sleep(0.01); cancel(future); end
            @test_throws CancellationError wait(future, timeout=1)
        end
        let future = Future()
            @test_throws TimeoutError wait(future, timeout=0.5)
        end
    end
end

end # module TestFutures


text:
using ModelGraphs
using JuMP
using GLPK

graph = ModelGraph()

function sub_model()
    m1 = Model()
    @variable(m1,x[1:2] >= 0)
    @constraint(m1,x[1] >= 3)
    @objective(m1,Min,x[1] + x[2])
    setsolver(m1,GLPKSolverLP())
    return m1
end

n2 = add_node!(graph,sub_model())
n3 = add_node!(graph,sub_model())

@linkvariable(graph,x[1:2] >= 0)
link_variables!(graph[:x][1],n2[:x][1])
link_variables!(graph[:x][2],n2[:x][2])
link_variables!(graph[:x][2],n3[:x][2])

optimize!(graph,with_optimizer(BendersOptimizer))


text:
# # Baroclinic adjustment
#
# In this example, we simulate the evolution and equilibration of a baroclinically
# unstable front.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, CairoMakie"
# ```

using Oceananigans
using Oceananigans.Units

# ## Grid
#
# We use a three-dimensional channel that is periodic in the `x` direction:

Lx = 1000kilometers # east-west extent [m]
Ly = 1000kilometers # north-south extent [m]
Lz = 1kilometers    # depth [m]

Nx = 64
Ny = 64
Nz = 40

grid = RectilinearGrid(CPU();
                       topology = (Periodic, Bounded, Bounded),
                       size = (Nx, Ny, Nz),
                       x = (0, Lx),
                       y = (-Ly/2, Ly/2),
                       z = (-Lz, 0),
                       halo = (3, 3, 3))

# ## Turbulence closures

# We prescribe the values of vertical viscocity and diffusivity according to the ratio
# of the vertical and lateral grid spacing.

Δx, Δz = Lx/Nx, Lz/Nz

𝒜 = Δz/Δx # Grid cell aspect ratio.

κh = 0.1    # [m² s⁻¹] horizontal diffusivity
νh = 0.1    # [m² s⁻¹] horizontal viscosity
κz = 𝒜 * κh # [m² s⁻¹] vertical diffusivity
νz = 𝒜 * νh # [m² s⁻¹] vertical viscosity

horizontal_diffusive_closure = HorizontalScalarDiffusivity(ν = νh, κ = κh)

vertical_diffusive_closure = VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretization();
                                                       ν = νz, κ = κz)
nothing #hide

# ## Model

# We built a `HydrostaticFreeSurfaceModel` with an `ImplicitFreeSurface` solver.
# Regarding Coriolis, we use a beta-plane centered at 45° South.

model = HydrostaticFreeSurfaceModel(; grid,
                                      coriolis = BetaPlane(latitude = -45),
                                      buoyancy = BuoyancyTracer(),
                                      tracers = :b,
                                      closure = (vertical_diffusive_closure, horizontal_diffusive_closure),
                                      momentum_advection = WENO5(),
                                      tracer_advection = WENO5(),
                                      free_surface = ImplicitFreeSurface())

# We want to initialize our model with a baroclinically unstable front plus some small-amplitude
# noise.

"""
    ramp(y, Δy)

Linear ramp from 0 to 1 between -Δy/2 and +Δy/2.

For example:
```
            y < -Δy/2 => ramp = 0
    -Δy/2 < y < -Δy/2 => ramp = y / Δy
            y >  Δy/2 => ramp = 1
```
"""
ramp(y, Δy) = min(max(0, y/Δy + 1/2), 1)
nothing #hide

# We then use `ramp(y, Δy)` to construct an initial buoyancy configuration of a baroclinically
# unstable front. The front has a buoyancy jump `Δb` over a latitudinal width `Δy`.

N² = 4e-6 # [s⁻²] buoyancy frequency / stratification
M² = 8e-8 # [s⁻²] horizontal buoyancy gradient

Δy = 50kilometers # width of the region of the front
Δb = Δy * M²      # buoyancy jump associated with the front
ϵb = 1e-2 * Δb    # noise amplitude

bᵢ(x, y, z) = N² * z + Δb * ramp(y, Δy) + ϵb * randn()

set!(model, b=bᵢ)

# Let's visualize the initial buoyancy distribution.

using CairoMakie

x, y, z = 1e-3 .* nodes((Center, Center, Center), grid) # convert m -> km

b = model.tracers.b

fig, ax, hm = heatmap(y, z, interior(b)[1, :, :],
                      colormap=:deep,
                      axis = (xlabel = "y [km]",
                              ylabel = "z [km]",
                              title = "b(x=0, y, z, t=0)",
                              titlesize = 24))

Colorbar(fig[1, 2], hm, label = "[m s⁻²]")

save("initial_buoyancy.svg", fig); nothing # hide

# ![](initial_buoyancy.svg)

# Now let's built a `Simulation`.

Δt₀ = 5minutes
stop_time = 40days

simulation = Simulation(model, Δt=Δt₀, stop_time=stop_time)

# We add a `TimeStepWizard` callback to adapt the siulation's time-step,

wizard = TimeStepWizard(cfl=0.2, max_change=1.1, max_Δt=20minutes)

simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(20))

# Also, we add a callback to print a message about how the simulation is going,

using Printf

wall_clock = [time_ns()]

function print_progress(sim)
    @printf("[%05.2f%%] i: %d, t: %s, wall time: %s, max(u): (%6.3e, %6.3e, %6.3e) m/s, next Δt: %s\n",
            100 * (sim.model.clock.time / sim.stop_time),
            sim.model.clock.iteration,
            prettytime(sim.model.clock.time),
            prettytime(1e-9 * (time_ns() - wall_clock[1])),
            maximum(abs, sim.model.velocities.u),
            maximum(abs, sim.model.velocities.v),
            maximum(abs, sim.model.velocities.w),
            prettytime(sim.Δt))

    wall_clock[1] = time_ns()

    return nothing
end

simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterval(100))

# ## Diagnostics/Output

# Add some diagnostics. Here, we save the buoyancy, ``b``, at the edges of our domain as well as
# the zonal (``x``) averages of buoyancy and zonal velocity ``u``.

u, v, w = model.velocities

B = Field(Average(b, dims=1))
U = Field(Average(u, dims=1))

filename = "baroclinic_adjustment"
save_fields_interval = 0.5day

slicers = (west = (1, :, :),
           east = (grid.Nx, :, :),
           south = (:, 1, :),
           north = (:, grid.Ny, :),
           bottom = (:, :, 1),
           top = (:, :, grid.Nz))

for side in keys(slicers)
    indices = slicers[side]

    simulation.output_writers[side] = JLD2OutputWriter(model, (; b, u);
                                                       filename = filename * "_$(side)_slice",
                                                       schedule = TimeInterval(save_fields_interval),
                                                       indices)
end

simulation.output_writers[:zonal] = JLD2OutputWriter(model, (b=B, u=U);
                                                     schedule = TimeInterval(save_fields_interval),
                                                     filename = filename * "_zonal_average")

# Now let's run!

@info "Running the simulation..."

run!(simulation)

@info "Simulation completed in " * prettytime(simulation.run_wall_time)

# ## Visualization

# Now we are ready to visualize our resutls! We use `CairoMakie` in this example.
# On a system with OpenGL `using GLMakie` is more convenient as figures will be
# displayed on the screen.

using CairoMakie

filename = "baroclinic_adjustment"

fig = Figure(resolution = (800, 500))
ax_b = fig[2, 1] = LScene(fig, show_axis=false)
nothing #hide

# We use Makie's `Observable` to animate the data. To dive into how `Observable`s work we
# refer to [Makie.jl's Documentation](https://makie.juliaplots.org/stable/documentation/nodes/index.html).

n = Observable(1)

# We load the saved buoyancy output on the top, bottom, and east surface as `FieldTimeSeries`es.

sides = keys(slicers)

slice_filenames = NamedTuple(side => filename * "_$(side)_slice.jld2" for side in sides)

b_timeserieses = (
      east = FieldTimeSeries(slice_filenames.east, "b"),
    bottom = FieldTimeSeries(slice_filenames.bottom, "b"),
       top = FieldTimeSeries(slice_filenames.top, "b")
)

b_avg_timeseries = FieldTimeSeries(filename * "_zonal_average.jld2", "b")

# We build the coordinates and we rescale the vertical coordinate for visualization purposes.

x, y, z = nodes(b_timeserieses[1])

yscale = 2.5
zscale = 600
z = z .* zscale
y = y .* yscale

zonal_slice_displacement = 1.5
nothing #hide

# Now let's make a 3D plot of the buoyancy and in front of it we'll use the zonally-averaged output
# to plot the instantaneous zonal-average of the buoyancy.

b_slices = (
      east = @lift(interior(b_timeserieses.east[$n], 1, :, :)),
    bottom = @lift(interior(b_timeserieses.bottom[$n], :, :, 1)),
       top = @lift(interior(b_timeserieses.top[$n], :, :, 1))
)

clims_b = @lift 1.1 .* extrema(b_timeserieses.top[$n][:])
kwargs_b = (colorrange = clims_b, colormap = :deep)

surface!(ax_b, y, z, b_slices.east; transformation = (:yz, x[end]), kwargs_b...)
surface!(ax_b, x, y, b_slices.bottom; transformation = (:xy, z[1]), kwargs_b...)
surface!(ax_b, x, y, b_slices.top; transformation = (:xy, z[end]), kwargs_b...)

b_avg = @lift interior(b_avg_timeseries[$n], 1, :, :)

surface!(ax_b, y, z, b_avg; transformation = (:yz, zonal_slice_displacement * x[end]),
         colorrange = clims_b,
         colormap = :deep)

contour!(ax_b, y, z, b_avg; levels = 15, transformation = (:yz, zonal_slice_displacement * x[end]),
         linewidth = 2,
         color = :black)

rotate_cam!(ax_b.scene, (π/20, -π/6, 0))

# Finally, we add a figure title with the time of the snapshot and then record a movie.

times = b_avg_timeseries.times

title = @lift "Buoyancy at t = " * string(round(times[$n] / day, digits=1)) * " days"

fig[1, 1] = Label(fig, title;
                  textsize = 24,
                  tellwidth = false,
                  padding = (0, 0, -120, 0))

frames = 1:length(times)

record(fig, filename * ".mp4", frames, framerate=8) do i
    @info "Plotting frame $i of $(frames[end])..."
    n[] = i
end
nothing #hide

# ![](baroclinic_adjustment.mp4)


text:
# Inter book operations
# TODO: Link on run time two books
# USAGE: You could have a private and a public book (both in different repositories)
# TODO: This might introduce cooperation capabilities

# Include system
# TODO: Redo include system (make it lazy)
# NOTES: Make it includes the whole book (even with the multiple passes maybe) only if currbook() === nothing
# NOTES: The rest of the time just include the modified files
# NOTES: The user can always force the whole loading at any time.
# TODO: redo the API

# Tags
# TODO: Include tag alias (maybe, not sure)
# TODO: Include tag types
# TODO: Include tag parameters
# TODO: Creates a tag recommendation system.
# TODO: Create a parametrized tag system. Like IL5.
# TODO: Make tags insensitive to space (non alpha numeric) and case. Create standardization tool (detect tag variants).

# Filetree
# TODO: Add folder sorting capabilities (book -> Readreport -> :year -> :author)
# TODO: Abstract the concept of a Structured Folder

# To Remember
# TODO: Create a 'To Remember' tag and a system for interacting with Anki (or any other to remember app)

# TODOS
# TODO: make a todo report

# Documments
# TODO: Revisit the doc stuctures (allow nested sections and sectionless documents)

module ResearchBooks

    import TOML
    import BibTeX
    import JSON
    import Random: randstring

    using ArgParse
    using Dates
    using OrderedCollections
    using EasyEvents

    # Re-exports
    using Reexport
    @reexport using FilesTreeTools
    @reexport using StringRepFilter

    include("Book/types.jl")
    include("Book/RBObject.jl")
    include("Book/RBDoc.jl")
    include("Book/RBSection.jl")
    include("Book/RBook.jl")
    include("Book/RBParagraph.jl")
    include("Book/RBQuote.jl")
    include("Book/RBNote.jl")
    include("Book/RBRef.jl")
    include("Book/RBTagMeta.jl")
    include("Book/RBRefList.jl")
    include("Book/utils.jl")

    export RBook, RBDoc, RBSection
    export RBObject
    export RBQuote, RBParagraph, RBNote
    export RBRef, RBRefList
    export RBTagMeta
    export refdict, getproperty
    export eachdoc, eachobj, eachtagmeta
    export tagmetas, objlabels
    export get_book, bookdir, get_label, localpath, localrelpath
    export get_doi, get_year, get_author, get_title, get_bibkey, get_tags
    export references

    include("FunctionalInterface/openbook.jl")
    include("FunctionalInterface/glob_state.jl")
    include("FunctionalInterface/getters.jl")
    include("FunctionalInterface/find.jl")
    # include("FunctionalInterface/markdown.jl")
    include("FunctionalInterface/replace_macros.jl")
    include("FunctionalInterface/genlabel.jl")
    include("FunctionalInterface/RBDoc_macros.jl")
    include("FunctionalInterface/RBSection_macros.jl")
    include("FunctionalInterface/RBTagMeta_macros.jl")
    include("FunctionalInterface/RBQuote_macros.jl")
    include("FunctionalInterface/RBNote_macros.jl")
    include("FunctionalInterface/RBParagraph_macros.jl")
    include("FunctionalInterface/newobj_macros.jl")
    include("FunctionalInterface/setmeta_macros.jl")
    include("FunctionalInterface/utils.jl")
    include("FunctionalInterface/show.jl")
    include("FunctionalInterface/filesys.jl")
    include("FunctionalInterface/tag_report.jl")

    export currdoc, currdoc!
    export currsec, currsec!
    export currobj, currobj!
    export currbook, currbook!
    export bookdir, srcline, srcfile
    export is_ipass0, is_ipass1, is_ipass2, get_ipass
    export genlabel, @genlabel!
    export openbook, @openbook
    export @new_document!, @new_section!
    export @new_paragraph!, @new_quote!, @new_note!
    export @setproperty!, @set_title!, @set_doi!, @set_text!
    export @set_author!, @set_year!, @set_bibkey!, @set_abstract!
    export @set_ctime!
    export @add_tag!, @add_todo!
    export @new_tagmeta!, @set_similars!, @set_antagonists!, @set_description!, @set_supertag!
    export show_tags
    export show_tag_report, @report
    # export @create_doc_readme

    export bookbib, findall_bookbib, findfirst_bookbib, filter_bookbib
    export crossrefs, findall_crossrefs, findfirst_crossrefs, filter_crossrefs
    export find_obj, filter_book

    include("Core/bookbib.jl")
    include("Core/crossref.jl")
    include("Core/git.jl")
    include("Core/open.jl")
    include("Core/configfile.jl")
    include("Core/rbfiles.jl")
    include("Core/utils.jl")

    include("ReadReport/dirtree.jl")
    include("ReadReport/create_read_report.jl")

    export create_read_report

    # include("Issues/dirtree.jl")
    # include("Issues/issues.jl")

end # end module


text:
# import Tmp # import module of slicerm
# commented out to use from REPL

using LinearAlgebra # to use dot fun

function g2(x, xa)
    return -0.5*dot(x,x)
end

res = Tmp.slicerm(g2, 2, [0.0 0.1], 0, N=500) # run a 500 sample MCMC

max_stars_repo_name: davidwcraig/slicer

max_issues_repo_name: davidwcraig/slicer

max_forks_repo_name: davidwcraig/slicer


text:
include("solver_preamble.jl")
import CSDP
factory = with_optimizer(CSDP.Optimizer, printlevel=0)
config = MOI.Test.TestConfig(atol=1e-4, rtol=1e-4, query=false)
@testset "Linear" begin
    Tests.linear_test(factory, config)
end
@testset "SDP" begin
    Tests.sd_test(factory, config)
end


text:
@doc raw"""
    solve_(prop:OptimalEstimation)

Rodgers-style linear optimal estimation, solves
```math
\begin{pmatrix} S_\epsilon^{-\frac{1}{2}}K \\ S_a^{-\frac{1}{2}} \end{pmatrix} x = \begin{pmatrix} S_\epsilon^{-\frac{1}{2}}y \\ 0 \end{pmatrix}
```
"""
function solve_stable(prob::OptimalEstimationProblem{FT}; errorAnalysis::Bool=true) where {FT}
    @unpack K, Sₑ, Sₐ, δy = prob
    n,m = size(K)
    # Compute S^{-0.5} using cholesky decomposition
    sqrtSₑ⁻¹ = inv(cholesky!(Sₑ).L);
    sqrtSₐ⁻¹ = inv(cholesky!(Sₐ).L);

    # A bit tediuous that way but reduces memory allocation:
    lhs = Array{FT}(undef,n+m,m)
    rhs = Array{FT}(undef,n+m)
    # K̃ = [sqrtSₑ⁻¹ * K; sqrtSₐ⁻¹]; ỹ = [sqrtSₑ⁻¹ * δy; 0]
    @views lhs[1:n,1:m]     = sqrtSₑ⁻¹ * K
    @views lhs[n+1:n+m,1:m] = sqrtSₐ⁻¹
    @views rhs[1:n]         = sqrtSₑ⁻¹ * δy
    @views rhs[n+1:end]    .= 0

    Kqr = qr!(lhs)

    # Solve like in Rodgers using QR decomposition, eq. 5.43,
    Δx̂ = Kqr \ rhs

    # Perform error analysis:
    if errorAnalysis
        Sₑ⁻¹ = sqrtSₑ⁻¹' * sqrtSₑ⁻¹
        # Compute Ŝ from K̃ = [sqrtSₑ⁻¹ * K; sqrtSₐ⁻¹]
        L = LinearAlgebra.inv!(UpperTriangular(Kqr.R));#UpperTriangular(Kqr.R)\I
        Ŝ = L * L'
        A = Ŝ * K' * Sₑ⁻¹ * K
        # residual r and reduced Χᵣ²
        r  = K*Δx̂ - δy
        # Reduced χ²
        χ² =  (r'* Sₑ⁻¹ *r)/(n-m)
        return OEsolution(Δx̂,Ŝ,A,r,χ²)
    end
    # Add xₐ here later?
    return Δx̂
end

@doc raw"""
    solve(prop:OptimalEstimation)

Rodgers-style linear optimal estimation, solves
\hat{x} = x_a + (K^TS_\epsilon^{-1} K + S_a^{-1})^{-1} K^T S_\epsilon^{-1}\delta y
"""
function solve(prob::OptimalEstimationProblem; errorAnalysis::Bool=true)
    @unpack K, Sₑ, Sₐ, δy = prob
    n,m = size(K)
    # Inverse(s)
    Sₑ⁻¹ = inv(Sₑ)
    Sₐ⁻¹ = inv(Sₐ)
    # Will be re-used:
    KᵀSₑ⁻¹K = K' * Sₑ⁻¹ * K

    # Posterior covariance matrix:
    Ŝ = inv(KᵀSₑ⁻¹K + Sₐ⁻¹)

    # Estimated state vector (eq 4.5):
    Δx̂ = Ŝ * K'Sₑ⁻¹ * δy
    if errorAnalysis
        # Averaging kernel matrix:
        A = Ŝ * KᵀSₑ⁻¹K
        # residual r and reduced Χᵣ²
        r  = K*Δx̂ - δy
        # Reduced χ²
        χ² =  (r'*Sₑ⁻¹*r)/(n-m)
        return OEsolution(Δx̂,Ŝ,A,r,χ²)
    end
    return Δx̂
end

# sqrtSₑ⁻¹ = sqrt.(inv(Sₑ));
# Solve like in Rodgers, eq. 5.43, supposedly more stable
#    lh = [sqrtSₑ⁻¹ * K; Sₐ]
#    rh = [sqrtSₑ⁻¹ * δy; zeros(nx)]
#    return lh \ rh

max_stars_repo_name: RadiativeTransfer/OptimalEstimation

max_issues_repo_name: RadiativeTransfer/OptimalEstimation

max_forks_repo_name: RadiativeTransfer/OptimalEstimation


text:
module TestExamples2DUnstructured

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "unstructured_2d_dgsem")

@testset "Unstructured Mesh for Euler" begin

  @trixi_testset "elixir_euler_periodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_periodic.jl"),
      l2   = [0.00010978828464875207, 0.00013010359527356914, 0.00013010359527326057, 0.0002987656724828824],
      linf = [0.00638626102818618, 0.009804042508242183, 0.009804042508253286, 0.02183139311614468])
  end

  @trixi_testset "elixir_euler_free_stream.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"),
      l2   = [3.357431396258123e-14, 2.439943089578555e-13, 1.4655386790023588e-13, 4.670410488845425e-13],
      linf = [1.0169198816356584e-11, 6.838458965763294e-11, 4.400946274074613e-11, 1.4071055431941204e-10],
      tspan = (0.0, 0.1),
      atol = 3.0e-13)
  end

  @trixi_testset "elixir_euler_wall_bc.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_wall_bc.jl"),
      l2   = [0.0401951408663404, 0.042562446022296446, 0.03734280979760256, 0.10058806400957201],
      linf = [0.24502910304890335, 0.298141188019722, 0.29446571031375074, 0.5937151600027155],
      tspan = (0.0, 0.25))
  end

  @trixi_testset "elixir_euler_basic.jl" begin
    @test_trixi_include(default_example_unstructured(),
      l2   = [0.0007258658867098887, 0.000676268065087451, 0.0006316238024054346, 0.0014729738442086392],
      linf = [0.004476908674416524, 0.0052614635050272085, 0.004926298866533951, 0.018058026023565432],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_euler_restart.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_restart.jl"),
      l2   = [0.0007258658867098887, 0.000676268065087451, 0.0006316238024054346, 0.0014729738442086392],
      linf = [0.004476908674416524, 0.0052614635050272085, 0.004926298866533951, 0.018058026023565432])
  end

  @trixi_testset "elixir_euler_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.06594600495903137, 0.10803914821786433, 0.10805946357846291, 0.1738171782368222],
      linf = [0.31880214280781305, 0.3468488554333352, 0.34592958184413264, 0.784555926860546],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_advection_basic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"),
      l2   = [0.00023219572238346008],
      linf = [0.0017401556568237275])
  end

  @trixi_testset "elixir_ape_gauss_wall.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_ape_gauss_wall.jl"),
      l2   = [0.029331247985489625, 0.02934616721732521, 0.03803253571320854, 0.0,
              7.465932985352019e-16, 1.4931865970704038e-15, 1.4931865970704038e-15],
      linf = [0.3626825396196784, 0.3684490307932018, 0.8477478712580901, 0.0,
              8.881784197001252e-16, 1.7763568394002505e-15, 1.7763568394002505e-15],
      tspan = (0.0, 5.0))
  end

  @trixi_testset "elixir_mhd_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
      l2   = [0.06418293357851637, 0.12085176618704108, 0.12085099342419513, 0.07743005602933221,
              0.1622218916638482, 0.04044434425257972, 0.04044440614962498, 0.05735896706356321,
              0.0020992340041681734],
      linf = [0.1417000509328017, 0.3210578460652491, 0.335041095545175, 0.22500796423572675,
              0.44230628074326406, 0.16743171716317784, 0.16745989278866702, 0.17700588224362557,
              0.02692320090677309],
      tspan = (0.0, 0.5))
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [5.377518922553881e-5, 0.09999999206243514, 0.09999999206243441, 0.1414213538550799,
              8.770450430886394e-6, 0.0999999926130084, 0.0999999926130088, 0.14142135396487032,
              1.1553833987291942e-5],
      linf = [0.00039334982566352483, 0.14144904937275282, 0.14144904937277897, 0.20003315928443416,
              6.826863293230012e-5, 0.14146512909995967, 0.14146512909994702, 0.20006706837452526,
              0.00013645610312810813],
      tspan = (0.0, 0.5))
  end
end

end # module


text:
using TensorCalculus
using Documenter

DocMeta.setdocmeta!(TensorCalculus, :DocTestSetup, :(using TensorCalculus); recursive=true)

makedocs(;
    modules=[TensorCalculus],
    authors="Elias Kempf <eli2008@gmx.de> and contributors",
    repo="https://github.com/icetube23/TensorCalculus.jl/blob/{commit}{path}#{line}",
    sitename="TensorCalculus.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://icetube23.github.io/TensorCalculus.jl",
        assets=String[],
        edit_link="main",
    ),
    pages=["Home" => "index.md"],
    strict=:doctest,
)

deploydocs(; repo="github.com/icetube23/TensorCalculus.jl", devbranch="main")


text:
using KernelDensity

"""
    initByMsmt(data, numSegments, numMeas, indMeas, selectMeas, selectSign, threshold, docond)

Condition data on measurement results on the same channel

numSegments = number of distinct segments in the sequence
numMeas = number of meauserements per segments
indMeas = list of indeces of the msmt's to compare for fidelity
selectMeas = indeces of msms't to postselect on. All the msm'ts in a
segment will be postselected on this one
selectSign = select the sign for postselection. 1 = keep greater than;
-1 = keep smaller than
threshold
docond. If 1, do postselection
"""
function initByMsmt(data,numSegments,numMeas,indMeas,selectMeas,selectSign,threshold,docond)

  ind0=indMeas[1,:]
  ind1=indMeas[2,:]
  numShots = Int64(length(data)/(numMeas*numSegments))
  data=real(data)'
  datamat = splitdata(data,numShots,numMeas*numSegments)
  bins = range(minimum(minimum(datamat)), stop=maximum(maximum(datamat)),length=500)

  PDFvec = zeros(length(bins),numMeas*numSegments)
  PDFvec_con = zeros(length(bins),numMeas*numSegments)

  for kk=1:numMeas*numSegments
    PDFvec[:,kk] = kde(datamat[:,kk], bins).density
  end

  opt_thr = zeros(length(ind0),1)
  fidvec_un = zeros(length(ind0),1)
  fidvec_con = zeros(length(ind0),1)
  err0 = zeros(length(ind0),1)
  err1 = zeros(length(ind0),1)
  data_con = zeros(numSegments, numMeas-length(selectMeas))

  for kk=1:size(PDFvec,2)
    PDFvec[:,kk]=PDFvec[:,kk]/sum(PDFvec[:,kk])
    PDFvec[:,kk]=PDFvec[:,kk]/sum(PDFvec[:,kk])
  end

  for kk=1:length(ind0)
    fidvec_un[kk] = 1-0.5*(1-0.5*sum(abs.(PDFvec[:,ind0[kk]]-PDFvec[:,ind1[kk]])))
    indmaximum = argmax(abs.(cumsum(PDFvec[:,ind0[kk]])-cumsum(PDFvec[:,ind1[kk]])))
    @printf("Optimum unconditioned fid. for segm. %d and %d = %.3f\n", ind0[kk], ind1[kk], fidvec_un[kk])
    tempvec0 = cumsum(abs.(PDFvec[:,ind0[kk]]))
    err0[kk] = tempvec0[indmaximum]
    @printf("Error for |0> for segm. %d and %d = %.3f\n", ind0[kk], ind1[kk], tempvec0[indmaximum])
    tempvec1 = 1-cumsum(abs.(PDFvec[:,ind1[kk]]))
    err1[kk] = tempvec1[indmaximum]
    @printf("Error for |1> for segm. %d and %d = %.3f\n", ind0[kk], ind1[kk], tempvec1[indmaximum])
    opt_thr[kk] = bins[indmaximum]
    @printf("Optimum threshold for segments %d and %d = %.4f\n", ind0[kk], ind1[kk], opt_thr[kk])
  end

  if(docond>0)
    for mm=1:numSegments
      for kk=1:length(selectMeas)
        ind = selectMeas[kk]
        for jj=1:numShots
          if(mm>4)
          end
          if selectSign == 1 && datamat[jj,numMeas*(mm-1)+ind] < threshold
            datamat[jj,numMeas*(mm-1)+1:numMeas*(mm-1)+numMeas]=NaN
          elseif selectSign == -1 && datamat[jj,numMeas*(mm-1)+ind] > threshold
            datamat[jj,numMeas*(mm-1)+1:numMeas*(mm-1)+numMeas]=NaN
          end
        end
      end
    end
  end

  for jj=1:numSegments
    thismeas=1;
    for kk=1:numMeas
      tempvec = filter!(vec->~isnan(vec) ,datamat[:,numMeas*(jj-1)+kk])
      PDFvec_con[:,(jj-1)*numMeas+kk] = kde(tempvec, bins).density

      #hist(~isnan(datamat[:,numMeas*(jj-1)+kk]).*datamat[:,numMeas*(jj-1)+kk], bins)[2];
      if(size(findall(selectMeas.==kk),1)==0)
        data_con[jj,thismeas] = nanmean(datamat[:,numMeas*(jj-1)+kk])
        thismeas=thismeas+1
      end
    end
  end

  #normalize
  for kk=1:size(PDFvec_con,2)
    PDFvec_con[:,kk]=PDFvec_con[:,kk]/sum(PDFvec_con[:,kk])
    PDFvec_con[:,kk]=PDFvec_con[:,kk]/sum(PDFvec_con[:,kk])
  end

  for kk=1:length(ind0)
    fidvec_con[kk] = 1-0.5*(1-0.5*sum(abs.(PDFvec_con[:,ind0[kk]]-PDFvec_con[:,ind1[kk]])))
    @printf("Optimum conditioned fid. for segm. %d and %d = %.3f\n", ind0[kk], ind1[kk], fidvec_con[kk])
  end
  return bins, PDFvec, PDFvec_con,data_con,opt_thr,datamat, fidvec_un, err0, err1, fidvec_con
end

function initByMsmt_2D(data,Anum::Vector{Int},numSegments,numMeas_A,indMeas,selectMeas,selectSign,docond,numCal;threshold::Vector{Float64}=NaN)
  #data = data qubits
  #data = single ancilla qubit. Data to be postselected on
  #Anum = index of ancilla measurement
  #numSegments = number of distinct segments in the sequence
  #numMeas_A = number of ancilla measurements per segments per channel (including those to be conditioned on)
  #indMeas = list of indeces of the msmt's to compare for fidelity. 1st dim = segments; 2nd dim: A channel. Note: removed comparison between multiple segments
  #selectMeas = indeces of msms't to postselect on. All the msm'ts in a
  #segment will be postselected on this one
  #selectSign = select the sign for postselection. 1 = keep greater than;
  #-1 = keep smaller than
  #docond. If 1, do postselection
  #numCal = number of calibration sequences. Currently, it assumes that they only have 1
  #measurement per qubit.
  #threshold. If not specified, uses the optimum to maximize fidelity

  numMeas_D=numMeas_A
  #the number of data measurements can be different from numMeas_A. For
  #instance, there may be multiple conditioning rounds before data
  #tomography. However, there is currently the constraint that all channels,
  #even if in different cards, acquire the same number of shots.

  ind0=indMeas[1,:]
  ind1=indMeas[2,:]
  nbins = 500

  data_A = zeros(length(data[1]), length(Anum))
  for (k, Aind) in enumerate(Anum)
    data_A[:, k] = real(data[Aind])
  end
  num_data_meas = length(data)-length(Anum)
  data_D = Dict()
  datamat = Dict()
  var_con = Dict()
  data_con = Dict()

  Dnum = setdiff(1:length(data), Anum)
  [data_D[kk] = data[Dnum[kk]] for kk = 1:num_data_meas]

  numShots_A = Int64(floor(length(data_A[:, 1])/(numMeas_A*(numSegments-numCal)+numCal))) #assume same number of meas for all As
  for kk = 1:num_data_meas
    datamat[kk] = splitdata(real(data_D[kk]),numShots_A,numMeas_D*(numSegments-numCal)+numCal)
  end

  datamat_A = zeros(numShots_A, numMeas_A*(numSegments-numCal)+numCal, length(Anum))
  bins = Array(Base.LinRange{Float64},length(Anum))
  PDFvec = zeros(nbins,numMeas_A*(numSegments-numCal)+numCal, length(Anum))

  for Aind = 1:length(Anum)

    datamat_A[:, :, Aind] = splitdata(data_A[:, Aind],numShots_A,numMeas_A*(numSegments-numCal)+numCal)
    bins[Aind] = range(minimum(minimum(datamat_A[:, :, Aind])), stop=maximum(maximum(datamat_A[:, :, Aind])),length=nbins)

    #PDFvec_con = zeros(length(bins),numMeas*numSegments);

    for kk=1:numMeas_A*(numSegments-numCal)+numCal
      PDFvec[:,kk,Aind] = kde(datamat_A[:,kk,Aind], bins[Aind]).density
    end
  end

  opt_thr = zeros(size(ind0))
  fidvec_un = zeros(size(ind0))
  #fidvec_con = zeros(length(ind0),1);
  err0 = zeros(size(ind0))
  err1 = zeros(size(ind0))
  datatemp = zeros(numSegments, numMeas_D)
  vartemp = zeros(numSegments, numMeas_D)
  #renormalize
  for Aind=1:length(Anum)
    for kk=1:size(PDFvec,2)
      PDFvec[:,kk,Aind]=PDFvec[:,kk,Aind]/sum(PDFvec[:,kk,Aind])
      PDFvec[:,kk,Aind]=PDFvec[:,kk,Aind]/sum(PDFvec[:,kk,Aind])
    end
  end

  for Aind = 1:length(Anum)
    @printf("Ancilla n.%d\n", Aind)
    #for kk=1:size(ind0,1)
    fidvec_un[Aind] = 1-0.5*(1-0.5*sum(abs.(PDFvec[:,ind0[Aind],Aind]-PDFvec[:,ind1[Aind],Aind])))
    indmaximum = argmax(abs.(cumsum(PDFvec[:,ind0[Aind],Aind])-cumsum(PDFvec[:,ind1[Aind],Aind])))
    @printf("Optimum unconditioned fid. for segm. %d and %d = %.3f\n", ind0[Aind], ind1[Aind], fidvec_un[Aind])
    tempvec0 = cumsum(abs.(PDFvec[:,ind0[Aind],Aind]))
    err0[Aind] = tempvec0[indmaximum]
    @printf("Error for |0> for segm. %d and %d = %.3f\n", ind0[Aind], ind1[Aind], tempvec0[indmaximum])
    tempvec1 = 1-cumsum(abs.(PDFvec[:,ind1[Aind],Aind]))
    err1[Aind] = tempvec1[indmaximum]
    @printf("Error for |1> for segm. %d and %d = %.3f\n", ind0[Aind], ind1[Aind], tempvec1[indmaximum])
    opt_thr[Aind] = bins[Aind][indmaximum]
    @printf("Optimum threshold for segments %d and %d = %.4f\n", ind0[Aind], ind1[Aind], opt_thr[Aind])
  end

  if docond>0
    if any(isnan(threshold))
      threshold  = opt_thr #uses one of the fidelity-maximizing values if unspecified
    end
    @printf("Length =%d\n", length(data_D))
    dataslice = datamat[1]
    for nn = 1:length(data_D)
      dataslice = datamat[nn]
      for Aind = 1:length(Anum)
        for mm=1:numSegments-numCal
          for kk=1:length(selectMeas)
            ind = selectMeas[kk]

            for jj=1:Int(numShots_A*numMeas_A/numMeas_D)
              if selectSign[kk, Aind] == 1 && datamat_A[jj,numMeas_A*(mm-1)+ind,Aind] < threshold[Aind]
                dataslice[jj,numMeas_A*(mm-1)+1:numMeas_A*(mm-1)+numMeas_A]=NaN
              elseif selectSign[kk, Aind] == -1 && datamat_A[jj,numMeas_A*(mm-1)+ind,Aind] > threshold[Aind]
                dataslice[jj,numMeas_A*(mm-1)+1:numMeas_A*(mm-1)+numMeas_A]=NaN
              end
            end

          end
        end
        datamat[nn] = dataslice
        frac = sum(sum(~isnan(dataslice[:,1:end-numCal])))/(size(dataslice,1)*(size(dataslice,2)-numCal))
        @printf("Fraction kept A%d = %.2f\n", Aind, frac)
      end

    end
  end

  for nn=1:length(data_D)
    dataslice = datamat[nn]
    for jj=1:numSegments
      thismeas=1
      for kk=1:numMeas_D
        if jj<numSegments-numCal+1
          datatemp[jj, thismeas] = nanmean(dataslice[:,numMeas_D*(jj-1)+kk])
          vartemp[jj, thismeas] = nanvar(dataslice[:, numMeas_D*(jj-1)+kk])
        else
          datatemp[jj, thismeas] = nanmean(dataslice[:,numMeas_D*(numSegments-numCal) + (jj-(numSegments-numCal))])
          vartemp[jj, thismeas] = nanvar(dataslice[:,numMeas_D*(numSegments-numCal) + (jj-(numSegments-numCal))]) #same cal. points for all data measurements
        end
        thismeas=thismeas+1
      end
    end

    data_con[nn] = datatemp[:]
    var_con[nn] = vartemp[:]

  end

  return (bins,PDFvec,data_con,var_con,opt_thr,datamat, fidvec_un, err0, err1)
end

function initByMsmt_2D(data,Anum::Int,numSegments,numMeas_A,indMeas,selectMeas,selectSign,docond,numCal;threshold::Float64=NaN)
  return initByMsmt_2D(data,[Anum],numSegments,numMeas_A,indMeas,selectMeas,selectSign,docond,numCal;threshold = fill(threshold,length(Anum)))
end

function splitdata(data, numseg, numsets)
  datamat = zeros(numseg, numsets);
  for kk=1:numsets
    datamat[:,kk] = data[kk:numsets:end];
  end
  return datamat
end


text:
using TriangularMatrices
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
using StaticArrays
using Compat, Compat.LinearAlgebra

x = @SMatrix randn(5,3)
y = @SMatrix randn(3,5)
SAS = x' * x
SAS2 = y * y'

TMS = xtx(x)
TMS2 = xxt(y)

# @show TMS
# @show SAS
# @show TMS2
# @show SAS2

@test all(TMS .== SAS)
@test all(TMS2 .== SAS2)

function check_approx_equality(A, R)
    @assert size(A) == size(R)
    m, n = size(A)
    all_true = true
    for i = 1:n, j = 1:m
        all_true = A[j,i] ≈ R[j,i]
        all_true || ( @show((j,i)); break )
    end
    all_true
end

# @time @testset for M = 1:128, N = round(Int, M/√2, RoundUp):round(Int, M*√2, RoundDown), P = max(round(Int, M/√2, RoundUp), round(Int, N/√2, RoundUp)):min(round(Int, M*√2, RoundDown),round(Int, N*√2, RoundDown))

@time @testset for M = 20:128, N = round(Int, M/√2, RoundUp):round(Int, M*√2, RoundDown), P = max(round(Int, M/√2, RoundUp), round(Int, N/√2, RoundUp)):min(round(Int, M*√2, RoundDown),round(Int, N*√2, RoundDown))
        d_r = randmat(M,P)
        a_r = randmat(M,N)
        x_r = randmat(N,P)
        TriangularMatrices.mul!(d_r, a_r, x_r)
        d_m = Matrix(d_r)
        if VERSION < v"0.7-"
            A_mul_B!(d_m, Matrix(a_r), Matrix(x_r))
        else
            mul!(d_m, Matrix(a_r), Matrix(x_r))
        end
        @test check_approx_equality(d_m, d_r)
    # end
end


text:
module PkgTestRoot

println("EVALUATING PkgTestRoot!")

fnX() = println("PkgTestRoot.fnX()")

function __init__()
	println("Initializing PkgTestRoot!")
end

end # module


text:
@cache mutable struct AitkenNevilleCache{uType,rateType,arrayType,dtType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  utilde::uType
  atmp::uNoUnitsType
  fsalfirst::rateType
  dtpropose::dtType
  T::arrayType
  cur_order::Int
  work::dtType
  A::Int
  step_no::Int
end

@cache mutable struct AitkenNevilleConstantCache{dtType,arrayType} <: OrdinaryDiffEqConstantCache
  dtpropose::dtType
  T::arrayType
  cur_order::Int
  work::dtType
  A::Int
  step_no::Int
end

function alg_cache(alg::AitkenNeville,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  utilde = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  cur_order = max(alg.init_order, alg.min_order)
  dtpropose = zero(dt)
  T = fill(zeros(eltype(u), size(u)), (alg.max_order, alg.max_order))
  work = zero(dt)
  A = one(Int)
  atmp = similar(u,uEltypeNoUnits)
  step_no = zero(Int)
  AitkenNevilleCache(u,uprev,tmp,k,utilde,atmp,fsalfirst,dtpropose,T,cur_order,work,A,step_no)
end

function alg_cache(alg::AitkenNeville,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  dtpropose = zero(dt)
  cur_order = max(alg.init_order, alg.min_order)
  T = fill(zero(eltype(u)), (alg.max_order, alg.max_order))
  work = zero(dt)
  A = one(Int)
  step_no = zero(Int)
  AitkenNevilleConstantCache(dtpropose,T,cur_order,work,A,step_no)
end


text:
using Test, TestSetExtensions, SafeTestsets, Logging, Pkg
using HTTP, Genie, Pluto

include("html.jl")

function testnotebook(input)
    session = Pluto.ServerSession();
    notebook = Pluto.SessionActions.open(session, input; run_async=false)
    html_contents = body(notebook)

    # split input string by '--', '.' and '_' to array of strings
    # e.g. '18--a_b_c.jl' -> ['18', 'a', 'b', 'c', 'jl']
    # html_input string
    first_input = last(split(input, "/"))
    new_input = replace(last(split(split(first_input, ".")[1], "--")), r"_" => s"-")

    write(joinpath(@__DIR__, "../build/guides/","$(new_input).html"), html_contents)

    errored=false
    for c in notebook.cells
        if c.errored
            errored=true
            @error "Error in  $(c.cell_id): $(c.output.body[:msg])\n $(c.code)"
        end
        @assert !c.errored
        @test !c.errored
    end
    !errored
end

@testset "testing guides" begin
    @info @__DIR__
    notebooks= [
    "0--Index/1--index.jl",
    "183--Simple_API_backend/183--Simple_API_backend.jl",
    "182--Genie_Plugins/182--Genie_Plugins.jl",
    "184--Interactive_environment/184--Interactive_environment.jl"
    ]
    for notebook in notebooks
        input=joinpath(@__DIR__,"PlutoNotebooks/",notebook)
        @show input
        @info "notebook: $(input)"
        @test testnotebook(input)
    end
end


text:
# ---
# title: 1282. Group the People Given the Group Size They Belong To
# id: problem1282
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/group-the-people-given-the-group-size-they-belong-to/description/>
# hidden: true
# ---
#
# There are `n` people that are split into some unknown number of groups. Each
# person is labeled with a  **unique ID**  from `0` to `n - 1`.
#
# You are given an integer array `groupSizes`, where `groupSizes[i]` is the size
# of the group that person `i` is in. For example, if `groupSizes[1] = 3`, then
# person `1` must be in a group of size `3`.
#
# Return  _a list of groups  such that each person `i` is in a group of size
# `groupSizes[i]`_.
#
# Each person should appear in  **exactly one group** , and every person must be
# in a group. If there are multiple answers, **return any of them**. It is
# **guaranteed** that there will be **at least one** valid solution for the
# given input.
#
#
#
# **Example 1:**
#
#
#
#     Input: groupSizes = [3,3,3,3,3,1,3]
#     Output: [[5],[0,1,2],[3,4,6]]
#     Explanation:
#     The first group is [5]. The size is 1, and groupSizes[5] = 1.
#     The second group is [0,1,2]. The size is 3, and groupSizes[0] = groupSizes[1] = groupSizes[2] = 3.
#     The third group is [3,4,6]. The size is 3, and groupSizes[3] = groupSizes[4] = groupSizes[6] = 3.
#     Other possible solutions are [[2,1,6],[5],[0,4,3]] and [[5],[0,6,2],[4,3,1]].
#
#
# **Example 2:**
#
#
#
#     Input: groupSizes = [2,1,3,3,3,2]
#     Output: [[1],[0,5],[2,3,4]]
#
#
#
#
# **Constraints:**
#
#   * `groupSizes.length == n`
#   * `1 <= n <= 500`
#   * `1 <= groupSizes[i] <= n`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
import QuadTreeMeshes
using Base.Test
import GeometryTypes

plot = true

if plot==true
    import Plots
    Plots.plotlyjs()
end

include("TestHelpers.jl")
include("CurrentTest.jl")

#include("MeshTests.jl")
#include("QuadTreeTests.jl")


text: # Place here configuration options that will be set for all environments

max_stars_repo_name: Fatalon/Dashboard


text:
abstract type AbstractViscosity end
abstract type Isotropic end
abstract type Anisotropic end

struct Isoviscous{T} <: AbstractViscosity
    val::Float64
end

struct TemperatureDependant{T} <: AbstractViscosity
    val::Vector{Float64}
end

struct IsoviscousPlastic{T} <: AbstractViscosity
    node::Float64
    ip::Matrix{Float64}
    τ_VonMises::Float64

    function IsoviscousPlastic(T, val, τ_VonMises, nel, nip)
        new{T}(val, zeros(nel, nip), τ_VonMises)
    end
end

struct TemperatureDependantPlastic{T} <: AbstractViscosity
    node::Vector{Float64}
    ip::Matrix{Float64}
    τ_VonMises::Float64

    function TemperatureDependantPlastic(T, val, τ_VonMises, nel, nip)
        new{T}(val, zeros(nel, nip), τ_VonMises)
    end
end

struct Mallard2016{T} <: AbstractViscosity
    node::Vector{Float64}
    ip::Matrix{Float64}
    τ_VonMises::Float64

    function TemperatureDependantPlastic(T, val, τ_VonMises, nel, nip)
        new{T}(val, zeros(nel, nip), τ_VonMises)
    end
end

Base.Val(x::Isoviscous{Isotropic})             = Val{Isotropic}()
Base.Val(x::Isoviscous{Anisotropic})           = Val{Anisotropic}()
Base.Val(x::TemperatureDependant{Isotropic})   = Val{Isotropic}()
Base.Val(x::TemperatureDependant{Anisotropic}) = Val{Anisotropic}()
Base.Val(x::TemperatureDependantPlastic{Isotropic}) = Val{Isotropic}()

function getviscosity(T, type, nel; η = 1, τ_VonMises = 5e3, nip = 7)

    if type == :IsoviscousIsotropic
        return Isoviscous{Isotropic}(float(η))

    elseif type == :IsoviscousAnisotropic
        return Isoviscous{Anisotropic}(float(η))

    elseif type == :TemperatureDependantIsotropic
        return TemperatureDependant{Isotropic}(@. exp(13.8156*(0.5-T)))

    elseif type == :TemperatureDependantIsotropicPlastic
        # return TemperatureDependantPlastic(Isotropic, @.(exp(13.8156*(0.5-T))), τ_VonMises, nel, nip)
        return TemperatureDependantPlastic(Isotropic, @.(exp(23.03/(1+T) -23.03/2)), τ_VonMises, nel, nip)

    elseif type == :VanHeckIsotropic
        return TemperatureDependant{Isotropic}(@. exp(23.03/(1+T) -23.03/2))

    elseif type == :TemperatureDependantAnisotropic
        return TemperatureDependant{Anisotropic}(@. exp(13.8156*(0.5-T)))

    elseif type == :TemperatureDependantAnisotropicPlastic
        # return TemperatureDependant(Anisotropic, @.(exp(13.8156*(0.5-T))), τ_VonMises, nel, nip)
        return TemperatureDependant(Anisotropic, @.(exp(23.03/(1+T) -23.03/2)), τ_VonMises, nel, nip)

    elseif type == :VanHeckAnisotropic
        return TemperatureDependant{Anisotropic}(@. exp(23.03/(1+T) -23.03/2))

    end

end

getviscosity!(η::Isoviscous, T) = nothing

function getviscosity!(η::TemperatureDependant, T)
    Threads.@threads for i in eachindex(η.val)
        # update viscosity
        @inbounds η.val[i] = exp(13.8156*(0.5-T[i]))
    end
end

function getviscosity!(η::TemperatureDependantPlastic, T, r)
    @inbounds for i in eachindex(η.node)
        depth = r[i] - 1.22
        # viscosity correction from Richards et al. [2001] and Tackley [2000b]
        m = T[i] < 0.6 + 2*(1-depth) ? 1 : 0.1
        # update viscosity
        η.node[i] = exp(23.03/(1+T[i]) -23.03/2) * m
    end
end

function getviscosity!(η::Mallard2016, T, r)

    # Parameters from Mallard et al 2016 - Nature
    a, B, d0, dstep = 1e6, 30, 0.276, 0.02

    # η(z,T) from Mallard et al 2016 - Nature
    @inbounds for i in eachindex(η)
        depth = 2.22 - r[i]
        # viscosity correction
        m = T[i] < 0.6 + 7.5*depth ? 1 : 0.1
        # depth dependent component
        ηz = a*exp(log(B)*(1-0.5*(1-tanh((d0-depth)/dstep))))
        # update viscosity
        η[i] = ηz*exp(0.064 - 30/(T[i]+1))*m
    end
end

# Temperature effect on density
# (α is thermal expansion coeff, ρ is defined at the nodes)
state_equation(α, T; ρ0 = 1) = @. ρ0 * (1-α*T);

state_equation!(ρ,α,T; ρ0 = 1) = ρ .= ones(length(T))

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
# Lar = LinearAlgebraicRepresentation

function interior_to_f(triangle,f,V,FV,EV,FE)
	v1,v2,v3 = triangle
	u = V[:,v2]-V[:,v1]
	v = V[:,v3]-V[:,v1]
	w = cross(u,v)
	T = [1. 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]; T[1:3,4] = -V[:,v1]
	R = [1. 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]; R[1:3,1:3] = [u v w]; R = R'
	mapping = R * T

	trianglepoints = [[V[:,v1] V[:,v2] V[:,v3]]; ones(3)']
	pts = mapping * trianglepoints
	points2D = [pts[r,:] for r = 1:size(pts,1)
				if !(all(pts[r,:].==0) || all(pts[r,:].==1) )]
	p2D = hcat(points2D...)'
	checkpoint = 0.4995 .* p2D[:,1] + 0.4995 .* p2D[:,2] + 0.001 .* p2D[:,3]

	cellpoints = [V[:,FV[f]]; ones(length(FV[f]))' ]
	points = mapping * cellpoints
	verts2D = [points[r,:] for r = 1:size(points,1)
				if !(all(points[r,:].==0) || all(points[r,:].==1) )]
	P2D = hcat(verts2D...)'

	vdict = Dict(collect(zip(FV[f], 1:length(FV[f]))))
	celledges = [[vdict[v] for v in EV[e]] for e in FE[f]]

	out = pointInPolygonClassification(P2D,celledges)(checkpoint)
	if out=="p_in"
		return true
	else
		return false
	end
end

function ordering(triangles,V)
	normals = []
	v1,v2,v3 = triangles[1]
	if v1>v2 v1,v2 = v2,v1 end
	e3 = normalize(V[:,v2]-V[:,v1])
	e1 = normalize(V[:,v3]-V[:,v1])
	e2 = normalize(cross(e1,e3))
	basis = [e1 e2 e3]
	transform = inv(basis)

	angles = []
	for (v1,v2,v3) in triangles
		w1 = normalize(V[:,v3]-V[:,v1])
		w2 = transform * w1
		w3 = cross([0,0,1],w2)
		push!(normals,w3)
	end
	for k=1:length(normals)
		angle = atan(normals[k][2],normals[k][1])
		push!(angles,angle)
	end
	#pairs = sort(collect(zip(angles,1:length(triangles))))
	pairs = sort([(angle,k) for (k,angle) in enumerate(angles)])
	order = [k for (angle,k) in pairs]
	return order
end

function ord(hinge::Int64, bd1::AbstractSparseVector{Int64,Int64}, V::Array{Float64,2},
FV::Array{Array{Int64,1},1}, EV::Array{Array{Int64,1},1}, FE::Array{Array{Int64,1},1})
	cells = SparseArrays.findnz(bd1)[1]
	triangles = []

	function area(v1,v2,v3)
		u = V[:,v2]-V[:,v1]
		v = V[:,v3]-V[:,v1]
		out = norm(cross(u,v)) # actually, to be divided by two
		return out
	end

	for f in cells
		v1,v2 = EV[hinge]
		index = findfirst(v -> (area(v1,v2,v)≠0), FV[f])
		v3 = FV[f][index]

		# test if [v1,v2,v3] interior to f
		while true
			if interior_to_f([v1, v2, v3],f,V,FV,EV,FE)
				push!(triangles, [v1,v2,v3])
				break
			else
				index = findnext(v -> (area(v1,v2,v)≠0), FV[f], index+1)
				v3 = FV[f][index]
			end
		end
	end
	order = ordering(triangles,V)
	return [cells[index] for index in order]
end

function mynext(cycle, pivot)
	len = length(cycle)
	ind = findall(x -> x==pivot, cycle)[1]
	nextIndex = ind==len ? 1 : ind+1
	return cycle[nextIndex][1]
end

function myprev(cycle, pivot)
	len = length(cycle)
	ind = findall(x->x==pivot, cycle)[1]
	nextIndex = ind==1 ? len : ind-1
	return cycle[nextIndex][1]
end

function build_copFC(rV, rcopEV, rcopFE)
#function build_copFC(V,FV,EV,copFE)

	# G&F -> Pao data structures
	V = convert(Lar.Points, rV')
	EV = cop2lar(rcopEV)
	fe = cop2lar(rcopFE)
	fv = [collect(Set(cat(EV[e] for e in fe[f]))) for f=1:length(fe)]
	FV = convert(Lar.Cells, fv)
	copFE = rcopFE

	copEF = copFE'
	FE = [SparseArrays.findnz(copFE[k,:])[1] for k=1:size(copFE,1)]
	# Initializations
	m,n = size(copEF)
	marks = zeros(Int,n);
	I = Int64[]; J = Int64[]; W = Int8[];
	jcol = 0
	choose(marks) = findfirst(x -> x<2, marks)

	# Main loop (adding one copFC's column stepwise)
	while sum(marks) < 2n
		# select a (d−1)-cell, "seed" of the column extraction
		σ = choose(marks)
		if marks[σ] == 0
			cd1 = sparsevec([σ], Int8[1], n)
		elseif marks[σ] == 1
			cd1 = sparsevec([σ], Int8[-1], n)
		end
		# compute boundary cd2 of seed cell
		cd2 = copEF * cd1
		# loop until (boundary) cd2 becomes empty
		while nnz(cd2)≠0
			corolla = sparsevec([], Int8[], n)
			# for each “hinge” τ cell
			for τ ∈ (.*)(SparseArrays.findnz(cd2)...)
				#compute the  coboundary
				tau = sparsevec([abs(τ)], [sign(τ)], m)
				bd1 = transpose(transpose(tau) * copEF)
				cells2D = SparseArrays.findnz(bd1)[1]
				# compute the  support
				pivot = intersect(cells2D, SparseArrays.findnz(cd1)[1])[1]
				# compute the new adj cell
				fan = ord(abs(τ),bd1,V,FV,EV,FE) # ord(pivot,bd1)
				if τ > 0
					adj = mynext(fan,pivot)
				elseif τ < 0
					adj = myprev(fan,pivot)
				end
				# orient adj
				if copEF[abs(τ),adj] ≠ copEF[abs(τ),pivot]
					corolla[adj] = cd1[pivot]
				else
					corolla[adj] = -(cd1[pivot])
				end
			end
			# insert corolla cells in current cd1
			for (k,val) in zip(SparseArrays.findnz(corolla)...)
				cd1[k] = val
			end
			# compute again the boundary of cd1
			cd2 = copEF * cd1
		end
		for σ ∈ SparseArrays.findnz(cd1)[1]
			# update the counters of used cells
			marks[σ] += 1
		end
		# append a new column to [∂d+]
		# copFC += cd1
		rows, vals = SparseArrays.findnz(cd1)
		jcol += 1
		append!(I,rows)
		append!(J,[ jcol for k=1:nnz(cd1) ])
		append!(W,vals)
	end
	copFC = sparse(I,J,W)
	return copFC
end

"""
    bbox(vertices::Points)

The axis aligned bounding box of the provided set of n-dim `vertices`.

The box is returned as the couple of `Points` of the two opposite corners of the box.
"""
function bbox(vertices::Points)
    minimum = mapslices(x->min(x...), vertices, dims=1)
    maximum = mapslices(x->max(x...), vertices, dims=1)
    minimum, maximum
end

"""
    bbox_contains(container, contained)

Check if the axis aligned bounding box `container` contains `contained`.

Each input box must be passed as the couple of `Points` standing on the opposite corners of the box.
"""
function bbox_contains(container, contained)
    b1_min, b1_max = container
    b2_min, b2_max = contained
    all(map((i,j,k,l)->i<=j<=k<=l, b1_min, b2_min, b2_max, b1_max))
end

"""
    face_area(V::Points, EV::Cells, face::Cell)

The area of `face` given a geometry `V` and an edge topology `EV`.
"""
function face_area(V::Points, EV::Cells, face::Cell)
    return face_area(V, build_copEV(EV), face)
end

function face_area(V::Points, EV::ChainOp, face::Cell)
    function triangle_area(triangle_points::Points)
        ret = ones(3,3)
        ret[:, 1:2] = triangle_points
        return .5*det(ret)
    end

    area = 0

    fv = buildFV(EV, face)

    verts_num = length(fv)
    v1 = fv[1]

    for i in 2:(verts_num-1)

        v2 = fv[i]
        v3 = fv[i+1]

        area += triangle_area(V[[v1, v2, v3], :])
    end

    return area
end

"""
    skel_merge(V1::Points, EV1::ChainOp, V2::Points, EV2::ChainOp)

Merge two **1-skeletons**
"""
function skel_merge(V1::Points, EV1::ChainOp, V2::Points, EV2::ChainOp)
	V = [V1; V2]
    EV = blockdiag(EV1,EV2)
    return V, EV
end

# function coppie_duplicate(A)
# 	t=()
# 	for i in 1:size(A,1)
# 		t= (t...,(A[i,1],A[i,2]))
# 	end
# 	return [y[i] for y in unique(t), i in 1:2]
# end
"""
    skel_merge(V1::Points, EV1::ChainOp, FE1::ChainOp, V2::Points, EV2::ChainOp, FE2::ChainOp)

Merge two **2-skeletons**
"""
function skel_merge(V1::Lar.Points, EV1::Lar.ChainOp, FE1::Lar.ChainOp,
					V2::Lar.Points, EV2::Lar.ChainOp, FE2::Lar.ChainOp)
    FE = blockdiag(FE1,FE2)
    V, EV = skel_merge(V1, EV1, V2, EV2)
    return V, EV, FE
end

"""
    delete_edges(todel, V::Points, EV::ChainOp)

Delete edges and remove unused vertices from a **2-skeleton**.

Loop over the `todel` edge index list and remove the marked edges from `EV`.
The vertices in `V` which remained unconnected after the edge deletion are deleted too.
"""
function delete_edges(todel, V::Points, EV::ChainOp)
    tokeep = setdiff(collect(1:EV.m), todel)
    EV = EV[tokeep, :]

    vertinds = 1:EV.n
    todel = Array{Int64, 1}()
    for i in vertinds
        if length(EV[:, i].nzind) == 0
            push!(todel, i)
        end
    end

    tokeep = setdiff(vertinds, todel)
    EV = EV[:, tokeep]
    V = V[tokeep, :]

    return V, EV
end

"""
    buildFV(EV::Cells, face::Cell)

The list of vertex indices that expresses the given `face`.

The returned list is made of the vertex indices ordered following the traversal order to keep a coherent face orientation.
The edges are need to understand the topology of the face.

In this method the input face must be expressed as a `Cell`(=`SparseVector{Int8, Int}`) and the edges as `Cells`.
"""
function buildFV(EV::Cells, face::Cell)
    return buildFV(build_copEV(EV), face)
end

"""
    buildFV(copEV::ChainOp, face::Cell)

The list of vertex indices that expresses the given `face`.

The returned list is made of the vertex indices ordered following the traversal order to keep a coherent face orientation.
The edges are need to understand the topology of the face.

In this method the input face must be expressed as a `Cell`(=`SparseVector{Int8, Int}`) and the edges as `ChainOp`.
"""
function buildFV(copEV::ChainOp, face::Cell)
    startv = -1
    nextv = 0
    edge = 0

    vs = Array{Int64, 1}()

    while startv != nextv
        if startv < 0
            edge = face.nzind[1]
            startv = copEV[edge,:].nzind[face[edge] < 0 ? 2 : 1]
            push!(vs, startv)
        else
            edge = setdiff(intersect(face.nzind, copEV[:, nextv].nzind), edge)[1]
        end
        nextv = copEV[edge,:].nzind[face[edge] < 0 ? 1 : 2]
        push!(vs, nextv)

    end

    return vs[1:end-1]
end

"""
    buildFV(copEV::ChainOp, face::Array{Int, 1})

The list of vertex indices that expresses the given `face`.

The returned list is made of the vertex indices ordered following the traversal order to keep a coherent face orientation.
The edges are need to understand the topology of the face.

In this method the input face must be expressed as a list of vertex indices and the edges as `ChainOp`.
"""
function buildFV(copEV::ChainOp, face::Array{Int, 1})
    startv = face[1]
    nextv = startv

    vs = []
    visited_edges = []

    while true
        curv = nextv
        push!(vs, curv)

        edge = 0

        for edgeEx in copEV[:, curv].nzind
            nextv = setdiff(copEV[edgeEx, :].nzind, curv)[1]
            if nextv in face && (nextv == startv || !(nextv in vs)) && !(edgeEx in visited_edges)
                edge = edgeEx
                break
            end
        end

        push!(visited_edges, edge)

        if nextv == startv
            break
        end
    end

    return vs
end

"""
    build_copFE(FV::Cells, EV::Cells)

The signed `ChainOp` from 1-cells (edges) to 2-cells (faces)
"""
function build_copFE(FV::Lar.Cells, EV::Lar.Cells)
	copFE = Lar.u_coboundary_1(FV, EV) # unsigned
	faceedges = [findnz(copFE[f,:])[1] for f=1:size(copFE,1)]

	f_edgepairs = Array{Array{Int64,1}}[]
	for f=1:size(copFE,1)
		edgepairs = Array{Int64,1}[]
		for v in FV[f]
			push!(edgepairs, [e for e in faceedges[f] if v in EV[e]])
		end
		push!(f_edgepairs, edgepairs)
	end
	for f=1:size(copFE,1)
		for (e1,e2) in f_edgepairs[f]
			v = intersect(EV[e1], EV[e2])[1]
			copFE[f,e1] = EV[e1][2]==v ? 1 : -1
			copFE[f,e2] = EV[e2][1]==v ? 1 : -1
		end
	end
	return copFE
end

"""
    build_copEV(EV::Cells, signed=true)

The signed (or not) `ChainOp` from 0-cells (vertices) to 1-cells (edges)
"""
function build_copEV(EV::Cells, signed=true)
    setValue = [-1, 1]
    if signed == false
        setValue = [1, 1]
    end

    maxv = max(map(x->max(x...), EV)...)
    copEV = spzeros(Int8, length(EV), maxv)

    for (i,e) in enumerate(EV)
        e = sort(collect(e))
        copEV[i, e] = setValue
    end

    return copEV
end

"""
    build_cops(edges::Cells, faces::Cells)

The vertices-edges and edges-faces chain operators (`copEV::ChainOp`, `copFE::ChainOp`)
"""
function build_cops(edges::Cells, faces::Cells)
    copEV = build_copEV(edges)
    FV = Cells(map(x->buildFV(copEV,x), faces))
    copFE = build_copFE(FV, edges)

    return [copEV, copFE]
end

"""
    vin(vertex, vertices_set)

Checks if `vertex` is one of the vertices inside `vertices_set`
"""
function vin(vertex, vertices_set)
    for v in vertices_set
        if vequals(vertex, v)
            return true
        end
    end
    return false
end

"""
    vequals(v1, v2)

Check the equality between vertex `v1` and vertex `v2`
"""
function vequals(v1, v2)
    err = 10e-8
    return length(v1) == length(v2) && all(map((x1, x2)->-err < x1-x2 < err, v1, v2))
end

function vcycle( copEV::Lar.ChainOp, copFE::Lar.ChainOp, f::Int64 )
	edges,signs = SparseArrays.findnz(copFE[f,:])
	vpairs = [s>0 ? SparseArrays.findnz(copEV[e,:])[1] :
					reverse(SparseArrays.findnz(copEV[e,:])[1])
				for (e,s) in zip(edges,signs)]
	a = [pair for pair in vpairs if length(pair)==2]
	function mycat(a::Lar.Cells)
		out=[]
		for cell in a append!(out,cell) end
		return out
	end
	vs = collect(Set(mycat(a)))
	vdict = Dict(zip(vs,1:length(vs)))
	edges = [[vdict[pair[1]], vdict[pair[2]]] for pair in vpairs if length(pair)==2]
	return vs, edges
end

"""
    triangulate(model::LARmodel)

Full constrained Delaunnay triangulation of the given 3-dimensional `LARmodel`
"""
function triangulate(model::LARmodel)
    V, topology = model
    cc = build_cops(topology...)
    return triangulate(V, cc)
end

"""
    triangulate(V::Points, cc::ChainComplex)

Full constrained Delaunnay triangulation of the given 3-dimensional model (given with topology as a `ChainComplex`)
"""
function triangulate(V::Lar.Points, cc::Lar.ChainComplex)
	copEV, copFE = cc[1:2]

    triangulated_faces = Array{Any, 1}(undef, copFE.m)

    for f in 1:copFE.m
        if f % 10 == 0
            print(".")
        end

        edges_idxs = copFE[f, :].nzind
        edge_num = length(edges_idxs)
        edges = zeros(Int64, edge_num, 2)

        #fv = Lar.buildFV(copEV, copFE[f, :])
        fv, edges = Lar.vcycle(copEV, copFE, f)
		if fv ≠ []
			vs = V[fv, :]
	        v1 = LinearAlgebra.normalize(vs[2, :] - vs[1, :])
	        v2 = [0, 0, 0]
	        v3 = [0, 0, 0]
	        err = 1e-8
	        i = 3
	        while -err < LinearAlgebra.norm(v3) < err
	            v2 = LinearAlgebra.normalize(vs[i, :] - vs[1, :])
	            v3 = LinearAlgebra.cross(v1, v2)
	            i = i % size(vs,1) + 1
	        end
	        M = reshape([v1; v2; v3], 3, 3)
	        vs = (vs*M)[:, 1:2]
			v = convert(Lar.Points, vs'[1:2,:])
			vmap = Dict(zip(fv,1:length(fv))) # vertex map
			mapv = Dict(zip(1:length(fv),fv)) # inverse vertex map

			trias = Lar.triangulate2d(v,edges)
			triangulated_faces[f] = [[mapv[v] for v in tria] for tria in trias]
		end
    end

    return triangulated_faces
end

"""
    point_in_face(point, V::Points, copEV::ChainOp)

Check if `point` is inside the area of the face bounded by the edges in `copEV`
"""
function point_in_face(point, V::Points, copEV::ChainOp)

    function pointInPolygonClassification(V,EV)

        function crossingTest(new, old, status, count)
        if status == 0
            status = new
            return status, (count + 0.5)
        else
            if status == old
                return 0, (count + 0.5)
            else
                return 0, (count - 0.5)
            end
        end
        end

        function setTile(box)
        tiles = [[9,1,5],[8,0,4],[10,2,6]]
        b1,b2,b3,b4 = box
        function tileCode(point)
            x,y = point
            code = 0
            if y>b1 code=code|1 end
            if y<b2 code=code|2 end
            if x>b3 code=code|4 end
            if x<b4 code=code|8 end
            return code
        end
        return tileCode
        end

        function pointInPolygonClassification0(pnt)
            x,y = pnt
            xmin,xmax,ymin,ymax = x,x,y,y
            tilecode = setTile([ymax,ymin,xmax,xmin])
            count,status = 0,0

            for k in 1:EV.m
                edge = EV[k,:]
                p1, p2 = V[edge.nzind[1], :], V[edge.nzind[2], :]
                (x1,y1),(x2,y2) = p1,p2
                c1,c2 = tilecode(p1),tilecode(p2)
                c_edge, c_un, c_int = xor(c1, c2), c1|c2, c1&c2

                if (c_edge == 0) & (c_un == 0) return "p_on"
                elseif (c_edge == 12) & (c_un == c_edge) return "p_on"
                elseif c_edge == 3
                    if c_int == 0 return "p_on"
                    elseif c_int == 4 count += 1 end
                elseif c_edge == 15
                    x_int = ((y-y2)*(x1-x2)/(y1-y2))+x2
                    if x_int > x count += 1
                    elseif x_int == x return "p_on" end
                elseif (c_edge == 13) & ((c1==4) | (c2==4))
                        status, count = crossingTest(1,2,status,count)
                elseif (c_edge == 14) & ((c1==4) | (c2==4))
                        status, count = crossingTest(2,1,status,count)
                elseif c_edge == 7 count += 1
                elseif c_edge == 11 count = count
                elseif c_edge == 1
                    if c_int == 0 return "p_on"
                    elseif c_int == 4
                        status, count = crossingTest(1,2,status,count)
                    end
                elseif c_edge == 2
                    if c_int == 0 return "p_on"
                    elseif c_int == 4
                        status, count = crossingTest(2,1,status,count)
                    end
                elseif (c_edge == 4) & (c_un == c_edge) return "p_on"
                elseif (c_edge == 8) & (c_un == c_edge) return "p_on"
                elseif c_edge == 5
                    if (c1==0) | (c2==0) return "p_on"
                    else
                        status, count = crossingTest(1,2,status,count)
                    end
                elseif c_edge == 6
                    if (c1==0) | (c2==0) return "p_on"
                    else
                        status, count = crossingTest(2,1,status,count)
                    end
                elseif (c_edge == 9) & ((c1==0) | (c2==0)) return "p_on"
                elseif (c_edge == 10) & ((c1==0) | (c2==0)) return "p_on"
                end
            end

            if (round(count)%2)==1
                return "p_in"
            else
                return "p_out"
            end
        end
        return pointInPolygonClassification0
    end

    return pointInPolygonClassification(V, copEV)(point) == "p_in"
end

"""
    lar2obj(V::Points, cc::ChainComplex)

Triangulated OBJ string representation of the model passed as input.

Use this function to export LAR models into OBJ

# Example

```julia
	julia> cube_1 = ([0 0 0 0 1 1 1 1; 0 0 1 1 0 0 1 1; 0 1 0 1 0 1 0 1],
	[[1,2,3,4],[5,6,7,8],[1,2,5,6],[3,4,7,8],[1,3,5,7],[2,4,6,8]],
	[[1,2],[3,4],[5,6],[7,8],[1,3],[2,4],[5,7],[6,8],[1,5],[2,6],[3,7],[4,8]] )

	julia> cube_2 = Lar.Struct([Lar.t(0,0,0.5), Lar.r(0,0,pi/3), cube_1])

	julia> V, FV, EV = Lar.struct2lar(Lar.Struct([ cube_1, cube_2 ]))

	julia> V, bases, coboundaries = Lar.chaincomplex(V,FV,EV)

	julia> (EV, FV, CV), (copEV, copFE, copCF) = bases, coboundaries

	julia> FV # bases[2]
	18-element Array{Array{Int64,1},1}:
	 [1, 3, 4, 6]
	 [2, 3, 5, 6]
	 [7, 8, 9, 10]
	 [1, 2, 3, 7, 8]
	 [4, 6, 9, 10, 11, 12]
	 [5, 6, 11, 12]
	 [1, 4, 7, 9]
	 [2, 5, 11, 13]
	 [2, 8, 10, 11, 13]
	 [2, 3, 14, 15, 16]
	 [11, 12, 13, 17]
	 [11, 12, 13, 18, 19, 20]
	 [2, 3, 13, 17]
	 [2, 13, 14, 18]
	 [15, 16, 19, 20]
	 [3, 6, 12, 15, 19]
	 [3, 6, 12, 17]
	 [14, 16, 18, 20]

	julia> CV # bases[3]
	3-element Array{Array{Int64,1},1}:
	 [2, 3, 5, 6, 11, 12, 13, 14, 15, 16, 18, 19, 20]
	 [2, 3, 5, 6, 11, 12, 13, 17]
	 [1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 17]

	julia> copEV # coboundaries[1]
	34×20 SparseMatrixCSC{Int8,Int64} with 68 stored entries: ...

	julia> copFE # coboundaries[2]
	18×34 SparseMatrixCSC{Int8,Int64} with 80 stored entries: ...

	julia> copCF # coboundaries[3]
	4×18 SparseMatrixCSC{Int8,Int64} with 36 stored entries: ...

	objs = Lar.lar2obj(V'::Lar.Points, [coboundaries...])

	open("./two_cubes.obj", "w") do f
    	write(f, objs)
	end
```
"""
function lar2obj(V::Lar.Points, cc::Lar.ChainComplex)
    copEV, copFE, copCF = cc
	if size(V,2) > 3
		V = convert(Lar.Points, V') # out V by rows
	end
    obj = ""
    for v in 1:size(V, 1)
        obj = string(obj, "v ",
    	round(V[v, 1], digits=6), " ",
    	round(V[v, 2], digits=6), " ",
    	round(V[v, 3], digits=6), "\n")
    end

    print("Triangulating")
    triangulated_faces = Lar.triangulate(V, cc[1:2])
    println("DONE")

    for c in 1:copCF.m
        obj = string(obj, "\ng cell", c, "\n")
        for f in copCF[c, :].nzind
            triangles = triangulated_faces[f]
            for tri in triangles
                #t = copCF[c, f] > 0 ? tri : tri[end:-1:1]
				t = tri
                obj = string(obj, "f ", t[1], " ", t[2], " ", t[3], "\n")
            end
        end
    end

    return obj
end

"""
	lar2obj(V::Lar.Points, TV::Lar.Cells)::String
```
```
"""
function lar2obj(V::Lar.Points, TV::Lar.Cells)
	obj = ""
    for v in 1:size(V, 1)
        obj = string(obj, "v ",
    	round(V[1,v], digits=6), " ",
    	round(V[2,v], digits=6), " ",
    	round(V[3,v], digits=6), "\n")
	end
	for tri in TV
		#t = copCF[c, f] > 0 ? tri : tri[end:-1:1]
		obj = string(obj, "f ", tri[1], " ", tri[2], " ", tri[3], "\n")
	end
end

"""
    obj2lar(path)

Read OBJ file at `path` and create a 2-skeleton as `Tuple{Points, ChainComplex}` from it.

This function does care about eventual internal grouping inside the OBJ file.
"""
function obj2lar(path)
    vs = Array{Float64, 2}(undef, 0, 3)
    edges = Array{Array{Array{Int, 1}, 1}, 1}()
    faces = Array{Array{Array{Int, 1}, 1}, 1}()
	push!(edges, Array{Array{Int, 1}, 1}[])
	push!(faces, Array{Array{Int, 1}, 1}[])
	g = 1

    open(path, "r") do fd
        for line in eachline(fd)
            elems = split(line)
            if length(elems) > 0
                if elems[1] == "v"
					# parse and store a vertex
                    x = parse(Float64, elems[2])
                    y = parse(Float64, elems[3])
                    z = parse(Float64, elems[4])
                    vs = [vs; x y z]

				elseif elems[1] == "f"
                    # Ignore the vertex tangents and normals
                    v1 = parse(Int, split(elems[2], "/")[1])
                    v2 = parse(Int, split(elems[3], "/")[1])
                    v3 = parse(Int, split(elems[4], "/")[1])

                    e1 = sort([v1, v2])
                    e2 = sort([v2, v3])
                    e3 = sort([v3, v1])

					push!(edges[g], e1)
					push!(edges[g], e2)
					push!(edges[g], e3)

                    push!(faces[g], sort([v1, v2, v3]))

				elseif elems[1] == "g"
					#start a new group of edges and faces
					#println(line)
					g += 1
					push!(edges, Array{Array{Int, 1}, 1}[])
					push!(faces, Array{Array{Int, 1}, 1}[])
				end
            end
        end
    end

	#return vs, build_cops(edges, faces)
	return convert(Lar.Points, vs'), edges[2:end], faces[2:end]
end

"""
    binaryRange(n)

Generate the first `n` binary numbers in string padded for max `2^n` length
"""
function binaryRange(n)
    return string.(range(0, length=2^n), base=2, pad=n)
end

function space_arrangement(V::Points, EV::ChainOp, FE::ChainOp, multiproc::Bool=false)

    fs_num = size(FE, 1)
    sp_idx = Lar.Arrangement.spatial_index(V, EV, FE)

    rV = Lar.Points(undef, 0,3)
    rEV = SparseArrays.spzeros(Int8,0,0)
    rFE = SparseArrays.spzeros(Int8,0,0)

    if (multiproc == true)
        in_chan = Distributed.RemoteChannel(()->Channel{Int64}(0))
        out_chan = Distributed.RemoteChannel(()->Channel{Tuple}(0))

        @async begin
            for sigma in 1:fs_num
                put!(in_chan, sigma)
            end
            for p in Distributed.workers()
                put!(in_chan, -1)
            end
        end

        for p in Distributed.workers()
            @async Base.remote_do(
                frag_face_channel, p, in_chan, out_chan, V, EV, FE, sp_idx)
        end

        for sigma in 1:fs_num
            rV, rEV, rFE = skel_merge(rV, rEV, rFE, take!(out_chan)...)
        end

    else

       for sigma in 1:fs_num
           # print(sigma, "/", fs_num, "\r")
           nV, nEV, nFE = Lar.Arrangement.frag_face(
           	V, EV, FE, sp_idx, sigma)
           a,b,c = Lar.skel_merge(
           	rV, rEV, rFE, nV, nEV, nFE)
           rV=a; rEV=b; rFE=c
       end

#		depot_V = Array{Array{Float64,2},1}(undef,fs_num)
#		depot_EV = Array{ChainOp,1}(undef,fs_num)
#		depot_FE = Array{ChainOp,1}(undef,fs_num)
#        for sigma in 1:fs_num
#            print(sigma, "/", fs_num, "\r")
#            nV, nEV, nFE = Arrangement.frag_face( V, EV, FE, sp_idx, sigma)
#            depot_V[sigma] = nV
#            depot_EV[sigma] = nEV
#            depot_FE[sigma] = nFE
#        end
#		rV = vcat(depot_V...)
#		rEV = SparseArrays.blockdiag(depot_EV...)
#		rFE = SparseArrays.blockdiag(depot_FE...)

    end

    rV, rEV, rFE = Lar.Arrangement.merge_vertices(rV, rEV, rFE)

    rCF = Arrangement.minimal_3cycles(rV, rEV, rFE)

    return rV, rEV, rFE, rCF
end

###  2D triangulation
Lar = LinearAlgebraicRepresentation
"""
	obj2lar2D(path::AbstractString)::Lar.LARmodel

Read a *triangulation* from file, given its `path`. Return a `LARmodel` object
"""
function obj2lar2D(path::AbstractString)::Lar.LARmodel
    vs = Array{Float64, 2}(undef, 0, 3)
    edges = Array{Array{Int, 1}, 1}()
    faces = Array{Array{Int, 1}, 1}()

    open(path, "r") do fd
		for line in eachline(fd)
			elems = split(line)
			if length(elems) > 0
				if elems[1] == "v"
					x = parse(Float64, elems[2])
					y = parse(Float64, elems[3])
					z = parse(Float64, elems[4])
					vs = [vs; x y z]
				elseif elems[1] == "f"
					# Ignore the vertex tangents and normals
					v1 = parse(Int, elems[2])
					v2 = parse(Int, elems[3])
					v3 = parse(Int, elems[4])
					append!(edges, map(sort,[[v1,v2],[v2,v3],[v3,v1]]))
					push!(faces, [v1, v2, v3])
				end
				edges = collect(Set(edges))
			end
		end
	end
    return (vs, [edges,faces])::Lar.LARmodel
end

"""
	lar2obj2D(V::Lar.Points,
			cc::Lar.ChainComplex)::String

Produce a *triangulation* from a `LARmodel`. Return a `String` object
"""
function lar2obj2D(V::Lar.Points, cc::Lar.ChainComplex)::String
    @assert length(cc) == 2
    copEV, copFE = cc
    V = [V zeros(size(V, 1))]

    obj = ""
    for v in 1:size(V, 1)
        	obj = string(obj, "v ",
        	round(V[v, 1]; digits=6), " ",
        	round(V[v, 2]; digits=6), " ",
        	round(V[v, 3]; digits=6), "\n")
    end

    #triangulated_faces = triangulate2D(V, cc)
    triangulated_faces = triangulate(V, cc)

	obj = string(obj, "\n")
	for f in 1:copFE.m
		triangles = triangulated_faces[f]
		for tri in triangles
			t = tri
			#t = copCF[c, f] > 0 ? tri : tri[end:-1:1]
			obj = string(obj, "f ", t[1], " ", t[2], " ", t[3], "\n")
		end
	end

    return obj
end

#TODO: finish by using a string as an IObuffer
#"""
#	lar2tria2lar(V::Lar.Points, cc::Lar.ChainComplex)::Lar.LARmodel
#
#Return a triangulated `LARmodel` starting from a stadard LARmodel.
#Useful for colour drawing a complex of non-convex cells.
#
#"""
#function lar2tria2lar(V::Lar.Points, cc::Lar.ChainComplex)::Lar.LARmodel
#	obj = Lar.lar2obj2D(V::Lar.Points, cc::Lar.ChainComplex)
#	vs, (edges,faces) = Lar.obj2lar2D(obj::AbstractString)::Lar.LARmodel
#	return (vs, [edges,faces])::Lar.LARmodel
#end

"""
	triangulate2D(V::Lar.Points,
			cc::Lar.ChainComplex)::Array{Any, 1}

Compute a *CDT* for each face of a `ChainComplex`. Return an `Array` of triangles.
"""
function triangulate2D(V::Lar.Points, cc::Lar.ChainComplex)::Array{Any, 1}
    copEV, copFE = cc
    triangulated_faces = Array{Any, 1}(undef, copFE.m)
    if size(V,2)==2
		V = [V zeros(size(V,1),1)]
	end

	polygons,edgecycles = faces2polygons(copEV, copFE) #new

    for f in 1:copFE.m
        edges_idxs = copFE[f, :].nzind
        edge_num = length(edges_idxs)
        edges = Array{Int64,1}[] #zeros(Int64, edge_num, 2)

		# fv = Lar.buildFV(copEV, copFE[f, :])
		fv = cat(polygons[f])
        vs = V[fv, :]
		edges = cat(edgecycles[f])
        edges = convert(Array{Int64,2}, hcat(edges...)')

		# triangulated_faces[f] = Triangle.constrained_triangulation(
        # 	vs, fv, edges, fill(true, edge_num))
		v = convert(Lar.Points, vs'[1:2,:])
		vmap = Dict(zip(fv,1:length(fv))) # vertex map
		mapv = Dict(zip(1:length(fv),fv)) # inverse vertex map
		ev = [[vmap[e] for e in edges[k,:]] for k=1:size(edges,1)]
		trias = Lar.triangulate2d(v,ev)
		triangulated_faces[f] = [[mapv[v] for v in tria] for tria in trias]

        tV = V[:, 1:2]

        area = Lar.face_area(tV, copEV, copFE[f, :])
        if area < 0
            for i in 1:length(triangulated_faces[f])
                triangulated_faces[f][i] = triangulated_faces[f][i][end:-1:1]
            end
        end
    end

    return triangulated_faces
end

"""
	lar2cop(CV::Lar.Cells)::Lar.ChainOp

Convert an array of array of integer indices to vertices into a sparse matrix.

# Examples

For a single 3D unit cube we get:

```
julia> V,(VV,EV,FV,CV) = Lar.cuboid([1,1,1],true);

julia> Matrix(Lar.lar2cop(EV))
12×8 Array{Int8,2}:
 1  1  0  0  0  0  0  0
 0  0  1  1  0  0  0  0
 0  0  0  0  1  1  0  0
 0  0  0  0  0  0  1  1
 1  0  1  0  0  0  0  0
 0  1  0  1  0  0  0  0
 0  0  0  0  1  0  1  0
 0  0  0  0  0  1  0  1
 1  0  0  0  1  0  0  0
 0  1  0  0  0  1  0  0
 0  0  1  0  0  0  1  0
 0  0  0  1  0  0  0  1

julia> Matrix(Lar.lar2cop(FV))
6×8 Array{Int8,2}:
 1  1  1  1  0  0  0  0
 0  0  0  0  1  1  1  1
 1  1  0  0  1  1  0  0
 0  0  1  1  0  0  1  1
 1  0  1  0  1  0  1  0
 0  1  0  1  0  1  0  1

julia> Matrix(Lar.lar2cop(CV))
1×8 Array{Int8,2}:
 1  1  1  1  1  1  1  1
```
"""
function lar2cop(CV::Lar.Cells)::Lar.ChainOp
	I = Int64[]; J = Int64[]; Value = Int8[];
	for k=1:size(CV,1)
		n = length(CV[k])
		append!(I, k * ones(Int64, n))
		append!(J, CV[k])
		append!(Value, ones(Int64, n))
	end
	return SparseArrays.sparse(I,J,Value)
end

"""
	cop2lar(cop::Lar.ChainOp)::Lar.Cells

Convert a sparse array of type `ChainOp` into an array of array of type `Cells`.

Notice that `cop2lar` is the inverse function of `lar2cop`. their composition is the identity function.

# Example

```
julia> V,(VV,EV,FV,CV) = Lar.cuboid([1,1,1],true);

julia> Lar.cop2lar(Lar.lar2cop(EV))
12-element Array{Array{Int64,1},1}:
 [1, 2]
 [3, 4]
   ...
 [2, 6]
 [3, 7]
 [4, 8]

julia> Lar.cop2lar(Lar.lar2cop(FV))
6-element Array{Array{Int64,1},1}:
 [1, 2, 3, 4]
 [5, 6, 7, 8]
 [1, 2, 5, 6]
 [3, 4, 7, 8]
 [1, 3, 5, 7]
 [2, 4, 6, 8]

julia> Lar.cop2lar(Lar.lar2cop(CV))
1-element Array{Array{Int64,1},1}:
 [1, 2, 3, 4, 5, 6, 7, 8]
```
"""
function cop2lar(cop::Lar.ChainOp)::Lar.Cells
	[findnz(cop[k,:])[1] for k=1:size(cop,1)]
end

function FV2EVs(copEV::Lar.ChainOp, copFE::Lar.ChainOp)
	EV = [findnz(copEV[k,:])[1] for k=1:size(copEV,1)]
	FE = [findnz(copFE[k,:])[1] for k=1:size(copFE,1)]
	EVs = [[EV[e] for e in fe] for fe in FE]
	return EVs
end

"""
	randomcuboids(n,scale=1.0

Generate the `LAR` model of a collection of `n` random cuboids.
Position, orientation and measure of sides are all random.
"""
function randomcuboids(n,scale=1.0)
	assembly = []
	for k=1:n
		corner = rand(Float64, 2)
		sizes = rand(Float64, 2)
		V,(_,EV,_) = Lar.cuboid(corner,true,corner+sizes)
		center = (corner + corner+sizes)/2
		angle = rand(Float64)*2*pi
		obj = Lar.Struct([ Lar.t(center...), Lar.r(angle),
				Lar.s(scale,scale), Lar.t(-center...), (V,EV) ])
		push!(assembly, obj)
	end
	Lar.struct2lar(Lar.Struct(assembly))
end

function randomcubes(n,scale=1.0)
	assembly = []
	for k=1:n
		corner = rand(Float64, 3)
		sizes = rand(Float64, 3)
		V,(_,EV,FV,_) = Lar.cuboid(corner,true,corner+sizes)
		center = (corner + corner+sizes)/2
		angle1 = rand(Float64)*2*pi
		angle2 = rand(Float64)*2*pi
		obj = Lar.Struct([ Lar.t(center...),
				Lar.r(0,angle2,0),Lar.r(0,0,angle1),
				Lar.s(scale,scale,scale), Lar.t(-center...), (V,EV,FV) ])
		push!(assembly, obj)
	end
	Lar.struct2lar(Lar.Struct(assembly))
end

"""
	compute_FV( copEV::Lar.ChainOp, copFE::Lar.ChainOp )::Lar.Cells

Compute the `FV` array of type `Lar.Cells` from two `Lar.ChainOp`, via
sparse array product.  To be generalized to open 2-manifolds.
"""
function compute_FV( copEV::Lar.ChainOp, copFE::Lar.ChainOp )
	# TODO: generalize for open 2-manifolds
	kFV = (x->div(x,2)).(abs.(copFE) * abs.(copEV)) # works only for closed surfaces
	FV = [SparseArrays.findnz(kFV[k,:])[1] for k=1:size(kFV,1)]
	return FV
end

"""
    triangulate2d(V, EV)

Contrained Delaunay Triangulation of LAR model (V,EV).
Discovery and removal of holes from triangulation, by comparing
original and generated edges. `V` is given by column.
"""
# function triangulate2d(V, EV)
#     # data for Constrained Delaunay Triangulation (CDT)
#     points = convert(Array{Float64,2}, V')
# 	points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
#     edges_list = convert(Array{Int64,2}, hcat(EV...)')
#     edge_boundary = [true for k=1:size(edges_list,1)]
#     triangles = Triangle.constrained_triangulation(points,points_map,edges_list)
#     # edges of the triangulation
#     ev = map(sort,cat([[[u,v], [v,w], [w,u]] for (u,v,w) in triangles]))
#     # remove duplicated edges from triangulation
#     ev_nodups = sort(collect(Set(ev)))
#     ##Plasm.view(Plasm.numbering(0.35)((V,[[[k] for k=1:size(V,2)], ev_nodups])))
#     # dictionary o original edges
# 	edge_dict = Dict(zip(EV,1:length(EV)))
# 	#edge_dict = Dict(zip(ev_nodups,1:length(ev_nodups)))
#
# 	triaedges = Array{Int64,1}(undef,0)
# 	for (u,v) in ev
# 		if haskey(edge_dict, [u,v])
# 			push!(triaedges, edge_dict[[u,v]])
# 		elseif haskey(edge_dict, [v,u])
# 			push!(triaedges, edge_dict[[v,u]])
# 		end
# 	end
#
#     # subdivide original edges between inner and outer
#     counters = zeros(size(edges_list,1))
#     for e in triaedges
#         counters[e]+=1
#     end
#     # compute inner triangles
#     inneredges = Array{Array{Int64,1},1}()
#     for (k,value) in enumerate(counters)
#        if value==2
#            push!(inneredges, EV[k], reverse(EV[k]))
#        end
#     end
#     # compute hole(s): wheater all (some?) triangle edges are inneredges
#     holes = Array{Array{Int64,1},1}()
#     for (k,(u,v,w)) in enumerate(triangles)
#         triangle = [[u,v],[v,w],[w,u]]
#         if setdiff(triangle,inneredges)==[]
#             push!(holes, triangles[k])
#         end
#     end
#     triangles = [tria for tria in triangles if !(tria in holes)]
#     return triangles
# end
function triangulate2d(V, EV)
    # data for Constrained Delaunay Triangulation (CDT)
    points = convert(Array{Float64,2}, V')
	points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
    edges_list = convert(Array{Int64,2}, hcat(EV...)')
    edge_boundary = [true for k=1:size(edges_list,1)] ## dead code !!
    trias = Triangle.constrained_triangulation(points,points_map,edges_list)
	innertriangles = Array{Int64,1}[]
	for (u,v,w) in trias
		point = (points[u,:]+points[v,:]+points[w,:])./3
		copEV = Lar.lar2cop(EV)
		inner = Lar.point_in_face(point, points::Lar.Points, copEV::Lar.ChainOp)
		if inner
			push!(innertriangles,[u,v,w])
		end
	end
    return innertriangles
end

function blockdiag_di_nardo(X::SparseMatrixCSC...)
    num = length(X)
    mX = Int[ size(x, 1) for x in X ]
    nX = Int[ size(x, 2) for x in X ]
    m = sum(mX)
    n = sum(nX)

    Tv = promote_type(map(x->eltype(x.nzval), X)...)
    Ti = isempty(X) ? Int : promote_type(map(x->eltype(x.rowval), X)...)

    colptr = Vector{Ti}(undef, n+1)
    nnzX = Int[ nnz(x) for x in X ]
    nnz_res = sum(nnzX)
    rowval = Vector{Ti}(undef, nnz_res)
    nzval = Vector{Tv}(undef, nnz_res)

    nnz_sofar = 0
    nX_sofar = 0
    mX_sofar = 0
    for i = 1 : num
        colptr[(1 : nX[i] + 1) .+ nX_sofar] = X[i].colptr .+ nnz_sofar
        rowval[(1 : nnzX[i]) .+ nnz_sofar] = X[i].rowval .+ mX_sofar
        nzval[(1 : nnzX[i]) .+ nnz_sofar] = X[i].nzval
        nnz_sofar += nnzX[i]
        nX_sofar += nX[i]
        mX_sofar += mX[i]
    end
    colptr[n+1] = nnz_sofar + 1

    SparseMatrixCSC(m, n, colptr, rowval, nzval)
end


text:
using Documenter, InvertibleNetworks

makedocs(sitename="Invertible Networks",
         doctest=false, clean=true,
         authors="Philipp Witte, Ali Siahkoohi, Gabbrio Rizzuti, Mathias Louboutin, Felix J. Herrmann",
         format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
         pages = Any[
             "Home" => "index.md",
             "Examples" => "examples.md",
             "API Reference" => "api.md",
             "LICENSE" => "LICENSE.md",
         ])

deploydocs(repo="github.com/slimgroup/InvertibleNetworks.jl")


text:
data_spec =
[
	(:name, AbstractString, "table"),
	(:format, VegaFormat, nothing),
	(:values, Union{Vector{Dict{Any, Any}}, Vector}, nothing),
	(:source, AbstractString, nothing),
	(:url, AbstractString, nothing),
	(:transform, Vector{VegaTransform}, nothing),
  (:modify, Vector{VegaStream}, nothing)
]

primitivefactory(:VegaData, data_spec)

max_stars_repo_name: jpthor/vega


text:
@testset "Hozz" begin
    @test 13 |> fizz |> buzz |> hozz == "Aho"
    @test 3 |> fizz |> buzz |> hozz == "FizzAho"
    @test 35 |> fizz |> buzz |> hozz == "BuzzAho"
    @test 30 |> fizz |> buzz |> hozz == "FizzBuzzAho"
    @test 30 |> hozz |> buzz |> fizz == "AhoBuzzFizz"
end


text:
using Distributions

"""
	function mcmc(data,r,total::Int)

Metropolis-Hastings MCMC of the SIR ODE model assuming a Poisson likelihood
for the cumulative death rate and a noninformative prior
data = cumulative deaths starting from the first death
r = array of 4 parameters used in the model
total = number of MCMC steps
output = (time series of loglikelihoods, time series of parameters, maximum loglikelihood, max loglikelihood parameters)

"""

function mcmc(data,r,total::Int)

	chi = loss(data,r)
	accept = 0
	rt = copy(r)
	chiml = chi
	rml = copy(r)

	chiout = Array{Float64,1}(undef,total)
	rout = Array{Float64,2}(undef,total,4)
	for step in 1:total

		sampleparams!(rt,r)
		chit = loss(data,rt)

		# mesa prior for psi: P(psi) = 1/50, psi < 50, 1/psi, psi >= 50
		if rand() < exp(chi - chit + log(max(r[4],50)) - log(max(rt[4],50)))
			chi = chit
			r = copy(rt)
			accept += 1
			if chi < chiml
				chiml = chi
				rml = copy(r)
			end
		end
		chiout[step] = chi
		rout[step,:] = r
	end

	return chiout, rout, accept/total,chiml,rml

end

# generate MCMC step guess from current parameter value r using logNormal distirbution
# with scale parameter manually adjusted so approximately a third of guesses are accepted
function sampleparams!(rt,r)
	for i in eachindex(r)
		d = Distributions.LogNormal(log(r[i]),.02)
		rt[i] = rand(d)
	end
end

# loglikelihood of data given parameters at r using a Poisson likelihood function
function loss(data,r)
	prediction = model(r)
	start = max(round(Int,r[4]),0) + 1
	chi = 0
	for i in eachindex(data)
		# d = Distributions.Normal(prediction[start+i-1],sqrt(prediction[start+i-1]))
		# println(max(prediction[start+i-1],eps(Float64)))
		d = Distributions.Poisson(max(prediction[start+i-1],eps(Float64)))
		chi -= loglikelihood(d,[data[i]])
	end
	return chi
end

# model prediction of model for cumulative deaths
function model(r)
	out = oxford(0.,1/6.67e7,round(Int,r[4])+17,.01,r[1],r[2],r[3]*6.67e7,round(Int,r[4]))
	return out[1]
end

function model2(r)
	out = oxfordrk2(0.,1/6.67e7,round(Int,r[4])+17,.01,r[1],r[2],r[3]*6.67e7,round(Int,r[4]))
	return out[1]
end

function modelconstrained(r)
	out = oxfordrk2(0.,1/6.67e7,round(Int,r[4])+17,.01,2.25*r[2],r[2],r[3]*6.67e7,round(Int,r[4]))
	return out[1]
end

# Lourenco  model, Runge-Kutta 2 stepper
function oxfordrk2(zinit,yinit,total,dt,beta,sigma,rtheta,psi)

	z = zinit
	y = yinit

	totalN = Int(total/dt)

	zout = Array{Float64,1}(undef,Int(total+1))
	yout = similar(zout)
	dout = similar(zout)
	tout = similar(zout)

	zout[1] = zinit
	yout[1] = yinit
	tout[1] = 0.
	interval = 1

	t = 1
	# Runge Kutta 2 method
	for n = 1:totalN
		y1 = dt*(beta*y*(1-z)-sigma*y)/2
		z1 = dt*beta*y*(1-z)/2
		y2 = dt*(beta*(y+y1)*(1-z-z1)-sigma*y)
		z2 = dt*beta*(y+y1)*(1-z-z1)
		y += y2
		z += z2
		if n*dt%interval == 0
			t += 1
			zout[t] = z
			yout[t] = y
			dout[t] = rtheta*zout[max(t-psi,0)+1]
			tout[t] = n*dt
		end
	end
	return dout,zout,yout,tout
end

# Lourenco  model, Euler stepper
function oxford(zinit,yinit,total,dt,beta,sigma,rtheta,psi)

	z = zinit
	y = yinit

	totalN = Int(total/dt)

	zout = Array{Float64,1}(undef,Int(total+1))
	yout = similar(zout)
	tout = similar(zout)
	dout = zeros(total+1)

	zout[1] = zinit
	yout[1] = yinit
	tout[1] = 0.
	interval = 1

	t = 1
	# Euler method
	for n = 1:totalN
		y += dt*(beta*y*(1-z)-sigma*y)
		z += dt*beta*y*(1-z)
		if n*dt%interval == 0
			t += 1
			zout[t] = z
			yout[t] = y
			dout[t] = rtheta*zout[max(t-psi,0)+1]
			tout[t] = n*dt
		end
	end
	return dout,zout,yout,tout
end

max_stars_repo_name: cannin/covid-sicr-test

max_issues_repo_name: cannin/covid-sicr-test

max_forks_repo_name: cannin/covid-sicr-test


text:
"""
`incidence(H::SimpleHypergraph)` returns the vertex-edge incidence
matrix of `H`in sparse form. Wrap in `Matrix` for full storage version.
"""
function incidence(H::SimpleHypergraph)
    n = NV(H)
    m = NE(H)
    M = spzeros(Int, n, m)

    VV = vlist(H)
    EE = elist(H)

    try
        sort!(VV)
    catch
    end

    for i = 1:n
        v = VV[i]
        for j = 1:m
            e = EE[j]
            if in(v, e)
                M[i, j] = 1
            end
        end
    end
    return M
end

# conversely, convert a matrix into a hypergraph
"""
`SimpleHypergraph(A)` where `A` is an `n`-by-`m` matrix creates a
hypergraph with `n` vertices and `m` edges determined by the nonzero entries
in the columns of `A`. This is a sort of inverse operation to `incidence`.
"""
function SimpleHypergraph(A::AbstractArray{T,2})::SimpleHypergraph{Int} where {T<:Number}
    n, m = size(A)
    H = IntHypergraph(n)

    for j = 1:m
        a = A[:, j]
        e = Set(findall(a .!= 0))
        add!(H, e)
    end
    return H
end