text:
function testAddFactor(client, context, factorLabels, factorTypes, factorVariables, factorData)
    resultIds = Task[]
    for (index, label) in enumerate(factorLabels)
        resultId = addFactor(client,context,Factor(label, factorTypes[index], factorVariables[index], factorData[index]))
        @test resultId != "Error"
        push!(resultIds, resultId)
    end

    waitForCompletion(client, resultIds, expectedStatuses=["Complete"])
    return resultIds
end

function testLsf(client, context, factorLabels, factorTypes)
    @test setdiff(factorLabels, fetch( lsf(client, context) )) == []
end

function testGetFactor(client, context, factorLabels, factorTypes)
    for i in 1:length(factorLabels)
        actualFactor = fetch( getFactor(client,context,factorLabels[i]) )
        @test actualFactor["label"] == factorLabels[i]
        @test actualFactor["fnctype"] == factorTypes[i]
    end
end

function testGetFactors(client, context, factorLabels, factorTypes)
    # Make a quick dictionary of the expected variable Types
    factorIdType = Dict(factorLabels .=> factorTypes)

    factors = fetch( getFactors(client, context, detail=FULL) )
    for f in factors
        @test f["fnctype"] == factorIdType[f["label"]]
    end
end

function runFactorTests(client, context)
    @testset "factor-testset" begin
        @info "Running factor-testset"

        # TODO: Refactor as a dictionary, or just do most of this in addFactor.
        factorLabels = ["x0x1f1","x0f1"]
        factorTypes = ["Pose2Pose2","PriorPose2"]
        factorVariables = [["x0", "x1"], ["x0"]]
        factorData = [Pose2Pose2Data(), PriorPose2Data()]

        @testset "Adding" begin testAddFactor(client, context, factorLabels, factorTypes, factorVariables, factorData) end
        @testset "Listing" begin testLsf(client, context, factorLabels, factorTypes) end
        @testset "Getting" begin testGetFactor(client, context, factorLabels, factorTypes) end
        @testset "Getting Lists" begin testGetFactors(client, context, factorLabels, factorTypes) end
    end
end


text:
include("../src/JOH.jl")
include("../MDMKP/MDMKP.jl")
include("SSIT.jl")

use_cplex = true
if use_cplex
	using CPLEX
	optimizer = CPLEX.Optimizer
	get_det_time = CPLEX.CPXgetdettime
else
	using Gurobi
	optimizer = Gurobi.Optimizer
	get_det_time = nothing
end

tight_problems = MDMKP.load_folder()
loose_problems = MDMKP.loosen.(tight_problems)
problems = vcat(tight_problems, loose_problems)
problems = MDMKP.set_decision.(problems)

ass = map(p->p.id.category, problems)
count(i->i=='A', ass)
count(i->i=='B', ass)
count(i->i=='C', ass)
count(i->i=='D', ass)

const ssit_methods = Dict(
	'A'=> JOH.Matheur.SSIT.make_SSIT_method(
			[.0001, .001, .005, .01, .02],
			[60*5,   60*3,  60*3, 60*3, 60*3],
			"A method", 1),
	'B'=> JOH.Matheur.SSIT.make_SSIT_method(
			[.001, .003, .005, .008, .01, .02],
			[60*3, 60*3, 60*3, 60*3, 60*5, 60*5],
			"B method", 1),
	'C'=> JOH.Matheur.SSIT.make_SSIT_method(
			[.005, .01,   .02,   .05],
			[60*3, 60*10, 60*10, 60*10],
			"C method", 1),
	'D'=> JOH.Matheur.SSIT.make_SSIT_method(
			[.005, .01,   .02,   .05],
			[60*3, 60*10, 60*10, 60*10],
			"D method", 1)
)

const base_case = JOH.Matheur.SSIT.make_SSIT_method(
			[.0001],
			[60*60],
			"Base method", 1)

"""
Configured SSIT run using the optimizer configured at the beginning of the file.
"""
function log_ssit_run(problem, method, dir;
		opt=optimizer,
		det_log=get_det_time,
		model_gen=MDMKP.create_MIPS_model)
	SE.log_ssit_run(model_gen(problem, opt), method, dir, opt, det_log)
end

"""
test the problems with the passed case on the corresponding method defined in the const ssit_methods
parameters:
	problems: list of MDMKP problems to test
	case: decision tree classification to test on
		accepts: A, B, C, or D
	loosen: if -1, do not loosen problem, if value from 0 to 1, loosen demand constraints of tested problems via
			demand * loosen
		default: -1
	use_base: should the base case be tested?
		default: true
	res_dir: directory to store results in
		default: results/decision_tree
	ssit_methods: dict of SSIT methods. Must contain entry for the passed case.
"""
function main(problems, case; res_dir="results/decision_tree/", use_base=true, loosen=-1,
		ssit_methods=ssit_methods, base_method=base_case)

	mkpath(res_dir)

	for problem in problems
		if problem.id.category != case #check case matches
			continue
		end

		if 0 < loosen < 1 #check if we should loosen the problem demand constraints
			problem = MDMKP.loosen(problem, percent=loosen, id_increment=0, label=case)
		end

		# make directory for resutls
		problem_dir = mkpath(joinpath(res_dir, "ssit", "$(problem.id.id)"))

		# run SSIT trial
		log_ssit_run(problem, ssit_methods[case], problem_dir)

		# run base test
		if use_base
			base_dir = mkpath(joinpath(res_dir, "base", "$(problem.id.id)"))
			log_ssit_run(problem, base_method, base_dir)
		end
	end
end

println("code loaded.")

#compilation run
main(problems[1:1], 'A', res_dir="results/decision_tree/A_v2/", use_base=false)
println("code compiled.")

main(problems, 'D', res_dir="results/decision_tree/D_95/", loosen=.95)

println("finished. Goodbye.")

ext: jl

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
################################################################################
#
#   (8,3)a
#
################################################################################

# Implementation
# - implementation 1
# - labels <: Any
function getUnitcell_8_3_a(
			unitcell_type  :: Type{U},
			implementation :: Val{1}
		) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

	# return a new Unitcell
	return newUnitcell(
		# Type of the unitcell
		U,
		# lattice vectors
		Vector{Float64}[
			Float64[1.0, 0.0, 0.0],
			Float64[-0.5, sqrt(3)/2., 0.0],
			Float64[0.0, 0.0, (3*sqrt(2))/5.],
		],
		# sites
		S[
			newSite(S, Float64[0.5, sqrt(3)/10., 0.0], getDefaultLabelN(LS,1)),
			newSite(S, Float64[3/5., sqrt(3)/5., (2*sqrt(2))/5.], getDefaultLabelN(LS,2)),
			newSite(S, Float64[0.1, (3*sqrt(3))/10., sqrt(2)/5.], getDefaultLabelN(LS,3)),
			newSite(S, Float64[0.4, sqrt(3)/5., sqrt(2)/5.], getDefaultLabelN(LS,4)),
			newSite(S, Float64[0.0, (2*sqrt(3))/5., 0.0], getDefaultLabelN(LS,5)),
			newSite(S, Float64[-0.1, (3*sqrt(3))/10., (2*sqrt(2))/5.], getDefaultLabelN(LS,6)),
		],
		# bonds
		B[
			newBond(B, 1, 4, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 4, 2, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 4, 3, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 5, 3, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 3, 6, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 4, 1, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 2, 4, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 3, 4, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 3, 5, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 6, 3, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 5, 1, getDefaultLabel(LB), (0, 1, 0)),
			newBond(B, 1, 5, getDefaultLabel(LB), (0, -1, 0)),
			newBond(B, 2, 6, getDefaultLabel(LB), (1, 0, 0)),
			newBond(B, 6, 2, getDefaultLabel(LB), (-1, 0, 0)),
			newBond(B, 1, 2, getDefaultLabel(LB), (0, 0, -1)),
			newBond(B, 2, 1, getDefaultLabel(LB), (0, 0, 1)),
			newBond(B, 5, 6, getDefaultLabel(LB), (0, 0, -1)),
			newBond(B, 6, 5, getDefaultLabel(LB), (0, 0, 1)),
		]
	)
end

# Implementation
# - implementation 4 (Kitaev)
# - labels <: Any
function getUnitcell_8_3_a(
			unitcell_type  :: Type{U},
			implementation :: Val{4}
		) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

	# return a new Unitcell
	return newUnitcell(
		# Type of the unitcell
		U,
		# lattice vectors
		Vector{Float64}[
			Float64[1.0, 0.0, 0.0],
			Float64[-0.5, sqrt(3)/2., 0.0],
			Float64[0.0, 0.0, (3*sqrt(2))/5.],
		],
		# sites
		S[
			newSite(S, Float64[0.5, sqrt(3)/10., 0.0], getDefaultLabelN(LS,1)),
			newSite(S, Float64[3/5., sqrt(3)/5., (2*sqrt(2))/5.], getDefaultLabelN(LS,2)),
			newSite(S, Float64[0.1, (3*sqrt(3))/10., sqrt(2)/5.], getDefaultLabelN(LS,3)),
			newSite(S, Float64[0.4, sqrt(3)/5., sqrt(2)/5.], getDefaultLabelN(LS,4)),
			newSite(S, Float64[0.0, (2*sqrt(3))/5., 0.0], getDefaultLabelN(LS,5)),
			newSite(S, Float64[-0.1, (3*sqrt(3))/10., (2*sqrt(2))/5.], getDefaultLabelN(LS,6)),
		],
		# bonds
		B[
			newBond(B, 1, 4, getDefaultLabelY(LB), (0, 0, 0)),
			newBond(B, 4, 2, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 4, 3, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 5, 3, getDefaultLabelY(LB), (0, 0, 0)),
			newBond(B, 3, 6, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 4, 1, getDefaultLabelY(LB), (0, 0, 0)),
			newBond(B, 2, 4, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 3, 4, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 3, 5, getDefaultLabelY(LB), (0, 0, 0)),
			newBond(B, 6, 3, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 5, 1, getDefaultLabelZ(LB), (0, 1, 0)),
			newBond(B, 1, 5, getDefaultLabelZ(LB), (0, -1, 0)),
			newBond(B, 2, 6, getDefaultLabelY(LB), (1, 0, 0)),
			newBond(B, 6, 2, getDefaultLabelY(LB), (-1, 0, 0)),
			newBond(B, 1, 2, getDefaultLabelX(LB), (0, 0, -1)),
			newBond(B, 2, 1, getDefaultLabelX(LB), (0, 0, 1)),
			newBond(B, 5, 6, getDefaultLabelX(LB), (0, 0, -1)),
			newBond(B, 6, 5, getDefaultLabelX(LB), (0, 0, 1)),
		]
	)
end


text:
#=
1. Create a regular $N^3$ grid T of random tetrahedra:
	(a) Generate four vertices for each tetrahedron randomly in a half-open fixed-size cube.
	(b) Let these cubes form a regular $N^3$ grid.
2. Create a regular $(N − 1)^3$ grid C of such cubes.
3. Align T and C such that the grid nodes of C are at the centers of the grid cells of T .
4. Measure time for T ∪ C.
=#

using ViewerGL; GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation;

function t1t2(i,j,k)
	v1 = rand(3,4).+ [2i,2j,2k]
	v2 = rand(3,4).+ [2i,2j,2k]
	ev = [[1,2],[1,3],[1,4],[2,3],[2,4],[3,4]]
	fv = [[2,3,4],[1,3,4],[1,2,4],[1,2,3]]
	t1 = [v1,fv,ev]
	t2 = [v2,fv,ev]
	return t1,t2
end

N = 7  # N ≥ 2
tets = []
for i=0:N-1, j=0:N-1, k=0.5:0.5
	push!(tets, t1t2(i,j,k)...)
end
V,FV,EV = Lar.struct2lar(Lar.Struct(tets))
GL.VIEW([ GL.GLGrid(V,FV, GL.COLORS[1],1), GL.GLFrame2 ]);

W,(_,EW,FW,_) = Lar.cuboidGrid([1,1,1],true);
U = (W .* [2(N),2(N),2]) .+ [-.5,-.5,-.5]
GL.VIEW([ GL.GLGrid(V,FV,GL.COLORS[1],0.5), GL.GLGrid(U,FW,GL.COLORS[2],0.5), GL.GLFrame2 ]);

V,FV,EV = Lar.struct2lar(Lar.Struct([Lar.Struct(tets), Lar.Struct([(U,FW,EW)]) ]))
GL.VIEW([ GL.GLGrid(V,EV, GL.COLORS[1],0.5), GL.GLFrame2 ]);

	cop_EV = convert(Lar.ChainOp, Lar.coboundary_0(EV::Lar.Cells));
	cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
	W = convert(Lar.Points, V');
	V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement( W, cop_EV, cop_FE);
	W = convert(Lar.Points, V');
	V,CVs,FVs,EVs = Lar.pols2tria(W, copEV, copFE, copCF);

GL.VIEW( GL.GLExplode(V,CVs[3:end],1.5,1.5,1.5,99,0.5) );
GL.VIEW( GL.GLExplode(V,CVs[3:end],1,1,1,99,0.5) );
GL.VIEW( GL.GLExplode(V,CVs[2:2],5,5,5,99,0.5) );
GL.VIEW( GL.GLExplode(V,CVs[12:12],5,5,5,99,0.5) );


text:
using Distributed
@time include("required.jl")

## Conditional arguments
# outcome = "raw" # desired outcome (required)
# outcome = "sdm" # desired outcome (required)
# save_figures = true # should figures be overwritten (optional)

# Make sure "outcome" is defined
if !(@isdefined outcome)
    @warn "'outcome' not defined, must be either 'raw' or 'sdm'"
elseif (outcome != "raw" && outcome != "sdm")
    @warn "'outcome' invalid, must be either 'raw' or 'sdm'"
else
    @info "'outcome' currently set to '$(outcome)'"
end

## Load distributions for all species
@load joinpath("data", "jld2", "$(outcome)-distributions.jld2") distributions spenames speindex
## Load matrix Y
Y = Ymatrix(distributions)

## Weighted endemism
# AOO (Area of occurrence): Species add 1/(number of sites occupied) to site scores

# Number of occupied sites per species
n_occ = map(x -> sum(filter(!isnothing, x)), eachcol(Y))
# Weight site score by number of occupied sites per species
Yweight = hcat([replace(Y[:, i], 1 => 1 / n_occ[i]) for i in 1:size(Y, 2)]...)
# Sum scores per site for all species
endemism_scores = fill(nothing, size(Yweight, 1)) |> Array{Union{Nothing,Float32}}
inds_obs = _indsobs(Yweight)
endemism_scores[inds_obs] = sum.(eachrow(Yweight[inds_obs]))
# Check stats
describe(filter(!isnothing, endemism_scores))
sort(filter(!isnothing, endemism_scores); rev=true) # possibly some outliers, scaling problems for visualization
# Reshape to grid format
endemism_grid = reshape(endemism_scores, size(distributions[1]))

# Create SimpleSDMLayer with endemism values
endemism = SimpleSDMResponse(Array(endemism_grid), distributions[1])

## Plot results
# Raw endemism scores
endemism_plot = plot_layer(
    endemism;
    c=:viridis,
    title="Endemism ($outcome distributions)",
    colorbar_title="Weigthed endemism (area of occurrence)",
)
# Quantile endemism scores
endemism_qplot = plot_layer(
    quantiles(endemism);
    c=:viridis,
    title="Endemism quantiles ($outcome distributions)",
    colorbar_title="Weighted endemism quantile (area of occurrence)",
)

## Save result
# save_figures = true # should figures be overwritten (optional)
if (@isdefined save_figures) && save_figures == true
    @info "Figures saved ($(outcome) lcbd)"
    savefig(endemism_plot, joinpath("fig", outcome, "07_$(outcome)_endemism.png"))
else
    @info "Figures not saved ($(outcome) lcbd)"
end
if (@isdefined save_figures) && save_figures == true
    @info "Figures saved ($(outcome) lcbd)"
    savefig(
        endemism_qplot, joinpath("fig", "quantiles", "07_$(outcome)_endemism_quantiles.png")
    )
else
    @info "Figures not saved ($(outcome) lcbd)"
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


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

using Markdown
using InteractiveUtils

# ╔═╡ 60c13530-6625-11eb-0b68-43bb98153c01
md"""
# Types

Types are a fundamental feature of computer programming that are found in some form in every programming language. They are a way for a programming language to describe the meaning of data stored in memory, describe the rules of how that data can be manipulated, enforce those rules and give the programmer tools to create their own systems of rules for the manipulation of data.

The terminology varies quite a lot between the worlds of mathematics and computer science, and even between different programming languages. The description below takes a largely Julia oriented explanation of types, as the terminology used in discussing Julia's type system is very close to the generic CS definition of the concepts, so makes a good foundation for knowledge of types. Despite the variation in terminology and in the features offered by different languages' type systems, a lot of the concepts transfer to any language.

## Concrete Types

A concrete type can be thought of as a way of interpreting a sequence of bits. For example, the bitstring `10111000` is `184` when viewed as an unsigned integer (`UInt8` in Julia), or `-74` when viewed as a signed integer (`Int8`). The pattern of ones and zeros in memory does not change, only the rules for deciding the meaning of the data changes.
"""

# ╔═╡ 4e34af60-6629-11eb-3df9-537ac7a37f07
bitstring(UInt8(184))

# ╔═╡ 40ef14a0-6627-11eb-27ca-8984ae94c309
bitstring(Int8(-72))

# ╔═╡ 878872c0-662d-11eb-1ee3-dfdc633d89d6
parse(Int, bitstring(Int8(-72)), base=2)

# ╔═╡ f5ad1f90-662c-11eb-1fb3-c96ee1d6be18
md"""
Data can change from being one type to another in several different ways. The most common way is by explicit type conversion, i.e. creating a new instantiation of a type by passing its constructor.

The Char `'A'` can be converted to an Int by passing it to the Int constructor.
"""

# ╔═╡ fe89a5c0-662c-11eb-21a9-1590e2af58cc
Int8('A')

# ╔═╡ 6868aae0-6632-11eb-22b3-ad8f1aa5654a
md"""
Likewise, Ints can be converted back to Chars with the Char constructor.
"""

# ╔═╡ cfd58da0-662e-11eb-0bb3-832ee64a0194
Char(65)

# ╔═╡ 0ab3be60-662f-11eb-3b12-ad38972ef557
md"""
In strongly typed languages, this is the only way to convert data from one type to another. In weakly typed languages, such as C, types can be changed by *casting* them from one type to another. This is where the bits in memory are treated as a different type by disregarding its previous type and treating that bitstring as if it is its new type.

```c
int i = (int)true; // i is now 1
i = (int)false; // i is now 0
```

This can backfire in situations where there are unexpected differences in the underlying hardware. For example, different systems can use different [endianness](https://en.wikipedia.org/wiki/Endianness). In general, the order of bits in a byte goes from most significant bit (MSB) to least significant bit (LSB). However, in types larger than a single byte, the order of bytes can differ from system to system and even between types.
"""

# ╔═╡ 6603f930-662d-11eb-304d-771cf3e9011f
bitstring(Char(65))

# ╔═╡ 68e40eb0-662d-11eb-2152-33ab3770c8b1
bitstring(Int32(65))

# ╔═╡ ec4d6a30-662d-11eb-11e8-711d50e4f58d
begin
	eachbyte(val) = join.([Iterators.partition(bitstring(val), 8)...])
	md"""
	The bitstring of the `Int32` is represented most significant byte first, or big-endian, and the `Char` is represented least significant byte first, or little-endian.
	"""
end

# ╔═╡ e38ae590-6636-11eb-35e0-8dd86e537f4f
Char(Int32(65)) == 'A'

# ╔═╡ 0b9de562-6631-11eb-077a-eb21ad4de957
eachbyte(Int32(65))

# ╔═╡ 6d3603b0-662d-11eb-02a1-d5dffdd90d1e
eachbyte('A')

# ╔═╡ dc468b90-6636-11eb-2c2d-efac00304e64
md"""
Even though they both represent the number 65, the underlying bit pattern in memory is different. A C style casting of the `Char` `'A'` to an `Int32` results in a number that is very much not `65`.
"""

# ╔═╡ 709b27b0-662d-11eb-2618-fdc4bfff6279
parse(Int32, bitstring('A'), base=2)

# ╔═╡ f90c2c40-6630-11eb-06d0-b3877f1e7653
md"""
Other languages, such as JavaScript, allow a loose comparison of values across types. There are two different types of comparison, `==` and `===` (which don't have names, but I call twoquals and threequals, so feel free to use those), which differ on whether the type information is used when determining equality.

```js
> "0" == 0
true

> "0" === 0
false

> false == 0
true

> false === 0
false
```

This is due to another form of type conversion, called implicit type conversion or duck typing. This is when the type of an object is helpfully changed to whatever the interpreter thinks your code probably means based on the types passed to it.

In JavaScript, this is most famously known from the confusion caused by implicit conversion between Numbers and Strings. In particular, the `+` operator means addition for Numbers and concatenation for Strings. However, due to the loose conversions between types, if a string can also represent a number, and the operator can only be used on numbers, such as `-`, `*` or `/`, both strings will be implicitly converted to numbers.

```js
> 1 + 1
2

> "1" + 1
"11"

> "1" - "1"
0

> "1" + "1" - "1"
10
```

This is either perfectly sensible or jaw-clenchingly infuriating, depending on how you feel about this sort of thing. I don't mind JavaScript's semantics, but I can see why other people hate it. Let's compare with Julia.
"""

# ╔═╡ ed740ee0-6638-11eb-2990-a38eccca8f31
md"""
Julia lets you increment the value of a `Char`, but adding an `Int` to a `String` is not allowed. Or rather, not implemented by default. Due to the magic of operator overloading and sadomasochism, it's possible to add JavaScript's `+` semantics to Julia.
"""

# ╔═╡ 82b00a80-663f-11eb-3f2d-891a54c4e511
begin
	import Base: +, -
	+(a::String, b::String) = a * b
	-(a::String, b::String) = parse(Int, a) - parse(Int, b)
	"1" + "1" - "1"
end

# ╔═╡ d9e5f460-6638-11eb-3aef-e517d47f7251
'1' + 1

# ╔═╡ db84fe60-6638-11eb-2828-bbb42b4b733a
"1" + 1

# ╔═╡ 8c46a0a0-6639-11eb-1af7-b583bfcbdc20
let	+(a::String, b::Int64) = a * string(b)
	"1" + 1
end

# ╔═╡ 901bd8a0-6646-11eb-1011-49da944c1584
md"""
Mmmmmm, *web scale* and *full stack*. 👌😎🔥🔥🔥🔥🔥🔥

While ridiculous and nauseating to many, this example helps demonstrate why anyone cares about types in the first place. Annotating the types of function parameters on the function definition tells the compiler which types that function definition (or *method*) should handle. In the example above, we defined a new method on the `+` function, so that `+` can now handle being passed `String` arguments.

Calling a different method on a function based on the types of the arguments is called *multiple dispatch*, and is very handy. It lets us add our super neato JS style string addition and subtraction without overwriting the boring normal way of adding and subtracting of numbers. The `+` functions can live happily side by side.
"""

# ╔═╡ 4231c040-6647-11eb-35d9-2901001b9f20
1 + 1        # your grampa's addition

# ╔═╡ 2bb893c0-6647-11eb-174c-d948da03a90a
"1" + "1"    # cool refreshing NEW Addition™

# ╔═╡ ccbb2d40-6639-11eb-3ac3-d1e78b34f12e
md"""
## Abstract Types

As mentioned before, the types we've been referring to have all been *concrete* types, that is, types which have a particular representation in memory. Due the practical limitation of only having a finite amount of memory available on any given computer, it's necessary to only use a certain number of bits to represent any given value. This means that any concrete type has a limited range of values it can represent.

`Int8` can only represent values in the range `-128:127`, because it only has eight bits available, so the maximum number of values it can represent is `2^8`. `Int64` can represent a considerably larger range of values.
"""

# ╔═╡ 662f9f40-663c-11eb-27dd-5f3a6253d396
typemin(Int8):typemax(Int8)

# ╔═╡ 608a9d5e-663c-11eb-39e0-f70eabad7b59
typemin(Int64):typemax(Int64)

# ╔═╡ fd3a0100-663c-11eb-2b00-2741aa087dc9
md"""
Even though these are technically different types, they represent the same concept, so it's useful to be able to refer to operations on the common value that these types represent, like with the `+` operator in the example above.

Now what if some people don't like our awesome new `+`? Well, we could make our own `addclassic` function that doesn't have the string methods, just like the `+` used to. This could be done using our old pal multiple dispatch again, something like:

```julia
addclassic(a::Int8, b::Int8) = a + b
addclassic(a::Int16, b::Int16) = a + b
addclassic(a::Int32, b::Int32) = a + b
addclassic(a::Int64, b::Int64) = a + b
```

This would work fine, but it would get tedious having to do this for every function that we wanted to work across multiple types, and this doesn't even cover what happens if we wanted to add an `Int64` to an `Int32`, or an `Int8` to an `Int16`, or any of the other combinations. Luckily there's an easier way.

> “If it walks like a duck, quacks like a duck, and subtypes `AbstractDuck`, it's a duck.”
> 	— *Dijkstra, probably.*

Abstract types are a way of arranging types in a hierarchy, so that methods can be written for a type further up the hierarchy (called a *supertype*), and that function can be called by any of the types that are below it in the hierarchy (*subtypes*).
"""

# ╔═╡ 2e3b8990-6651-11eb-2ffa-b1600faf3e74
supertype.((Int8, Int16, Int32, Int64))

# ╔═╡ 2fbfea2e-6652-11eb-0aa2-0b07461c90b0
md"""
When each of those types was defined, it was assigned the supertype `Signed`, which is the supertype which represents all signed integer types.

```julia
type Int8 <: Signed end
```

Since each of the types we wrote `add` methods for share a supertype, so we could instead write the following code and have the same functionality.

```julia
add(a::Signed, b::Signed) = a + b
```

`Signed` itself also has further supertypes, so let's have a look at those.
"""

# ╔═╡ 9e84d930-6652-11eb-2a7a-d543e31df803
supertype(Signed)

# ╔═╡ 0d1571c0-6653-11eb-1b4c-ebe1e40a2c09
supertype(Integer)

# ╔═╡ 0fe581ae-6653-11eb-0fae-dba6c6847d8b
supertype(Real)

# ╔═╡ 132e9230-6653-11eb-2ae6-ff6393016405
supertype(Number)

# ╔═╡ 15c0aec0-6653-11eb-139a-c7c7bc9d4496
supertype(Any)

# ╔═╡ 19443030-6653-11eb-2ab4-593b361d76f3
md"""
As you can see, the type hierarchy goes quite high, all the way up to `Any`, the supertype of all supertypes that's even its own supertype. All the types form a big tree with `Any` sitting at the top (in the mathematical sense of a tree, it's not my fault mathematicians don't know which way up trees go). `Any` is the default type, and variables and parameters which don't have a type annotation default to `Any`.

Our `add` function could target `Any`, with either of the following equivalent definitions.

```julia
addclassic(a::Any, b::Any) = a + b
addclassic(a, b) = a + b
```

This leaves us back with our JS style string concatenation that we were defining our special strict `add` function to get away from. We want to add numbers and nothing else, so let's target `Number`.
"""

# ╔═╡ aba57c40-6653-11eb-0fa6-01eae675683e
addclassic(a::Number, b::Number) = a + b

# ╔═╡ 83505ca2-6654-11eb-3a24-8b83a78ec479
addclassic(1, 1)

# ╔═╡ 860f7ca0-6654-11eb-2536-27ebe0cda50d
addclassic(Int8(1), Int8(1))

# ╔═╡ 91285e40-6654-11eb-3960-b3226865c76d
addclassic("1", "1")

# ╔═╡ 95b11ab0-6654-11eb-1dfc-7d7103cb1a5e
md"""
Success! `addclassic` works just like `+` used to before we upgraded it. Now nobody can complain.
"""

# ╔═╡ Cell order:
# ╟─60c13530-6625-11eb-0b68-43bb98153c01
# ╠═4e34af60-6629-11eb-3df9-537ac7a37f07
# ╠═40ef14a0-6627-11eb-27ca-8984ae94c309
# ╠═878872c0-662d-11eb-1ee3-dfdc633d89d6
# ╟─f5ad1f90-662c-11eb-1fb3-c96ee1d6be18
# ╠═fe89a5c0-662c-11eb-21a9-1590e2af58cc
# ╟─6868aae0-6632-11eb-22b3-ad8f1aa5654a
# ╠═cfd58da0-662e-11eb-0bb3-832ee64a0194
# ╟─0ab3be60-662f-11eb-3b12-ad38972ef557
# ╠═6603f930-662d-11eb-304d-771cf3e9011f
# ╠═68e40eb0-662d-11eb-2152-33ab3770c8b1
# ╟─ec4d6a30-662d-11eb-11e8-711d50e4f58d
# ╠═e38ae590-6636-11eb-35e0-8dd86e537f4f
# ╠═0b9de562-6631-11eb-077a-eb21ad4de957
# ╠═6d3603b0-662d-11eb-02a1-d5dffdd90d1e
# ╟─dc468b90-6636-11eb-2c2d-efac00304e64
# ╠═709b27b0-662d-11eb-2618-fdc4bfff6279
# ╟─f90c2c40-6630-11eb-06d0-b3877f1e7653
# ╠═d9e5f460-6638-11eb-3aef-e517d47f7251
# ╠═db84fe60-6638-11eb-2828-bbb42b4b733a
# ╟─ed740ee0-6638-11eb-2990-a38eccca8f31
# ╠═8c46a0a0-6639-11eb-1af7-b583bfcbdc20
# ╠═82b00a80-663f-11eb-3f2d-891a54c4e511
# ╟─901bd8a0-6646-11eb-1011-49da944c1584
# ╠═4231c040-6647-11eb-35d9-2901001b9f20
# ╠═2bb893c0-6647-11eb-174c-d948da03a90a
# ╟─ccbb2d40-6639-11eb-3ac3-d1e78b34f12e
# ╠═662f9f40-663c-11eb-27dd-5f3a6253d396
# ╠═608a9d5e-663c-11eb-39e0-f70eabad7b59
# ╟─fd3a0100-663c-11eb-2b00-2741aa087dc9
# ╠═2e3b8990-6651-11eb-2ffa-b1600faf3e74
# ╟─2fbfea2e-6652-11eb-0aa2-0b07461c90b0
# ╠═9e84d930-6652-11eb-2a7a-d543e31df803
# ╠═0d1571c0-6653-11eb-1b4c-ebe1e40a2c09
# ╠═0fe581ae-6653-11eb-0fae-dba6c6847d8b
# ╠═132e9230-6653-11eb-2ae6-ff6393016405
# ╠═15c0aec0-6653-11eb-139a-c7c7bc9d4496
# ╟─19443030-6653-11eb-2ab4-593b361d76f3
# ╠═aba57c40-6653-11eb-0fa6-01eae675683e
# ╠═83505ca2-6654-11eb-3a24-8b83a78ec479
# ╠═860f7ca0-6654-11eb-2536-27ebe0cda50d
# ╠═91285e40-6654-11eb-3960-b3226865c76d
# ╟─95b11ab0-6654-11eb-1dfc-7d7103cb1a5e


text:
module URIs

import Base.==

include("urlparser.jl")

export URI, URL,
       hasscheme, scheme,
       hashostname, hostname,
       haspath, path,
       hasquery, query,
       hasfragment, fragment,
       hasuserinfo, userinfo,
       hasport, port,
       resource, host,
       escape, unescape,
       splitpath

"""
    HTTP.URL(host; userinfo="", path="", query="", fragment="", isconnect=false)
    HTTP.URI(; scheme="", hostname="", port="", ...)
    HTTP.URI(str; isconnect=false)
    parse(HTTP.URI, str::String; isconnect=false)

A type representing a valid uri. Can be constructed from distinct parts using the various
supported keyword arguments. With a raw, already-encoded uri string, use `parse(HTTP.URI, str)`
to parse the `HTTP.URI` directly. The `HTTP.URI` constructors will automatically escape any provided
`query` arguments, typically provided as `"key"=>"value"::Pair` or `Dict("key"=>"value")`.
Note that multiple values for a single query key can provided like `Dict("key"=>["value1", "value2"])`.

For efficiency, the internal representation is stored as a set of offsets and lengths to the various uri components.
To access and return these components as strings, use the various accessor methods:
  * `HTTP.scheme`: returns the scheme (if any) associated with the uri
  * `HTTP.userinfo`: returns the userinfo (if any) associated with the uri
  * `HTTP.hostname`: returns the hostname only of the uri
  * `HTTP.port`: returns the port of the uri; will return "80" or "443" by default if the scheme is "http" or "https", respectively
  * `HTTP.host`: returns the "hostname:port" combination
  * `HTTP.path`: returns the path for a uri
  * `HTTP.query`: returns the query for a uri
  * `HTTP.fragment`: returns the fragment for a uri
  * `HTTP.resource`: returns the path-query-fragment combination
"""
struct URI
    data::Vector{UInt8}
    offsets::NTuple{7, Offset}
end

function URI(;hostname::String="", path::String="",
            scheme::String="", userinfo::String="",
            port::Union{Integer,String}="", query="",
            fragment::String="", isconnect::Bool=false)
    hostname != "" && scheme == "" && !isconnect && (scheme = "http")
    io = IOBuffer()
    printuri(io, scheme, userinfo, hostname, string(port), path, escape(query), fragment)
    return Base.parse(URI, String(take!(io)); isconnect=isconnect)
end

# we assume `str` is at least hostname & port
# if all others keywords are empty, assume CONNECT
# can include path, userinfo, query, & fragment
function URL(str::String; userinfo::String="", path::String="",
                          query="", fragment::String="",
                          isconnect::Bool=false)
    if str != ""
        if startswith(str, "http") || startswith(str, "https")
            str = string(str, path, ifelse(query == "", "", "?" * escape(query)),
                         ifelse(fragment == "", "", "#$fragment"))
        else
            if startswith(str, "/") || str == "*"
                # relative uri like "/" or "*", leave it alone
            elseif path == "" && userinfo == "" && query == "" && fragment == "" && ':' in str
                isconnect = true
            else
                str = string("http://", userinfo == "" ? "" : "$userinfo@",
                             str, path, ifelse(query == "", "", "?" * escape(query)),
                             ifelse(fragment == "", "", "#$fragment"))
            end
        end
    end
    return Base.parse(URI, str; isconnect=isconnect)
end
URI(str::String; isconnect::Bool=false) = Base.parse(URI, str; isconnect=isconnect)
Base.parse(::Type{URI}, str::String; isconnect::Bool=false) = http_parser_parse_url(Vector{UInt8}(str), 1, sizeof(str), isconnect)

==(a::URI,b::URI) = scheme(a)   == scheme(b)    &&
                    hostname(a) == hostname(b)  &&
                    path(a)     == path(b)      &&
                    query(a)    == query(b)     &&
                    fragment(a) == fragment(b)  &&
                    userinfo(a) == userinfo(b)  &&
                    ((!hasport(a) || !hasport(b)) || (port(a) == port(b)))

# accessors
for uf in instances(http_parser_url_fields)
    uf == UF_MAX && break
    nm = lowercase(string(uf)[4:end])
    has = Symbol(string("has", nm))
    @eval $has(uri::URI) = uri.offsets[Int($uf)].len > 0
    uf == UF_PORT && continue
    @eval $(Symbol(nm))(uri::URI) = String(uri.data[uri.offsets[Int($uf)]])
end

# special def for port
function port(uri::URI)
    if hasport(uri)
        return String(uri.data[uri.offsets[Int(UF_PORT)]])
    else
        sch = scheme(uri)
        return sch == "http" ? "80" : sch == "https" ? "443" : ""
    end
end

resource(uri::URI; isconnect::Bool=false) = isconnect ? host(uri) : path(uri) * (isempty(query(uri)) ? "" : "?$(query(uri))") * (isempty(fragment(uri)) ? "" : "#$(fragment(uri))")
host(uri::URI) = hostname(uri) * (isempty(port(uri)) ? "" : ":$(port(uri))")

Base.show(io::IO, uri::URI) = print(io, "HTTP.URI(\"", uri, "\")")

Base.print(io::IO, u::URI) = printuri(io, scheme(u), userinfo(u), hostname(u), port(u), path(u), query(u), fragment(u))
function printuri(io::IO, sch::String, userinfo::String, hostname::String, port::String, path::String, query::String, fragment::String)
    if sch in uses_authority
        print(io, sch, "://")
        !isempty(userinfo) && print(io, userinfo, "@")
        print(io, ':' in hostname ? "[$hostname]" : hostname)
        print(io, ((sch == "http" && port == "80") ||
                   (sch == "https" && port == "443") || isempty(port)) ? "" : ":$port")
    elseif path != "" && path != "*" && sch != ""
        print(io, sch, ":")
    elseif hostname != "" && port != "" # CONNECT
        print(io, hostname, ":", port)
    end
    if (isempty(hostname) || hostname[end] != '/') &&
       (isempty(path) || path[1] != '/') &&
       (!isempty(fragment) || !isempty(path))
        path = (!isempty(sch) && sch == "http" || sch == "https") ? string("/", path) : path
    end
    print(io, path, isempty(query) ? "" : "?$query", isempty(fragment) ? "" : "#$fragment")
end

# Validate known URI formats
const uses_authority = ["hdfs", "ftp", "http", "gopher", "nntp", "telnet", "imap", "wais", "file", "mms", "https", "shttp", "snews", "prospero", "rtsp", "rtspu", "rsync", "svn", "svn+ssh", "sftp" ,"nfs", "git", "git+ssh", "ldap", "s3"]
const uses_params = ["ftp", "hdl", "prospero", "http", "imap", "https", "shttp", "rtsp", "rtspu", "sip", "sips", "mms", "sftp", "tel"]
const non_hierarchical = ["gopher", "hdl", "mailto", "news", "telnet", "wais", "imap", "snews", "sip", "sips"]
const uses_query = ["http", "wais", "imap", "https", "shttp", "mms", "gopher", "rtsp", "rtspu", "sip", "sips", "ldap"]
const uses_fragment = ["hdfs", "ftp", "hdl", "http", "gopher", "news", "nntp", "wais", "https", "shttp", "snews", "file", "prospero"]

"checks if a `HTTP.URI` is valid"
function Base.isvalid(uri::URI)
    sch = scheme(uri)
    isempty(sch) && throw(ArgumentError("can not validate relative URI"))
    if ((sch in non_hierarchical) && (search(path(uri), '/') > 1)) ||       # path hierarchy not allowed
       (!(sch in uses_query) && !isempty(query(uri))) ||                    # query component not allowed
       (!(sch in uses_fragment) && !isempty(fragment(uri))) ||              # fragment identifier component not allowed
       (!(sch in uses_authority) && (!isempty(hostname(uri)) || ("" != port(uri)) || !isempty(userinfo(uri)))) # authority component not allowed
        return false
    end
    return true
end

# RFC3986 Unreserved Characters (and '~' Unsafe per RFC1738).
@inline issafe(c::Char) = c == '-' ||
                          c == '.' ||
                          c == '_' ||
                          (isascii(c) && isalnum(c))

utf8_chars(str::AbstractString) = (Char(c) for c in Vector{UInt8}(str))

"percent-encode a string, dict, or pair for a uri"
function escape end

escape(c::Char) = string('%', uppercase(hex(c,2)))
escape(str::AbstractString, safe::Function=issafe) =
    join(safe(c) ? c : escape(c) for c in utf8_chars(str))

escape(bytes::Vector{UInt8}) = bytes
escape(v::Number) = escape(string(v))
escape(v::Nullable) = Base.isnull(v) ? "" : escape(get(v))

escape(key, value) = string(escape(key), "=", escape(value))
escape(key, values::Vector) = escape(key => v for v in values)
escape(query) = join((escape(k, v) for (k,v) in query), "&")

"unescape a percent-encoded uri/url"
function unescape(str)
    out = IOBuffer()
    i = 1
    while !done(str, i)
        c, i = next(str, i)
        if c == '%'
            c1, i = next(str, i)
            c, i = next(str, i)
            write(out, Base.parse(UInt8, string(c1, c), 16))
        else
            write(out, c)
        end
    end
    return String(take!(out))
end

"""
Splits the path into components
See: http://tools.ietf.org/html/rfc3986#section-3.3
"""
function splitpath(uri::URI, starting=2)
    return splitpath(path(uri), starting)
end

function splitpath(p::String, starting=2)
    elems = String[]
    len = length(p)
    len > 1 || return elems
    start_ind = i = ifelse(p[1] == '/', 2, 1)
    while true
        c = p[i]
        if c == '/'
            push!(elems, p[start_ind:i-1])
            start_ind = i + 1
        elseif i == len
            push!(elems, p[start_ind:i])
        end
        i += 1
        (i > len || c in ('?', '#')) && break
    end
    return elems
end

end # module


text:
using Test

using DEC2D
using LinearAlgebra

@testset "Geometry algorithms" begin
    @testset "Mid vertices" begin
        @test DEC2D.mid_point([2,2], [4,4]) == [3,3]
        @test DEC2D.mid_point([-2,-2], [4,4]) == [1,1]
    end

    @testset "Orientation" begin
        @test DEC2D.orientation([-1,0], [0,1], [1,0]) == DEC2D.Clockwise
        @test DEC2D.orientation([1,0], [0,1], [-1,0]) == DEC2D.CounterClockwise
        @test DEC2D.orientation([1,0], [0,0], [-1,0]) == DEC2D.Colinear
        @test DEC2D.orientation([0,0], [0,0], [-1,0]) == DEC2D.Colinear
    end

    @testset "Circumcenter" begin
        @test DEC2D.circumcenter([0,0], [0,3], [4,0]) ≈ [2.0, 1.5]
        @test DEC2D.circumcenter([0,0], [0,-3], [-4,0]) ≈ -[2.0, 1.5]
        @test DEC2D.circumcenter([0,1], [-cos(π/3),-sin(π/3)], [cos(π/3),-sin(π/3)]) ≈ [0.0, 0.0]
    end

    @testset "Length" begin
        @test DEC2D.length([0,0], [0, 1.0]) ≈ 1.0
        @test DEC2D.length([0,0], [1, 1.0]) ≈ √2.0
    end

    @testset "Area" begin
        @test DEC2D.area([[0,0], [0,1], [1, 0]]) ≈ 0.5
        @test DEC2D.area([[0,0], [0,1], [1, 1], [1, 0]]) ≈ 1.0

        areapolyn(n) = n * sin(π/n) * cos(π/n)
        polyn(n) = [[cos(i * 2π / n), sin(i * 2π / n)] for i in 1:n]
        @test DEC2D.area(polyn(3)) ≈ areapolyn(3)
        @test DEC2D.area(polyn(4)) ≈ areapolyn(4)
        @test DEC2D.area(polyn(7)) ≈ areapolyn(7)

        @test DEC2D.area([p .+ 1 for p in polyn(3)]) ≈ areapolyn(3)
        @test DEC2D.area([p .+ 1 for p in polyn(4)]) ≈ areapolyn(4)
        @test DEC2D.area([p .+ 1 for p in polyn(7)]) ≈ areapolyn(7)
    end
end

@testset "Triangulation and Dual" begin
    @testset "Parted triangle" begin
        vertices = vcat([[0,0]], [[cos(i * 2π/3), sin(i * 2π/3)] for i in 0:2])
        edges = [[3, 4], [1, 3], [2, 4], [2, 3], [1, 4], [1, 2]]
        triangles = [[1, 2, 3], [1, 3, 4], [1, 4, 2]]

        mesh = DEC2D.DECMesh(vertices, edges, triangles)

        @show mesh.tri.vertices
        @show mesh.tri.edges
        @show mesh.tri.triangles

        println("")

        @show mesh.dtri.n
        @show mesh.dtri.dtriangles
        @show mesh.dtri.dedges
        @show mesh.dtri.dvertices
    end
end


text:
"""
    PolicyWrapper

Flexible utility wrapper for a policy designed for collecting statistics about planning.

Carries a function, a policy, and optionally a payload (that can be any type).

The function should typically be defined with the do syntax. Each time `action` is called on the wrapper, this function will be called.

If there is no payload, it will be called with two argments: the policy and the state/belief. If there is a payload, it will be called with three arguments: the policy, the payload, and the current state or belief. The function should return an appropriate action. The idea is that, in this function, `action(policy, s)` should be called, statistics from the policy/planner should be collected and saved in the payload, exceptions can be handled, and the action should be returned.

# Example

    using POMDPModels
    using POMDPToolbox

    mdp = GridWorld()
    policy = RandomPolicy(mdp)
    counts = Dict(a=>0 for a in iterator(actions(mdp)))

    # with a payload
    statswrapper = PolicyWrapper(policy, payload=counts) do policy, counts, s
        a = action(policy, s)
        counts[a] += 1
        return a
    end

    h = simulate(HistoryRecorder(max_steps=100), mdp, statswrapper)
    for (a, count) in payload(statswrapper)
        println("policy chose action \$a \$count of \$(n_steps(h)) times.")
    end

    # without a payload
    errwrapper = PolicyWrapper(policy) do policy, s
        try
            a = action(policy, s)
        catch ex
            warn("Caught error in policy; using default")
            a = :left
        end
        return a
    end

    h = simulate(HistoryRecorder(max_steps=100), mdp, errwrapper)
"""
mutable struct PolicyWrapper{P<:Policy, F<:Function, PL} <: Policy
    f::F
    policy::P
    payload::Nullable{PL}
end

function PolicyWrapper(f::Function, policy::Policy; payload=Nullable())
    return PolicyWrapper(f, policy, convert(Nullable, payload))
end

function PolicyWrapper(policy::Policy; payload=Nullable())
    return PolicyWrapper((p,s)->action(p.policy,s), policy, convert(Nullable, payload))
end

function action(p::PolicyWrapper, s)
    if isnull(p.payload)
        return p.f(p.policy, s)
    else
        return p.f(p.policy, get(p.payload), s)
    end
end

updater(p::PolicyWrapper) = updater(p.policy)

payload(p::PolicyWrapper) = get(p.payload, nothing)

Base.srand(p::PolicyWrapper, seed) = srand(p.policy, seed)


text:
export discretizepolygon

"""
Discretize the polygon given by the (complex) entries of Z.

# Arguments
- `z`: Polygon vertices given as complex values. Pass an empty array for the unit disk.
  Pass a single value for a disk centered at that value. Pass two values for Chebyshev points.
- `include_interior_points`: Whether to return interior points.
- `npts`: Number of boundary points to return.
- `nptsint`: Minimum number of interior points (more than this may be returned).

# Return values
- `boundary_points`: Boundary points (`npts` of these), followed by the input
  polygon vertices, with an additional final vertex equal to the first vertex.
- `interior_points`: Interior points, if specified to be returned, otherwise an
  empty vector. If specified, there will be at least `nptsint` of these.
"""
function discretizepolygon(
    z::Vector{CT} = [],
    include_interior_points::Bool = false,
    npts::Int = 10000,
    nptsint::Int = 5) where CT<:Complex{<:Real}

    if isempty(z) # return unit disk
        z = [CT(0)]
    end

    if length(z) == 1 # return disk centered at z
        zz = z[1] .+ cis.(2 * pi * (1:npts)/npts)
    elseif length(z) == 2 # return Chebyshev points
        zz = (z[2]-z[1])/2 * (cos.(pi*(npts-1:-1:0)/(npts-1)).+1) .+ z[1]
    else
        z = [z[:]; z[1]] # close polygon
        L = sum(abs.(diff(z))) # length of polygon

        # 0<=alph<=1 is position of current point between z[ind] and z[ind+1]
        ind = 1
        alph = 0.0
        zz = [z[ind]]
        remL = L/npts
        while length(zz) < npts
            d = abs(z[ind+1] - z[ind])
            if (1-alph)*d < remL # can go to end of edge
                ind += 1
                remL -= (1-alph) * d
                alph = 0.0
            else # next point zz is between z[ind] and z[ind+1]
                alph += remL / d
                remL = L / npts
                push!(zz, z[ind] + alph * (z[ind+1] - z[ind]))
            end
        end
    end

    append!(zz, z)

    Z = Vector{CT}()
    if include_interior_points
        if length(z) == 2 # interval case
            xnr = 2*nptsint
            if (xnr&1) == 0 xnr += 1 end
            xpts = range(z[1], stop = z[2], length = xnr)
            Z = collect(xpts[2:2:end])
            return zz, Z
        end

        points = length(z) == 1 ? zz : z
        realz = real(points)
        imagz = imag(points)
        real_min, real_max = extrema(realz)
        imag_min, imag_max = extrema(imagz)

        iter = 0
        spacing = (real_max - real_min) / 2.0001 / sqrt(nptsint)
        while length(Z) < nptsint
            iter += 1
            if iter > 10
                error("Failed to find interior polygon points. Polygon too narrow? (Note that intervals should be given by their two endpoints only.)")
            end

            xnr = trunc(Int, (real_max - real_min) / (2*spacing))
            ynr = trunc(Int, (imag_max - imag_min) / (2*spacing))
            spacing /= sqrt(sqrt(2))
            if xnr <= 1 || ynr <= 1
                continue
            end

            xpts = range(real_min, stop = real_max, length = xnr)
            xpts = xpts[2:2:end]

            ypts = range(imag_min-eps(), stop = imag_max+eps(), length = ynr)
            ypts = ypts[2:2:end]

            Z = [x + y*im for x in xpts for y in ypts]::Vector{CT}
            Z = filter(p -> inpolygon(real(p), imag(p), realz, imagz), Z)
        end
    end

    return zz, Z
end


text:
isabspath(ARGS[1]) || error("Please use an absolute path to the output folder.")
println("Reading experiments from: '$(ARGS[1])'")
data_output_root = ARGS[1]

using OneClassActiveLearning, SVDD
using DataFrames, CSV
using Unmarshal, JSON
using Logging, Statistics
using StatsBase

#### Simplify naming

ALIASES = Dict{Symbol, Dict{String, Union{String, Missing}}}(
    :init_strategy_C => Dict("BoundedTaxErrorEstimate" => "TaxBound",
                             "TaxErrorEstimate" => "Tax",
                             "BinarySearchCStrategy" => "BinarySearch",
                             "FixedCStrategy" => "Fix"),
    :init_strategy_gamma  => Dict("FixedGammaStrategy" => "Fix",
                                  "WangGammaStrategy" =>  "Wang",
                                  "WangCombinedInitializationStrategy" => "WangTax",
                                  "RuleOfThumbSilverman" => "Silverman",
                                  "RuleOfThumbScott" => "Scott")
)

function simplify(str, key)
   for (pattern, alias) in ALIASES[key]
       occursin(pattern, str) && return alias
    end
    return str
end

##### Build CSV

array_stats(x) = minimum(x), mean(x), maximum(x)
prefix_named_tuples(x, prefix) = merge(NamedTuple(), [Symbol("$(prefix)_$k")=>v for (k, v) in pairs(x)])

function get_time_mem_stats(r::OneClassActiveLearning.Result)
    if !haskey(r.al_history, :time_fit) || !haskey(r.al_history, :mem_fit) ||
       !haskey(r.al_history, :time_qs) || !haskey(r.al_history, :mem_qs)
        return (time_fit_min=missing, time_fit_mean=missing, time_fit_max=missing,
                mem_fit_min=missing, mem_fit_mean=missing, mem_fit_max=missing,
                time_qs_min=missing, time_qs_mean=missing, time_qs_max=missing,
                mem_qs_min=missing, mem_qs_mean=missing, mem_qs_max=missing)
    end
    total_run_time = r.al_summary[:runtime][:time_exp]
    time_fit_min, time_fit_mean, time_fit_max = array_stats(values(r.al_history, :time_fit))
    mem_fit_min, mem_fit_mean, mem_fit_max = array_stats(values(r.al_history, :mem_fit))
    time_qs_min, time_qs_mean, time_qs_max = array_stats(values(r.al_history, :time_qs))
    mem_qs_min, mem_qs_mean, mem_qs_max = array_stats(values(r.al_history, :mem_qs))
    return (total_run_time=total_run_time,
            time_fit_min=time_fit_min, time_fit_mean=time_fit_mean, time_fit_max=time_fit_max,
            mem_fit_min=mem_fit_min, mem_fit_mean=mem_fit_mean, mem_fit_max=mem_fit_max,
            time_qs_min=time_qs_min, time_qs_mean=time_qs_mean, time_qs_max=time_qs_max,
            mem_qs_min=mem_qs_min, mem_qs_mean=mem_qs_mean, mem_qs_max=mem_qs_max)
end

function get_al_summary(r::OneClassActiveLearning.Result, metric::Symbol, metric_short_name::String, n=5, stability_window=10)
    if !haskey(r.al_summary, metric)
        summary = (start_quality=missing, end_quality=missing, maximum=missing, ramp_up=missing, quality_range=missing, total_quality_range=missing,
                   average_end_quality=missing, average_quality_change=missing, average_gain=missing, average_loss=missing,
                   learning_stability=missing, ratio_of_outlier_queries=missing,
                   aulc=missing, reyes_paulc=missing, reyes_naulc=missing, reyes_tpr=missing, reyes_tnr=missing, reyes_tp=missing)
    else
        cur_data = r.al_summary[metric]
        start_quality = cur_data[:start_quality]
        end_quality = cur_data[:end_quality]
        maximum = cur_data[:maximum]
        ramp_up = cur_data[:ramp_up][n]
        quality_range = cur_data[:quality_range]
        total_quality_range = cur_data[:total_quality_range]
        average_end_quality = if length(cur_data[:average_end_quality]) - n + 1 > 0
                                    cur_data[:average_end_quality][end-n+1]
                              else missing end
        average_quality_change = cur_data[:average_quality_change]
        average_gain = cur_data[:average_gain]
        average_loss = cur_data[:average_loss]
        learning_stability = if length(cur_data[:learning_stability]) - stability_window + 1 > 0
                                    cur_data[:learning_stability][end-stability_window+1]
                              else missing end
        ratio_of_outlier_queries = cur_data[:ratio_of_outlier_queries]
        aulc_based_stats = if :aulc in keys(cur_data)
            aulc = cur_data[:aulc]
            reyes_paulc = cur_data[:reyes_paulc]
            reyes_naulc = cur_data[:reyes_naulc]
            reyes_tpr = cur_data[:reyes_tpr]
            reyes_tnr = cur_data[:reyes_tnr]
            reyes_tp = cur_data[:reyes_tp]
            (aulc=aulc, reyes_paulc=reyes_paulc, reyes_naulc=reyes_naulc, reyes_tpr=reyes_tpr, reyes_tnr=reyes_tnr, reyes_tp=reyes_tp)
        else
            (aulc=missing, reyes_paulc=missing, reyes_naulc=missing, reyes_tpr=missing, reyes_tnr=missing, reyes_tp=missing)
        end
        summary = (start_quality=start_quality, end_quality=end_quality, maximum=maximum, ramp_up=ramp_up, quality_range=quality_range, total_quality_range=total_quality_range,
                   average_end_quality=average_end_quality, average_quality_change=average_quality_change, average_gain=average_gain, average_loss=average_loss,
                   learning_stability=learning_stability, ratio_of_outlier_queries=ratio_of_outlier_queries, aulc_based_stats...)
    end
    return prefix_named_tuples(summary, metric_short_name)
end

function calc_stats(r::OneClassActiveLearning.Result)
    id = string(r.id)
    dir, file_name = splitdir(normpath(r.experiment[:data_file]))
    data_set = r.experiment[:data_set_name]
    split_strategy = r.experiment[:split_strategy_name][1]
    initial_pool_strategy = r.experiment[:initial_pool_strategy_name][1]
    initial_pool_resample_version = r.experiment[:param][:initial_pool_resample_version]
    initial_pool_num_labeled = sum(r.experiment[:param][:initial_pools] .!= "U")
    model = r.experiment[:model][:type]
    if model == "SSAD"
        kappa = get(r.experiment[:model][:param], :κ, 1.0)
        model = "$(model)_$(kappa)"
    end
    init_strategy_C = simplify(r.experiment[:model][:init_strategy], :init_strategy_C)
    init_strategy_gamma = simplify(r.experiment[:model][:init_strategy], :init_strategy_gamma)
    init_strategy_gamma_fitted_value = if haskey(r.experiment[:model], :fitted) && haskey(r.experiment[:model][:fitted], :kernel)
            first(match(r"\((.*)\)", r.experiment[:model][:fitted][:kernel]).captures)
        else missing end

    qs = split(string(r.experiment[:query_strategy][:type]), ".")[end]

    num_points = r.data_stats.num_observations
    num_dimensions = r.data_stats.num_dimensions
    num_al_iterations = r.experiment[:param][:num_al_iterations]

    exit_code = string(r.status[:exit_code])
    res = (id=id, file_name=file_name, data_set=data_set, split_strategy=split_strategy, initial_pool_strategy=initial_pool_strategy,
            initial_pool_resample_version=initial_pool_resample_version, initial_pool_num_labeled=initial_pool_num_labeled, model=model,
            init_strategy_C=init_strategy_C, init_strategy_gamma=init_strategy_gamma,
            init_strategy_gamma_fitted_value=init_strategy_gamma_fitted_value,
            qs=qs,
            num_points=num_points, num_dimensions=num_dimensions, num_al_iterations=num_al_iterations, exit_code=exit_code)
    stats = reduce(merge, [get_time_mem_stats(r),
            get_al_summary(r, :matthews_corr, "m"),
            get_al_summary(r, :cohens_kappa, "c"),
            get_al_summary(r, :auc, "auc"),
            get_al_summary(r, Symbol("auc_fpr_normalized_0.05"), "pauc")])
    return merge(res, stats)
end

function reduce_results(input_path)
    df = Any[]
    for scenario in readdir(input_path)
        if occursin(".csv", scenario)
            continue
        end
        for ds in readdir(joinpath(input_path, scenario, "results"))
            @info ds
            result_path = joinpath(input_path, scenario, "results", ds)
            result_files = joinpath.(result_path, filter(x -> !startswith(x, "."), readdir(result_path)))
            for r in result_files
                @info r
                res = Unmarshal.unmarshal(OneClassActiveLearning.Result, JSON.parsefile(r))
                push!(df,  merge((scenario=scenario,), calc_stats(res)))
            end
        end
    end
    df_types = mapreduce(x -> [typeof(t) for t in x], (d1, d2) -> [Union{t1,t2} for (t1,t2) in zip(d1, d2)], df)
    df_names = keys(df[1])
    res_df = DataFrame(df_types, [df_names...], 0)
    for row in df
        push!(res_df, row)
    end
    output_file = joinpath(input_path, "$(split(input_path, "/")[end]).csv")
    @info "Writing result to '$(output_file)'"
    CSV.write(output_file, res_df)
    return res_df
end

reduce_results(data_output_root)

max_stars_repo_name: kit-dbis/ocal-evaluation

max_issues_repo_name: kit-dbis/ocal-evaluation

max_forks_repo_name: kit-dbis/ocal-evaluation


text:
function preprocess!(aff::AffineSets, conic_sets::ConicSets)
    c_orig = zeros(1)
    if length(conic_sets.sdpcone) >= 1 || length(conic_sets.socone) >= 1
        all_cone_vars = Int[]
        for (idx, sdp) in enumerate(conic_sets.sdpcone)
            M = zeros(Int, sdp.sq_side, sdp.sq_side)
            iv = conic_sets.sdpcone[idx].vec_i
            im = conic_sets.sdpcone[idx].mat_i
            for i in eachindex(iv)
                M[im[i]] = iv[i]
            end
            X = LinearAlgebra.Symmetric(M, :L)

            n = size(X)[1] # columns or line
            cont = 1
            sdp_vars = zeros(Int, div(sdp.sq_side*(sdp.sq_side+1), 2))
            for j in 1:n, i in j:n
                sdp_vars[cont] = X[i, j]
                cont += 1
            end
            append!(all_cone_vars, sdp_vars)
        end
        for (idx, soc) in enumerate(conic_sets.socone)
            soc_vars = copy(soc.idx)
            append!(all_cone_vars, soc_vars)
        end

        totvars = aff.n
        extra_vars = sort(collect(setdiff(Set(collect(1:totvars)),Set(all_cone_vars))))
        ord = vcat(all_cone_vars, extra_vars)
    else
        ord = collect(1:aff.n)
    end

    c_orig = copy(aff.c)

    aff.A, aff.G, aff.c = aff.A[:, ord], aff.G[:, ord], aff.c[ord]
    return c_orig[ord], sortperm(ord)
end

function norm_scaling(affine_sets::AffineSets, cones::ConicSets)
    cte = (sqrt(2.) / 2.)
    rows = SparseArrays.rowvals(affine_sets.A)
    cont = 1
    for sdp in cones.sdpcone, j in 1:sdp.sq_side, i in 1:j
        if i != j
            for line in SparseArrays.nzrange(affine_sets.A, cont)
                affine_sets.A[rows[line], cont] *= cte
            end
        end
        cont += 1
    end
    rows = SparseArrays.rowvals(affine_sets.G)
    cont = 1
    for sdp in cones.sdpcone, j in 1:sdp.sq_side, i in 1:j
        if i != j
            for line in SparseArrays.nzrange(affine_sets.G, cont)
                affine_sets.G[rows[line], cont] *= cte
            end
        end
        cont += 1
    end
    cont = 1
    @inbounds for sdp in cones.sdpcone, j in 1:sdp.sq_side, i in 1:j
        if i != j
            affine_sets.c[cont] *= cte
        end
        cont += 1
    end
    return nothing
end


text:
# utility functions for the CUSPARSE wrappers
#
# TODO: move raw ccall wrappers to libcusparse.jl

"""
convert `SparseChar` {`N`,`T`,`C`} to `cusparseOperation_t`
"""
function cusparseop(trans::SparseChar)
    if trans == 'N'
        return CUSPARSE_OPERATION_NON_TRANSPOSE
    end
    if trans == 'T'
        return CUSPARSE_OPERATION_TRANSPOSE
    end
    if trans == 'C'
        return CUSPARSE_OPERATION_CONJUGATE_TRANSPOSE
    end
    throw(ArgumentError("unknown cusparse operation $trans"))
end

"""
convert `SparseChar` {`G`,`S`,`H`,`T`} to `cusparseMatrixType_t`
"""
function cusparsetype(mattype::SparseChar)
    if mattype == 'G'
        return CUSPARSE_MATRIX_TYPE_GENERAL
    end
    if mattype == 'T'
        return CUSPARSE_MATRIX_TYPE_TRIANGULAR
    end
    if mattype == 'S'
        return CUSPARSE_MATRIX_TYPE_SYMMETRIC
    end
    if mattype == 'H'
        return CUSPARSE_MATRIX_TYPE_HERMITIAN
    end
    throw(ArgumentError("unknown cusparse matrix type $mattype"))
end

"""
convert `SparseChar` {`U`,`L`} to `cusparseFillMode_t`
"""
function cusparsefill(uplo::SparseChar)
    if uplo == 'U'
        return CUSPARSE_FILL_MODE_UPPER
    end
    if uplo == 'L'
        return CUSPARSE_FILL_MODE_LOWER
    end
    throw(ArgumentError("unknown cusparse fill mode $uplo"))
end

"""
convert `SparseChar` {`U`,`N`} to `cusparseDiagType_t`
"""
function cusparsediag(diag::SparseChar)
    if diag == 'U'
        return CUSPARSE_DIAG_TYPE_UNIT
    end
    if diag == 'N'
        return CUSPARSE_DIAG_TYPE_NON_UNIT
    end
    throw(ArgumentError("unknown cusparse diag mode $diag"))
end

"""
convert `SparseChar` {`Z`,`O`} to `cusparseIndexBase_t`
"""
function cusparseindex(index::SparseChar)
    if index == 'Z'
        return CUSPARSE_INDEX_BASE_ZERO
    end
    if index == 'O'
        return CUSPARSE_INDEX_BASE_ONE
    end
    throw(ArgumentError("unknown cusparse index base"))
end

"""
convert `SparseChar` {`R`,`C`} to `cusparseDirection_t`
"""
function cusparsedir(dir::SparseChar)
    if dir == 'R'
        return CUSPARSE_DIRECTION_ROW
    end
    if dir == 'C'
        return CUSPARSE_DIRECTION_COL
    end
    throw(ArgumentError("unknown cusparse direction $dir"))
end

"""
check that the dimensions of matrix `X` and vector `Y` make sense for a multiplication
"""
function chkmvdims(X, n, Y, m)
    if length(X) != n
        throw(DimensionMismatch("X must have length $n, but has length $(length(X))"))
    elseif length(Y) != m
        throw(DimensionMismatch("Y must have length $m, but has length $(length(Y))"))
    end
end

"""
check that the dimensions of matrices `X` and `Y` make sense for a multiplication
"""
function chkmmdims( B, C, k, l, m, n )
    if size(B) != (k,l)
        throw(DimensionMismatch("B has dimensions $(size(B)) but needs ($k,$l)"))
    elseif size(C) != (m,n)
        throw(DimensionMismatch("C has dimensions $(size(C)) but needs ($m,$n)"))
    end
end

"""
form a `cusparseMatDescr_t` from a `CuSparseMatrix`
"""
function getDescr( A::CuSparseMatrix, index::SparseChar )
    cuind = cusparseindex(index)
    typ   = CUSPARSE_MATRIX_TYPE_GENERAL
    fill  = CUSPARSE_FILL_MODE_LOWER
    if ishermitian(A)
        typ  = CUSPARSE_MATRIX_TYPE_HERMITIAN
        fill = cusparsefill(A.uplo)
    elseif issymmetric(A)
        typ  = CUSPARSE_MATRIX_TYPE_SYMMETRIC
        fill = cusparsefill(A.uplo)
    end
    cudesc = cusparseMatDescr_t(typ, fill,CUSPARSE_DIAG_TYPE_NON_UNIT, cuind)
end

function getDescr( A::Symmetric, index::SparseChar )
    cuind = cusparseindex(index)
    typ  = CUSPARSE_MATRIX_TYPE_SYMMETRIC
    fill = cusparsefill(A.uplo)
    cudesc = cusparseMatDescr_t(typ, fill,CUSPARSE_DIAG_TYPE_NON_UNIT, cuind)
end

function getDescr( A::Hermitian, index::SparseChar )
    cuind = cusparseindex(index)
    typ  = CUSPARSE_MATRIX_TYPE_HERMITIAN
    fill = cusparsefill(A.uplo)
    cudesc = cusparseMatDescr_t(typ, fill,CUSPARSE_DIAG_TYPE_NON_UNIT, cuind)
end

# type conversion
for (fname,elty) in ((:cusparseScsr2csc, :Float32),
                     (:cusparseDcsr2csc, :Float64),
                     (:cusparseCcsr2csc, :ComplexF32),
                     (:cusparseZcsr2csc, :ComplexF64))
    @eval begin
        function switch2csc(csr::CuSparseMatrixCSR{$elty},inda::SparseChar='O')
            cuind = cusparseindex(inda)
            m,n = csr.dims
            colPtr = cuzeros(Cint, n+1)
            rowVal = cuzeros(Cint, csr.nnz)
            nzVal = cuzeros($elty, csr.nnz)
            @check ccall(($(string(fname)),libcusparse), cusparseStatus_t,
                              (cusparseHandle_t, Cint, Cint, Cint, CuPtr{$elty},
                               CuPtr{Cint}, CuPtr{Cint}, CuPtr{$elty}, CuPtr{Cint},
                               CuPtr{Cint}, cusparseAction_t, cusparseIndexBase_t),
                               handle(), m, n, csr.nnz, csr.nzVal,
                               csr.rowPtr, csr.colVal, nzVal, rowVal,
                               colPtr, CUSPARSE_ACTION_NUMERIC, cuind)
            csc = CuSparseMatrixCSC(colPtr,rowVal,nzVal,csr.nnz,csr.dims)
            csc
        end
        function switch2csr(csc::CuSparseMatrixCSC{$elty},inda::SparseChar='O')
            cuind = cusparseindex(inda)
            m,n = csc.dims
            rowPtr = cuzeros(Cint,m+1)
            colVal = cuzeros(Cint,csc.nnz)
            nzVal = cuzeros($elty,csc.nnz)
            @check ccall(($(string(fname)),libcusparse), cusparseStatus_t,
                              (cusparseHandle_t, Cint, Cint, Cint, CuPtr{$elty},
                               CuPtr{Cint}, CuPtr{Cint}, CuPtr{$elty}, CuPtr{Cint},
                               CuPtr{Cint}, cusparseAction_t, cusparseIndexBase_t),
                               handle(), n, m, csc.nnz, csc.nzVal,
                               csc.colPtr, csc.rowVal, nzVal, colVal,
                               rowPtr, CUSPARSE_ACTION_NUMERIC, cuind)
            csr = CuSparseMatrixCSR(rowPtr,colVal,nzVal,csc.nnz,csc.dims)
            csr
        end
    end
end

for (fname,elty) in ((:cusparseScsr2bsr, :Float32),
                     (:cusparseDcsr2bsr, :Float64),
                     (:cusparseCcsr2bsr, :ComplexF32),
                     (:cusparseZcsr2bsr, :ComplexF64))
    @eval begin
        function switch2bsr(csr::CuSparseMatrixCSR{$elty},
                            blockDim::Cint,
                            dir::SparseChar='R',
                            inda::SparseChar='O',
                            indc::SparseChar='O')
            cudir = cusparsedir(dir)
            cuinda = cusparseindex(inda)
            cuindc = cusparseindex(indc)
            m,n = csr.dims
            nnz = Ref{Cint}(1)
            mb = div((m + blockDim - 1),blockDim)
            nb = div((n + blockDim - 1),blockDim)
            bsrRowPtr = cuzeros(Cint,mb + 1)
            cudesca = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
            cudescc = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuindc)
            @check ccall((:cusparseXcsr2bsrNnz,libcusparse), cusparseStatus_t,
                              (cusparseHandle_t, cusparseDirection_t, Cint, Cint,
                               Ptr{cusparseMatDescr_t}, CuPtr{Cint},
                               CuPtr{Cint}, Cint, Ptr{cusparseMatDescr_t},
                               CuPtr{Cint}, Ptr{Cint}),
                              handle(), cudir, m, n, Ref(cudesca), csr.rowPtr,
                              csr.colVal, blockDim, Ref(cudescc), bsrRowPtr, nnz)
            bsrNzVal = cuzeros($elty, nnz[] * blockDim * blockDim )
            bsrColInd = cuzeros(Cint, nnz[])
            @check ccall(($(string(fname)),libcusparse), cusparseStatus_t,
                              (cusparseHandle_t, cusparseDirection_t, Cint,
                               Cint, Ptr{cusparseMatDescr_t}, CuPtr{$elty},
                               CuPtr{Cint}, CuPtr{Cint}, Cint,
                               Ptr{cusparseMatDescr_t}, CuPtr{$elty}, CuPtr{Cint},
                               CuPtr{Cint}), handle(), cudir, m, n,
                               Ref(cudesca), csr.nzVal, csr.rowPtr, csr.colVal,
                               blockDim, Ref(cudescc), bsrNzVal, bsrRowPtr,
                               bsrColInd)
            CuSparseMatrixBSR{$elty}(bsrRowPtr, bsrColInd, bsrNzVal, csr.dims, blockDim, dir, nnz[])
        end
        function switch2bsr(csc::CuSparseMatrixCSC{$elty},
                            blockDim::Cint,
                            dir::SparseChar='R',
                            inda::SparseChar='O',
                            indc::SparseChar='O')
                switch2bsr(switch2csr(csc),blockDim,dir,inda,indc)
        end
    end
end

for (fname,elty) in ((:cusparseSbsr2csr, :Float32),
                     (:cusparseDbsr2csr, :Float64),
                     (:cusparseCbsr2csr, :ComplexF32),
                     (:cusparseZbsr2csr, :ComplexF64))
    @eval begin
        function switch2csr(bsr::CuSparseMatrixBSR{$elty},
                            inda::SparseChar='O',
                            indc::SparseChar='O')
            cudir = cusparsedir(bsr.dir)
            cuinda = cusparseindex(inda)
            cuindc = cusparseindex(indc)
            m,n = bsr.dims
            mb = div(m,bsr.blockDim)
            nb = div(n,bsr.blockDim)
            nnz = bsr.nnz * bsr.blockDim * bsr.blockDim
            cudesca = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
            cudescc = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuindc)
            csrRowPtr = cuzeros(Cint, m + 1)
            csrColInd = cuzeros(Cint, nnz)
            csrNzVal  = cuzeros($elty, nnz)
            @check ccall(($(string(fname)),libcusparse), cusparseStatus_t,
                              (cusparseHandle_t, cusparseDirection_t, Cint,
                               Cint, Ptr{cusparseMatDescr_t}, CuPtr{$elty},
                               CuPtr{Cint}, CuPtr{Cint}, Cint,
                               Ptr{cusparseMatDescr_t}, CuPtr{$elty}, CuPtr{Cint},
                               CuPtr{Cint}), handle(), cudir, mb, nb,
                               Ref(cudesca), bsr.nzVal, bsr.rowPtr, bsr.colVal,
                               bsr.blockDim, Ref(cudescc), csrNzVal, csrRowPtr,
                               csrColInd)
            CuSparseMatrixCSR(csrRowPtr, csrColInd, csrNzVal, convert(Cint,nnz), bsr.dims)
        end
        function switch2csc(bsr::CuSparseMatrixBSR{$elty},
                            inda::SparseChar='O',
                            indc::SparseChar='O')
            switch2csc(switch2csr(bsr,inda,indc))
        end
    end
end

for (cname,rname,elty) in ((:cusparseScsc2dense, :cusparseScsr2dense, :Float32),
                           (:cusparseDcsc2dense, :cusparseDcsr2dense, :Float64),
                           (:cusparseCcsc2dense, :cusparseCcsr2dense, :ComplexF32),
                           (:cusparseZcsc2dense, :cusparseZcsr2dense, :ComplexF64))
    @eval begin
        function Base.Array(csr::CuSparseMatrixCSR{$elty},ind::SparseChar='O')
            cuind = cusparseindex(ind)
            m,n = csr.dims
            denseA = cuzeros($elty,m,n)
            lda = max(1,stride(denseA,2))
            cudesc = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuind)
            @check ccall(($(string(rname)),libcusparse), cusparseStatus_t,
                              (cusparseHandle_t, Cint, Cint,
                               Ptr{cusparseMatDescr_t}, CuPtr{$elty},
                               CuPtr{Cint}, CuPtr{Cint}, CuPtr{$elty}, Cint),
                               handle(), m, n, Ref(cudesc), csr.nzVal,
                               csr.rowPtr, csr.colVal, denseA, lda)
            denseA
        end
        function Base.Array(csc::CuSparseMatrixCSC{$elty},ind::SparseChar='O')
            cuind = cusparseindex(ind)
            m,n = csc.dims
            denseA = cuzeros($elty,m,n)
            lda = max(1,stride(denseA,2))
            cudesc = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuind)
            @check ccall(($(string(cname)),libcusparse), cusparseStatus_t,
                              (cusparseHandle_t, Cint, Cint,
                               Ptr{cusparseMatDescr_t}, CuPtr{$elty},
                               CuPtr{Cint}, CuPtr{Cint}, CuPtr{$elty}, Cint),
                               handle(), m, n, Ref(cudesc), csc.nzVal,
                               csc.rowVal, csc.colPtr, denseA, lda)
            denseA
        end
        function Base.Array(hyb::CuSparseMatrixHYB{$elty},ind::SparseChar='O')
            Array(switch2csr(hyb,ind))
        end
        function Base.Array(bsr::CuSparseMatrixBSR{$elty},ind::SparseChar='O')
            Array(switch2csr(bsr,ind))
        end
    end
end

for (nname,cname,rname,hname,elty) in ((:cusparseSnnz, :cusparseSdense2csc, :cusparseSdense2csr, :cusparseSdense2hyb, :Float32),
                                       (:cusparseDnnz, :cusparseDdense2csc, :cusparseDdense2csr, :cusparseDdense2hyb, :Float64),
                                       (:cusparseCnnz, :cusparseCdense2csc, :cusparseCdense2csr, :cusparseCdense2hyb, :ComplexF32),
                                       (:cusparseZnnz, :cusparseZdense2csc, :cusparseZdense2csr, :cusparseZdense2hyb, :ComplexF64))
    @eval begin
        function sparse(A::CuMatrix{$elty},fmt::SparseChar='R',ind::SparseChar='O')
            cuind = cusparseindex(ind)
            cudir = cusparsedir('R')
            if( fmt == 'C' )
                cudir = cusparsedir(fmt)
            end
            m,n = size(A)
            lda = max(1,stride(A,2))
            cudesc = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuind)
            nnzRowCol = cuzeros(Cint, fmt == 'R' ? m : n)
            nnzTotal = Ref{Cint}(1)
            @check ccall(($(string(nname)),libcusparse), cusparseStatus_t,
                              (cusparseHandle_t, cusparseDirection_t,
                               Cint, Cint, Ptr{cusparseMatDescr_t}, CuPtr{$elty},
                               Cint, CuPtr{Cint}, Ptr{Cint}), handle(),
                               cudir, m, n, Ref(cudesc), A, lda, nnzRowCol,
                               nnzTotal)
            nzVal = cuzeros($elty,nnzTotal[])
            if(fmt == 'R')
                rowPtr = cuzeros(Cint,m+1)
                colInd = cuzeros(Cint,nnzTotal[])
                @check ccall(($(string(rname)),libcusparse), cusparseStatus_t,
                                  (cusparseHandle_t, Cint, Cint,
                                   Ptr{cusparseMatDescr_t}, CuPtr{$elty},
                                   Cint, CuPtr{Cint}, CuPtr{$elty}, CuPtr{Cint},
                                   CuPtr{Cint}), handle(), m, n, Ref(cudesc), A,
                                   lda, nnzRowCol, nzVal, rowPtr, colInd)
                return CuSparseMatrixCSR(rowPtr,colInd,nzVal,nnzTotal[],size(A))
            end
            if(fmt == 'C')
                colPtr = cuzeros(Cint,n+1)
                rowInd = cuzeros(Cint,nnzTotal[])
                @check ccall(($(string(cname)),libcusparse), cusparseStatus_t,
                                  (cusparseHandle_t, Cint, Cint,
                                   Ptr{cusparseMatDescr_t}, CuPtr{$elty},
                                   Cint, CuPtr{Cint}, CuPtr{$elty}, CuPtr{Cint},
                                   CuPtr{Cint}), handle(), m, n, Ref(cudesc), A,
                                   lda, nnzRowCol, nzVal, rowInd, colPtr)
                return CuSparseMatrixCSC(colPtr,rowInd,nzVal,nnzTotal[],size(A))
            end
            if(fmt == 'B')
                return switch2bsr(sparse(A,'R',ind),convert(Cint,gcd(m,n)))
            end
            if(fmt == 'H')
                hyb = cusparseHybMat_t[0]
                @check ccall((:cusparseCreateHybMat,libcusparse), cusparseStatus_t,
                                  (Ptr{cusparseHybMat_t},), hyb)
                @check ccall(($(string(hname)),libcusparse), cusparseStatus_t,
                                  (cusparseHandle_t, Cint, Cint,
                                   Ptr{cusparseMatDescr_t}, CuPtr{$elty},
                                   Cint, CuPtr{Cint}, cusparseHybMat_t,
                                   Cint, cusparseHybPartition_t),
                                  handle(), m, n, Ref(cudesc), A, lda, nnzRowCol,
                                  hyb[1], 0, CUSPARSE_HYB_PARTITION_AUTO)
                return CuSparseMatrixHYB{$elty}(hyb[1],size(A),nnzTotal[])
            end
        end
    end
end

"""
    switch2hyb(csr::CuSparseMatrixCSR, inda::SparseChar='O')

Convert a `CuSparseMatrixCSR` to the hybrid CUDA storage format.
"""
switch2hyb(csr::CuSparseMatrixCSR, inda::SparseChar='O')

"""
    switch2hyb(csc::CuSparseMatrixCSC, inda::SparseChar='O')

Convert a `CuSparseMatrixCSC` to the hybrid CUDA storage format.
"""
switch2hyb(csc::CuSparseMatrixCSC, inda::SparseChar='O')

for (rname,cname,elty) in ((:cusparseScsr2hyb, :cusparseScsc2hyb, :Float32),
                           (:cusparseDcsr2hyb, :cusparseDcsc2hyb, :Float64),
                           (:cusparseCcsr2hyb, :cusparseCcsc2hyb, :ComplexF32),
                           (:cusparseZcsr2hyb, :cusparseZcsc2hyb, :ComplexF64))
    @eval begin
        function switch2hyb(csr::CuSparseMatrixCSR{$elty},
                            inda::SparseChar='O')
            cuinda = cusparseindex(inda)
            m,n = csr.dims
            cudesca = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
            hyb = cusparseHybMat_t[0]
            @check ccall((:cusparseCreateHybMat,libcusparse), cusparseStatus_t,
                              (Ptr{cusparseHybMat_t},), hyb)
            @check ccall(($(string(rname)),libcusparse), cusparseStatus_t,
                              (cusparseHandle_t, Cint, Cint,
                               Ptr{cusparseMatDescr_t}, CuPtr{$elty},
                               CuPtr{Cint}, CuPtr{Cint}, cusparseHybMat_t,
                               Cint, cusparseHybPartition_t), handle(),
                               m, n, Ref(cudesca), csr.nzVal, csr.rowPtr, csr.colVal,
                               hyb[1], 0, CUSPARSE_HYB_PARTITION_AUTO)
            CuSparseMatrixHYB{$elty}(hyb[1], csr.dims, csr.nnz)
        end
        function switch2hyb(csc::CuSparseMatrixCSC{$elty},
                            inda::SparseChar='O')
            switch2hyb(switch2csr(csc,inda),inda)
        end
    end
end

for (rname,cname,elty) in ((:cusparseShyb2csr, :cusparseShyb2csc, :Float32),
                           (:cusparseDhyb2csr, :cusparseDhyb2csc, :Float64),
                           (:cusparseChyb2csr, :cusparseChyb2csc, :ComplexF32),
                           (:cusparseZhyb2csr, :cusparseZhyb2csc, :ComplexF64))
    @eval begin
        function switch2csr(hyb::CuSparseMatrixHYB{$elty},
                            inda::SparseChar='O')
            cuinda = cusparseindex(inda)
            m,n = hyb.dims
            cudesca = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
            csrRowPtr = cuzeros(Cint, m + 1)
            csrColInd = cuzeros(Cint, hyb.nnz)
            csrNzVal = cuzeros($elty, hyb.nnz)
            @check ccall(($(string(rname)),libcusparse), cusparseStatus_t,
                              (cusparseHandle_t, Ptr{cusparseMatDescr_t},
                               cusparseHybMat_t, CuPtr{$elty}, CuPtr{Cint},
                               CuPtr{Cint}), handle(), Ref(cudesca),
                               hyb.Mat, csrNzVal, csrRowPtr, csrColInd)
            CuSparseMatrixCSR(csrRowPtr, csrColInd, csrNzVal, hyb.nnz, hyb.dims)
        end
        function switch2csc(hyb::CuSparseMatrixHYB{$elty},
                            inda::SparseChar='O')
            switch2csc(switch2csr(hyb,inda),inda)
        end
    end
end

"""
    switch2csr(csr::CuSparseMatrixCSR, inda::SparseChar='O')

Convert a `CuSparseMatrixCSR` to the compressed sparse column format.
"""
function switch2csc(csr::CuSparseMatrixCSR, inda::SparseChar='O') end

"""
    switch2csr(csc::CuSparseMatrixCSC, inda::SparseChar='O')

Convert a `CuSparseMatrixCSC` to the compressed sparse row format.
"""
function switch2csr(csc::CuSparseMatrixCSC, inda::SparseChar='O') end

"""
    switch2bsr(csr::CuSparseMatrixCSR, blockDim::Cint, dir::SparseChar='R', inda::SparseChar='O', indc::SparseChar='O')

Convert a `CuSparseMatrixCSR` to the compressed block sparse row format. `blockDim` sets the block dimension of the compressed sparse blocks and `indc` determines whether the new matrix will be one- or zero-indexed.
"""
function switch2bsr(csr::CuSparseMatrixCSR, blockDim::Cint, dir::SparseChar='R', inda::SparseChar='O', indc::SparseChar='O') end


text:
using IterTools

cd(dirname(@__FILE__))
lines = open("six.txt", "r") do f
    readlines(f)
end

n1 = (groupby(!isempty, lines)
    |> l -> map(join, l)
    |> l -> filter(!isempty, l)
    |> l -> map(Set, l)
    |> s -> map(length, s)
    |> sum)

println(n1)

n2 = (groupby(!isempty, lines)
    |> l -> collect(l)
    |> l -> filter(!=([""]), l)
    |> l -> map(x -> [Set(s) for s in x], l)
    |> l -> map(s -> intersect(s...), l)
    |> s -> map(length, s)
    |> sum)

println(n2)


text:
using Base.Test
using Brim

anyerrors = false

my_tests = ["recursive_brim.jl",
  "no_empty_modules.jl",
  "reorder_by_module.jl",
  "network_roles.jl",
  "partition_single.jl",
  "partition_lp.jl",
  "permutations.jl",
  "Qr.jl"]

for my_test in my_tests
  try
    include(my_test)
    println("\033[1m\033[32mPASSED\033[0m\t$(my_test)")
  catch e
    anyerrors = true
    println("\033[1m\033[31mFAILED\033[0m\t$(my_test)")
    showerror(STDOUT, e, backtrace())
    println()
  end
end

if anyerrors
  throw("Tests failed")
end


text:
struct ManifoldInterpolationData{F,uType,tType,kType,cacheType,TM} <:
       OrdinaryDiffEq.OrdinaryDiffEqInterpolation{cacheType}
    f::F
    timeseries::uType
    ts::tType
    ks::kType
    dense::Bool
    cache::cacheType
    manifold::TM
end

function (interp::ManifoldInterpolationData)(t, idxs, deriv, p, continuity)
    return ode_interpolation(t, interp, idxs, deriv, p, continuity)
end

function manifold_linear_interpolation(
    M::AbstractManifold,
    Θ,
    dt,
    p,
    q,
    idxs::Nothing,
    T::Type{Val{0}},
)
    return shortest_geodesic(M, p, q, Θ)
end

function manifold_linear_interpolation(
    M::AbstractManifold,
    Θ,
    dt,
    ps,
    qs,
    idxs,
    T::Type{Val{0}},
)
    return map(i -> shortest_geodesic(M, ps[i], qs[i], Θ), idxs)
end

function ode_interpolation(
    tval::Number,
    id::ManifoldInterpolationData,
    idxs,
    deriv,
    p,
    continuity::Symbol = :left,
)
    # implmenented based on `ode_interpolation` from OrdinaryDiffEq.jl
    @unpack ts, timeseries, ks, f, cache = id
    @inbounds tdir = sign(ts[end] - ts[1])

    if continuity === :left
        # we have i₋ = i₊ = 1 if tval = ts[1], i₊ = i₋ + 1 = lastindex(ts) if tval > ts[end],
        # and otherwise i₋ and i₊ satisfy ts[i₋] < tval ≤ ts[i₊]
        i₊ = min(lastindex(ts), OrdinaryDiffEq._searchsortedfirst(ts, tval, 2, tdir > 0))
        i₋ = i₊ > 1 ? i₊ - 1 : i₊
    else
        # we have i₋ = i₊ - 1 = 1 if tval < ts[1], i₊ = i₋ = lastindex(ts) if tval = ts[end],
        # and otherwise i₋ and i₊ satisfy ts[i₋] ≤ tval < ts[i₊]
        i₋ = max(1, OrdinaryDiffEq._searchsortedlast(ts, tval, 1, tdir > 0))
        i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
    end

    begin
        dt = ts[i₊] - ts[i₋]
        Θ = iszero(dt) ? oneunit(tval) / oneunit(dt) : (tval - ts[i₋]) / dt

        # TODO: use manifold Hermite interpolation
        val = manifold_linear_interpolation(
            id.manifold,
            Θ,
            dt,
            timeseries[i₋],
            timeseries[i₊],
            idxs,
            deriv,
        )
    end

    return val
end


text:
import BasisFunctions: discretemeasure, measure

"A generic routine to process a few arguments given to the Fun interface."
function ap_generic_options(ap::ApproximationProblem;
        samplingstyle = -1,
        solverstyle = -1,
        normalizedsampling = false,
        options...)

    if samplingstyle == -1
        samplingstyle = SamplingStyle(platform(ap))
    end
    if normalizedsampling && (samplingstyle isa DiscreteStyle)
        samplingstyle = NormalizedSampling(samplingstyle)
    end
    if solverstyle == -1
        solverstyle = SolverStyle(platform(ap), samplingstyle)
    end
    cache!(ap, SamplingStyle, samplingstyle)
    cache!(ap, SolverStyle, solverstyle)
    options
end

ap_process_options(ap::AdaptiveApproximation; options...) =
    ap_generic_options(ap; options...)

function ap_process_options(ap::PlatformApproximation; options...)
    opt = ap_generic_options(ap; options...)
    # compute and cache the sampling parameter
    samplingparameter(ap; opt...)
    opt
end

"""
    @ap_interface `function`

Generate an interface for the given function in terms of an approximation
problem. The approximation problem is constructed from the arguments of the
function.
"""
macro ap_interface(ex)
    opt = Meta.parse("options_"*string(ex))
    ret = quote
        # generate an approximation problem and pass to options parsing
        $(ex)(dict::Dictionary, args...; options...) =
            $(opt)(approximationproblem(dict, args...); options...)
        $(ex)(platform::Platform, args...; options...) =
            $(opt)(approximationproblem(platform, args...); options...)
        # Accept a function as the first argument
        $(ex)(fun, dict::Dictionary, args...; options...) =
            $(opt)(approximationproblem(fun, dict, args...); options...)
        $(ex)(fun, platform::Platform, args...; options...) =
            $(opt)(approximationproblem(fun, platform, args...); options...)
        # parse a few options and invoke the function again
        $(opt)(ap::ApproximationProblem; options...) =
            $(ex)(ap; ap_process_options(ap; options...)...)
    end
    esc(ret)
end

"Compute a property and store it in the ap cache."
macro ap_property(property, default)
    ret = quote
        $(property)(ap::ApproximationProblem; options...) = ap_hasproperty(ap, $(property)) ?
            cache(ap, $(property)) : cache!(ap, $(property), $(default))
    end
    esc(ret)
end
macro ap_property(property)
    compute = Meta.parse("compute_"*string(property))
    ret = quote
        $(property)(ap::ApproximationProblem; options...) = ap_hasproperty(ap, $(property)) ?
            cache(ap, $(property)) : cache!(ap, $(property), $(compute)(ap; options...))
    end
    esc(ret)
end

"""
Return the cached property of an approximation problem, where `property` is a
function.

If the result is not cached yet, compute it by invoking:
property(ap; options...) =
    property(SamplingStyle(ap), platform(ap), platformparameter(ap), samplingparameter(ap); options...)

Then cache the result.
"""
macro ap_sampling_property_to_platform(property)
    ret = quote
        $(property)(ap::ApproximationProblem; options...) = ap_hasproperty(ap, $(property)) ?
            cache(ap, $(property)) : cache!(ap, $(property), $(property)(SamplingStyle(ap), ap; options...))
        $(property)(sstyle::SamplingStyle, ap::ApproximationProblem; options...) =
            $(property)(sstyle, platform(ap), platformparameter(ap), samplingparameter(ap); options...)
    end
    esc(ret)
end

macro ap_sampling_property(property)
    ret = quote
        $(property)(ap::ApproximationProblem; options...) = ap_hasproperty(ap, $(property)) ?
            cache(ap, $(property)) : cache!(ap, $(property), $(property)(SamplingStyle(ap), ap; options...))
        function $(property)(ss::ProductSamplingStyle, ap::ApproximationProblem; options...)
            @assert nfactors(ss) == nfactors(platform(ap))
            to_product( map((st,pl,p_par,s_par)->$(property)(pl, p_par, s_par; samplingstyle=st, options...),
                factors(ss), factors(platform(ap)), productparameter(ap), samplingparameter(ap))...)
        end
    end
    esc(ret)
end

to_product(args...) = tensorproduct(args...)

## Basic properties

# - SamplingStyle and SolverStyle
@ap_property SamplingStyle SamplingStyle(platform(ap))
@ap_property SolverStyle SolverStyle(platform(ap), SamplingStyle(ap))

# - sampling parameter
export samplingparameter
@ap_interface samplingparameter

function samplingparameter(ap::ApproximationProblem; options...)
    ap_hasproperty(ap, samplingparameter) ? cache(ap, samplingparameter) : begin
        smpl_par = samplingparameter(SamplingStyle(ap), platform(ap), platformparameter(ap); options...)
        if platform(ap) isa ProductPlatform
            @assert length(smpl_par) == length(productparameter(ap))
        end
        cache!(ap, samplingparameter, smpl_par)
    end
end

## Other properties are implemented in separate files
include("samplingparameter.jl")
include("grids.jl")
include("samplingoperator.jl")

export discretemeasure
@ap_interface discretemeasure
@ap_property discretemeasure discretemeasure(sampling_grid(ap; options...), sampling_weights(ap; options...))

export discretization, dualdiscretization
@ap_interface discretization
@ap_property discretization apply(sampling_operator(ap), dictionary(ap); options...)

@ap_interface dualdiscretization
@ap_property dualdiscretization apply(dual_sampling_operator(ap), dualdictionary(ap); options...)

export sample_data
@ap_interface sample_data
@ap_property sample_data apply(sampling_operator(ap), ap_fun(ap); options...)

export full_discretization
@ap_interface full_discretization
@ap_property full_discretization (discretization(ap; options...), sample_data(ap; options...))

export solver
@ap_interface solver
include("solver.jl")

include("az.jl")
include("approximate.jl")

export smoothingoperator
smoothingoperator(ap::ApproximationProblem; options...) =
    WeightedSmoothingOperator(dictionary(ap); options...)


text:
 #-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*---------------------------------------------------------------------
#* Author: Estefane Horn, N�bia do Carmo Ferreira
#*$Id$
#*-------------------------------------------------------------------
#-------------------------------------------------------------------
#* Model of a valve (simplified)
#*--------------------------------------------------------------------
#*
#* Author: Paula B. Staudt
#*--------------------------------------------------------------------
type valve_simplified
	valve_simplified()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Type=>"PP"
			)),
			DanaInteger(),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream",
				:PosX=>0,
				:PosY=>0.7365,
				:Symbol=>"_{in}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>1,
				:PosY=>0.7365,
				:Symbol=>"_{out}"
			)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Plug Position"
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Fluid Density",
				:Default=>1e3
			)),
			vol_mol(Dict{Symbol,Any}(
				:Brief=>"Specific volume",
				:Default=>1e3
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop",
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Pressure Ratio",
				:Symbol=>"P_{ratio}"
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Reference Density",
				:Default=>1e4
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Valve Constant",
				:Unit=>"gal/min/psi^0.5"
			)),
			[
				:(Inlet.F = Outlet.F),
				:(Inlet.z = Outlet.z),
				:(Inlet.h = Outlet.h),
				:(Outlet.P = Inlet.P - Pdrop),
				:(Outlet.P = Inlet.P * Pratio),
				:(rho = Inlet.v*PP.VapourDensity((Inlet.T+Outlet.T)/2, (Inlet.P+Outlet.P)/2, Outlet.z) + (1-Inlet.v)*PP.LiquidDensity((Inlet.T+Outlet.T)/2, (Inlet.P+Outlet.P)/2, Outlet.z)),
				:(v = Inlet.v*PP.VapourVolume((Inlet.T+Outlet.T)/2, (Inlet.P+Outlet.P)/2, Outlet.z) + (1-Inlet.v)*PP.LiquidVolume((Inlet.T+Outlet.T)/2, (Inlet.P+Outlet.P)/2, Outlet.z)),
				:(Outlet.F * v = k*x*sqrt(Pdrop * rho_ref / rho )),
				:(Outlet.F = 0 * "kmol/h"),
			],
			[
				"Overall Molar Balance","Componente Molar Balance","Energy Balance","Pressure Drop","Pressure Ratio","Density","Volume","Flow","Closed",
			],
			[:PP,:NComp,:rho_ref,:k,],
			[:Inlet,:Outlet,:x,:rho,:v,:Pdrop,:Pratio,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Inlet::stream
	Outlet::streamPH
	x::fraction
	rho::dens_mass
	v::vol_mol
	Pdrop::press_delta
	Pratio::positive
	rho_ref::dens_mass
	k::DanaReal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export valve_simplified
function setEquationFlow(in::valve_simplified)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	if Pdrop > 0 * "atm"
		addEquation(8)
	else
		addEquation(9)
	end
end
function atributes(in::valve_simplified,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Valve"
	fields[:Brief]="Model of a very simple valve - used in distillation column models."
	fields[:Info]="== Assumptions ==
* no flashing liquid in the valve;
* the flow in the valve is adiabatic;
* dynamics in the valve are neglected;
* linear flow type.

== Specify ==
* the inlet stream
* the plug position (x) OR outlet temperature (Outlet.T) OR outlet pressure (Outlet.P)

	OR

* the inlet stream excluding its flow (Inlet.F)
* the outlet pressure (Outlet.P) OR outlet flow (Outlet.F)
* the plug position (x)
"
	drive!(fields,_)
	return fields
end
valve_simplified(_::Dict{Symbol,Any})=begin
	newModel=valve_simplified()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
import SpecialFunctions
using BenchmarkTools
import FunctionZeros: besselj_zero, besselj_zero_asymptotic

## Times reported below are all measured on the same machine
## Intel(R) Core(TM) i7-4712HQ CPU @ 2.30GHz

# This is an asymptotic form with more terms than
# the series used in ../src/. It gives no
# advantage as the initial point for root finding
# since the root finding takes far more time.
#
# This could be used for large enough n without
# root finding, especially if the user specifies
# a desired precision.
## time nu=1, n=20: 125ns
function besselj_zero_asymptotic_long0(nu::Real, n::Integer)
    beta = MathConstants.pi * (n + nu / 2 - 1//4)
    delta = 8 * beta
    mu = 4 * nu^2
    t1 = 1
    t2 = 4 * (7 * mu - 31) / (3 * delta^2)
    t3 = 32 * (84 * mu^2 - 982 * mu + 3779) / (15 * delta^4)
    t4 = 64 * (6949 * mu^3 - 153855 * mu^2 + 1585743 * mu - 6277237) /
        (105 * delta^6)
    return beta - (mu - 1) / delta * (t1 + t2 + t3 + t4)
end

# Here we compute the powers more efficiently. A Horner formula might be better.
# This is about 4.5 times faster than the version above.
## time nu=1, n=20: 29ns
function besselj_zero_asymptotic_long(nu::Real, n::Integer)
    beta = MathConstants.pi * (n + nu / 2 - 1//4)
    delta = 8 * beta
    mu = 4 * nu^2
    mup2 = mu * mu
    mup3 = mup2 * mu
    deltap2 = delta * delta
    deltap3 = deltap2 * delta
    deltap4 = deltap2 * deltap2
    deltap6 = deltap3 * deltap3
    t1 = 1
    t2 = 4 * (7 * mu - 31) / (3 * deltap2)
    t3 = 32 * (84 * mup2 - 982 * mu + 3779) / (15 * deltap4)
    t4 = 64 * (6949 * mup3 - 153855 * mup2 + 1585743 * mu - 6277237) /
        (105 * deltap6)
    zero_asymp = beta - (mu - 1)  / delta * (t1 + t2 + t3 + t4)
    return zero_asymp
end

# Here, we convert prefactors in the numerator and denominator to a single rational
## time nu=1, n=20: 316ns
function besselj_zero_asymptotic_long1(nu::Real, n::Integer)
    beta = MathConstants.pi * (n + nu / 2 - 1//4)
    delta = 8 * beta
    mu = 4 * nu^2
    mup2 = mu * mu
    mup3 = mup2 * mu
    deltap2 = delta * delta
    deltap3 = deltap2 * delta
    deltap4 = deltap2 * deltap2
    deltap6 = deltap3 * deltap3
    t1 = 1
    t2 = 4//3 * (7 * mu - 31) / deltap2
    t3 = 32//15  * (84 * mup2 - 982 * mu + 3779) /  deltap4
    t4 = 64//105 * (6949 * mup3 - 153855 * mup2 + 1585743 * mu - 6277237) /
        deltap6
    return beta - (mu - 1) / delta * (t1 + t2 + t3 + t4)
end

function test_besselj_zero_asymptotic_long(nu::Real, n::Integer)
    z = besselj_zero_asymptotic_long(nu, n)
    return SpecialFunctions.besselj(nu, z)
end

function test_besselj_zero_asymptotic(nu::Real, n::Integer)
    z = besselj_zero_asymptotic(nu, n)
    return SpecialFunctions.besselj(nu, z)
end

function test_besselj_zero(nu::Real, n::Integer)
    z = besselj_zero(nu, n)
    return SpecialFunctions.besselj(nu, z)
end


text:
using Scalpels
using LinearAlgebra
import Polynomials

function calc_vertex_of_quadratic_fit(x::AbstractVector{T1}, y::AbstractVector{T2}) where {T1<:Real, T2<:Real}
    pfit = Polynomials.fit(x, y, 2)
    @assert length(Polynomials.coeffs(pfit)) >= 3   # just in case fails to fit a quadratic
    c, b, a = Polynomials.coeffs(pfit)
    v_at_min_of_quadratic = -b/(2*a)
    return v_at_min_of_quadratic
end

function est_rv_from_test_ccf(v_grid::AbstractVector{T1}, ccf::AbstractArray{T2,1}) where {T1<:Real, T2<:Real}
    idx_at_min = findmin(ccf)[2]
    @assert 1 <= idx_at_min <= length(v_grid)
    Δ_idx_to_include_in_fit = floor(Int,length(v_grid)//8)
    idx_around_min = (idx_at_min-Δ_idx_to_include_in_fit):(idx_at_min+Δ_idx_to_include_in_fit)
    @assert 1 <= minimum(idx_around_min) < maximum(idx_around_min) <= length(v_grid)
    (min_ccf,max_ccf) = extrema(ccf)
    ys = log.((ccf./max_ccf))
    v = calc_vertex_of_quadratic_fit(v_grid,ys)
    return v
end

# Hermite functions (physicists)
H0(x) = exp(-x^2)
H1(x) = 2*x*exp(-x^2)
H2(x) = (4*x^2-2)*exp(-x^2)
H3(x) = x*(8*x^2-12)*exp(-x^2)
H4(x) = (16*x^4-48*x^2+12)*exp(-x^2)

function make_rotation_matrix(angle::Real; ndims::Integer = 2, x::Integer = 1, y::Integer =2 )
    @assert( 1<=x<=ndims)
    @assert( 1<=y<=ndims)
    @assert( x!=y)
    matrix = diagm(ones(ndims))
    matrix[x,x] = cos(angle)
    matrix[y,y] = cos(angle)
    matrix[y,x] = sin(angle)
    matrix[x,y] = -matrix[y,x]
    return matrix
end

function make_simulated_rvs_and_ccfs_to_test_scalpels(; num_rvs::Integer = 40, num_vels::Integer=200, v_max::Real = 20000,
                σ_ccf::Real = 0, depth::Real = 0.5, σv::Real = 7000, periods = [1/2, 1/3, 1/4], amplitudes = fill(0.001,length(periods)) )
    @assert length(periods) == length(amplitudes)
    @assert 1 <= length(amplitudes) <= 5
    times = range(0.0,stop=1,length=num_rvs)
    coeffs = amplitudes'.*cos.(2π.*times./periods')
    rotated_coeffs = coeffs

    if length(amplitudes)>1
        rot_angle = π/6
        rot_matrix = make_rotation_matrix(rot_angle, ndims=length(amplitudes), x=1, y=2)
        for i in 2:length(amplitudes)-1
             rot_matrix *= make_rotation_matrix(rot_angle, ndims=length(amplitudes), x=i, y=i+1)
        end
        rotated_coeffs *= rot_matrix
    end

    #println(rotated_coeffs)
    v_grid = range(-v_max,stop=v_max,length=num_vels)
    hermite_functions = [H1,H2,H3,H4]
    basis = mapreduce(h->h.(v_grid./σv), hcat, hermite_functions[1:length(amplitudes)] )

    ccfs = one(depth) .- depth .* (H0.(v_grid./σv) .+ mapreduce(i->basis[:,i].*rotated_coeffs[:,i]', +, 1:size(basis,2) ) )
    if σ_ccf > zero(σ_ccf)
        ccfs .*= one(σ_ccf) .+ σ_ccf .* randn(size(ccfs))
    end
    rvs_dirty = map(t->est_rv_from_test_ccf(v_grid,ccfs[:,t]),1:num_rvs)
    return (rvs=rvs_dirty, ccfs=ccfs)
end


text:
@generator NMCExampleGen begin
    start() = mult(a)
    a() = choose(Int, -2, 3)
    a() = choose(Float64, 9.1, 9.2)
    a() = choose(Bool)
end

@generator NMCXListGen begin
	start() = "x" * (choose(Bool) ? start() : "")
end

@generator NMCXYListGen begin
	start() = xy() * (choose(Bool) ? start() : "")
	xy() = "x"
	xy() = "y"
end

@testset "NMCS choice model" begin

	@testset "constructor (one level)" begin

	    gn = NMCExampleGen()

		setchoicemodel!(gn, SamplerChoiceModel(gn))
		setchoicemodel!(gn, NMCSChoiceModel(gn, x->1.0))
		cm = choicemodel(gn)

	    @test typeof(cm) == DataGenerators.NMCSChoiceModel

		setchoicemodel!(gn, NMCSChoiceModel(gn, x->1.0))
		cm2 = choicemodel(gn)

	    @test typeof(cm2) == DataGenerators.NMCSChoiceModel

	end

	@testset "set/get parameters and ranges" begin

	    gn = NMCExampleGen()

		setchoicemodel!(gn, SamplerChoiceModel(gn))
		scm = deepcopy(choicemodel(gn))
		setchoicemodel!(gn, NMCSChoiceModel(gn, x->1.0))
		cm = choicemodel(gn)

		@testset "paramranges" begin
		    ranges = paramranges(cm)
		    @test typeof(ranges) <: Vector{Tuple{Float64,Float64}}
		    @test ranges == paramranges(scm)
		end

		@testset "getparams" begin
		    params = getparams(cm)
		    @test typeof(params) <: Vector{Float64}
		    @test params == getparams(scm)
		end

		@testset "setparams" begin
		    newparams = [0.7, 0.4, 0.7, 0.4, 0.7]
			setparams!(scm, newparams)
		    setparams!(cm, newparams)
		    @test getparams(cm) == getparams(scm)
		end

	end

	@testset "always successful with sufficient samplesize" begin

		gn = NMCXListGen()
		setchoicemodel!(gn, SamplerChoiceModel(gn))
		fitness(x) = abs(length(x)-20)
		setchoicemodel!(gn, NMCSChoiceModel(gn, fitness, 10))

		@mtestset "length is always 20" begin
			y = choose(gn)
			@test length(y) == 20
		end

	end

	@testset "sometimes unsuccessful with insufficient samplesize" begin

		gn = NMCXListGen()
		setchoicemodel!(gn, SamplerChoiceModel(gn))
		fitness(x) = abs(length(x)-20)
		setchoicemodel!(gn, NMCSChoiceModel(gn, fitness, 1))

		@mtestset "length is sometimes not 20" begin
			y = choose(gn)
			@mtest_that_sometimes length(y) != 20
		end

	end

	@testset "sometimes successful with inbetween samplesize" begin

		gn = NMCXListGen()
		setchoicemodel!(gn, SamplerChoiceModel(gn))
		fitness(x) = abs(length(x)-20)
		setchoicemodel!(gn, NMCSChoiceModel(gn, fitness, 2))

		@mtestset "length is sometimes 20" begin
			y = choose(gn)
			@mtest_that_sometimes length(y) == 20
		end

	end

	@testset "state is reset on each choice" begin

		gn = NMCXYListGen()
		setchoicemodel!(gn, SamplerChoiceModel(gn))
		fitness(x) = abs(length(x)-20)
		setchoicemodel!(gn, NMCSChoiceModel(gn, fitness, 10))

		@mtestset "length is always 40 but values are different" begin
			y = choose(gn)
			@test length(y) == 20
			@mtest_values_vary y
		end

	end

end