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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ ae625e40-7f56-11ec-3cd5-d1244e23a3ba
begin
using PlutoUI, HTTP, Gumbo, Cascadia, CSV, JLD2

md"""# Fridge.jl
by *Ward Van Belle*
"""

end

# ╔═╡ 06db3cc2-1680-49b9-a054-281a6b35205a
md"""$(Resource("https://github.com/wardvanbelle/Fridge.jl/raw/master/FridgeLogo.png", :width => 150, :align => "right"))

A package that optimizes your fridge use while reducing your waste pile!!

This package tries to find the best recipes for you based on a recipe database. In our eyes (and the eyes of the objective function), the best recipes are the ones that use as much ingredients from your fridge as possible and that don't need extra ingredients from the grocery store."""

# ╔═╡ 3ad1af72-7096-4106-a3c7-47326695e2cb
md"""
## Checking The Ingredients
If you want to optimize fridge usage. It's best to know which ingredients are in the fridge. By checking all the ingredients and offering possible alternatives, it is easier for the algorithm to give better results. For instance, cheese may be replaced by swiss cheese.\
\
It takes two inputs:
- The fridge list: A list containing the different foods in your fridge as a string.
- The ingredient list: A list containing all the different ingredients that are used in the recipe database.

After searching for the different ingredients, it outputs an adapted fridge list.
"""

# ╔═╡ 8e58f807-1db9-495e-8679-430f953292ea
function checkIngredients(fridgeList,ingredientList)
    print("Checking if the food in your fridge is found in our database.\n\n")

    # check if the ingredients in the fridge are found in the database
    for food in fridgeList
        if food in ingredientList
            print("Found $food in the ingredient database.\n")
        else
            print("Did not find $food in the ingredient database.\n\n")
            if any(occursin.(food,ingredientList))
                print("Possible alternatives in database:\n")
                alternatives = ingredientList[occursin.(food,ingredientList)]
                for (indexNum, alternative) in enumerate(alternatives)
                    print("[$indexNum] $alternative\n")
                end
                print("If you want to take an alternative type its number, else type no.\n")
                answer = readline()
                if answer != "no"
                    correctInput = false
                    while !correctInput
                        try
                            answer = parse(Int64,answer)
                            correctInput = true
                        catch
                            print("Please only type the number of the alternative, eg. 1\n")
                            answer = readline()
                        end
                    end
                    fridgeList[fridgeList .== food] .= alternatives[answer]
                    print("Replaced $food with $(alternatives[answer])\n")
                end
            else
                print("Did not find any alternative in the ingredient database.\n")
            end
        end
    end
    print("Done checking ingredients.\n\n")

    return fridgeList
end

# ╔═╡ 264b1eb5-c421-4a69-a55c-8fbbbb7715d1
md"""
## Objective Functions
Of course, we want to be able to score how good our objective function will be. This objective function is based on the following formula where ``u`` is the number of ingredients used from the fridge, ``r`` is the number of ingredients remaining in the fridge and ``e`` are the extra ingredients needed. ``w_1``, ``w_2`` and ``w_3`` are a set of weights.

``score = w_1 * u + w_2 * r + w_3*e``
"""

# ╔═╡ 1752889a-513a-4620-add4-a9ca7fa92c8e
md"""
## The Search Algorithm
Let's now take a look at the different search algorithms that you can use in the Fridge.jl module.
"""

# ╔═╡ a65c0925-7a67-4da3-9a6d-2ade8d3529e9
md"""### Greedy Algorithm
This algorithm works like the one seen in class.
It ranks all recipes based on the previously discussed objective function. Picks the best scoring recipe. Adapts the recipe options and picks the new best recipe.
It keeps repeating this until the amount of chosen recipes is equal to `numRecipes` or if there are no recipes left to pick.\
\
As an input, it takes the earlier mentioned `fridgeList`, a dictionary containing all the recipes and their ingredients (`recipeDict`) and, the max number of recipes that a combination should have (`numRecipes`). \
\
After searching, it outputs the best-found combination.
"""

# ╔═╡ 954fe2b3-171b-43ad-a1df-cbb3414988de
md"""
### Simulated Annealing
To improve our earlier result, we can use simulated annealing. The algorithm below is based on the one seen in class. The biggest difference is the fact that we also use a Tabu list to block recipes for a certain number of cycles.

To start searching, it needs the earlier mentioned `fridgeList`, `recipeDict` and `numRecipes` but also the following things:
- `initSolution`: The initial solution from which it should start looking.
- `tabuList`: A list of recipes that should not be used in the found neighbour.
- `randRecipe`: A Boolean `true` or `false` value. When `true`, random recipes are used to create the neighbour.

After searching it also outputs the best-found combination of recipes.
"""

# ╔═╡ eee53924-fd20-44d8-a3eb-a2af9f9112b9
md"""
To optimize this searching method, there are a lot of parameters that can be tweaked.
The following options are available:
- `kT`: repetitions per temperature
- `r`: cooling rate
- `Tmax`: maximal temperature to start
- `Tmin`: minimal temperature to end
- `tabuLength`: Number of cycles that recipe needs to be blocked

You can experiment with them further on below!
"""

# ╔═╡ 7315f676-0d0e-43ae-8231-b6274dc5c8de
md"""
### Neighbours
There are two types of neighbours to use in the module. One uses a random combination of recipes, the other is based on the Greedy Search algorithm. Both are implemented in the `Neighbour()` function.
"""

# ╔═╡ 134a2dd7-69d2-449e-82ab-5bfff90ede1b
md"""A neighbour is created by first removing a recipe from the current combination `curSolution`. Next, the Neighbour function looks for other recipes that are compatible with the other recipes that are in the current combination."""

# ╔═╡ b6205bd3-41da-45e0-81f1-a480a3f0e737
md"""
## Overview Function
Finally, we want to bring all the functions together in one clean overview function.
The only things that need to be provided are:
- A list containing the different foods in your fridge as a string. `fridgeList`
- A relative or absolute path to a .csv or .jld2 file containing the recipe database. `dataPath`

Optionally you can also adapt the max amount of recipes in a combination (`numRecipes`) and wether or not random recipes should be used to create a neighbour (`randRecipe`).

"""

# ╔═╡ 6af2a6dd-ba7e-470e-8915-12ae6eff8357
md"""## Example Time
To start of we'll need to get our hands on a recipe database. Don't fear if you don't have one laying around. By using the `scrapeRecipe()` function provided in this module, you can create a .csv file containing recipes from [the cosylab recipe database](https://cosylab.iiitd.edu.in/recipedb/).

You just need to provide the following things:
- a starting recipe number
- a stopping recipe number
- the path to where the .csv file should be saved (you can copy this from your file explorer)

Next al that is left to do is to check the `Download database now` box. Since this will take a few minutes we provide you with a 🍪 and some ☕, enjoy!
"""

# ╔═╡ 656f211d-d1e8-4325-a12b-b89eca0351da
md"""
#### recipe database selection menu
The standard inputs are 106 541 to 106 866 since this downloads all the Belgian recipes. If you would like to download recipes from another country, please take a look at the explanation provided in the readme on the [GitHub repository](https://github.com/wardvanbelle/Fridge.jl).

Starting recipe number: $(@bind startNumber TextField(default="106541"))\
Stopping recipe number: $(@bind endNumber TextField(default="106866"))\
Filepath: $(@bind filePath TextField(default="../data/BelgianRecipeDB.csv"))\
Download database now: $(@bind downloadDB CheckBox())

It is recommended to only use this once. Otherwise, you will have copies of some recipes in your database. So you better uncheck the box above!
"""

# ╔═╡ bf2725c1-55e6-4c6a-9b83-53f5984cae75
md"Next it is time to look at the things we have left in our kitchen:"

# ╔═╡ fd02d76a-0ae6-4869-96bb-516c536fb91d
fridgeList = ["chocolate", "potato", "salt", "cheese"]

# ╔═╡ 34544e10-e4a1-4cfe-b0bc-43f8927bb19c
md"""
Now that we have a small database and the content of our fridge, we can start looking for a recipe combination. Keep an eye on your terminal since it is possible that you need to give some input! If you are ready, check the box below!
"""

# ╔═╡ 312fabf7-7d11-4bac-b52b-14de82d13e73
md"Ready to start searching? $(@bind readyToSearch CheckBox())"

# ╔═╡ 15a97461-7096-4ea0-acbc-edbc5db79cde
md"""
Congratulations! You made it through the notebook!
You can now continue to the **experimental section** below to play with the different parameters used in the algorithms. Enjoy the food! 🍔
"""

# ╔═╡ d49cf674-088f-47a1-9795-7641b3f01cd1
md"""## Experimental Section
We will now load the .csv file and save it in a dictionary so that it doesn't need to load every time you run the simulated annealing algorithm.
"""

# ╔═╡ 86f03d01-4f97-4487-9683-38f8dc9d826a
md"""
### Parameter Options:
**General**\
`numRecipes` : $(@bind numRecipesExp Slider(1:10, default=3, show_value=true))\
`tabuLength` : $(@bind tabuLengthExp Slider(1:15, default=5, show_value=true))\
`randRecipe` : $(@bind randRecipeExp CheckBox())\
\
**Simulated Annealing specific**\
`kT` : $(@bind kTExp Slider(20:10:150, default=100, show_value=true))\
`r` : $(@bind rExp Slider(0.1:0.05:0.95, default=0.75, show_value=true))\
`Tmax` : $(@bind TmaxExp Slider(2:10, default=4, show_value=true))\
`Tmin` : $(@bind TminExp Slider(1:9, default=1, show_value=true))\
\
**Objective specific**\
``w_1``: $(@bind w_1 Slider(0:10, default=1, show_value=true))\
``w_2``: $(@bind w_2 Slider(0:10, default=6, show_value=true))\
``w_3``: $(@bind w_3 Slider(0:10, default=2, show_value=true))\
"""

# ╔═╡ 70ffadd9-80d4-4eb9-92d9-e3d5ffa50074
md"Start experimenting: $(@bind startExperimenting CheckBox())"

# ╔═╡ c51f2b76-fe4c-4c29-87dc-c7784e5200c2
md"""
## Supporting Functions
Below you can find all the supporting functions used in this notebook. These can also be found in the src folder in the [GitHub repository](https://github.com/wardvanbelle/Fridge.jl). Please do not adapt them.
"""

# ╔═╡ 09d6393c-9b14-4401-babe-708a88ef0a16
function scrapeRecipe(scrapeBegin,scrapeEnd,csvPath)
"""
    scrapeRecipe(scrapeBegin,scrapeEnd,csvPath)

Download recipetitles and their corresponding ingredients from the recipe database of cosylab as a dictionary.
The recipes are downloaded based on their recipe number. This number can range from 2610 to 149191.
To get a recipenumber one should look at the last number of the url of a certain recipe.
For example, the recipe for 'Speculoosbavarois' is number 106585.
The recipes get automatically saved in a csv file where the first column is the recipetitle and the second column
is the list of ingredients.

## Input:
- scrapeBegin: The recipe number where the iteration should begin.
- scrapeEnd: The recipe number after which the iteration should and.
- csvPath: The path where the csv file is stored.

## Examples:

The example below downloads the recipes 2700 to 2702 and
stores them in the csv file 'recipedb.csv' in the current folder.

```julia-repl
julia> scrapeRecipe(2700,2702,"./recipedb.csv")
```
"""
    for i = scrapeBegin:scrapeEnd
        # get the webpage
        htmlText = HTTP.request("GET","https://cosylab.iiitd.edu.in/recipedb/search_recipeInfo/$i")
        htmlBody = parsehtml(String(htmlText.body))

        # get the needed elements out of the webpage
        recipeTitle = eachmatch(Selector("h3"),htmlBody.root)[1].children[1].text
        print("recipe number $i : $recipeTitle\n")
        ingredientTab = eachmatch(Selector("#ingredient_nutri"),htmlBody.root)
        ingredientTabel = eachmatch(Selector("#myTable"),ingredientTab[1])
        ingredientLinks = eachmatch(Selector("td a"),ingredientTabel[1])
        ingredientList = []
        for ingredientLink in ingredientLinks
            try
                append!(ingredientList, [ingredientLink.children[1].text])
            catch y
            end
        end

        # write the data to the DB
        if isfile(csvPath)
            CSV.write(csvPath, Dict(recipeTitle => ingredientList), append = true)
        else
            CSV.write(csvPath, Dict(recipeTitle => ingredientList), append = false)
        end

    end
    print("done")
end

# ╔═╡ cd5458a5-2136-4559-b4e5-2d91abea870c
if downloadDB
	scrapeRecipe(parse(Int64,startNumber), parse(Int64,endNumber), filePath)
end

# ╔═╡ 87b72838-260f-4ffb-ac79-03860165afa1
function loadRecipeDBCSV(csvPath)
    # read the dictionary from the csv file
    print("Loading the recipe database.\n")
    tempDict = CSV.File(csvPath) |> Dict
    recipeDict = Dict()

    # parse the ingredient list to the right format
    print("Parsing ingredients to right format.\n\n")
    for recipe in keys(tempDict)
        ingredientList = tempDict[recipe]
        recipeDict[recipe] = eval(Meta.parse(ingredientList))
    end
    return recipeDict
end

# ╔═╡ 11b589ff-3a09-4395-921a-872feb190022
if isfile(filePath)
	recipeDictExp = loadRecipeDBCSV(filePath)
end

# ╔═╡ adde7105-4a4a-453b-a7a8-9c9471dc29d3
function createIngredientDatabase(recipeDict)
    # create an unique vector of all ingredients used in the recipe database
    ingredientList = []
    for ingredients in values(recipeDict)
        append!(ingredientList,ingredients)
    end
    ingredientList = unique(ingredientList)

    return ingredientList
end

# ╔═╡ 89e58385-4db9-4ae3-9495-d8a2e63f2d72
function recipeToNumVector(fridgeList,ingredientList)
"""
    recipeToNumVector(fridgeList,ingredientList)

This function changes an ingredients list of a recipe to a numeric vector. This vector has a length equal
to the number of foods in your fridge plus one. If a certain food from the fridge is used in the recipe, that position in the vector is a 1.
If not it is a 0. The last index of the vector contains the amount of extra ingredients needed.

## Input:
- fridgeList: A list containing the different foods in your fridge as a string.
- ingredientList: A list containing all the ingredients used in a specific recipe.

## Output:
- numVector: This vector has a length equal to the number of foods in your fridge plus one. If a certain food from the fridge is used in the recipe,
that position in the vector is a 1. If not it is a 0. The last index of the vector contains the amount of extra ingredients needed.

"""
    numVector = zeros(Int64,length(fridgeList)+1)
    for i in 1:length(fridgeList)
        numVector[i] = fridgeList[i] in ingredientList ? 1 : 0
    end
    numVector[end] = length(ingredientList) - sum(numVector)
    return numVector
end

# ╔═╡ 3ed830c9-e3f7-4632-9ff3-cdd48c00f315
function randomCombo(fridgeList, recipeDict, numRecipes)

    randCombo = Dict()
    ingredientsArray = []
    namesArray = []

    for (name,ingredients) in recipeDict
        push!(namesArray,name)
        push!(ingredientsArray, recipeToNumVector(fridgeList, ingredients))
    end

    for i = 1:numRecipes

        randIndex = rand(1:length(namesArray))
        tempRecipeName = namesArray[randIndex]
        randCombo[tempRecipeName] = ingredientsArray[randIndex]

        if all(isone.(sum(values(randCombo))))
            break
        end

        namesArray = [name for (name, ingredientList) in zip(namesArray, ingredientsArray) if sum(ingredientList[1:end-1] .& randCombo[tempRecipeName][1:end-1] ) == 0]
        ingredientsArray = [ingredientList for ingredientList in ingredientsArray if sum(ingredientList[1:end-1] .& randCombo[tempRecipeName][1:end-1] ) == 0]

        if isempty(namesArray)
            break
        end

    end

    return randCombo
end

# ╔═╡ a6e77624-606a-49ea-941d-8180cbd94a27
compatible(x) = !any(sum(x)[1:end-1] .>= 2) # checks if two recipes use a same ingredient

# ╔═╡ 644c2809-5b2b-4426-8deb-f952372a8438
begin

# objective function used to score a single recipe
fridgeObjective(x::Array{Int64}) = w_1 * sum(x[1:end-1]) + w_2 *sum(x[1:end-1] .== 0) + w_3 *x[end]

# objective function used to score a combination of recipes
fridgeObjective(x) = compatible(x) ? sum(sum(x) .== 1)*w_1 + sum(sum(x) .== 0)*w_2 + sum(x)[end]*w_3 : Inf

end

# ╔═╡ ace6ecfa-319b-48f3-a675-51b4bc26abd0
function greedyFindCombo(fridgeList, recipeDict, numRecipes)

    bestCombo = Dict()
    ingredientsArray = []
    namesArray = []

    for (name,ingredients) in recipeDict
        push!(namesArray,name)
        push!(ingredientsArray, recipeToNumVector(fridgeList, ingredients))
    end

    for i = 1:numRecipes

        bestOrder = sortperm(ingredientsArray, by=fridgeObjective)
        ingredientsArray = ingredientsArray[bestOrder]
        namesArray = namesArray[bestOrder]

        # break if all recipes that are left, don't use anything from the fridge
        if all(ingredientsArray[1][1:end-1] .== 0)
            break
        end

        tempRecipeName = namesArray[1]
        bestCombo[tempRecipeName] = ingredientsArray[1]

        # break if all ingredients are used
        if all(isone.(sum(values(bestCombo))))
            break
        end

        namesArray = [name for (name, ingredientList) in zip(namesArray, ingredientsArray) if sum(ingredientList[1:end-1] .& bestCombo[tempRecipeName][1:end-1] ) == 0]
        ingredientsArray = [ingredientList for ingredientList in ingredientsArray if sum(ingredientList[1:end-1] .& bestCombo[tempRecipeName][1:end-1] ) == 0]

        # break if there are no recipes left
        if isempty(namesArray)
            break
        end

    end

    return bestCombo
end

# ╔═╡ 7365d57b-fd33-46af-9dcd-2c07908f64a3
function Neighbour(curSolution, fridgeList, recipeDict, numRecipes, tabuList, randRecipe)

    toRemove = rand(curSolution)[1]
    # adapt the fridgeList so that only ingredients from the removed ingredient are available
    tempFridgeList = copy(fridgeList)
    for recipe in keys(curSolution)
        if recipe != toRemove
            tempFridgeList = [i for i in fridgeList if !in(i,recipeDict[recipe])]
            numRecipes -= 1
        end
    end

    # adapt the recipeDict and use greedy search to find a new solution
    tempRecipeDict = copy(recipeDict)
    for recipe in keys(curSolution)
        delete!(tempRecipeDict,recipe)
    end

    for recipe in tabuList
        try delete!(tempRecipeDict,recipe)
        catch e
        end
    end

    if randRecipe
        neighbour = randomCombo(tempFridgeList, tempRecipeDict, numRecipes)
    else
        neighbour = greedyFindCombo(tempFridgeList, tempRecipeDict, numRecipes)
    end

    # correct recipe vectors
    for recipe in keys(neighbour)
        neighbour[recipe] = recipeToNumVector(fridgeList,recipeDict[recipe])
    end

    # here combine the two dictionaries
    tempCurSolution = copy(curSolution)
    delete!(tempCurSolution,toRemove)
    neighbour = merge(neighbour,tempCurSolution)

    return neighbour
end

# ╔═╡ 6f1d6e44-30a1-42ee-b888-4cbb105b5348
function SAFindCombo(initSolution,  fridgeList, recipeDict, numRecipes, randRecipe;
    kT=100, # repetitions per temperature
    r=0.75, # cooling rate
    Tmax=4, # maximal temperature to start
    Tmin=1, # minimal temperature to end
    tabuLength=3) # number of cycli that recipe needs to be blocked

    @assert 0 < Tmin < Tmax "Temperatures should be positive"
	@assert 0 < r < 1 "cooling rate is between 0 and 1"
	solution = initSolution
	obj = fridgeObjective([i for i in values(solution)])
    tabuList = String[i for i in keys(initSolution)]

	# current temperature
	T = Tmax
	while T > Tmin
		# repeat kT times
		for i in 1:kT
			sn = Neighbour(solution, fridgeList, recipeDict, numRecipes, tabuList, randRecipe)  # random neighbor
			obj_sn = fridgeObjective([i for i in values(sn)])
			# if the neighbor improves the solution, keep it
			# otherwise accept with a probability determined by the
			# Metropolis heuristic
			if obj_sn < obj || rand() < exp(-(obj_sn-obj)/T)
				solution = sn
				obj = obj_sn
			end

            for recipe in keys(sn)
                if !in(recipe, tabuList)
                    if length(tabuList) < tabuLength
                        pushfirst!(tabuList,recipe)
                    else
                        pop!(tabuList)
                        pushfirst!(tabuList,recipe)
                    end
                end
            end
		end

		# decay temperature
		T *= r
	end

    return solution
end

# ╔═╡ 9bfc9f3d-bc05-4c15-8b4c-2cfce3507afc
function findBestRecipe(fridgeList, dataPath; numRecipes=3, randRecipe=false)

    # load the recipe dictionary from the db file
    recipeDict = dataPath[end-3:end] == ".csv" ? loadRecipeDBCSV(dataPath) : load(dataPath)

    # create a list of all ingredients in your database
    ingredientList = createIngredientDatabase(recipeDict)

    # check for every food in your fridge if it's in the database. If not check if there are alternatives.
    fridgeList = checkIngredients(fridgeList, ingredientList)

    # find the best greedy recipe
    greedySolution = greedyFindCombo(fridgeList, recipeDict, numRecipes)

    # find the best recipe with SA
    SASolution = SAFindCombo(greedySolution,  fridgeList, recipeDict, numRecipes, randRecipe)

    # print the solution
    print("The best recipes to make are:\n\n")
    for recipeName in keys(SASolution)
        print("$(recipeName) : $(recipeDict[recipeName])\n")
    end

    return SASolution
end

# ╔═╡ 3569550e-0157-4d7a-8cda-4df327c92b81
if readyToSearch
	findBestRecipe(fridgeList, filePath)
end

# ╔═╡ e10c7ddf-aef0-470a-be22-ac250d84fb7e
begin
if startExperimenting
	greedySolutionExp = greedyFindCombo(fridgeList, recipeDictExp, numRecipesExp)
	SASolutionExp = SAFindCombo(greedySolutionExp,  fridgeList, recipeDictExp, numRecipesExp, randRecipeExp, kT=kTExp, r=rExp, Tmax=TmaxExp, Tmin=TminExp,  tabuLength=tabuLengthExp)

	md"""
	greedy output: $([recipe * ", " for recipe in keys(greedySolutionExp)])\
	SA output: $([recipe * ", " for recipe in keys(SASolutionExp)])
	"""
else
	md"Check the `Start experimenting` checkbox to see the output of your parameter combination."
end
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
Cascadia = "54eefc05-d75b-58de-a785-1a3403f0919f"
Gumbo = "708ec375-b3d6-5a57-a7ce-8257bf98657a"
HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
CSV = "~0.10.2"
Cascadia = "~1.0.1"
Gumbo = "~0.8.0"
HTTP = "~0.9.17"
JLD2 = "~0.4.18"
PlutoUI = "~0.7.32"
"""

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

[[AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

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

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

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

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

[[CSV]]
deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings"]
git-tree-sha1 = "9519274b50500b8029973d241d32cfbf0b127d97"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.10.2"

[[Cascadia]]
deps = ["AbstractTrees", "Gumbo"]
git-tree-sha1 = "95629728197821d21a41778d0e0a49bc2d58ab9b"
uuid = "54eefc05-d75b-58de-a785-1a3403f0919f"
version = "1.0.1"

[[CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

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

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.41.0"

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

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

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

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

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

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

[[FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "67551df041955cc6ee2ed098718c8fcd7fc7aebe"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.12.0"

[[FilePathsBase]]
deps = ["Compat", "Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "04d13bfa8ef11720c24e4d840c0033d145537df7"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.17"

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

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[Gumbo]]
deps = ["AbstractTrees", "Gumbo_jll", "Libdl"]
git-tree-sha1 = "e711d08d896018037d6ff0ad4ebe675ca67119d4"
uuid = "708ec375-b3d6-5a57-a7ce-8257bf98657a"
version = "0.8.0"

[[Gumbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "29070dee9df18d9565276d68a596854b1764aa38"
uuid = "528830af-5a63-567c-a44a-034ed33b8444"
version = "0.10.2+0"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "0fa77022fe4b511826b39c894c90daf5fce3334a"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.17"

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

[[HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

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

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[InlineStrings]]
deps = ["Parsers"]
git-tree-sha1 = "8d70835a3759cdd75881426fced1508bb7b7e1b6"
uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
version = "1.1.1"

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

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

[[JLD2]]
deps = ["DataStructures", "FileIO", "MacroTools", "Mmap", "Pkg", "Printf", "Reexport", "TranscodingStreams", "UUIDs"]
git-tree-sha1 = "39f22411266cdd1621092c762a3f0648dbdc8433"
uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
version = "0.4.18"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "22df5b96feef82434b07327e2d3c770a9b21e023"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.0"

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

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

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

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

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

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

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

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

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

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

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

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

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

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

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

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

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

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

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

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "db3a23166af8aebf4db5ef87ac5b00d36eb771e2"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.0"

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

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

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

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

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

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

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

[[SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "15dfe6b103c2a993be24404124b8791a09460983"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.11"

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

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

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

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

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

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

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

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

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

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

[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

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

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

[[WeakRefStrings]]
deps = ["DataAPI", "InlineStrings", "Parsers"]
git-tree-sha1 = "c69f9da3ff2f4f02e811c3323c22e5dfcb584cfa"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.4.1"

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

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

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

# ╔═╡ Cell order:
# ╟─ae625e40-7f56-11ec-3cd5-d1244e23a3ba
# ╟─06db3cc2-1680-49b9-a054-281a6b35205a
# ╟─3ad1af72-7096-4106-a3c7-47326695e2cb
# ╠═8e58f807-1db9-495e-8679-430f953292ea
# ╟─264b1eb5-c421-4a69-a55c-8fbbbb7715d1
# ╠═644c2809-5b2b-4426-8deb-f952372a8438
# ╟─1752889a-513a-4620-add4-a9ca7fa92c8e
# ╟─a65c0925-7a67-4da3-9a6d-2ade8d3529e9
# ╠═ace6ecfa-319b-48f3-a675-51b4bc26abd0
# ╟─954fe2b3-171b-43ad-a1df-cbb3414988de
# ╠═6f1d6e44-30a1-42ee-b888-4cbb105b5348
# ╟─eee53924-fd20-44d8-a3eb-a2af9f9112b9
# ╟─7315f676-0d0e-43ae-8231-b6274dc5c8de
# ╠═3ed830c9-e3f7-4632-9ff3-cdd48c00f315
# ╟─134a2dd7-69d2-449e-82ab-5bfff90ede1b
# ╠═7365d57b-fd33-46af-9dcd-2c07908f64a3
# ╟─b6205bd3-41da-45e0-81f1-a480a3f0e737
# ╠═9bfc9f3d-bc05-4c15-8b4c-2cfce3507afc
# ╟─6af2a6dd-ba7e-470e-8915-12ae6eff8357
# ╟─656f211d-d1e8-4325-a12b-b89eca0351da
# ╠═cd5458a5-2136-4559-b4e5-2d91abea870c
# ╟─bf2725c1-55e6-4c6a-9b83-53f5984cae75
# ╠═fd02d76a-0ae6-4869-96bb-516c536fb91d
# ╟─34544e10-e4a1-4cfe-b0bc-43f8927bb19c
# ╟─312fabf7-7d11-4bac-b52b-14de82d13e73
# ╠═3569550e-0157-4d7a-8cda-4df327c92b81
# ╟─15a97461-7096-4ea0-acbc-edbc5db79cde
# ╟─d49cf674-088f-47a1-9795-7641b3f01cd1
# ╠═11b589ff-3a09-4395-921a-872feb190022
# ╟─86f03d01-4f97-4487-9683-38f8dc9d826a
# ╟─70ffadd9-80d4-4eb9-92d9-e3d5ffa50074
# ╟─e10c7ddf-aef0-470a-be22-ac250d84fb7e
# ╟─c51f2b76-fe4c-4c29-87dc-c7784e5200c2
# ╟─09d6393c-9b14-4401-babe-708a88ef0a16
# ╟─87b72838-260f-4ffb-ac79-03860165afa1
# ╟─adde7105-4a4a-453b-a7a8-9c9471dc29d3
# ╟─89e58385-4db9-4ae3-9495-d8a2e63f2d72
# ╟─a6e77624-606a-49ea-941d-8180cbd94a27
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using Gadfly

set_default_plot_size(12inch, 6inch)

independent = [1,2,3,1,2,3]
dependent = [4,5,6,1,2,3]
group = [1,1,1,2,2,2]
vpos = plot(x=independent, y=  dependent, color=group, Geom.bar(position=:stack, orientation=:vertical));
hpos = plot(x=  dependent, y=independent, color=group, Geom.bar(position=:stack, orientation=:horizontal));
vneg = plot(x=independent, y= -dependent, color=group, Geom.bar(position=:stack, orientation=:vertical));
hneg = plot(x= -dependent, y=independent, color=group, Geom.bar(position=:stack, orientation=:horizontal));
gridstack([vpos hpos; vneg hneg])

ext: jl


text:
using Plots
include("src/FluidQueues.jl")
include("model_def.jl")
h(t) = (t>0.0) ? FluidQueues.hitting_times_cdf(model,t,1.0,2.0) : zeros(2,4)
tvec = 0.0:0.1:3
cdfs = zeros(length(tvec),4)
for (c,t) in enumerate(tvec)
    h_mat = h(t)
    cdfs[c,:] = h_mat[3:6]
end
plot(tvec,cdfs)


text:
"""
    synthesis(f0::Vector{Float64}, spectrogram::Matrix{Float64}, aperiodicity::Matrix{Float64}, fs::Int, frame_period::Float64 = DEFAULT_FRAME_PERIOD)

Synthesize an audio signal `y::Vector{Float64}` using WORLD from and `f0` track, a spectrogram, and an aperiodicity map, at the sampling rate `fs` per second.
Optionally specify a `frame_period` in milliseconds.
"""
function synthesis(f0::Vector{Float64}, spectrogram::Matrix{Float64}, aperiodicity::Matrix{Float64}, fs::Real, frame_period::Float64 = DEFAULT_FRAME_PERIOD)

    fs = convert(Int, fs)

    if size(spectrogram) != size(aperiodicity)
        error("Sizes of spectrogram $(size(spectrogram)) and aperiodicity $(size(aperiodicity)) don't match.")
    end

    nfreqs = size(spectrogram, 1)

    aperiodicity_pointer_vector = [Ref(aperiodicity, offset)
        for offset in 1:nfreqs:length(aperiodicity)]

    spectrogram_pointer_vector = [Ref(spectrogram, offset)
        for offset in 1:nfreqs:length(spectrogram)]

    fft_size = (nfreqs - 1) * 2

    f0_length = length(f0)
    y_length = trunc(Int, f0_length * frame_period * fs / 1000)
    y = Vector{Float64}(undef, y_length)

    @ccall libworld.Synthesis(
        f0::Ptr{Cdouble},
        f0_length::Cint,
        spectrogram_pointer_vector::Ptr{Ptr{Cdouble}},
        aperiodicity_pointer_vector::Ptr{Ptr{Cdouble}},
        fft_size::Cint,
        frame_period::Cdouble,
        fs::Cint,
        y_length::Cint,
        y::Ptr{Cdouble}
    )::Cvoid

    y
end


text:
function GHmsjPolynomials(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jint)
    return GHmsjPolynomials((jdouble, jdouble, jdouble, jdouble, jint), arg0, arg1, arg2, arg3, arg4)
end

function get_gmsj(obj::GHmsjPolynomials, arg0::jint, arg1::jint, arg2::jint)
    return jcall(obj, "getGmsj", jdouble, (jint, jint, jint), arg0, arg1, arg2)
end

function get_hmsj(obj::GHmsjPolynomials, arg0::jint, arg1::jint, arg2::jint)
    return jcall(obj, "getHmsj", jdouble, (jint, jint, jint), arg0, arg1, arg2)
end

function getd_gmsd_alpha(obj::GHmsjPolynomials, arg0::jint, arg1::jint, arg2::jint)
    return jcall(obj, "getdGmsdAlpha", jdouble, (jint, jint, jint), arg0, arg1, arg2)
end

function getd_gmsd_beta(obj::GHmsjPolynomials, arg0::jint, arg1::jint, arg2::jint)
    return jcall(obj, "getdGmsdBeta", jdouble, (jint, jint, jint), arg0, arg1, arg2)
end

function getd_gmsdh(obj::GHmsjPolynomials, arg0::jint, arg1::jint, arg2::jint)
    return jcall(obj, "getdGmsdh", jdouble, (jint, jint, jint), arg0, arg1, arg2)
end

function getd_gmsdk(obj::GHmsjPolynomials, arg0::jint, arg1::jint, arg2::jint)
    return jcall(obj, "getdGmsdk", jdouble, (jint, jint, jint), arg0, arg1, arg2)
end

function getd_hmsd_alpha(obj::GHmsjPolynomials, arg0::jint, arg1::jint, arg2::jint)
    return jcall(obj, "getdHmsdAlpha", jdouble, (jint, jint, jint), arg0, arg1, arg2)
end

function getd_hmsd_beta(obj::GHmsjPolynomials, arg0::jint, arg1::jint, arg2::jint)
    return jcall(obj, "getdHmsdBeta", jdouble, (jint, jint, jint), arg0, arg1, arg2)
end

function getd_hmsdh(obj::GHmsjPolynomials, arg0::jint, arg1::jint, arg2::jint)
    return jcall(obj, "getdHmsdh", jdouble, (jint, jint, jint), arg0, arg1, arg2)
end

function getd_hmsdk(obj::GHmsjPolynomials, arg0::jint, arg1::jint, arg2::jint)
    return jcall(obj, "getdHmsdk", jdouble, (jint, jint, jint), arg0, arg1, arg2)
end


text:
using SnpParser

using Test

# Single-port Z parameterexample
snp = TouchstoneSnP("example/one_port_z_params")

# Test SnP parameters
@test snp.frequency_unit == "MHz"
@test snp.number_of_ports == 1
@test snp.parameter_type == "Z"
@test snp.number_of_frequencies == 5
@test snp.reference_resistance == 20.0

# Test network data
@test snp.network_data[100.0] == [74.25 -4.0]
@test snp.network_data[200.0] == [60.0 -22.0]
@test snp.network_data[300.0] == [53.025 -45.0]
@test snp.network_data[400.0] == [30.0 -62.0]
@test snp.network_data[500.0] == [0.75 -89.0]


text:
include("Utils.jl")

struct Vertex
    a
    p # x of the vertex
    q # y of the vertex
end

function toexpanded(F::Vertex)
    _validitycheck

    b = -2F.a * F.p
    c = F.a * F.p^2 + F.q

    return Expanded(F.a, b, c)
end

function tofactored(F::Vertex)
    _validitycheck

    under_root = -F.q / F.a

    if (under_root < 0)
        under_root = Complex(under_root)
    end

    x1 = F.p - sqrt(under_root)
    x2 = F.p + sqrt(under_root)

    return Factored(F.a, x1, x2)
end

function calcvalue(F::Vertex, x)
    return F.a * (x - F.p)^2 + F.q
end


text:
using MAVs, NLOptControl, BenchmarkTools, JLD, FileIO, DataFrames

#'''
#This function creates the benchmark tests
#'''
function bench(;baseBench::Bool=false)
  if baseBench
    name="master.jld"
  else
    name="pr.jld"
    master=load("master.jld", "results")
    m1=median(master["psMethods"])
  end
  # create a DataFrame
  fname = joinpath(dirname(@__FILE__), "data.jld")
  df = DataFrames.DataFrame()

  # Define a parent BenchmarkGroup to contain our suite
  const suite = BenchmarkGroup()

  # Add some child groups to our benchmark suite.
  suite["psMethods"] = BenchmarkGroup(["integrationScheme","Nck"])
  for scheme in (:lgrExplicit,:lgrImplicit)
    for Nck in ([10,8,6],[12,10,8,6])
      c=defineCase(;(:mode=>:autoBench));
      setMisc!(c;integrationScheme=scheme,Nck=Nck,max_cpu_time=20.0)
      n=initializeAutonomousControl(c);
      suite["psMethods"][scheme,Nck] = @benchmarkable optimize!($n) samples=10
    end
  end

  paramspath = joinpath(dirname(@__FILE__), "params.jld")
  if isfile(paramspath)
      loadparams!(suite, BenchmarkTools.load(paramspath, "suite"), :evals);
  else
      tune!(suite)
      BenchmarkTools.save(paramspath, "suite", params(suite));
  end
  results=run(suite,verbose=true);              # run all benchmarks
  BenchmarkTools.save(name, "results", results) # save results to JLD file

  m2=median(results["psMethods"]);

  # save the data
  file = jldopen(fname, "w")
  write(file, "df", df)
  close(file)

  if !baseBench
    return m2, judge(m2,m1)
  else
    return m2
  end
end

paramspath = joinpath(Pkg.dir("MAVs/examples/Benchmarks/"), "params.jld")

loadparams!(suite, BenchmarkTools.load(paramspath, "suite"), :evals);

master=load("master.jld", "results")

###########
file = joinpath(Pkg.dir("MAVs/examples/Benchmarks/"), "temp.jld")
res=master;
writeresults(file,res)
#https://github.com/JuliaCI/PkgBenchmark.jl/blob/master/src/runbenchmark.jl
function writeresults(file, res)
    save(File(format"JLD", file), "time", time(), "trials", res)
end

#https://github.com/JuliaCI/PkgBenchmark.jl/blob/master/src/runbenchmark.jl
function readresults(file)
    JLD.jldopen(file,"r") do f
        read(f, "trials")
    end
end

#=
if saveresults
               tosave = if promptsave
                   print("File results of this run? (commit=$(sha[1:6]), resultsdir=$resultsdir) (Y/n) ")
                   response = readline() |> strip
                   response == "" || lowercase(response) == "y"
               else true end
               if tosave
                   !isdir(resultsdir) && mkpath(resultsdir)
                   resfile = joinpath(resultsdir, sha*".jld")
                   writeresults(resfile, res)
                   info("Results of the benchmark were written to $resfile")
end
=#
#=
m=median(run(suite,verbose=true,seconds=10));
save("m2.jld","m",params(m));

# compare
tempDict=load("m1.jld")
m_temp=get(tempDict,"m1",0)

judge(m_temp,m)

##################OLD
#'''
# schemes=[:lgrExplicit,:lgrImplicit]
# m1=benchMark(schemes[1])
# m2=benchMark(schemes[2])
# judge(m1,m2)
#'''
function benchMark(scheme)
  c=defineCase(;(:mode=>:autoBench));
  setMisc!(c;integrationScheme=scheme)
  n=initializeAutonomousControl(c);
  b=@benchmark optimize!($n) samples=10
  return median(b)
end
=#


text:
""" Rewriting for GAT expressions.

The current content of this module is just a stopgap until I can implement
a generic term rewriting system.
"""
module Rewrite
export associate, associate_unit, distribute_unary, involute

using Compat
using ..Syntax

""" Simplify associative binary operation.

Maintains the normal form `op(e1,e2,...)` where `e1`,`e2`,... are expressions
that are *not* applications of `op()`
"""
function associate(expr::E)::E where E <: GATExpr
  op, e1, e2 = head(expr), first(expr), last(expr)
  args1 = head(e1) == op ? args(e1) : [e1]
  args2 = head(e2) == op ? args(e2) : [e2]
  E([args1; args2], type_args(expr))
end

""" Simplify associative binary operation with unit.

Reduces a freely generated (typed) monoid to normal form.
"""
function associate_unit(expr::GATExpr, unit::Function)::GATExpr
  e1, e2 = first(expr), last(expr)
  if (head(e1) == nameof(unit)) e2
  elseif (head(e2) == nameof(unit)) e1
  else associate(expr) end
end

""" Distribute unary operation over binary operation.
"""
function distribute_unary(expr::GATExpr, unary::Function, binary::Function;
                          unit::Union{Function,Nothing}=nothing,
                          contravariant::Bool=false)::GATExpr
  if (head(expr) != nameof(unary)) return expr end
  @assert length(args(expr)) == 1
  arg = first(expr)
  if head(arg) == nameof(binary)
    binary(map(unary, (contravariant ? reverse : identity)(args(arg))))
  elseif !isnothing(unit) && head(arg) == nameof(unit)
    arg
  else
    expr
  end
end

""" Simplify involutive unary operation.
"""
function involute(expr::GATExpr)
  @assert length(args(expr)) == 1
  arg = first(expr)
  head(expr) == head(arg) ? first(arg) : expr
end

end


text:
using Oceananigans.AbstractOperations: KernelFunctionOperation
using Oceananigans.Fields: Field, ZeroField

buoyancy(model) = buoyancy(model.buoyancy, model)

buoyancy(::Nothing, model) = ZeroField()
buoyancy(::BuoyancyTracer, model) = model.tracers.b
buoyancy(b, model) = KernelFunctionOperation{Center, Center, Center}(buoyancy_perturbation, model.grid, computed_dependencies=(b.model, model.tracers))

BuoyancyField(model) = Field(buoyancy(model))


text:
include("runtests.jl")

@time @testset "New Dummy Testing" begin

    Random.seed!(2019)

    # Dynamics
    dt=0.1
    function dynamics(x,u)
        v,ϕ = u
        xp = x
        xp += dt*([v*cos(x[3]);v*sin(x[3]); ϕ])
        xp[3] = xp[3]
        return xp
    end

    W = dt*0.05*Matrix{Float64}(I,3,3)

    dmodel = NonlinearDynamicsModel(dynamics,W)

    # Measurement
    observe(x,u) = [norm(x[1:2])]

    V = 0.1*Matrix{Float64}(I,1,1)

    omodel = NonlinearObservationModel(observe,V)

    # Build UKF
    ukf = UnscentedKalmanFilter(dmodel,omodel);

    # Simulation
    times = 0:dt:5
    action_sequence = [[1.0, 2*sin(pi/2*t)] for t in times]

    b0 = GaussianBelief([0.0,0.0,0.0], 0.01*Matrix{Float64}(I,3,3))

    sim_states, sim_measurements = simulation(b0,action_sequence,ukf);

    # Filter
    filtered_beliefs = run_filter(b0, action_sequence, sim_measurements, ukf)
    μ, Σ = unpack(filtered_beliefs)

    # Tests
    @test ukf.d isa NonlinearDynamicsModel
    @test ukf.o isa NonlinearObservationModel
    @test length(sim_measurements) == length(action_sequence)

    first_update = update(b0, action_sequence[1], sim_measurements[1], ukf)
    first_predict = predict(b0, action_sequence[1], ukf)
    first_predict_measure = measure(first_predict, sim_measurements[1], ukf,
                            u = action_sequence[1])
    array_isapprox(first_update.μ, first_predict_measure.μ)
    array_isapprox(first_update.Σ, first_predict_measure.Σ)

    @test length(filtered_beliefs[end].μ) == length(b0.μ)
    @test size(filtered_beliefs[end].Σ) == size(b0.Σ)
    @test size(μ) == (length(filtered_beliefs), length(filtered_beliefs[1].μ))
    @test size(Σ) == (length(filtered_beliefs), size(filtered_beliefs[1].Σ)...)

    # change with algorithm / rng changes

    array_isapprox(sim_states[end], [0.16459, 3.34923, 2.65721];atol=0.001)
    array_isapprox(sim_measurements[end], [3.17086];atol=0.001)
    array_isapprox(filtered_beliefs[end].μ, [1.83418, 2.78963, 1.37579];
                    atol=0.001)
    array_isapprox(filtered_beliefs[end].Σ, [0.739291 -0.472135 -0.282821;
                                            -0.472135 0.342115 0.197056;
                                            -0.282821 0.197056 0.185713]; atol=0.001)

end

#=
@show sim_states[end]
@show sim_measurements[end]
@show filtered_beliefs[end].μ
@show filtered_beliefs[end].Σ
=#


text:
using SingleCrystal
using Test
using LinearAlgebra
using DataFrames
using JSON

@testset "fcc/bcc unit cell->super cell->vacancy" begin
    element = "Fe"
    mass = 42.
    el2atom_map = Dict(element => SingleCrystal.Atom(name=element, mass=mass))
    a = 2.9
    nx, ny, nz = 2, 3, 4
    for b in ["fcc", "bcc"]
        if b == "fcc"
            unitcell = SingleCrystal.make_fcc_unitcell(element, a)
            num = 4
        elseif b == "bcc"
            unitcell = SingleCrystal.make_bcc_unitcell(element, a)
            num = 2
        end
        @testset "$(b)" begin

            @testset "unitcell: box" begin
                @test unitcell.box.M ≈ a*I(3)
                @test unitcell.edge_lengths ≈ a*ones(3)
            end
            @testset "unitcell: positions & atoms" begin
                @test length(unitcell.positions) == length(unitcell.atoms)
                @test length(unitcell.positions) == num
            end

            supercell = SingleCrystal.make_supercell(unitcell, nx=nx, ny=ny, nz=nz)
            @testset "supercell: box" begin
                @test supercell.box.M[1,1] ≈ nx * unitcell.box.M[1,1]
                @test supercell.box.M[2,2] ≈ ny * unitcell.box.M[2,2]
                @test supercell.box.M[3,3] ≈ nz * unitcell.box.M[3,3]
            end
            @testset "supercell: positions & atoms" begin
                @test length(supercell.positions) == length(supercell.atoms)
                @test length(supercell.positions) == num*2*3*4
            end

            @testset "vacancy" begin
                vac = SingleCrystal.add_vacancies(supercell, ixs=[1])
                @test length(vac.positions) == length(supercell.positions) - 1
            end
        end
    end
end

struct NeighbourFinder{T}
    nb_matrix::BitArray{2} # defines which atom pairs we'll be happy to check at all
    dist_cutoff::T
    rcut2::T
end

NeighbourFinder(nb_matrix, dist_cutoff) = NeighbourFinder{typeof(dist_cutoff)}(nb_matrix, dist_cutoff, dist_cutoff^2)

function find_neighbours(cell::SingleCrystal.Cell, nf::NeighbourFinder)
    neighbours = []
    rs = []
    for i in 1:length(cell.positions)
        ci = cell.positions[i]
        for j in 1:length(cell.positions)
            if i==j
                continue
            end

            r2 = sum(abs2, SingleCrystal.vector(ci, cell.positions[j], cell.edge_lengths))
            if r2 <= nf.rcut2 && nf.nb_matrix[j,i]
                push!(neighbours, (i,j))
                push!(rs, sqrt(r2))
            end
        end
    end
    return neighbours, rs
end

function get_distance_df(cell::SingleCrystal.Cell; dist_cutoff::Real=2)
    n_atoms = length(cell.atoms)
    nb_matrix = trues(n_atoms,n_atoms)
    nf = NeighbourFinder(nb_matrix, dist_cutoff)
    idxs, rs = find_neighbours(cell, nf)
    rs_df = sort(combine(groupby(DataFrame("distances"=>rs),[:distances]),
 nrow=>:count), [:distances])
    return rs_df, rs
end

@testset "Neighbours" begin

    d = 2.
    el = "Fe"
    el2atom_map = Dict(el => SingleCrystal.Atom(name=el, mass=1))

    @testset "bcc" begin
        unitcell = SingleCrystal.make_bcc_unitcell(el, 1.,)
        supercell = SingleCrystal.make_supercell(unitcell, nx=3, ny=3, nz=3)
        rs_df, rs = get_distance_df(supercell, dist_cutoff=d)
        @test rs_df.distances[1] ≈ sqrt((sqrt(2)/2)^2 + 1/2^2)
        @test rs_df.distances[2] ≈ 1
        @test rs_df.distances[3] ≈ sqrt(2)
        @test rs_df.distances[4] ≈ sqrt((sqrt(2)/2)^2 + (3/2)^2)
        @test rs_df.distances[5] ≈ sqrt(sqrt(2)^2 + 1^2)
    end

    @testset "fcc" begin
        unitcell = SingleCrystal.make_fcc_unitcell(el, 1.)
        supercell = SingleCrystal.make_supercell(unitcell, nx=3, ny=3, nz=3)
        rs_df, rs = get_distance_df(supercell, dist_cutoff=d)
        @test rs_df.distances[1] ≈ sqrt(1^2+1^2)/2
        @test rs_df.distances[2] ≈ 1
        @test rs_df.distances[3] ≈ sqrt(1^2+(sqrt(2)/2)^2)
        @test rs_df.distances[4] ≈ sqrt(1^2+1^2)
        @test rs_df.distances[5] ≈ sqrt(3^2+1^2)/2
    end
end

crystal_specs = Dict(
    # NaCl structure
    "NaCl" => Dict(
        "symbols" => ["Na", "Cl"],
        "basis" => [[0. 0. 0.], [.5 .5 .5]], # scaled coordinates
        "nr" => 225,
        "setting" => 1,
        "a_direction" => [1.; 0.; 0.],
        "ab_normal" => [0.; 0.; 1.],
        "cellpar" => [5.64, 5.64, 5.64, 90, 90, 90]
    ),
    # Al fcc structure
    "Al_fcc" => Dict(
        "symbols" => ["Al"],
        "basis" => [[0. 0. 0.],], # scaled coordinates
        "nr" => 225,
        "setting" => 1,
        "a_direction" => [1.; 0.; 0.],
        "ab_normal" => [0.; 0.; 1.],
        "cellpar" => [4.05, 4.05, 4.05, 90, 90, 90]
    ),
    # Fe bcc structure
    "Fe_bcc" => Dict(
        "symbols" => ["Fe"],
        "basis" => [[0. 0. 0.],], # scaled coordinates
        "nr" => 229,
        "setting" => 1,
        "a_direction" => [1.; 0.; 0.],
        "ab_normal" => [0.; 0.; 1.],
        "cellpar" => [2.87, 2.87, 2.87, 90, 90, 90]
    ),
    # Mg hcp structure
    "Mg_hcp" => Dict(
        "symbols" => ["Mg"],
        "basis" => [[1/3 2/3 3/4],], # scaled coordinates
        "nr" => 194,
        "setting" => 1,
        "a_direction" => [1.; 0.; 0.],
        "ab_normal" => [0.; 0.; 1.],
        "cellpar" => [3.21, 3.21, 5.21, 90, 90, 120]
    ),
    # Diamond structure
    "Diamond" => Dict(
        "symbols" => ["C"],
        "basis" => [[0. 0. 0.],], # scaled coordinates
        "nr" => 227,
        "setting" => 1,
        "a_direction" => [1.; 0.; 0.],
        "ab_normal" => [0.; 0.; 1.],
        "cellpar" => [3.57, 3.57, 3.57, 90, 90, 90]
    ),
    # Rutile structure
    "Rutile" => Dict(
        "symbols" => ["Ti", "O"],
        "basis" => [[0. 0. 0.], [.3 .3 0.]], # scaled coordinates
        "nr" => 136,
        "setting" => 1,
        "a_direction" => [1.; 0.; 0.],
        "ab_normal" => [0.; 0.; 1.],
        "cellpar" => [4.6, 4.6, 2.95, 90, 90, 90]
    ),
    # CoSb3 skudderudite
    "Skudderudite" => Dict(
        "symbols" => ["Co", "Sb"],
        "basis" => [[.25 .25 .25], [0. .335 .158]], # scaled coordinates
        "nr" => 204,
        "setting" => 1,
        "a_direction" => [1.; 0.; 0.],
        "ab_normal" => [0.; 0.; 1.],
        "cellpar" => [9.04, 9.04, 9.04, 90, 90, 90]
    )
)

function positions_match_ase(crystal::SingleCrystal.Cell, ase_crystal::Dict)
    return all([p0 ≈ p1 for (p0,p1) in zip(ase_crystal["positions"],crystal.positions)])
end

function cell_matches_ase(crystal::SingleCrystal.Cell, ase_crystal::Dict)
    return crystal.box.M ≈ ase_crystal["cell"]
end

@testset "various single crystals" begin

    json_ase_crystals = SingleCrystal.load_refs()

    for name in keys(crystal_specs)
        @testset "$(name)" begin
            nr = crystal_specs[name]["nr"]
            setting = crystal_specs[name]["setting"]
            basis = crystal_specs[name]["basis"]
            symbols = crystal_specs[name]["symbols"]
            a_direction = crystal_specs[name]["a_direction"]
            ab_normal = crystal_specs[name]["ab_normal"]
            cellpar = crystal_specs[name]["cellpar"]

            crystal = SingleCrystal.make_unitcell(basis, symbols, nr, setting, cellpar,
                                   a_direction=a_direction, ab_normal=ab_normal)

            ase_crystal = SingleCrystal.parse_json_crystal(json_ase_crystals[name])

            @testset "positions match" begin
                @test positions_match_ase(crystal, ase_crystal)
            end
            @testset "cell match" begin
                @test cell_matches_ase(crystal, ase_crystal)
            end
        end
    end
end


text:
using Random, Distributions
using DocStringExtensions

"
A mutable structure for models.
Contains model-level parameters and state variables. Intitialize with values from `global_dict`.
$(TYPEDFIELDS)
"
mutable struct Model
    remove_juveniles::Array{Int64,1} # weekdays on which juveniles are removed from the experiment
    renew_medium::Array{Int64,1} # weekdays on which medium is renewed
    frct_medium_renewal::Float64 # fraction of mefudium that is replaced at each renewal
    timestep::Int64
    tstep::Int64
    tday::Float64
    food_level::Float64
    food_density::Float64
    age_range_init::Int64
    unique_id_count::Int64
    fleas::Array
    fleas_dict::Dict
    stressor_conc::Array{Float64,1}
end

"""
A `toxicity_particle` contains all parameters related to toxic effects. \n
`k_d` is a one-dimenstional Array: There is only one dominant rate constant per stressor. \n
Hazard-related parameters are given in one-dimensional Arrays (one for each stressor). \n
Stress-related parameters are given in two dimensional Arrays (one for each stressor, one for each pMoA). \n
$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct toxicity_particle
    k_d::Array{Float64,1}
    S_max::Array{Float64,2}
    S_ED50::Array{Float64,2}
    S_beta::Array{Float64,2}
    h_max::Array{Float64,1}
    h_ED50::Array{Float64,1}
    h_beta::Array{Float64,1}
end

"
A mutable structure for individuals. \n
Contains individual-level parameters and state variables. Initialize with values from `fleas_dict` using function `make_flea`.
$(TYPEDEF)
"
mutable struct Flea
    unique_id::Int64
    species::String
    cv::Float64

    # correction factors
    T_fact::Float64

    # scaled DEB parameters
    K::Float64
    J_XAm_rate::Float64
    U_Hb::Float64
    U_Hp::Float64
    k_M_rate::Float64
    k_J_rate::Float64
    g::Float64
    v_rate::Float64
    kap::Float64
    kap_R::Float64
    shape_factor::Float64
    mortality_constant::Float64
    starvation_hazard_rate::Float64
    sG::Float64
    h_a::Float64
    L0::Float64

    # embryonal parameters
    egg_weight::Float64
    L_embryo::Float64
    U_E_embryo::Float64

    # auxiliary
    time_between_molts::Float64

    crit_mass::Float64

    # toxicity parameters
    tox_particle::toxicity_particle

    # state variables
    die::Bool
    cause_of_death::String
    juvenile::Bool
    adult::Bool
    offspring_number::Int64
    cum_repro::Int64
    final_offspring::Int64
    molt_time::Float64
    develop_time::Float64
    U_E::Float64
    dU_E::Float64
    U_H::Float64
    dU_H::Float64
    U_R::Float64
    dU_R::Float64
    h_rate::Float64
    dh_rate::Float64
    q_acceleration::Float64
    dq_acceleration::Float64
    L::Float64
    dL::Float64
    L_max::Float64
    Xi::Float64
    functresp::Float64
    food_taken::Float64
    S_A::Float64
    S_C::Float64
    e_scaled::Float64

    # _int values for everything that can be changed by a stressor
    J_XAm_rate_int::Float64
    k_J_rate_int::Float64
    g_int::Float64
    k_M_rate_int::Float64
    kap_R_int::Float64

    # toxicity-related state variables per stressor
    D::Array{Float64,1} # damage
    h_z::Array{Float64,1} # hazard rate
    p_z::Array{Float64,1} # survival probability
    stress::Array{Float64,2} # stress

    # toxicity-related state variables per pMoA
    f_stress::Float64
    g_stress::Float64
    k_M_stress_1::Float64
    k_M_stress_2::Float64
    k_J_stress::Float64
    repro_stress_1::Float64
    repro_stress_2::Float64

    # individual modifier
    scatter_multiplier::Float64
end

function make_tox_particle(tox_particle_dict::Dict)
    tox_particle = toxicity_particle(
    tox_particle_dict["k_d"],
    tox_particle_dict["S_max"],
    tox_particle_dict["S_ED50"],
    tox_particle_dict["S_beta"],
    tox_particle_dict["h_max"],
    tox_particle_dict["h_ED50"],
    tox_particle_dict["h_beta"]
    )
    return tox_particle
end

"
A function to initialize a Flea. \n
Choose start `juvenile=1` or start `adult=1` for initial individuals.
$(SIGNATURES)
"
function make_flea(
        m::Model,
        start_juvenile::Bool,
        start_adult::Bool
    )

    fleas_dict = m.fleas_dict

    # generate Flea object
    f = Flea(
        m.unique_id_count,
        fleas_dict["species"],
        fleas_dict["cv"],
        # correction factors
        1.,

        # individual-level parameters
        fleas_dict["K"],
        fleas_dict["J_XAm_rate"],
        fleas_dict["U_Hb"],
        fleas_dict["U_Hp"],
        fleas_dict["k_M_rate"],
        fleas_dict["k_J_rate"],
        fleas_dict["g"],
        fleas_dict["v_rate"],
        fleas_dict["kap"],
        fleas_dict["kap_R"],
        fleas_dict["shape_factor"],
        fleas_dict["mortality_constant"],
        fleas_dict["starvation_hazard_rate"],
        fleas_dict["sG"],
        fleas_dict["h_a"],
        fleas_dict["L0"],

        # embryonal parameters
        fleas_dict["egg_weight"],
        fleas_dict["L_embryo"],
        fleas_dict["U_E_embryo"],

        # other
        fleas_dict["time_between_molts"],
        fleas_dict["crit_mass"],

        # toxicity parameters
        make_tox_particle(fleas_dict["toxicity_particle"]),

        # state variables
        false, # die?
        "None",
        0, # juvenile
        0, # adults
        0, # offspring_number
        0, # final_offspring
        0, # cum_repro
        0., # molt_time
        0., # develop_time
        0., # U_E
        0., # dU_E
        0., # U_H
        0., # dU_H
        0., # U_R
        0., # dU_R
        0., # h_rate
        0., # dh_rate
        0., # q_acceleration
        0., # dq_acceleration
        0., # L
        0., # dL
        0., # L_max
        0., # Xi
        0., # functresp
        0., # food_taken
        0., # S_A
        0., # S_C
        0., # e_scaled

        # _int values
        fleas_dict["J_XAm_rate"],
        fleas_dict["k_J_rate"],
        fleas_dict["g"],
        fleas_dict["k_M_rate"],
        fleas_dict["kap_R"],

        # toxicity-related state variables per stressor
        [0., 0., 0.], # damage
        [0., 0., 0.], # hazard rate
        [1., 1., 1.], # survival probability
        repeat([0., 0., 1., 0.], 1, 3), # stress

        # toxicity-related state variables per pMoA
        1., # f_stress
        0., # g_stress
        0., # k_M_stress_1
        0., # k_M_stress_2
        0., # k_J_stress
        0., # repro_stress_1
        0., # repro_stress_2

        # individual variability
        0. # scatter_multiplier
    )
    # TODO at some point: implement temperature correction
    f.T_fact = 1

    f = individual_variability(f)

    if start_juvenile
        f = init_random_age(f, m)
    elseif start_adult
        f = init_adulthood(f, m)
    end

    return f
end

"""
Initialize a juvenile with age <= age _ range _ init. Age is given in units of model timesteps.
Initial juveniles are typically 0-24 hours old.
$(SIGNATURES)
"""
function init_random_age(f::Flea, m::Model)
    f.develop_time = -1
    f.juvenile = false

    age = rand(1:m.age_range_init)

    for a in 1:age
        f.functresp = 1
        f, m = feed(f, m, false)
        f = dU_E(f)
        f = dU_H(f)
        f = dU_R(f)
        f = dL(f)
        f = update(f, f.cv, m.timestep)
    end
    return f
end

"
Simulate the live of an adult before the start of an experiment.
It is assumed that the individual could always feed ad libitum, therefore the scaled functional response
is fixed to 1.
Investment in reproduction is accounted for and individuals have the potential to produce offspring immediately after the onset of the experiment.

Usage:
```
f = init_adulthood(f, m)
```
"
function init_adulthood(f::Flea, m::Model)
    init_time = 0
    while (f.shape_factor*10<3) & (init_time < m.timestep*20)
        init_time +=1
        f,m = feed(f, m, false)
        f = dU_E(f)
        f = dU_H(f)
        f = dU_R(f)
        f = dL(f)
        f = update(f, f.cv, m.timestep)

        if f.adult
            f.molt_time = f.molt_time + (f.T_fact/m.timestep)

            # emulate molting + reproduction
            if f.molt_time > f.time_between_molts
                f.molt_time = 0
                f.final_offspring = Int64(floor((f.U_R/f.kap_R)/(f.egg_weight/f.T_fact)))

                # TODO:
                # this line is not included in netlogo model
                # but it makes sense?
                # think about whether to include this
                f.U_R = f.U_R - ((floor(f.U_R/(f.egg_weight/f.T_fact))) * (f.egg_weight/f.T_fact))
            end
        end

    # TODO
    # figure out whether this is really necessary (relict of netlogo model)
    # annotate and remove magic numbers
    f.molt_time =  -1*rand(1:4)
    f.U_R = (f.egg_weight/f.T_fact) * (15 * f.molt_time) / f.time_between_molts

    end
    return f
end

"""
Feeding of individuals.
Returns updated model-level (food density) and individual-level (conditional food density, functional response) state variables.
$(SIGNATURES)
"""
function feed(;
        U_H::Float64,
        U_Hb::Float64,
        K::Float64,
        L::Float64,
        T_fact::Float64,
        J_XAm_rate::Float64,
        f_stress::Float64,
        timestep::Int64,
        food_density::Float64,
        update_food_density::Bool
    )
    # Embryos don't feed
    if U_H <= U_Hb
        functresp = 0
        Xi = 0
    # Juveniles and adults do
    else
        Xi = food_density
        # functional response
        functresp = Xi/(K + Xi)
    end

    # to convert scaled functional response to amount of ingested food
    # - multiply with surface area: the result is identical with assimilation flux S_A,
    # - multiply with the surface area-specific maximum ingestion rate
    # - multiply with temperature correction factor
    food_taken = (functresp * L^2 * T_fact * J_XAm_rate) / timestep

    if update_food_density
        # subtract ingested food from environmental food density
        # max() is applied just to make sure that floating point errors don't mess things up
        # (if food_density ever becomes negative, the error will propagate almost irreversibly)
        food_density = max(0, food_density - food_taken)
    end

    return Xi, functresp, food_taken, food_density
end

function feed(f::Flea, m::Model, update_food_density::Bool)
    f.Xi, f.functresp, f.food_taken, m.food_density = feed(
            U_H          = f.U_H,
            U_Hb         = f.U_Hb,
            K            = f.K,
            L            = f.L,
            T_fact       = f.T_fact,
            J_XAm_rate   = f.J_XAm_rate,
            f_stress     = f.f_stress,
            timestep     = m.timestep,
            food_density = m.food_density,
            update_food_density = update_food_density)
    return f, m
end

"""
Induction of individual variability, based on scatter_multiplier.
"""
function individual_variability(f::Flea)
    f.scatter_multiplier = exp(rand(Normal(0, f.cv)))

    f.g    = f.g / f.scatter_multiplier
    f.U_Hb = f.U_Hb / f.scatter_multiplier
    f.U_Hp = f.U_Hp / f.scatter_multiplier
    f.K    = f.K * f.scatter_multiplier
    f.J_XAm_rate = f.J_XAm_rate * f.scatter_multiplier
    f.L   = f.L_embryo * f.scatter_multiplier
    f.U_E = f.U_E_embryo * f.scatter_multiplier
    f.juvenile = true
    f.adult    = false
    f.offspring_number = 0
    f.U_R  = 0
    f.dU_R = 0
    f.U_H  = f.U_Hb
    f.dU_H = 0
    f.h_rate  = 0
    f.dh_rate = 0
    f.q_acceleration  = 0
    f.dq_acceleration = 0

    # store _int values for everything that can be changed by a stressor
    f.k_J_rate_int = f.k_J_rate
    f.g_int        = f.g
    f.k_M_rate_int = f.k_M_rate
    f.kap_R_int    = f.kap_R

    return f
end

"""
Calculate change in scaled maturity.
"""
function dU_H(f::Flea)
    if f.U_H < f.U_Hp/f.T_fact
        f.dU_H = (1-f.kap) * f.S_C - (f.T_fact * f.k_J_rate) * f.U_H
    else
        f.dU_H = 0
    end
    return f
end

"""
Calculate change in scaled reserves.
Also: Assimilation flux, mobilization flux, scaled reserve density.
"""
function dU_E(f::Flea)
    if f.juvenile
        f.S_A      = f.functresp * f.L^2 * f.f_stress
        f.e_scaled = f.T_fact * f.v_rate * (f.U_E/(f.L^3))
        f.S_C      = f.L^2 * (f.g*f.e_scaled/(f.g+f.e_scaled)) * (1+(f.L/(f.g*(f.T_fact*f.v_rate/(f.g*f.T_fact*f.k_M_rate)))))
        f.dU_E     = f.S_A - f.S_C
    end
    return f
end

"""
Calculate change in reproduction buffer
"""
function dU_R(f::Flea)
    if f.U_H > f.U_Hp
        f.dU_R = (1-f.kap) * f.S_C - (f.T_fact * f.k_J_rate) * (f.U_Hp / f.T_fact)
    else
        f.dU_R = 0
    end
    return f
end

"""
Calculate change in scaled size.
"""
function dL(f::Flea)

    f.dL = (1/3) *
    (((f.T_fact * f.v_rate)/(f.g*(f.L^2)) * f.S_C) - (
            f.T_fact * f.k_M_rate) * f.L)

    # starvation if growth is negative
    """
    if f.e_scaled < f.L / (f.v_rate / (f.g * f.k_M_rate))
        f.dL = 0

        # juveniles divert maturation energy into maintenance
        if f.dU_H < (f.U_Hp / f.T_fact)
            f.dU_H = (1-f.kap) * f.e_scaled * f.L^2 - f.k_J_rate * f.U_Hp - f.kap * f.L^2 * (f.L / (f.v_rate / (f.g * f.k_M_rate)) - f.e_scaled)
        # adults divert reproduction energy into maintenance
        else
            f.dU_R = (1-f.kap) * f.e_scaled * f.L^2 - f.k_J_rate * f.U_Hp - f.kap * f.L^2 * (f.L / (v.v_rate / (f.g * f.k_M_rate)) - f.e_scaled)
        end

        f.dU_E = f.S_A - f.e_scaled * f.L ^ 2

        if f.U_H < (f.U_Hp / f.T_fact)
            if f.dU_H < 0
                f.die = true
                end
        else
            if f.dU_R < 0
                f.die = true
            end
        end
    end
    """

    return f
end

"""
Calculate ageing acceleration
"""
function dq_acceleration(f::Flea)
    f.dq_acceleration = (f.q_acceleration * ((f.L^3)/(((f.T_fact*f.v_rate)/(f.g*(f.T_fact*f.k_M_rate)))^3))*f.sG+f.h_a)*f.e_scaled*((f.T_fact*f.v_rate/f.L)-((3/f.L)*f.dL))-((3/f.L)*f.dL)*f.q_acceleration
    return f
end

function dh_rate(f::Flea)
    f.dh_rate = f.q_acceleration - ((3/f.L)*f.dL)*f.h_rate
    return f
end

function update(f::Flea, cv::Float64, timestep::Int64)
    # update main state variables
    f.U_E = f.U_E + f.dU_E/timestep
    f.U_H = f.U_H + f.dU_H/timestep
    f.U_R = max(f.U_R + f.dU_R/timestep, 0)
      f.L = max(0, f.L + f.dL/timestep)
    f.L_max = max(f.L, f.L_max)
    # update ageing and hazard rate
    if f.U_H > f.U_Hb
        f.q_acceleration = f.q_acceleration + f.dq_acceleration/timestep
        f.h_rate = f.h_rate + f.dh_rate/timestep
    end
    # transition from embryo to juvenile
    if f.U_H >= f.U_Hb/f.T_fact
        f.juvenile = true
    end
    # transition from juvenile to adult
    if f.U_H > f.U_Hp/f.T_fact
        f.adult = true
    end
    # check development status
    if f.juvenile==false
        f.develop_time = f.develop_time - f.T_fact
        # introduce individual variability in embryos
        if f.develop_time <= 0
            f = individual_variability(f)
        end
    end
    return f
end

function reset(f::Flea)
    f.adult    = 0
    f.juvenile = 0
    f.offspring_number = 0
    f.L    = 0
    f.U_E  = 0
    f.U_H  = 0
    f.U_R  = 0
    f.dU_R = 0
    f.h_rate  = 0
    f.dh_rate = 0
    f.q_acceleration  = 0
    f.dq_acceleration = 0
    f.S_A = 0
    f.S_C = 0
    f.dL  = 0
    f.L_max = 0 # maximum structural size
    f.develop_time = f.T_fact - 1
    f.die = false
    return f
end

"""
Inheritance attributes.
$(SIGNATURES)
"""
function inherit(parent::Flea, child::Flea)
    child.g    = parent.g
    child.U_Hb = parent.U_Hb
    child.U_Hp = parent.U_Hp
    child.K    = parent.K
    child.J_XAm_rate = parent.J_XAm_rate
    child.L_embryo   = parent.L_embryo
    child.U_E_embryo = parent.U_E_embryo

    return child
end

function reproduce(f::Flea, m::Model)
    # molting controls timing of reproduction
    f.molt_time = f.molt_time + (f.T_fact/m.timestep)
    if f.molt_time > f.time_between_molts
        # determine no. of eggs based on reproduction buffer and reproduction efficiency
        num_eggs = (f.U_R*f.kap_R)/(f.egg_weight/f.T_fact)
        # convert to integer
        if isnan(num_eggs)
            f.offspring_number = 0
        else
            f.offspring_number = Int64(floor(num_eggs))
            f.cum_repro = f.cum_repro + Int64(floor(num_eggs))
        end
        # reset molting
        f.molt_time = 0
        # the subtracted amount of reserve does not take kap_R into account
        # (energy is spent, but not effectively turned into offspring)
        f.U_R = f.U_R - ((floor(f.U_R/(f.egg_weight/f.T_fact))) * (f.egg_weight/f.T_fact))
        # create new offspring
        for i in 1:f.offspring_number
            NewFlea = make_flea(m, false, false)
            NewFlea = inherit(f, NewFlea)
            NewFlea = reset(NewFlea)
            NewFlea = individual_variability(NewFlea)
            m.unique_id_count += 1
            # add new flea to fleas array
            push!(m.fleas, NewFlea)
        end
    end
    return f, m
end

function death(f::Flea, m::Model)
    # death by starvation
    condition_1a = (f.L^3) < f.crit_mass * (f.L_max^3)
    condition_1b = rand() < 1 - ((1-f.starvation_hazard_rate)^(f.T_fact/m.timestep))

    if condition_1a & condition_1b
        f.cause_of_death = "starvation"
        f.die = true
    end

    # death by ageing
    f.h_rate = min(f.h_rate, 1)
    condition_2 = rand() < 1-((1-f.h_rate)^(1/m.timestep))

    if condition_2
        f.cause_of_death = "age"
        f.die = true
    end

    # experimental removal of animals
    condition_3a = true in [floor(m.tday)%7==x for x in m.remove_juveniles]
    condition_3b = f.unique_id != 0

    if condition_3a & condition_3b
        f.cause_of_death = "removal"
        f.die = true
    end

    return f
end

"""
Hazard rate dependent on scaled reserve density.
$(SIGNATURES)
"""
function density_dependent_mortality(f::Flea, m::Model)
    if (f.juvenile==true) & (f.adult==false)
        mortality_probability = (1 - f.e_scaled) * (1-(1-f.mortality_constant)^(1/m.timestep))
        if rand() < mortality_probability
            f.cause_of_death = "density"
            f.die = true
        end
    end
    return f
end

"""
change in damage
$(SIGNATURES)
"""
function dD(D, k_d, stressor_conc)
    return D + k_d * (stressor_conc - D)
end

"""
log-logistic response
$(SIGNATURES)
"""
function log_logistic(x, upper_limit, EC50, slope)
    return upper_limit * (1/(1+(x/EC50)^slope))
end

"""
Lethal and sublethal effects of an arbitrary number of stressors. All stressors can have all pMoAs simultaneously.
$(SIGNATURES)
"""
function stressor_toxicity(f::Flea, m::Model)

    #### Step 1: Reset physiological stress ####
    f.f_stress = 1
    f.g_stress = 0
    f.k_M_stress_1 = 0
    f.k_J_stress = 0
    f.repro_stress_1 = 0
    f.repro_stress_2 = 0

    #### Step 2: update scaled damage ####
    f.D = max.(0, dD.(f.D, f.tox_particle.k_d, m.stressor_conc))

    #### Step 3: lethal effects ####
    ## Step 3.1: calculate hazard rate ##
    f.h_z = log_logistic.(
        f.D,
        f.tox_particle.h_max,
        f.tox_particle.h_ED50,
        .-f.tox_particle.h_beta
    )

    ## Step 3.2: convert to survival probability ##
    f.p_z = exp.(.- f.h_z ./m.timestep)

    ## Step 3.3: stochastic death ##
    dice = rand(length(f.p_z))
    if sum(dice .> f.p_z) .>=1
        f.die=true
        f.cause_of_death="stressor"
    end

    # apply log-logistic function for every stressor
    f.stress = log_logistic.(
        f.D',
        f.tox_particle.S_max,
        f.tox_particle.S_ED50,
        .-f.tox_particle.S_beta
    )

    #### Step 4: calculate mixture stress for every pMoA ####
    # Growth
    f.g_stress     = prod(1 .+ f.stress[1,:])
    f.k_M_stress_2 = prod(1 ./ (1 .+ f.stress[1,:]))

    # Maintenance
    f.k_J_stress   = prod(1 .+ f.stress[2,:])
    f.k_M_stress_1 = prod(1 .+ f.stress[2,:])

    # Assimilation
    f.f_stress = prod(1 .- f.stress[3,:])

    # Reproduction
    f.repro_stress_1 = prod(1 ./ (1 .+ f.stress[4,:]))

    #### Step 5: apply stress to physiological parameters ####
    # note that f_stress ins applied in `dU_E()`
    f.g = f.g_int * f.g_stress
    f.k_M_rate = f.k_M_rate_int * f.k_M_stress_1 * f.k_M_stress_2
    f.k_J_rate = f.k_J_rate_int * f.k_J_stress
    f.kap_R = f.kap_R_int * f.repro_stress_1

    return f
end

"""
Update environmental state variables. \n
Includes renewal of medium (=removal of food) and addition of food.
$(SIGNATURES)
"""
function update_environment(m::Model)
    # add food daily
    if m.tstep % m.timestep == 0
        # renew media
        if true in [floor(m.tday)%7==x for x in m.renew_medium]
             m.food_density = m.food_density * (1 - m.frct_medium_renewal)
        end
        m.food_density = m.food_density + m.food_level * exp(rand(Normal(0.0, 0.2)))
    end
    return m
end

"""
Execute a single model step.
"""
function Step(m::Model)
    # randomize order in which individuals are called
    shuffle!(m.fleas)
    # execute individual schedule
    for f in m.fleas

        f = stressor_toxicity(f, m)
        f, m = feed(f, m, true)
        if f.juvenile
            f = f |>
            dU_E |>
            dU_H |>
            dU_R |>
            dL |>
            dq_acceleration |>
            dh_rate
        end
        f = density_dependent_mortality(f, m)
        f = update(f, f.cv, m.timestep)
        if f.adult
            f, m = reproduce(f, m)
        end
        f = death(f, m)
    end
    # add food, renew media etc
    m = update_environment(m)
    # remove dead individuals
    filter!(f->f.die==false, m.fleas)
    m.tstep +=1
    m.tday = m.tstep/m.timestep
    return m
end


text:
include("autotuning_utils.jl")
include("zygote_rules.jl")
include("forwarddiff_rules.jl")

function update_hyperparameters!(m::GP, state, x, y)
    μ₀ = pr_mean(m.f) # Get prior means
    k = kernel(m.f) # Get kernels
    if !isnothing(opt(m.f))
        if ADBACKEND[] == :Zygote
            hp_state = state.hyperopt_state
            # Compute gradients
            Δμ₀, Δk = Zygote.gradient(μ₀, k) do μ₀, k # Compute gradients for the whole model
                ELBO(m, x, y, μ₀, k)
            end

            # Optimize prior mean
            hp_state = state.hyperopt_state
            if !isnothing(Δμ₀)
                hp_state = update!(μ₀, Δμ₀, hp_state)
            end

            # Optimize kernel parameters
            if isnothing(Δk)
                @warn "Kernel gradients are equal to zero" maxlog = 1
            else
                state_k = update_kernel!(opt(m.f), kernel(m.f), Δk, hp_state)
                hp_state = merge(hp_state, (; state_k))
            end
            state = merge(state, (; hyperopt_state=hp_state))
        elseif ADBACKEND[] == :ForwardDiff
            θ, re = destructure((μ₀, k))
            Δ = ForwardDiff.gradient(θ) do θ
                ELBO(m, re(θ)...)
            end
        end
    end
    # end
    return state
end

# @traitfn function update_hyperparameters!(
#     m::TGP,
# ) where {TGP <: AbstractGPModel; IsFull{TGP}}
#     update_hyperparameters!.(m.f, [xview(m)])
#     setHPupdated!(m.inference, true)
# end

@traitfn function update_hyperparameters!(
    m::TGP, state, x, y
) where {TGP <: AbstractGPModel; IsFull{TGP}}
    if any((!) ∘ isnothing ∘ opt, m.f) # Check there is a least one optimiser
        hp_state = state.hyperopt_state
        μ₀ = pr_means(m) # Get prior means
        ks = kernels(m) # Get kernels
        if ADBACKEND[] == :Zygote
            Δμ₀, Δk = Zygote.gradient(μ₀, ks) do μ₀, ks # Compute gradients for the whole model
                ELBO(m, x, y, μ₀, ks, state)
            end
            # Optimize prior mean
            if !isnothing(Δμ₀)
                hp_state = update!.(μ₀, Δμ₀, hp_state)
            end
            if isnothing(Δk)
                @warn "Kernel gradients are equal to zero" maxlog = 1
            else
                hp_state = map(m.f, Δk, hp_state) do gp, Δ, hp_st
                    if isnothing(opt(gp))
                        return hp_st
                    else
                        state_k = update_kernel!(opt(gp), kernel(gp), Δ, hp_st.state_k)
                        return merge(hp_st, (; state_k))
                    end
                end
            end
        elseif ADBACKEND[] == :ForwardDiff
            θ, re = destructure((μ₀, ks))
            Δ = ForwardDiff.gradient(θ) do θ
                ELBO(m, re(θ)...)
            end
        end
        state = merge(state, (; hyperopt_state=hp_state))
    end
    return state
end

@traitfn function update_hyperparameters!(
    m::TGP, state, x, y
) where {TGP <: AbstractGPModel; !IsFull{TGP}}
    # Check that here is least one optimiser
    if any((!) ∘ isnothing ∘ opt, m.f) || any((!) ∘ isnothing ∘ Zopt, m.f)
        hp_state = state.hyperopt_state
        μ₀ = pr_means(m)
        ks = kernels(m)
        Zs = Zviews(m)
        if ADBACKEND[] == :Zygote
            Δμ₀, Δk, ΔZ = Zygote.gradient(μ₀, ks, Zs) do μ₀, ks, Zs
                ELBO(m, x, y, μ₀, ks, Zs, state)
            end
            # Optimize prior mean
            if !isnothing(Δμ₀)
                hp_state = update!.(μ₀, Δμ₀, hp_state)
            end

            # Optimize kernel parameters
            if isnothing(Δk)
                @warn "Kernel gradients are equal to zero" maxlog = 1
            else
                hp_state = map(m.f, Δk, hp_state) do gp, Δ, hp_st
                    if !isnothing(opt(gp))
                        state_k = update_kernel!(opt(gp), kernel(gp), Δ, hp_st.state_k)
                        return merge(hp_st, (; state_k))
                    else
                        return hp_st
                    end
                end
            end

            # Optimize inducing point locations
            if isnothing(ΔZ)
                @warn "Inducing point locations gradients are equal to zero" maxlog = 1
            else
                hp_state = map(m.f, ΔZ, hp_state) do gp, Δ, hp_st
                    if !isnothing(Zopt(gp))
                        state_Z = update_Z!(Zopt(gp), Zview(gp), Δ, hp_st.state_Z)
                        return merge(hp_st, (; state_Z))
                    else
                        return hp_st
                    end
                end
            end
        elseif ADBACKEND[] == :ForwardDiff
            θ, re = destructure((μ₀, ks, Zs))
            Δ = ForwardDiff.gradient(θ) do θ
                ELBO(m, re(θ)...)
            end
        end
        state = merge(state, (; hyperopt_state=hp_state))
    end
    return state
end

## Update all hyperparameters for the full batch GP models ##
function update_hyperparameters!(gp::AbstractLatent, X::AbstractVector)
    if !isnothing(gp.opt)
        f_l, f_μ₀ = hyperparameter_gradient_function(gp, X)
        ad_backend = K_ADBACKEND[] == :auto ? ADBACKEND[] : K_ADBACKEND[]
        Δμ₀ = f_μ₀()
        Δk = if ad_backend == :forward
            ∇L_ρ_forward(f_l, gp, X)
        elseif ad_backend == :zygote
            ∇L_ρ_zygote(f_l, gp, X)
        end
        apply_Δk!(gp.opt, kernel(gp), Δk) # Apply gradients to the kernel parameters
        apply_gradients_mean_prior!(pr_mean(gp), Δμ₀, X)
    end
end

## Update all hyperparameters for the sparse variational GP models ##
function update_hyperparameters!(
    gp::AbstractLatent,
    X::AbstractVector,
    ∇E_μ::AbstractVector,
    ∇E_Σ::AbstractVector,
    i::AbstractInference,
    vi_opt::InferenceOptimizer,
)
    Δμ₀, Δk = if !isnothing(gp.opt)
        f_ρ, f_Z, f_μ₀ = hyperparameter_gradient_function(gp)
        ad_backend = K_ADBACKEND[] == :auto ? ADBACKEND[] : K_ADBACKEND[]
        Δμ₀ = f_μ₀()
        Δk = if ad_backend == :forward
            ∇L_ρ_forward(f_ρ, gp, X, ∇E_μ, ∇E_Σ, i, vi_opt)
        elseif ad_backend == :zygote
            ∇L_ρ_zygote(f_ρ, gp, X, ∇E_μ, ∇E_Σ, i, vi_opt)
        end
        (Δμ₀, Δk)
    else
        nothing, nothing
    end
    if !isnothing(Zopt(gp))
        ad_backend = Z_ADBACKEND[] == :auto ? ADBACKEND[] : Z_ADBACKEND[]
        Z_grads = if ad_backend == :forward
            Z_gradient_forward(gp, f_Z, X, ∇E_μ, ∇E_Σ, i, vi_opt) #Compute the gradient given the inducing points location
        elseif ad_backend == :zygote
            Z_gradient_zygote(gp, f_Z, X, ∇E_μ, ∇E_Σ, i, vi_opt)
        end
        update_Z!(opt(gp.Z), gp.Z, Z_grads) #Apply the gradients on the location
    end
    if !all([isnothing(Δk), isnothing(Δμ₀)])
        apply_Δk!(gp.opt, kernel(gp), Δk) # Apply gradients to the kernel parameters
        apply_gradients_mean_prior!(pr_mean(gp), Δμ₀, X)
    end
end

function update_hyperparameters!(
    gp::AbstractLatent,
    l::AbstractLikelihood,
    i::AbstractInference,
    X::AbstractVector,
    Y::AbstractVector,
)
    Δμ₀, Δk = if !isnothing(gp.opt)
        ad_backend = K_ADBACKEND[] == :auto ? ADBACKEND[] : K_ADBACKEND[]
        (Δμ₀, Δk) = if ad_backend == :forward
            ∇L_ρ_forward(f_ρ, gp, X, ∇E_μ, ∇E_Σ, i, vi_opt)
        elseif ad_backend == :zygote
            ∇L_ρ_zygote(gp, l, i, X, Y)
        end
    else
        nothing, nothing
    end
    if !isnothing(opt(gp.Z))
        ad_backend = Z_ADBACKEND[] == :auto ? ADBACKEND[] : Z_ADBACKEND[]
        global Z_grads = if ad_backend == :forward
            Z_gradient_forward(gp, f_Z, X, ∇E_μ, ∇E_Σ, i, vi_opt) #Compute the gradient given the inducing points location
        elseif ad_backend == :zygote
            Z_gradient_zygote(gp, l, i, X, Y)
        end
        update!(opt(gp.Z), gp.Z.Z, Z_grads) #Apply the gradients on the location
    end
    if !all([isnothing(Δk), isnothing(Δμ₀)])
        apply!(kernel(gp), Δk, gp.opt) # Apply gradients to the kernel parameters
        update!(pr_mean(gp), Δμ₀, X)
    end
end

## Return the derivative of the KL divergence between the posterior and the GP prior ##
function hyperparameter_KL_gradient(J::AbstractMatrix, A::AbstractMatrix)
    return trace_ABt(J, A) / 2
end

function hyperparameter_gradient_function(gp::LatentGP{T}, ::AbstractVector) where {T}
    A = (inv(cov(gp)) - mean(gp) * transpose(mean(gp))) # μ = inv(K+σ²)*(y-μ₀)
    return (function (Jnn)
        return -hyperparameter_KL_gradient(Jnn, A)
    end, function ()
        return -mean(gp)
    end)
end

## Return functions computing gradients of the ELBO given the kernel hyperparameters for a non-sparse model ##
function hyperparameter_gradient_function(
    gp::VarLatent{T}, X::AbstractVector
) where {T<:Real}
    μ₀ = pr_mean(gp, X)
    A =
        (I - pr_cov(gp) \ (cov(gp) + (mean(gp) - μ₀) * transpose(mean(gp) - μ₀))) /
        pr_cov(gp)
    return (
        function (Jnn)
            return -hyperparameter_KL_gradient(Jnn, A)
        end,
        function ()
            return pr_cov(gp) \ (mean(gp) - μ₀)
        end,
    )
end

## Return functions computing gradients of the ELBO given the kernel hyperparameters for a non-sparse latent GP ##
function hyperparameter_gradient_function(gp::SparseVarLatent{T}) where {T<:Real}
    μ₀ = pr_mean(gp, gp.Z)
    A =
        (I(dim(gp)) - pr_cov(gp) \ (cov(gp) + (mean(gp) - μ₀) * transpose(mean(gp) - μ₀))) /
        cov(gp)
    κΣ = gp.κ * cov(gp)
    return (
        function (Jmm, Jnm, Jnn, ∇E_μ, ∇E_Σ, i, opt)
            return (
                hyperparameter_expec_gradient(gp, ∇E_μ, ∇E_Σ, i, opt, κΣ, Jmm, Jnm, Jnn) -
                hyperparameter_KL_gradient(Jmm, A)
            )
        end,
        function (Jmm, Jnm, ∇E_μ, ∇E_Σ, i, opt)
            return hyperparameter_expec_gradient(
                gp, ∇E_μ, ∇E_Σ, i, opt, κΣ, Jmm, Jnm, zero(gp.K̃)
            ) - hyperparameter_KL_gradient(Jmm, A)
        end,
        function ()
            return pr_cov(gp) \ (mean(gp) - μ₀)
        end,
    )
end

function hyperparameter_gradient_function(
    gp::TVarLatent{T}, X::AbstractVector
) where {T<:Real}
    μ₀ = pr_mean(gp, X)
    A =
        (
            I(dim(gp)) .-
            (pr_cov(gp)) \ (cov(gp) .+ (mean(gp) - μ₀) * transpose(mean(gp) - μ₀))
        ) / (pr_cov(gp))
    return (
        function (Jnn)
            return -hyperparameter_KL_gradient(prior(gp).χ * Jnn, A)
        end,
        function ()
            return (pr_cov(gp)) \ (mean(gp) - μ₀)
        end,
    )
end

function hyperparameter_gradient_function(gp::OnlineVarLatent{T}) where {T<:Real}
    μ₀ = pr_mean(gp, gp.Z)
    A =
        (I(dim(gp)) - pr_cov(gp) \ (cov(gp) + (mean(gp) - μ₀) * transpose(mean(gp) - μ₀))) /
        pr_cov(gp)
    κΣ = gp.κ * cov(gp)
    κₐΣ = gp.κₐ * cov(gp)
    return (
        function (Jmm, Jnm, Jnn, Jab, Jaa, ∇E_μ, ∇E_Σ, i, opt)
            ∇E = hyperparameter_expec_gradient(gp, ∇E_μ, ∇E_Σ, i, opt, κΣ, Jmm, Jnm, Jnn)
            ∇KLₐ = hyperparameter_online_gradient(gp, κₐΣ, Jmm, Jab, Jaa)
            ∇KL = hyperparameter_KL_gradient(Jmm, A)
            return ∇E + ∇KLₐ - ∇KL
        end, # Function gradient given kernel parameters
        function (Jmm, Jnm, Jab, ∇E_μ, ∇E_Σ, i, opt)
            return hyperparameter_expec_gradient(
                gp, ∇E_μ, ∇E_Σ, i, opt, κΣ, Jmm, Jnm, zero(gp.K̃)
            ) + hyperparameter_online_gradient(
                gp, κₐΣ, Jmm, Jab, zeros(T, length(gp.Zₐ), length(gp.Zₐ))
            ) - hyperparameter_KL_gradient(Jmm, A)
        end, # Function gradient given inducing points locations
        function ()
            return -(pr_cov(gp) \ (μ₀ - mean(gp)))
        end,
    ) # Function gradient given mean prior
end

## Gradient with respect to hyperparameter for analytical VI ##
function hyperparameter_expec_gradient(
    gp::AbstractLatent,
    ∇E_μ::AbstractVector{<:Real},
    ∇E_Σ::AbstractVector{<:Real},
    i::AnalyticVI,
    opt::AVIOptimizer,
    κΣ::AbstractMatrix{<:Real},
    Jmm::AbstractMatrix{<:Real},
    Jnm::AbstractMatrix{<:Real},
    Jnn::AbstractVector{<:Real},
)
    ι = (Jnm - gp.κ * Jmm) / pr_cov(gp)
    J̃ = Jnn - (diag_ABt(ι, gp.Knm) + diag_ABt(gp.κ, Jnm))
    dμ = dot(∇E_μ, ι * mean(gp))
    dΣ = -dot(∇E_Σ, J̃)
    dΣ += -dot(∇E_Σ, 2.0 * (diag_ABt(ι, κΣ)))
    dΣ += -dot(∇E_Σ, 2.0 * (ι * mean(gp)) .* (gp.κ * mean(gp)))
    return ρ(i) * (dμ + dΣ)
end

## Gradient with respect to hyperparameters for numerical VI ##
function hyperparameter_expec_gradient(
    gp::AbstractLatent,
    ∇E_μ::AbstractVector{<:Real},
    ∇E_Σ::AbstractVector{<:Real},
    i::NumericalVI,
    opt::NVIOptimizer,
    κΣ::AbstractMatrix{<:Real},
    Jmm::AbstractMatrix{<:Real},
    Jnm::AbstractMatrix{<:Real},
    Jnn::AbstractVector{<:Real},
)
    ι = (Jnm - gp.κ * Jmm) / pr_cov(gp)
    J̃ = Jnn - (diag_ABt(ι, gp.Knm) + diag_ABt(gp.κ, Jnm))
    dμ = dot(∇E_μ, ι * mean(gp))
    dΣ = dot(∇E_Σ, J̃ + 2.0 * diag_ABt(ι, κΣ))
    return ρ(i) * (dμ + dΣ)
end

function hyperparameter_online_gradient(
    gp::AbstractLatent,
    κₐΣ::Matrix{<:Real},
    Jmm::AbstractMatrix{<:Real},
    Jab::AbstractMatrix{<:Real},
    Jaa::AbstractMatrix{<:Real},
)
    ιₐ = (Jab - gp.κₐ * Jmm) / pr_cov(gp)
    trace_term =
        -sum(
            trace_ABt.(
                [gp.invDₐ],
                [
                    Jaa,
                    2 * ιₐ * transpose(κₐΣ),
                    -ιₐ * transpose(gp.Kab),
                    -gp.κₐ * transpose(Jab),
                ],
            ),
        ) / 2
    term_1 = dot(gp.prevη₁, ιₐ * mean(gp))
    term_2 = -dot(ιₐ * mean(gp), gp.invDₐ * gp.κₐ * mean(gp))
    return trace_term + term_1 + term_2
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 8183b318-5ebb-11eb-1cd8-a96e8704a378
using Pkg, DrWatson

# ╔═╡ 686dac30-5ebb-11eb-00f1-434980dba906
begin
	@quickactivate "StatsModelComparisons"
	using StatsModelComparisons
	using StanSample, StatsFuns, StatsPlots
	using DataFrames, CSV, JSON
end

# ╔═╡ 923212a8-630d-11eb-390f-75d21be80011
begin
	ProjDir = joinpath(mc_path, "..", "examples", "arsenic")
	include(joinpath(ProjDir, "cvit.jl"))
end;

# ╔═╡ ca0d916e-5ebe-11eb-3af1-0993bf5b82c8
md" ### Psis_loo."

# ╔═╡ 2ee71b16-5ebd-11eb-3a7d-4553f2f22274
md" ##### Data."

# ╔═╡ 37f06bc2-5ebd-11eb-31f6-65bfa74ea9cf
begin
	data = JSON.parsefile(joinpath(ProjDir, "wells.data.json"))
	y = Float64.(data["switched"])
	x = Float64[data["arsenic"]  data["dist"]]
	n, m = size(x)
end

# ╔═╡ 37f0b87a-5ebd-11eb-092a-3d9b9e8e427b
md" ##### Model."

# ╔═╡ 37f14af6-5ebd-11eb-0aae-d57e8b3d951a
begin
	model_str = read(open(joinpath(ProjDir, "arsenic_logistic.stan")), String)
	tmpdir = joinpath(ProjDir, "tmp")
	sm1 = SampleModel("arsenic_logistic", model_str; tmpdir)
	data1 = (p = m, N = n, y = Int.(y), x = x)
	rc1 = stan_sample(sm1; data=data1)
end

# ╔═╡ 38128be4-5ebd-11eb-36da-a96a68ca1316
md" ##### Fit the model in Stan."

# ╔═╡ 38135bc8-5ebd-11eb-0644-511bb08b3e59
if success(rc1)
	nt1 = read_samples(sm1)

	# Compute LOO and standard error

	log_lik1 = nt1.log_lik'
	loo1, loos1, pk1 = psisloo(log_lik1)
	elpd_loo1 = sum(loos1)
	se_elpd_loo1 = std(loos1) * sqrt(n)
	(elpd_loo = elpd_loo1, se_elpd_loo = se_elpd_loo1)
end

# ╔═╡ 44457a20-5ebd-11eb-00a1-855bfbb45e22
if success(rc1)

	# Check the shape parameter k of the generalized Pareto distribution

	pk_qualify(pk1)
end

# ╔═╡ 5a76e4aa-5ebd-11eb-2e15-6d5808ead825
md" ##### Fit a second model, using log(arsenic) instead of arsenic."

# ╔═╡ 5a771bfa-5ebd-11eb-1431-a5bc45055a5a
begin
	x2 = Float64[log.(data["arsenic"])  data["dist"]]
	data2 = (p = m, N = n, y = Int.(y), x = x2)
	rc2 = stan_sample(sm1; data=data2)
end;

# ╔═╡ 5a77d188-5ebd-11eb-0090-8d80d495f1b1
if success(rc2)
    nt2 = read_samples(sm1)

    # Compute LOO and standard error

    log_lik2 = nt2.log_lik'
    loo2, loos2, pk2 = psisloo(log_lik2)
    elpd_loo2 = sum(loos2)
    se_elpd_loo2 = std(loos2) * sqrt(n)
    (elpd_loo = elpd_loo2, se_elpd_loo = se_elpd_loo2)
end

# ╔═╡ b88657e2-5ebd-11eb-343e-9d666eef86cc
if success(rc2)

	# Check the shape parameter k of the generalized Pareto distribution

	pk_qualify(pk2)
end

# ╔═╡ dd28d430-5ebd-11eb-1854-ab4c13e82c34
if success(rc1) && success(rc2)

    ## Compare the models

    loodiff = loos1 - loos2
    (elpd_diff = sum(loodiff), se_elpd_diff = std(loodiff) * sqrt(n))
end

# ╔═╡ 98d7208e-5ebd-11eb-3a64-e568c1ceebc0
md" ### k-fold-CV."

# ╔═╡ e86ee35c-5ebd-11eb-1932-ffbde638f11e
md" ##### k-fold-CV should be used if several pk>0.5, in this case it is not needed, but provided as an example."

# ╔═╡ e86f0ab2-5ebd-11eb-1a0b-6958db3a7435
begin
	model_str_2 = read(open(joinpath(ProjDir, "arsenic_logistic_t.stan")), String)
	sm3 = SampleModel("arsenic_logistic_t", model_str_2)
end;

# ╔═╡ d83916ba-6309-11eb-309b-474efacd9f85
md" ##### Be patient... 10 cross Validations."

# ╔═╡ e70f73aa-5eb8-11eb-1960-bf6731681898
begin
	cvitr, cvitst = cvit(n, 10, true)
	kfcvs = similar(loos1)
	for cvi in 1:10
		standatacv = (p = m, N = length(cvitr[cvi]), Nt = length(cvitst[cvi]),
			x = x[cvitr[cvi],:], y = Int.(y[cvitr[cvi]]),
			xt = x[cvitst[cvi],:], yt = Int.(y[cvitst[cvi]]))

		# Fit the model in Stan
		rc3 = stan_sample(sm3; data=standatacv)
		if success(rc3)
			nt3 = read_samples(sm3)
			# Compute LOO and standard error
			log_likt = nt3.log_likt'
        	local n_sam, n_obs = size(log_likt)
        	kfcvs[cvitst[cvi]] .=
            	reshape(logsumexp(log_likt .- log(n_sam), dims=1), n_obs)
		end
	end
end

# ╔═╡ 1bf8d5a2-5ebe-11eb-12b3-2501e2f9d288
begin

	# compare PSIS-LOO and k-fold-CV

	plot([-3.5, 0], [-3.5, 0], color=:red)
	scatter!(loos1, kfcvs, xlab = "PSIS-LOO", ylab = "10-fold-CV",
		leg=false, color=:darkblue)
end

# ╔═╡ Cell order:
# ╠═8183b318-5ebb-11eb-1cd8-a96e8704a378
# ╠═686dac30-5ebb-11eb-00f1-434980dba906
# ╠═923212a8-630d-11eb-390f-75d21be80011
# ╟─ca0d916e-5ebe-11eb-3af1-0993bf5b82c8
# ╟─2ee71b16-5ebd-11eb-3a7d-4553f2f22274
# ╠═37f06bc2-5ebd-11eb-31f6-65bfa74ea9cf
# ╟─37f0b87a-5ebd-11eb-092a-3d9b9e8e427b
# ╠═37f14af6-5ebd-11eb-0aae-d57e8b3d951a
# ╟─38128be4-5ebd-11eb-36da-a96a68ca1316
# ╠═38135bc8-5ebd-11eb-0644-511bb08b3e59
# ╠═44457a20-5ebd-11eb-00a1-855bfbb45e22
# ╟─5a76e4aa-5ebd-11eb-2e15-6d5808ead825
# ╠═5a771bfa-5ebd-11eb-1431-a5bc45055a5a
# ╠═5a77d188-5ebd-11eb-0090-8d80d495f1b1
# ╠═b88657e2-5ebd-11eb-343e-9d666eef86cc
# ╠═dd28d430-5ebd-11eb-1854-ab4c13e82c34
# ╟─98d7208e-5ebd-11eb-3a64-e568c1ceebc0
# ╟─e86ee35c-5ebd-11eb-1932-ffbde638f11e
# ╠═e86f0ab2-5ebd-11eb-1a0b-6958db3a7435
# ╟─d83916ba-6309-11eb-309b-474efacd9f85
# ╠═e70f73aa-5eb8-11eb-1960-bf6731681898
# ╠═1bf8d5a2-5ebe-11eb-12b3-2501e2f9d288


text:
using ProbNumDiffEq
using Test
using LinearAlgebra
using FiniteDiff, ForwardDiff, Zygote, ReverseDiff

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems;
importodeproblems();
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_fitzhughnagumo

prob = ProbNumDiffEq.remake_prob_with_jac(prob_ode_fitzhughnagumo)

function param_to_loss(p)
    sol = solve(remake(prob, p=p), EK1(order=3), sensealg=SensitivityADPassThrough())
    return norm(sol.u[end])  # Dummy loss
end
function startval_to_loss(u0)
    sol = solve(remake(prob, u0=u0), EK1(order=3), sensealg=SensitivityADPassThrough())
    return norm(sol.u[end])  # Dummy loss
end

dldp = FiniteDiff.finite_difference_gradient(param_to_loss, prob.p)
dldu0 = FiniteDiff.finite_difference_gradient(startval_to_loss, prob.u0)

@testset "ForwardDiff.jl" begin
    @test ForwardDiff.gradient(param_to_loss, prob.p) ≈ dldp rtol = 1e-4
    @test ForwardDiff.gradient(startval_to_loss, prob.u0) ≈ dldu0 rtol = 5e-4
end

@testset "ReverseDiff.jl" begin
    @test_broken ReverseDiff.gradient(param_to_loss, prob.p) ≈ dldp
    @test_broken ReverseDiff.gradient(startval_to_loss, prob.u0) ≈ dldu0
end

@testset "Zygote.jl" begin
    @test_broken Zygote.gradient(param_to_loss, prob.p) ≈ dldp
    @test_broken Zygote.gradient(startval_to_loss, prob.u0) ≈ dldu0
end


text:
export
    TraceMetricExtractor,

    get_score,
    reset!,
    extract!,

    extract_log_likelihood,
    extract_sum_square_jerk,

    RootWeightedSquareError,
    SumSquareJerk,
    EmergentKLDivergence

include("metrics.jl")