text:
"""
```
abstract type AbstractProgression end
struct LinearProgression <: AbstractProgression end
struct DoubleProgression <: AbstractProgression end
struct PeriodProgression <: AbstractProgression end
struct BlockProgression <: AbstractProgression end
```
Progression types. Use these to make specialised functions and add new functionality.
"""
abstract type AbstractProgression end
struct LinearProgression <: AbstractProgression end
struct DoubleProgression <: AbstractProgression end
struct PeriodProgression <: AbstractProgression end
struct BlockProgression <: AbstractProgression end

"""
```
abstract type AbstractProgramme end
```
Programme types.
"""
abstract type AbstractProgramme end

import Base: length, round, getindex, iterate, @_inline_meta
Base.round(x::Real, y::Real, mode::Function = floor) = mode(x / y) * y
Base.length(x::Function) = 1

"""
```
calcIntensity(reps::Integer, rpe::Real)
```
This function calculates a set's intensity as a function of reps and RPE. It does a disconcertedly good job of reproducing RPE charts. It is defined as follows:

``z = \\dfrac{1}{a + b (x + 10 - y) + (x - 1) \\left(\\dfrac{c}{x} + \\dfrac{d}{y}\\right)},``

where ``z \\equiv`` intensity, ``x \\equiv`` reps, ``y \\equiv`` RPE, ``a = 0.995,~ b = 0.0333,~c = 0.0025,~d = 0.1``.

We cap `rpe` to 10 and `intensity` to 1.
"""
function calcIntensity(reps::Integer, rpe::Real)
    a = 0.995
    b = 0.0333
    c = 0.0025
    d = 0.1
    rpe = minimum((rpe, 10))
    return minimum((1 / (a + b * (reps + 10 - rpe) + (reps - 1) * (c / reps + d / rpe)), 1))
end

"""
```
calcRPE(reps::Integer, intensity::Real)
```
Solving [`calcIntensity`](@ref) for `rpe` yields the following function:

``x = \\dfrac{\\sqrt{\\left(a y z - b y^2 z + 10 b y z + c y z - d z - y\\right)^2 + 4 c y z (b y z + d z)} - a y z + b y^2 z - 10 b y z - c y z + d z + y}{2 z (b y + d)},``

where the variables are the same as [`calcIntensity`](@ref).

We cap `intensity` to 1 and the `rpe` to 10.
"""
function calcRPE(reps::Integer, intensity::Real)
    a = 0.995
    b = 0.0333
    c = 0.0025
    d = 0.1
    intensity = minimum((intensity, 1))
    rpe =
        (
            sqrt(
                (
                    a * reps * intensity +
                    b * reps^2 * intensity +
                    10 * b * reps * intensity +
                    c * reps * intensity - c * intensity - reps
                )^2 +
                4 * b * reps * intensity * (d * reps^2 * intensity - d * reps * intensity),
            ) +
            a * reps * intensity +
            b * reps^2 * intensity +
            10 * b * reps * intensity +
            c * reps * intensity - c * intensity - reps
        ) / (2 * b * reps * intensity)

    return minimum((rpe, 10))
end

"""
```
calcReps(intensity::Real, rpe::Real)
```
Solving [`calcIntensity`](@ref) for `reps` yields the following function:

``y = \\dfrac{\\sqrt{\\left(a x z + b x^2 z + 10 b x z + c x z - c z - x\\right)^2 + 4 b x z \\left(d x^2 z - d x z\\right)} + a x z + b x^2 z + 10 b x z + c x z - c z - x}{2 b x z},``

where the variables are the same as [`calcIntensity`](@ref).

We cap `intensity` to 1 and `rpe` to 10.
"""
function calcReps(intensity::Real, rpe::Real)
    a = 0.995
    b = 0.0333
    c = 0.0025
    d = 0.1
    intensity = minimum((intensity, 1))
    rpe = minimum((rpe, 10))
    reps =
        (
            sqrt(
                (
                    a * rpe * intensity - b * rpe^2 * intensity +
                    10 * b * rpe * intensity +
                    c * rpe * intensity - d * intensity - rpe
                )^2 + 4 * c * rpe * intensity * (b * rpe * intensity + d * intensity),
            ) - a * rpe * intensity + b * rpe^2 * intensity - 10 * b * rpe * intensity -
            c * rpe * intensity +
            d * intensity +
            rpe
        ) / (2 * intensity * (b * rpe + d))

    return Int(round(reps, digits = 0))
end

"""
```
calcIntensityRatio(
    actualReps::Integer,
    actualRPE::Real,
    targetReps::Integer,
    targetRPE::Real,
)
```
Calculate the ratio between `targetIntensity/actualIntensity`.

All RPE values are capped to 10.
"""
function calcIntensityRatio(
    actualReps::Integer,
    actualRPE::Real,
    targetReps::Integer,
    targetRPE::Real,
)
    actualRPE = minimum((actualRPE, 10))
    targetRPE = minimum((targetRPE, 10))

    actualIntensity = calcIntensity(actualReps, actualRPE)
    targetIntensity = calcIntensity(targetReps, targetRPE)
    return targetIntensity / actualIntensity
end

"""
```
calcRPERatio(
    actualReps::Integer,
    actualIntensity::Real,
    targetReps::Integer,
    targetIntensity::Real,
)
```
Calculate the ratio between `targetRPE/actualRPE`.

All intensity values are capped to 1.
"""
function calcRPERatio(
    actualReps::Integer,
    actualIntensity::Real,
    targetReps::Integer,
    targetIntensity::Real,
)
    actualIntensity = minimum((actualIntensity, 1))
    targetIntensity = minimum((targetIntensity, 1))

    actualRPE = calcRPE(actualReps, actualIntensity)
    targetRPE = calcRPE(targetReps, targetIntensity)
    return targetRPE / actualRPE
end

"""
```
calcRepRatio(
    actualIntensity::Real,
    actualRPE::Real,
    targetIntensity::Real,
    targetRPE::Real,
)
```
Calculate the ratio between `targetReps/actualReps`.

All intensity values are capped to 1 and all rpe values to 10.
"""
function calcRepRatio(
    actualIntensity::Real,
    actualRPE::Real,
    targetIntensity::Real,
    targetRPE::Real,
)
    actualIntensity = minimum((actualIntensity, 1.0))
    targetIntensity = minimum((targetIntensity, 1.0))
    actualRPE = minimum((actualRPE, 10.0))
    targetRPE = minimum((targetRPE, 10.0))

    actualReps = calcReps(actualIntensity, actualRPE)
    targetReps = calcReps(targetIntensity, targetRPE)
    return targetReps / actualReps
end

"""
```
calcRepMax(
    weight::Real,
    actualReps::Integer,
    actualRPE::Real,
    targetReps::Integer,
    targetRPE::Real,
)
```
Calculates the rep weight for a target number of reps, `targetReps`, at a target rpe, `targetRPE`, given an actual number of reps, `actualReps`, and an actual rpe, `actualRPE` and `weight`. It's just the weight multiplied by the intensity ratio.
"""
function calcRepMax(
    weight::Real,
    actualReps::Integer,
    actualRPE::Real,
    targetReps::Integer,
    targetRPE::Real,
)

    actualRPE = minimum((actualRPE, 10.0))
    targetRPE = minimum((targetRPE, 10.0))

    intensity = calcIntensityRatio(actualReps, actualRPE, targetReps, targetRPE)
    return weight * intensity
end

"""
```
intensityArb(var::Integer)
```
Calculates intensity given an arbitrary variable. I've seen this used as a proxy for reps at a given RPE. I've seen `var` be anything between the number of reps to the number of reps + 2, 4, 6, 8 depending on a target rpe. This makes for a good rough guide. It is defined as

``z = \\dfrac{1}{a + b x},``

where ``x \\equiv`` var and the constants are the same as [`calcIntensity`](@ref).

!!! note
    `calcIntensity` works over a *much* wider range of RPE and rep combinations.
"""
function intensityArb(var::Integer)
    return 1 / (0.995 + 0.0333 * var)
end

"""
```
mutable struct SetScheme{
    T1 <: Union{<:AbstractString, Vector{<:AbstractString}},
    T2 <: Union{<:Integer, Vector{<:Integer}},
    T3 <: Union{<:Real, Vector{<:Real}},
    T4 <: Union{<:Function, Vector{<:Function}},
    T5 <: Bool,
}
    type::T1
    sets::T2
    reps::T2
    intensity::T3
    rpe::T3
    addWeight::T3
    roundMode::T4
    wght::T3
    rpeMode::T5
end
```
This structure contains user-created set schemes. Its parameters are as follows:

- `type` defines the set type,
- `sets` the number of sets of the corresponding type,
- `reps` is the number of reps in the set,
- `intensity` is the set intensity,
- `rpe` is the set's target RPE,
- `addWeight` is extra weight on top of the normally calculated weight (useful for chains, bands or microplates),
- `wght` is the set weight, which is calculated later when pairing progression schemes to exercises,
- `rpeMode` whether the set scheme is defined via RPE or percentage intensity, if true, the values in `intensity` will be assumed to be RPE and the percentage intensity will be calculated using RPE and the number of reps in the set, else the values in `intensity` will be used as provided.

It uses the following constructor, which takes care of any required computations.
```
SetScheme(;
    type::T1 = "Default",
    sets::T2 = 5,
    reps::T2 = 5,
    intensity::T3 = 0.75,
    addWeight::T4 = 0,
    roundMode::T5 = floor,
    rpeMode::Bool = false,
) where {
    T1 <: Union{<:AbstractString, Vector{<:AbstractString}},
    T2 <: Union{<:Integer, Vector{<:Integer}},
    T3 <: Union{<:Real, Vector{<:Real}},
    T4 <: Union{<:Real, Vector{<:Real}},
    T5 <: Union{<:Function, Vector{<:Function}},
}
```
Every parameter has defaults, so users have the ability to provide only the ones they want to modify.

If `rpeMode == true`, intensity is based on RPE, which lies in the interval ``x ∈ (0, 10]``. If `rpeMode == false`, then intensity is relative to training max and therefore lies in the interval ``x ∈ [0, 1]``. The lenghts of all input vectors must be equal.

## Example
Here we define a simple set scheme. Since we are not providing any added weight via `addWeight`, that parameter becomes a zero vector.
```
julia> SampleScheme = SetScheme(;
           type = ["Long Rest", "Longer Rest", "Longest Rest", "Optional Forced Reps"],
           sets = [1, 2, 1, 1],
           reps = [12, 14, 10, 5],
           intensity = [9.5, 10, 10, 10],
           roundMode = [floor, floor, ceil, ceil],
           rpeMode = true,
       )
SetScheme{Array{String,1},Array{Int64,1},Array{Float64,1},Array{Function,1},Bool}(["Long Rest", "Longer Rest", "Longest Rest", "Optional Forced Reps"], [1, 2, 1, 1], [12, 14, 10, 5], [0.6538806237677003, 0.6275409806672554, 0.7041013906002465, 0.8309098462816784], [9.5, 10.0, 10.0, 10.0], [0.0, 0.0, 0.0, 0.0], Function[floor, floor, ceil, ceil], [0.0, 0.0, 0.0, 0.0], true)
```
"""
mutable struct SetScheme{
    T1<:Union{<:AbstractString,Vector{<:AbstractString}},
    T2<:Union{<:Integer,Vector{<:Integer}},
    T3<:Union{<:Real,Vector{<:Real}},
    T4<:Union{<:Function,Vector{<:Function}},
    T5<:Bool,
}
    type::T1
    sets::T2
    reps::T2
    intensity::T3
    rpe::T3
    addWeight::T3
    roundMode::T4
    wght::T3
    rpeMode::T5

    function SetScheme(;
        type::T1 = "Default",
        sets::T2 = 5,
        reps::T2 = 5,
        intensity::T3 = 0.75,
        addWeight::T4 = 0,
        roundMode::T5 = floor,
        rpeMode::Bool = false,
    ) where {
        T1<:Union{<:AbstractString,Vector{<:AbstractString}},
        T2<:Union{<:Integer,Vector{<:Integer}},
        T3<:Union{<:Real,Vector{<:Real}},
        T4<:Union{<:Real,Vector{<:Real}},
        T5<:Union{<:Function,Vector{<:Function}},
    }
        difSets = length(sets)
        difSets > 1 && length(intensity) == 1 ? intensity = fill(0.75, difSets) : nothing
        difSets > 1 && length(addWeight) == 1 ? addWeight = zeros(difSets) : nothing
        difSets > 1 && length(roundMode) == 1 ? roundMode = fill(floor, difSets) : nothing

        @assert length(sets) ==
                length(reps) ==
                length(intensity) ==
                length(addWeight) ==
                length(roundMode) "lengths of sets $(length(sets)), reps $(length(reps)), intensity $(length(intensity)), addWeight $(length(addWeight)) and roundMode $(length(roundMode)) must be equal."
        rpe = 0.0
        wght = 0.0
        if difSets > 1
            rpe = zeros(difSets)
            wght = zeros(difSets)
        end
        if rpeMode
            rpe = intensity
            intensity = calcIntensity.(reps, intensity)
        else
            rpe = calcRPE.(reps, intensity)
        end
        return new{
            typeof(type),
            typeof(sets),
            typeof(intensity),
            typeof(roundMode),
            typeof(rpeMode),
        }(
            type,
            sets,
            reps,
            intensity,
            rpe,
            convert.(eltype(intensity), addWeight),
            roundMode,
            wght,
            rpeMode,
        )
    end
end
getindex(x::SetScheme, i::Integer) = i == 1 ? x : throw(BoundsError)
length(x::SetScheme) = 1
iterate(A::SetScheme, i = 1) =
    (@_inline_meta; (i % UInt) - 1 < length(A) ? (A[1], i + 1) : nothing)

"""
```
struct Progression{
    T1 <: AbstractProgression,
    T2 <: AbstractString,
    T3 <: Integer,
    T4 <: Union{<:SetScheme, Vector{<:SetScheme}},
} <: AbstractProgression
    type::T1
    name::T2
    sessions::T3
    period::T3
    setScheme::T4
end
```
This structure holds user-defined progressions. Its parameters are as follows:
- `type` is the progression type which is a subtype of [`AbstractProgression`](@ref),
- `name` is the progression name,
- `sessions` is the number of different sessions in a week,
- `period` is the progression period,
- `setScheme` is a previously defined set scheme or vector of set schemes, the length must be equal to `sessions * period`.

This structure lets users create different types of progressions, from simple day in day out progressions to complex undulating periodisation and block progressions.

Progressions are created with the following constructor function,
```
function Progression(;
    type::T1,
    name::T2,
    sessions::T3 = 1,
    period::T3 = 1,
    setScheme::T4,
) where {
    T1 <: AbstractProgression,
    T2 <: AbstractString,
    T3 <: Integer,
    T4 <: Union{<:SetScheme, Vector{<:SetScheme}},
}
```

## Example

Assuming we've run the example in [`SetScheme`](@ref) we can create a simple progression.
```
julia> SampleProgression = Progression(;
           type = LinearProgression(),
           name = "Progression Name",
           setScheme = SampleScheme,
       )
Progression{LinearProgression,String,Int64,SetScheme{Array{String,1},Array{Int64,1},Array{Float64,1},Array{Function,1},Bool}}(LinearProgression(), "Progression Name", 1, 1, SetScheme{Array{String,1},Array{Int64,1},Array{Float64,1},Array{Function,1},Bool}(["Long Rest", "Longer Rest", "Longest Rest", "Optional Forced Reps"], [1, 2, 1, 1], [12, 14, 10, 5], [0.6538806237677003, 0.6275409806672554, 0.7041013906002465, 0.8309098462816784], [9.5, 10.0, 10.0, 10.0], [0.0, 0.0, 0.0, 0.0], Function[floor, floor, ceil, ceil], [0.0, 0.0, 0.0, 0.0], true))
```
"""
struct Progression{
    T1<:AbstractProgression,
    T2<:AbstractString,
    T3<:Integer,
    T4<:Union{<:SetScheme,Vector{<:SetScheme}},
} <: AbstractProgression
    type::T1
    name::T2
    sessions::T3
    period::T3
    setScheme::T4

    function Progression(;
        type::T1,
        name::T2,
        sessions::T3 = 1,
        period::T3 = 1,
        setScheme::T4,
    ) where {
        T1<:AbstractProgression,
        T2<:AbstractString,
        T3<:Integer,
        T4<:Union{<:SetScheme,Vector{<:SetScheme}},
    }
        @assert length(setScheme) == sessions * period "length of setScheme, $(length(setScheme)), must be equal to sessions * period, $(sessions*period)."

        return new{typeof(type),typeof(name),typeof(sessions),typeof(setScheme)}(
            type,
            name,
            sessions,
            period,
            setScheme,
        )
    end
end
getindex(x::Progression, i::Integer) = i == 1 ? x : throw(BoundsError)
length(x::Progression) = 1
iterate(A::Progression, i = 1) =
    (@_inline_meta; (i % UInt) - 1 < length(A) ? (A[1], i + 1) : nothing)

"""
```
mutable struct Exercise{
    T1 <: AbstractString,
    T2 <: Union{AbstractString, Vector{<:AbstractString}},
    T3 <: Union{AbstractString, Vector{<:AbstractString}},
    T4 <: Union{AbstractString, Vector{<:AbstractString}},
    T5 <: Union{AbstractString, Vector{<:AbstractString}},
    T6 <: Real,
    T7 <: Real,
    T8 <: Function,
}
    name::T1
    equipment::T2
    modality::T3
    size::T4
    muscles::T5
    trainingMax::T6
    roundBase::T7
    roundMode::T8
end
```
User created exercises are encapsulated in this structure. Its parameters are as follows:
- `name` is the exercise name,
- `equipment` is the equipment required,
- `modality` are any modifiers to the exercise (for example blocks),
- `size` are the size/resistance of said modifiers (for example 2 inch),
- `muscles` are the muscles targeted,
- `trainingMax` is the training maximum used to calculate set weights when pairing with a progression,
- `roundBase` is the base to which set weights are rounded,
- `roundMode` is the rounding function used.

This uses the following constructor which provides default values for anything users don't need/want to track,
```
function Exercise(;
    name::T1,
    equipment::T2 = "NA",
    modality::T3 = "Default",
    muscles::T4 = "NA",
    trainingMax::T6 = 0,
    size::T5 = "NA",
    roundBase::T7 = 2.5,
    roundMode::T8 = floor,
) where {
    T1 <: AbstractString,
    T2 <: Union{AbstractString, Vector{<:AbstractString}},
    T3 <: Union{AbstractString, Vector{<:AbstractString}},
    T4 <: Union{AbstractString, Vector{<:AbstractString}},
    T5 <: Union{AbstractString, Vector{<:AbstractString}},
    T6 <: Real,
    T7 <: Real,
    T8 <: Function,
}
end
```

## Examples

Say we want to bench press off a 1 inch block and a semi supinated grip on a swiss bar. We can easily create this exercise
```
julia> BenchAccessory = Exercise(;
           name = "Swiss Bar Bench",
           modality = ["Block" , "Semi-Supinated"],
           equipment = "Swiss Bar",
           size = "1 inch",
           muscles = ["Pecs", "Triceps", "Front Delts"],
           trainingMax = 100,
           roundBase = 2.5,
           roundMode = floor,
       )
Exercise{String,String,Array{String,1},String,Array{String,1},Float64,Float64,typeof(floor)}("Swiss Bar Bench", "Swiss Bar", ["Block", "Semi-Supinated"], "1 inch", ["Pecs", "Triceps", "Front Delts"], 100.0, 2.5, floor)
```
"""
mutable struct Exercise{
    T1<:AbstractString,
    T2<:Union{AbstractString,Vector{<:AbstractString}},
    T3<:Union{AbstractString,Vector{<:AbstractString}},
    T4<:Union{AbstractString,Vector{<:AbstractString}},
    T5<:Union{AbstractString,Vector{<:AbstractString}},
    T6<:Real,
    T7<:Real,
    T8<:Function,
}
    name::T1
    equipment::T2
    modality::T3
    size::T4
    muscles::T5
    trainingMax::T6
    roundBase::T7
    roundMode::T8

    function Exercise(;
        name::T1,
        equipment::T2 = "NA",
        modality::T3 = "Default",
        muscles::T4 = "NA",
        trainingMax::T6 = 0,
        size::T5 = "NA",
        roundBase::T7 = 2.5,
        roundMode::T8 = floor,
    ) where {
        T1<:AbstractString,
        T2<:Union{AbstractString,Vector{<:AbstractString}},
        T3<:Union{AbstractString,Vector{<:AbstractString}},
        T4<:Union{AbstractString,Vector{<:AbstractString}},
        T5<:Union{AbstractString,Vector{<:AbstractString}},
        T6<:Real,
        T7<:Real,
        T8<:Function,
    }
        trainingMaxRound = round(trainingMax, roundBase, roundMode)

        return new{
            typeof(name),
            typeof(equipment),
            typeof(modality),
            typeof(size),
            typeof(muscles),
            typeof(trainingMaxRound),
            typeof(roundBase),
            typeof(roundMode),
        }(
            name,
            equipment,
            modality,
            size,
            muscles,
            trainingMaxRound,
            roundBase,
            roundMode,
        )
    end
end
getindex(x::Exercise, i::Integer) = i == 1 ? x : throw(BoundsError)
length(x::Exercise) = 1
iterate(A::Exercise, i = 1) =
    (@_inline_meta; (i % UInt) - 1 < length(A) ? (A[1], i + 1) : nothing)

"""
```
calcWeights(exercise::Exercise, setScheme::SetScheme)
```
This function calculates and updates the weights for each set according to the training max specified in `exercise.trainingMax` and the intensities of the sets in `setScheme`.
"""
function calcWeights(exercise::Exercise, setScheme::SetScheme)
    trainingMax = exercise.trainingMax
    roundBase = exercise.roundBase
    reps = setScheme.reps
    intensity = setScheme.intensity
    addWeight = setScheme.addWeight
    roundMode = setScheme.roundMode

    # Calculate wieghts.
    setScheme.wght = round.(trainingMax * intensity + addWeight, roundBase, roundMode)

    # Calculate target minimum RPE for a set.
    intense = setScheme.wght / trainingMax
    setScheme.rpe = round.(calcRPE.(reps, intense), digits = 2)

    return setScheme
end
function calcWeights(exercise::Exercise, prog::Progression)
    return calcWeights.(exercise, prog.setScheme)
end

"""
```
struct Programme{
    T1 <: AbstractProgramme,
    T2 <: AbstractString,
    T3 <: Dict{Any, Any},
    T4 <: Vector{Any},
}
    type::T1
    name::T2
    exerProg::T3
    days::T4
end
```
User created programmes are made with this structure. The parameters are as follows:
- `type` is the exercise type and is a subtype of [`AbstractProgramme`](@ref) (we recommend users define custom concrete types for their programmes),
- `name` is the programme name (we recommend the name is the same as the programme type),
- `exerProg` is a dictionary which pairs exercises and progressions with the key-value pair `name => (exercise, progression)`.
- `days` is the programme schedule where each entry in the vector has all the exercises for the day.

Each programme is different so the recommendation is to create a typed [`makeDays`](@ref) function for the programme type.

!!! warning
    [`makeDays`](@ref) must be explicitly imported either with `import Lifting: makeDays` or `Lifting.makeDays` for it to be automatically used in [`updateMaxes!`](@ref).

## Examples

Assuming we are using the previously defined [`Exercise`](@ref) and [`Progression`](@ref), we can create a programme. There are a few ways to do so but this is recommended.
```julia
# Define the programme's concrete type. Can be used to extend `makeDays`.
struct SampleProgramme <: AbstractProgramme end

# Define a dictionary that will contain an Exercise and its corresponding
# progression as a tuple or vector.
exerProg = Dict()

# Push exercises and progressions to dictionary. Using named tuples is
# recommended because it enables dot syntax as well as regular indexing.
# Each exercise is paired with its progression.
push!(
    exerProg,
    "BenchAccessory" => (exercise = BenchAccessory, progression = SampleProgression),
)

# Define the a function that creates the schedule. This lets you automatically
# update training maxes and makes changing the programme easier. Here we also
# define rest days but you can leave them out if need be. `makeDays` will be
# called for whatever type is given as its first argument. Different programme
# types will necessarily have and call different `makeDays` functions.
import Lifting: makeDays
function makeDays(::SampleProgramme, exerProg)

    # One sub array per day you want to track.
    week = [[], [], [], [], [], [], []]

    # You can 'unroll' the dictionary contents or push them directly to your
    # schedule. If using named tuples dot syntax and normal numeric indexing
    # both work
    # exerProg["BenchAccessory"].exercise == exerProg["BenchAccessory"][1]
    BenchAccessory = exerProg["BenchAccessory"].exercise
    BenchAccessoryProg = exerProg["BenchAccessory"].progression

    # Push the exercise and progression to whichever day of the week you want.
    # If the progression has multiple set schemes, the corresponding scheme can
    # be accessed by adding the index as the last argument.
    # push!(week[1], BenchAccessory, BenchAccessoryProg, 2)
    # Which would fail because there is only one set scheme, but if
    # BenchAccessoryProg were to have more entries that would take the second.
    # push! is overloaded so it also automatically calculates the required
    # weights. You can do this manually but this is easier and more reliable.
    push!(week[1], BenchAccessory, BenchAccessoryProg)   # Day 1
    push!(week[2], "Rest")  # Day 2
    push!(week[3], BenchAccessory, BenchAccessoryProg)   # Day 3
    push!(week[4], "Rest")  # Day 4
    push!(week[5], BenchAccessory, BenchAccessoryProg)   # Day 5
    push!(week[5], "Rest")  # Day 6
    push!(week[7], "Rest")  # Day 7

    return week
end

# Create weekly schedule.
week = makeDays(SampleProgramme(), exerProg)

# Create programme.
sampleProgramme = Programme(
                            SampleProgramme(),
                            "SampleProgramme",
                            exerProg,
                            week
                        )
```

It's recommended users either import the default dictionary [`Lifting_Programmes`](@ref) and store their programmes there with key equal to the programme name and type. They can also create their own dictionary for their own creations.
"""
struct Programme{T1<:AbstractProgramme,T2<:AbstractString,T3<:Dict{Any,Any},T4<:Any}

    type::T1
    name::T2
    exerProg::T3
    days::T4

    function Programme(
        type::T1,
        name::T2,
        exerProg::T3,
        days::T4,
    ) where {T1<:AbstractProgramme,T2<:AbstractString,T3<:Dict{Any,Any},T4<:Any}

        return new{typeof(type),typeof(name),typeof(exerProg),typeof(days)}(
            type,
            name,
            exerProg,
            days,
        )
    end
end

getindex(p::Programme, idx...) = [p.days[i] for i in idx]
getindex(p::Programme, i) = p.days[i]
length(p::Programme) = length(p.days)

import Base: push!, insert!
"""
```
push!(
    A::AbstractArray{T, 1} where {T},
    exercise::Exercise,
    progression::Progression,
    i::Integer = 1,
)
```
Calculates set weights using `exercise.trainingMax` and `progression.setScheme`, then unrolls the relevant data into a named tuple and pushes it to the end of vector A. See [`insert!`](@ref) which does the same but at a specified index of A.
"""
function push!(
    A::AbstractArray{T,1} where {T},
    exercise::Exercise,
    progression::Progression,
    i::Integer = 1,
)
    calcWeights.(exercise, progression.setScheme[i])
    if typeof(progression.setScheme[i].type) == String
        type = progression.setScheme[i].type
    else
        type = Tuple(progression.setScheme[i].type)
    end
    return if exercise.modality == "Default"
        push!(
            A,
            (
                name = exercise.name,
                type = type,
                sets = Tuple(progression.setScheme[i].sets),
                reps = Tuple(progression.setScheme[i].reps),
                wght = Tuple(progression.setScheme[i].wght),
                rpe = Tuple(progression.setScheme[i].rpe),
                intensity = Tuple(progression.setScheme[i].intensity),
            ),
        )
    else
        push!(
            A,
            (
                name = exercise.name,
                modality = exercise.modality,
                type = type,
                sets = Tuple(progression.setScheme[i].sets),
                reps = Tuple(progression.setScheme[i].reps),
                wght = Tuple(progression.setScheme[i].wght),
                rpe = Tuple(progression.setScheme[i].rpe),
                intensity = Tuple(progression.setScheme[i].intensity),
            ),
        )
    end
end

"""
```
insert!(
    A::AbstractArray{T, 1} where {T},
    index::Integer,
    exercise::Exercise,
    progression::Progression,
    i::Integer = 1,
)
```
Calculates set weights using `exercise.trainingMax` and `progression.setScheme`, then unrolls the relevant data into a named tuple and pushes it to index `index` of vector A. See [`push!`](@ref) which does the same but in the last index of A.
"""
function insert!(
    A::AbstractArray{T,1} where {T},
    index::Integer,
    exercise::Exercise,
    progression::Progression,
    i::Integer = 1,
)
    calcWeights.(exercise, progression.setScheme[i])
    if typeof(progression.setScheme[i].type) == String
        type = progression.setScheme[i].type
    else
        type = Tuple(progression.setScheme[i].type)
    end
    return if exercise.modality == "Default"
        insert!(
            A,
            index,
            (
                name = exercise.name,
                type = type,
                sets = Tuple(progression.setScheme[i].sets),
                reps = Tuple(progression.setScheme[i].reps),
                wght = Tuple(progression.setScheme[i].wght),
                rpe = Tuple(progression.setScheme[i].rpe),
                intensity = Tuple(progression.setScheme[i].intensity),
            ),
        )
    else
        insert!(
            A,
            index,
            (
                name = exercise.name,
                modality = exercise.modality,
                type = type,
                sets = Tuple(progression.setScheme[i].sets),
                reps = Tuple(progression.setScheme[i].reps),
                wght = Tuple(progression.setScheme[i].wght),
                rpe = Tuple(progression.setScheme[i].rpe),
                intensity = Tuple(progression.setScheme[i].intensity),
            ),
        )
    end
end

"""
```
adjustMaxes!(
    name::AbstractString,
    dict::Dict{Any, Any},
    actualReps::Integer;
    weight = missing,
)
```
Adjusts a programme's training maxes using what is expected by the programme and the actual reps, `actualReps`, performed at the programme's highest intensity for the exercise contained in `dict[name]`. However, if `weight` is provided, `updateMaxes!` will use that to calculate the new training max instead of whatever the programme prescribes.

!!! note
    Adjustments made to training max are calculated based on what the programme expects vs actual performance. So if you spectacularly fail a set, the adjustment *down* will be relatively aggressive. Conversely if you absolutely smash it, the adjustment *up* will be relatively aggressive. Downward adjustments are more aggressive than upward ones.
"""
function adjustMaxes!(
    name::AbstractString,
    dict::Dict{Any,Any},
    actualReps::Integer;
    weight = missing,
)
    entry = dict[name]
    exercise = entry[1]
    prog = entry[2]
    setScheme = prog.setScheme
    numSets = length(setScheme)

    idx1 = 0
    idx2 = 0
    tmp = 0.0
    old = 0.0
    maxWght = 0
    for i = 1:numSets
        old = tmp
        wght = setScheme[i].wght
        tmp = maximum(wght)
        tmp = maximum((tmp, old))
        if tmp > old || maxWght == 0
            maxWght = tmp
            idx1 = i
            idx3 = findfirst(x -> x == maxWght, wght)
            !isnothing(idx3) ? idx2 = idx3 : idx2 = idx1
        end
    end
    targetReps = setScheme[idx1].reps[idx2]
    targetIntensity = setScheme[idx1].intensity[idx2]
    targetRPE = calcRPE(targetReps, targetIntensity)
    actualRPE = calcRPE(actualReps, targetIntensity)
    metTarget = targetRPE <= actualRPE
    actualRPE = 2 * targetRPE - actualRPE

    # If we gave a specific weight we want to use it instead of using the maximum weight.
    if ismissing(weight)
        trainingMax = calcRepMax(maxWght, actualReps, actualRPE, targetReps, targetRPE)
    else
        trainingMax = calcRepMax(weight, actualReps, actualRPE, targetReps, targetRPE)
    end

    roundBase = exercise.roundBase
    roundMode = exercise.roundMode
    change = 0
    if metTarget
        if trainingMax < exercise.trainingMax
            change = roundBase
        else
            change = round(trainingMax - exercise.trainingMax, roundBase, roundMode)
            change = maximum((change, roundBase))
        end
    else
        if trainingMax < exercise.trainingMax
            change = round(exercise.trainingMax - trainingMax, roundBase, roundMode)
            change = -maximum((change, roundBase))
        else
            change = round(trainingMax - exercise.trainingMax, roundBase, roundMode)
            change = maximum((change, roundBase))
        end
    end

    exercise.trainingMax += change
    return calcWeights(exercise, prog)
end

"""
```
adjustMaxes(
    name::AbstractString,
    dict::Dict{Any, Any},
    actualReps::Integer;
    weight = missing,
)
```
Does the same as [`adjustMaxes!`](@ref) but without updating the programme's training maxes. It returns a tuple with the new training max and the would-be change to the programme's training max, `(trainingMax, change)``.
"""
function adjustMaxes(
    name::AbstractString,
    dict::Dict{Any,Any},
    actualReps::Integer;
    weight = missing,
)
    entry = dict[name]
    exercise = entry[1]
    prog = entry[2]
    setScheme = prog.setScheme
    numSets = length(setScheme)

    idx1 = 0
    idx2 = 0
    tmp = 0.0
    old = 0.0
    maxWght = 0
    for i = 1:numSets
        old = tmp
        wght = setScheme[i].wght
        tmp = maximum(wght)
        tmp = maximum((tmp, old))
        if tmp > old || maxWght == 0
            maxWght = tmp
            idx1 = i
            idx3 = findfirst(x -> x == maxWght, wght)
            !isnothing(idx3) ? idx2 = idx3 : idx2 = idx1
        end
    end
    targetReps = setScheme[idx1].reps[idx2]
    targetIntensity = setScheme[idx1].intensity[idx2]
    targetRPE = calcRPE(targetReps, targetIntensity)
    actualRPE = calcRPE(actualReps, targetIntensity)
    metTarget = targetRPE <= actualRPE
    actualRPE = 2 * targetRPE - actualRPE
    # If we gave a specific weight we want to use it instead of using the maximum weight.
    if ismissing(weight)
        trainingMax = calcRepMax(maxWght, actualReps, actualRPE, targetReps, targetRPE)
    else
        trainingMax = calcRepMax(weight, actualReps, actualRPE, targetReps, targetRPE)
    end

    roundBase = exercise.roundBase
    roundMode = exercise.roundMode
    change = 0
    if metTarget
        if trainingMax < exercise.trainingMax
            change = roundBase
        else
            change = round(trainingMax - exercise.trainingMax, roundBase, roundMode)
            change = maximum((change, roundBase))
        end
    else
        if trainingMax < exercise.trainingMax
            change = round(exercise.trainingMax - trainingMax, roundBase, roundMode)
            change = -maximum((change, roundBase))
        else
            change = round(trainingMax - exercise.trainingMax, roundBase, roundMode)
            change = maximum((change, roundBase))
        end
    end

    return trainingMax, change
end

"""
```
makeDays()
```
Generic function that can be explicitly typed by the user to create their programme's daily schedule. Typed `makeDays` functions are used by [`updateMaxes!`](@ref) to update the training days with the new training maxes calculated by [`adjustMaxes!`](@ref). See the example in [`Programme`](@ref) to see how to define typed `makeDays`.

!!! warning
    [`makeDays`](@ref) must be explicitly imported either with `import Lifting: makeDays` or `Lifting.makeDays` for it to be automatically used in [`updateMaxes!`](@ref).
"""
function makeDays() end

"""
```
updateMaxes!(prog::Programme, names, reps; idx = missing)
```
`names` are the keys of the exercises you want to update, `reps` contains the reps performed at the corresponding exercise's maximum intensity set of the cycle, `idx` is optional in case each entry of reps is a vector and you would like to specifically use an index instead of the last entry.

Updates the programme's training maxes and adjusts the set weights for the next training cycle. It uses [`adjustMaxes`](@ref) to update training maxes and [`makeDays`](@ref) to calculate new set weights according to what is expected by the programme vs actual performance. All this is done on a per exercise basis.

See [`updateMaxes`](@ref) for a non-mutating version of this function.
"""
function updateMaxes!(prog::Programme, names, reps; idx = missing)
    typeof(names) == String ? names = [names] : nothing

    exerProg = prog.exerProg
    for (i, name) in enumerate(names)
        isempty(reps[name]) ? continue : nothing
        ismissing(idx) ? numReps = reps[name][end] : numReps = reps[name][idx]
        numReps < 0 ? continue : nothing
        adjustMaxes!(name, exerProg, numReps)
    end

    prog.days .= makeDays(prog.type, exerProg)
    return prog
end

"""
```
updateMaxes(prog::Programme, names, reps; idx = missing)
```
The arguments are the same as [`updateMaxes!`](@ref).

Calculates the would-be new training maxes and changes to the old training maxes for the next cycle according to what is expected by the programme vs actual performance but does **not** update the programme in any way. Returns a tuple `(trainingMaxes, change)` where each entry in the tuple is an array whose entries correspond 1 to 1 to the exercises in `names`. The new value for `trainingMax` is not necessarily going to equal to the old value of `trainingMax + change`, because `change` is calculated based on how far off the new value for `trainingMax` is from the old value of `trainingMax`, and the result is rounded according to the exercise's `roundBase` and `roundMode`.

This is essentially calls [`adjustMaxes`](@ref) for all the arguments in `names`. In order to see what a programme's training days would look like under these new training maxes, the easiest thing to do is to create a deepcopy() of the programme and run [`updateMaxes!`](@ref) on it.
"""
function updateMaxes(prog::Programme, names, reps; idx = missing)
    typeof(names) == String ? names = [names] : nothing

    trainingMaxes = zeros(length(names))
    change = zeros(length(names))
    exerProg = prog.exerProg
    for (i, name) in enumerate(names)
        isempty(reps[name]) ? continue : nothing
        ismissing(idx) ? numReps = reps[name][end] : numReps = reps[name][idx]
        numReps < 0 ? continue : nothing
        trainingMaxes[i], change[i] = adjustMaxes(name, exerProg, numReps)
    end
    return trainingMaxes, change

end

"""
```
calcTrainingMaxLogs(prog::Programme, names, reps, weight)
```
This function does something similar to [`updateMaxes`](@ref) but for the whole history of entries instead of a single point in time. This is useful for analysing progress.
"""
function calcTrainingMaxLogs(prog::Programme, names, reps, weight)
    typeof(names) == String ? names = [names] : nothing

    trainingMax = deepcopy(weight)
    change = deepcopy(weight)
    for name in names
        isempty(reps[name]) ? continue : nothing
        for j = 1:length(reps[name])
            change[name][j] = 0.0
            trainingMax[name][j] = 0.0
            trainingMax[name][j], change[name][j] =
                adjustMaxes(name, prog.exerProg, reps[name][j]; weight = weight[name][j])
        end
    end
    return trainingMax, change
end

max_stars_repo_name: dcelisgarza/Lifting


text:
# --------------------------------------------------------------------------
# --------------------------------------------------------------------------
# The following gradients were taken from https://github.com/BIDS/colormap/blob/master/colormaps.py
# Here is the licensing note which accompanied this:

    # New matplotlib colormaps by Nathaniel J. Smith, Stefan van der Walt,
    # and (in the case of viridis) Eric Firing.
    #
    # This file and the colormaps in it are released under the CC0 license /
    # public domain dedication. We would appreciate credit if you use or
    # redistribute these colormaps, but do not impose any legal restrictions.
    #
    # To the extent possible under law, the persons who associated CC0 with
    # mpl-colormaps have waived all copyright and related or neighboring rights
    # to mpl-colormaps.
    #
    # You should have received a copy of the CC0 legalcode along with this
    # work.  If not, see <http://creativecommons.org/publicdomain/zero/1.0/>.

function sample_evenly(v::AVec, n::Integer = length(v))
    idx = Int[round(Int, x) for x in linspace(1, length(v), n)]
    v[idx]
end

# note: to use the full arrays, just take out the second param (n) from the sample_evenly calls

register_gradient_colors(:magma, sample_evenly([
    RGB(0.001462, 0.000466, 0.013866),
    RGB(0.002258, 0.001295, 0.018331),
    RGB(0.003279, 0.002305, 0.023708),
    RGB(0.004512, 0.003490, 0.029965),
    RGB(0.005950, 0.004843, 0.037130),
    RGB(0.007588, 0.006356, 0.044973),
    RGB(0.009426, 0.008022, 0.052844),
    RGB(0.011465, 0.009828, 0.060750),
    RGB(0.013708, 0.011771, 0.068667),
    RGB(0.016156, 0.013840, 0.076603),
    RGB(0.018815, 0.016026, 0.084584),
    RGB(0.021692, 0.018320, 0.092610),
    RGB(0.024792, 0.020715, 0.100676),
    RGB(0.028123, 0.023201, 0.108787),
    RGB(0.031696, 0.025765, 0.116965),
    RGB(0.035520, 0.028397, 0.125209),
    RGB(0.039608, 0.031090, 0.133515),
    RGB(0.043830, 0.033830, 0.141886),
    RGB(0.048062, 0.036607, 0.150327),
    RGB(0.052320, 0.039407, 0.158841),
    RGB(0.056615, 0.042160, 0.167446),
    RGB(0.060949, 0.044794, 0.176129),
    RGB(0.065330, 0.047318, 0.184892),
    RGB(0.069764, 0.049726, 0.193735),
    RGB(0.074257, 0.052017, 0.202660),
    RGB(0.078815, 0.054184, 0.211667),
    RGB(0.083446, 0.056225, 0.220755),
    RGB(0.088155, 0.058133, 0.229922),
    RGB(0.092949, 0.059904, 0.239164),
    RGB(0.097833, 0.061531, 0.248477),
    RGB(0.102815, 0.063010, 0.257854),
    RGB(0.107899, 0.064335, 0.267289),
    RGB(0.113094, 0.065492, 0.276784),
    RGB(0.118405, 0.066479, 0.286321),
    RGB(0.123833, 0.067295, 0.295879),
    RGB(0.129380, 0.067935, 0.305443),
    RGB(0.135053, 0.068391, 0.315000),
    RGB(0.140858, 0.068654, 0.324538),
    RGB(0.146785, 0.068738, 0.334011),
    RGB(0.152839, 0.068637, 0.343404),
    RGB(0.159018, 0.068354, 0.352688),
    RGB(0.165308, 0.067911, 0.361816),
    RGB(0.171713, 0.067305, 0.370771),
    RGB(0.178212, 0.066576, 0.379497),
    RGB(0.184801, 0.065732, 0.387973),
    RGB(0.191460, 0.064818, 0.396152),
    RGB(0.198177, 0.063862, 0.404009),
    RGB(0.204935, 0.062907, 0.411514),
    RGB(0.211718, 0.061992, 0.418647),
    RGB(0.218512, 0.061158, 0.425392),
    RGB(0.225302, 0.060445, 0.431742),
    RGB(0.232077, 0.059889, 0.437695),
    RGB(0.238826, 0.059517, 0.443256),
    RGB(0.245543, 0.059352, 0.448436),
    RGB(0.252220, 0.059415, 0.453248),
    RGB(0.258857, 0.059706, 0.457710),
    RGB(0.265447, 0.060237, 0.461840),
    RGB(0.271994, 0.060994, 0.465660),
    RGB(0.278493, 0.061978, 0.469190),
    RGB(0.284951, 0.063168, 0.472451),
    RGB(0.291366, 0.064553, 0.475462),
    RGB(0.297740, 0.066117, 0.478243),
    RGB(0.304081, 0.067835, 0.480812),
    RGB(0.310382, 0.069702, 0.483186),
    RGB(0.316654, 0.071690, 0.485380),
    RGB(0.322899, 0.073782, 0.487408),
    RGB(0.329114, 0.075972, 0.489287),
    RGB(0.335308, 0.078236, 0.491024),
    RGB(0.341482, 0.080564, 0.492631),
    RGB(0.347636, 0.082946, 0.494121),
    RGB(0.353773, 0.085373, 0.495501),
    RGB(0.359898, 0.087831, 0.496778),
    RGB(0.366012, 0.090314, 0.497960),
    RGB(0.372116, 0.092816, 0.499053),
    RGB(0.378211, 0.095332, 0.500067),
    RGB(0.384299, 0.097855, 0.501002),
    RGB(0.390384, 0.100379, 0.501864),
    RGB(0.396467, 0.102902, 0.502658),
    RGB(0.402548, 0.105420, 0.503386),
    RGB(0.408629, 0.107930, 0.504052),
    RGB(0.414709, 0.110431, 0.504662),
    RGB(0.420791, 0.112920, 0.505215),
    RGB(0.426877, 0.115395, 0.505714),
    RGB(0.432967, 0.117855, 0.506160),
    RGB(0.439062, 0.120298, 0.506555),
    RGB(0.445163, 0.122724, 0.506901),
    RGB(0.451271, 0.125132, 0.507198),
    RGB(0.457386, 0.127522, 0.507448),
    RGB(0.463508, 0.129893, 0.507652),
    RGB(0.469640, 0.132245, 0.507809),
    RGB(0.475780, 0.134577, 0.507921),
    RGB(0.481929, 0.136891, 0.507989),
    RGB(0.488088, 0.139186, 0.508011),
    RGB(0.494258, 0.141462, 0.507988),
    RGB(0.500438, 0.143719, 0.507920),
    RGB(0.506629, 0.145958, 0.507806),
    RGB(0.512831, 0.148179, 0.507648),
    RGB(0.519045, 0.150383, 0.507443),
    RGB(0.525270, 0.152569, 0.507192),
    RGB(0.531507, 0.154739, 0.506895),
    RGB(0.537755, 0.156894, 0.506551),
    RGB(0.544015, 0.159033, 0.506159),
    RGB(0.550287, 0.161158, 0.505719),
    RGB(0.556571, 0.163269, 0.505230),
    RGB(0.562866, 0.165368, 0.504692),
    RGB(0.569172, 0.167454, 0.504105),
    RGB(0.575490, 0.169530, 0.503466),
    RGB(0.581819, 0.171596, 0.502777),
    RGB(0.588158, 0.173652, 0.502035),
    RGB(0.594508, 0.175701, 0.501241),
    RGB(0.600868, 0.177743, 0.500394),
    RGB(0.607238, 0.179779, 0.499492),
    RGB(0.613617, 0.181811, 0.498536),
    RGB(0.620005, 0.183840, 0.497524),
    RGB(0.626401, 0.185867, 0.496456),
    RGB(0.632805, 0.187893, 0.495332),
    RGB(0.639216, 0.189921, 0.494150),
    RGB(0.645633, 0.191952, 0.492910),
    RGB(0.652056, 0.193986, 0.491611),
    RGB(0.658483, 0.196027, 0.490253),
    RGB(0.664915, 0.198075, 0.488836),
    RGB(0.671349, 0.200133, 0.487358),
    RGB(0.677786, 0.202203, 0.485819),
    RGB(0.684224, 0.204286, 0.484219),
    RGB(0.690661, 0.206384, 0.482558),
    RGB(0.697098, 0.208501, 0.480835),
    RGB(0.703532, 0.210638, 0.479049),
    RGB(0.709962, 0.212797, 0.477201),
    RGB(0.716387, 0.214982, 0.475290),
    RGB(0.722805, 0.217194, 0.473316),
    RGB(0.729216, 0.219437, 0.471279),
    RGB(0.735616, 0.221713, 0.469180),
    RGB(0.742004, 0.224025, 0.467018),
    RGB(0.748378, 0.226377, 0.464794),
    RGB(0.754737, 0.228772, 0.462509),
    RGB(0.761077, 0.231214, 0.460162),
    RGB(0.767398, 0.233705, 0.457755),
    RGB(0.773695, 0.236249, 0.455289),
    RGB(0.779968, 0.238851, 0.452765),
    RGB(0.786212, 0.241514, 0.450184),
    RGB(0.792427, 0.244242, 0.447543),
    RGB(0.798608, 0.247040, 0.444848),
    RGB(0.804752, 0.249911, 0.442102),
    RGB(0.810855, 0.252861, 0.439305),
    RGB(0.816914, 0.255895, 0.436461),
    RGB(0.822926, 0.259016, 0.433573),
    RGB(0.828886, 0.262229, 0.430644),
    RGB(0.834791, 0.265540, 0.427671),
    RGB(0.840636, 0.268953, 0.424666),
    RGB(0.846416, 0.272473, 0.421631),
    RGB(0.852126, 0.276106, 0.418573),
    RGB(0.857763, 0.279857, 0.415496),
    RGB(0.863320, 0.283729, 0.412403),
    RGB(0.868793, 0.287728, 0.409303),
    RGB(0.874176, 0.291859, 0.406205),
    RGB(0.879464, 0.296125, 0.403118),
    RGB(0.884651, 0.300530, 0.400047),
    RGB(0.889731, 0.305079, 0.397002),
    RGB(0.894700, 0.309773, 0.393995),
    RGB(0.899552, 0.314616, 0.391037),
    RGB(0.904281, 0.319610, 0.388137),
    RGB(0.908884, 0.324755, 0.385308),
    RGB(0.913354, 0.330052, 0.382563),
    RGB(0.917689, 0.335500, 0.379915),
    RGB(0.921884, 0.341098, 0.377376),
    RGB(0.925937, 0.346844, 0.374959),
    RGB(0.929845, 0.352734, 0.372677),
    RGB(0.933606, 0.358764, 0.370541),
    RGB(0.937221, 0.364929, 0.368567),
    RGB(0.940687, 0.371224, 0.366762),
    RGB(0.944006, 0.377643, 0.365136),
    RGB(0.947180, 0.384178, 0.363701),
    RGB(0.950210, 0.390820, 0.362468),
    RGB(0.953099, 0.397563, 0.361438),
    RGB(0.955849, 0.404400, 0.360619),
    RGB(0.958464, 0.411324, 0.360014),
    RGB(0.960949, 0.418323, 0.359630),
    RGB(0.963310, 0.425390, 0.359469),
    RGB(0.965549, 0.432519, 0.359529),
    RGB(0.967671, 0.439703, 0.359810),
    RGB(0.969680, 0.446936, 0.360311),
    RGB(0.971582, 0.454210, 0.361030),
    RGB(0.973381, 0.461520, 0.361965),
    RGB(0.975082, 0.468861, 0.363111),
    RGB(0.976690, 0.476226, 0.364466),
    RGB(0.978210, 0.483612, 0.366025),
    RGB(0.979645, 0.491014, 0.367783),
    RGB(0.981000, 0.498428, 0.369734),
    RGB(0.982279, 0.505851, 0.371874),
    RGB(0.983485, 0.513280, 0.374198),
    RGB(0.984622, 0.520713, 0.376698),
    RGB(0.985693, 0.528148, 0.379371),
    RGB(0.986700, 0.535582, 0.382210),
    RGB(0.987646, 0.543015, 0.385210),
    RGB(0.988533, 0.550446, 0.388365),
    RGB(0.989363, 0.557873, 0.391671),
    RGB(0.990138, 0.565296, 0.395122),
    RGB(0.990871, 0.572706, 0.398714),
    RGB(0.991558, 0.580107, 0.402441),
    RGB(0.992196, 0.587502, 0.406299),
    RGB(0.992785, 0.594891, 0.410283),
    RGB(0.993326, 0.602275, 0.414390),
    RGB(0.993834, 0.609644, 0.418613),
    RGB(0.994309, 0.616999, 0.422950),
    RGB(0.994738, 0.624350, 0.427397),
    RGB(0.995122, 0.631696, 0.431951),
    RGB(0.995480, 0.639027, 0.436607),
    RGB(0.995810, 0.646344, 0.441361),
    RGB(0.996096, 0.653659, 0.446213),
    RGB(0.996341, 0.660969, 0.451160),
    RGB(0.996580, 0.668256, 0.456192),
    RGB(0.996775, 0.675541, 0.461314),
    RGB(0.996925, 0.682828, 0.466526),
    RGB(0.997077, 0.690088, 0.471811),
    RGB(0.997186, 0.697349, 0.477182),
    RGB(0.997254, 0.704611, 0.482635),
    RGB(0.997325, 0.711848, 0.488154),
    RGB(0.997351, 0.719089, 0.493755),
    RGB(0.997351, 0.726324, 0.499428),
    RGB(0.997341, 0.733545, 0.505167),
    RGB(0.997285, 0.740772, 0.510983),
    RGB(0.997228, 0.747981, 0.516859),
    RGB(0.997138, 0.755190, 0.522806),
    RGB(0.997019, 0.762398, 0.528821),
    RGB(0.996898, 0.769591, 0.534892),
    RGB(0.996727, 0.776795, 0.541039),
    RGB(0.996571, 0.783977, 0.547233),
    RGB(0.996369, 0.791167, 0.553499),
    RGB(0.996162, 0.798348, 0.559820),
    RGB(0.995932, 0.805527, 0.566202),
    RGB(0.995680, 0.812706, 0.572645),
    RGB(0.995424, 0.819875, 0.579140),
    RGB(0.995131, 0.827052, 0.585701),
    RGB(0.994851, 0.834213, 0.592307),
    RGB(0.994524, 0.841387, 0.598983),
    RGB(0.994222, 0.848540, 0.605696),
    RGB(0.993866, 0.855711, 0.612482),
    RGB(0.993545, 0.862859, 0.619299),
    RGB(0.993170, 0.870024, 0.626189),
    RGB(0.992831, 0.877168, 0.633109),
    RGB(0.992440, 0.884330, 0.640099),
    RGB(0.992089, 0.891470, 0.647116),
    RGB(0.991688, 0.898627, 0.654202),
    RGB(0.991332, 0.905763, 0.661309),
    RGB(0.990930, 0.912915, 0.668481),
    RGB(0.990570, 0.920049, 0.675675),
    RGB(0.990175, 0.927196, 0.682926),
    RGB(0.989815, 0.934329, 0.690198),
    RGB(0.989434, 0.941470, 0.697519),
    RGB(0.989077, 0.948604, 0.704863),
    RGB(0.988717, 0.955742, 0.712242),
    RGB(0.988367, 0.962878, 0.719649),
    RGB(0.988033, 0.970012, 0.727077),
    RGB(0.987691, 0.977154, 0.734536),
    RGB(0.987387, 0.984288, 0.742002),
    RGB(0.987053, 0.991438, 0.749504)
], 30))

register_gradient_colors(:inferno, sample_evenly([
    RGB(0.001462, 0.000466, 0.013866),
    RGB(0.002267, 0.001270, 0.018570),
    RGB(0.003299, 0.002249, 0.024239),
    RGB(0.004547, 0.003392, 0.030909),
    RGB(0.006006, 0.004692, 0.038558),
    RGB(0.007676, 0.006136, 0.046836),
    RGB(0.009561, 0.007713, 0.055143),
    RGB(0.011663, 0.009417, 0.063460),
    RGB(0.013995, 0.011225, 0.071862),
    RGB(0.016561, 0.013136, 0.080282),
    RGB(0.019373, 0.015133, 0.088767),
    RGB(0.022447, 0.017199, 0.097327),
    RGB(0.025793, 0.019331, 0.105930),
    RGB(0.029432, 0.021503, 0.114621),
    RGB(0.033385, 0.023702, 0.123397),
    RGB(0.037668, 0.025921, 0.132232),
    RGB(0.042253, 0.028139, 0.141141),
    RGB(0.046915, 0.030324, 0.150164),
    RGB(0.051644, 0.032474, 0.159254),
    RGB(0.056449, 0.034569, 0.168414),
    RGB(0.061340, 0.036590, 0.177642),
    RGB(0.066331, 0.038504, 0.186962),
    RGB(0.071429, 0.040294, 0.196354),
    RGB(0.076637, 0.041905, 0.205799),
    RGB(0.081962, 0.043328, 0.215289),
    RGB(0.087411, 0.044556, 0.224813),
    RGB(0.092990, 0.045583, 0.234358),
    RGB(0.098702, 0.046402, 0.243904),
    RGB(0.104551, 0.047008, 0.253430),
    RGB(0.110536, 0.047399, 0.262912),
    RGB(0.116656, 0.047574, 0.272321),
    RGB(0.122908, 0.047536, 0.281624),
    RGB(0.129285, 0.047293, 0.290788),
    RGB(0.135778, 0.046856, 0.299776),
    RGB(0.142378, 0.046242, 0.308553),
    RGB(0.149073, 0.045468, 0.317085),
    RGB(0.155850, 0.044559, 0.325338),
    RGB(0.162689, 0.043554, 0.333277),
    RGB(0.169575, 0.042489, 0.340874),
    RGB(0.176493, 0.041402, 0.348111),
    RGB(0.183429, 0.040329, 0.354971),
    RGB(0.190367, 0.039309, 0.361447),
    RGB(0.197297, 0.038400, 0.367535),
    RGB(0.204209, 0.037632, 0.373238),
    RGB(0.211095, 0.037030, 0.378563),
    RGB(0.217949, 0.036615, 0.383522),
    RGB(0.224763, 0.036405, 0.388129),
    RGB(0.231538, 0.036405, 0.392400),
    RGB(0.238273, 0.036621, 0.396353),
    RGB(0.244967, 0.037055, 0.400007),
    RGB(0.251620, 0.037705, 0.403378),
    RGB(0.258234, 0.038571, 0.406485),
    RGB(0.264810, 0.039647, 0.409345),
    RGB(0.271347, 0.040922, 0.411976),
    RGB(0.277850, 0.042353, 0.414392),
    RGB(0.284321, 0.043933, 0.416608),
    RGB(0.290763, 0.045644, 0.418637),
    RGB(0.297178, 0.047470, 0.420491),
    RGB(0.303568, 0.049396, 0.422182),
    RGB(0.309935, 0.051407, 0.423721),
    RGB(0.316282, 0.053490, 0.425116),
    RGB(0.322610, 0.055634, 0.426377),
    RGB(0.328921, 0.057827, 0.427511),
    RGB(0.335217, 0.060060, 0.428524),
    RGB(0.341500, 0.062325, 0.429425),
    RGB(0.347771, 0.064616, 0.430217),
    RGB(0.354032, 0.066925, 0.430906),
    RGB(0.360284, 0.069247, 0.431497),
    RGB(0.366529, 0.071579, 0.431994),
    RGB(0.372768, 0.073915, 0.432400),
    RGB(0.379001, 0.076253, 0.432719),
    RGB(0.385228, 0.078591, 0.432955),
    RGB(0.391453, 0.080927, 0.433109),
    RGB(0.397674, 0.083257, 0.433183),
    RGB(0.403894, 0.085580, 0.433179),
    RGB(0.410113, 0.087896, 0.433098),
    RGB(0.416331, 0.090203, 0.432943),
    RGB(0.422549, 0.092501, 0.432714),
    RGB(0.428768, 0.094790, 0.432412),
    RGB(0.434987, 0.097069, 0.432039),
    RGB(0.441207, 0.099338, 0.431594),
    RGB(0.447428, 0.101597, 0.431080),
    RGB(0.453651, 0.103848, 0.430498),
    RGB(0.459875, 0.106089, 0.429846),
    RGB(0.466100, 0.108322, 0.429125),
    RGB(0.472328, 0.110547, 0.428334),
    RGB(0.478558, 0.112764, 0.427475),
    RGB(0.484789, 0.114974, 0.426548),
    RGB(0.491022, 0.117179, 0.425552),
    RGB(0.497257, 0.119379, 0.424488),
    RGB(0.503493, 0.121575, 0.423356),
    RGB(0.509730, 0.123769, 0.422156),
    RGB(0.515967, 0.125960, 0.420887),
    RGB(0.522206, 0.128150, 0.419549),
    RGB(0.528444, 0.130341, 0.418142),
    RGB(0.534683, 0.132534, 0.416667),
    RGB(0.540920, 0.134729, 0.415123),
    RGB(0.547157, 0.136929, 0.413511),
    RGB(0.553392, 0.139134, 0.411829),
    RGB(0.559624, 0.141346, 0.410078),
    RGB(0.565854, 0.143567, 0.408258),
    RGB(0.572081, 0.145797, 0.406369),
    RGB(0.578304, 0.148039, 0.404411),
    RGB(0.584521, 0.150294, 0.402385),
    RGB(0.590734, 0.152563, 0.400290),
    RGB(0.596940, 0.154848, 0.398125),
    RGB(0.603139, 0.157151, 0.395891),
    RGB(0.609330, 0.159474, 0.393589),
    RGB(0.615513, 0.161817, 0.391219),
    RGB(0.621685, 0.164184, 0.388781),
    RGB(0.627847, 0.166575, 0.386276),
    RGB(0.633998, 0.168992, 0.383704),
    RGB(0.640135, 0.171438, 0.381065),
    RGB(0.646260, 0.173914, 0.378359),
    RGB(0.652369, 0.176421, 0.375586),
    RGB(0.658463, 0.178962, 0.372748),
    RGB(0.664540, 0.181539, 0.369846),
    RGB(0.670599, 0.184153, 0.366879),
    RGB(0.676638, 0.186807, 0.363849),
    RGB(0.682656, 0.189501, 0.360757),
    RGB(0.688653, 0.192239, 0.357603),
    RGB(0.694627, 0.195021, 0.354388),
    RGB(0.700576, 0.197851, 0.351113),
    RGB(0.706500, 0.200728, 0.347777),
    RGB(0.712396, 0.203656, 0.344383),
    RGB(0.718264, 0.206636, 0.340931),
    RGB(0.724103, 0.209670, 0.337424),
    RGB(0.729909, 0.212759, 0.333861),
    RGB(0.735683, 0.215906, 0.330245),
    RGB(0.741423, 0.219112, 0.326576),
    RGB(0.747127, 0.222378, 0.322856),
    RGB(0.752794, 0.225706, 0.319085),
    RGB(0.758422, 0.229097, 0.315266),
    RGB(0.764010, 0.232554, 0.311399),
    RGB(0.769556, 0.236077, 0.307485),
    RGB(0.775059, 0.239667, 0.303526),
    RGB(0.780517, 0.243327, 0.299523),
    RGB(0.785929, 0.247056, 0.295477),
    RGB(0.791293, 0.250856, 0.291390),
    RGB(0.796607, 0.254728, 0.287264),
    RGB(0.801871, 0.258674, 0.283099),
    RGB(0.807082, 0.262692, 0.278898),
    RGB(0.812239, 0.266786, 0.274661),
    RGB(0.817341, 0.270954, 0.270390),
    RGB(0.822386, 0.275197, 0.266085),
    RGB(0.827372, 0.279517, 0.261750),
    RGB(0.832299, 0.283913, 0.257383),
    RGB(0.837165, 0.288385, 0.252988),
    RGB(0.841969, 0.292933, 0.248564),
    RGB(0.846709, 0.297559, 0.244113),
    RGB(0.851384, 0.302260, 0.239636),
    RGB(0.855992, 0.307038, 0.235133),
    RGB(0.860533, 0.311892, 0.230606),
    RGB(0.865006, 0.316822, 0.226055),
    RGB(0.869409, 0.321827, 0.221482),
    RGB(0.873741, 0.326906, 0.216886),
    RGB(0.878001, 0.332060, 0.212268),
    RGB(0.882188, 0.337287, 0.207628),
    RGB(0.886302, 0.342586, 0.202968),
    RGB(0.890341, 0.347957, 0.198286),
    RGB(0.894305, 0.353399, 0.193584),
    RGB(0.898192, 0.358911, 0.188860),
    RGB(0.902003, 0.364492, 0.184116),
    RGB(0.905735, 0.370140, 0.179350),
    RGB(0.909390, 0.375856, 0.174563),
    RGB(0.912966, 0.381636, 0.169755),
    RGB(0.916462, 0.387481, 0.164924),
    RGB(0.919879, 0.393389, 0.160070),
    RGB(0.923215, 0.399359, 0.155193),
    RGB(0.926470, 0.405389, 0.150292),
    RGB(0.929644, 0.411479, 0.145367),
    RGB(0.932737, 0.417627, 0.140417),
    RGB(0.935747, 0.423831, 0.135440),
    RGB(0.938675, 0.430091, 0.130438),
    RGB(0.941521, 0.436405, 0.125409),
    RGB(0.944285, 0.442772, 0.120354),
    RGB(0.946965, 0.449191, 0.115272),
    RGB(0.949562, 0.455660, 0.110164),
    RGB(0.952075, 0.462178, 0.105031),
    RGB(0.954506, 0.468744, 0.099874),
    RGB(0.956852, 0.475356, 0.094695),
    RGB(0.959114, 0.482014, 0.089499),
    RGB(0.961293, 0.488716, 0.084289),
    RGB(0.963387, 0.495462, 0.079073),
    RGB(0.965397, 0.502249, 0.073859),
    RGB(0.967322, 0.509078, 0.068659),
    RGB(0.969163, 0.515946, 0.063488),
    RGB(0.970919, 0.522853, 0.058367),
    RGB(0.972590, 0.529798, 0.053324),
    RGB(0.974176, 0.536780, 0.048392),
    RGB(0.975677, 0.543798, 0.043618),
    RGB(0.977092, 0.550850, 0.039050),
    RGB(0.978422, 0.557937, 0.034931),
    RGB(0.979666, 0.565057, 0.031409),
    RGB(0.980824, 0.572209, 0.028508),
    RGB(0.981895, 0.579392, 0.026250),
    RGB(0.982881, 0.586606, 0.024661),
    RGB(0.983779, 0.593849, 0.023770),
    RGB(0.984591, 0.601122, 0.023606),
    RGB(0.985315, 0.608422, 0.024202),
    RGB(0.985952, 0.615750, 0.025592),
    RGB(0.986502, 0.623105, 0.027814),
    RGB(0.986964, 0.630485, 0.030908),
    RGB(0.987337, 0.637890, 0.034916),
    RGB(0.987622, 0.645320, 0.039886),
    RGB(0.987819, 0.652773, 0.045581),
    RGB(0.987926, 0.660250, 0.051750),
    RGB(0.987945, 0.667748, 0.058329),
    RGB(0.987874, 0.675267, 0.065257),
    RGB(0.987714, 0.682807, 0.072489),
    RGB(0.987464, 0.690366, 0.079990),
    RGB(0.987124, 0.697944, 0.087731),
    RGB(0.986694, 0.705540, 0.095694),
    RGB(0.986175, 0.713153, 0.103863),
    RGB(0.985566, 0.720782, 0.112229),
    RGB(0.984865, 0.728427, 0.120785),
    RGB(0.984075, 0.736087, 0.129527),
    RGB(0.983196, 0.743758, 0.138453),
    RGB(0.982228, 0.751442, 0.147565),
    RGB(0.981173, 0.759135, 0.156863),
    RGB(0.980032, 0.766837, 0.166353),
    RGB(0.978806, 0.774545, 0.176037),
    RGB(0.977497, 0.782258, 0.185923),
    RGB(0.976108, 0.789974, 0.196018),
    RGB(0.974638, 0.797692, 0.206332),
    RGB(0.973088, 0.805409, 0.216877),
    RGB(0.971468, 0.813122, 0.227658),
    RGB(0.969783, 0.820825, 0.238686),
    RGB(0.968041, 0.828515, 0.249972),
    RGB(0.966243, 0.836191, 0.261534),
    RGB(0.964394, 0.843848, 0.273391),
    RGB(0.962517, 0.851476, 0.285546),
    RGB(0.960626, 0.859069, 0.298010),
    RGB(0.958720, 0.866624, 0.310820),
    RGB(0.956834, 0.874129, 0.323974),
    RGB(0.954997, 0.881569, 0.337475),
    RGB(0.953215, 0.888942, 0.351369),
    RGB(0.951546, 0.896226, 0.365627),
    RGB(0.950018, 0.903409, 0.380271),
    RGB(0.948683, 0.910473, 0.395289),
    RGB(0.947594, 0.917399, 0.410665),
    RGB(0.946809, 0.924168, 0.426373),
    RGB(0.946392, 0.930761, 0.442367),
    RGB(0.946403, 0.937159, 0.458592),
    RGB(0.946903, 0.943348, 0.474970),
    RGB(0.947937, 0.949318, 0.491426),
    RGB(0.949545, 0.955063, 0.507860),
    RGB(0.951740, 0.960587, 0.524203),
    RGB(0.954529, 0.965896, 0.540361),
    RGB(0.957896, 0.971003, 0.556275),
    RGB(0.961812, 0.975924, 0.571925),
    RGB(0.966249, 0.980678, 0.587206),
    RGB(0.971162, 0.985282, 0.602154),
    RGB(0.976511, 0.989753, 0.616760),
    RGB(0.982257, 0.994109, 0.631017),
    RGB(0.988362, 0.998364, 0.644924)
], 30))

register_gradient_colors(:plasma, sample_evenly([
    RGB(0.050383, 0.029803, 0.527975),
    RGB(0.063536, 0.028426, 0.533124),
    RGB(0.075353, 0.027206, 0.538007),
    RGB(0.086222, 0.026125, 0.542658),
    RGB(0.096379, 0.025165, 0.547103),
    RGB(0.105980, 0.024309, 0.551368),
    RGB(0.115124, 0.023556, 0.555468),
    RGB(0.123903, 0.022878, 0.559423),
    RGB(0.132381, 0.022258, 0.563250),
    RGB(0.140603, 0.021687, 0.566959),
    RGB(0.148607, 0.021154, 0.570562),
    RGB(0.156421, 0.020651, 0.574065),
    RGB(0.164070, 0.020171, 0.577478),
    RGB(0.171574, 0.019706, 0.580806),
    RGB(0.178950, 0.019252, 0.584054),
    RGB(0.186213, 0.018803, 0.587228),
    RGB(0.193374, 0.018354, 0.590330),
    RGB(0.200445, 0.017902, 0.593364),
    RGB(0.207435, 0.017442, 0.596333),
    RGB(0.214350, 0.016973, 0.599239),
    RGB(0.221197, 0.016497, 0.602083),
    RGB(0.227983, 0.016007, 0.604867),
    RGB(0.234715, 0.015502, 0.607592),
    RGB(0.241396, 0.014979, 0.610259),
    RGB(0.248032, 0.014439, 0.612868),
    RGB(0.254627, 0.013882, 0.615419),
    RGB(0.261183, 0.013308, 0.617911),
    RGB(0.267703, 0.012716, 0.620346),
    RGB(0.274191, 0.012109, 0.622722),
    RGB(0.280648, 0.011488, 0.625038),
    RGB(0.287076, 0.010855, 0.627295),
    RGB(0.293478, 0.010213, 0.629490),
    RGB(0.299855, 0.009561, 0.631624),
    RGB(0.306210, 0.008902, 0.633694),
    RGB(0.312543, 0.008239, 0.635700),
    RGB(0.318856, 0.007576, 0.637640),
    RGB(0.325150, 0.006915, 0.639512),
    RGB(0.331426, 0.006261, 0.641316),
    RGB(0.337683, 0.005618, 0.643049),
    RGB(0.343925, 0.004991, 0.644710),
    RGB(0.350150, 0.004382, 0.646298),
    RGB(0.356359, 0.003798, 0.647810),
    RGB(0.362553, 0.003243, 0.649245),
    RGB(0.368733, 0.002724, 0.650601),
    RGB(0.374897, 0.002245, 0.651876),
    RGB(0.381047, 0.001814, 0.653068),
    RGB(0.387183, 0.001434, 0.654177),
    RGB(0.393304, 0.001114, 0.655199),
    RGB(0.399411, 0.000859, 0.656133),
    RGB(0.405503, 0.000678, 0.656977),
    RGB(0.411580, 0.000577, 0.657730),
    RGB(0.417642, 0.000564, 0.658390),
    RGB(0.423689, 0.000646, 0.658956),
    RGB(0.429719, 0.000831, 0.659425),
    RGB(0.435734, 0.001127, 0.659797),
    RGB(0.441732, 0.001540, 0.660069),
    RGB(0.447714, 0.002080, 0.660240),
    RGB(0.453677, 0.002755, 0.660310),
    RGB(0.459623, 0.003574, 0.660277),
    RGB(0.465550, 0.004545, 0.660139),
    RGB(0.471457, 0.005678, 0.659897),
    RGB(0.477344, 0.006980, 0.659549),
    RGB(0.483210, 0.008460, 0.659095),
    RGB(0.489055, 0.010127, 0.658534),
    RGB(0.494877, 0.011990, 0.657865),
    RGB(0.500678, 0.014055, 0.657088),
    RGB(0.506454, 0.016333, 0.656202),
    RGB(0.512206, 0.018833, 0.655209),
    RGB(0.517933, 0.021563, 0.654109),
    RGB(0.523633, 0.024532, 0.652901),
    RGB(0.529306, 0.027747, 0.651586),
    RGB(0.534952, 0.031217, 0.650165),
    RGB(0.540570, 0.034950, 0.648640),
    RGB(0.546157, 0.038954, 0.647010),
    RGB(0.551715, 0.043136, 0.645277),
    RGB(0.557243, 0.047331, 0.643443),
    RGB(0.562738, 0.051545, 0.641509),
    RGB(0.568201, 0.055778, 0.639477),
    RGB(0.573632, 0.060028, 0.637349),
    RGB(0.579029, 0.064296, 0.635126),
    RGB(0.584391, 0.068579, 0.632812),
    RGB(0.589719, 0.072878, 0.630408),
    RGB(0.595011, 0.077190, 0.627917),
    RGB(0.600266, 0.081516, 0.625342),
    RGB(0.605485, 0.085854, 0.622686),
    RGB(0.610667, 0.090204, 0.619951),
    RGB(0.615812, 0.094564, 0.617140),
    RGB(0.620919, 0.098934, 0.614257),
    RGB(0.625987, 0.103312, 0.611305),
    RGB(0.631017, 0.107699, 0.608287),
    RGB(0.636008, 0.112092, 0.605205),
    RGB(0.640959, 0.116492, 0.602065),
    RGB(0.645872, 0.120898, 0.598867),
    RGB(0.650746, 0.125309, 0.595617),
    RGB(0.655580, 0.129725, 0.592317),
    RGB(0.660374, 0.134144, 0.588971),
    RGB(0.665129, 0.138566, 0.585582),
    RGB(0.669845, 0.142992, 0.582154),
    RGB(0.674522, 0.147419, 0.578688),
    RGB(0.679160, 0.151848, 0.575189),
    RGB(0.683758, 0.156278, 0.571660),
    RGB(0.688318, 0.160709, 0.568103),
    RGB(0.692840, 0.165141, 0.564522),
    RGB(0.697324, 0.169573, 0.560919),
    RGB(0.701769, 0.174005, 0.557296),
    RGB(0.706178, 0.178437, 0.553657),
    RGB(0.710549, 0.182868, 0.550004),
    RGB(0.714883, 0.187299, 0.546338),
    RGB(0.719181, 0.191729, 0.542663),
    RGB(0.723444, 0.196158, 0.538981),
    RGB(0.727670, 0.200586, 0.535293),
    RGB(0.731862, 0.205013, 0.531601),
    RGB(0.736019, 0.209439, 0.527908),
    RGB(0.740143, 0.213864, 0.524216),
    RGB(0.744232, 0.218288, 0.520524),
    RGB(0.748289, 0.222711, 0.516834),
    RGB(0.752312, 0.227133, 0.513149),
    RGB(0.756304, 0.231555, 0.509468),
    RGB(0.760264, 0.235976, 0.505794),
    RGB(0.764193, 0.240396, 0.502126),
    RGB(0.768090, 0.244817, 0.498465),
    RGB(0.771958, 0.249237, 0.494813),
    RGB(0.775796, 0.253658, 0.491171),
    RGB(0.779604, 0.258078, 0.487539),
    RGB(0.783383, 0.262500, 0.483918),
    RGB(0.787133, 0.266922, 0.480307),
    RGB(0.790855, 0.271345, 0.476706),
    RGB(0.794549, 0.275770, 0.473117),
    RGB(0.798216, 0.280197, 0.469538),
    RGB(0.801855, 0.284626, 0.465971),
    RGB(0.805467, 0.289057, 0.462415),
    RGB(0.809052, 0.293491, 0.458870),
    RGB(0.812612, 0.297928, 0.455338),
    RGB(0.816144, 0.302368, 0.451816),
    RGB(0.819651, 0.306812, 0.448306),
    RGB(0.823132, 0.311261, 0.444806),
    RGB(0.826588, 0.315714, 0.441316),
    RGB(0.830018, 0.320172, 0.437836),
    RGB(0.833422, 0.324635, 0.434366),
    RGB(0.836801, 0.329105, 0.430905),
    RGB(0.840155, 0.333580, 0.427455),
    RGB(0.843484, 0.338062, 0.424013),
    RGB(0.846788, 0.342551, 0.420579),
    RGB(0.850066, 0.347048, 0.417153),
    RGB(0.853319, 0.351553, 0.413734),
    RGB(0.856547, 0.356066, 0.410322),
    RGB(0.859750, 0.360588, 0.406917),
    RGB(0.862927, 0.365119, 0.403519),
    RGB(0.866078, 0.369660, 0.400126),
    RGB(0.869203, 0.374212, 0.396738),
    RGB(0.872303, 0.378774, 0.393355),
    RGB(0.875376, 0.383347, 0.389976),
    RGB(0.878423, 0.387932, 0.386600),
    RGB(0.881443, 0.392529, 0.383229),
    RGB(0.884436, 0.397139, 0.379860),
    RGB(0.887402, 0.401762, 0.376494),
    RGB(0.890340, 0.406398, 0.373130),
    RGB(0.893250, 0.411048, 0.369768),
    RGB(0.896131, 0.415712, 0.366407),
    RGB(0.898984, 0.420392, 0.363047),
    RGB(0.901807, 0.425087, 0.359688),
    RGB(0.904601, 0.429797, 0.356329),
    RGB(0.907365, 0.434524, 0.352970),
    RGB(0.910098, 0.439268, 0.349610),
    RGB(0.912800, 0.444029, 0.346251),
    RGB(0.915471, 0.448807, 0.342890),
    RGB(0.918109, 0.453603, 0.339529),
    RGB(0.920714, 0.458417, 0.336166),
    RGB(0.923287, 0.463251, 0.332801),
    RGB(0.925825, 0.468103, 0.329435),
    RGB(0.928329, 0.472975, 0.326067),
    RGB(0.930798, 0.477867, 0.322697),
    RGB(0.933232, 0.482780, 0.319325),
    RGB(0.935630, 0.487712, 0.315952),
    RGB(0.937990, 0.492667, 0.312575),
    RGB(0.940313, 0.497642, 0.309197),
    RGB(0.942598, 0.502639, 0.305816),
    RGB(0.944844, 0.507658, 0.302433),
    RGB(0.947051, 0.512699, 0.299049),
    RGB(0.949217, 0.517763, 0.295662),
    RGB(0.951344, 0.522850, 0.292275),
    RGB(0.953428, 0.527960, 0.288883),
    RGB(0.955470, 0.533093, 0.285490),
    RGB(0.957469, 0.538250, 0.282096),
    RGB(0.959424, 0.543431, 0.278701),
    RGB(0.961336, 0.548636, 0.275305),
    RGB(0.963203, 0.553865, 0.271909),
    RGB(0.965024, 0.559118, 0.268513),
    RGB(0.966798, 0.564396, 0.265118),
    RGB(0.968526, 0.569700, 0.261721),
    RGB(0.970205, 0.575028, 0.258325),
    RGB(0.971835, 0.580382, 0.254931),
    RGB(0.973416, 0.585761, 0.251540),
    RGB(0.974947, 0.591165, 0.248151),
    RGB(0.976428, 0.596595, 0.244767),
    RGB(0.977856, 0.602051, 0.241387),
    RGB(0.979233, 0.607532, 0.238013),
    RGB(0.980556, 0.613039, 0.234646),
    RGB(0.981826, 0.618572, 0.231287),
    RGB(0.983041, 0.624131, 0.227937),
    RGB(0.984199, 0.629718, 0.224595),
    RGB(0.985301, 0.635330, 0.221265),
    RGB(0.986345, 0.640969, 0.217948),
    RGB(0.987332, 0.646633, 0.214648),
    RGB(0.988260, 0.652325, 0.211364),
    RGB(0.989128, 0.658043, 0.208100),
    RGB(0.989935, 0.663787, 0.204859),
    RGB(0.990681, 0.669558, 0.201642),
    RGB(0.991365, 0.675355, 0.198453),
    RGB(0.991985, 0.681179, 0.195295),
    RGB(0.992541, 0.687030, 0.192170),
    RGB(0.993032, 0.692907, 0.189084),
    RGB(0.993456, 0.698810, 0.186041),
    RGB(0.993814, 0.704741, 0.183043),
    RGB(0.994103, 0.710698, 0.180097),
    RGB(0.994324, 0.716681, 0.177208),
    RGB(0.994474, 0.722691, 0.174381),
    RGB(0.994553, 0.728728, 0.171622),
    RGB(0.994561, 0.734791, 0.168938),
    RGB(0.994495, 0.740880, 0.166335),
    RGB(0.994355, 0.746995, 0.163821),
    RGB(0.994141, 0.753137, 0.161404),
    RGB(0.993851, 0.759304, 0.159092),
    RGB(0.993482, 0.765499, 0.156891),
    RGB(0.993033, 0.771720, 0.154808),
    RGB(0.992505, 0.777967, 0.152855),
    RGB(0.991897, 0.784239, 0.151042),
    RGB(0.991209, 0.790537, 0.149377),
    RGB(0.990439, 0.796859, 0.147870),
    RGB(0.989587, 0.803205, 0.146529),
    RGB(0.988648, 0.809579, 0.145357),
    RGB(0.987621, 0.815978, 0.144363),
    RGB(0.986509, 0.822401, 0.143557),
    RGB(0.985314, 0.828846, 0.142945),
    RGB(0.984031, 0.835315, 0.142528),
    RGB(0.982653, 0.841812, 0.142303),
    RGB(0.981190, 0.848329, 0.142279),
    RGB(0.979644, 0.854866, 0.142453),
    RGB(0.977995, 0.861432, 0.142808),
    RGB(0.976265, 0.868016, 0.143351),
    RGB(0.974443, 0.874622, 0.144061),
    RGB(0.972530, 0.881250, 0.144923),
    RGB(0.970533, 0.887896, 0.145919),
    RGB(0.968443, 0.894564, 0.147014),
    RGB(0.966271, 0.901249, 0.148180),
    RGB(0.964021, 0.907950, 0.149370),
    RGB(0.961681, 0.914672, 0.150520),
    RGB(0.959276, 0.921407, 0.151566),
    RGB(0.956808, 0.928152, 0.152409),
    RGB(0.954287, 0.934908, 0.152921),
    RGB(0.951726, 0.941671, 0.152925),
    RGB(0.949151, 0.948435, 0.152178),
    RGB(0.946602, 0.955190, 0.150328),
    RGB(0.944152, 0.961916, 0.146861),
    RGB(0.941896, 0.968590, 0.140956),
    RGB(0.940015, 0.975158, 0.131326)
], 30))

register_gradient_colors(:viridis, sample_evenly([
    RGB(0.267004, 0.004874, 0.329415),
    RGB(0.268510, 0.009605, 0.335427),
    RGB(0.269944, 0.014625, 0.341379),
    RGB(0.271305, 0.019942, 0.347269),
    RGB(0.272594, 0.025563, 0.353093),
    RGB(0.273809, 0.031497, 0.358853),
    RGB(0.274952, 0.037752, 0.364543),
    RGB(0.276022, 0.044167, 0.370164),
    RGB(0.277018, 0.050344, 0.375715),
    RGB(0.277941, 0.056324, 0.381191),
    RGB(0.278791, 0.062145, 0.386592),
    RGB(0.279566, 0.067836, 0.391917),
    RGB(0.280267, 0.073417, 0.397163),
    RGB(0.280894, 0.078907, 0.402329),
    RGB(0.281446, 0.084320, 0.407414),
    RGB(0.281924, 0.089666, 0.412415),
    RGB(0.282327, 0.094955, 0.417331),
    RGB(0.282656, 0.100196, 0.422160),
    RGB(0.282910, 0.105393, 0.426902),
    RGB(0.283091, 0.110553, 0.431554),
    RGB(0.283197, 0.115680, 0.436115),
    RGB(0.283229, 0.120777, 0.440584),
    RGB(0.283187, 0.125848, 0.444960),
    RGB(0.283072, 0.130895, 0.449241),
    RGB(0.282884, 0.135920, 0.453427),
    RGB(0.282623, 0.140926, 0.457517),
    RGB(0.282290, 0.145912, 0.461510),
    RGB(0.281887, 0.150881, 0.465405),
    RGB(0.281412, 0.155834, 0.469201),
    RGB(0.280868, 0.160771, 0.472899),
    RGB(0.280255, 0.165693, 0.476498),
    RGB(0.279574, 0.170599, 0.479997),
    RGB(0.278826, 0.175490, 0.483397),
    RGB(0.278012, 0.180367, 0.486697),
    RGB(0.277134, 0.185228, 0.489898),
    RGB(0.276194, 0.190074, 0.493001),
    RGB(0.275191, 0.194905, 0.496005),
    RGB(0.274128, 0.199721, 0.498911),
    RGB(0.273006, 0.204520, 0.501721),
    RGB(0.271828, 0.209303, 0.504434),
    RGB(0.270595, 0.214069, 0.507052),
    RGB(0.269308, 0.218818, 0.509577),
    RGB(0.267968, 0.223549, 0.512008),
    RGB(0.266580, 0.228262, 0.514349),
    RGB(0.265145, 0.232956, 0.516599),
    RGB(0.263663, 0.237631, 0.518762),
    RGB(0.262138, 0.242286, 0.520837),
    RGB(0.260571, 0.246922, 0.522828),
    RGB(0.258965, 0.251537, 0.524736),
    RGB(0.257322, 0.256130, 0.526563),
    RGB(0.255645, 0.260703, 0.528312),
    RGB(0.253935, 0.265254, 0.529983),
    RGB(0.252194, 0.269783, 0.531579),
    RGB(0.250425, 0.274290, 0.533103),
    RGB(0.248629, 0.278775, 0.534556),
    RGB(0.246811, 0.283237, 0.535941),
    RGB(0.244972, 0.287675, 0.537260),
    RGB(0.243113, 0.292092, 0.538516),
    RGB(0.241237, 0.296485, 0.539709),
    RGB(0.239346, 0.300855, 0.540844),
    RGB(0.237441, 0.305202, 0.541921),
    RGB(0.235526, 0.309527, 0.542944),
    RGB(0.233603, 0.313828, 0.543914),
    RGB(0.231674, 0.318106, 0.544834),
    RGB(0.229739, 0.322361, 0.545706),
    RGB(0.227802, 0.326594, 0.546532),
    RGB(0.225863, 0.330805, 0.547314),
    RGB(0.223925, 0.334994, 0.548053),
    RGB(0.221989, 0.339161, 0.548752),
    RGB(0.220057, 0.343307, 0.549413),
    RGB(0.218130, 0.347432, 0.550038),
    RGB(0.216210, 0.351535, 0.550627),
    RGB(0.214298, 0.355619, 0.551184),
    RGB(0.212395, 0.359683, 0.551710),
    RGB(0.210503, 0.363727, 0.552206),
    RGB(0.208623, 0.367752, 0.552675),
    RGB(0.206756, 0.371758, 0.553117),
    RGB(0.204903, 0.375746, 0.553533),
    RGB(0.203063, 0.379716, 0.553925),
    RGB(0.201239, 0.383670, 0.554294),
    RGB(0.199430, 0.387607, 0.554642),
    RGB(0.197636, 0.391528, 0.554969),
    RGB(0.195860, 0.395433, 0.555276),
    RGB(0.194100, 0.399323, 0.555565),
    RGB(0.192357, 0.403199, 0.555836),
    RGB(0.190631, 0.407061, 0.556089),
    RGB(0.188923, 0.410910, 0.556326),
    RGB(0.187231, 0.414746, 0.556547),
    RGB(0.185556, 0.418570, 0.556753),
    RGB(0.183898, 0.422383, 0.556944),
    RGB(0.182256, 0.426184, 0.557120),
    RGB(0.180629, 0.429975, 0.557282),
    RGB(0.179019, 0.433756, 0.557430),
    RGB(0.177423, 0.437527, 0.557565),
    RGB(0.175841, 0.441290, 0.557685),
    RGB(0.174274, 0.445044, 0.557792),
    RGB(0.172719, 0.448791, 0.557885),
    RGB(0.171176, 0.452530, 0.557965),
    RGB(0.169646, 0.456262, 0.558030),
    RGB(0.168126, 0.459988, 0.558082),
    RGB(0.166617, 0.463708, 0.558119),
    RGB(0.165117, 0.467423, 0.558141),
    RGB(0.163625, 0.471133, 0.558148),
    RGB(0.162142, 0.474838, 0.558140),
    RGB(0.160665, 0.478540, 0.558115),
    RGB(0.159194, 0.482237, 0.558073),
    RGB(0.157729, 0.485932, 0.558013),
    RGB(0.156270, 0.489624, 0.557936),
    RGB(0.154815, 0.493313, 0.557840),
    RGB(0.153364, 0.497000, 0.557724),
    RGB(0.151918, 0.500685, 0.557587),
    RGB(0.150476, 0.504369, 0.557430),
    RGB(0.149039, 0.508051, 0.557250),
    RGB(0.147607, 0.511733, 0.557049),
    RGB(0.146180, 0.515413, 0.556823),
    RGB(0.144759, 0.519093, 0.556572),
    RGB(0.143343, 0.522773, 0.556295),
    RGB(0.141935, 0.526453, 0.555991),
    RGB(0.140536, 0.530132, 0.555659),
    RGB(0.139147, 0.533812, 0.555298),
    RGB(0.137770, 0.537492, 0.554906),
    RGB(0.136408, 0.541173, 0.554483),
    RGB(0.135066, 0.544853, 0.554029),
    RGB(0.133743, 0.548535, 0.553541),
    RGB(0.132444, 0.552216, 0.553018),
    RGB(0.131172, 0.555899, 0.552459),
    RGB(0.129933, 0.559582, 0.551864),
    RGB(0.128729, 0.563265, 0.551229),
    RGB(0.127568, 0.566949, 0.550556),
    RGB(0.126453, 0.570633, 0.549841),
    RGB(0.125394, 0.574318, 0.549086),
    RGB(0.124395, 0.578002, 0.548287),
    RGB(0.123463, 0.581687, 0.547445),
    RGB(0.122606, 0.585371, 0.546557),
    RGB(0.121831, 0.589055, 0.545623),
    RGB(0.121148, 0.592739, 0.544641),
    RGB(0.120565, 0.596422, 0.543611),
    RGB(0.120092, 0.600104, 0.542530),
    RGB(0.119738, 0.603785, 0.541400),
    RGB(0.119512, 0.607464, 0.540218),
    RGB(0.119423, 0.611141, 0.538982),
    RGB(0.119483, 0.614817, 0.537692),
    RGB(0.119699, 0.618490, 0.536347),
    RGB(0.120081, 0.622161, 0.534946),
    RGB(0.120638, 0.625828, 0.533488),
    RGB(0.121380, 0.629492, 0.531973),
    RGB(0.122312, 0.633153, 0.530398),
    RGB(0.123444, 0.636809, 0.528763),
    RGB(0.124780, 0.640461, 0.527068),
    RGB(0.126326, 0.644107, 0.525311),
    RGB(0.128087, 0.647749, 0.523491),
    RGB(0.130067, 0.651384, 0.521608),
    RGB(0.132268, 0.655014, 0.519661),
    RGB(0.134692, 0.658636, 0.517649),
    RGB(0.137339, 0.662252, 0.515571),
    RGB(0.140210, 0.665859, 0.513427),
    RGB(0.143303, 0.669459, 0.511215),
    RGB(0.146616, 0.673050, 0.508936),
    RGB(0.150148, 0.676631, 0.506589),
    RGB(0.153894, 0.680203, 0.504172),
    RGB(0.157851, 0.683765, 0.501686),
    RGB(0.162016, 0.687316, 0.499129),
    RGB(0.166383, 0.690856, 0.496502),
    RGB(0.170948, 0.694384, 0.493803),
    RGB(0.175707, 0.697900, 0.491033),
    RGB(0.180653, 0.701402, 0.488189),
    RGB(0.185783, 0.704891, 0.485273),
    RGB(0.191090, 0.708366, 0.482284),
    RGB(0.196571, 0.711827, 0.479221),
    RGB(0.202219, 0.715272, 0.476084),
    RGB(0.208030, 0.718701, 0.472873),
    RGB(0.214000, 0.722114, 0.469588),
    RGB(0.220124, 0.725509, 0.466226),
    RGB(0.226397, 0.728888, 0.462789),
    RGB(0.232815, 0.732247, 0.459277),
    RGB(0.239374, 0.735588, 0.455688),
    RGB(0.246070, 0.738910, 0.452024),
    RGB(0.252899, 0.742211, 0.448284),
    RGB(0.259857, 0.745492, 0.444467),
    RGB(0.266941, 0.748751, 0.440573),
    RGB(0.274149, 0.751988, 0.436601),
    RGB(0.281477, 0.755203, 0.432552),
    RGB(0.288921, 0.758394, 0.428426),
    RGB(0.296479, 0.761561, 0.424223),
    RGB(0.304148, 0.764704, 0.419943),
    RGB(0.311925, 0.767822, 0.415586),
    RGB(0.319809, 0.770914, 0.411152),
    RGB(0.327796, 0.773980, 0.406640),
    RGB(0.335885, 0.777018, 0.402049),
    RGB(0.344074, 0.780029, 0.397381),
    RGB(0.352360, 0.783011, 0.392636),
    RGB(0.360741, 0.785964, 0.387814),
    RGB(0.369214, 0.788888, 0.382914),
    RGB(0.377779, 0.791781, 0.377939),
    RGB(0.386433, 0.794644, 0.372886),
    RGB(0.395174, 0.797475, 0.367757),
    RGB(0.404001, 0.800275, 0.362552),
    RGB(0.412913, 0.803041, 0.357269),
    RGB(0.421908, 0.805774, 0.351910),
    RGB(0.430983, 0.808473, 0.346476),
    RGB(0.440137, 0.811138, 0.340967),
    RGB(0.449368, 0.813768, 0.335384),
    RGB(0.458674, 0.816363, 0.329727),
    RGB(0.468053, 0.818921, 0.323998),
    RGB(0.477504, 0.821444, 0.318195),
    RGB(0.487026, 0.823929, 0.312321),
    RGB(0.496615, 0.826376, 0.306377),
    RGB(0.506271, 0.828786, 0.300362),
    RGB(0.515992, 0.831158, 0.294279),
    RGB(0.525776, 0.833491, 0.288127),
    RGB(0.535621, 0.835785, 0.281908),
    RGB(0.545524, 0.838039, 0.275626),
    RGB(0.555484, 0.840254, 0.269281),
    RGB(0.565498, 0.842430, 0.262877),
    RGB(0.575563, 0.844566, 0.256415),
    RGB(0.585678, 0.846661, 0.249897),
    RGB(0.595839, 0.848717, 0.243329),
    RGB(0.606045, 0.850733, 0.236712),
    RGB(0.616293, 0.852709, 0.230052),
    RGB(0.626579, 0.854645, 0.223353),
    RGB(0.636902, 0.856542, 0.216620),
    RGB(0.647257, 0.858400, 0.209861),
    RGB(0.657642, 0.860219, 0.203082),
    RGB(0.668054, 0.861999, 0.196293),
    RGB(0.678489, 0.863742, 0.189503),
    RGB(0.688944, 0.865448, 0.182725),
    RGB(0.699415, 0.867117, 0.175971),
    RGB(0.709898, 0.868751, 0.169257),
    RGB(0.720391, 0.870350, 0.162603),
    RGB(0.730889, 0.871916, 0.156029),
    RGB(0.741388, 0.873449, 0.149561),
    RGB(0.751884, 0.874951, 0.143228),
    RGB(0.762373, 0.876424, 0.137064),
    RGB(0.772852, 0.877868, 0.131109),
    RGB(0.783315, 0.879285, 0.125405),
    RGB(0.793760, 0.880678, 0.120005),
    RGB(0.804182, 0.882046, 0.114965),
    RGB(0.814576, 0.883393, 0.110347),
    RGB(0.824940, 0.884720, 0.106217),
    RGB(0.835270, 0.886029, 0.102646),
    RGB(0.845561, 0.887322, 0.099702),
    RGB(0.855810, 0.888601, 0.097452),
    RGB(0.866013, 0.889868, 0.095953),
    RGB(0.876168, 0.891125, 0.095250),
    RGB(0.886271, 0.892374, 0.095374),
    RGB(0.896320, 0.893616, 0.096335),
    RGB(0.906311, 0.894855, 0.098125),
    RGB(0.916242, 0.896091, 0.100717),
    RGB(0.926106, 0.897330, 0.104071),
    RGB(0.935904, 0.898570, 0.108131),
    RGB(0.945636, 0.899815, 0.112838),
    RGB(0.955300, 0.901065, 0.118128),
    RGB(0.964894, 0.902323, 0.123941),
    RGB(0.974417, 0.903590, 0.130215),
    RGB(0.983868, 0.904867, 0.136897),
    RGB(0.993248, 0.906157, 0.143936)
], 30))

# end of matplotlib colormaps
# ----------------------------------------------------------------------
# ----------------------------------------------------------------------

ext: jl


text:
function calculate_reaction_rate(problem::Dict{String,Any}, Weights::Array{Float64,1}, k_activator::Float64, n_activator::Float64, x::Float64)::Float64

    # Incorporate problem dictionary data
    kcat = problem["Catalytic_rate"] * 3600
    conc_F6P = problem["Fructose_6_Phosphate_concentration"]
    k_F6P = problem["F6P_binding_constant"]
    conc_ATP = problem["Adenosine_Tri_Phosphate_concentration"]
    k_ATP = problem["ATP_binding_constant"]
    E1 = problem["Phosphofructokinase_concentration"]

    # Calculate maximum reaction rate
    r1 = kcat * E1 * (conc_F6P / (k_F6P + conc_F6P)) * (conc_ATP / (k_ATP + conc_ATP))

    # Calculate activation control variable, v
    v1 = calculate_transcriptional_control_variable(Weights, k_activator, n_activator, x)

    # Calculate regulated reaction rate
    r̂ = r1 * v1
    return r̂
end


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using CSV#
using DataFrames
using Turing
using Dates
using StatsBase
using Distributions
using CovidSurvey
# using Gadfly
using StatsPlots
using LaTeXStrings
using NamedTupleTools
using StatsFuns
using Random
using TransformVariables, LogDensityProblems
using Optim
using ReverseDiff
using BSON
using Underscores
import TimeSeriesDecompositions as TSD
plotlyjs()

# using DataFramesMeta
# plotlyjs()
const dateform = "yyyymmdd"
include( projectdir("scripts/survey/rawsurvey_utils.jl") )

## ============================================================================
# util functions

function plot_contacts_per_day(df)
    gp = groupby(df[:,[:date]], :date, sort=true)
    cb = combine(gp, nrow)
    return plot(cb.date, cb.nrow, ylab="# contacts", xticks=:native)
end

function support_index(stepsize, support_position)
    support_position == :start && return 1
    support_position == :end && return stepsize
    iseven(stepsize) && return stepsize ÷ 2
    return stepsize ÷ 2 + 1
end

function sampledown(dates, stepsize=1, support_position=:middle)
    @assert issorted(dates)
    @assert all(unique(dates) .== first(dates):Day(1):last(dates))

    s = support_index(stepsize, support_position)
    n = length(dates)

    new_dates = Vector{Date}(undef, n)
    date = first(dates) + Day(s - 1)
    i = 1
    d_old = first(dates)
    for (j,d) in enumerate(dates)
        new_dates[j] = date
        if d == (d_old + Day(1))
            i += 1
            i % stepsize == 0 && ( date += Day(stepsize))
            d_old = date
        end
    end
    new_dates
end

## =============================================================================
# zero inflated analysis

function main(aggregate=1)
    plotlyjs()
    for region in 1:5
        @info "load data region $region"
        rawcontacts = readcontacts( normpath( homedir(), "data/covidsurvey/rawcontacts.csv"), select_region=region )
        plot_contacts_per_day(rawcontacts) |> display
        threshold = (family=50, colleagues=100, friends=100, strangers=1000)
        contacts = filteroutliers(rawcontacts, threshold) #NOTE National: p = 99.9
        contacts[!,:date] = sampledown(contacts.date, aggregate)
        contacts[!,:total] .= 0
        for c in ctype
            contacts[!,:total] += contacts[:,c]
        end
        # contacts = contacts[:, [:total, :strangers, :date]]

        @info "zero inflated regression region $region"
        surveyproblem = SurveyInferenceProblem(
            ZINegBinomialProblem,
            datacols(contacts),
            contacts
        )

        results = optimize(surveyproblem, LBFGS(),
            Optim.Options(show_trace=true, iterations=2000);
        )

        contactrates = DataFrame(results)
        save(projectdir("data","raw_contact_rates_region=$(region)_aggregation=$(aggregate).csv"), contactrates)

        let
            fname = projectdir("scripts/survey/","survey_inference_results_method=$(results.problemtype)_region=$(region)_aggregation=$(aggregate).bson")
            save(fname, θ = results)
        end

        plot_confidence_timeseries(results) |> display
        # pyplot()
        # p = plot_survey_histograms(results, day=nothing)
        # savefig(p, projectdir("figures/histogram_survey_data_vs_fake_data_region=$region.png"))
        # plotlyjs()

        # @info "decomposition region $region"
        # slowtrends = TSD.decompose(contactrates, Optim.Options(show_trace=true, iterations=2000); fdir=projectdir("figures/"))
        # save(projectdir("data","smoothed_contact_rates_region=$(region)_aggregation=$(aggregate).csv"), slowtrends)
    end
end

main(7)

## =========================================================================================================
#  Simple data cleaning with mean instead of zero inflated Negative Binomial
function rollingmean(xs::AbstractVector, window)
    n  = length(xs)-window+1
    ys = Vector{Float64}(undef, n)
    for i in 1:n
        x = xs[i:i+window-1]
        ys[i] = mean(skipmissing(x))
        @assert isfinite(ys[i]) "infinite at i=$i"
    end
    return ys
end

function rollingmean(df, window)
    @assert isodd(window)
    Δ = window ÷ 2
    smoothed = DataFrame( :date => df.date[1+Δ:end-Δ] )
    for col in names(df, Not(:date))
        @info col
        smoothed[!,col] = rollingmean(df[:,col], window)
        @assert all( isfinite.(smoothed[!,col]) ) "$(findall(x-> !isfinite(x), smoothed[!,col]))"
    end
    smoothed
end

function main(aggregate=1)
    plotlyjs()
    threshold = (family=50, colleagues=100, friends=100, strangers=1000)
    @progress for region in 1:5
        @info "load data region $region"
        contacts = @_ normpath( homedir(), "data/covidsurvey/rawcontacts.csv") |>
            readcontacts(__ , select_region=region ) |>
            filteroutliers(__, threshold) #NOTE National: p = 99.9

        contacts[!,:date] = sampledown(contacts.date, aggregate)
        contacts[!,:total] .= 0
        for c in ctype
            contacts[!,:total] += contacts[:,c]
        end
        # contacts = contacts[:, [:total, :strangers, :date]]

        contactrates = @_ contacts |>
             groupby(__, :date; sort=true) |>
             combine(__, names(__, Not("date")) .=> mean; renamecols=false) |>
             rollingmean(__, 7) |>
             rename(__, "date"=>"dates")

        save(projectdir("data/contacts/dk","averaged_contact_rates_region=$(region).csv"), contactrates)

        # @info "decomposition region $region"
        # slowtrends = TSD.decompose(contactrates, Optim.Options(show_trace=true, iterations=2000); fdir=projectdir("figures/"))
        # save(projectdir("data","smoothed_mean_contact_rates_region=$region.csv"), slowtrends)
    end
end

main(1)

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
"""
`Energy` - Compute the energy

E = -0.5 * J ∑ xᵢxⱼ

over a periodic N×N square lattice.  Note that this counts each bond twice,
hence the factor of 0.5

### Fields
* `x` - State of the lattice
### Optional Fields
* 'J=1' - Coupling constant>0
"""
function Energy(x::Tx; J::TF=1.0) where {TI<:Integer, Tx<:AbstractArray{TI}, TF<:AbstractFloat}
    E = 0.0;

    L= size(x)[1];
    for i in 1:L, j in 1:L
        # sum of energy contributions
        k, l = mod1(i+1,L), j;
        E += - 0.5 * J * x[i,j] * x[k,l];

        k,l = mod1(i-1,L), j;
        E += - 0.5 * J * x[i,j] * x[k,l];

        k, l = i, mod1(j+1,L);;
        E += - 0.5 * J * x[i,j] * x[k,l];

        k, l = i, mod1(j-1,L);
        E += - 0.5 * J * x[i,j] * x[k,l];
    end
    return E
end

"""
`Magnetization` - Compute the magnetization

M = ∑ xᵢ

over a periodic N×N square lattice.

### Fields
* `x` - State of the lattice
"""
function Magnetization(x::Tx) where {TI<:Integer, Tx<:AbstractArray{TI}}
    return sum(x)
end

"""
`ClusterNeighbors!` - Mark on the cluster map, `c`, where the current cluster
continues.

### Fields
* `c` - Cluster map
* `x` - State of the lattice
* `i,j` - Indices of the current site
* `D` - Direction of the current cluster (+/-1)
* `id` - Cluster id
* `L` - Lattice size, L × L
"""
function ClusterNeighbors!(c, x, i, j, D, id, L)

    k,l = mod1(i+1,L), j;
    if (x[k,l] == D && c[k,l]==0)
        c[k,l] = id;
        ClusterNeighbors!(c, x, k, l, D, id, L)
    end

    k,l = mod1(i-1,L), j;
    if (x[k,l]==D && c[k,l]==0)
        c[k,l] = id;
        ClusterNeighbors!(c, x, k, l, D, id, L)
    end

    k,l = i, mod1(j+1,L);
    if (x[k,l]==D && c[k,l]==0)
        c[k,l] = id;
        ClusterNeighbors!(c, x, k, l, D, id, L)
    end

    k,l = i, mod1(j-1,L);
    if (x[k,l]==D && c[k,l]==0)
        c[k,l] = id;
        ClusterNeighbors!(c, x, k, l, D, id, L)
    end
    c
end

"""
`Cluster` - Find all the clusters in the lattice.  Note that the numbering is
arbitrary.

### Fields
* `x` - State of the lattice
"""
function Clusters(x::Tx) where {TI<:Integer, Tx<:AbstractArray{TI}}
    L= size(x)[1];
    id = 0;
    c = zeros(Int, L, L);

    for i in 1:L, j in 1:L
        if(c[i,j]==0)
            id+=1;
            c[i,j] = id;
            D = x[i,j];
            ClusterNeighbors!(c, x, i, j, D, id, L)
        end
    end

    return c
end


text:
@testset "Frequency pitches" begin
    @testset "types ok and compact" begin
        @test isbitstype(FreqInterval)
        @test isbitstype(FreqIC)
        @test isbitstype(Pitch{FreqInterval})
        @test isbitstype(Pitch{FreqIC})
    end

    @testset "constructors" begin
        @test logfreqi(log(3.8)) == freqi(3.8)
        @test logfreqic(log(3.8)) == freqic(3.8)
        @test logfreqic(log(3.8) + log(2)) ≈ freqic(3.8)
        @test logfreqp(log(441)) == freqp(441)
        @test logfreqpc(log(441)) == freqpc(441)
        @test logfreqpc(log(441) + log(2)) ≈ freqpc(441)
    end

    @testset "printing" begin
        #@test string(freqi(3.8))  == "fr3.8"
        #@test string(freqic(3.8)) == "fr[1.9]"
        #@test string(freqp(441)) == "441.0Hz"
        #@test string(freqpc(441)) == "[1.7226562500000004]Hz"
    end

    @testset "hashing" begin
        @test hash(logfreqi(log(3.8))) == hash(freqi(3.8))
        @test hash(logfreqic(log(3.8))) == hash(freqic(3.8))
        @test hash(logfreqi(log(441))) == hash(freqi(441))
        @test hash(logfreqic(log(441))) == hash(freqic(441))
    end

    @testset "interval interface" begin
        @test freqi(3/2) + freqi(4/3) ≈ freqi(2)
        @test freqi(2/3) + freqi(4/3) ≈ freqi(8/9)
        @test (-freqi(3/2)) + freqi(4/3) ≈ freqi(8/9)
        @test freqi(4/3) - freqi(3/2) ≈ -freqi(9/8)
        @test -freqi(3/2) ≈ freqi(2/3)
        @test -freqi(4/5) ≈ unison(FreqInterval) - freqi(4/5)

        @test unison(FreqInterval) == freqi(1)
        @test unison(freqi(0.42)) == freqi(1)

        @test freqi(3/2) * 2 == freqi(9/4)
        @test freqi(3/2) * -2 == freqi(4/9)
        @test 2 * freqi(3/2) == freqi(9/4)
        @test -2 * freqi(3/2) == freqi(4/9)

        @test sign(freqi(3/2)) == 1
        @test sign(freqi(2/3)) == -1
        @test sign(freqi(1)) == 0
        @test abs(freqi(3/2)) == freqi(3/2)
        @test abs(freqi(2/3)) ≈ freqi(3/2)
        @test abs(freqi(1)) == freqi(1)

        @test octave(FreqInterval) == freqi(2)
        @test chromsemi(FreqInterval) == logfreqi(log(2)/12)

        @test ic(freqi(3)) ≈ freqic(3)
        @test ic(freqi(3)) ≈ freqic(3/2)
        @test ic(freqi(3/2)) ≈ freqic(3)
        @test embed(freqi(3/2)) == freqi(3/2)
        @test intervaltype(FreqInterval) == FreqInterval
        @test intervalclasstype(FreqInterval) == FreqIC
    end

    @testset "interval class interface" begin
        @test freqic(3/2) + freqic(4/3) ≈ freqic(2)
        @test freqic(2/3) + freqic(4/3) ≈ freqic(8/9)
        @test (-freqic(3/2)) + freqic(4/3) ≈ freqic(8/9)
        @test freqic(4/3) - freqic(3/2) ≈ -freqic(9/8)
        @test -freqic(3/2) ≈ freqic(2/3)
        @test -freqic(4/5) ≈ unison(FreqIC) - freqic(4/5)

        @test unison(FreqIC) == freqic(1)
        @test unison(freqic(0.42)) == freqic(1)

        @test freqic(3/2) * 2 == freqic(9/4)
        @test freqic(3/2) * -2 == freqic(4/9)
        @test 2 * freqic(3/2) == freqic(9/4)
        @test -2 * freqic(3/2) == freqic(4/9)

        @test sign(freqic(3/2)) == -1
        @test sign(freqic(2/3)) == 1
        @test sign(freqic(1)) == 0
        @test abs(freqic(3/2)) == freqic(2/3)
        @test abs(freqic(2/3)) ≈ freqic(2/3)
        @test abs(freqic(1)) == freqic(1)

        @test octave(FreqIC) == freqic(2)
        @test chromsemi(FreqIC) == logfreqic(log(2)/12)

        @test ic(freqic(3)) ≈ freqic(3)
        @test ic(freqic(3)) ≈ freqic(3/2)
        @test ic(freqic(3/2)) ≈ freqic(3)
        @test embed(freqic(3/2)) == freqi(3/2)
        @test intervaltype(FreqIC) == FreqInterval
        @test intervalclasstype(FreqIC) == FreqIC
    end

    @testset "pitch interface" begin
        @test topitch(freqi(441)) == freqp(441)
        @test tointerval(freqp(441)) == freqi(441)

        @test freqp(441) + freqi(2) ≈ freqp(882)
        @test freqp(441) + freqi(1/2) ≈ freqp(220.5)
        @test freqp(441) - freqi(2) ≈ freqp(220.5)
        @test freqp(441) - freqp(431.5) ≈ freqi(441/431.5)

        @test pc(freqp(441)) == freqpc(441)
        @test embed(freqp(441)) == freqp(441)
    end

    @testset "pitch class interface" begin
        @test topitch(freqic(441)) == freqpc(441)
        @test tointerval(freqpc(441)) == freqic(441)

        @test freqpc(441) + freqic(3/2) ≈ freqpc(661.5)
        @test freqpc(441) + freqic(2/3) ≈ freqpc(294)
        @test freqpc(441) - freqic(3/2) ≈ freqpc(294)
        @test freqpc(441) - freqpc(431.5) ≈ freqic(441/431.5)

        @test pc(freqpc(441)) == freqpc(441)
        @test embed(freqpc(441)) == logfreqp(mod(log(441), log(2)))
        @test embed(freqpc(441), 8) == freqp(441)
    end
end


text:
@testset "test And" begin
    @test And()(TRUE(), TRUE()) == TRUE()
    @test And()(TRUE(), FALSE()) == FALSE()
    @test And()(FALSE(), TRUE()) == FALSE()
    @test And()(FALSE(), FALSE()) == FALSE()
end

@testset "test Or" begin
    @test Or()(TRUE(), TRUE()) == TRUE()
    @test Or()(TRUE(), FALSE()) == TRUE()
    @test Or()(FALSE(), TRUE()) == TRUE()
    @test Or()(FALSE(), FALSE()) == FALSE()
end

@testset "test Not" begin
    @test Not()(TRUE()) == FALSE()
    @test Not()(FALSE()) == TRUE()
end

@testset "test Eq" begin
    @test Eq()(TRUE(), TRUE()) == TRUE()
    @test Eq()(TRUE(), FALSE()) == FALSE()
    @test Eq()(FALSE(), TRUE()) == FALSE()
    @test Eq()(FALSE(), FALSE()) == TRUE()
end

@testset "test Xor" begin
    @test Xor()(TRUE(), TRUE()) == FALSE()
    @test Xor()(TRUE(), FALSE()) == TRUE()
    @test Xor()(FALSE(), TRUE()) == TRUE()
    @test Xor()(FALSE(), FALSE()) == FALSE()
end


text:
"""
    generate random particle-hole-symmetric positive-definite hamiltonian.
"""
function test_random_particle_hole_positive_definite_hamiltonian(N::Int)
    @assert mod(N,2)==0
    S = rand(N).+0.5
    q1,r1 = qr(rand(N,N))
    q2,r2 = qr(rand(N,N))
    S = diagm(vcat(S,S))
    q = [q1 q2;q2 q1]
    return transpose(q)*S*q
end

"""
    test if the transform V is canonical
"""
function test_canonicality(V::Matrix,ϵ::Float64)
    @assert size(V)[1]==size(V)[2]
    N = size(V)[1]÷2
    η = diagm(vcat([1.0 for i=1:N],[-1.0 for i=1:N]))
    @test norm(transpose(V)*η*V-η)<ϵ
end

"""
    test if the transform V is canonical with zeromodes
"""
function test_canonicality_zeromode(V::Matrix,ϵ::Float64)
    @assert size(V)[1]==size(V)[2]
    N = size(V)[1]÷2
    η = diagm(vcat([1.0 for i=1:N],[-1.0 for i=1:N]))
    test_η = transpose(V)*η*V
    for i = 1:2N
        if abs(test_η[i,i])<1E-10
            test_η[i,i] = i>N ? -1.0 : 1.0
        end
    end
    @test norm(test_η-η)<ϵ
end

"""
    test if the matrix is diagonalized
"""
function test_diag(s::Matrix,ϵ::Float64=1E-10)
    @test norm(s-diagm(diag(s)))<ϵ
end

"""
    test if S = transpose(V)*A*V
    this will hold if use normal and abnormal deal zeromodes
"""
function test_right(s::Matrix,V::Matrix,A::Matrix,ϵ::Float64)
    N = size(s)[1]÷2
    @test norm(transpose(V)*A*V-s)<ϵ
end

"""
    test if my_inverse(V^T)*S*my_inverse(V) = A
    this will not hold if use normal and abnormal deal zeromodes
"""
function test_full_right(s::Matrix,V::Matrix,A::Matrix,ϵ::Float64)
    N = size(s)[1]÷2
    @test norm(transpose(V)*A*V-s)<ϵ
    Vs = V[1:N,1:N]
    Us = V[1:N,1+N:2N]
    iV = [transpose(Vs) -transpose(Us);-transpose(Us) transpose(Vs)]
    @test norm(A-transpose(iV)*s*iV)<ϵ
end

"""
    OBC Schwinger Boson hamiltonian
"""

function schwinger_boson_hamiltonian(N::Int,bc::String)
    A = 1
    λ = 1
    begin
        Q = zeros(N,N)
        for id in CartesianIndices((1:N,1:N))
            if id[1]-id[2]==1
                Q[id] = 1.0*A
            elseif id[2] - id[1] ==1
                Q[id] = -1.0*A
            end
        end

        if bc == "PBC"
            Q[N,1] = -1.0*A
            Q[1,N] =  1.0*A
        end
        Q = Q/2.0
        Q
    end
    ϵ = 5*1E-10
    Lambda = diagm([λ for i = 1:N])
    M = [Lambda Q;transpose(Q) Lambda]
    return M
end