text:
using ReversibleTuringMachine
using Test
using NiLang

@testset "inst" begin
    tape = [BLANK, 1, 1, 0, 0, 1]
    q = 3
    rule = Quadruple(3, 1, 0, 2)
    @instr ReversibleTuringMachine.instr(q, tape[3], rule, 3)
    @test q == 2
    @test tape[3] == 0
end

@testset "rtm" begin
    prog = quadruples([
        1 BLANK BLANK 2;
        2 SLASH RIGHT 3;
        3 0     1     4;
        3 1     0     2;
        3 BLANK BLANK 4;
        4 SLASH LEFT  5;
        5 0     0     4;
        5 BLANK BLANK 6;
    ])

    @test prog isa Vector{Quadruple{Int,Int}}

    rtm = RTM(1,6,prog)
    @test is_deterministic(rtm)
    @test is_reversible(rtm)

    tape = [BLANK, 1, 1, 0, 0, 1, BLANK]
    loc = 1
    pc = 1
    @instr run!(rtm, tape, loc, pc)
    @test tape == [BLANK, 0, 0, 1, 0, 1, BLANK]
    @instr (~run!)(rtm, tape, loc, pc)
    @test tape == [BLANK, 1, 1, 0, 0, 1, BLANK]
end


text:
module CSV

# stdlib
using Mmap, Dates, Unicode
# Parsers.jl is used for core type parsing from byte buffers
# and all other parsing options (quoted fields, delimiters, dateformats etc.)
using Parsers
# Tables.jl allows integration with all other table/data file formats
using Tables
# PooledArrays.jl is used for materializing pooled columns
using PooledArrays
# FilePathBase allows passing FilePaths instead of just strings for the file name
using FilePathsBase
# WeakRefStrings allows for more efficient materializing of string columns via StringVector
using WeakRefStrings

using CategoricalArrays, DataFrames

struct Error <: Exception
    msg::String
end

Base.showerror(io::IO, e::Error) = println(io, e.msg)

include("utils.jl")
include("detection.jl")
include("header.jl")
include("file.jl")
include("iteration.jl")
include("tables.jl")
include("rows.jl")
include("write.jl")

"""
`CSV.read(source; copycols::Bool=false, kwargs...)` => `DataFrame`

Parses a delimited file into a `DataFrame`. `copycols` determines whether a copy of columns should be made when creating the DataFrame; by default, no copy is made, and the DataFrame is built with immutable, read-only `CSV.Column` vectors. If mutable operations are needed on the DataFrame columns, set `copycols=true`.

`CSV.read` supports the same keyword arguments as [`CSV.File`](@ref).
"""
read(source; copycols::Bool=false, kwargs...) = DataFrame(CSV.File(source; kwargs...), copycols=copycols)

DataFrames.DataFrame(f::CSV.File; copycols::Bool=true) = DataFrame(getcolumns(f), getnames(f); copycols=copycols)

end # module

ext: jl


text:
"""
    Parameters

Abstract type representing a generic set of paramters for finding optimal
two-stage designs.
"""
abstract type Parameters end

# make parameters iterable for automatic broadcasting
Base.size(::Parameters) = ()

Base.length(par::Parameters) = 1

Base.getindex(par::Parameters, i) = par

"""
    label(par::Parameters)

Return the label of a [`Parameters`](@ref)-object.
"""
label(par::Parameters) = try par.label catch error("not implemented") end

"""
    null(par::Parameters)

Return the response rate under the null hypothesis of a [`Parameters`](@ref)-object.
"""
null(par::Parameters) = try par.p0 catch error("not implemented") end

"""
    mtoer(par::Parameters)

Return the maximal type one error rate of a [`Parameters`](@ref)-object.
"""
mtoer(par::Parameters) = try par.mtoer catch error("not implemented") end

"""
    mcrv(par::Parameters)

Return the minimal clinically relevant response rate of a [`Parameters`](@ref)-object.
"""
mcrv(par::Parameters) = try par.pmcrv catch error("not implemented") end

"""
    samplespace(par::Parameters)

Return the [`SampleSpace`](@ref)-object of a [`Parameters`](@ref)-object.
"""
samplespace(par::Parameters) = try par.samplespace catch error("not implemeted") end

"""
    maxsamplesize(par::Parameters)

Return the maximal overall sample size of a [`Parameters`](@ref)-object.
"""
maxsamplesize(par::Parameters) = error("not implemented")

"""
    isgroupsequential(par::Parameters)::Bool

Returns `true` if the undrlying [`SampleSpace`](@ref) is constraint to
group-sequential designs.
"""
isgroupsequential(par::Parameters)::Bool = error("not implemented")

isunimodal(par::Parameters)::Bool = true

"""
    allowsstoppingforfutility(par::Parameters)::Bool

Returns `true` if the undrlying [`SampleSpace`](@ref) allows stopping-for-futlity.
"""
allowsstoppingforfutility(par::Parameters) = error("not implemented")

"""
    NoParameters <: Parameters

This dummy object is used is a [`Design`](@ref) is created manually without
providing information on the respective optimality criterion.
"""
struct NoParameters <: Parameters end

null(par::NoParameters) = error("NoParameters do not have null hypothesis")

mtoer(par::NoParameters) = error("NoParameters do not have mtoer")

samplespace(par::NoParameters) = error("NoParameters do not have sample space")

maxsamplesize(par::NoParameters) = error("NoParameters do not have maximal sample size")

Base.show(io::IO, object::NoParameters) = print("NoParameters")

"""
    PointAlternative <: Parameters

Abstract type representing a generic set of paramters for finding optimal
two-stage designs which are characterized by a point alternative.
"""
abstract type PointAlternative <: Parameters end

alternative(par::PointAlternative) = try par.p1 catch error("not implemented") end

mcrv(par::PointAlternative) = try par.p1 catch error("not implemented") end

mtter(par::PointAlternative) = try par.beta catch error("not implemented") end

"""
    VagueAlternative <: Parameters

Abstract type representing a generic set of paramters for finding optimal
two-stage designs which are characterized by a prior distribution over the
response probability.
"""
abstract type VagueAlternative <: Parameters end

prior(par::VagueAlternative, p::T) where {T<:Real} = try par.prior(p) catch error("not implemented") end

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
using CPUSummary
using Documenter

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

makedocs(;
  modules = [CPUSummary],
  authors = "Chris Elrod <elrodc@gmail.com> and contributors",
  repo = "https://github.com/JuliaSIMD/CPUSummary.jl/blob/{commit}{path}#{line}",
  sitename = "CPUSummary.jl",
  format = Documenter.HTML(;
    prettyurls = get(ENV, "CI", "false") == "true",
    canonical = "https://JuliaSIMD.github.io/CPUSummary.jl",
    assets = String[],
  ),
  pages = ["Home" => "index.md"],
)

deploydocs(; repo = "github.com/JuliaSIMD/CPUSummary.jl")


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

using Markdown
using InteractiveUtils

# ╔═╡ fa4e46d8-b3b0-451c-9a25-a36dfc5b1dfe
begin
	#import Pkg
	#Pkg.activate(mktempdir())
	#using Pkg
	#Pkg.add(["Plots", "PlutoUI", "PyCall"])	# instalação de packages
	using PyCall   # chamar linguagem Python
	using PlutoUI # user-interface do Pluto.jl
	using Plots  # Julia package para gráficos
end

# ╔═╡ 4614a560-9b71-11eb-091e-4f3ed9bc35e3
md"""
# *Notebook*: pmeNotebook1.jl [![forthebadge](https://forthebadge.com/images/badges/made-with-julia.svg)](https://julialang.org/)
##### *Packages*: [Plots](http://docs.juliaplots.org/latest/), [PlutoUI](https://juliahub.com/docs/PlutoUI/abXFp/0.7.6/), [PyCall](https://github.com/JuliaPy/PyCall.jl)
##### Copiar o URL do *notebook* em: [![](https://img.shields.io/badge/PME-GitHub-green.svg)](https://github.com/Ricardo-Luis/PME), para o abrir no seu Julia\Pluto
"""

# ╔═╡ 43efcca9-24ee-49eb-9bbc-1448e86af0e4
html"<button onclick='present()'>present</button>"

# ╔═╡ 98fa8d3c-9632-4cd4-81a4-f4aa107faae6
md"""
# *Notebooks* em *Julia\Pluto*
"""

# ╔═╡ 5570fe10-a878-4686-8c70-bf40741fbafc
md"""
**_Pluto_** é uma interface para *notebook* baseada em navegador *web* (recomendado: Mozilla Firefox ou Google Chrome) para a linguagem de programação *Julia*. Ele permite combinar o código *Julia* e os resultados computacionais com texto, expressões matemáticas, imagens, etc., de uma forma interligada e interactiva com o utilizador.
"""

# ╔═╡ 95173c4c-176d-4790-b236-27de3dfe9cea
md"""
## Porquê *Julia*?
"""

# ╔═╡ 581f1d64-29ee-4851-9330-fe3d4a98cad7
md"""
- linguagem de computação científica _open-source_
- fácil adopção, principalmente para quem conhece MATLAB: ver comparativo entre [MATLAB _vs._ Python _vs._ Julia](https://cheatsheets.quantecon.org/)
- artigo de opinião: [10 Reasons Why You Should Learn Julia](https://blog.goodaudience.com/10-reasons-why-you-should-learn-julia-d786ac29c6ca)
"""

# ╔═╡ 464d6e10-3cca-47b8-aa9b-64e24d1ca14f
md"""
## Porquê *Pluto*?
"""

# ╔═╡ af9761e5-14fb-4cc0-9627-3b251a78da62
md"""
- permite programação literária e reactiva. [Vídeo de apresentação do Pluto, 24min.](https://youtu.be/IAF8DjrQSSk)
- os _notebooks_ criados no *Pluto* são ficheiros **_Julia_**, ou seja, os mesmos podem ser utilizados em editores de código ou em ambiente de desenvolvimento integrado (IDE) como Visual Studio, Atom...
- permite fazer ciência reproduzível!! \
![](https://i.imgur.com/shIpaY1.png) \
Fonte: Xing Su, "Reproducible Research Course Notes", Data Science Specialization, Johns Hopkins University, Fev. 2016.

- artigo de opinião: [A Guide to Building Reactive Notebooks for Scientific Computing With Julia and Pluto.jl](https://medium.com/swlh/a-guide-to-building-reactive-notebooks-for-scientific-computing-with-julia-and-pluto-jl-1a2c0c455d51)
"""

# ╔═╡ d750080e-69a3-4d17-a274-a938be51cd70
md"""
### Limitações do *Pluto*
- o facto de executar tudo... e poder haver partes de longa execução que não seriam necessárias. Mas estão a trabalhar nisso:
- **PlutoCon 2021**: [Visualizing cell dependencies and execution barriers](https://youtu.be/UJRRauTR9Qg?t=827), Benjamin Lungwitz

**Resumo**:
Reactivity for notebooks is great! But in practice, cell dependencies causing re-evaluations are not always obvious (and sometimes even accidential). In addition, there may be long-running cells which should not be updated automatically for each upstream change.  \
In this presentation I show you the current state of cell dependency visualizations and execution barriers in Pluto.

- construir slides para uma apresentação parece ainda ter alguns problemas e exige alguns truques: [Building slides with Pluto.jl](https://andreaskroepelin.de/blog/plutoslides/), Andreas Kröpelin.
"""

# ╔═╡ a2b16233-1a5a-47df-87b8-be8724d7a770
md"""
## Instalação do Julia e Pluto
"""

# ╔═╡ b649e7cf-3153-4929-a64e-2ca64de0eee7
md"""
Versão vídeo e versão escrita de instalação de Julia e Pluto: [First-time setup: Install Julia & Pluto](https://computationalthinking.mit.edu/Spring21/installation/)
"""

# ╔═╡ c86c0119-487c-429d-a7a8-1dc022c794ca
md"""
# FEMM
"""

# ╔═╡ c255a6d2-e66e-41ef-bd9f-ca1f9e248122
md"""
## Exemplo de modelo com CAD paramétrico
"""

# ╔═╡ d46b14ee-7185-42e1-b590-59f85ca9619f
md"""
**Objectivo**: desenhar um núcleo de transformador com a geometria da figura, mas de dimensões variáveis (CAD paramétrico). \
![](https://i.imgur.com/kyvAs6W.jpg)
"""

# ╔═╡ 937a7902-4553-46f5-8e25-3a438719a6a4
md"""
### Importação da biblioteca de comandos do FEMM a partir de *Python*
"""

# ╔═╡ 0adbd32d-c418-4272-a92c-976d6ae95de6
@pyimport femm as f

# ╔═╡ 83337f39-ac9a-49df-8028-5c31ac278fd3
md"""
### Definição de problema no FEMM
"""

# ╔═╡ 37cf28fc-618e-4694-a892-d3427a6f6339
md"""
### Parâmetros CAD
"""

# ╔═╡ a63683ea-29e2-4360-bcbd-0d8a8597e88a
begin
	low=44.1 	#outer legs width, mm
	lcw=42.5 	#center legs width, mm
	hc=111.6 	#half core height, mm
	tc=43.6 	#top/botton core height, mm
	ww=44.5 	#core window width, mm
	cl=2low+lcw+2ww 	#core lenght, mm
	wh=hc-tc 	#half core window height, mm
	cd=44.1 	#core deph, mm
end;

# ╔═╡ 10181196-8cfe-45ac-a1bb-bd58a6ca9b33
begin
	f.openfemm()
	f.newdocument(0)
	f.mi_setgrid(5,"cart")
    f.mi_probdef(0,"millimeters","planar",1.E-8,cd,30,0)
end

# ╔═╡ d30e9a75-d225-4a04-96a7-56b7f7802193
md"""
### Desenho do núcleo
"""

# ╔═╡ 4b90e489-0a38-4f6d-9c46-34f00595e125
begin
	#half core drawing
	f.mi_addnode(0,0)
	f.mi_addnode(low,0)
	f.mi_addnode(low,wh)
	f.mi_addnode(0,wh) #4
	f.mi_addnode(low,hc)
	f.mi_addnode(low+ww,wh)
	f.mi_addnode(low+ww,0)
	f.mi_addnode(low+ww+lcw,0)
	f.mi_addnode(low+ww+lcw,wh) #9
	f.mi_addnode(low+2ww+lcw,wh)
	f.mi_addnode(low+2ww+lcw,0)
	f.mi_addnode(2low+2ww+lcw,0)
	f.mi_addnode(2low+2ww+lcw,wh)
	f.mi_addnode(low+2ww+lcw,hc)
	f.mi_addsegment(0,0,low,0)
	f.mi_addsegment(low,0,low,wh)
	f.mi_addsegment(0,0,0,wh)
	f.mi_addsegment(low,wh,low+ww,wh)
	f.mi_addsegment(low+ww,wh,low+ww,0)  #5
	f.mi_addsegment(low+ww,0,low+ww+lcw,0)
	f.mi_addsegment(low+ww+lcw,0,low+ww+lcw,wh)
	f.mi_addsegment(low+ww+lcw,wh,low+2ww+lcw,wh)
	f.mi_addsegment(low+2ww+lcw,wh,low+2ww+lcw,0)
	f.mi_addsegment(low+2ww+lcw,0,2low+2ww+lcw,0)
	f.mi_addsegment(2low+2ww+lcw,0,2low+2ww+lcw,wh)
	f.mi_addsegment(low,hc,low+2ww+lcw,hc)
	f.mi_addarc(low,hc,0,wh,90,1)
	f.mi_addarc(2low+2ww+lcw,wh,low+2ww+lcw,hc,90,1)
	#2nd half core:
	f.mi_seteditmode("group")
	f.mi_selectgroup(0)
	f.mi_mirror(0,0,2low+2ww+lcw,0)
	f.mi_selectgroup(0)
	f.mi_setgroup(1)
	f.mi_clearselected()
end

# ╔═╡ 9057bc6b-56fb-4ea2-82b3-2603a06e6dca
md"""
### *zoom ; save as, mesh*
"""

# ╔═╡ 1d8380bd-e958-4352-978e-9f9688d2ab93
begin
	f.mi_zoomnatural()
	f.mi_saveas("1phShell.fem")
	f.mi_createmesh()
end;

# ╔═╡ ae36a7f8-6347-4f11-ad47-6c8e4dbf62b4
md"""
### Adicionar materiais
"""

# ╔═╡ d0451715-efdf-4418-bbdc-194e8106f324
begin
	Cduct=10
	Lam_thck=0.025
	f.mi_addmaterial("Metglas2605SA1.025mm", 1, 1, 0, 0, Cduct, Lam_thck, 0, 0.98, 0, 0, 0, 0, 0);
    f.mi_addbhpoint("Metglas2605SA1.025mm",0.0,0.0)
	f.mi_addbhpoint("Metglas2605SA1.025mm",0.111,0.5)
	f.mi_addbhpoint("Metglas2605SA1.025mm",0.259,0.9)
	f.mi_addbhpoint("Metglas2605SA1.025mm",0.424,1.3)
	f.mi_addbhpoint("Metglas2605SA1.025mm",0.576,1.6)
    f.mi_addbhpoint("Metglas2605SA1.025mm",0.689,1.8)
	f.mi_addbhpoint("Metglas2605SA1.025mm",0.758,1.9)
	f.mi_addbhpoint("Metglas2605SA1.025mm",0.858,2.0)
	f.mi_addbhpoint("Metglas2605SA1.025mm",0.915,2.1)
	f.mi_addbhpoint("Metglas2605SA1.025mm",0.959,2.2)
	f.mi_addbhpoint("Metglas2605SA1.025mm",1.059,2.4)
	f.mi_addbhpoint("Metglas2605SA1.025mm",1.110,2.5)
	f.mi_addbhpoint("Metglas2605SA1.025mm",1.148,2.7)
	f.mi_addbhpoint("Metglas2605SA1.025mm",1.185,3.0)
	f.mi_addbhpoint("Metglas2605SA1.025mm",1.219,3.3)
	f.mi_addbhpoint("Metglas2605SA1.025mm",1.251,3.6)
	f.mi_addbhpoint("Metglas2605SA1.025mm",1.292,4.3)
	f.mi_addbhpoint("Metglas2605SA1.025mm",1.321,5.0)
	f.mi_addbhpoint("Metglas2605SA1.025mm",1.350,5.8)
	f.mi_addbhpoint("Metglas2605SA1.025mm",1.387,7.4)
	f.mi_addbhpoint("Metglas2605SA1.025mm",1.387,7.4)
	f.mi_addbhpoint("Metglas2605SA1.025mm",1.415,10.1)
	f.mi_addbhpoint("Metglas2605SA1.025mm",1.442,22.0)
	f.mi_addbhpoint("Metglas2605SA1.025mm",1.449,33.4)
end

# ╔═╡ 035137fa-e80c-4ad8-8173-220394a5b481
f.mi_saveas("1phShell.fem")

# ╔═╡ bda74f3c-8871-49d7-b79f-359ec8522e1d
md"""
.\
.\
.\
... este projecto há-de continuar!
"""

# ╔═╡ e43c564a-2e9e-4aa4-bbda-0c1f1d7c7171
md"""
# TOC
"""

# ╔═╡ 5b468d3e-4650-42fb-afd6-778130b04fb6
TableOfContents()

# ╔═╡ ccf8217f-3a9c-43e1-bde8-536a14710c13
md"""
# Setup
"""

# ╔═╡ f49ecfdc-ca4a-402f-a3e5-0ea541e8b365
md"""
O *package* **PyCall** serve para chamar instruções de programação *Python*, por conseguinte, é necessário ter previamente [*Python*](https://www.python.org/downloads/) instalado.
"""

# ╔═╡ 32a6ddd0-d7f1-42fe-af59-4cfd17e69a56

# ╔═╡ 4a91269d-e022-4039-9b1f-cc1a310160cb
version=VERSION;

# ╔═╡ 74f1620e-ad77-4a7a-88d1-0a789762853f
md"""
PME2\
*Notebook* realizado em linguagem de programação *Julia* versão: $(version)  \
**Ricardo Luís** (Professor Adjunto) \
ISEL, 12/Abr/2021
"""

# ╔═╡ Cell order:
# ╟─4614a560-9b71-11eb-091e-4f3ed9bc35e3
# ╟─43efcca9-24ee-49eb-9bbc-1448e86af0e4
# ╟─98fa8d3c-9632-4cd4-81a4-f4aa107faae6
# ╟─5570fe10-a878-4686-8c70-bf40741fbafc
# ╟─95173c4c-176d-4790-b236-27de3dfe9cea
# ╟─581f1d64-29ee-4851-9330-fe3d4a98cad7
# ╟─464d6e10-3cca-47b8-aa9b-64e24d1ca14f
# ╟─af9761e5-14fb-4cc0-9627-3b251a78da62
# ╟─d750080e-69a3-4d17-a274-a938be51cd70
# ╟─a2b16233-1a5a-47df-87b8-be8724d7a770
# ╟─b649e7cf-3153-4929-a64e-2ca64de0eee7
# ╟─c86c0119-487c-429d-a7a8-1dc022c794ca
# ╟─c255a6d2-e66e-41ef-bd9f-ca1f9e248122
# ╟─d46b14ee-7185-42e1-b590-59f85ca9619f
# ╟─937a7902-4553-46f5-8e25-3a438719a6a4
# ╠═0adbd32d-c418-4272-a92c-976d6ae95de6
# ╟─83337f39-ac9a-49df-8028-5c31ac278fd3
# ╠═10181196-8cfe-45ac-a1bb-bd58a6ca9b33
# ╟─37cf28fc-618e-4694-a892-d3427a6f6339
# ╠═a63683ea-29e2-4360-bcbd-0d8a8597e88a
# ╟─d30e9a75-d225-4a04-96a7-56b7f7802193
# ╠═4b90e489-0a38-4f6d-9c46-34f00595e125
# ╟─9057bc6b-56fb-4ea2-82b3-2603a06e6dca
# ╠═1d8380bd-e958-4352-978e-9f9688d2ab93
# ╟─ae36a7f8-6347-4f11-ad47-6c8e4dbf62b4
# ╠═d0451715-efdf-4418-bbdc-194e8106f324
# ╠═035137fa-e80c-4ad8-8173-220394a5b481
# ╟─bda74f3c-8871-49d7-b79f-359ec8522e1d
# ╟─e43c564a-2e9e-4aa4-bbda-0c1f1d7c7171
# ╠═5b468d3e-4650-42fb-afd6-778130b04fb6
# ╟─ccf8217f-3a9c-43e1-bde8-536a14710c13
# ╠═fa4e46d8-b3b0-451c-9a25-a36dfc5b1dfe
# ╟─f49ecfdc-ca4a-402f-a3e5-0ea541e8b365
# ╟─32a6ddd0-d7f1-42fe-af59-4cfd17e69a56
# ╟─74f1620e-ad77-4a7a-88d1-0a789762853f
# ╟─4a91269d-e022-4039-9b1f-cc1a310160cb

max_stars_repo_name: Ricardo-Luis/PME

max_issues_repo_name: Ricardo-Luis/PME

max_forks_repo_name: Ricardo-Luis/PME


text:
function array2list(cells)
	return PyObject([Any[cell[h] for h=1:length(cell)] for cell in cells])
end

function doublefirst(cells)
	return p.AL([cells[1],cells])
end


text:
# This file was formerly a part of Julia. License is MIT: https://julialang.org/license

## Sparse matrix performance
include("../perfutil.jl")

include("getindex.jl")
sparse_getindex_perf()

include("getindex_skinny.jl")
sparse_getindex_skinny_perf()

include("fem.jl")
fem_perf()

max_stars_repo_name: JuliaTagBot/Benchmarks

max_issues_repo_name: JuliaArchive/Benchmarks

max_forks_repo_name: JuliaArchive/Benchmarks


text:
using InstaRound
using Test

@testset "No Names" begin
    @test round(IGRound, 999; names=false) == "999"
    @test round(IGRound, BigInt(10) ^ 3; names=false) == "1.0K"
    @test round(IGRound, BigInt(10) ^ 4; names=false) == "10.0K"
    @test round(IGRound, BigInt(10) ^ 5; names=false) == "100.0K"
    @test round(IGRound, BigInt(10) ^ 6; names=false) == "1.0M"
    @test round(IGRound, BigInt(10) ^ 7; names=false) == "10.0M"
    @test round(IGRound, BigInt(10) ^ 8; names=false) == "100.0M"
    @test round(IGRound, BigInt(10) ^ 9; names=false) == "1.0B"
    @test round(IGRound, BigInt(10) ^ 10; names=false) == "10.0B"
    @test round(IGRound, BigInt(10) ^ 11; names=false) == "100.0B"
    @test round(IGRound, BigInt(10) ^ 12; names=false) == "1.0t"
    @test round(IGRound, BigInt(10) ^ 13; names=false) == "10.0t"
    @test round(IGRound, BigInt(10) ^ 14; names=false) == "100.0t"
    @test round(IGRound, BigInt(10) ^ 15; names=false) == "1.0q"
    @test round(IGRound, BigInt(10) ^ 16; names=false) == "10.0q"
    @test round(IGRound, BigInt(10) ^ 17; names=false) == "100.0q"
    @test round(IGRound, BigInt(10) ^ 18; names=false) == "1.0Q"
    @test round(IGRound, BigInt(10) ^ 19; names=false) == "10.0Q"
    @test round(IGRound, BigInt(10) ^ 20; names=false) == "100.0Q"
    @test round(IGRound, BigInt(10) ^ 21; names=false) == "1.0s"
    @test round(IGRound, BigInt(10) ^ 22; names=false) == "10.0s"
    @test round(IGRound, BigInt(10) ^ 23; names=false) == "100.0s"
    @test round(IGRound, BigInt(10) ^ 24; names=false) == "1.0S"
    @test round(IGRound, BigInt(10) ^ 25; names=false) == "10.0S"
    @test round(IGRound, BigInt(10) ^ 26; names=false) == "100.0S"
    @test round(IGRound, BigInt(10) ^ 27; names=false) == "1.0o"
    @test round(IGRound, BigInt(10) ^ 28; names=false) == "10.0o"
    @test round(IGRound, BigInt(10) ^ 29; names=false) == "100.0o"
    @test round(IGRound, BigInt(10) ^ 30; names=false) == "1.0n"
    @test round(IGRound, BigInt(10) ^ 31; names=false) == "10.0n"
    @test round(IGRound, BigInt(10) ^ 32; names=false) == "100.0n"
    @test round(IGRound, BigInt(10) ^ 33; names=false) == "1.0d"
    @test round(IGRound, BigInt(10) ^ 34; names=false) == "10.0d"
    @test round(IGRound, BigInt(10) ^ 35; names=false) == "100.0d"
    @test round(IGRound, BigInt(10) ^ 36; names=false) == "1.0U"
    @test round(IGRound, BigInt(10) ^ 37; names=false) == "10.0U"
    @test round(IGRound, BigInt(10) ^ 38; names=false) == "100.0U"
    @test round(IGRound, BigInt(10) ^ 39; names=false) == "1.0D"
    @test round(IGRound, BigInt(10) ^ 40; names=false) == "10.0D"
    @test round(IGRound, BigInt(10) ^ 41; names=false) == "100.0D"
    @test round(IGRound, BigInt(10) ^ 42; names=false) == "1.0T"
    @test round(IGRound, BigInt(10) ^ 43; names=false) == "10.0T"
    @test round(IGRound, BigInt(10) ^ 44; names=false) == "100.0T"
    @test round(IGRound, BigInt(10) ^ 45; names=false) == "1.0Qt"
    @test round(IGRound, BigInt(10) ^ 46; names=false) == "10.0Qt"
    @test round(IGRound, BigInt(10) ^ 47; names=false) == "100.0Qt"
    @test round(IGRound, BigInt(10) ^ 48; names=false) == "1.0Qd"
    @test round(IGRound, BigInt(10) ^ 49; names=false) == "10.0Qd"
    @test round(IGRound, BigInt(10) ^ 50; names=false) == "100.0Qd"
    @test round(IGRound, BigInt(10) ^ 51; names=false) == "1.0Sd"
    @test round(IGRound, BigInt(10) ^ 52; names=false) == "10.0Sd"
    @test round(IGRound, BigInt(10) ^ 53; names=false) == "100.0Sd"
    @test round(IGRound, BigInt(10) ^ 54; names=false) == "1.0St"
    @test round(IGRound, BigInt(10) ^ 55; names=false) == "10.0St"
    @test round(IGRound, BigInt(10) ^ 56; names=false) == "100.0St"
    @test round(IGRound, BigInt(10) ^ 57; names=false) == "1.0O"
    @test round(IGRound, BigInt(10) ^ 58; names=false) == "10.0O"
    @test round(IGRound, BigInt(10) ^ 59; names=false) == "100.0O"
    @test round(IGRound, BigInt(10) ^ 60; names=false) == "1.0N"
    @test round(IGRound, BigInt(10) ^ 61; names=false) == "10.0N"
    @test round(IGRound, BigInt(10) ^ 62; names=false) == "100.0N"
    @test round(IGRound, BigInt(10) ^ 63; names=false) == "1.0v"
    @test round(IGRound, BigInt(10) ^ 64; names=false) == "10.0v"
    @test round(IGRound, BigInt(10) ^ 65; names=false) == "100.0v"
    @test round(IGRound, BigInt(10) ^ 66; names=false) == "1.0c"
    @test round(IGRound, BigInt(10) ^ 67; names=false) == "10.0c"
    @test round(IGRound, BigInt(10) ^ 68; names=false) == "100.0c"
end

@testset "With Names" begin
    @test round(IGRound, 999; names=true) == "999"
    @test round(IGRound, BigInt(10) ^ 3; names=true) == "1.0 Thousand"
    @test round(IGRound, BigInt(10) ^ 4; names=true) == "10.0 Thousand"
    @test round(IGRound, BigInt(10) ^ 5; names=true) == "100.0 Thousand"
    @test round(IGRound, BigInt(10) ^ 6; names=true) == "1.0 Million"
    @test round(IGRound, BigInt(10) ^ 7; names=true) == "10.0 Million"
    @test round(IGRound, BigInt(10) ^ 8; names=true) == "100.0 Million"
    @test round(IGRound, BigInt(10) ^ 9; names=true) == "1.0 Billion"
    @test round(IGRound, BigInt(10) ^ 10; names=true) == "10.0 Billion"
    @test round(IGRound, BigInt(10) ^ 11; names=true) == "100.0 Billion"
    @test round(IGRound, BigInt(10) ^ 12; names=true) == "1.0 Trillion"
    @test round(IGRound, BigInt(10) ^ 13; names=true) == "10.0 Trillion"
    @test round(IGRound, BigInt(10) ^ 14; names=true) == "100.0 Trillion"
    @test round(IGRound, BigInt(10) ^ 15; names=true) == "1.0 Quadrillion"
    @test round(IGRound, BigInt(10) ^ 16; names=true) == "10.0 Quadrillion"
    @test round(IGRound, BigInt(10) ^ 17; names=true) == "100.0 Quadrillion"
    @test round(IGRound, BigInt(10) ^ 18; names=true) == "1.0 Quintillion"
    @test round(IGRound, BigInt(10) ^ 19; names=true) == "10.0 Quintillion"
    @test round(IGRound, BigInt(10) ^ 20; names=true) == "100.0 Quintillion"
    @test round(IGRound, BigInt(10) ^ 21; names=true) == "1.0 Sextillion"
    @test round(IGRound, BigInt(10) ^ 22; names=true) == "10.0 Sextillion"
    @test round(IGRound, BigInt(10) ^ 23; names=true) == "100.0 Sextillion"
    @test round(IGRound, BigInt(10) ^ 24; names=true) == "1.0 Septillion"
    @test round(IGRound, BigInt(10) ^ 25; names=true) == "10.0 Septillion"
    @test round(IGRound, BigInt(10) ^ 26; names=true) == "100.0 Septillion"
    @test round(IGRound, BigInt(10) ^ 27; names=true) == "1.0 Octillion"
    @test round(IGRound, BigInt(10) ^ 28; names=true) == "10.0 Octillion"
    @test round(IGRound, BigInt(10) ^ 29; names=true) == "100.0 Octillion"
    @test round(IGRound, BigInt(10) ^ 30; names=true) == "1.0 Nonillion"
    @test round(IGRound, BigInt(10) ^ 31; names=true) == "10.0 Nonillion"
    @test round(IGRound, BigInt(10) ^ 32; names=true) == "100.0 Nonillion"
    @test round(IGRound, BigInt(10) ^ 33; names=true) == "1.0 Decillion"
    @test round(IGRound, BigInt(10) ^ 34; names=true) == "10.0 Decillion"
    @test round(IGRound, BigInt(10) ^ 35; names=true) == "100.0 Decillion"
    @test round(IGRound, BigInt(10) ^ 36; names=true) == "1.0 Undecillion"
    @test round(IGRound, BigInt(10) ^ 37; names=true) == "10.0 Undecillion"
    @test round(IGRound, BigInt(10) ^ 38; names=true) == "100.0 Undecillion"
    @test round(IGRound, BigInt(10) ^ 39; names=true) == "1.0 Duodecillion"
    @test round(IGRound, BigInt(10) ^ 40; names=true) == "10.0 Duodecillion"
    @test round(IGRound, BigInt(10) ^ 41; names=true) == "100.0 Duodecillion"
    @test round(IGRound, BigInt(10) ^ 42; names=true) == "1.0 Tredecillion"
    @test round(IGRound, BigInt(10) ^ 43; names=true) == "10.0 Tredecillion"
    @test round(IGRound, BigInt(10) ^ 44; names=true) == "100.0 Tredecillion"
    @test round(IGRound, BigInt(10) ^ 45; names=true) == "1.0 Quattuordecillion"
    @test round(IGRound, BigInt(10) ^ 46; names=true) == "10.0 Quattuordecillion"
    @test round(IGRound, BigInt(10) ^ 47; names=true) == "100.0 Quattuordecillion"
    @test round(IGRound, BigInt(10) ^ 48; names=true) == "1.0 Quindecillion"
    @test round(IGRound, BigInt(10) ^ 49; names=true) == "10.0 Quindecillion"
    @test round(IGRound, BigInt(10) ^ 50; names=true) == "100.0 Quindecillion"
    @test round(IGRound, BigInt(10) ^ 51; names=true) == "1.0 Sexdecillion"
    @test round(IGRound, BigInt(10) ^ 52; names=true) == "10.0 Sexdecillion"
    @test round(IGRound, BigInt(10) ^ 53; names=true) == "100.0 Sexdecillion"
    @test round(IGRound, BigInt(10) ^ 54; names=true) == "1.0 Septendecillion"
    @test round(IGRound, BigInt(10) ^ 55; names=true) == "10.0 Septendecillion"
    @test round(IGRound, BigInt(10) ^ 56; names=true) == "100.0 Septendecillion"
    @test round(IGRound, BigInt(10) ^ 57; names=true) == "1.0 Octodecillion"
    @test round(IGRound, BigInt(10) ^ 58; names=true) == "10.0 Octodecillion"
    @test round(IGRound, BigInt(10) ^ 59; names=true) == "100.0 Octodecillion"
    @test round(IGRound, BigInt(10) ^ 60; names=true) == "1.0 Novemdecillion"
    @test round(IGRound, BigInt(10) ^ 61; names=true) == "10.0 Novemdecillion"
    @test round(IGRound, BigInt(10) ^ 62; names=true) == "100.0 Novemdecillion"
    @test round(IGRound, BigInt(10) ^ 63; names=true) == "1.0 Vigintillion"
    @test round(IGRound, BigInt(10) ^ 64; names=true) == "10.0 Vigintillion"
    @test round(IGRound, BigInt(10) ^ 65; names=true) == "100.0 Vigintillion"
    @test round(IGRound, BigInt(10) ^ 66; names=true) == "1.0 Unvigintillion"
    @test round(IGRound, BigInt(10) ^ 67; names=true) == "10.0 Unvigintillion"
    @test round(IGRound, BigInt(10) ^ 68; names=true) == "100.0 Unvigintillion"
end

@testset "With floating points" begin
    @test round(IGRound, 4_123_456) == "4.12M"
end


text:
using ExSup
using Test

# @testset "ExSu.jl" begin
#     # Write your own tests here.
# end

include("experiment_support_test.jl")


text:
const spectroscopic = "spdfghiklmnoqrtuvwxyz"
function spectroscopic_label(ℓ)
    ℓ >= 0 || throw(DomainError(ℓ, "Negative ℓ values not allowed"))
    ℓ + 1 ≤ length(spectroscopic) ? spectroscopic[ℓ+1] : "[$(ℓ)]"
end


text:
using Construct
using Test

@testset "Construct.jl" begin
    # Write your own tests here.
end


text:
# DESCRIPTION
# Runge-Kutta 4 integration scheme.
#
# INPUT
# - fdyn: name of system dynamics function
# - finp: name of control input function
# - t: time
# - dt: time step
# - x_curr: current state vector
# - u0: initial control vector
#
# OUTPUT
# - x_new: new state vector after integration
#
# ------------------------------------------------------------------------------

function rk4(fdyn,finp,t,dt,x_curr,u0)
    #
    u=finp(t)+u0
    xd, y = fdyn(x_curr,u)
    k1=dt*xd
    #
    u=finp(t+0.5*dt)+u0
    xd, y = fdyn(x_curr+0.5*k1,u)
    k2=dt*xd
    #
    xd, y = fdyn(x_curr+0.5*k2,u)
    k3=dt*xd
    #
    u=finp(t+dt)+u0
    xd, y = fdyn(x_curr+k3,u)
    k4=dt*xd
    # weighted average
    x_new=x_curr+k1/6+k2/3+k3/3+k4/6
    #
    return x_new
end

max_issues_repo_name: umbertosaetti/J-SimpleHel

max_forks_repo_name: umbertosaetti/J-SimpleHel


text:
function bonferroni(pValues)
  @check_pValues pValues
  min(1, pValues * length(pValues))
end


text: Dict{String,Any}("zero" => Dict{String,Any}("value" => "0","type" => "integer"),"posanswer" => Dict{String,Any}("value" => "42","type" => "integer"),"answer" => Dict{String,Any}("value" => "42","type" => "integer"),"neganswer" => Dict{String,Any}("value" => "-42","type" => "integer"))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using Documenter, PropertyFiles

makedocs(
    modules = [PropertyFiles],
    format = Documenter.HTML(),
    checkdocs = :exports,
    sitename = "PropertyFiles.jl",
    pages = ["Home" => "index.md",
            "Interface" => "lib/interface.md"]
)

deploydocs(
    repo = "github.com/cserteGT3/PropertyFiles.jl.git",
)


text:
@doc raw"""
    ∂log(Z::Rotor)

Return the gradient of `log(Z)` with respect to the components of `Z`.

The result includes "off-shell" components of the gradient, meaning that even
though change of `Z` in a direction that changes its norm would not be allowed
for a `Rotor`, we measure the gradient in that direction anyway.  That is, the
elements of the returned vector of quaternions is
```math
\begin{aligned}
  \left[
    \frac{\partial} {\partial Z_w} \log(Z),
    \frac{\partial} {\partial Z_x} \log(Z),
    \frac{\partial} {\partial Z_y} \log(Z),
    \frac{\partial} {\partial Z_z} \log(Z)
  \right].
\end{aligned}
```

Note that, even though `log(::Rotor)` is a `QuatVec`, the derivative (and
therefore each element of the result) is a general `Quaternion`.

See also [`∂exp`](@ref) for a similar function, as well as [`log∂log`](@ref)
for a function to compute the value along with the gradient.

# Examples
```julia
julia> ∂log∂w, ∂log∂x, ∂log∂y, ∂log∂z = ∂log(randn(QuatVecF64));

```
"""
function ∂log(Z::Rotor)
    # log(Z::Rotor) = log(abs2) / 2 + f(Z) * Z.vec
    # f(Z) = atan(absvec / Z.w) / absvec = acos(Z.w) / absvec
    a2 = abs2vec(Z)
    a = sqrt(a2)
    if a < 2eps(typeof(a))
        return Quaternion.([Quaternion(one(a)), zero(a)*imx, zero(a)*imy, one(a)*imz])
    end
    f = acos(Z.w) / a
    fprime = (Z.w - f) / a2
    V = 1 + fprime * QuatVec(Z)
    [
        Quaternion(conj(Z)),
        Z.x * V + f * imx,
        Z.y * V + f * imy,
        Z.z * V + f * imz,
    ]
end

"""
    log∂log(Z::Rotor)

Return the value and gradient of `log(Z)` with respect to the components of `Z`.

See [`∂log`](@ref) for more explanation of the components of the gradient.

# Examples
```julia
julia> l, ∂l = log∂log(randn(RotorF64));

```
"""
function log∂log(Z::Rotor)
    # log(Z::Rotor) = log(abs2) / 2 + f(Z) * Z.vec
    # f(Z) = atan(absvec / Z.w) / absvec = acos(Z.w) / absvec
    a2 = abs2vec(Z)
    a = sqrt(a2)
    if a < 2eps(typeof(a))
        return (
            QuatVec(zero(a)),
            Quaternion.([Quaternion(one(a)), zero(a)*imx, zero(a)*imy, one(a)*imz])
        )
    end
    f = acos(Z.w) / a
    fprime = (Z.w - f) / a2
    V = 1 + fprime * QuatVec(Z)
    ∂log = [
        Quaternion(conj(Z)),
        Z.x * V + f * imx,
        Z.y * V + f * imy,
        Z.z * V + f * imz,
    ]
    (f * QuatVec(Z), ∂log)
end

@doc raw"""
    ∂exp(Z::QuatVec)

Return the gradient of `exp(Z)` with respect to the components of `Z`.

The result includes "off-shell" components of the gradient, meaning that even
though a scalar component of `Z` would not be allowed for a `QuatVec`, we
measure the gradient in that direction anyway.  That is, the first element
of the returned vector of quaternions is
```math
\begin{aligned}
  \left.\frac{\partial} {\partial Z_w} \exp(Z) \right|_{Z_w=0}.
\end{aligned}
```

Note that, even though `exp(::QuatVec)` is a `Rotor`, the derivative (and
therefore each element of the result) is a general `Quaternion`.

See also [`∂log`](@ref) for a similar function, as well as [`exp∂exp`](@ref)
for a function to compute the value along with the gradient.

# Examples
```julia
julia> ∂exp∂w, ∂exp∂x, ∂exp∂y, ∂exp∂z = ∂exp(randn(QuatVecF64));

```
"""
function ∂exp(Z::QuatVec)
    # exp(Z::Quaternion) = exp(Z.w) * cos(absvec) + g(Z) * Z.vec
    # g(Z) = exp(Z.w) * sin(absvec) / absvec
    a2 = abs2vec(Z)
    a = sqrt(a2)
    if a < 2eps(typeof(a))
        return Quaternion.([Quaternion(one(a)), one(a)*imx, one(a)*imy, one(a)*imz])
    end
    s, c = sincos(a)
    sinca = s / a
    g = sinca
    gprime = (c - g) / a2
    V = -sinca + gprime * Z
    [
        c + sinca * Z,
        Z.x * V + g * imx,
        Z.y * V + g * imy,
        Z.z * V + g * imz,
    ]
end

"""
    exp∂exp(Z::QuatVec)

Return the value and gradient of `exp(Z)` with respect to the components of `Z`.

See [`∂exp`](@ref) for more explanation of the components of the gradient.

# Examples
```julia
julia> e, ∂e = exp∂exp(randn(QuatVecF64));

```
"""
function exp∂exp(Z::QuatVec)
    # exp(Z::Quaternion) = exp(Z.w) * cos(absvec) + g(Z) * Z.vec
    # g(Z) = exp(Z.w) * sin(absvec) / absvec
    a2 = abs2vec(Z)
    a = sqrt(a2)
    if a < 2eps(typeof(a))
        return (
            Quaternion(one(a)),
            Quaternion.([Quaternion(one(a)), one(a)*imx, one(a)*imy, one(a)*imz])
        )
    end
    s, c = sincos(a)
    sinca = s / a
    g = sinca
    gprime = (c - g) / a2
    V = -sinca + gprime * Z
    ∂exp = [
        c + sinca * Z,
        Z.x * V + g * imx,
        Z.y * V + g * imy,
        Z.z * V + g * imz,
    ]
    (c+g*Z, ∂exp)
end

@doc raw"""
    slerp∂slerp(q₁, q₂, τ)

Return the value and gradient of `slerp`.

The gradient is with respect to each of the input arguments in turn, with each
quaternion regarded as a series of four arguments.  That is, a total of 10
quaternions will be returned:
```math
\begin{aligned}
  &\big[\\
    &\quad \mathrm{slerp}(q₁, q₂, τ), \\
    &\quad \frac{\partial}{\partial q₁.w} \mathrm{slerp}(q₁, q₂, τ), \\
    &\quad \frac{\partial}{\partial q₁.x} \mathrm{slerp}(q₁, q₂, τ), \\
    &\quad \frac{\partial}{\partial q₁.y} \mathrm{slerp}(q₁, q₂, τ), \\
    &\quad \frac{\partial}{\partial q₁.z} \mathrm{slerp}(q₁, q₂, τ), \\
    &\quad \frac{\partial}{\partial q₂.w} \mathrm{slerp}(q₁, q₂, τ), \\
    &\quad \frac{\partial}{\partial q₂.x} \mathrm{slerp}(q₁, q₂, τ), \\
    &\quad \frac{\partial}{\partial q₂.y} \mathrm{slerp}(q₁, q₂, τ), \\
    &\quad \frac{\partial}{\partial q₂.z} \mathrm{slerp}(q₁, q₂, τ), \\
    &\quad \frac{\partial}{\partial \tau} \mathrm{slerp}(q₁, q₂, τ), \\
  &\big]
\end{aligned}
```
For convenience, this will be a 4-tuple with the `slerp` as the first element,
the first four components of the derivative, followed by the next four
components of the derivative, followed by the last component of the derivative.

See also [`slerp`](@ref) for just the value, and [`slerp∂slerp∂τ`](@ref) for
just the value and the derivative with respect to `τ`.

# Examples
```julia
julia> (q1, q2), τ = randn(RotorF64, 2), rand();

julia> s, ∂s∂q₁, ∂s∂q₂, ∂s∂τ = slerp∂slerp(q₁, q₂, τ);

```
"""
function slerp∂slerp(q₁::Rotor, q₂::Rotor, τ::Real)
    r = q₂ / q₁
    basis = typeof(r)[1, imx, imy, imz]
    l, ∂l = log∂log(r)
    e, ∂e = exp∂exp(τ*l)
    s = e * q₁
    (
        s,
        [
            sum(
                ∂e[b] * τ * ∂l[c][b] * (
                    q₂ * conj(basis[a]) #/ abs2(q₁)
                    - 2q₁[a] * (q₂/q₁) #/ abs2(q₁)
                )[c] * q₁
                for b in 1:4
                for c in 1:4
            ) + e * basis[a]
            for a in 1:4
        ],
        [
            sum(
                ∂e[b] * τ * ∂l[c][b] * (
                    basis[a] / q₁
                )[c] * q₁
                for b in 1:4
                for c in 1:4
            )
            for a in 1:4
        ],
        l * s
    )
end

"""
    slerp∂slerp∂τ(q₁, q₂, τ)

Return the value and time-derivative of `slerp`.

See also [`slerp∂slerp`](@ref), which returns the value and *all* of the
derivatives of `slerp`.

"""
function slerp∂slerp∂τ(q₁::Rotor, q₂::Rotor, τ::Real)
    l = log(q₂ / q₁)
    e = exp(τ*l)
    s = e * q₁
    ∂s∂t = l * s
    (s, ∂s∂t)
end

"""
    squad∂squad∂t(qᵢ, A, B, qᵢ₊₁, ta, tb, t)

Compute the value and time-derivative of [`squad`](@ref).

This is primarily an internal helper function, taking various parameters
computed within the `squad` function.  This will be used to compute the
derivative of `squad` when the angular velocity is also requested.  To actually
obtain the derivative, simply pass the relevant keyword to the `squad`
function.

"""
function squad∂squad∂t(qᵢ, A, B, qᵢ₊₁, ta, tb, t)
    # squad = slerp(slerp(qᵢ, qᵢ₊₁, τ), slerp(A, B, τ), 2τ*(1-τ))
    # τ = (t - ta) / (tb - ta)
    # ∂τ∂t = 1 / (tb - ta)
    # ∂(2τ*(1-τ))/∂t = ∂τ/∂t * (2 - 4τ)

    τ = (t - ta) / (tb - ta)
    ∂τ∂t = 1 / (tb - ta)

    X, ∂X∂τ = slerp∂slerp∂τ(qᵢ, qᵢ₊₁, τ)
    ∂X∂t = ∂τ∂t * ∂X∂τ

    Y, ∂Y∂τ = slerp∂slerp∂τ(A, B, τ)
    ∂Y∂t = ∂τ∂t * ∂Y∂τ

    s, ∂s∂X, ∂s∂Y, ∂s∂τparabola = slerp∂slerp(X, Y, 2τ*(1-τ))
    ∂s∂t = ∂τ∂t * (2 - 4τ) * ∂s∂τparabola

    (
        s,
        sum(∂X∂t[b] * ∂s∂X[b] for b in 1:4)
        + sum(∂Y∂t[b] * ∂s∂Y[b] for b in 1:4)
        + ∂s∂t
    )
end


text:
################################################################################
#  Copyright 2021, Marta Vanin                                                 #
################################################################################
function plot_result_caseA(path_to_caseA_csv::String, which_plot::String, ymin::Float64, ymax::Float64)

    casea = CSV.read(path_to_caseA_csv)

    casea_wls = casea |> Query.@filter(_.criterion .== "wls") |> DataFrames.DataFrame
    casea_rwls = casea |> Query.@filter(_.criterion .== "rwls") |> DataFrames.DataFrame
    casea_rwlav = casea |> Query.@filter(_.criterion .== "rwlav") |> DataFrames.DataFrame

    ACP_df_wls = casea_wls |> Query.@filter(_.eq_model .== "rACP") |> DataFrames.DataFrame
    ACR_df_wls = casea_wls |> Query.@filter(_.eq_model .== "rACR") |> DataFrames.DataFrame
    IVR_df_wls = casea_wls |> Query.@filter(_.eq_model .== "rIVR") |> DataFrames.DataFrame

    ACP_df_rwls = casea_rwls |> Query.@filter(_.eq_model .== "rACP") |> DataFrames.DataFrame
    ACR_df_rwls = casea_rwls |> Query.@filter(_.eq_model .== "rACR") |> DataFrames.DataFrame
    IVR_df_rwls = casea_rwls |> Query.@filter(_.eq_model .== "rIVR") |> DataFrames.DataFrame

    ACP_df_rwlav = casea_rwlav |> Query.@filter(_.eq_model .== "rACP") |> DataFrames.DataFrame
    ACR_df_rwlav = casea_rwlav |> Query.@filter(_.eq_model .== "rACR") |> DataFrames.DataFrame
    IVR_df_rwlav = casea_rwlav |> Query.@filter(_.eq_model .== "rIVR") |> DataFrames.DataFrame

    if which_plot == "time"
        p1 = Plots.scatter([ACP_df_wls.n_bus, ACR_df_wls.n_bus, IVR_df_wls.n_bus], [ACP_df_wls.solve_time, ACR_df_wls.solve_time, IVR_df_wls.solve_time], markershape=[:circle :rect :utriangle], label=["ACP" "ACR" "IVR"], ylabel="Solve time [s]",
                yscale=:log, legend=false, ylims=(ymin, ymax), title="WLS")
        p2 = Plots.scatter([ACP_df_rwls.n_bus, ACR_df_rwls.n_bus, IVR_df_rwls.n_bus], [ACP_df_rwls.solve_time, ACR_df_rwls.solve_time, IVR_df_rwls.solve_time], markershape=[:circle :rect :utriangle], label=["ACP" "ACR" "IVR"],
                legend=:bottomright, yscale=:log, xlabel="number of buses [-]", ylims=(ymin, ymax), yaxis=nothing, title="rWLS")
        p3 = Plots.scatter([ACP_df_rwlav.n_bus, ACR_df_rwlav.n_bus, IVR_df_rwlav.n_bus], [ACP_df_rwlav.solve_time, ACR_df_rwlav.solve_time, IVR_df_rwlav.solve_time], markershape=[:circle :rect :utriangle],
                legend=false, yscale=:log, ylims=(ymin, ymax), yaxis=nothing, title="rWLAV")
        Plots.plot(p1, p2, p3, layout = (1,3))
    elseif which_plot == "error_ph1"
        Plots.scatter([ACR_df_rwlav.n_bus, ACR_df_rwlav.n_bus], [ACR_df_rwlav.err_max_1, ACR_df_rwlav.err_avg_1], markershape=[:circle :utriangle], label=["max. abs. error" "avg. abs. error"], ylabel="Absolute error ϵ [p.u.]", xlabel="Number of buses [-]",
                legend=:topright, title="Error plot for case study A - Phase 1", ylims = (ymin, ymax))
    elseif which_plot == "error_ph2"
        Plots.scatter([ACR_df_rwlav.n_bus, ACR_df_rwlav.n_bus], [ACR_df_rwlav.err_max_2, ACR_df_rwlav.err_avg_2], markershape=[:circle :utriangle], label=["max. abs. error" "avg. abs. error"], ylabel="Absolute error ϵ [p.u.]", xlabel="Number of buses [-]",
                legend=:topright, title="Error plot for case study A - Phase 2", ylims = (ymin, ymax))
    elseif which_plot == "error_ph3"
        Plots.scatter([ACR_df_rwlav.n_bus, ACR_df_rwlav.n_bus], [ACR_df_rwlav.err_max_3, ACR_df_rwlav.err_avg_3], markershape=[:circle :utriangle], label=["max. abs. error" "avg. abs. error"], ylabel="Absolute error ϵ [p.u.]", xlabel="Number of buses [-]",
                legend=:topright, title="Error plot for case study A - Phase 3", ylims = (ymin, ymax))
    else
        display("ERROR: plot type $which_plot in argument `which_plot` not recognized. Possibilities are: \"time\", \"error_ph1\", \"error_ph2\", \"error_ph3\"")
    end
end

function plot_result_caseB(path_to_caseB_csv::String, which_plot::String)
    caseb = CSV.read(path_to_caseB_csv)

    LDF_df = caseb |> Query.@filter(_.eq_model .== "LD3F") |> DataFrames.DataFrame
    IVR_df = caseb |> Query.@filter(_.eq_model .== "rIVR") |> DataFrames.DataFrame

    if which_plot == "time"
        Plots.scatter([IVR_df.n_bus, LDF_df.n_bus], [IVR_df.solve_time, LDF_df.solve_time], markershape=[:utriangle :diamond], label=["IVR" "LD3F"], ylabel="Solve time [s]", xlabel="Number of buses [-]",
                legend=:bottomright, title="$which_plot plot for case study B", yscale=:log)
    elseif which_plot == "error_ph1"
        Plots.scatter([LDF_df.n_bus, LDF_df.n_bus], [LDF_df.err_max_1, LDF_df.err_avg_1], markershape=[:circle :utriangle], label=["max. abs. error" "avg. abs. error"], ylabel="Absolute error ϵ [p.u.]", xlabel="Number of buses [-]",
                legend=:topright, title="Error plot for case study B - Phase 1")
    elseif which_plot == "error_ph2"
        Plots.scatter([LDF_df.n_bus, LDF_df.n_bus], [LDF_df.err_max_2, LDF_df.err_avg_2], markershape=[:circle :utriangle], label=["max. abs. error" "avg. abs. error"], ylabel="Absolute error ϵ [p.u.]", xlabel="Number of buses [-]",
                legend=:topright, title="Error plot for case study B - Phase 2")
    elseif which_plot == "error_ph3"
        Plots.scatter([LDF_df.n_bus, LDF_df.n_bus], [LDF_df.err_max_3, LDF_df.err_avg_3], markershape=[:circle :utriangle], label=["max. abs. error" "avg. abs. error"], ylabel="Absolute error ϵ [p.u.]", xlabel="Number of buses [-]",
                legend=:topright, title="Error plot for case study B - Phase 3")
    else
        display("ERROR: plot type $which_plot in argument `which_plot` not recognized. Possibilities are: \"time\", \"error_ph1\", \"error_ph2\", \"error_ph3\"")
    end
end

function plot_result_caseC(path_to_caseC_csv::String)
    casec = CSV.read(path_to_caseC_csv)

    linIVR_df = casec |> Query.@filter(_.eq_model .== "linIVR") |> DataFrames.DataFrame
    rIVR_df = casec |> Query.@filter(_.eq_model .== "rIVR") |> DataFrames.DataFrame

    Plots.scatter([linIVR_df.n_bus, rIVR_df.n_bus], [linIVR_df.solve_time, rIVR_df.solve_time], markershape=[:circle :utriangle], label=["IVR - linear" "IVR - nonlinear"], ylabel="Solve time [s]", xlabel="Number of buses [-]",
             legend=:bottomright, title="Plot for case study C", yscale=:log)
end

function plot_result_caseD(path_to_caseD_csv::String, which_plot::String)

    caseD = CSV.read(path_to_caseD_csv)

    x = (caseD.n_meas.-3)/(3*55)*100

    if which_plot == "error_ph1"
        Plots.scatter([x, x], [caseD.err_max_1, caseD.err_avg_1], markershape=[:circle :utriangle], label=["max. abs. error" "avg. abs. error"], ylabel="Absolute error ϵ [p.u.]", xlabel="Measured users [%]",
                legend=:bottomleft, title="Error plot for case study D - Phase 1", yscale=:log)
    elseif which_plot == "error_ph2"
        Plots.scatter([x, x], [caseD.err_max_2, caseD.err_avg_2], markershape=[:circle :utriangle], label=["max. abs. error" "avg. abs. error"], ylabel="Absolute error ϵ [p.u.]", xlabel="Measured users [%]",
                legend=:bottomleft, title="Error plot for case study D - Phase 2", yscale=:log)
    elseif which_plot == "error_ph3"
        Plots.scatter([x, x], [caseD.err_max_3, caseD.err_avg_3], markershape=[:circle :utriangle], label=["max. abs. error" "avg. abs. error"], ylabel="Absolute error ϵ [p.u.]", xlabel="Measured users [%]",
                legend=:bottomleft, title="Error plot for case study D - Phase 3", yscale=:log)
    else
        display("ERROR: plot type $which_plot in argument `which_plot` not recognized. Possibilities are: \"error_ph1\", \"error_ph2\", \"error_ph3\"")
    end
end


text:
using Revise
using NLS_Fit
using DelimitedFiles
using LinearAlgebra

rootDir = joinpath(dirname(pathof(NLS_Fit)), "..")
dataDir = joinpath(rootDir,"data")
dataFile = joinpath(dataDir,"simple_gaussian.txt")

# Getting started example ****************
#

# Data ================
#
XY=readdlm(dataFile);
X=XY[:,1];
Y=XY[:,2];

# Fitting ================
#
model = Gaussian_Peak()

θ_init = Float64[1,10,5]

nls = NLS_ForwardDiff_From_Model2Fit(model,X,Y);
conf = NLS_Solver.LevenbergMarquardt_Conf()
result = NLS_Solver.solve(nls,θ_init,conf)

# Fit result ================
#
@assert NLS_Solver.converged(result)

θ_true = Float64[2,25,3] # cheating...
θ_solver = NLS_Solver.solution(result)

println("Found $θ_solver, error is $(norm(θ_solver-θ_true,Inf))")
println("(caveat: error is mainly due to noisy data, not due to solver accuracy)")

# Check solver accuracy by using a noise free input ****************
#

# Generate noise free input ================
#
n = 50
model = Gaussian_Peak()
θ = Float64[2,n/2,3]

X=Float64[1:n;];
Y=eval_y(model,X,θ);

# Solve and print new noise free results ================
#
nls = NLS_ForwardDiff_From_Model2Fit(model,X,Y);
result = NLS_Solver.solve(nls,θ_init,conf)
θ_solver = NLS_Solver.solution(result)

println("Found $θ_solver, error is $(norm(θ_solver-θ_true,Inf))")


text:
import LinearAlgebra: I, Hermitian, tr, inv
import ForneyLab: unsafeCov, unsafeMean, unsafePrecision, VariateType,
				  collectNaiveVariationalNodeInbounds, assembleClamp!, ultimatePartner
include("util.jl")

export ruleVariationalNARXOutNPPPPP,
       ruleVariationalNARXIn1PNPPPP,
       ruleVariationalNARXIn2PPNPPP,
       ruleVariationalNARXIn3PPPNPP,
	   ruleVariationalNARXIn4PPPPNP,
	   ruleVariationalNARXIn5PPPPPN

function ruleVariationalNARXOutNPPPPP(g :: Function,
									  marg_y :: Nothing,
									  marg_θ :: ProbabilityDistribution{Multivariate},
									  marg_x :: ProbabilityDistribution{Multivariate},
									  marg_u :: ProbabilityDistribution{Univariate},
                                      marg_z :: ProbabilityDistribution{Multivariate},
                                      marg_τ :: ProbabilityDistribution{Univariate})

	# Extract moments of beliefs
	mθ = unsafeMean(marg_θ)
	mx = unsafeMean(marg_x)
	mz = unsafeMean(marg_z)
	mu = unsafeMean(marg_u)
	mτ = unsafeMean(marg_τ)

	# Evaluate f at mθ
	fθ = mθ'*g([mx; mu; mz])

	# Set outgoing message
	return Message(Univariate, GaussianMeanPrecision, m=fθ, w=mτ)
end

function ruleVariationalNARXIn1PNPPPP(g :: Function,
									  marg_y :: ProbabilityDistribution{Univariate},
                                      marg_θ :: Nothing,
									  marg_x :: ProbabilityDistribution{Multivariate},
									  marg_u :: ProbabilityDistribution{Univariate},
                                      marg_z :: ProbabilityDistribution{Multivariate},
                                      marg_τ :: ProbabilityDistribution{Univariate})

    # Extract moments of beliefs
	my = unsafeMean(marg_y)
	mx = unsafeMean(marg_x)
	mu = unsafeMean(marg_u)
	mz = unsafeMean(marg_z)
	mτ = unsafeMean(marg_τ)

	# Jacobian of f w.r.t. θ
	Jθ = g([mx; mu; mz])

	# Update parameters
	Φ = mτ*Jθ*Jθ'
	ϕ = mτ*my*Jθ

	# Set message
    return Message(Multivariate, GaussianWeightedMeanPrecision, xi=ϕ, w=Φ)
end

function ruleVariationalNARXIn2PPNPPP(g :: Function,
									  marg_y :: ProbabilityDistribution{Univariate},
									  marg_θ :: ProbabilityDistribution{Multivariate},
									  marg_x :: Nothing,
									  marg_u :: ProbabilityDistribution{Univariate},
							  	      marg_z :: ProbabilityDistribution{Multivariate},
                                      marg_τ :: ProbabilityDistribution{Univariate})

	error("Output history vector should be observed.")
    return Message(vague(GaussianWeightedMeanPrecision, 2))
end

function ruleVariationalNARXIn3PPPNPP(g :: Function,
									  marg_y :: ProbabilityDistribution{Univariate},
									  marg_θ :: ProbabilityDistribution{Multivariate},
									  marg_x :: ProbabilityDistribution{Multivariate},
									  marg_u :: Nothing,
							  	      marg_z :: ProbabilityDistribution{Multivariate},
                                      marg_τ :: ProbabilityDistribution{Univariate})

	error("Current input should be observed.")
    return Message(vague(GaussianWeightedMeanPrecision, 2))
end

function ruleVariationalNARXIn4PPPPNP(g :: Function,
									  marg_y :: ProbabilityDistribution{Univariate},
									  marg_θ :: ProbabilityDistribution{Multivariate},
									  marg_x :: ProbabilityDistribution{Multivariate},
									  marg_u :: ProbabilityDistribution{Univariate},
							  	      marg_z :: Nothing,
                                      marg_τ :: ProbabilityDistribution{Univariate})

	error("Input history vector should be observed.")
    return Message(vague(GaussianWeightedMeanPrecision))
end

function ruleVariationalNARXIn5PPPPPN(g :: Function,
									  marg_y :: ProbabilityDistribution{Univariate},
									  marg_θ :: ProbabilityDistribution{Multivariate},
									  marg_x :: ProbabilityDistribution{Multivariate},
									  marg_u :: ProbabilityDistribution{Univariate},
							  	      marg_z :: ProbabilityDistribution{Multivariate},
                                      marg_τ :: Nothing)

    # Extract moments of beliefs
	my = unsafeMean(marg_y)
	mθ = unsafeMean(marg_θ)
	mx = unsafeMean(marg_x)
	mz = unsafeMean(marg_z)
	mu = unsafeMean(marg_u)
	Vθ = unsafeCov(marg_θ)

	# Jacobian of f w.r.t. θ
	Jθ = g([mx; mu; mz])

	# Auto-regression function
	fθ = mθ'*Jθ

	# Update parameters
	a = 3/2.
	b = (my^2 - 2*my*fθ + fθ^2 + Jθ'*Vθ*Jθ)/2.

	# Set message
    return Message(Univariate, Gamma, a=a, b=b)
end

function collectNaiveVariationalNodeInbounds(node::NAutoregressiveX, entry::ScheduleEntry)
	inbounds = Any[]

	# Push function to calling signature (g needs to be defined in user scope)
	push!(inbounds, Dict{Symbol, Any}(:g => node.g, :keyword => false))

    target_to_marginal_entry = currentInferenceAlgorithm().target_to_marginal_entry

    for node_interface in entry.interface.node.interfaces
        inbound_interface = ultimatePartner(node_interface)
        if node_interface === entry.interface
            # Ignore marginal of outbound edge
            push!(inbounds, nothing)
        elseif (inbound_interface != nothing) && isa(inbound_interface.node, Clamp)
            # Hard-code marginal of constant node in schedule
            push!(inbounds, assembleClamp!(inbound_interface.node, ProbabilityDistribution))
        else
            # Collect entry from marginal schedule
            push!(inbounds, target_to_marginal_entry[node_interface.edge.variable])
        end
    end

    return inbounds
end

max_stars_repo_name: wmkouw/NARX

max_issues_repo_name: wmkouw/NARX

max_forks_repo_name: wmkouw/NARX


text:
struct RegInfo
    creg_ssa_map::Dict{Symbol, Int}
    ssa_creg_map::Dict{Int, Symbol}
    creg_size::Dict{Symbol, Int}
    qreg_to_locs::Dict{Int, Vector{Int}}
    locs_to_qreg::Dict{Int, Int}
    locs_to_addr::Dict{Int, Int}
end

RegInfo(ci::CodeInfo) = RegInfo(OpenQASMTarget(), ci)

function RegInfo(target::OpenQASMTarget, ci::CodeInfo)
    creg_ssa_map, ssa_creg_map = extract_creg_map(ci)
    creg_size = extract_creg_size(ci, creg_ssa_map)
    locs_to_qreg = extract_qreg(target, ci)
    qreg_to_locs, locs_to_qreg, locs_to_addr = extract_locs_address(locs_to_qreg)
    RegInfo(creg_ssa_map, ssa_creg_map, creg_size, qreg_to_locs, locs_to_qreg, locs_to_addr)
end

function Base.getindex(ri::RegInfo, loc::Locations{Int})
    plain_loc = plain(loc)
    qreg = ri.locs_to_qreg[plain_loc]
    addr = ri.locs_to_addr[plain_loc]
    if length(ri.qreg_to_locs[qreg]) == 1 && addr == 0
        return Bit(qreg_name(qreg))
    else
        return Bit(qreg_name(qreg), addr)
    end
end

function Base.getindex(ri::RegInfo, loc::CtrlLocations)
    length(loc) == 1 || error("multiple control is not supported by QASM")
    loc.flags[1] || error("inverse control is not supported by QASM")
    return ri[loc.storage]
end

function Base.show(io::IO, ::MIME"text/plain", info::RegInfo)
    println(io, "RegInfo")
    tab = "  "
    isempty(info.creg_size) || printstyled(io, "creg:\n"; color=:light_cyan)
    for (name, size) in info.creg_size
        printstyled(io, tab, name; color=:green)
        printstyled(io, "[", size, "]"; color=:light_blue)
        println(io)
    end

    isempty(info.qreg_to_locs) || printstyled(io, "qreg:\n"; color=:light_cyan)
    for (k, (qreg, locs)) in enumerate(info.qreg_to_locs)
        printstyled(io, tab, "%", qreg; color=:light_magenta)
        printstyled(io, "[", length(locs), "]"; color=:light_blue)
        print(io, " => ", locs)
        if k != length(info.qreg_to_locs)
            println(io)
        end
    end
end

function extract_creg_map(ci::CodeInfo)
    creg_ssa_map = Dict{Symbol, Int}()
    ssa_creg_map = Dict{Int, Symbol}()

    # find alias name for return NamedTuple
    for (v, stmt) in enumerate(ci.code)
        ret = @match stmt begin
            ReturnNode(SSAValue(id)) => ci.code[id]
            _ => nothing
        end

        isnothing(ret) && continue

        @switch ret begin
            @case Expr(:new, type::Type{<:NamedTuple}, args...)
                cregs = type.parameters[1]
                ctypes = type.parameters[2]
                ctypes <: NTuple{N, <:MeasureResult} where N || continue
                for (name, val) in zip(cregs, args)
                    val::SSAValue
                    creg_ssa_map[name] = val.id
                    ssa_creg_map[val.id] = name
                end
            @case _
                nothing
        end
    end
    return creg_ssa_map, ssa_creg_map
end

function extract_creg_size(ci::CodeInfo, creg_ssa_map::Dict{Symbol, Int})
    creg_size = Dict{Symbol, Int}()
    for (key, id) in creg_ssa_map
        stmt = ci.code[id]
        @switch stmt begin
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :measure), _, QuoteNode(locs::Locations))
                creg_size[key] = length(locs)
            @case _
                error("expect measurement statement, got $stmt")
        end
    end
    return creg_size
end

function extract_locs_address(locs_to_qreg::Dict{Int, Int})
    qreg_to_locs = Dict{Int, Vector{Int}}()
    for (k, r) in locs_to_qreg
        locs = get!(qreg_to_locs, r, Int[])
        push!(locs, k)
    end

    # loc => qreg, addr
    locs_to_qreg = Dict{Int, Int}()
    locs_to_addr = Dict{Int, Int}()

    for (r, locs) in qreg_to_locs
        sort!(locs)
        for (k, loc) in enumerate(locs)
            locs_to_qreg[loc] = r
            # QASM starts from 0
            locs_to_addr[loc] = k - 1
        end
    end
    return qreg_to_locs, locs_to_qreg, locs_to_addr
end

function extract_qreg(target::OpenQASMTarget, ci::CodeInfo)
    locs_to_qreg = Dict{Int, Int}()

    for (v, stmt) in enumerate(ci.code)
        @switch stmt begin
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :apply), _, _, QuoteNode(locs::Locations))
                record_locations!(locs_to_qreg, target, locs)
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :apply), _, _, QuoteNode(locs::Locations), QuoteNode(ctrl::CtrlLocations))
                record_locations!(locs_to_qreg, target, locs)
                record_locations!(locs_to_qreg, target, ctrl.storage)
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :measure), _, QuoteNode(locs::Locations))
                allocate_new_qreg!(locs_to_qreg, locs)
                record_locations!(locs_to_qreg, target, locs)
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :barrier), _, QuoteNode(locs::Locations))
                record_locations!(locs_to_qreg, target, locs)
            # skip other statement
            @case _
                nothing
        end
    end
    return locs_to_qreg
end

function record_locations!(locs_to_qreg::Dict{Int, Int}, target::OpenQASMTarget, locs::Locations)
    for loc in locs
        plain_loc = plain(loc)
        if target.toplevel
            get!(locs_to_qreg, plain_loc, 1)
        else
            # force allocate new register for qasm gate
            if !haskey(locs_to_qreg, plain_loc)
                locs_to_qreg[plain_loc] = length(keys(locs_to_qreg)) + 1
            end
        end
    end
    return locs_to_qreg
end

function allocate_new_qreg!(locs_to_qreg::Dict{Int, Int}, locs::Locations)
    max_qreg_addr = maximum(locs) do loc
        # we will add 1 later
        get(locs_to_qreg, plain(loc), 0)
    end

    for loc in locs
        locs_to_qreg[plain(loc)] = max_qreg_addr + 1
    end
    return locs_to_qreg
end


text:
@testset "SampleArray" begin
    function cmparrays(a, b;
        rate_=rate(a), domain_=domain(a), nchannels_=nchannels(a),
        nframes_=nframes(a), names_=names(a), eltype_=eltype(a), typeof_=typeof(a), data_=data(a))
        isnothing(rate_) || rate(b) == rate_ || error("rate $(rate(b)) != $(rate_)")
        isnothing(domain_) || domain(b) == domain_ || error("domain $(domain(b)) != $(domain_)")
        isnothing(nchannels_) || nchannels(b) == nchannels_ || error("nchannels $(nchannels(b)) != $(nchannels_)")
        isnothing(nframes_) || nframes(b) == nframes_ || error("nframes $(nframes(b)) != $(nframes_)")
        isnothing(names_) || names(b) == names_ || error("names $(names(b)) != $(names_)")
        isnothing(eltype_) || eltype(b) == eltype_ || error("eltype $(eltype(b)) != $(eltype_)")
        isnothing(typeof_) || typeof(b) == typeof_ || error("typeof $(typeof(b)) != $(typeof_)")
        isnothing(data_) || data(b) == data_ || error("data not same!")
        true
    end

    sa = SampleArray(sin.(0:0.1:4pi), 2, [:center]) # from vector
    sb = SampleArray(hcat(cos.(0:0.1:4pi), -cos.(0:0.1:4pi)), 2Hz, [:left, :Right]) # from matrix
    sc = copy(sb)
    sd = SampleArray(hcat(cos.(0:0.1:4pi), -cos.(0:0.1:4pi)), 2Hz)
    se = SampleArray([float(j) for i in 1:100, j in 1:8], 44100Hz, [
            :front_left, :front_right, :rear_left, :rear_right,
            :front_center, :lfe, :side_left, :side_right])

    @test_throws DimensionMismatch SampleArray(rand(3, 2), 1Hz, [:center])
    @test_throws ArgumentError SampleArray(rand(3, 2), 1Hz, [:center, :center])

    @test cmparrays(sb, sc)
    @test sc !== sb
    sd = SampleArray(hcat(cos.(0:0.1:4pi), -cos.(0:0.1:4pi)), 2Hz)

    @testset "basic functions" begin
        @test domain(sa) == domain(sb)
        @test nchannels(sa) == 1
        @test nchannels(sb) == 2
        @test nchannels(se) == 8
        @test nframes(se) == 100
        @test cmparrays(sa, sb; nchannels_=2, names_=[:left, :Right], data_=nothing)
        @test getduration(sa) == 63s
        @test size(sb) == (126, 2)
        @test eltype(sa) == Float64

        @test names(sa)[1] == :center
        @test names(sd) == [Symbol(1), Symbol(2)]
        @test names(sb) == [:left, :Right]
        @test names(se) == [:front_left, :front_right, :rear_left, :rear_right,
            :front_center, :lfe, :side_left, :side_right]

        sb2 = copy(sb)
        @test_throws ArgumentError names!(sb2, :Right, 1) # non unique name
        names!(sb2, :Left, 1)
        @test cmparrays(sb, sb2; names_=[:Left, :Right])
        names!(sb2, [:L, :R], 1:2)
        @test cmparrays(sb, sb2; names_=[:L, :R])
        names!(sb2, [:left, :right])
        @test_throws ArgumentError names!(sb2, [:same, :same]) # non unique name
        @test cmparrays(sb, sb2; names_=[:left, :right])
        names!(sc, [:left, :right])
        @test cmparrays(sb2, sc)
        names!(sb2, [:RIGHT, :LEFT], [:right, :left])
        @test cmparrays(sb, sb2; names_=[:LEFT, :RIGHT])
    end

    @testset "equality" begin
        dta = rand(3, 2)
        sf = SampleArray(dta, 1Hz)
        sg = copy(sf)
        sh = SampleArray(dta, 2Hz)
        @test sf == sg
        @test sf != sh && cmparrays(sf, sh; rate_=nothing, domain_=nothing)
        names!(sg, [:a, :b])
        @test sf != sg && cmparrays(sf, sg; names_=nothing)
    end

    @testset "indexing" begin
        # linear index
        @test se[1] == data(se)[1]
        @test typeof(se[1]) == eltype(se)

        # frame addressing, select all channels
        @test cmparrays(se, se[1:10]; domain_=nothing, nframes_=10, data_=data(se)[1:10, :])

        # channel addressing
        @test cmparrays(se[1:10, :], se[1:10])

        @test cmparrays(se, se[:, 2]; nchannels_=1, names_=[:front_right], data_=data(se)[:, 2:2])
        @test cmparrays(se[:, 2], se[:, 2:2])
        @test cmparrays(se[:, [3, 4, 5, 6, 7]], se[:, 3:7])
        @test cmparrays(se[:, 2], se[:, [2]])
        @test cmparrays(se[:, 2], se[:, :front_right])
        @test cmparrays(se[:, :front_right], se[:, [:front_right]])
        @test_throws ArgumentError se[:, :wrong_channel_name]
        @test cmparrays(se[:, :front_right], se[:, [:front_right]])
        @test cmparrays(se[:, :rear_left..:side_left], se[:, 3:7])
        @test cmparrays(se[:, [name for name in names(se) if occursin("left", String(name))]], se[:, [1, 3, 7]])

        # getting integer indices
        @test [SampleArrays.toindex(sa, t) for t in [0s, 0.8s, 1s, 10s]] == [1, 3, 3, 21]
        @test [SampleArrays.toindex(sa, t) for t in [0, 0.8, 1, 10]] == [1, 3, 3, 21]

        @test [SampleArrays.toindexdelta(sa, t) for t in [1s, 2s, 5s]] == [2, 4, 10]
        @test [SampleArrays.toindexdelta(sa, t) for t in [1.0, 2, 5]] == [2, 4, 10]
        @test [SampleArrays.toindexdelta(sa, t) for t in [1, 2, 5]] == [1, 2, 5]

        @test_throws ArgumentError SampleArrays.toframeidx(sa, 0Hz..1Hz) == 1:3
        @test SampleArrays.toframeidx(sa, 0s..1s) == 1:3
        @test SampleArrays.toframeidx(sa, 0s..0.74s) == 1:2
        @test SampleArrays.toframeidx(sa, 0s..0.75s) == 1:3

        # frame adressing by time
        @test cmparrays(sa[0s..0.75s, 1:1], sa[SampleArrays.toframeidx(sa, 0s..0.75s), 1:1])
        @test cmparrays(sa[0s..0.75s, 1], sa[0s..0.75s, 1:1])
        @test cmparrays(sb[0s..0.75s], sb[0s..0.75s, :])

        # Array indexing
        sb2 = sb[[1,3,5], [1]]
        @test data(sb2) == data(sb)[[1,3,5], [1]]

        # Bool[] and BitArray indexing
        sb2 = sb[data(sa)[:] .> 0, [true, false]]
        @test data(sb2) == data(sb)[data(sa)[:] .> 0, [true, false]]

        # zero channels
        sb2 = sb[:, []]
        @test nchannels(sb2) == 0
        sb3 = sb[:, [false, false]]
        @test nchannels(sb3) == 0
        @test cmparrays(sb2, sb3)

        # setindex!
        sc[1] = -2
        @test data(sc)[1] == -2

        @test_throws ArgumentError sc[1:5] = -2
        sc[1:5] .= -2
        @test all(data(sc)[1:5, :] .== -2)

        sc[1:10,:] .= -5
        @test all(data(sc)[1:10, 1:2] .== -5)

        sc[0s..0.75s, 1:2] .= -13
        @test all(data(sc[0s..0.75s, 1:2]) .== -13)

        sc[0s..0.75s, :left] .= -14
        @test all(data(sc[0s..0.75s, 1:1]) .== -14)

        sc[0s..0.75s, [:left, :right]] .= -15
        @test all(data(sc[0s..0.75s, 1:2]) .== -15)

        sc[[1,3,5], [1]] .= -20
        @test all(data(sc)[[1,3,5], [1]] .== -20)

        mask = data(sc)[:, 1] .> 0
        sc[mask, [false, true]] .= -50
        @test all(data(sc)[mask, [false, true]] .== -50)
    end

    @testset "similar" begin
        @test cmparrays(sb, similar(sb); data_=nothing)
        @test cmparrays(sb, similar(sb, Float32); data_=nothing, eltype_=Float32, typeof_=SampleArray{Float32})
        @test cmparrays(sb, similar(sb, 10, 2); domain_=nothing, nframes_=10, data_=nothing)
        @test cmparrays(sb, similar(sb, Float32, 10, 2); domain_=nothing, nframes_=10,
            data_=nothing, eltype_=Float32, typeof_=SampleArray{Float32})
        @test cmparrays(sb, similar(sb, 10, 1); domain_=nothing, nframes_=10, data_=nothing, nchannels_=1, names_=[:left])
        @test cmparrays(sb, similar(sb, 10, 3); domain_=nothing, nframes_=10, data_=nothing,
            nchannels_=3, names_=SampleArrays._default_channel_names(3))
    end

    @testset "view" begin
        sc = copy(sb)
        vsc = view(sc, 1:100, 1:1)
        typeof(vsc) <: SubArray

        vsc[:] .= -5
        @test all(sc[1:100, 1:1] .== -5)
    end

    @testset "broadcasting" begin
        @test all((sa.^2 + sb[:, 1:1].^2) .≈ 1)
        @test sa + sa == 2*sa
        @test data(sa .+ 5) == data(sa) .+ 5

        sa32 = Float32.(sa)
        ssum = sa + sa32
        @test ssum ≈ 2sa

        @test cmparrays(ssum, sa32 .+ sb; nchannels_=2, names_=nothing, data_=nothing)
    end

    @testset "resampling" begin
        @test all(abs.(resample(resample(sa, 4Hz), 2) .- sa) .< 0.01) # should not this be more precise?
    end

    @testset "concatenation" begin
        @test data(hcat(sa, sb)) == hcat(data(sa), data(sb))
        sc2 = hcat(sc, sc)
        @test cmparrays(sc, sc2[:, 1:2])
        @test cmparrays(sc, sc2[:, 3:4]; names_=[:left_2, :Right_2])
        sc2 = hcat(sc, Float32.(sc))
        @test cmparrays(sc, sc2[:, 1:2])
        @test cmparrays(sc, sc2[:, 3:4]; names_=[:left_2, :Right_2], data_=nothing)
        @test_throws ArgumentError hcat(sa, se) # different rate
        sc2 = vcat(sc, Float32.(sc))
        @test data(sc2) ≈ vcat(data(sc), data(sc))
        @test_throws ArgumentError vcat(sa, sc)
        @test_throws ErrorException cat(sa, sa; dims=1)
    end
end


text:
"mask with all false columns should produce an error"
let
    mask = transpose([false false false; true false true; false true false])
    @test_throws ErrorException MDPs.Examples.random(size(mask, 1), 2, mask)
end

"mask sizes that do not match the states and actions arguments should throw error"
let
    mask = transpose([true false false; true false true; false true false])
    @test_throws ErrorException MDPs.Examples.random(size(mask, 1) + 1, 2, mask)
end

"random transition with a 2d mask argument should have zeros everywhere the mask is false"
let
    mask = transpose([true false false; true false true; false true false])
    transition, reward = MDPs.Examples.random(3, 2, mask)
    @test (transition[:, :, 2] .!= 0.0) == mask
    @test (transition[:, :, 2] .!= 0.0) == mask
end

"random transition with a 3d mask argument should have zeros everywhere the mask is false"
let
    mask = [x==y for x in 1:3, y in 1:3, z in 1:2]
    transition, reward = MDPs.Examples.random(3, 2, mask)
    @test (transition .!= 0.0) == mask
end

"random should accept a Nullable{Array{Bool}} as a mask and the "
"returned array types should be Float64"
let
    transition, reward = MDPs.Examples.random(
        10,
        2,
        Nullable{Array{Bool,2}}(),
    )
    @test typeof(transition) == Array{Float64,3}
    @test typeof(reward) == Array{Float64,2}
end

"random with just two numbers should word and return Float64 arrays"
let
    transition, reward = MDPs.Examples.random(10, 2)
    @test typeof(transition) == Array{Float64,3}
    @test typeof(reward) == Array{Float64,2}
end

"random with different types for transition and reward should return arrays of those types"
let
    transition, reward = MDPs.Examples.random(Float64, Float32, 10, 2)
    @test typeof(transition) == Array{Float64,3}
    @test typeof(reward) == Array{Float32,2}
end

"random should allow Integer reward types"
let
    transition, reward = MDPs.Examples.random(Float64, Int8, 10, 2)
    @test typeof(transition) == Array{Float64,3}
    @test typeof(reward) == Array{Int8,2}
end


text:
using ITensors
using ITensorInfiniteMPS

##############################################################################
# VUMPS parameters
#

maxdim = 64 # Maximum bond dimension
cutoff = 1e-6 # Singular value cutoff when increasing the bond dimension
max_vumps_iters = 100 # Maximum number of iterations of the VUMPS algorithm at each bond dimension
vumps_tol = 1e-6
outer_iters = 10 # Number of times to increase the bond dimension
eager = true

##############################################################################
# CODE BELOW HERE DOES NOT NEED TO BE MODIFIED
#
N = 3# Number of sites in the unit cell
J = -1.0
J₂ = -0.2
h = 1.0;

function space_shifted(q̃sz)
  return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1]
end

space_ = fill(space_shifted(1), N);
s = infsiteinds("S=1/2", N; space=space_)
initstate(n) = "↑"
ψ = InfMPS(s, initstate);

model = Model("ising_extended");
H = InfiniteSum{MPO}(model, s; J=J, J₂=J₂, h=h);

@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...));

vumps_kwargs = (tol=vumps_tol, maxiter=max_vumps_iters, eager)
subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
ψ_0 = @time vumps(H, ψ; vumps_kwargs...)

@time for j in 1:outer_iters
  println("\nIncrease bond dimension")
  ψ_1 = @time subspace_expansion(ψ_0, H; subspace_expansion_kwargs...)
  println("Run VUMPS with new bond dimension")
  global ψ_0 = @time vumps(H, ψ_1; vumps_kwargs...)
end

Sz = [expect(ψ_0, "Sz", n) for n in 1:N]
energy_infinite = expect(ψ_0, H)

Nfinite = 100
sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
Hfinite = MPO(model, sfinite; J=J, J₂=J₂, h=h)
ψfinite = randomMPS(sfinite, initstate; linkdims=10)
@show flux(ψfinite)
sweeps = Sweeps(10)
setmaxdim!(sweeps, maxdim)
setcutoff!(sweeps, cutoff)
energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps)

nfinite = Nfinite ÷ 2
Sz_finite = expect(ψfinite, "Sz")[nfinite:(nfinite + N - 1)]

@show energy_infinite
@show energy_finite_total / Nfinite
@show reference(model, Observable("energy"); J=J, h=h, J₂=J₂)

@show Sz
@show Sz_finite


text:
#
#	Types
#
mutable struct Cloud
	N::Vector{Int64}
	points::Array{Float64,2}
end

mutable struct SubCloud
	cl::Cloud
	N::Vector{Int64}
	verts::Vector{Int64}
end

# complex scalar field
mutable struct Cscalar
	cl::Cloud
	values::Vector{ComplexF64}
end

mutable struct SubCscalar
	cs::Cscalar
	N::Vector{Int64}
	verts::Vector{Int64}
end

#
#	Functions
#

function toint(A::Vector{Int}, N::Vector{Int})
	if length(N)==1
		return A[1]+1
	else
		return A[end]*prod(N[1:end-1]) + toint(A[1:end-1],N[1:end-1])
	end
end

function tovec(i::Int, N::Vector{Int})
	if length(N)==1
		return [rem(i-1,N[1])]
	else
		return vcat(tovec(i,N[1:end-1]),[rem(div(i-1,prod(N[1:end-1])),N[end])])
	end
end

function getverts(cl::Cloud, axis::Integer)

	axis <= dim(cl) || throw(ArgumentError("your cloud is $(dim(cl))-dimensional but you are asking to slice along axis $(axis)"))

	A = div.(cl.N,2)
	A[axis] = 0
	verts = [toint(A,cl.N)]

	for i=1:cl.N[axis]-1
		A[axis] += 1
		verts = vcat(verts,toint(A,cl.N))
	end

	N = ones(Int, dim(cl))
	N[axis] = cl.N[axis]

	return N, verts
end

function slice(cl::Cloud, axis::Integer)

	N, verts = getverts(cl,axis)

	SubCloud(cl,N,verts)
end

function slice(cs::Cscalar, axis::Integer)

	N, verts = getverts(cs.cl, axis)

	SubCscalar(cs, N, verts)
end

# Constructors
function Cloud(
	d::Integer,
	N::Vector{T} where T<:Integer)

	Cloud(N,Array{Float64,2}(undef,Npoints(N),d))
end

function Cscalar(cl::Cloud)
	Cscalar(cl, Vector{ComplexF64}(undef,Npoints(cl)))
end

dim(N::Vector{T} where T<: Integer) = length(N)
dim(cl::Cloud)   	= dim(cl.N)
dim(scl::SubCloud) 	= dim(scl.cl.N)
dim(cs::Cscalar) 	= dim(cs.cl.N)
dim(scs::SubCscalar)= dim(scs.cs.cl.N)

Npoints(N::Vector{T} where T<:Integer) = prod(N)
Npoints(cl::Cloud)		= Npoints(cl.N)
Npoints(cs::Cscalar)	= Npoints(cs.cl)
Npoints(scl::SubCloud)	= Npoints(scl.N)
Npoints(scs::SubCscalar)= Npoints(scs.N)

function init!(
	cl::Cloud,
	L::Vector{T} where T<:AbstractFloat)

	for i=1:Npoints(cl)
		A = tovec(i, cl.N)
		@. cl.points[i,1:end] = L*(-1 + 2*A/cl.N)
	end
end

# sample a function, f, on a cloud as a complex scalar field
function sample!(cs::Cscalar, f::Function)
	for i=1:Npoints(cs.cl)
		cs.values[i] = f(cs.cl.points[i,1:end])
	end
end

# printers
function tovtk(cl::Cloud, filename::AbstractString)
	f = open(filename, "w")
	writepreamble(f, "Point cloud")
	writepoints(f, dim(cl), Npoints(cl), cl.points)
	writeverts(f,Npoints(cl))
	close(f)
end

function tovtk(scl::SubCloud, filename::AbstractString)
	f = open(filename, "w")
	writepreamble(f, "Point (sub) cloud")
	writepoints(f, dim(scl), Npoints(scl), scl.cl.points[scl.verts,1:end])
	writeverts(f,Npoints(scl))
	close(f)
end

function tovtk(cs::Cscalar, filename::AbstractString)
	f = open(filename, "w")
	writepreamble(f, "Complex scalar field")
	writepoints(f, dim(cs), Npoints(cs.cl), cs.cl.points)
	writeverts(f,Npoints(cs.cl))
	print(f, "CELL_DATA "*string(Npoints(cs.cl))*"\n")
	writescalars(f,"real",Npoints(cs.cl),real.(cs.values))
	writescalars(f,"imag",Npoints(cs.cl),imag.(cs.values))
	close(f)
end

function tovtk(scs::SubCscalar, filename::AbstractString)
	f = open(filename, "w")
	writepreamble(f, "Complex (sub) scalar field")
	writepoints(f, dim(scs), Npoints(scs), scs.cs.cl.points[scs.verts,1:end])
	writeverts(f,Npoints(scs))
	print(f, "CELL_DATA "*string(Npoints(scs))*"\n")
	writescalars(f,"real",Npoints(scs),real.(scs.cs.values[scs.verts]))
	writescalars(f,"imag",Npoints(scs),imag.(scs.cs.values[scs.verts]))
	close(f)
end

function writepreamble(f::IO, label::AbstractString)
	preamble = """
	# vtk DataFile Version 4.0
	"""*label*"""\n
	ASCII
	DATASET UNSTRUCTURED_GRID
	"""
	print(f,preamble)
end

function writepoints(
	f::IO,
	DIMENSION::Integer,
	N::Integer,
	points::Array{Float64,2})

	print(f, "POINTS $N float\n")
	for i=1:N
        for j=1:DIMENSION
            print(f, "$(points[i,j]) ")
        end
        for j=DIMENSION+1:3
            print(f, "0.0 ")
        end
        print(f, "\n")
    end
end

function writeverts(f::IO, N::Integer)
    print(f,"CELLS $N $(2N)\n")
    for i=1:N
        print(f, "1 $(i-1)\n")
    end
    print(f,"CELL_TYPES $N\n")
    for i=1:N
        print(f, "1\n")
    end
end

function writescalars(
        f::IO,
        name::AbstractString,
        N::Integer,
        scalars::Vector{Float64}
    )
    print(f, "SCALARS "*name*" float\nLOOKUP_TABLE default\n")
    for i=1:N
        print(f, "$(scalars[i])\n")
    end
end

function tocsv(scs::SubCscalar, filename)
	f = open(filename, "w")
	if dim(scs)==3
		print(f, "x,y,z,Re,Im\n")
		for i=1:Npoints(scs)
			print(f,
				"$(scs.cs.cl.points[scs.verts[i],1]),"*
				"$(scs.cs.cl.points[scs.verts[i],2]),"*
				"$(scs.cs.cl.points[scs.verts[i],3]),"*
				"$(real(scs.cs.values[scs.verts[i]])),"*
				"$(imag(scs.cs.values[scs.verts[i]]))\n")
		end
	elseif dim(scs)==2
		print(f, "x,y,Re,Im\n")
		for i=1:Npoints(scs)
			print(f,
				"$(scs.cs.cl.points[scs.verts[i],1]),"*
				"$(scs.cs.cl.points[scs.verts[i],2]),"*
				"$(real(scs.cs.values[scs.verts[i]])),"*
				"$(imag(scs.cs.values[scs.verts[i]]))\n")
		end
	elseif dim(scs)==1
		print(f, "x,Re,Im\n")
		for i=1:Npoints(scs)
			print(f,
				"$(scs.cs.cl.points[scs.verts[i],1]),"*
				"$(real(scs.cs.values[scs.verts[i]])),"*
				"$(imag(scs.cs.values[scs.verts[i]]))\n")
		end
	else
		throw(ArgumentError("invalid number of dimensions specified."))
	end
	close(f)
end

# warning: this function creates a big memory block;
# not yet optimised for performance.
function surface(cs::Cscalar)

	dim(cs) == 2 || throw(ArgumentError("input Cscalar needs to be dimension 2"))

	N = [cs.cl.N[1], cs.cl.N[2], 1]

	# allocate memory
	points = Array{Float64,2}(undef, Npoints(N), dim(N))

	@. points[1:end,1:2] = cs.cl.points[1:end,1:2]
	@. points[1:end,3]   = real(cs.values[1:end])

	Cloud(N,points)
end


text:
"""
    BaB(optimizer, ϵ::Float64)

BaB uses branch and bound to estimate the range of the output node.

# Problem requirement
1. Network: any depth, ReLU activation, single output
2. Input: hyperrectangle
3. Output: hyperrectangle (1d interval)

# Return
`CounterExampleResult` or `ReachabilityResult`

# Method
Branch and bound.
For branch, it uses iterative interval refinement.
For bound, it computes concrete bounds by sampling, approximated bound by optimization.
- `optimizer` default `GLPKSolverMIP()`
- `ϵ` is the desired accurancy for termination, default `0.1`.

# Property
Sound and complete.

# Reference
[R. Bunel, I. Turkaslan, P. H. Torr, P. Kohli, and M. P. Kumar,
"A Unified View of Piecewise Linear Neural Network Verification,"
*ArXiv Preprint ArXiv:1711.00455*, 2017.](https://arxiv.org/abs/1711.00455)
"""
@with_kw struct BaB <: Solver
    optimizer = GLPK.Optimizer
    ϵ::Float64 = 0.1
end

function solve(solver::BaB, problem::Problem)
    (u_approx, u, x_u) = output_bound(solver, problem, :max)
    (l_approx, l, x_l) = output_bound(solver, problem, :min)
    bound = Hyperrectangle(low = [l], high = [u])
    reach = Hyperrectangle(low = [l_approx], high = [u_approx])

    output = problem.output
    if reach ⊆ output
        return ReachabilityResult(:holds, [reach])
    end
    high(bound) > high(output) && return CounterExampleResult(:violated, x_u)
    low(bound)  < low(output)  && return CounterExampleResult(:violated, x_l)
    return ReachabilityResult(:unknown, [reach])
end

function output_bound(solver::BaB, problem::Problem, type::Symbol)
    nnet = problem.network
    global_concrete, x_star = concrete_bound(nnet, problem.input, type)
    global_approx = approx_bound(nnet, problem.input, solver.optimizer, type)
    doms = Tuple{Float64, Hyperrectangle}[(global_approx, problem.input)]
    index = ifelse(type == :max, 1, -1)
    while index * (global_approx - global_concrete) > solver.ϵ
        dom, doms = pick_out(doms) # pick_out implements the search strategy
        subdoms = split_dom(dom[2]) # split implements the branching rule
        for i in 1:length(subdoms)
            dom_concrete, x = concrete_bound(nnet, subdoms[i], type) # Sample
            dom_approx = approx_bound(nnet, subdoms[i], solver.optimizer, type)
            if index * (dom_concrete - global_concrete) > 0
                (global_concrete, x_star) = (dom_concrete, x)
            end
            if index * (dom_approx - global_concrete) > 0
                add_domain!(doms, (dom_approx, subdoms[i]), type)
            end
        end
        length(doms) == 0 && return (global_approx, global_concrete, x_star)
        global_approx = doms[1][1]
    end
    return (global_approx, global_concrete, x_star)
end

# Always pick the first dom
function pick_out(doms)
    return (doms[1], doms[2:end])
end

function add_domain!(doms::Vector{Tuple{Float64, Hyperrectangle}}, new::Tuple{Float64, Hyperrectangle}, type::Symbol)
    rank = length(doms) + 1
    index = ifelse(type == :max, 1, -1)
    for i in 1:length(doms)
        if index * (new[1] - doms[i][1]) >= 0
            rank = i
            break
        end
    end
    insert!(doms, rank, new)
end

# Always split the longest input dimension
function split_dom(dom::Hyperrectangle)
    max_value, index_to_split = findmax(dom.radius)
    return split_interval(dom, index_to_split)
end

# For simplicity
function concrete_bound(nnet::Network, subdom::Hyperrectangle, type::Symbol)
    points = [center(subdom), low(subdom), high(subdom)]
    vals = [sum(compute_output(nnet, p)) for p in points]
    i = (type == :min) ? argmin(vals) : argmax(vals)
    return vals[i], points[i]
end

function approx_bound(nnet::Network, dom::Hyperrectangle, optimizer, type::Symbol)
    model = Model(optimizer)
    model[:bounds] = get_bounds(nnet, dom, false)
    model[:before_act] = true
    z = init_vars(model, nnet, :z, with_input=true)

    add_set_constraint!(model, dom, first(z))
    encode_network!(model, nnet, TriangularRelaxedLP())
    index = ifelse(type == :max, 1, -1)
    o = sum(last(z))
    @objective(model, Max, index * o)
    optimize!(model)
    termination_status(model) == OPTIMAL && return value(o)
    error("Could not find bound for dom: ", dom)
end


text:
export MeshPartitionTree
export tree, leaf_meshes
# A data structure to hold a hierarchical partition of a mesh.
# Since the mesh is partitioned, we only need to store the leaves of the partition
# tree to reconstruct the mesh.
#
# Each node has the name of the partition and an ID.
# If the node is an internal node, it has id = 0.
# If the node is a leaf node, it has a non-zero id, corresponding to the index by
# which the mesh may be found in the MeshPartitionTree type's leaf_meshes.
struct MeshPartitionTree{M <: AbstractMesh}
    partition_tree::Tree{Tuple{Int64,String}}
    leaf_meshes::Vector{M}
end

tree(mpt::MeshPartitionTree) = mpt.partition_tree
leaf_meshes(mpt::MeshPartitionTree) = mpt.leaf_meshes

# Partitions a mesh based upon the names of its groups
# For a hierarchical partition, all partitions in the hierarchy must contain face sets
# of the form "<name>_L<N>" where name is a string, and N is an integer.
# N is the level of the node in the tree: 1, 2, 3, etc...
# Example: "Grid_L1_triangle", or "Partition_L3"
function MeshPartitionTree(mesh::AbstractMesh; by::String="MPACT")
    @info "Partitioning mesh: "*name(mesh)
    # Extract the names of all face sets that contain 'by' (the variable)
    partition_names = _get_partition_names(mesh, by)
    # Create a tree to store the partition hierarchy.
    root = _create_mesh_partition_tree(mesh, partition_names)
    # Construct the leaf meshes
    leaf_meshes = _create_leaf_meshes(mesh, root)
    return MeshPartitionTree(root, leaf_meshes)
end

function _create_leaf_meshes(mesh::AbstractMesh, root::Tree)
    leaf_nodes = sort!(leaves(root), by=x->x.data)
    leaf_meshes = Vector{typeof(mesh)}(undef, length(leaf_nodes))
    leaf_ctr = 1
    node_children = children(root)
    if !isnothing(node_children)
        for child ∈ sort(node_children, by=x->x.data)
            child_mesh = submesh(mesh, data(child)[2])
            leaf_ctr = _create_leaf_meshes!(child_mesh, child, leaf_meshes, leaf_ctr)
        end
    else
        # Remove any Lattice, Module, or Coarse_Cell groups, since this info should now
        # be encoded in the tree
        mpact_groups = filter(x->isa_MPACT_partition_name(x), keys(groups(mesh)))
        for grp in mpact_groups
            pop!(groups(mesh), grp)
        end
        name = data(root)[2]
        root.data = (leaf_ctr, name)
        leaf_meshes[leaf_ctr] = mesh
    end
    return leaf_meshes
end

function _create_leaf_meshes!(mesh::AbstractMesh,
                              node::Tree,
                              leaf_meshes,
                              leaf_ctr::Int64=1
                             )
    node_children = children(node)
    if !isnothing(node_children)
        for child ∈ sort(node_children, by=x->x.data)
            child_mesh = submesh(mesh, data(child)[2])
            leaf_ctr = _create_leaf_meshes!(child_mesh, child, leaf_meshes, leaf_ctr)
        end
    else
        # Remove any Lattice, Module, or Coarse_Cell groups, since this info should now
        # be encoded in the tree
        mpact_groups = filter(x->isa_MPACT_partition_name(x), keys(groups(mesh)))
        for grp in mpact_groups
            pop!(mesh.groups, grp)
        end
        leaf_meshes[leaf_ctr] = mesh
        name = data(node)[2]
        node.data = (leaf_ctr, name)
        leaf_ctr += 1
    end
    return leaf_ctr
end

# Create a tree to store grid relationships.
function _create_mesh_partition_tree(mesh::AbstractMesh, partition_names::Vector{String})
    root = Tree((0, name(mesh)))
    parents = [root]
    tree_type = typeof(root)
    next_parents = tree_type[]
    remaining_names = copy(partition_names)
    # Extract the sets that are not a subset of any other set.
    this_level = Int64[]
    mesh_groups = groups(mesh)
    while length(remaining_names) > 0
        for i in eachindex(remaining_names)
            name_i = remaining_names[i]
            i_isa_lattice = startswith(name_i, "Lattice")
            i_isa_module = startswith(name_i, "Module")
            isa_subset = false
            for j in eachindex(remaining_names)
                if i === j
                    continue
                end
                name_j = remaining_names[j]
                j_isa_coarsecell = startswith(name_j, "Coarse")
                if j_isa_coarsecell && (i_isa_module || i_isa_lattice)
                    continue
                end
                j_isa_module = startswith(name_j, "Module")
                if j_isa_module && i_isa_lattice
                    continue
                end
                if mesh_groups[name_i] ⊆ mesh_groups[name_j]
                    isa_subset = true
                    break
                end
            end
            if !isa_subset
                push!(this_level, i)
            end
        end
        # Add the groups that are not a subset to the tree
        # and next_parents
        for id in this_level
            name = remaining_names[id]
            for parent in parents
                if isroot(parent)
                    node = Tree((0, name), parent)
                    push!(next_parents, node)
                else
                    parent_name = data(parent)[2]
                    if mesh_groups[name] ⊆ mesh_groups[parent_name]
                        node = Tree((0, name), parent)
                        push!(next_parents, node)
                        break
                    end
                end
            end
        end
        parents = next_parents
        next_parents = tree_type[]
        # Remove the groups added to the tree from the remaining names
        deleteat!(remaining_names, this_level)
        this_level = Int64[]
    end
    return root
end

function isa_MPACT_partition_name(x::String)
    return startswith(x, "Coarse_Cell_(") ||
           startswith(x, "Module_("     ) ||
           startswith(x, "Lattice_("    )
end
# Extract partition names
function _get_partition_names(mesh::AbstractMesh, by::String)
    # If by === MPACT, then we will partition by Lattice, Module, and Coarse Cell
    # Otherwise, partition by the string contained in by
    partition_names = collect(keys(groups(mesh)))
    if length(partition_names) === 0
        error("The mesh does not have any groups.")
    end
    if by === "MPACT"
        filter!(x->isa_MPACT_partition_name(x), partition_names)
        if length(partition_names) === 0
            error("The mesh does not have any MPACT grid hierarchy groups.")
        end
    else
        filter!(x->occursin(by, x), partition_names)
        if length(partition_names) === 0
            error("The mesh does not have any groups containing '", by, "'.")
        end
    end
    return sort!(partition_names)
end

function Base.show(io::IO, mpt::MeshPartitionTree{M}) where {M}
    println("MeshPartitionTree{",M,"}")
    println(tree(mpt))
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
# Model
include_model("biped_pinned")
model = free_time_model(additive_noise_model(model))

# # Visualize
# include(joinpath(pwd(), "models/visualize.jl"))
# vis = Visualizer()
# open(vis)
#
# urdf = joinpath(pwd(), "models/biped/urdf/biped_left_pinned.urdf")
# mechanism = parse_urdf(urdf, floating=false)
# mvis = MechanismVisualizer(mechanism,
#     URDFVisuals(urdf, package_path=[dirname(dirname(urdf))]), vis)

ϵ = 1.0e-8
θ = 12.5 * pi / 180
h = model.l2 + model.l1 * cos(θ)
ψ = acos(h / (model.l1 + model.l2))
stride = sin(θ) * model.l1 + sin(ψ) * (model.l1 + model.l2)
x1 = [π - θ, π + ψ, θ, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
xT = [π + ψ, π - θ - ϵ, 0.0, θ, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
kinematics(model, x1)[1]
kinematics(model, xT)[1] * 2.0

kinematics(model, x1)[2]
kinematics(model, xT)[2]

q1 = transformation_to_urdf_left_pinned(model, x1[1:5])
# set_configuration!(mvis, q1)

qT = transformation_to_urdf_left_pinned(model, xT[1:5])
# set_configuration!(mvis, qT)

# Horizon
T = 21

tf0 = 2.0
h0 = tf0 / (T-1)

# Bounds
ul, uu = control_bounds(model, T,
	[-5.0 * ones(model.m - 1); 0.0 * h0],
	[5.0 * ones(model.m - 1); h0])

xl, xu = state_bounds(model, T, x1 = [x1[1:5]; Inf * ones(5)])

# Objective
include_objective("nonlinear_stage")
qq = 1.0 * [1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0]
Q = [Diagonal(qq) for t = 1:T]
R = [Diagonal([1.0e-3 * ones(model.m - 1); h0]) for t = 1:T-1]

obj_track = quadratic_time_tracking_objective(
		Q,
		R,
    	[xT for t = 1:T],
		[zeros(model.m) for t = 1:T-1],
		1.0)

l_stage_fh(x, u, t) = 100.0 * (kinematics(model, view(x, 1:5))[2] - 0.25)^2.0
l_terminal_fh(x) = 0.0
obj_fh = nonlinear_stage_objective(l_stage_fh, l_terminal_fh)

obj_multi = MultiObjective([obj_track, obj_fh])

# Constraints]
include_constraints(["loop_delta", "free_time"])
con_loop = loop_delta_constraints(model, (1:model.n), 1, T)
con_free_time = free_time_constraints(T)
con = multiple_constraints([con_loop, con_free_time])#, con_pinned_foot])

# Problem
prob = trajectory_optimization_problem(model,
           obj_multi,
           T,
           h = h,
           xl = xl,
           xu = xu,
           ul = ul,
           uu = uu,
		   con = con)

# Trajectory initialization
x0 = linear_interpolation(x1, xT, T) # linear interpolation on state
u0 = [ones(model.m) for t = 1:T-1]

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

# Solve
if false
	include_snopt()
	@time z̄, info = solve(prob, copy(z0),
		nlp = :SNOPT7,
		time_limit = 60 * 5)
	@save joinpath(@__DIR__, "sol_to.jld2") z̄
else
	println("Loading solution...")
	@load joinpath(@__DIR__, "sol_to.jld2") z̄
end

# Unpack solutions
x̄, ū = unpack(z̄, prob)
@show tf = sum([ū[t][end] for t = 1:T-1])
# visualize!(mvis, model, x̄, Δt = ū[1][end])


text:
using UUID

# v1 -------

u1 = UUID.v1()
u1::UUID.Uuid

@assert UUID.get_version(u1) == 1
@assert string(u1) != nothing
@assert int(u1) != nothing
@assert hex(u1) != nothing

@assert UUID.v1() != UUID.v1()  # just checking :)

# v4 -------

u4 = UUID.v4()
u4::UUID.Uuid

@assert UUID.get_version(u4) == 4
@assert string(u4) != nothing
@assert int(u4) != nothing
@assert uint(u4) != nothing
@assert hex(u4) != nothing

@assert UUID.v4() != UUID.v4()


text:
"""
	A global agents control method
"""
module Airspace_Control_module

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

using Constants_module, Stats_module
using Distributions, CAS_module, Aircraft_module

include("check_status.jl")

export Airspace_Controller, Constant_Number_Controller, Constant_Spawn_Rate_Controller
export CAS_Property_Controller, Uniform_CAS, Random_CAS
export Sensor_Property_Controller, Uniform_Sensor, Random_Sensor
export airspace_initialization, airspace_control

abstract type Airspace_Controller end
mutable struct Constant_Number_Controller <: Airspace_Controller
	ac_number::Int64
end
mutable struct Constant_Spawn_Rate_Controller <: Airspace_Controller
	ac_number::Int64 # initial value
	spawn_rate::Float64 # num new ac per time step
end

abstract type CAS_Property_Controller end
mutable struct Uniform_CAS <: CAS_Property_Controller
	the_cas::CAS
end
struct Random_CAS <: CAS_Property_Controller
	the_cas::Void
	function Random_CAS()
		return new(nothing)
	end
end

abstract type Sensor_Property_Controller end
mutable struct Uniform_Sensor <: Sensor_Property_Controller
	the_sensor::Sensor
end
struct Random_Sensor <: Sensor_Property_Controller end

"""
	Initialize airspace
"""
function airspace_initialization(stats::Stats, init_num_ac::Int64, cas_ctrl::CAS_Property_Controller,
	sensor_ctrl::Sensor_Property_Controller, t::Int64)
	All_AC = Array{Aircraft, 1}(0);
	for i = 1 : init_num_ac
		ac = create_ac(cas_ctrl, sensor_ctrl, t)
		push!(All_AC, ac)
	end
	return All_AC
end

"""
	Control the airspace traffic density using constant ac number
"""
# Use airspace_control() after dynamics updated for All_AC
function airspace_control(stats::Stats, All_AC::Array{Aircraft, 1}, as_ctrl::Constant_Number_Controller,
	cas_ctrl::CAS_Property_Controller, sensor_ctrl::Sensor_Property_Controller, t::Int64)
	while length(All_AC) < as_ctrl.ac_number
		ac = create_ac(cas_ctrl, sensor_ctrl, t)
		push!(All_AC, ac)
	end
	for i = 1:length(All_AC)
		All_AC[i].desti = check_desti(All_AC[i])
		if All_AC[i].desti == true
			stats.total_nominal_route_len += All_AC[i].nominal_route_len
			stats.total_route_len += (All_AC[i].actual_route_len + dest_range)
			deleteat!(All_AC, i)
			ac = create_ac(cas_ctrl, sensor_ctrl, t)
			push!(All_AC, ac)
		end
	end
	return All_AC
end

"""
	Control the airspace traffic density using constant spawn rate
	Analogying Poisson process
"""
function airspace_control(stats::Stats, All_AC::Array{Aircraft, 1}, as_ctrl::Constant_Spawn_Rate_Controller,
	cas_ctrl::CAS_Property_Controller, sensor_ctrl::Sensor_Property_Controller, t::Int64)
	# Correction weight control
	# if correction_type == :closest && as_ctrl.spawn_rate > 5900./3600
	# 	policy.env.correction_weight = correction_weight
	# elseif correction_type == :sector && as_ctrl.spawn_rate > 7900./3600
	# 	policy.env.correction_weight = correction_weight
	# end

	# deactivate unqualified ac
	i = 1
	while i <= length(All_AC)
		All_AC[i].desti = check_desti(All_AC[i])
		if All_AC[i].desti == true
			stats.total_nominal_route_len += All_AC[i].nominal_route_len
			stats.total_route_len += (All_AC[i].actual_route_len + dest_range)
			deleteat!(All_AC, i)
		else
			i += 1
		end
	end
	# add new ac following Poisson process
	num_new_ac = rand(Poisson(as_ctrl.spawn_rate))
	for i = 1 : num_new_ac
		ac = create_ac(cas_ctrl, sensor_ctrl, t)
		# while check_colli(ac, All_AC) == true
		# 	ac = create_ac(cas_ctrl, sensor_ctrl, t)
		# end
		push!(All_AC, ac)
	end
	return All_AC
end

##################################################################################
"""
	Helpers:
"""
function init_ac_helper()
	id = Base.Random.uuid1()
	x_0 = rand() * airspace_dim
	y_0 = rand() * airspace_dim
	heading_0 = rand() * 2 * pi - pi
	turn_rate_0 = 0.
	x_dest_0 = rand() * airspace_dim
	y_dest_0 = rand() * airspace_dim
	colli_0 = false
	desti_0 = false
	advisory_0 = :COC # COC
	nominal_route_len = sqrt((x_dest_0 - x_0)^2 + (y_dest_0 - y_0)^2)
	actual_route_len = 0.
	t_to_nominal_t_0 = 0.
	d_to_nominal_d_0 = 1.
	cumulative_num_ac_encountered = 0
	num_ac_encountered_arr = Int64[]
	dist_to_dest = nominal_route_len
	return (id, x_0, y_0, heading_0, turn_rate_0, x_dest_0, y_dest_0,
		colli_0, desti_0, advisory_0, nominal_route_len, actual_route_len,
		t_to_nominal_t_0, d_to_nominal_d_0, cumulative_num_ac_encountered, num_ac_encountered_arr, dist_to_dest)
end

function init_random_cas()
	num_cas = length(subtypes(CAS))
	return random_cas = CAS_DICT[subtypes(CAS)[rand(1:num_cas)]]
end

function init_random_sensor()
	num_sensor = length(subtypes(Sensor))
	return random_sensor = subtypes(Sensor)[rand(1:num_sensor)](sensing_range)
end

function create_ac(cas_ctrl::Uniform_CAS, sensor_ctrl::Uniform_Sensor, t::Int64)
	(id, x_0, y_0, heading_0, turn_rate_0, x_dest_0, y_dest_0,
		colli_0, desti_0, advisory_0, nominal_route_len, actual_route_len,
		t_to_nominal_t_0, d_to_nominal_d_0, cumulative_num_ac_encountered,
		num_ac_encountered_arr, dist_to_dest) = init_ac_helper()
	v_0 = rand() * (max_speed - min_speed) + min_speed
	dyn_0 = Dynamics(x_0, y_0, dt, v_0, heading_0, turn_rate_0)
	return Aircraft(id, dyn_0, x_0, y_0, x_dest_0, y_dest_0, colli_0, desti_0,
		advisory_0, cas_ctrl.the_cas, sensor_ctrl.the_sensor, t, nominal_route_len, actual_route_len, dist_to_dest,
		t_to_nominal_t_0, d_to_nominal_d_0, cumulative_num_ac_encountered, num_ac_encountered_arr)
end

function create_ac(::Random_CAS, sensor_ctrl::Uniform_Sensor, t::Int64)
	(id, x_0, y_0, heading_0, turn_rate_0, x_dest_0, y_dest_0,
		colli_0, desti_0, advisory_0, nominal_route_len, actual_route_len,
		t_to_nominal_t_0, d_to_nominal_d_0, cumulative_num_ac_encountered,
		num_ac_encountered_arr, dist_to_dest) = init_ac_helper()
	rand_cas = init_random_cas()
	v_0 = rand() * (max_speed - min_speed) + min_speed
	dyn_0 = Dynamics(x_0, y_0, dt, v_0, heading_0, turn_rate_0)
	return Aircraft(id, dyn_0, x_0, y_0, x_dest_0, y_dest_0, colli_0, desti_0,
		advisory_0, rand_cas, sensor_ctrl.the_sensor, t, nominal_route_len, actual_route_len, dist_to_dest,
		t_to_nominal_t_0, d_to_nominal_d_0, cumulative_num_ac_encountered, num_ac_encountered_arr)
end

function create_ac(cas_ctrl::Uniform_CAS, ::Random_Sensor, t::Int64)
	(id, x_0, y_0, heading_0, turn_rate_0, x_dest_0, y_dest_0,
		colli_0, desti_0, advisory_0, nominal_route_len, actual_route_len,
		t_to_nominal_t_0, d_to_nominal_d_0, cumulative_num_ac_encountered,
		num_ac_encountered_arr, dist_to_dest) = init_ac_helper()
	rand_sensor = init_random_sensor()
	v_0 = rand() * (max_speed - min_speed) + min_speed
	dyn_0 = Dynamics(x_0, y_0, dt, v_0, heading_0, turn_rate_0)
	return Aircraft(id, dyn_0, x_0, y_0, x_dest_0, y_dest_0, colli_0, desti_0,
		in_airspace_0, advisory_0, rand_cas, rand_sensor, t, nominal_route_len, actual_route_len, dist_to_dest,
		t_to_nominal_t_0, d_to_nominal_d_0, cumulative_num_ac_encountered, num_ac_encountered_arr)
end

function create_ac(::Random_CAS, ::Random_Sensor, t::Int64)
	(id, x_0, y_0, heading_0, turn_rate_0, x_dest_0, y_dest_0,
		colli_0, desti_0, advisory_0, nominal_route_len, actual_route_len,
		t_to_nominal_t_0, d_to_nominal_d_0, cumulative_num_ac_encountered,
		num_ac_encountered_arr, dist_to_dest) = init_ac_helper()
	rand_cas = init_random_cas()
	rand_sensor = init_random_sensor()
	v_0 = rand() * (max_speed - min_speed) + min_speed
	dyn_0 = Dynamics(x_0, y_0, dt, v_0, heading_0, turn_rate_0)
	return Aircraft(id, dyn_0, x_0, y_0, x_dest_0, y_dest_0, colli_0, desti_0,
		in_airspace_0, advisory_0, rand_cas, rand_sensor, t, nominal_route_len, actual_route_len, dist_to_dest,
		t_to_nominal_t_0, d_to_nominal_d_0, cumulative_num_ac_encountered, num_ac_encountered_arr)
end

end # module

max_stars_repo_name: parachutel/DeepCorrectionCAS