text:
################################################################################
# Breast Cancer diagnosis problem
#
# In this problem we are given a dataset containing real world characteristics of observed breast cancer (size, pattern,..) together with the associated diagnosis in terms of malignity or benignity of the cancer.
# Our task is to build a linear classifier using the perceptron algorithm that we studied and train it in order to make diagnosis based on the cancer characteristics.

# ### Environment set-up and data loading

# 1) Start by setting the working directory to the directory of this file and activate it. If you have the provided `Manifest.toml` file in the directory, just run `Pkg.instantiate()`, otherwise manually add the packages Pipe, HTTP, StatsPlots and BetaML.
# Also, seed the random number generator with the integer `123`.

# 2) Load the packages/modules Statistics, DelimitedFiles, LinearAlgebra, Pipe, HTTP, StatsPlots, BetaML

# 3) Load from internet or from localfile the input data and shuffle its rows (records)

dataURL = "https://raw.githubusercontent.com/sylvaticus/IntroSPMLJuliaCourse/main/lessonsMaterial/03_ML1/BreastCancerDiagnosisWithPerceptron/data/wdbc.data.csv"

# Source: [Breast Cancer Wisconsin (Diagnostic) Data Set, UCI Machine Learning Repository](https://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+%28Diagnostic%29)

# 4) Map the data to (X,y)
# The data you have loaded contains the actual diagnosis for the cancer in the second column, coded with a string "B" for "Benign" and "M" for "Malign", and the characteristics of the cancer foir the next 30 columns.
# Save the diagnosis to the vector `y`, coding malign cancers with `1` and benign cancers with `-1`
# Save the characteristics to the feature matrix `X` (and be sure it is made of Float64)

# 5) (this task is provided) Plot the first 2 attributes of training points and define the function plot2DClassifierWithData()
colors = [y == 1 ? "red" : "green" for y in y]
labels = [y == 1 ? "malign" : "benign" for y in y]
scatter(X[:,1],X[:,2], colour=colors, title="Classified tumors",xlabel="Tumor Radius", ylabel="Tumor Texture", group=labels)

function plot2DClassifierWithData(X,y,θ;d1=1,d2=2,origin=false,xlabel="Dimx: $(d1)",ylabel="Dimy: $(d2)")
    nR    = size(X,1)
    X     = hcat(ones(nR),X)
    X     = scale(X) # for visualisation
    d1    += 1
    d2    += 1
    colors = [y == 1 ? "red" : "green" for y in y]
    labels = [y == 1 ? "malign" : "benign" for y in y]
    minD1,maxD1 = extrema(X[:,d1])
    minD2,maxD2 = extrema(X[:,d2])
    myplot = scatter(X[:,d1],X[:,d2], colour=colors, title="Linear classifier in 2D",xlabel=xlabel, ylabel=ylabel, group=labels)
    d2Class(x) = -θ[1]/θ[d2] -x * θ[d1]/θ[d2]
    if θ[d2] == 0
        vline!([0], color= "blue",label="",linewidth=5)
    else
        plot!(d2Class,minD1,maxD1, color= "blue",label="",linewidth=5)
    end
    display(myplot)
end

# 6) (provided) Define the Model and training options structures

abstract type SupervisedModel end
abstract type TrainingOptions end

mutable struct Perceptron <: SupervisedModel
    θ::Vector{Float64}
end

mutable struct PerceptronTrainingOptions <: TrainingOptions
    epochs::Int64
    verbose::Bool
    shuffle::Bool
    function PerceptronTrainingOptions(;epochs=1,verbose=false,shuffle=false)
        return new(epochs,verbose,shuffle)
    end
end

# 7) (provided) Implement the functions `predict()` and `update()`
function predict(model::Perceptron,x::AbstractVector)
    x = vcat(1.0,x)
    x' * model.θ > eps() ? (return 1) : (return -1)
end

function predict(model::Perceptron,X::AbstractMatrix)
    return [predict(model,r) for r in eachrow(X)]
end

function update!(model::Perceptron,X::Vector,y)
    X       = vcat(1.0,X)
    model.θ = model.θ .+ y .* X
    return model.θ
end

# 8) Implement the function `train!(model::Perceptron,X,y,ops=PerceptronTrainingOptions()::TrainingOptions)`
# Compared to the function we saw in the 0302-perceptron.jl file, add, if you wish, a counter to eventually return early if there are no more errors in an epoch (i.e., all points are correctly classified)
function train!(model::Perceptron,X,y,ops=PerceptronTrainingOptions()::TrainingOptions)
    #...
    for t in 1:epochs
        #...
        if ops.shuffle
           #...
        end
        for n in 1:nR
           #...
        end
        #...
    end
    #...
    return model.θ
end

# 9) Instanziate a `Perceptron` object with a parameter vector of nD+1 zeros and a `PerceptronTrainingOption` object with 5 epochs and shuffling, use the options to train the model on the whole dataset, compute the model predictions and the accuracy relative to the whole sample.

m   = Perceptron(zeros(size(X,2)+1))
ops = #...
train!(m,X,y,ops)
plot2DClassifierWithData(X,y,m.θ,d1=1,d2=2,xlabel="Tumor Radius", ylabel="Tumor Texture")
ŷ           = #...
inSampleAcc = accuracy(#= ... =#) # 0.91

# 10) Partition the data in `(xtrain,xtest)` and `(ytrain,ytest)` keeping 65% of the data for training and reserving 35% for testing

((xtrain,xtest),(ytrain,ytest)) = partition(#=...=#)

# 11) Using a 10-folds cross-validation strategy, find the best hyperparameters within the following ranges :

sampler = KFold(#=...=#)

epochsSet  = 1:5:150
shuffleSet = [false,true]

bestE       = 0
bestShuffle = false
bestAcc     = 0.0

accuraciesNonShuffle = []
accuraciesShuffle    = []

for e in epochsSet, s in shuffleSet
    global bestE, bestShuffle, bestAcc, accuraciesNonShuffle, accuraciesShuffle
    local acc
    local ops  = PerceptronTrainingOptions(#=...=#)
    (acc,_)    = crossValidation([xtrain,ytrain],sampler) do trainData,valData,rng
                    (xtrain,ytrain) = trainData; (xval,yval) = valData
                    m               = Perceptron(zeros(size(xtrain,2)+1))
                    train!(#=...=#)
                    ŷval            = predict(#=...=#)
                    valAccuracy     = accuracy(#=...=#)
                    return valAccuracy
                end
    if acc > bestAcc
        bestAcc     = acc
        bestE       = e
        bestShuffle = s
    end
    if s
        push!(accuraciesShuffle,acc)
    else
        push!(accuraciesNonShuffle,acc)
    end
end

bestAcc # 0.91
bestE
bestShuffle

plot(epochsSet,accuraciesNonShuffle,label="Val accuracy without shuffling", legend=:bottomright)
plot!(epochsSet,accuraciesShuffle, label="Val accuracy with shuffling")

# 12) Using the "best" hyperparameters found in the previous step, instantiate a new model and options, train the model using `(xtrain,ytrain)`, make your predicitons for the testing features (`xtest`) and compute your output accuracy compared with those of the true `ytest` (use the BetaML function `accuracy`)

ops = PerceptronTrainingOptions(#=...=#)
m   = Perceptron(zeros(size(xtest,2)+1))
train!(#=...=#)
ŷtest           = predict(#=...=#)
testAccuracy    = accuracy(#=...=#) # 0.89

plot2DClassifierWithData(xtest,ytest,m.θ,xlabel="Tumor Radius", ylabel="Tumor Texture")
plot2DClassifierWithData(xtest,ytest,m.θ,d1=3,d2=4)

# 13) (optional) Add a scaling passage to the workflow and test it with cross-validation if it improves the accuracy

epochsSet  = 1:10:150
shuffleSet = [false,true]
scalingSet = [false,true]

bestE       = 0
bestShuffle = false
bestScaling = false
bestAcc     = 0.0

accNShuffleNSc = Float64[]
accNShuffleSc  = Float64[]
accShuffleNSc  = Float64[]
accShuffleSc   = Float64[]

scfactors = getScaleFactors(xtrain)

for e in epochsSet, s in shuffleSet, sc in scalingSet
    global bestE, bestShuffle, bestAcc, accNShuffleNSc, accNShuffleSc, accShuffleNSc, accShuffleSc
    local acc
    local ops  = PerceptronTrainingOptions(#=...=#)
    xtrainsc= copy(xtrain)
    if(sc)
        xtrainsc =scale(xtrain,scfactors)
    end
    (acc,_)    = crossValidation([xtrainsc,ytrain],sampler) do trainData,valData,rng
                    #...
                    return valAccuracy
                end
    if acc > bestAcc
        bestAcc     = acc
        bestE       = e
        bestShuffle = s
        bestScaling = sc
    end
    if s && sc
        push!(accShuffleSc,acc)
    elseif  s && !sc
        push!(accShuffleNSc,acc)
    elseif  !s && sc
        push!(accNShuffleSc,acc)
    elseif  !s && !sc
        push!(accNShuffleNSc,acc)
    else
        @error "Something wrong here"
    end
end

bestAcc # 0.97
bestE
bestShuffle
bestScaling

plot(epochsSet,accShuffleSc,label="Val accuracy shuffling scaling", legend=:bottomright)
plot!(epochsSet,accNShuffleSc,label="Val accuracy NON shuffling scaling", legend=:bottomright)
plot!(epochsSet,accShuffleNSc,label="Val accuracy shuffling NON scaling", legend=:bottomright)
plot!(epochsSet,accNShuffleNSc,label="Val accuracy NON shuffling NON scaling", legend=:bottomright)

ops = PerceptronTrainingOptions(#=...=#)
m   = Perceptron(#=...=#)
if bestScaling
    train!(m,scale(xtrain),ytrain,ops)
    ŷtest  = predict(m,scale(xtest))
else
    train!(m,xtrain,ytrain,ops)
    ŷtest  = predict(m,xtest)
end
testAccuracy    = accuracy(#=...=#) # 0.96

plot2DClassifierWithData(xtest,ytest,m.θ,xlabel="Tumor Radius", ylabel="Tumor Texture")

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
# layout-related functions

ContViews{T,N} =  Union{ContiguousView{T,N},UnsafeContiguousView{T,N}}
NonContViews{T,N,M} =  Union{StridedView{T,N,M}, UnsafeStridedView{T,N,M}}

ContiguousArray{T,N} =  Union{Array{T,N}, ContiguousView{T,N}}
ContiguousVector{T} =  ContiguousArray{T,1}
ContiguousMatrix{T} =  ContiguousArray{T,2}

import Base: IndexStyle, IndexLinear, IndexCartesian
IndexStyle(a::ContViews) = IndexLinear()
IndexStyle(a::NonContViews) = IndexCartesian()
IndexStyle(a::NonContViews{T,1}) where {T} = IndexLinear()

## strides method

strides(a::ContViews{T,1}) where {T} = (1,)
strides(a::ContViews{T,2}) where {T} = (1, a.shp[1])
strides(a::ContViews{T,3}) where {T} = ((d1, d2, d3) = size(a); (1, d1, d1 * d2))

strides(a::ContViews{T,4}) where {T} = ((d1, d2, d3, d4) = size(a); (1, d1, d1 * d2, d1 * d2 * d3))

function strides(a::ContViews{T,5}) where T
    (d1, d2, d3, d4, d5) = a.shp
    s3 = d1 * d2
    s4 = s3 * d3
    s5 = s4 * d4
    (1, d1, s3, s4, s5)
end

strides(a::NonContViews) = a.strides

## stride method

stride(a::ContiguousView{T,1}, d::Integer) where {T} =
    (d > 0 || error("dimension out of range.");
     d == 1 ? 1 : length(a))::Int

stride(a::ContiguousView{T,2}, d::Integer) where {T} =
    (d > 0 || error("dimension out of range.");
     d == 1 ? 1 :
     d == 2 ? a.shp[1] : length(a))::Int

stride(a::ContiguousView{T,3}, d::Integer) where {T} =
    (d > 0 || error("dimension out of range.");
     d == 1 ? 1 :
     d == 2 ? a.shp[1] :
     d == 3 ? a.shp[1] * a.shp[2] : length(a))::Int

 stride(a::ContiguousView{T,4}, d::Integer) where {T} =
    (d > 0 || error("dimension out of range.");
     d == 1 ? 1 :
     d == 2 ? a.shp[1] :
     d == 3 ? a.shp[1] * a.shp[2] :
     d == 4 ? a.shp[1] * a.shp[2] * a.shp[3] : length(a))::Int

stride(a::ContiguousView{T,5}, d::Integer) where {T} =
    (d > 0 || error("dimension out of range.");
     d == 1 ? 1 :
     d == 2 ? a.shp[1] :
     d == 3 ? a.shp[1] * a.shp[2] :
     d == 4 ? a.shp[1] * a.shp[2] * a.shp[3] :
     d == 5 ? a.shp[1] * a.shp[2] * a.shp[3] * a.shp[4] : length(a))::Int

stride(a::NonContViews{T,N}, d::Integer) where {T,N} =
    (d > 0 || error("dimension out of range.");
     d <= N ? a.strides[d] : a.strides[N])

 # test contiguousness

iscontiguous(a::ContViews) = true

iscontiguous(a::NonContViews{T,N}) where {T,N} = _iscontiguous(size(a), strides(a))

_iscontiguous(shp::NTuple{0,Int}, strides::NTuple{0,Int}) = true
_iscontiguous(shp::NTuple{1,Int}, strides::NTuple{1,Int}) = (strides[1] == 1)
_iscontiguous(shp::NTuple{2,Int}, strides::NTuple{2,Int}) = (strides[1] == 1 && strides[2] == shp[1])
_iscontiguous(shp::NTuple{3,Int}, strides::NTuple{3,Int}) =
    (strides[1] == 1 &&
     strides[2] == shp[1] &&
     strides[3] == shp[1] * shp[2])

function _iscontiguous(shp::NTuple{4,Int}, strides::NTuple{4,Int})
    s2 = shp[1]
    s3 = s2 * shp[2]
    s4 = s3 * shp[3]
    strides[1] == 1 && strides[2] == s2 && strides[3] == s3 && strides[4] == s4
end

function _iscontiguous(shp::NTuple{5,Int}, strides::NTuple{5,Int})
    s2 = shp[1]
    s3 = s2 * shp[2]
    s4 = s3 * shp[3]
    s5 = s4 * shp[4]
    strides[1] == 1 && strides[2] == s2 && strides[3] == s3 && strides[4] == s4 && strides[5] == s5
end

ext: jl


text:
using Plots,ApproxFun

x=Fun()
D=Derivative()

L=D^10 + cosh(x)*D^8 + x^3*D^6 + x^4*D^4 + cos(x)*D^2 + x^2
d=Interval()
B=[dirichlet(d);neumann(d);[Evaluation(Interval(),false,k) for k=2:4];
           [Evaluation(Interval(),true,k) for k=2:4]]

u=[B;L]\[0.;0.;1.;1.;zeros(6);exp(x)]

ApproxFun.plot(u)


text:
module LandscapeChange

using Reexport

@reexport using DynamicGrids

using DynamicGrids.Setfield
using DynamicGrids.StaticArrays

export LandCoverCount, LandCoverPotential, LandCoverSuitability, LandCoverChange

export NeutralLUCModel, RandomConstraintMatch, GrowingClusters

export SuitabilityScore, MaxProbability, ChangeProbabilities, RawSuitabilities

export PreallocatedUnorderedList, NamedVector

export sim

include("unordered_list.jl")
include("named_vector.jl")
include("rules.jl")
include("neutral.jl")
include("scores.jl")

end


text:
include("header.jl")

#=
# Isentropic quasi-1D steady compressible flow
This notebook demonstrates the use of the compressible flow tools for computing
states and processes in isentropic quasi-1d steady compressible flow.
=#

# ### Set up the module
#-
using MAE103
#-
using Plots

#=
### Using isentropic relations
Let us apply some basic isentropic relations.
=#

#=
#### Example 1
Let's suppose two states of the flow in air, 1 and 2, are connected with each
other by an isentropic process. We know the pressure p1 is 101 kPa, the temperature
T1 is 20 degrees C, and the pressure p2 is 80 kPa. What is temperature T2?

To answer this, we will use the relationship

$$\dfrac{T_2}{T_1} = \left( \dfrac{p_2}{p_1}\right)^{(\gamma-1)/\gamma}$$
=#
p1 = Pressure(101u"kPa")
T1 = Temperature(20u"°C")
p2 = Pressure(80u"kPa")

#=
First, let's set the pressure ratio:
=#
p2_over_p1 = PressureRatio(p2/p1)

#=
Now find the temperature ratio. Note below that we specify the argument
`Isentropic` to make sure it is clear that we are using the isentropic relation.
We only need this argument when it is needed for clarity.
=#
T2_over_T1 = TemperatureRatio(p2_over_p1,Isentropic)

#=
It is important to understand that the tools "know" what formula you want to use
(the one listed above), based on the fact that (a) you supplied it with a pressure
ratio (the purpose of the `PressureRatio` line above), and (b) you told it that
the process is `Isentropic`. It figures out the rest.
=#
#=
Finally, calculate $T_2 = T_1 (T_2/T_1)$:
=#
T2 = Temperature(T1*T2_over_T1)

#=
or, in Celsius, if desired
=#
value(T2,u"°C")

#=
We could also do all of this in one line, though it is a bit harder to debug
if something goes wrong:
=#
T2 = Temperature(T1*TemperatureRatio(PressureRatio(p2/p1),Isentropic))

#=
#### Example 2
If the temperature ratio $T/T_0$ is 0.2381, what is the Mach number?
=#
MachNumber(TemperatureRatio(0.2381),Isentropic)

#=
#### Example 3
If the Mach number is 4.4 and stagnation pressure is 800 kPa, what is the pressure?
=#
Pressure(StagnationPressure(800u"kPa"),MachNumber(4.4),Isentropic)

#=
### Mach - area relations
A big part of isentropic quasi-1D flow deals with changes of the flow in variable-area
ducts. For these calculations, we make use of the *sonic area* $A_*$ as a reference area.
Remember that, for any ratio of $A/A_*$, there are two possible Mach numbers,
corresponding to a **subsonic flow** and a **supersonic flow**. Let us see that by
plotting $A/A_*$ versus Mach number $M$. If, for example, $A/A_* = 2$, then
note where the dashed line crosses the plot:
=#
Mrange = range(0,6,length=601)
Aratio = []
for M in Mrange
    push!(Aratio,value(AOverAStar(MachNumber(M),Isentropic)))
end

plot(Mrange,Aratio,xlim=(0,4),ylim=(0,12),yticks=0:1:12,xlabel="Mach number",ylabel=L"A/A_*",legend=false)
scatter!([1],[1])
plot!(Mrange,2*ones(length(Mrange)),style=:dash)

#=
#### Example 4
What are the subsonic and supersonic Mach numbers associated with, for example,
an area ratio $A/A_*$ of 2 (the dashed line in the plot above)? We simply input
the desired area ratio and find the two solutions:
=#
M1, M2 = MachNumber(AreaRatio(2),Isentropic);
#-
M1
#-
M2

#=
A related question: What is the local sonic reference area $A_*$ when the Mach
number is 7.1 and the local area is 50 sq cm?

We first compute $A/A_*$ from $M = 7.1$, then compute $A_* = A/(A/A_*)$:
=#
A = Area(50u"cm^2")
M = MachNumber(7.1)
A_over_Astar = AOverAStar(M,Isentropic)
Astar = Area(A/A_over_Astar)
value(Astar,u"cm^2")
#=
So the throat would have to be 0.45 sq cm, much smaller than 50 sq cm!
=#
#=
Note that there is a convenience function to do those steps all in one:
=#
value(AStar(A,M),u"cm^2")

#=
#### Example 5
Consider the flow of air through a converging-diverging nozzle, leaving a stagnant
reservoir at pressure $p_0$ = 700 kPa and temperature $T_0 = 30$ degrees C. The Mach
number at a location (1) in the converging section with area 50 sq cm is equal to 0.4.
The exit of the nozzle has area 60 sq cm.

(a) What are the possible Mach numbers at the exit in choked isentropic conditions?
What are the exit (i.e., "back") pressures $p_2$ associated with these two Mach numbers?

(b) What is the mass flow rate through the nozzle in these conditions?
=#

# First, we set the known values
p0 = StagnationPressure(700u"kPa")
T0 = StagnationTemperature(30u"°C")
A1 = Area(50u"cm^2")
A2 = Area(60u"cm^2")
M1 = MachNumber(0.4)

#=
Now, we compute $A_1/A_*$ from $M_1$. Use this to calculate $A_*$ from $A_1/(A_1/A_*)$.
Then calculate $A_2/A_*$:
=#
A1_over_Astar = AOverAStar(M1,Isentropic,gas=Air)
Astar = Area(A1/A1_over_Astar)
A2_over_Astar = AreaRatio(A2/Astar)

#=
Now calculate the Mach numbers at location 2 (the nozzle exit):
=#
M2sub = SubsonicMachNumber(A2_over_Astar,Isentropic,gas=Air);
#-
M2sup = SupersonicMachNumber(A2_over_Astar,Isentropic,gas=Air);
#=
Actually, all of the last few steps can be done in *one step* with a different
version of the function `MachNumber`:
=#
M2sub, M2sup = MachNumber(M1,A1,A2,Isentropic,gas=Air);

#=
Now let's determine the exit pressures (location 2) corresponding to these two Mach numbers:
=#
p0_over_p2sub = P0OverP(M2sub,Isentropic,gas=Air)
p2sub = Pressure(p0/p0_over_p2sub)
value(p2sub,u"kPa")
#-
p0_over_p2sup = P0OverP(M2sup,Isentropic,gas=Air)
p2sup = Pressure(p0/p0_over_p2sup)
value(p2sup,u"kPa")

#=
So if the exit pressure is 651 kPa, then the flow will remain choked and **subsonic**
throughout, and if the exit pressure is 71.5 kPa, then the flow will remain choked
and **supersonic** throughout.

Note that, because the flow is choked, the mass flow rate is the same for both of
these cases. Let's calculate that mass flow rate (using $\rho_1 u_1 A_1$). We need
$\rho_1$ and $u_1$. First, calculate the stagnation density in this nozzle:
=#
ρ0 = StagnationDensity(p0,T0,gas=Air) # this uses the perfect gas law

# Using $M_1$, find $\rho_0/\rho_1$:
ρ0_over_ρ1 = ρ0Overρ(M1,Isentropic,gas=Air)

# So we can get $\rho_1$ from $\rho_0/(\rho_0/\rho_1)$:
ρ1 = Density(ρ0/ρ0_over_ρ1)

# Now get $u_1$ from $M_1 c_1$, where $c_1$ is the speed of sound. For that, we need $T_1$:
T0_over_T1 = T0OverT(M1,Isentropic,gas=Air)
T1 = Temperature(T0/T0_over_T1)
c1 = SoundSpeed(T1)
#-
u1 = Velocity(M1*c1)

# And now we can put it all together:
mdot = MassFlowRate(ρ1*u1*A1)

# So the mass flow rate is stuck at 5.1 kg/s


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

using Markdown
using InteractiveUtils

# ╔═╡ 6ec070b0-7214-11eb-04d3-a54e83e79e71
using GMT

# ╔═╡ 411f1300-7214-11eb-0cf5-7b1317346330
md"# Plotyy example
Likewise *contourf* GMT does not have a *plotyy* module. A `plotyy` plot is a one where we plot Y1 versus X1 with y-axis labeling on the left and plots Y2 versus X2 with y-axis labeling on the right. So it's basically a two plots overlain but the details to make it nice can be a bit boring and long.
"

# ╔═╡ 73ff5dc0-7214-11eb-21db-232604c38428
begin
	x = 0:0.01:20;
	y1 = 200 * exp.(-0.05x) .* sin.(x);
	y2 = 0.8 * exp.(-0.5x)  .* sin.(10x);
end;

# ╔═╡ 967ba4d0-7214-11eb-10a9-23ca5705e4ea
plotyy(y1, y2, title="Vibrating dishes", ylabel=:Knifes, xlabel=:Forks, seclabel=:Spoons, show=true)

# ╔═╡ 571727a0-721f-11eb-0d87-dbe6673f8e03
md"Note that to make the command shorter and nicer to read we have used a less known option in GMT. The *secondary* label of an axes. In this example we also didn't set the *xx* coordinates so the program plotted from 1 to numbers of points.

In the general case the data has *xx* coordinates and they don't even need to be the same for Y1 and Y2 (but they need to have a shared interval). In that case we **should** set the plot limits because otherwise the guessing done from *xx*,Y1 risk to not capture the total Y1+Y2 extent.
"

# ╔═╡ 320ffa82-7220-11eb-1c74-b1f9e7d3e8a0
plotyy([x[:] y1[:]], [x[:] y2[:]], title="Vibrating dishes", ylabel=:Knifes, xlabel="2000 Forks", seclabel=:Spoons, show=true)

# ╔═╡ Cell order:
# ╟─411f1300-7214-11eb-0cf5-7b1317346330
# ╠═6ec070b0-7214-11eb-04d3-a54e83e79e71
# ╠═73ff5dc0-7214-11eb-21db-232604c38428
# ╠═967ba4d0-7214-11eb-10a9-23ca5705e4ea
# ╟─571727a0-721f-11eb-0d87-dbe6673f8e03
# ╠═320ffa82-7220-11eb-1c74-b1f9e7d3e8a0


text:
importall MathProgBase.SolverInterface

export Presolver

type Presolver <: AbstractMathProgSolver
    realsolver :: AbstractMathProgSolver
    # Presolver() gives an error. says get_defualt_solver undefined, why?
    function Presolver()
        realsolver = get_defualt_solver()
        new(realsolver)
    end
    function Presolver(solver::AbstractMathProgSolver)
        realsolver = solver
        new(realsolver)
    end
end

type PresolveMathProgModel <: AbstractMathProgModel
    m::Int
    n::Int
    A::SparseMatrixCSC{Float64,Int}
    b::Vector{Float64}
    c::Vector{Float64}
    p::Presolve_Problem
    orig_sense::Symbol
    solve_stat::Symbol
    obj_val::Float64
    primal_sol::Array{Float64,1}
    dual_sol::Array{Float64,1}
    innermodel :: AbstractMathProgModel
    realsolver #gives an error if I specify the type. Can't converr from SCSSolver to AbstractMathProgSolver
    function PresolveMathProgModel()
        realsolver = get_default_solver()
        innermodel = LinearQuadraticModel(realsolver)
        p = Presolve_Problem()
        new(0,0,spzeros(0,0),Int[],Int[],:Min,p,:NotSolved,0.0,Int[],model,realsolver)
    end
    function PresolveMathProgModel(solver::AbstractMathProgSolver)
        realsolver=solver
        innermodel = LinearQuadraticModel(realsolver)
        p = Presolve_Problem()
        new(0,0,spzeros(0,0),Int[],Int[],:Min,p,:NotSolved,0.0,Int[],model,realsolver)
    end
end

function ConicModel(s::Presolver)
    return PresolveMathProgModel(s.realsolver)
end

LinearQuadraticModel(s::Presolver) = ConicToLPQPBridge(ConicModel(s))

status(model::PresolveMathProgModel) = model.solve_stat
getobjval(model::PresolveMathProgModel) = model.obj_val
getsolution(model::PresolveMathProgModel) = copy(model.primal_sol)
numvar(model::PresolveMathProgModel) = model.m  # my m doesn't change. in future it might.
numconstr(model::PresolveMathProgModel) = model.n
supportedcones(s::Presolver) = [:Free, :Zero, :NonNeg, :NonPos]

loadproblem!(model::PresolveMathProgModel,c,A,b,constr_cones,var_cones)
    = loadproblem!(model,c,sparse(A),b,constr_cones,var_cones)

function loadproblem!(model::PresolveMathProgModel, c, A::SparseMatrixCSC, b, constr_cones, var_cones)
    bad_cones = [:SOC,:SOCRotated,:SDP,:ExpPrimal,:ExpDual]
    for cone_vars in constr_cones
        cone_vars[1] in bad_cones && error("Cone type $(cone_vars[1]) is currently not supported")
        cone_vars[1] == :Free && error("There is a free constraint")
    end
    for cone_vars in var_cones
        cone_vars[1] in bad_cones && error("Cone type $(cone_vars[1]) is currently not supported")
    end

    # going to be used later. c_cones, and v_cones?
    cone_dict = Dict{Symbol,Int64}()
    cone_dict[:Zero] = 1
    cone_dict[:NonNeg] = 2
    cone_dict[:NonPos] = 3
    cone_dict[:Free] = 4

    c_cones = (Array{Int,1}(),Array{Int,1}(),Array{Int,1}(),Array{Int,1}())
    for (cone,id) in constr_cones
        if(isa(id,Number))
           push!(c_cones[cone_dict[c]],id)
        else
           append!(c_cones[cone_dict[c]],id)
        end
        push!(c_cones[cone_dict[cone]],id)
    end

    v_cones = (Array{Int,1}(),Array{Int,1}(),Array{Int,1}(),Array{Int,1}())
    for (cone,id) in var_cones
        if(isa(id,Number))
           push!(v_cones[cone_dict[c]],id)
        else
           append!(v_cones[cone_dict[c]],id)
        end
        push!(v_cones[cone_dict[cone]],id)
    end

    # Column or Variable bounds
    collb = fill(-Inf, length(c))
    colub = fill(Inf, length(c))
    for (cone,idxs) in var_cones
        l = (cone == :Free || cone == :NonPos) ? -Inf : 0.0
        u = (cone == :Free || cone == :NonNeg) ?  Inf : 0.0
        for idx in idxs
            collb[idx] = l
            colub[idx] = u
        end
    end

    # Row or Constraint bounds
    rowlb = Array(Float64,length(b))
    rowub = Array(Float64,length(b))
    for (cone,idxs) in constr_cones
        for idx in idxs
            rowlb[idx] = (cone == :Zero || cone == :NonPos) ? b[idx] : -Inf
            rowub[idx] = (cone == :Zero || cone == :NonNeg) ? b[idx] :  Inf
        end
    end

    model.m,model.n = size(A)
    model.p = Presolve_Problem(false,model.m,model.n)
    model.A = A
    model.b = b
    model.c = c

    #TODO.. change the variable orders in presolve routines.

    make_presolve!(false,model.p,A,collb,colub,c,rowlb,rowub)

    #TODO.. add code for Dual Variables.
    newA,newcollb,newcolub,newc,newrowlb,newrowub = presolver!(false,model.p,A,collb,colub,c,rowlb,rowub)

    loadproblem!(model.innermodel,A,collb,colub,c,rowlb,rowub,model.orig_sense)

    return model
end

function optimize!(m::PresolveMathProgModel)
    constr_primalsol = Array{Float64,1}()
    constr_dualsol = Array{Float64,1}()
    var_primalsol = Array{Float64,1}()
    var_dualsol = Array{Float64,1}()

    if(length(find(m.p.independentvar))!=0)
        optimize!(m.innermodel)
        if(m.innermodel.solve_stat!= :Optimal)
            error("Problem status is not optimal, solution might be inaccurate")
        end
        m.solve_stat = m.innermodel.solve_stat
        constr_primalsol = m.innermodel.getconstrsolution()
        constr_dualsol = m.innermodel.getconstrduals()
        var_primalsol = m.innermodel.getsolution()
        var_dualsol = m.innermodel.getreducedcosts()

        #TODO.. need to account for DUAL VARIABLES HERE.
    end

    sol = return_postsolved(constr_primalsol,constr_dualsol,var_primalsol,var_dualsol,m.p.independentvar,m.p.pstack)

    m.primal_sol = sol[1]
    # WHAT exactly is the dual sol in the conic problem format.
    # figure that out by reading the cblib or mosek manual. Might need to convert lp dual to conic dual.
    #TODO
    m.dual_sol = sol[2]
    m.obj_val = dot(m.c, m.primal_sol)

    # you  should get the dual solution and then convert it to the type required for Conic Models.
end

#getsolution        for x ,             col primal or var_primalsol
#getconstrsolution  for y (Ax),         row primal or constr_primalsol
#getreducedcosts    for gamma-delta,    col dual or var_dualsol
#getconstrduals     for alpha-beta,     row dual or constr_dualsol

# TODO: getdual and getvardaul, can be done once add the dualsolution support for all the presolve routines.
function getdual(m::PresolveMathProgModel)
    dual = m.dual_sol[m.row_map_ind]
    # flip sign for NonPos since it's treated as NonNeg by SCS
    for i in 1:length(m.row_map_type)
        if m.row_map_type[i] == :NonPos
            dual[i] = -dual[i]
        end
    end
    return dual
end

function getvardual(m::PresolveMathProgModel)
    dual = zeros(length(m.col_map_ind))
    for i in 1:length(m.col_map_type)
        if m.col_map_type[i] == :Free
            continue # dual is zero
        elseif m.col_map_type[i] == :NonPos
            # flip sign for NonPos since it's treated as NonNeg by SCS
            dual[i] = -m.dual_sol[m.col_map_ind[i]]
        else
            dual[i] = m.dual_sol[m.col_map_ind[i]]
        end
    end
    return dual
end


text:
using GraphMatFun
dirlist=["exp" ; "sqrt_plus_one"]
for dir in dirlist
    data_dir=joinpath("..","data",dir);
    gen_code_dir=joinpath("..","data","generated",dir);
    # Generated code
    langlist=[LangJulia(),LangMatlab(), LangC_MKL(), LangC_OpenBLAS()];
    extlist=[".jl",".m","_MKL.c","_OpenBLAS.c"];
    for fname in readdir(data_dir);

        if (endswith(fname,".cgr"))
            # Only process cgr-files
            basename=split(fname,".")[1];
            println("Generating code for $fname: ");
            graph=import_compgraph(joinpath(data_dir,fname));
            compress_graph!(graph);

            if (eltype(graph) == BigFloat)
                graph=Compgraph(Float64,graph);
            end
            if (eltype(graph) == Complex{BigFloat})
                graph=Compgraph(ComplexF64,graph);
            end

            for (i,lang) in enumerate(langlist)
                ext=extlist[i];

                gen_code(joinpath(gen_code_dir,"$(basename)$(ext)"),
                         graph,
                         funname=basename,lang=lang);

            end

        end

    end
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using Test
using Profile
using gcSolver
using ForwardDiff
using LinearAlgebra

rxntfR = ones(gcSolver.Nparams) * 0.1

tps = [0.0, 0.1, 1.0, 10.0, 100.0, 1000.0, 10000.0, 100000.0]

include("model.jl")
include("varprop.jl")


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["earcut"], :earcut),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/SimonDanisch/EarCutBuilder/releases/download/v1.0.0"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/EarCut.v1.0.0.aarch64-linux-gnu.tar.gz", "5c9b3bbbefca40653a0f8c76b4fb81c1b3ffe38a595d4fea8ceb4042e61d076c"),
    Linux(:aarch64, :musl) => ("$bin_prefix/EarCut.v1.0.0.aarch64-linux-musl.tar.gz", "33804140cbc0548e5b7b6019dde699f00887874bb20d5455f72773c6c941e096"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/EarCut.v1.0.0.arm-linux-gnueabihf.tar.gz", "5530d6c7e745721df2eebb16b4434e6bc376f4f6e954cf684209bc10a7d80590"),
    Linux(:armv7l, :musl, :eabihf) => ("$bin_prefix/EarCut.v1.0.0.arm-linux-musleabihf.tar.gz", "6ad7fc447cc12a3c78520ee6f74263a7ddbc5a22efebad177e044478e9f1c009"),
    Linux(:i686, :glibc) => ("$bin_prefix/EarCut.v1.0.0.i686-linux-gnu.tar.gz", "af984da94d01608ed4e56500d98043622f0a6d8e9c89be7f871c9952dabda39e"),
    Linux(:i686, :musl) => ("$bin_prefix/EarCut.v1.0.0.i686-linux-musl.tar.gz", "fedf73674acc407fb619b838dd92624654c4e4af2589c7aa9d162544ca9dc0b4"),
    Windows(:i686) => ("$bin_prefix/EarCut.v1.0.0.i686-w64-mingw32.tar.gz", "d7e090ca7bff14ba1fcb0063445ebe29aaf67efbcefe05438169fab65c27b11e"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/EarCut.v1.0.0.powerpc64le-linux-gnu.tar.gz", "2fc25b6c4aedba11ebb6da708046d1d798e9dd3ba6a2cf7e522f08689f616de3"),
    MacOS(:x86_64) => ("$bin_prefix/EarCut.v1.0.0.x86_64-apple-darwin14.tar.gz", "fe39007cbc391b9021b37ea16cedb478f10ecae270f3950ea87fa394acee8bd0"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/EarCut.v1.0.0.x86_64-linux-gnu.tar.gz", "12624a769191b22ae1040f8171ae8d94b5ee8cf06288c58f216337d846d1258d"),
    Linux(:x86_64, :musl) => ("$bin_prefix/EarCut.v1.0.0.x86_64-linux-musl.tar.gz", "025fa1d07a067cb6b0eed1465f549ce88aac9732a39cded2f8a547eeb80a89ec"),
    FreeBSD(:x86_64) => ("$bin_prefix/EarCut.v1.0.0.x86_64-unknown-freebsd11.1.tar.gz", "6dea8e7ddea50a4e4426fd7dafe25807a6b27b3ebbbb3c821fc03e9e2f6e23e1"),
    Windows(:x86_64) => ("$bin_prefix/EarCut.v1.0.0.x86_64-w64-mingw32.tar.gz", "e0bc2089108bc2aabfb1e11d0a7ddf40adeec501fa8d1414871c0056b927762a"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform $(triplet(platform_key())) is not supported by this package!")
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products)


text:
## Kullback-Leibler Discrepancy as a function of the model order
# for white noise
using EstimateSDE
using Serialization
using Pipe: @pipe
using Plots
pgfplotsx()
using DataFrames
using Statistics
using StatsPlots
using REPL
using REPL.TerminalMenus

## Settings
name = "experiment01"
ssize = (300,200)

## Load batch settings and results
results_dir = settings_resultsdir(name)
S = deserialize(joinpath(results_dir,"settings.jld"))
df_eval = deserialize(joinpath(results_dir,"model-evaluation.jld")).df
df_data = deserialize(joinpath(results_dir,"data.jld")).df

df_eval = @pipe(
  innerjoin(df_eval,df_data,on=:index) |>
  rename(_,:label=>:kld_label,:T=>:kld_T))

## Generate graphs
figure_dir = joinpath(results_dir,"figures")
if !isdir(figure_dir) mkpath(figure_dir) end

@info "Figures are saved to $figure_dir"

## Count of selected model order
colfix = (engine=SdeRootsEngine,proclabel=:ICMLA_W,method="MLE pmax",kld_label="same")
# colfix = (engine=SdeRootsEngine,proclabel=:ICMLA_W,method="MLE pmax",kld_label="regular",kld_T=1.0)
df = @pipe(
  filter(df_eval) do x
    for (k,v) in pairs(colfix)
      if x[k]≠v return false end
    end
    return true
  end |>
  select(_,[:pmethod,:kld,:p])
)

## KLD as a function of model order
D = @pipe(
  groupby(df,:pmethod) |>
  combine(_) do df
    x = df.kld
    c = length(x)
    μ = mean(x); σ = std(x)/√c
    (kld=μ,kld_std=σ,count=c)
  end
)

## Theoretical result
pmax = maximum(D.pmethod)
ps = 0:pmax
kld_theory = 0.5*(ps.+1) # parameters: p SDE coefficients plus one for σϵ

## KLD vs order
plt = plot(xlabel="p",ylabel="Do",fg_legend=nothing,bg_legend=nothing,
  legend=:topleft,size=ssize,ylim=(0,5))
plot!(ps,kld_theory,linecolor=:gray,markercolor=:gray,markershape=:square,
  markerstrokecolor=:auto,markersize=2,label="theory")
plot!(D.pmethod,D.kld,yerror=D.kld_std,markerstrokecolor=:auto,label="SDE(p)",
  markershape=:square,markersize=2)
display(plt)

## Save figure
for ext in ["tex","pdf"]
  f = "whitenoise.$ext"
  savefig(plt,joinpath(figure_dir,f))
  @info "Saved $f"
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMI
using Flux
using DifferentialEquations: Tsit5

import Random
Random.seed!(1234);

FMUPath = joinpath(dirname(@__FILE__), "..", "model", "SpringPendulumExtForce1D.fmu")

t_start = 0.0
t_step = 0.01
t_stop = 5.0
tData = t_start:t_step:t_stop

# generate traing data
myFMU = fmiLoad(FMUPath)
fmiInstantiate!(myFMU; loggingOn=false)
fmiSetupExperiment(myFMU, t_start, t_stop)
fmiSetReal(myFMU, "mass_s0", 1.3)
fmiEnterInitializationMode(myFMU)
fmiExitInitializationMode(myFMU)
_, realSimData = fmiSimulateCS(myFMU, t_start, t_stop; recordValues=["mass.a"], setup=false, reset=false, saveat=tData)
fmiUnload(myFMU)

# setup traing data
function extForce(t)
    return [sin(t), cos(t)]
end
accData = collect(data[1] for data in realSimData.saveval)

# loss function for training
function losssum()
    solution = problem(extForce, t_step)

    accNet = collect(data[1] for data in solution)

    Flux.Losses.mse(accNet, accData)
end

# callback function for training
global iterCB = 0
global lastLoss = 0.0
function callb()
    global iterCB += 1
    global lastLoss

    if iterCB == 1
        lastLoss = losssum()
    end

    if iterCB % 50 == 0
        loss = losssum()
        @info "Loss: $loss"
        @test loss < lastLoss
        lastLoss = loss
    end
end

# Load FMUs
fmus = []
for i in 1:2 # how many FMUs do you want?
    fmu = fmiLoad(FMUPath)
    fmiInstantiate!(fmu; loggingOn=false)
    #fmiSetupExperiment(fmu, t_start, t_stop)
    #fmiEnterInitializationMode(fmu)
    #fmiExitInitializationMode(fmu)
    push!(fmus, fmu)
end

# NeuralFMU setup
total_fmu_outdim = sum(map(x->length(x.modelDescription.outputValueReferences), fmus))

net = Chain(
    Parallel(
        vcat,
        inputs -> fmi2InputDoStepCSOutput(fmus[1], t_step, inputs[1:1]),
        inputs -> fmi2InputDoStepCSOutput(fmus[2], t_step, inputs[2:2])
    ),
    Dense(total_fmu_outdim, 16, tanh),
    Dense(16, 16, tanh),
    Dense(16, length(fmus[1].modelDescription.outputValueReferences)),
)

problem = CS_NeuralFMU(fmus, net, (t_start, t_stop); saveat=tData)
@test problem != nothing

solutionBefore = problem(extForce, t_step)

# train it ...
p_net = Flux.params(problem)

optim = ADAM()
Flux.train!(losssum, p_net, Iterators.repeated((), 100), optim; cb=callb)

# check results
solutionAfter = problem(extForce, t_step)

fmiUnload(myFMU)


text:
function converge_ensemble!(e::IPM_Ensemble, instruction::Permute_Instruct; max_iterates=typemax(Int64), backup::Tuple{Bool,Integer}=(false,0), clean::Tuple{Bool,Integer,Integer}=(false,0,0), verbose::Bool=false, converge_criterion::String="standard", converge_factor::AbstractFloat=.001, progargs...)
    N = length(e.models); curr_it=length(e.log_Li)

    curr_it>1 && isfile(e.path*"/tuner") ? (tuner=deserialize(e.path*"/tuner")) : (tuner = Permute_Tuner(instruction)) #restore tuner from saved if any
    wk_mon = Worker_Monitor([1])
    meter = ProgressNS(e, wk_mon, tuner, 0.; start_it=curr_it, progargs...)

    converge_check = get_convfunc(converge_criterion)
    while !converge_check(e, converge_factor) && (curr_it <= max_iterates)
        warn,step_report = nested_step!(e, instruction) #step the ensemble
        warn == 1 && #"1" passed for warn code means no workers persist; all have hit the permute limit
                (@error "Failed to find new models, aborting at current iterate."; return e) #if there is a warning, iust return the ensemble and print info
        curr_it += 1

        tune_weights!(tuner, step_report)
        instruction = tuner.inst

        backup[1] && curr_it%backup[2] == 0 && e_backup(e,tuner) #every backup interval, serialise the ensemble and instruction
        clean[1] &&  !e.sample_posterior && curr_it%clean[2] == 0 && clean_ensemble_dir(e,clean[3]) #every clean interval, remove old discarded models

        update!(meter, converge_check(e,converge_factor,vals=true)...)
    end

    if converge_check(e,converge_factor)
        final_logZ = measure_evidence(e)
        @info "Job done, sampled to convergence. Final logZ $(final_logZ.val) ± $(final_logZ.err)"

        e_backup(e,tuner)
        clean[1] && !e.sample_posterior && clean_ensemble_dir(e,0) #final clean
        return final_logZ
    elseif curr_it==max_iterates
        @info "Job done, sampled to maximum iterate $max_iterates. Convergence criterion not obtained."

        e_backup(e,tuner)
        clean[1] && !e.sample_posterior && clean_ensemble_dir(e,0) #final clean
        return e.log_Zi[end]
    end
end

function converge_ensemble!(e::IPM_Ensemble, instruction::Permute_Instruct, wk_pool::Vector{Int64}; max_iterates=typemax(Int64), backup::Tuple{Bool,Integer}=(false,0), clean::Tuple{Bool,Integer,Integer}=(false,0,0), verbose::Bool=false, converge_criterion::String="standard", converge_factor::AbstractFloat=.001, progargs...)
    N = length(e.models)

    converge_check = get_convfunc(converge_criterion)
    model_chan= RemoteChannel(()->Channel{Tuple{Union{ICA_PWM_Model,String},Integer, AbstractVector{<:Tuple}}}(10*length(wk_pool))) #channel to take EM iterates off of
    job_chan = RemoteChannel(()->Channel{Tuple{<:AbstractVector{<:Model_Record}, Float64, Union{Permute_Instruct,String}}}(1))
    put!(job_chan,(e.models, e.contour, instruction))

    if !converge_check(e,converge_factor) #sequence workers only if not already converged
        @async sequence_workers(wk_pool, permute_IPM, e, job_chan, model_chan)
    end

    curr_it=length(e.log_Li)
    wk_mon=Worker_Monitor(wk_pool)
    curr_it>1 && isfile(e.path*"/tuner") ? (tuner=deserialize(e.path*"/tuner")) : (tuner = Permute_Tuner(instruction)) #restore tuner from saved if any
    meter = ProgressNS(e, wk_mon, tuner, 0.; start_it=curr_it, progargs...)

    while !converge_check(e, converge_factor) && (curr_it <= max_iterates)
        warn, step_report = nested_step!(e, model_chan, wk_mon) #step the ensemble
        warn == 1 && #"1" passed for warn code means no workers persist; all have hit the permute limit
                (@error "All workers failed to find new models, aborting at current iterate."; return e) #if there is a warning, iust return the ensemble and print info
        curr_it += 1
        tune_weights!(tuner, step_report)
        instruction = tuner.inst
        take!(job_chan); put!(job_chan,(e.models,e.contour,instruction))
        backup[1] && curr_it%backup[2] == 0 && e_backup(e,tuner) #every backup interval, serialise the ensemble and instruction
        clean[1] && !e.sample_posterior && curr_it%clean[2] == 0 && clean_ensemble_dir(e,clean[3]) #every clean interval, remove old discarded models

        update!(meter, converge_check(e,converge_factor,vals=true)...)
    end

    take!(job_chan); put!(job_chan, (e.models, e.contour, "stop")) #stop instruction terminates worker functions

    if converge_check(e,converge_factor)
        final_logZ = measure_evidence(e)
        @info "Job done, sampled to convergence. Final logZ $(final_logZ.val) ± $(final_logZ.err)"

        e_backup(e,tuner)
        clean[1] && !e.sample_posterior && clean_ensemble_dir(e,0) #final clean
        return final_logZ
    elseif curr_it==max_iterates
        @info "Job done, sampled to maximum iterate $max_iterates. Convergence criterion not obtained."

        e_backup(e,tuner)
        clean[1] && !e.sample_posterior && clean_ensemble_dir(e,0) #final clean
        return e.log_Zi[end]
    end
end

                function evidence_converge(e, evidence_fraction; vals=false)
                    val=lps(findmax([model.log_Li for model in e.models])[1],  e.log_Xi[end])
                    thresh=lps(log(evidence_fraction),e.log_Zi[end])
                    vals ? (return val, thresh) : (return val<thresh)
                end

                function compress_converge(e, compression_ratio; vals=false)
                    val=findmax([model.log_Li for model in e.models])[1]-e.contour
                    thresh=compression_ratio
                    vals ? (return val, thresh) : (return val<thresh)
                end

                function get_convfunc(criterion)
                    if criterion == "standard"
                        return evidence_converge
                    elseif criterion == "compression"
                        return compress_converge
                    else
                        throw(ArgumentError("Convergence criterion $criterion not supported! Try \"standard\" or \"compression\"."))
                    end
                end


text:
using Courgette
using Test

import Random
import GLPK
import LinearAlgebra
using JuMP
using SparseArrays: sprand

@testset "CGLP cuts fractional point" begin
    Random.seed!(33)
    A = round.(Int, 10 .* sprand(100, 50, 0.5))
    b = round.(Int, rand(100) .* 20)
    Ah = vcat(A, -Int.(LinearAlgebra.Diagonal(LinearAlgebra.I, 50)), Int.(LinearAlgebra.Diagonal(LinearAlgebra.I, 50)))
    bh = vcat(b, -ones(Int, 50), zeros(Int, 50))
    primal_model = Model(
        with_optimizer(GLPK.Optimizer, msg_lev = GLPK.OFF)
    )
    @variable(primal_model, 0 <= x[1:50] <= 1)
    @constraint(primal_model, A * x .>= b)
    @objective(primal_model, Min, sum(x))
    optimize!(primal_model)
    @test termination_status(primal_model) == MOI.OPTIMAL
    xvalue = JuMP.value.(x)
    jfrac = 1
    while jfrac <= length(xvalue)
        if xvalue[jfrac] > 5 * eps(Float64) && xvalue[jfrac] < 1 - 5 * eps(Float64)
            break
        end
        jfrac += 1
    end
    π = eachindex(xvalue) .== jfrac
    π0 = 0

    (cglp, variables) = Courgette.cut_generating_cone(with_optimizer(GLPK.Optimizer, msg_lev = GLPK.OFF), Ah, xvalue, bh, π, π0)
    Courgette.normalize_cglp!(Courgette.StandardNormalization(), cglp::JuMP.Model, variables)
    optimize!(cglp)
    @test termination_status(cglp) == MOI.OPTIMAL
    @test LinearAlgebra.dot(JuMP.value.(variables[:γ]), xvalue) - JuMP.value(variables[:γ0]) < 0
end


text:
@testset "DOAs" begin
    @testset "Static DOAs" begin
        doa = @inferred GNSSSimulator.propagate(
            SVector(1.0, 2.0, 3.0),
            1s,
            Random.GLOBAL_RNG
        )
        @test doa == SVector(1.0, 2.0, 3.0)
        @test get_doa(SVector(1.0, 2.0, 3.0)) == SVector(1.0, 2.0, 3.0)

        doa = @inferred GNSSSimulator.propagate(
            Spherical(1.0, 0.0, π / 2),
            1s,
            Random.GLOBAL_RNG
        )
        @test doa == Spherical(1.0, 0.0, π / 2)
        @test get_doa(Spherical(1.0, 0.0, π / 2)) ≈ SVector(0.0, 0.0, 1.0)
    end

    @testset "Dynamic DOAs" begin
        doas = [SVector(0.0, 1.0+0.1i, 2.0+0.1i) for i = 0:10]
        dynamic_doa = @inferred DynamicDOA(doas, 0.0s, 1.0Hz)
        next_dynamic_doa = @inferred GNSSSimulator.propagate(
            dynamic_doa,
            1s,
            Random.GLOBAL_RNG
        )
        @test next_dynamic_doa == DynamicDOA(doas, 1.0s, 1.0Hz)
        @test get_doa(next_dynamic_doa) == SVector(0.0, 1.1, 2.1)
        @test get_doa(
            GNSSSimulator.propagate(dynamic_doa, 100s, Random.GLOBAL_RNG)
        ) ≈ SVector(0.0, 2.0, 3.0)

        linear_dynamic_doa = @inferred LinearDynamicDOA(
            Spherical(1.0, 0.0, 0.0),
            0.1rad/s,
            0.0rad/s
        )
        next_linear_dynamic_doa = @inferred GNSSSimulator.propagate(
            linear_dynamic_doa,
            1s,
            Random.GLOBAL_RNG
        )
        @test get_doa(next_linear_dynamic_doa) ≈ Spherical(1.0, 0.1, 0.0)
    end
end


text:
# Benchmark iterative integration scheme against ground truth.
function benchmark_iterative(integration_scheme::AbstractIntegrationScheme, iterations::AbstractArray, noise_levels::AbstractArray, shape::AbstractSyntheticShape = SynthSphere(), r::Real = 50; log_scale::Bool = true)
    # setup ground truth
    Zₜ = ground_truth(shape, radius = r)
    Zₜ = Zₜ ./ maximum(Zₜ).*35

    # setup plotting
    scene = Scene()
    runs = first(size(iterations))

    # For each noise level bechmark the scheme
    for i in eachindex(noise_levels)

        # setup test
        p, q = synthetic_gradient(shape, radius = r, noise = noise_levels[i])
        RMSE = zeros(Float64, 0)
        scene2 = Scene()
        Zout = zeros(Float64, size(p))

        for j in eachindex(iterations)
            println("Part $i: $j/$runs")

            # Run tests
            Z = convert_gradient(integration_scheme, p, q, iterations[j])
            Z = Z./maximum(Z).*35
            Zout = Z

            # Calculate error
            E = 0.0
            for i in CartesianIndices(Z)
                E = E + (Zₜ[i] - Z[i])^2
            end
            push!(RMSE, sqrt(E / (first(size(Z)) * last(size(Z))))/(maximum(Zₜ) - minimum(Zₜ)))
        end

        # Change to log scale
        iterations_display = copy(iterations)
        if log_scale == true
            iterations_display = log10.(iterations)
        end

        # Display reconstructed surface
        range = 1.0:1:151
        scene2 = surface(range, range, Zout)

        # Plot results
        c = RGB{Float64}(rand(), rand(), rand())
        scatter!(iterations_display, RMSE, markersize = 0.25, color = c)
        lines!(scene, iterations_display, RMSE, color = c)
        display(Makie.vbox(scene, scene2))
        display(AbstractPlotting.PlotDisplay(), Makie.vbox(scene, scene2))
    end
end

# Benchmark non-iterative integration scheme against ground truth.
function benchmark_noniterative(integration_scheme::AbstractIntegrationScheme, noise_levels::AbstractArray, shape::AbstractSyntheticShape = SynthSphere(), r::Real = 50; trials::Int = 50)
    # setup ground truth
    Zₜ = ground_truth(shape, radius = r)
    Zₜ = Zₜ ./ maximum(Zₜ)

    # setup plotting
    RMSE = zeros(Float64, 0)
    error_bars = Array{Pair{Point{2,Float32},Point{2,Float32}}, 1}(undef, 0)
    scene = Scene()

    # For each noise level bechmark the scheme
    for i in eachindex(noise_levels)

        # Average across multiple trials
        error_list = Array{Float64}(undef, 0)
        for j = 1:trials
            println("Part $i: $j/$trials")

            # Setup and Run tests
            p, q = synthetic_gradient(shape, radius = r, noise = noise_levels[i])
            Z = convert_gradient(integration_scheme,p,q)
            Z = Z./maximum(Z)

            # Calculate error
            E = 0.0
            for i in CartesianIndices(Z)
                E = E + (Zₜ[i] - Z[i])^2
            end
            push!(error_list, sqrt(E / (first(size(Z)) * last(size(Z)))))
        end

        # Average error and calculate error bars.
        temp_error = mean(error_list)/(maximum(Zₜ) - minimum(Zₜ))
        push!(RMSE, temp_error)
        push!(error_bars, Point2f0(noise_levels[i], temp_error - std(error_list./(maximum(Zₜ) - minimum(Zₜ)))) => Point2f0(noise_levels[i], temp_error + std(error_list./(maximum(Zₜ) - minimum(Zₜ)))))
    end

    # Plot results
    scatter!(scene, noise_levels, RMSE, markersize = 0.15, color = :black)
    linesegments!(scene, error_bars, color = :red)
    display(scene)
end

# Benchmark multiple integration schemes (using defult settings) against ground truth for comparison.
function compare_benchmark(integration_schemes::Array{AbstractIntegrationScheme,1}, noise_levels::AbstractArray, shape = SynthSphere(), r::Real = 50; trials::Int = 50)
    # setup ground truth
    Zₜ = ground_truth(shape, radius = r)
    Zₜ = Zₜ./maximum(Zₜ)

    # setup plotting
    scene = Scene()
    runs = first(size(noise_levels))

    # Benchmark each scheme
    for i in eachindex(integration_schemes)

        # Setup error_list
        RMSE = zeros(Float64, 0)
        error_bars = Array{Pair{Point{2,Float32},Point{2,Float32}}, 1}(undef, 0)

        # For each noise level bechmark the scheme
        for j in eachindex(noise_levels)
            error_list = Array{Float64}(undef, 0)
            println("Part $i: $j/$runs")

            # Average across multiple trials
            for k = 1:trials
                # Setup and run test
                p, q = synthetic_gradient(shape, radius = r, noise = noise_levels[j])
                Z = convert_gradient(integration_schemes[i],p,q)
                Z = Z./maximum(Z)

                # Calculate error
                E = 0.0
                for i in CartesianIndices(Z)
                    E = E + (Zₜ[i] - Z[i])^2
                end
                push!(error_list, sqrt(E / (first(size(Z)) * last(size(Z)))))
            end

            # Average error and calculate error bars.
            temp_error = mean(error_list)/(maximum(Zₜ) - minimum(Zₜ))
            push!(RMSE, temp_error)
            push!(error_bars, Point2f0(noise_levels[j], temp_error - std(error_list./(maximum(Zₜ) - minimum(Zₜ)))) => Point2f0(noise_levels[j], temp_error + std(error_list./(maximum(Zₜ) - minimum(Zₜ)))))
        end

        #plot result
        c = RGB{Float64}(rand(),rand(),rand())
        scatter!(scene, noise_levels, RMSE, markersize = 0.15, color = c)
        linesegments!(scene, error_bars, color = c)
        display(scene)
    end
end

# Benchmark reconsrtucted shape against ground truth and return error map
function benchmark_shape(Z::AbstractArray, ground::AbstractArray)
    # Setup test
    z = copy(Z)
    g = copy(ground)
    z = z./maximum(z)
    g = g./maximum(g)
    a = zeros(size(z))

    # Calculate error
    E = 0.0
    for i in CartesianIndices(Z)
        E = E + (g[i] - z[i])^2
        a[i] = sqrt((g[i] - z[i])^2)
    end
    error = sqrt(E / (first(size(Z)) * last(size(Z))))
    return error/(maximum(g) - minimum(g)), a./(maximum(g) - minimum(g))
end


text:
# Some functions in this file are inspired from Polyhedra.jl
#
export set_rtol, set_ztol, set_atol

"""
    Tolerance{N<:Number}

Type that represents the tolerances for a given numeric type.

### Fields

- `rtol` -- relative tolerance
- `ztol` -- zero tolerance or absolute tolerance for comparison against zero
- `atol` -- absolute tolerance

### Notes

The type `Tolerance`, parametric in the numeric type `N`, is used to store default
values for numeric comparisons. It is mutable and setting the value of a field
affects the getter functions hence it can be used to fix the tolerance globally
in `LazySets`.

Default values are defined for the most commonly used numeric types, and for those
cases when other numeric types are needed one can extend the default values
as explained next.

The cases `Float64` and `Rational` are special in the sense that they are the most
commonly used types in applications. Getting and setting default tolerances
is achieved with the functions `_rtol` and `set_rtol` (and similarly for the other
tolerances); the implementation creates an instance of `Tolerance{Float64}`
(resp. `Tolerance{Rational}`) and sets some default values. Again since `Tolerance`
is mutable, setting a value is possible e.g. `set_rtol(Type{Float64}, ε)` for some
floating-point `ε`.

For all other cases, a dictionary mapping numeric types to instances of `Tolerance`
for that numeric type is used. For floating-point types, a default value has been
defined through `default_tolerance` as follows:

```julia
default_tolerance(N::Type{<:AbstractFloat}) = Tolerance(Base.rtoldefault(N), sqrt(eps(N)), zero(N))
```
Hence to set a single tolerance (either `rtol`, `ztol` or `atol`) for a given
floating-point type, use the corresponding `set_rtol` function, while the values
which have not been set will be pulled from `default_tolerance`. If you would like
to define the three default values at once, or are computing with a non floating-point
numeric type, you can just extend `default_tolerance(N::Type{<:Number})`.
"""
mutable struct Tolerance{N<:Number}
    rtol::N
    ztol::N
    atol::N
end

# global Float64 tolerances
const _TOL_F64 = Tolerance(Base.rtoldefault(Float64), Float64(10)*sqrt(eps(Float64)), zero(Float64))

_rtol(N::Type{Float64}) = _TOL_F64.rtol
_ztol(N::Type{Float64}) = _TOL_F64.ztol
_atol(N::Type{Float64}) = _TOL_F64.atol

set_rtol(N::Type{Float64}, ε::Float64) = _TOL_F64.rtol = ε
set_ztol(N::Type{Float64}, ε::Float64) = _TOL_F64.ztol = ε
set_atol(N::Type{Float64}, ε::Float64) = _TOL_F64.atol = ε

# global rational tolerances
const _TOL_RAT = Tolerance(zero(Rational), zero(Rational), zero(Rational))

_rtol(N::Type{<:Rational}) = _TOL_RAT.rtol
_ztol(N::Type{<:Rational}) = _TOL_RAT.ztol
_atol(N::Type{<:Rational}) = _TOL_RAT.atol

set_rtol(N::Type{<:Rational}, ε::Rational) = _TOL_RAT.rtol = ε
set_ztol(N::Type{<:Rational}, ε::Rational) = _TOL_RAT.ztol = ε
set_atol(N::Type{<:Rational}, ε::Rational) = _TOL_RAT.atol = ε

# global default tolerances for other numeric types
TOL_N = Dict{Type{<:Number}, Tolerance}()
_rtol(N::Type{<:Number}) = get!(TOL_N, N, default_tolerance(N)).rtol
_ztol(N::Type{<:Number}) = get!(TOL_N, N, default_tolerance(N)).ztol
_atol(N::Type{<:Number}) = get!(TOL_N, N, default_tolerance(N)).atol

set_rtol(N::Type{NT}, ε::NT) where {NT<:Number} = begin
    if N ∉ keys(TOL_N)
        TOL_N[N] = default_tolerance(N)
    end
    TOL_N[N].rtol = ε
end

set_ztol(N::Type{NT}, ε::NT) where {NT<:Number} = begin
    if N ∉ keys(TOL_N)
        TOL_N[N] = default_tolerance(N)
    end
    TOL_N[N].ztol = ε
end

set_atol(N::Type{NT}, ε::NT) where {NT<:Number} = begin
    if N ∉ keys(TOL_N)
        TOL_N[N] = default_tolerance(N)
    end
    TOL_N[N].atol = ε
end

default_tolerance(N::Type{<:Number}) = error("default tolerance for numeric type $N is not defined")
default_tolerance(N::Type{<:Rational}) = Tolerance(zero(N), zero(N), zero(N))
default_tolerance(N::Type{<:Integer}) = Tolerance(zero(N), zero(N), zero(N))
default_tolerance(N::Type{<:AbstractFloat}) = Tolerance(Base.rtoldefault(N), sqrt(eps(N)), zero(N))

"""
    _leq(x::N, y::N; [kwargs...]) where {N<:Real}

Determine if `x` is smaller than or equal to `y`.

### Input

- `x`      -- number
- `y`      -- another number (of the same numeric type as `x`)
- `kwargs` -- not used

### Output

A boolean that is `true` iff `x <= y`.

### Algorithm

This is a fallback implementation for numbers of type `Real`. If the arguments
are floating point numbers, see `_leq(x::AbstractFloat, y::AbstractFloat)`.
"""
_leq(x::N, y::N; kwargs...) where {N<:Real} = x <= y

"""
    _leq(x::N, y::M; [kwargs...]) where {N<:Real, M<:Real}

Determine if `x` is smaller than or equal to `y`.

### Input

- `x`      -- number
- `y`      -- another number (of possibly different numeric type than `x`)
- `kwargs` -- optional arguments; see `?_leq` for the available options

### Output

A boolean that is `true` iff `x <= y`.

### Algorithm

This implementation calls Julia's `promote(x, y)` function, which converts all
arguments to a common numeric type, returning them as a tuple. The conversion
is such that the common type to which the values are converted can represent
them as faithfully as possible.
"""
function _leq(x::N, y::M; kwargs...) where {N<:Real, M<:Real}
    return _leq(promote(x, y)...; kwargs...)
end

"""
    _geq(x::Real, y::Real; [kwargs...])

Determine if `x` is greater than or equal to `y`.

### Input

- `x` -- number
- `y` -- another number (of possibly different numeric type than `x`)

### Output

A boolean that is `true` iff `x >= y`.

### Algorithm

This function falls back to `_leq(y, x)`, with type promotion if needed. See the
documentation of `_leq` for further details.
"""
_geq(x::Real, y::Real; kwargs...) = _leq(y, x; kwargs...)

"""
    isapproxzero(x::N; ztol::Real=_ztol(N)) where {N<:Real}

Determine if `x` is approximately zero.

### Input

- `x`    -- number
- `ztol` -- (optional, default: `_ztol(N)`) tolerance against zero

### Output

A boolean that is `true` iff `x ≈ 0`.

### Algorithm

It is considered that `x ≈ 0` whenever `x` (in absolute value) is smaller than
the tolerance for zero, `ztol`.
"""
function isapproxzero(x::N; ztol::Real=_ztol(N)) where {N<:Real}
    return abs(x) <= ztol
end

"""
    _isapprox(x::N, y::N;
              rtol::Real=_rtol(N),
              ztol::Real=_ztol(N),
              atol::Real=_atol(N)) where {N<:Real}

Determine if `x` is approximately equal to `y`.

### Input

- `x`    -- number
- `y`    -- another number (of the same numeric type as `x`)
- `rtol` -- (optional, default: `_rtol(N)`) relative tolerance
- `ztol` -- (optional, default: `_ztol(N)`) absolute tolerance for comparison
            against zero
- `atol` -- (optional, default: `_atol(N)`) absolute tolerance

### Output

A boolean that is `true` iff `x ≈ y`.

### Algorithm

We first check if `x` and `y` are both approximately zero, using
`isapproxzero(x, y)`.
If that fails, we check if `x ≈ y`, using Julia's `isapprox(x, y)`.
In the latter check we use `atol` absolute tolerance and `rtol` relative
tolerance.

Comparing to zero with default tolerances is a special case in Julia's
`isapprox`, see the last paragraph in `?isapprox`. This function tries to
combine `isapprox` with its default values and a branch for `x ≈ y ≈ 0` which
includes `x == y == 0` but also admits a tolerance `ztol`.

Note that if `x = ztol` and `y = -ztol`, then `|x-y| = 2*ztol` and still
`_isapprox` returns `true`.
"""
function _isapprox(x::N, y::N;
                   rtol::Real=_rtol(N),
                   ztol::Real=_ztol(N),
                   atol::Real=_atol(N)) where {N<:Real}
    if isapproxzero(x, ztol=ztol) && isapproxzero(y, ztol=ztol)
        return true
    else
        return isapprox(x, y, rtol=rtol, atol=atol)
    end
end

# different numeric types with promotion
function _isapprox(x::N, y::M; kwargs...) where {N<:Real, M<:Real}
    return _isapprox(promote(x, y)...; kwargs...)
end

# generic "dense"
function _isapprox(x::AbstractVector{N}, y::AbstractVector{N};
                   rtol::Real=_rtol(N),
                   ztol::Real=_ztol(N),
                   atol::Real=_atol(N)) where {N<:Real}
    n = length(x)
    if length(x) != length(y)
        return false
    end
    @inbounds for i in 1:n
        if !_isapprox(x[i], y[i], rtol=rtol, ztol=ztol, atol=atol)
            return false
        end
    end
    return true
end

# sparse
function _isapprox(x::SparseVector{N}, y::SparseVector{N};
                   rtol::Real=_rtol(N),
                   ztol::Real=_ztol(N),
                   atol::Real=_atol(N)) where {N<:Real}
    @assert length(x) == length(y)
    return x.nzind == y.nzind && _isapprox(x.nzval, y.nzval, rtol=rtol, ztol=ztol, atol=atol)
end

# different numeric types with promotion
function _isapprox(x::AbstractVector{N}, y::AbstractVector{M};
                   kwargs...) where {N<:Real, M<:Real}
    _isapprox(promote(x, y)...; kwargs...)
end

# matrices
function _isapprox(A::AbstractMatrix{N}, B::AbstractMatrix{N};
                   rtol::Real=_rtol(N),
                   ztol::Real=_ztol(N),
                   atol::Real=_atol(N)) where {N<:Real}
    if size(A) != size(B)
        return false
    end
    @inbounds for i in eachindex(A)
        if !_isapprox(A[i], B[i], rtol=rtol, ztol=ztol, atol=atol)
            return false
        end
    end
    return true
end

# matrices with different numeric types with promotion
function _isapprox(A::AbstractMatrix{N}, B::AbstractMatrix{M};
                   kwargs...) where {N<:Real, M<:Real}
    _isapprox(promote(A, B)...; kwargs...)
end

"""
    ispermutation(u::AbstractVector{T}, v::AbstractVector) where {T}

Check that two vectors contain the same elements up to reordering.

### Input

- `u` -- first vector
- `v` -- second vector

### Output

`true` iff the vectors are identical up to reordering.

### Examples

```jldoctest
julia> LazySets.ispermutation([1, 2, 2], [2, 2, 1])
true

julia> LazySets.ispermutation([1, 2, 2], [1, 1, 2])
false
```

### Notes

Containment check is performed using `LazySets._in(e, v)`, so in the case of
floating point numbers, the precision to which the check is made is determined
by the type of elements in `v`. See `_in` and `_isapprox` for more information.
"""
function ispermutation(u::AbstractVector{T}, v::AbstractVector) where {T}
    if length(u) != length(v)
        return false
    end
    occurrence_map = Dict{T, Int}()
    has_duplicates = false
    for e in u
        if !_in(e, v)
            return false
        end
        if haskey(occurrence_map, e)
            occurrence_map[e] += 1
            has_duplicates = true
        else
            occurrence_map[e] = 1
        end
    end
    if has_duplicates
        for e in v
            if !haskey(occurrence_map, e) || occurrence_map[e] == 0
                return false
            end
            occurrence_map[e] -= 1
        end
    end
    return true
end

# alias for Julia's containment check
function _in(x, itr)
    return x ∈ itr
end

# approximate containment check for numbers in floating point
function _in(x::T, itr) where {T<:AbstractFloat}
    return any(y -> _isapprox(x, y), itr)
end

# approximate containment check for a vector in floating point
function _in(x::AbstractVector{T}, itr) where {T<:AbstractFloat}
    return any(y -> _isapprox(x, y), itr)
end

"""
    _leq(x::N, y::N;
         rtol::Real=_rtol(N),
         ztol::Real=_ztol(N),
         atol::Real=_atol(N)) where {N<:AbstractFloat}

Determine if `x` is smaller than or equal to `y`.

### Input

- `x`    -- number
- `y`    -- another number (of the same numeric type as `x`)
- `rtol` -- (optional, default: `_rtol(N)`) relative tolerance
- `ztol` -- (optional, default: `_ztol(N)`) absolute tolerance for comparison
            against zero
- `atol` -- (optional, default: `_atol(N)`) absolute tolerance

### Output

A boolean that is `true` iff `x <= y`.

### Algorithm

The `x <= y` comparison is split into `x < y` or `x ≈ y`; the latter is
implemented by extending Juila's built-in `isapprox(x, y)` with an absolute
tolerance that is used to compare against zero.
"""
function _leq(x::N, y::N;
              rtol::Real=_rtol(N),
              ztol::Real=_ztol(N),
              atol::Real=_atol(N)) where {N<:AbstractFloat}
    return x <= y || _isapprox(x, y, rtol=rtol, ztol=ztol, atol=atol)
end

"""
    minmax(A, B, C)

Compute the minimum and maximum of three numbers A, B, C.

### Input

- `A` -- first number
- `B` -- second number
- `C` -- third number

### Output

The minimum and maximum of three given numbers.

### Examples

```jldoctest
julia> LazySets.minmax(1.4, 52.4, -5.2)
(-5.2, 52.4)
```
"""
function minmax(A, B, C)
    if A > B
        min, max = B, A
    else
        min, max = A, B
    end
    if C > max
        max = C
    elseif C < min
        min = C
    end
    return min, max
end

"""
    arg_minmax(A, B, C)

Compute the index (1, 2, 3) of the minimum and maximum of three numbers A, B, C.

### Input

- `A` -- first number
- `B` -- second number
- `C` -- third number

### Output

The index of the minimum and maximum of the three given numbers.

### Examples

```jldoctest
julia> LazySets.arg_minmax(1.4, 52.4, -5.2)
(3, 2)
```
"""
function arg_minmax(A, B, C)
    if A > B
        min, max = B, A
        imin, imax = 2, 1
    else
        min, max = A, B
        imin, imax = 1, 2
    end
    if C > max
        imax = 3
    elseif C < min
        imin = 3
    end
    return imin, imax
end


text:
function get_node(obj::IAUPole, arg0::AbsoluteDate)
    return jcall(obj, "getNode", Vector3D, (AbsoluteDate,), arg0)
end

function get_node(obj::IAUPole, arg0::FieldAbsoluteDate)
    return jcall(obj, "getNode", FieldVector3D, (FieldAbsoluteDate,), arg0)
end

function get_pole(obj::IAUPole, arg0::AbsoluteDate)
    return jcall(obj, "getPole", Vector3D, (AbsoluteDate,), arg0)
end

function get_pole(obj::IAUPole, arg0::FieldAbsoluteDate)
    return jcall(obj, "getPole", FieldVector3D, (FieldAbsoluteDate,), arg0)
end

function get_prime_meridian_angle(obj::IAUPole, arg0::AbsoluteDate)
    return jcall(obj, "getPrimeMeridianAngle", jdouble, (AbsoluteDate,), arg0)
end

function get_prime_meridian_angle(obj::IAUPole, arg0::FieldAbsoluteDate)
    return jcall(obj, "getPrimeMeridianAngle", RealFieldElement, (FieldAbsoluteDate,), arg0)
end


text:
# This is file provides the ability for the user to override the default prior.

# This code is `include()`'ed inside of `/scripts/venus.jl`, so it will **literarly replace** the
# function lnprior inside `src/model.jl` (all of it!). Be careful to make sure it accepts the same arguments!

using DiskJockey.model
import DiskJockey.model.lnprior_base
import DiskJockey.constants.ModelException

# for example, a replacement for the standard model might look like
function lnprior(pars::ParametersStandard, grid::Grid)
    lnp = lnprior_base(pars)

    # Modify this code here to enforce your prior!

    r_out = grid.Rs[end]/AU # [AU]
    # A somewhat arbitrary cutoff regarding the gridsize to prevent the disk from being too large
    # to fit on the model grid.
    if (3 * pars.r_c) > r_out
        throw(ModelException("User prior violated."))
    end

    return lnp
end

# while an example for the cavity model would look like this
function lnprior(pars::ParametersCavity, grid::Grid)
    lnp = lnprior_base(pars)

    # Modify this code here to enforce your prior!

    r_in = grid.Rs[1]/AU # [AU]
    r_out = grid.Rs[end]/AU # [AU]
    # A somewhat arbitrary cutoff regarding the gridsize to prevent the disk from being too large
    # to fit on the model grid.

    # Also check to make sure that r_cav is less than r_c but larger than r_in.
    if (3 * pars.r_c) > r_out || pars.r_cav < r_in || pars.r_cav > pars.r_c
        throw(ModelException("User prior violated."))
    end

    return lnp

end

function lnprior(pars::ParametersVertical, grid::Grid)
    # Create a giant short-circuit or loop to test for sensical parameter values.
    if pars.M_star <= 0.0 || pars.ksi <= 0. || pars.T_10a <= 0. || pars.T_10m <= 0. || pars.r_c <= 0.0  || pars.T_10a > 1500. || pars.q_m < 0. || pars.q_a < 0. || pars.q_m > 1.0 || pars.q_a > 1.0 || pars.incl < 0. || pars.incl > 180. || pars.PA < -180. || pars.PA > 520. || pars.X_freeze > 1.0 || pars.sigma_s < 0.0
        throw(ModelException("Parameters outside of prior range."))
    end

    # Check to see that the temperatures make sense
    if pars.T_10m > pars.T_10a
        throw(ModelException("Atmosphere is cooler than midplane."))
    end

    # If we've passed all the hard-cut offs by this point, return the geometrical inclination prior.
    lnp = log(0.5 * sind(pars.incl))

    r_out = grid.Rs[end]/AU # [AU]
    # A somewhat arbitrary cutoff regarding the gridsize to prevent the disk from being too large
    # to fit on the model grid.
    if (3 * pars.r_c) > r_out
        throw(ModelException("Model radius too large for grid size."))
    else
        return lnp
    end

end


text:
##
ensureDir(dest) = mkpath(dirname(dest))
export ensureDir
##
pathExtract(path::AbstractString)::RegexMatch = match(r"(.*/)?(([^/]*?)(?:\.([^.]+))?)$", path)
export pathExtract

fileNameNoExt(path::AbstractString) = pathExtract(path)[3]
export fileNameNoExt
##


text:
export Pade21

"""
    Pade21()

This uses a Padé-approximation of type (2,1) for prediction.
"""
struct Pade21 <: AbstractPredictor end
struct Pade21Cache{T,AV<:AbstractVector{T}} <: AbstractPredictorCache
    x²::Vector{T}
    x³::Vector{T}
    x_h::AV

    u::Vector{T}
    u₁::Vector{T}
    u₂::Vector{T}
end

function cache(::Pade21, H, x, ẋ, t)
    x², x³, x_h = copy(ẋ), copy(ẋ), copy(x)
    u = evaluate(H, x, t)
    u₁, u₂ = copy(u), copy(u)
    Pade21Cache(x², x³, x_h, u, u₁, u₂)
end

@inline function g!(u, H, x, ẋ, t, h, x_h)
    @inbounds for i in eachindex(x)
        x_h[i] = x[i] + h * ẋ[i]
    end
    evaluate!(u, H, x_h, t + h)
end

@inline function g!(u, H, x, ẋ, ẍ, t, h, x_h)
    h² = 0.5 * h * h
    @inbounds for i in eachindex(x)
        x_h[i] = x[i] + h * ẋ[i] + h² * ẍ[i]
    end
    evaluate!(u, H, x_h, t + h)
end

function update!(cache::Pade21Cache, H, x, ẋ, t, Jac::JacobianMonitor, eval_err::Float64)
    # unpack stuff to make the rest easier to read
    @unpack u, u₁, u₂, x_h, x², x³ = cache
    ψ = 10eval_err
    h₂ = nthroot(ψ, 4)
    g!(u₁, H, x, ẋ, t, h₂, x_h)
    g!(u₂, H, x, ẋ, t, -h₂, x_h)

    h² = h₂^2
    @inbounds for i in eachindex(u₁)
        u[i] = -(u₁[i] + u₂[i]) / h²
    end
    LA.ldiv!(x², Jac, u)

    h₃ = nthroot(ψ, 5)
    g!(u₁, H, x, ẋ, x², t, h₃, x_h)
    g!(u₂, H, x, ẋ, x², t, -h₃, x_h)

    h³ = h₃^3
    @inbounds for i in eachindex(u₁)
        u[i] = -3 * (u₁[i] - u₂[i]) / h³
    end
    LA.ldiv!(x³, Jac, u)

    cache
end

function predict!(x̂, cache::Pade21Cache, H::HomotopyWithCache, x, t, Δt, ẋ, J)
    @unpack x², x³ = cache
    @inbounds for i in eachindex(x)
        δ = 1 - Δt * x³[i] / (3 * x²[i])
        if isnan(δ) || iszero(δ)
            x̂[i] = x[i] + Δt * ẋ[i] + 0.5 * Δt^2 * x²[i]
        else
            x̂[i] = x[i] + Δt * ẋ[i] + (0.5 * Δt^2 * x²[i]) / δ
        end
    end
    nothing
end

order(::Pade21Cache) = 4
@inline highest_derivative(cache::Pade21Cache) = (cache.x³, 3)
second_derivative(cache::Pade21Cache) = cache.x²
third_derivative(cache::Pade21Cache) = cache.x³


text:
# NOTE: As of https://github.com/JuliaLang/julia/pull/21746, all macro calls have
# a LineNumberNode. Since our functions that do expression transformation don't have
# access to the source code locations, we don't insert a LineNumberNode, which means
# that a transformed expression won't compare equal to an expression literal due to
# line information. Thus for testing purposes we define a function that removes line
# information from expressions as well as a custom equality operator that uses this.
# Note that it's only necessary for comparisons that deal with macro calls.
function skip_line_info(ex::Expr)
    map!(arg->arg isa LineNumberNode ? nothing : skip_line_info(arg), ex.args, ex.args)
    ex
end
skip_line_info(ex) = ex
≃(a, b) = skip_line_info(a) == skip_line_info(b)

@testset "code_transformation/differentiable" begin

    import Nabla.Nabla: unionise_type, unionise_arg, unionise_subtype, unionise_eval,
        unionise_macro_eval, unionise_sig, unionise, unionise_struct

    # Test Nabla.unionise_arg. This depends upon Nabla.unionise_type, so we express
    # our tests in terms of this.
    @test unionise_arg(:x) == :x
    @test unionise_arg(:foo) == :foo
    @test unionise_arg(:(::x)) == :(::$(unionise_type(:x)))
    @test unionise_arg(:(::T{V})) == :(::$(unionise_type(:(T{V}))))
    @test unionise_arg(:(x::T)) == :(x::$(unionise_type(:T)))
    @test unionise_arg(:(x::T{V})) == :(x::$(unionise_type(:(T{V}))))
    @test unionise_arg(:(::T{V} where V)) == :(::$(unionise_type(:(T{V} where V))))
    @test unionise_arg(:(::T{V} where V<:Q)) == :(::$(unionise_type(:(T{V} where V<:Q))))
    @test unionise_arg(:(x::T{V} where V<:Q)) == :(x::$(unionise_type(:(T{V} where V<:Q))))

    # Test unionise_arg for varargs case.
    let dots = Symbol("...")
        @test unionise_arg(Expr(dots, :x)) == Expr(dots, :x)
        @test unionise_arg(Expr(dots, :(x::T))) == Expr(dots, unionise_arg(:(x::T)))
    end

    # Test unionise_subtype.
    @test unionise_subtype(:T) == :T
    @test unionise_subtype(:(T<:V)) == :(T<:$(unionise_type(:V)))

    # Test Nabla.unionise_eval.
    @test unionise_eval(:(eval(:foo))) ≃ :(eval(:(@unionise foo)))
    @test unionise_eval(:(eval(DiffBase, :foo))) ≃ :(eval(DiffBase, :(@unionise foo)))
    @test unionise_eval(:(eval(:(println("foo"))))) ≃ :(eval(:(@unionise println("foo"))))
    @test unionise_eval(:(eval(DiffBase, :(println("foo"))))) ≃
        :(eval(DiffBase, :(@unionise println("foo"))))

    # Test Nabla.unionise_macro_eval.
    @test unionise_macro_eval(:(@eval foo)) ≃ :(@eval @unionise foo)
    @test unionise_macro_eval(:(@eval DiffBase foo)) ≃ :(@eval DiffBase @unionise foo)
    @test unionise_macro_eval(:(@eval println("foo"))) ≃ :(@eval @unionise println("foo"))
    @test unionise_macro_eval(:(@eval DiffBase println("foo"))) ≃
        :(@eval DiffBase @unionise println("foo"))

    # Test Nabla.unionise. This depends upon Nabla.unionise_arg, so we express
    # our tests in terms of this function.
    @test unionise_sig(:((x,))) == :((x,))
    @test unionise_sig(:((x, y))) == :((x, y))
    @test unionise_sig(:(foo(x))) == :(foo(x))
    @test unionise_sig(:(foo(x, y))) == :(foo(x, y))
    @test unionise_sig(:((x::T,))) == :(($(unionise_arg(:(x::T))),))
    @test unionise_sig(:(foo(x::T))) == :(foo($(unionise_arg(:(x::T)))))
    @test unionise_sig(:(foo(x::T) where T)) == :(foo($(unionise_arg(:(x::T)))) where T)

    # Test Nabla.unionise_struct. Written in terms of Nabla.unionise_arg.
    @test unionise_struct(:(struct Foo end)) == :(struct Foo end)
    @test unionise_struct(:(struct Foo{T} end)) ==
        :(struct Foo{$(unionise_subtype(:T))} end)
    @test unionise_struct(:(struct Foo{T<:Real} end)) ==
        :(struct Foo{$(unionise_subtype(:(T<:Real)))} end)
    @test unionise_struct(:(struct Foo{T<:V, Q<:U} end)) ==
        :(struct Foo{$(unionise_subtype(:(T<:V))), $(unionise_subtype(:(Q<:U)))} end)
    @test unionise_struct(:(struct Foo <: Bar end)) == :(struct Foo <: Bar end)
    @test unionise_struct(:(struct Foo{T} <: Bar end)) == :(struct Foo{T} <: Bar end)
    @test unionise_struct(:(struct Foo{T<:V} <: Bar end)) ==
        :(struct Foo{$(unionise_subtype(:(T<:V)))} <: Bar end)
    @test unionise_struct(:(struct Foo{T<:V, Q<:U} <: Bar end)) ==
        :(struct Foo{$(unionise_subtype(:(T<:V))), $(unionise_subtype(:(Q<:U)))} <: Bar end)

    # Test Nabla.make_accept_nodes. Heavily depends upon Nabla.unionise_sig, so we
    # express tests in terms of this function.
    @test unionise(:hi) == :hi
    @test unionise(:(N = 5)) == :(N = 5)
    @test unionise(:(N, M = 5, 4)) == :(N, M = 5, 4)
    @test unionise(:((N, M) = (5, 4))) == :((N, M) = (5, 4))
    @test unionise(:((N, M) = 5, 4)) == :((N, M) = 5, 4)
    @test unionise(:(N, M = (5, 4))) == :(N, M = (5, 4))
    @test unionise(Expr(Symbol("->"), :x, :x)) == Expr(Symbol("->"), :x, :x)
    @test unionise(Expr(Symbol("->"), :(x::Real), :(5x))) ==
        Expr(Symbol("->"), unionise_sig(:(x::Real)), :(5x))
    @test unionise(Expr(Symbol("->"), :((x::Real,)), :(5x))) ==
        Expr(Symbol("->"), unionise_sig(:((x::Real,))), :(5x))
    @test unionise(Expr(Symbol("->"), :((x::Real, y::Real)), :(5x))) ==
        Expr(Symbol("->"), unionise_sig(:((x::Real, y::Real))), :(5x))
    @test unionise(Expr(:function, :((x,)), :(return x))) ==
        Expr(:function, unionise_sig(:((x,))), :(return x))
    @test unionise(Expr(:function, :((x::Real,)), :(return x))) ==
        Expr(:function, unionise_sig(:((x::Real,))), :(return x))
    @test unionise(Expr(:function, :(foo(x)), :(return x))) ==
        Expr(:function, unionise_sig(:(foo(x))), :(return x))
    @test unionise(Expr(:function, :(foo(x::Real)), :(return x))) ==
        Expr(:function, unionise_sig(:(foo(x::Real))), :(return x))
    @test unionise(Expr(Symbol("="), :(foo(x)), :x)) ==
        Expr(Symbol("="), unionise_sig(:(foo(x))), :x)
    @test unionise(Expr(Symbol("="), :(foo(x::Real)), :x)) ==
        Expr(Symbol("="), unionise_sig(:(foo(x::Real))), :x)
    @test unionise(Expr(Symbol("="), :(foo(x::T, y::T) where T), :x)) ==
        Expr(Symbol("="), unionise_sig(:(foo(x::T, y::T) where T)), :x)
    @test unionise(:(eval(:foo))) == unionise_eval(:(eval(:foo)))
    @test unionise(:(eval(DiffBase, :foo))) == unionise_eval(:(eval(DiffBase, :foo)))
    @test unionise(:(@eval foo)) ≃ unionise_macro_eval(:(@eval foo))
    @test unionise(:(@eval DiffBase foo)) ≃ unionise_macro_eval(:(@eval DiffBase foo))
    @test unionise(:(struct Foo{T<:V} end)) == unionise_struct(:(struct Foo{T<:V} end))

    # @unionise with default values and keywords
    UT = unionise_type(:T)
    raw = unionise(:(f(x::T, y::T=2; z::T=4) = x + y + z))
    new = :(f(x::$UT, y::$UT=2; z::T=4) = x + y + z)
    @test raw ≃ new

    # @unionise error conditions
    @test_throws ArgumentError unionise(:(f(@nospecialize x) = x))
end


text:
using StatsBase
using Test

a = [1.0, 2.0, 2.0, 3.0, 4.0, 4.0, 4.0, 5.0]
x = [3.0, 1.0, 2.0, 4.0, 4.0, 2.0, 5.0, 4.0]  # x is a permutated version of a
xm = [3.0, 1.0, 2.0, 4.0, 4.0, 2.0, 5.0, 4.0, missing]
s = ["c", "a", "b", "d", "d", "b", "e", "d"] # s is a vector of strings ordered like x

@test ordinalrank(a) == [1, 2, 3, 4, 5, 6, 7, 8]
@test ordinalrank(x) == [4, 1, 2, 5, 6, 3, 8, 7]
@test isequal(ordinalrank(xm), [4, 1, 2, 5, 6, 3, 8, 7, missing])
@test isequal(ordinalrank([missing, missing]), [missing, missing])
@test ordinalrank(s) == ordinalrank(x)
@test ordinalrank(x, rev = true) == ordinalrank(-x)
@test ordinalrank(x, lt = (x, y) -> isless(y, x)) == ordinalrank(-x)

@test competerank(a) == [1, 2, 2, 4, 5, 5, 5, 8]
@test competerank(x) == [4, 1, 2, 5, 5, 2, 8, 5]
@test isequal(competerank(xm), [4, 1, 2, 5, 5, 2, 8, 5, missing])
@test isequal(competerank([missing, missing]), [missing, missing])
@test competerank(s) == competerank(x)
@test competerank(x, rev = true) == competerank(-x)
@test competerank(x, lt = (x, y) -> isless(y, x)) == competerank(-x)

@test denserank(a) == [1, 2, 2, 3, 4, 4, 4, 5]
@test denserank(x) == [3, 1, 2, 4, 4, 2, 5, 4]
@test isequal(denserank(xm), [3, 1, 2, 4, 4, 2, 5, 4, missing])
@test isequal(denserank([missing, missing]), [missing, missing])
@test denserank(s) == denserank(x)
@test denserank(x, rev = true) == denserank(-x)
@test denserank(x, lt = (x, y) -> isless(y, x)) == denserank(-x)

@test tiedrank(a) == [1.0, 2.5, 2.5, 4.0, 6.0, 6.0, 6.0, 8.0]
@test tiedrank(x) == [4.0, 1.0, 2.5, 6.0, 6.0, 2.5, 8.0, 6.0]
@test isequal(tiedrank(xm), [4.0, 1.0, 2.5, 6.0, 6.0, 2.5, 8.0, 6.0, missing])
@test isequal(tiedrank([missing, missing]), [missing, missing])
@test tiedrank(s) == tiedrank(x)
@test tiedrank(x, rev = true) == tiedrank(-x)
@test tiedrank(x, lt = (x, y) -> isless(y, x)) == tiedrank(-x)


text:
module TestUnits

using Modia
@usingModiaPlot

UnitTest = Model(

    T = 0.2u"s",
    u0 = 8.5u"kg",
    k = 1u"kg^-1",
    m = 10u"kg",
    F = 20u"N",

    x = Var(init = 5.5u"kg"),    # 5500.0u"g"),
    v = Var(init = 1u"m/s"),
    equations = :[
        T*der(x) + x = u0
        y = x
        z = x*u"kg^-1"

        F = m * a
        der(v) = a ]
)

model = @instantiateModel(UnitTest)

simulate!(model, Tsit5(), stopTime = 1.0, requiredFinalStates = [8.479786015016273, 2.999999999999999])

#plot(model, ["T", "x", "der(x)", "y", "a", "der(v)"])
plot(model, ["x", "der(x)", "y", "a", "der(v)"])

end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _geo2geodetic{T<:AbstractFloat}(lat::T, long::T, alt::T, eqrad::T, polrad::T)
    e = sqrt(eqrad^2 - polrad^2)/eqrad
    lat = deg2rad(lat)
    long_rad = deg2rad(long)
    x = (eqrad + alt)*cos(lat)*cos(long_rad)
    y = (eqrad + alt)*cos(lat)*sin(long_rad)
    z = (eqrad + alt)*sin(lat)
    r = hypot(x, y)
    s    = hypot(r, z)*(1.0 - eqrad*sqrt((1.0 - e^2)/((1.0 - e^2)*r^2 + z^2)))
    t0   = 1.0 + s*sqrt(1.0 - (e*z)^2/(r^2 + z^2))/eqrad
    dzeta1 = z*t0
    xi1  = r*(t0 - e^2)
    rho1 = hypot(xi1, dzeta1)
    c1   = xi1/rho1
    s1   = dzeta1/rho1
    b1   = eqrad/sqrt(1.0 - (e*s1)^2)
    u1   = b1*c1
    w1   = b1*s1*(1.0 - e^2)
    return rad2deg(atan2(s1, c1)), long, hypot(r - u1, z - w1)
end

"""
    geo2geodetic(latitude, longitude, altitude) -> latitude, longitude, altitude
    geo2geodetic(latitude, longitude, altitude, planet) -> latitude, longitude, altitude
    geo2geodetic(latitude, longitude, altitude, equatorial_radius, polar_radius) -> latitude, longitude, altitude

### Purpose ###

Convert from geographic (or planetographic) to geodetic coordinates.

### Explanation ###

Converts from geographic (latitude, longitude, altitude) to geodetic (latitude,
longitude, altitude).  In geographic coordinates, the Earth is assumed a perfect
sphere with a radius equal to its equatorial radius.  The geodetic (or
ellipsoidal) coordinate system takes into account the Earth's oblateness.

Geographic and geodetic longitudes are identical.  Geodetic latitude is the
angle between local zenith and the equatorial plane.  Geographic and geodetic
altitudes are both the closest distance between the satellite and the ground.

### Arguments ###

The function has two base methods.  The arguments common to all methods and
always mandatory are `latitude`, `longitude`, and `altitude`:

* `latitude`: geographic latitude, in degrees.
* `longitude`: geographic longitude, in degrees.
* `altitude`: geographic altitude, in kilometers.

In order to convert to geodetic coordinates, you can either provide custom
equatorial and polar radii of the planet or use the values of one of the planets
of Solar System (Pluto included).

If you want to use the method with explicit equatorial and polar radii the
additional mandatory arguments are:

* `equatorial_radius`: value of the equatorial radius of the body, in
  kilometers.
* `polar_radius`: value of the polar radius of the body, in kilometers.

Instead, if you want to use the method with the selection of a planet, the only
additional argument is the planet name:

* `planet` (optional string argument): string with the name of the Solar System
  planet, from "Mercury" to "Pluto".  If omitted (so, when only `latitude`,
  `longitude`, and `altitude` are provided), the default is "Earth".

In all cases, the three coordinates can be passed as a 3-tuple `(latitude,
longitude, altitude)`.  In addition, geographical `latitude`, `longitude`, and
`altitude` can be given as arrays of the same length.

### Output ###

The 3-tuple `(latitude, longitude, altitude)` in geodetic coordinates, for the
body with specified equatorial and polar radii (Earth by default).

If geographical coordinates are given as arrays, a 3-tuple of arrays of the same
length is returned.

### Method ###

Stephen P.  Keeler and Yves Nievergelt, "Computing geodetic coordinates", SIAM
Rev. Vol. 40, No. 2, pp. 300-309, June 1998
(DOI:[10.1137/S0036144597323921](http://dx.doi.org/10.1137/S0036144597323921)).

Planetary constants are from Planetary Fact Sheet
(http://nssdc.gsfc.nasa.gov/planetary/factsheet/index.html).

### Example ###

Locate the Earth geographic North pole (latitude: 90°, longitude: 0°, altitude 0
km), in geodetic coordinates:

``` julia
geo2geodetic(90, 0, 0)
# => (90.0,0.0,21.38499999999931)
```

The same for Jupiter:

``` julia
geo2geodetic(90, 0, 0, "Jupiter")
# => (90.0,0.0,4355.443799999994)
```

Find geodetic coordinates for point of geographic coordinates (latitude,
longitude, altitude) = (43.16°, -24.32°, 3.87 km) on a planet with equatorial
radius 8724.32 km and polar radius 8619.19 km:

``` julia
geo2geodetic(43.16, -24.32, 3.87, 8724.32, 8619.19)
# => (43.849399515234516,-24.32,53.53354478670836)
```

### Notes ###

Whereas the conversion from geodetic to geographic coordinates is given by an
exact, analytical formula, the conversion from geographic to geodetic isn't.
Approximative iterations (as used here) exist, but tend to become less good with
increasing eccentricity and altitude.  The formula used in this routine should
give correct results within six digits for all spatial locations, for an
ellipsoid (planet) with an eccentricity similar to or less than Earth's.  More
accurate results can be obtained via calculus, needing a non-determined amount
of iterations.

In any case, the function `geodetic2geo`, which converts from geodetic (or
planetodetic) to geographic coordinates, can be used to estimate the accuracy of
`geo2geodetic`.

``` julia
collect(geodetic2geo(geo2geodetic(67.2, 13.4, 1.2))) - [67.2, 13.4, 1.2]
# => 3-element Array{Float64,1}:
#     -3.56724e-9
#      0.0
#      9.47512e-10
```

Code of this function is based on IDL Astronomy User's Library.
"""
geo2geodetic(lat::Real, long::Real, alt::Real, eq::Real, pol::Real) =
    _geo2geodetic(promote(float(lat), float(long), float(alt),
                          float(eq), float(pol))...)

geo2geodetic(lla::Tuple{Real, Real, Real}, eq::Real, pol::Real) =
    geo2geodetic(lla..., eq, pol)

function geo2geodetic{LA<:Real, LO<:Real, AL<:Real}(lat::AbstractArray{LA},
                                                    long::AbstractArray{LO},
                                                    alt::AbstractArray{AL},
                                                    eq::Real, pol::Real)
    @assert length(lat) == length(long) == length(alt)
    typela  = typeof(float(one(LA)))
    outlat  = similar(lat, typela)
    outlong = similar(lat, typela)
    outalt  = similar(lat, typela)
    for i in eachindex(lat)
        outlat[i], outlong[i], outalt[i] =
            geo2geodetic(lat[i], long[i], alt[i], eq, pol)
    end
    return outlat, outlong, outalt
end

##### Select a planet.
geo2geodetic(lat::Real, long::Real, alt::Real, planet::AbstractString="earth") =
    (planet = lowercase(strip(planet));
     geo2geodetic(lat, long, alt, planets[planet].eqradius*1e-3,
                  planets[planet].polradius*1e-3))

geo2geodetic(lla::Tuple{Real, Real, Real}, planet::AbstractString="earth") =
    geo2geodetic(lla..., planet)

function geo2geodetic{LA<:Real, LO<:Real, AL<:Real}(lat::AbstractArray{LA},
                                                    long::AbstractArray{LO},
                                                    alt::AbstractArray{AL},
                                                    planet::AbstractString="earth")
    @assert length(lat) == length(long) == length(alt)
    typela  = typeof(float(one(LA)))
    outlat  = similar(lat, typela)
    outlong = similar(lat, typela)
    outalt  = similar(lat, typela)
    for i in eachindex(lat)
        outlat[i], outlong[i], outalt[i] =
            geo2geodetic(lat[i], long[i], alt[i], planet)
    end
    return outlat, outlong, outalt
end