text:
struct FiniteDifference{T<:Real,S<:SVector,LBC,RBC} <: AbstractDerivativeOperator{T}
    derivative_order    :: Int
    approximation_order :: Int
    dx                  :: Vector{T}
    dimension           :: Int
    stencil_length      :: Int
    stencil_coefs       :: Vector{S}
    boundary_point_count:: Tuple{Int,Int}
    boundary_length     :: Tuple{Int,Int}
    low_boundary_coefs  :: Ref{Vector{Vector{T}}}
    high_boundary_coefs :: Ref{Vector{Vector{T}}}
    boundary_condition  :: Ref{Tuple{Tuple{T,T,Any},Tuple{T,T,Any}}}
    t                   :: Ref{Int}

    function FiniteDifference{T,S,LBC,RBC}(derivative_order::Int,
                                           approximation_order::Int,
                                           dx::Vector{T},
                                           dimension::Int, BC) where
                                           {T<:Real,S<:SVector,LBC,RBC}
        # dimension            = dimension
        # dx                   = dx
        stencil_length       = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2
        bl                   = derivative_order + approximation_order
        boundary_length      = (bl,bl)
        bpc                  = stencil_length - div(stencil_length,2) + 1
        bpc_array            = [bpc,bpc]
        grid_step            = dx #unnecessary

        if any(x->x<zero(T),dx)
            error("All grid steps must be greater than 0.0!")
        end
        x                    = [zero(T); cumsum(dx)]
        low_boundary_coefs   = Vector{T}[]
        high_boundary_coefs  = Vector{T}[]

        stl_2 = div(stencil_length,2)

        stencil_coefs        =[convert(SVector{stencil_length, T}, calculate_weights(derivative_order, zero(T),
                               x[i-stl_2 : i+stl_2].-x[i])) for i in stl_2+1:dimension-stl_2]

        left_bndry = initialize_left_boundary!(Val{:FD},low_boundary_coefs,stencil_coefs[1:bl],
                                                BC,derivative_order,grid_step[1:bl],bl,bpc_array,LBC)

        right_bndry = initialize_right_boundary!(Val{:FD},high_boundary_coefs,stencil_coefs[end-bl+1:end],
                                                 BC,derivative_order,grid_step[end-bl+1:end],bl,bpc_array,RBC)

        boundary_condition = (left_bndry, right_bndry)
        boundary_point_count = (bpc_array[1],bpc_array[2])

        t = 0

        new(derivative_order, approximation_order, dx, dimension, stencil_length,
            stencil_coefs,
            boundary_point_count,
            boundary_length,
            low_boundary_coefs,
            high_boundary_coefs,
            boundary_condition,
            t
            )
    end
    FiniteDifference{T}(dorder::Int,aorder::Int,dx::Vector{T},dim::Int,LBC::Symbol,RBC::Symbol;BC=(zero(T),zero(T))) where {T<:Real} =
        FiniteDifference{T, SVector{dorder+aorder-1+(dorder+aorder)%2,T}, LBC, RBC}(dorder, aorder, dx, dim, BC)
end

function LinearAlgebra.mul!(x_temp::AbstractVector{T}, A::FiniteDifference{T}, x::AbstractVector{T}) where T<:Real
    convolve_BC_left!(x_temp, x, A)
    convolve_interior!(x_temp, x, A)
    convolve_BC_right!(x_temp, x, A)
    A.t[] += 1 # incrementing the internal time stamp
end

#=
    This function is used to update the boundary conditions especially if they evolve with
    time.
=#
function DiffEqBase.update_coefficients!(A::FiniteDifference{T,S,LBC,RBC};BC=nothing) where {T<:Real,S<:SVector,RBC,LBC}
    if BC != nothing
        LBC == :Robin ? (length(BC[1])==3 || error("Enter the new left boundary condition as a 1-tuple")) :
                        (length(BC[1])==1 || error("Robin BC needs a 3-tuple for left boundary condition"))

        RBC == :Robin ? length(BC[2])==3 || error("Enter the new right boundary condition as a 1-tuple") :
                        length(BC[2])==1 || error("Robin BC needs a 3-tuple for right boundary condition")

        left_bndry = initialize_left_boundary!(Val{:FD}, A.low_boundary_coefs[],A.stencil_coefs[1:A.boundary_length],BC,
                                               A.derivative_order,A.dx[1:A.boundary_length],A.boundary_length[1],LBC)

        right_bndry = initialize_right_boundary!(A.high_boundary_coefs[],A.stencil_coefs[end-A.boundary_length:end],BC,
                                                 A.derivative_order,A.dx[end-A.boundary_length:end],A.boundary_length[2],RBC)

        boundary_condition = (left_bndry, right_bndry)
        A.boundary_condition[] = boundary_condition
    end
end

#################################################################################################

#=
    This function sets us up to apply the boundary condition correctly. It returns a
    3-tuple which is basically the coefficients in the equation
                            a*u + b*du/dt = c(t)
    The RHS ie. 'c' can be a function of time as well and therefore it is implemented
    as an anonymous function.
=#
function initialize_left_boundary!(::Type{Val{:FD}},low_boundary_coefs,stencil_coefs,BC,derivative_order,grid_step::Vector{T},
                                   boundary_length,boundary_point_count,LBC) where T
    stencil_length = length(stencil_coefs[1])

    if LBC == :None
        #=
            Val{:FD} is type parametrization on symbols. There are two different right_None_BC!
            functions, one for LinearOperator and one for UpwindOperator. So to select the correct
            function without changing the name, this trick has been applied.
        =#
        boundary_point_count[1] = div(stencil_length,2)
        return (zero(T),zero(T),left_None_BC!(Val{:FD},low_boundary_coefs,stencil_length,derivative_order,
                                              grid_step,boundary_length)*BC[1])
    elseif LBC == :Neumann
        return (zero(T),one(T),left_Neumann_BC!(Val{:FD},low_boundary_coefs,stencil_length,derivative_order,
                                                grid_step,boundary_length)*BC[1])
    elseif LBC == :Robin
        error("LBC Robin not yet supported for irregular grid.")
        # return (BC[1][1],-BC[1][2],left_Robin_BC!(Val{:FD},low_boundary_coefs,stencil_length,
                                                   # BC[1],derivative_order,grid_step,
                                                   # boundary_length,dx)*BC[1][3])
    elseif LBC == :Dirichlet0
        return (one(T),zero(T),zero(T)*BC[1])

    elseif LBC == :Dirichlet
        typeof(BC[1]) <: Real ? ret = t->BC[1] : ret = BC[1]
        return (one(T),zero(T),ret)

    elseif LBC == :Neumann0
        return (zero(T),one(T),zero(T))

    elseif LBC == :periodic
        return (zero(T),zero(T),zero(T))

    else
        error("Unrecognized Boundary Type!")
    end
end

#=
    This function sets us up to apply the boundary condition correctly. It returns a
    3-tuple which is basically the coefficients in the equation
                            a*u + b*du/dt = c(t)
    The RHS ie. 'c' can be a function of time as well and therefore it is implemented
    as an anonymous function.
=#
function initialize_right_boundary!(::Type{Val{:FD}},high_boundary_coefs,stencil_coefs,BC,derivative_order,grid_step::Vector{T},
                                    boundary_length,boundary_point_count,RBC) where T
    stencil_length = length(stencil_coefs[1])

    if RBC == :None
        boundary_point_count[2] = div(stencil_length,2)
        #=
            Val{:FD} is type parametrization on symbols. There are two different right_None_BC!
            functions, one for LinearOperator and one for UpwindOperator. So to select the correct
            function without changing the name, this trick has been applied.
        =#
        return (zero(T),zero(T),right_None_BC!(Val{:FD},high_boundary_coefs,stencil_length,derivative_order,
                               grid_step,boundary_length)*BC[2])
    elseif RBC == :Neumann
        error("RBC Robin not yet supported for irregular grid.")
        # return (zero(T),one(T),right_Neumann_BC!(Val{:FD},high_boundary_coefs,stencil_length,derivative_order,
        #                           grid_step,boundary_length)*BC[2])
    elseif RBC == :Robin
        error("RBC Robin not yet supported for irregular grid.")
        # return (BC[2][1],BC[2][2],right_Robin_BC!(Val{:FD},high_boundary_coefs,stencil_length,
                                                    # BC[2],derivative_order,grid_step,
                                                    # boundary_length,dx)*BC[2][3])
    elseif RBC == :Dirichlet0
        return (one(T),zero(T),zero(T)*BC[2])

    elseif RBC == :Dirichlet
        typeof(BC[2]) <: Real ? ret = t->BC[2] : ret = BC[2]
        return (one(T),zero(T),ret)

    elseif RBC == :Neumann0
        return (zero(T),one(T),zero(T))

    elseif RBC == :periodic
        return (zero(T),zero(T),zero(T))

    else
        error("Unrecognized Boundary Type!")
    end
end

function left_None_BC!(::Type{Val{:FD}},low_boundary_coefs,stencil_length,derivative_order,
                       grid_step::Vector{T},boundary_length) where T
    # Fixes the problem excessive boundary points
    boundary_point_count = div(stencil_length,2)
    l_diff               = zero(T)
    mid                  = div(stencil_length,2)
    x = [0.0; cumsum(grid_step)]
    for i in 1 : boundary_point_count
        # One-sided stencils require more points for same approximation order
        # TODO: I don't know if this is the correct stencil length for i > 1?
        push!(low_boundary_coefs, calculate_weights(derivative_order, zero(T), x[1:boundary_length].-x[i]))
    end
    return l_diff
end

function right_None_BC!(::Type{Val{:FD}},high_boundary_coefs,stencil_length,derivative_order,
                        grid_step::Vector{T},boundary_length) where T
    boundary_point_count = div(stencil_length,2)
    high_temp            = zeros(T,boundary_length)
    flag                 = derivative_order*boundary_point_count%2
    aorder               = boundary_length - 1
    r_diff               = zero(T)
    x                    = [0.0; cumsum(grid_step)]

    for i in length(x) : -1 : length(x) - boundary_point_count + 1
        # One-sided stencils require more points for same approximation order
        push!(high_boundary_coefs, calculate_weights(derivative_order,x[i],x[end - boundary_length + 1 : end]))
    end
    return r_diff
end

function left_Neumann_BC!(::Type{Val{:FD}},low_boundary_coefs,stencil_length,derivative_order,
                          grid_step::Vector{T},boundary_length) where T
    boundary_point_count = stencil_length - div(stencil_length,2) + 1
    # first_order_coeffs   = zeros(T,boundary_length)
    # original_coeffs      = zeros(T,boundary_length)
    l_diff               = one(T)
    mid                  = div(stencil_length,2)+1
    x = [0.0;cumsum(grid_step)]

    first_order_coeffs = calculate_weights(1, zero(T), x[1:boundary_length])
    original_coeffs =  calculate_weights(derivative_order, zero(T), x[1:boundary_length])
    l_diff = original_coeffs[end]/first_order_coeffs[end]
    rmul!(first_order_coeffs, original_coeffs[end]/first_order_coeffs[end])
    # rmul!(original_coeffs, first_order_coeffs[end]/original_coeffs[end])
    @. original_coeffs = original_coeffs - first_order_coeffs
    # copyto!(first_order_coeffs, first_order_coeffs[1:end-1])
    push!(low_boundary_coefs, original_coeffs[1:end-1])

    for i in 2 : boundary_point_count
        #=  this means that a stencil will suffice ie. we dont't need to worry about the boundary point
            being considered in the low_boundary_coefs
        =#

        if i > mid
            pos=i
            push!(low_boundary_coefs, calculate_weights(derivative_order, x[pos], x[1:boundary_length]))
        else
            pos=i-1
            push!(low_boundary_coefs, append!([zero(T)],calculate_weights(derivative_order, x[pos], x[1:boundary_length])))
        end
    end

    return l_diff
end

#TODO:::
# function right_Neumann_BC!(::Type{Val{:FD}},high_boundary_coefs,stencil_length,derivative_order,
#                            grid_step::Vector{T},boundary_length) where T
#     boundary_point_count = stencil_length - div(stencil_length,2) + 1
#     flag                 = derivative_order*boundary_point_count%2
#     # original_coeffs      = zeros(T,boundary_length)
#     r_diff               = one(T)
#     mid                  = div(stencil_length,2)+1
#     x                    = [0.0; cumsum(grid_step)]
#     @show flag
#
#     # this part is to incorporate the value of first derivative at the right boundary
#     first_order_coeffs = calculate_weights(1, x[end-boundary_point_count+1], x[end - boundary_length + 1 : end])
#     original_coeffs = calculate_weights(derivative_order, x[end-boundary_point_count+1], x[end - boundary_length + 1 : end])
#
#     r_diff = original_coeffs[1]/first_order_coeffs[1]
#     rmul!(first_order_coeffs, original_coeffs[1]/first_order_coeffs[1])
#     # rmul!(original_coeffs, first_order_coeffs[1]/original_coeffs[1])
#     @. original_coeffs = original_coeffs - first_order_coeffs
#     # copyto!(first_order_coeffs, first_order_coeffs[1:end-1])
#     push!(high_boundary_coefs, original_coeffs[2:end])
#     #TODO doesnt work... i don't get it.
#     for i in length(x) - boundary_point_count+2 : length(x)
#         # push!(high_boundary_coefs,ones(T,boundary_length)*i)
#         # @show i
#         # @show length(x)-i+2
#         if i > length(x)-mid+1
#             pos=i
#             push!(high_boundary_coefs, calculate_weights(derivative_order, x[pos],x[end - boundary_length + 1 : end]))
#         else
#             pos=i-1
#             push!(high_boundary_coefs, append!([zero(T)], calculate_weights(derivative_order, x[pos],x[end - boundary_length  : end-1])))
#         end
#     end
#
#     reverse!(high_boundary_coefs)
#
#
#     return r_diff
# end

# function left_Robin_BC!(::Type{Val{:FD}},low_boundary_coefs,stencil_length,params,
#                         derivative_order,grid_step::Vector{T},boundary_length,dx) where T
#     boundary_point_count = stencil_length - div(stencil_length,2) + 1
#     first_order_coeffs   = zeros(T,boundary_length)
#     original_coeffs      = zeros(T,boundary_length)
#     l_diff               = one(T)
#     mid                  = div(stencil_length,2)+1
#
#     # in Robin BC the left boundary has opposite sign by convention
#     first_order_coeffs = -params[2]*calculate_weights(1, (0)*grid_step, collect(zero(T) : grid_step : (boundary_length-1)* grid_step))
#     first_order_coeffs[1] += dx*params[1]
#     original_coeffs =  calculate_weights(derivative_order, (0)*grid_step, collect(zero(T) : grid_step : (boundary_length-1) * grid_step))
#
#     l_diff = original_coeffs[end]/first_order_coeffs[end]
#     rmul!(first_order_coeffs, original_coeffs[end]/first_order_coeffs[end])
#     # rmul!(original_coeffs, first_order_coeffs[end]/original_coeffs[end])
#     @. original_coeffs = original_coeffs - first_order_coeffs
#     # copyto!(first_order_coeffs, first_order_coeffs[1:end-1])
#     push!(low_boundary_coefs, original_coeffs[1:end-1])
#
#     for i in 2 : boundary_point_count
#         #=
#             this means that a stencil will suffice ie. we dont't need to worry about the boundary point
#             being considered in the low_boundary_coefs
#         =#
#         if i > mid
#             pos=i-1
#             push!(low_boundary_coefs, calculate_weights(derivative_order, pos*grid_step, collect(zero(T) : grid_step : (boundary_length-1) * grid_step)))
#         else
#             pos=i-2
#             push!(low_boundary_coefs, append!([zero(T)],calculate_weights(derivative_order, pos*grid_step, collect(zero(T) : grid_step : (boundary_length-1) * grid_step))))
#         end
#     end
#
#     return l_diff
# end

# function right_Robin_BC!(::Type{Val{:FD}},high_boundary_coefs,stencil_length,params,
#                         derivative_order,grid_step::Vector{T},boundary_length,dx) where T
#     boundary_point_count = stencil_length - div(stencil_length,2) + 1
#     flag                 = derivative_order*boundary_point_count%2
#     original_coeffs      = zeros(T,boundary_length)
#     r_diff               = one(T)
#     mid                  = div(stencil_length,2)+1
#
#     first_order_coeffs = params[2]*calculate_weights(1, (boundary_point_count-1)*grid_step,
#                                                      collect(zero(T) : grid_step : (boundary_length-1) * grid_step))
#     first_order_coeffs[end] += dx*params[1]
#     reverse!(first_order_coeffs)
#     isodd(flag) ? negate!(first_order_coeffs) : nothing
#
#     copyto!(original_coeffs, calculate_weights(derivative_order, (boundary_point_count-1)*grid_step,
#                                              collect(zero(T) : grid_step : (boundary_length-1) * grid_step)))
#     reverse!(original_coeffs)
#     isodd(flag) ? negate!(original_coeffs) : nothing
#
#     r_diff = original_coeffs[1]/first_order_coeffs[1]
#     rmul!(first_order_coeffs, original_coeffs[1]/first_order_coeffs[1])
#     # rmul!(original_coeffs, first_order_coeffs[1]/original_coeffs[1])
#     @. original_coeffs = original_coeffs - first_order_coeffs
#     # copyto!(first_order_coeffs, first_order_coeffs[1:end-1])
#
#     for i in 2 : boundary_point_count
#         #=
#             this means that a stencil will suffice ie. we dont't need to worry about the boundary point
#             being considered in the high_boundary_coefs. Same code for low_boundary_coefs but reversed
#             at the end
#         =#
#         if i > mid
#             pos=i-1
#             push!(high_boundary_coefs, calculate_weights(derivative_order, pos*grid_step,
#                                                          collect(zero(T) : grid_step : (boundary_length-1) * grid_step)))
#         else
#             pos=i-2
#             push!(high_boundary_coefs, append!([zero(T)],
#                                                calculate_weights(derivative_order,pos*grid_step,
#                                                                  collect(zero(T) : grid_step : (boundary_length-1) * grid_step))))
#         end
#     end
#     if flag == 1
#         negate!(high_boundary_coefs)
#     end
#     reverse!(high_boundary_coefs)
#     push!(high_boundary_coefs, original_coeffs[2:end])
#     return r_diff
# end

#################################################################################################

(L::FiniteDifference)(u,p,t) = L*u
(L::FiniteDifference)(du,u,p,t) = mul!(du,L,u)
# get_LBC(::FiniteDifference{A,B,C,D}) where {A,B,C,D} = C
# get_RBC(::FiniteDifference{A,B,C,D}) where {A,B,C,D} = D


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

using Markdown
using InteractiveUtils

# ╔═╡ f502a2cc-ea34-11eb-0eb3-93d1cb0a9ad2
# In version 0.15.1 and later, Pluto.jl will automatically load these packages and their dependencies.

using CSV     # docs: https://csv.juliadata.org/stable/

# ╔═╡ e78d0da4-69b2-450a-b9c3-a7087a698618
using Plots   # docs: http://docs.juliaplots.org/latest/

# ╔═╡ 65c4e11f-7c12-46ef-8c9c-e97d139a80d5
using Dates   # docs: https://docs.julialang.org/en/v1/stdlib/Dates/

# ╔═╡ 1bd45843-5bdb-4f19-9f04-23d3917650db
using FFTW    # Fourier transforms; docs: http://www.fftw.org/doc/index.html#Top

# ╔═╡ 12032637-059c-4227-840a-55c32a5c402c
md"""
The code in this notebook accompanies the essay "Riding the Covid Coaster," by Brian Hayes, published 17 August 2021 at [http://bit-player.org/2021/riding-the-covid-coaster](http://bit-player.org/2021/riding-the-covid-coaster). It explores various approaches to understanding the trajectory of the Covid-19 epidemic in the United States during its first 18 months.

This is one of three notebooks with code for the project:

* covidcurve.jl covers the reconstruction of the Covid-19 incidence curve and attempts to apply Fourier methods, looking for periodic signals in the data. Also a look at time derivatives. Includes code for figures 1, 3, 4, 5, 6, 7, 8, 9, 12.

* SIR_model.jl deals with exponential growth in epidemic processes, and models based on that idea. Figures 10, 11, 13, 14, 16, 17, 18, 19, 20, 21, 23, 24.

* us_states_data.jl breaks down the national Covid data set into records for inividual states and regions. Code for Figure 22.

The underlying data for most of this work was collected and curated by _The New York Times_. The archive is available at [https://github.com/nytimes/covid-19-data](https://github.com/nytimes/covid-19-data). The version I used in the code below, including data from 21 January 2020 through 20 July 2021, is commit c3ab8c1beba1f4728d284c7b1e58d7074254aff8. You can access the identical set of files at [https://github.com/nytimes/covid-19-data/tree/c3ab8c1beba1f4728d284c7b1e58d7074254aff8](https://github.com/nytimes/covid-19-data/tree/c3ab8c1beba1f4728d284c7b1e58d7074254aff8). But you should probably get the most up-to-date version, with the latest data.

Within the _Times_ archive I have used the files within a folder named "rolling-averages," which has daily new cases rather than cumulative cases, and has precomputed the seven-day trailing average of the data sequences.

The code is written in the [Julia programming language](https://julialang.org/) and organized in [Pluto.jl notebooks](https://github.com/fonsp/Pluto.jl).

Instructions:

1. You'll need a working installation of Julia with the Pluto.jl package installed. Version 0.15.1 or later is needed to automatically find, load, and install packages.

2. In a terminal window, cd to the directory where you have stored this file. Start Julia, then at the "julia>" prompt, type "using Pluto; Pluto.run('covidcurve.jl')".

3. A browser window will open with this notebook displayed. Pluto will then download all the necessary packages and their dependencies, and will evaluate the cells in the proper order. This can take several minutes, and it can rev up your CPU fan. But then you're done.

"""

# ╔═╡ 47bc3ba1-7d7c-435d-8066-87fb07baa293
# Assume: the Times data archive is stored in the same directory
# as this notebook, in a subdirectory named "covid-19-data".

# load the full US rolling-average file, just to have a look

tbl_us = CSV.File("covid-19-data/rolling-averages/us.csv")

# ╔═╡ a7f2bdba-91ea-4225-8a82-a0d7bc2ed1d2
# We have no use for most of the columns in the data set, so now load a
# leaner selection.

tbl_us_cases = CSV.File("covid-19-data/rolling-averages/us.csv"; select=[:date, :cases, :cases_avg])

# ╔═╡ 0ad0cad5-1a38-485a-89af-3b9d4d8b1b08
# And as a further convenience, collect each of the three table columns into
# an array.

dates = collect(tbl_us_cases.date)

# ╔═╡ 7c75d0ce-58fc-4849-90db-8a66cefc2334
# Raw (unsmoothed) daily counts.

cases = collect(tbl_us_cases.cases)

# ╔═╡ 8f84f330-3158-4420-9e48-32133dfad8cc
# Seven-day trailing average

cases_7avg = collect(tbl_us_cases.cases_avg)

# ╔═╡ d7ff514d-cdce-432a-b046-92ccc3205cfe
# Sanity check: This should yield something that looks roughly like the
# Times Covid curve.

plot(dates, cases_7avg)

# ╔═╡ f3f76fab-a1f4-445f-ab84-4dae6630c2c2
# Also get a glimpse of the raw (i.e., unsmoothed) data

plot(dates, cases,
	seriestype = :sticks, linewidth = 0.4, linecolor=RGBA(0.8,0.1,0.1,0.5))

# ╔═╡ a5ec6c2c-5c68-42c8-8549-12c380adb0fb
# savefig(us_cases_avg_full_population, "us_cases_avg_full_population.svg")

# ╔═╡ f75da0f6-9a4f-461b-8503-4c36d89ac268
# Generate a list of first-day-of-the-month dates

monthlist = [Dates.value(d) for d in Date(2020,1,1):Dates.Month(1):Date(2021,7,1)]

# ╔═╡ e5f7be60-31a6-42db-8149-bf9a153d41c2
# Label the first of each month with a single-letter abbreviation

monthly_xticks = (monthlist, ["J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D", "J", "F", "M", "A", "M", "J", "J"])

# ╔═╡ 1ade71f3-2a7c-4a08-94de-4ab00c6eaea6
# Figure 1

us_cases =
plot(tbl_us_cases.date, [tbl_us_cases.cases, tbl_us_cases.cases_avg],
		size = (640, 400),
		seriestype = [:sticks :line],
		linewidth = [0.4 2.0],
		linecolor = [RGBA(0.2,0.1,0.9,0.4) "#4100eb"],
		xticks = monthly_xticks,
		yticks = ([0, 50000, 100000, 150000, 200000, 250000, 300000], ["0", "50", "100", "150", "200", "250", "300"]),
		xlabel = "       2020                                                     2021",
		ylabel = "daily new Covid cases (thousands)",
		guidefontsize = 9,
		titlefontsize = 9,
		xlims = [Dates.value(Date(2020,1,1)), Dates.value(Date(2021,7,31))],
		ylims = [-20000, 300000],
		label = ["daily reports" "7-day rolling average"],
		legend = (0.2, 0.7))

# ╔═╡ 1f058d22-aa90-4cb2-9719-3cdfbb10bb9e
# Figure 8

nov_jan_details_graph =
plot(tbl_us_cases.date[305:366],
	[Array(tbl_us_cases.cases)[305:366],
		Array(tbl_us_cases.cases_avg)[305:366]],
	size = [640, 400],
	seriestype = [:bar :line],
	linewidth = [0.4 3.0],
	linecolor = [:black "#4100eb"],
	fillcolor = RGBA(0.2,0.1,0.9,0.4),
	fontfamily = "helvetica",
	ylabel = "daily new Covid cases",
	xticks = ([	Date(2020,11,20),
				Date(2020,11,26),
				Date(2020,12,25),
				Date(2021,1,1),
				Date(2021,1,20)],
			  [	Dates.format(Date(2020,11,20), "u d"),
				Dates.format(Date(2020,11,26), "u d"),
				Dates.format(Date(2020,12,25), "u d"),
				Dates.format(Date(2021,1,1), "u d"),
				Dates.format(Date(2021,1,20), "u d")]),
	yticks = ([0, 100000, 200000, 300000], ["0", "100,000", "200,000", "300,000"]),
	guidefontsize = 9,
	legend = false)

# ╔═╡ 6f2b3f2a-a580-4005-93e4-f39f90305ae6
# savefig(nov_jan_details_graph, "nov_jan_details_graph.svg")

# ╔═╡ cdde1708-43d3-4903-826c-1b9ef1af24f7
# Figure 12 -- The Covid curve plotted with a y scale that extends not
# to 300,000 but to 330,000,000.

us_cases_avg_full_population =
plot(dates, cases_7avg,
		size = (320, 200),
		linewidth = 2.0,
		linecolor = "#4100eb",
		xticks = monthly_xticks,
		yticks = ([0, 100000000, 200000000, 300000000], ["0", "100", "200", "300"]),
		xlabel = "     2020                    2021",
		ylabel = "millions",
		titlefontsize = 8,
		guidefontsize = 9,
		xlims = [Dates.value(Date(2020,1,1)), Dates.value(Date(2021,7,31))],
		ylims = [0, 330000000],
		label = "U.S. Covid cases,\n 7-day rolling average",
		legend = false)

# ╔═╡ a9b12d92-2cbd-4d8b-a7da-38b03a8bd242
md"""
#### Fourier analysis of the Covid curve

This is done with Julia bindings to the FFTW package ("Fastest Fourier Transform in the West"). After experimenting with plain FFT and the Hartley transform, I wound
up using the discrete cosine transform for the published figures.

"""

# ╔═╡ 43666ed8-122f-4a40-b413-27a386e3913a
# a compound sinusoid, sum of cos(x) and cos(3x)

coscurve = [cos(x) + cos(3x) for x in 0:(π/128):(4π)]

# ╔═╡ cf6fdd02-6a0a-48ce-8c91-a6a36bdd71d5
# Figure 3

coscurve_fig =
plot(coscurve,
	size = (640, 200),
	seriestype = :line,
	linewidth = 2.0,
	linecolor = "#4100eb",
	xticks = false,
	yticks = false,
	# xlabel = false,
	ylabel = false,
	guidefontsize = 9,
	titlefontsize = 9,
	label = false,
	legend = false)

# ╔═╡ 368ff9d5-f29d-4975-b119-bcb4c5b86150
# savefig(coscurve_fig, "cos(x)+cos(3x)curve.svg")

# ╔═╡ 77b8c0ab-6687-4238-8c3c-68d503eb16fd
# Compute the discrete cosine transform of our sinusoidal curve

dct_coscurve = FFTW.dct(coscurve)

# ╔═╡ 9ac69ef5-1c27-4eea-b58f-83e7172b708e
# Figure 4: Graph the first 21 dct coefs for coscurve

coscurve_fourier_coefs =
plot([0:20], dct_coscurve[1:21],
	size = (500, 200),
	bar_edges = false,
	seriestype = :bar,
	fillcolor = "#4100eb",
	yticks = [],
	xticks = (0:1:20, ["0", "", "1", "", "2", "", "3", "", "4", "", "5", "", "6", "", "7", "", "8", "", "9", "", "10"]),
	xlabel = "wavenumber",
	fontfamily = "helvetica",
	guidefontsize = 9,
	titlefontsize = 9,
	legend = false)

# ╔═╡ 5713564a-fc5c-4e6d-a4cf-c137dcf6cffc
# savefig(coscurve_fourier_coefs, "DCT_coefficiencts_of_coscurve.svg")

# ╔═╡ 1ad022ea-7eed-4f10-b18a-172692e39921
# Create a new vector of the same length as F (the vector of Fourier coefs),
# with zeros in all entries except those specified by 'indices'.

function selectFourierCoefs(F, indices)
	cF = zeros(length(F))
	for i in indices
		cF[i] = F[i]
	end
	return cF
end

# ╔═╡ b6586485-cd85-45c4-9a55-3e6b7050006d
# Figure not published. Just to check that the two selected coefs (at
# wavenumbers 2 and 6) suffice to reconstruct the input curve. Note that
# the coefs are numbers 5 and 13 in the Julia array numbering.

reconstructed_coscurve =
plot([1:length(coscurve)], FFTW.idct(selectFourierCoefs(dct_coscurve, [5, 13])),
	size = (640, 200))

# ╔═╡ c6e8a34b-2904-45cd-acf9-e73cab338e3a
# Now apply the discrete cosine transform to the Covid data.

dct_cases_7avg = FFTW.dct(cases_7avg)

# ╔═╡ 5d526be1-3d3e-4d30-8e7d-ee8a8743cd3a
# Figure 5; Graph the first 21 coefs of Covid data DCT

dct_coef_graph_for_cases_7avg =
plot([0:20], 1/547 * dct_cases_7avg[1:21],
	seriestype = :bar,
	size = [500, 250],
	fillcolor = "#4100eb",
	yticks = [-1000, 0, 1000, 2000],
	xticks = (0:1:20, ["0", "", "1", "", "2", "", "3", "", "4", "", "5", "", "6", "", "7", "", "8", "", "9", "", "10"]),
	xlabel = "wavenumber",
	guidefontsize = 9,
	titlefontsize = 9,
	legend = false)

# ╔═╡ fdbd5aee-188e-44a7-bedb-375e8dbafb83
# savefig(dct_coef_graph_for_cases_7avg, "dct_coef_graph_for_cases_7avg.svg")

# ╔═╡ b3f318d5-83ec-44bf-907f-58706b386018
# Run the inverse DCT to reconstruct the original curve
# from selected Fourier coefs

function idct_curve(coefs, selectors)
	return FFTW.idct(selectFourierCoefs(coefs, selectors))
end

# ╔═╡ 740e3274-a134-4c17-9a0c-1256fb16464f
# Figure 6 -- Plot each of the first 13 component curves from the Covid DCT

fourier_components_from_dct_cases_7avg =
plot(dates, [idct_curve(dct_cases_7avg, [n]) for n in 1:13],
		seriestype = :line,
		size = (500, 400),
		linewidth = 2.0,
		linecolor = "#4100eb",
		fontfamily = "helvetica",
		xticks = monthly_xticks,
		yticks = false,
		xlabel = "     2020                                       2021",
		guidefontsize = 10,
		titlefontsize = 10,
		xlims = [Dates.value(Date(2020,1,1)), Dates.value(Date(2021,7,31))],
		legend = false)

# ╔═╡ 8e7af716-ce34-46b4-97ad-6d9d19236ae7
# savefig(fourier_components_from_dct_cases_7avg, "fourier_components_from_dct_cases_7avg.svg")

# ╔═╡ 72b56571-dcae-4998-bad9-b22bf64f42c9
# Figure 7 -- Plot the sum of the 13 DCT components against a background
# of the 7avg curve

recon_from_dct_cases_7avg =
plot(tbl_us_cases.date, [cases_7avg,
		idct_curve(dct_cases_7avg, [1:13])],
		seriestype = [:line :line],
		fillcolor = [RGBA(0.8,0.1,0.1,0.2) nothing],
		fillrange = 0,
		size = (640, 400),
		linewidth = [1.5 2.5],
		linecolor = [RGBA(0.8,0.1,0.1,0.2) "#4100eb"],
		fontfamily = "helvetica",
		xticks = monthly_xticks,
		yticks = ([0, 50000, 100000, 150000, 200000, 250000], ["0", "50,000", "100,000", "150,000", "200,000", "250,000"]),
		xlabel = "     2020                                       2021",
		ylabel = "daily new Covid cases",
		guidefontsize = 10,
		titlefontsize = 10,
		xlims = [Dates.value(Date(2020,1,1)), Dates.value(Date(2021,7,31))],
		label = ["7-day average case count" "first 13 Fourier coefficients"],
		legend = (0.15, 0.7))

# ╔═╡ 7cc40b50-80cc-496f-96ac-4c6fb502e638
# savefig(recon_from_dct_cases_7avg, "reconstructed_covid_curve_from_14_coefs.svg")

# ╔═╡ 3fd9fec9-80d2-43ee-98c5-1392a9f77c65
md"""
#### Time derivative of the Covid data
"""

# ╔═╡ 4d76574b-df99-4e37-ba4b-d8a39f60efe1
first_diffs_cases = diff(cases)

# ╔═╡ a706728a-b042-46bf-b6a6-eb865427ab71
first_diffs_cases_7avg = diff(cases_7avg)

# ╔═╡ 74d23fa6-71c0-4200-a8a2-6eb5b693b33b
first_diffs_idct_curve_1_13 = diff(idct_curve(dct_cases_7avg, [1:13]))

# ╔═╡ df4639f4-e766-40c2-b8ea-d878032ff44d
# Figure 9 -- Plots of time derivatives

first_derivatives_plot =
plot(dates[1:546],
	[first_diffs_cases,
		180000 .+ (2 .* first_diffs_cases_7avg),
		240000 .+ (4 .* first_diffs_idct_curve_1_13)],
	size = [640, 550],
	xlims = [Dates.value(Date(2020,1,1)), Dates.value(Date(2021,7,31))],
	yticks = ([0, 180000, 240000], ["0", "0", "0"]),
	xticks = monthly_xticks,
	xlabel = "       2020                                                     2021",
	fontfamily = "helvetica",
	linecolor = ["#4100eb" "#dd4411" :green],
	linewidth = 1.5,
	annotationfontsize = 10,
	annotationhalign = :left,
	annotations =
		[(Dates.value(Date(2020,1,1)), 17000, "time derivative of raw data"),
		 (Dates.value(Date(2020,1,1)), 195000, "time derivative of averaged data (2x)"),
		 (Dates.value(Date(2020,1,1)), 255000, "time derivative of FFT reconstructed data (4x)")],
	legend = false
	)

# ╔═╡ 8556c18b-169a-4b28-bd4c-93baa60d80c3
# savefig(first_derivatives_plot, "first_derivatives.svg")

# ╔═╡ fb525a7e-3b0b-4886-bffe-ee28110ec8f9
# Count the zero crossings in a vector. Note that we count only
# actual crossings; merely touching zero doesn't count

function count_zero_crossings(data)
	crossings = map(sign, data)
	thecount = 0
	i = 1
	while crossings[i] == 0
		i += 1
	end
	previous = crossings[i]
	for j in i+1:length(crossings)
		if crossings[j] == -previous
			thecount += 1
			previous = -previous
		end
	end
	return thecount
end

# ╔═╡ 26db9d45-8d83-4d29-b89a-6ae0c9dcf4df
# Get counts for raw data, 7-day smoothed, and the DCT reconstruction

map(count_zero_crossings, [first_diffs_cases, first_diffs_cases_7avg, first_diffs_idct_curve_1_13])

# ╔═╡ d2634e23-e155-4614-aaa9-55b0bfe90c89
# Figure 10

exp_log_plot =
begin
	xs = 0:0.05:20
	plot(xs,
		[[0.8^x for x in xs],
		 [1.0^x for x in xs],
		 [1.25^x for x in xs],
		 [2.0^x for x in xs],
		 [3.0^x for x in xs]],
		ylims = [0.01, 1e10],
		yscale = :log10,
		yticks = [0.01, 1.0, 1e2, 1e4, 1e6, 1e8, 1e10],
		linewidth = 2.0,
		xlabel = "serial passage times (SPT)",
		ylabel = "number of new cases",
		fontfamily = "helvetica",
		legend = false)
end

# ╔═╡ e78cb3a9-7319-459e-b772-401f96f89cfb
# savefig(exp_log_plot, "exp_log_plot.svg")

# ╔═╡ 6dca6932-d9fc-435a-972b-7c8caa797512
# In the text I talk about summing many sinusoidal waves that
# have the same frequency f but random phases: The sum is always
# anothr sinusoid of frequency f. I could not make that statement
# without checking it computationally.

function sinesum(N, xmax)
	samples = 400
	sumcurve = zeros(samples)
	for i in 1:N
		phase = rand() * 2π
		for x in 1:samples
			sumcurve[x] += sin(((x * xmax) / samples) + phase)
		end
	end
	return sumcurve
end

# ╔═╡ 1aab13b2-3e7b-4fc3-9287-1676318501bd
# The proof of the pudding. For any N and any xmax, the result
# will be a sinusoidal curve

plot(sinesum(6, 10π))

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

[compat]
CSV = "~0.8.5"
FFTW = "~1.4.3"
Plots = "~1.20.0"
"""

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

[[AbstractFFTs]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "485ee0867925449198280d4af84bdb46a2a404d0"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.0.1"

[[Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

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

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

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

[[Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c3598e525718abcc440f69cc6d5f60dda0a1b61e"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.6+5"

[[CSV]]
deps = ["Dates", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode"]
git-tree-sha1 = "b83aa3f513be680454437a0eee21001607e5d983"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.8.5"

[[Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "e2f47f6d8337369411569fd45ae5753ca10394c6"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.0+6"

[[ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random", "StaticArrays"]
git-tree-sha1 = "ed268efe58512df8c7e224d2e170afd76dd6a417"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.13.0"

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

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

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

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[Contour]]
deps = ["StaticArrays"]
git-tree-sha1 = "9f02045d934dc030edad45944ea80dbd1f0ebea7"
uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
version = "0.5.7"

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

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

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

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

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

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

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

[[EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "92d8f9f208637e8d2d28c664051a00569c01493d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.1.5+1"

[[Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b3bfd02e98aedfa5cf885665493c5598c350cd2f"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.2.10+0"

[[FFMPEG]]
deps = ["FFMPEG_jll"]
git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
version = "0.4.1"

[[FFMPEG_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "LibVPX_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
git-tree-sha1 = "3cc57ad0a213808473eafef4845a74766242e05f"
uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
version = "4.3.1+4"

[[FFTW]]
deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
git-tree-sha1 = "f985af3b9f4e278b1d24434cbb546d6092fca661"
uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
version = "1.4.3"

[[FFTW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3676abafff7e4ff07bbd2c42b3d8201f31653dcc"
uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
version = "3.3.9+8"

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

[[Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "35895cf184ceaab11fd778b4590144034a167a2f"
uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
version = "2.13.1+14"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "cbd58c9deb1d304f5a245a0b7eb841a2560cfec6"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.10.1+5"

[[FriBidi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91"
uuid = "559328eb-81f9-559d-9380-de523a88c83c"
version = "1.0.10+0"

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

[[GLFW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"]
git-tree-sha1 = "dba1e8614e98949abfa60480b13653813d8f0157"
uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
version = "3.3.5+0"

[[GR]]
deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"]
git-tree-sha1 = "182da592436e287758ded5be6e32c406de3a2e47"
uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
version = "0.58.1"

[[GR_jll]]
deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "d59e8320c2747553788e4fc42231489cc602fa50"
uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
version = "0.58.1+0"

[[GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "58bcdf5ebc057b085e58d95c138725628dd7453c"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.4.1"

[[Gettext_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
version = "0.21.0+0"

[[Glib_jll]]
deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "7bf67e9a481712b3dbe9cb3dac852dc4b1162e02"
uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
version = "2.68.3+0"

[[Grisu]]
git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
version = "1.0.2"

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

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

[[IntelOpenMP_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d979e54b71da82f3a65b62553da4fc3d18c9004c"
uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
version = "2018.0.3+2"

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

[[IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

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

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

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

[[JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d735490ac75c5cb9f1b00d8b5509c11984dc6943"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.0+0"

[[LAME_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
version = "3.100.1+0"

[[LZO_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
version = "2.10.1+0"

[[LaTeXStrings]]
git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.2.1"

[[Latexify]]
deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
git-tree-sha1 = "a4b12a1bd2ebade87891ab7e36fdbce582301a92"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.15.6"

[[LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

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

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

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

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

[[LibVPX_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "12ee7e23fa4d18361e7c2cde8f8337d4c3101bc7"
uuid = "dd192d2f-8180-539f-9fb4-cc70b1dcf69a"
version = "1.10.0+0"

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

[[Libffi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "761a393aeccd6aa92ec3515e428c26bf99575b3b"
uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
version = "3.2.2+0"

[[Libgcrypt_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae"
uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
version = "1.8.7+0"

[[Libglvnd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"]
git-tree-sha1 = "7739f837d6447403596a75d19ed01fd08d6f56bf"
uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29"
version = "1.3.0+3"

[[Libgpg_error_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
version = "1.42.0+0"

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[Libmount_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
version = "2.35.0+0"

[[Libtiff_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "340e257aada13f95f98ee352d316c3bed37c8ab9"
uuid = "89763e89-9b03-5906-acba-b20f662cd828"
version = "4.3.0+0"

[[Libuuid_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066"
uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
version = "2.36.0+0"

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

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

[[MKL_jll]]
deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
git-tree-sha1 = "c253236b0ed414624b083e6b72bfe891fbd2c7af"
uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
version = "2021.1.1+1"

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

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

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

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

[[Measures]]
git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f"
uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
version = "0.3.1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "4ea90bd5d3985ae1f9a908bd4500ae88921c5ce7"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.0"

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

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

[[NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

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

[[Ogg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7937eda4681660b4d6aeeecc2f7e1c81c8ee4e2f"
uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
version = "1.3.5+0"

[[OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "15003dcb7d8db3c6c857fda14891a539a8f2705a"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.10+0"

[[Opus_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720"
uuid = "91d4177d-7536-5919-b921-800302f37372"
version = "1.3.2+0"

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

[[PCRE_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b2a7af664e098055a7529ad1a900ded962bca488"
uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc"
version = "8.44.0+0"

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

[[Pixman_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29"
uuid = "30392449-352a-5448-841d-b1acce4e97dc"
version = "0.40.1+0"

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

[[PlotThemes]]
deps = ["PlotUtils", "Requires", "Statistics"]
git-tree-sha1 = "a3a964ce9dc7898193536002a6dd892b1b5a6f1d"
uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
version = "2.0.1"

[[PlotUtils]]
deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
git-tree-sha1 = "501c20a63a34ac1d015d5304da0e645f42d91c9f"
uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
version = "1.0.11"

[[Plots]]
deps = ["Base64", "Contour", "Dates", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs"]
git-tree-sha1 = "e39bea10478c6aff5495ab522517fae5134b40e3"
uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
version = "1.20.0"

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

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

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

[[Qt5Base_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"]
git-tree-sha1 = "ad368663a5e20dbb8d6dc2fddeefe4dae0781ae8"
uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1"
version = "5.15.3+0"

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

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

[[RecipesBase]]
git-tree-sha1 = "b3fb709f3c97bfc6e948be68beeecb55a0b340ae"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.1.1"

[[RecipesPipeline]]
deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"]
git-tree-sha1 = "2a7a2469ed5d94a98dea0e85c46fa653d76be0cd"
uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
version = "0.3.4"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

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

[[Scratch]]
deps = ["Dates"]
git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda"
uuid = "6c6a2e73-6563-6170-7368-637461726353"
version = "1.1.0"

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

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

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

[[Showoff]]
deps = ["Dates", "Grisu"]
git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de"
uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
version = "1.0.3"

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

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

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

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "b28f39450421d07d89ab5d126fd15e5246350e8a"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.11"

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

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "fed1ec1e65749c4d96fc20dd13bea72b55457e62"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.9"

[[StructArrays]]
deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
git-tree-sha1 = "000e168f5cc9aded17b6999a560b7c11dda69095"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.0"

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

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

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

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

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

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

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

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

[[Wayland_jll]]
deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "3e61f0b86f90dacb0bc0e73a0c5a83f6a8636e23"
uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89"
version = "1.19.0+0"

[[Wayland_protocols_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll"]
git-tree-sha1 = "2839f1c1296940218e35df0bbb220f2a79686670"
uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91"
version = "1.18.0+4"

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[XSLT_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
version = "1.1.34+0"

[[Xorg_libX11_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527"
uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
version = "1.6.9+4"

[[Xorg_libXau_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e"
uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
version = "1.0.9+4"

[[Xorg_libXcursor_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "12e0eb3bc634fa2080c1c37fccf56f7c22989afd"
uuid = "935fb764-8cf2-53bf-bb30-45bb1f8bf724"
version = "1.2.0+4"

[[Xorg_libXdmcp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4"
uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
version = "1.1.3+4"

[[Xorg_libXext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
version = "1.3.4+4"

[[Xorg_libXfixes_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "0e0dc7431e7a0587559f9294aeec269471c991a4"
uuid = "d091e8ba-531a-589c-9de9-94069b037ed8"
version = "5.0.3+4"

[[Xorg_libXi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXfixes_jll"]
git-tree-sha1 = "89b52bc2160aadc84d707093930ef0bffa641246"
uuid = "a51aa0fd-4e3c-5386-b890-e753decda492"
version = "1.7.10+4"

[[Xorg_libXinerama_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll"]
git-tree-sha1 = "26be8b1c342929259317d8b9f7b53bf2bb73b123"
uuid = "d1454406-59df-5ea1-beac-c340f2130bc3"
version = "1.1.4+4"

[[Xorg_libXrandr_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "34cea83cb726fb58f325887bf0612c6b3fb17631"
uuid = "ec84b674-ba8e-5d96-8ba1-2a689ba10484"
version = "1.5.2+4"

[[Xorg_libXrender_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
version = "0.9.10+4"

[[Xorg_libpthread_stubs_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb"
uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
version = "0.1.0+3"

[[Xorg_libxcb_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6"
uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
version = "1.13.0+3"

[[Xorg_libxkbfile_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "926af861744212db0eb001d9e40b5d16292080b2"
uuid = "cc61e674-0454-545c-8b26-ed2c68acab7a"
version = "1.1.0+4"

[[Xorg_xcb_util_image_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "0fab0a40349ba1cba2c1da699243396ff8e94b97"
uuid = "12413925-8142-5f55-bb0e-6d7ca50bb09b"
version = "0.4.0+1"

[[Xorg_xcb_util_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll"]
git-tree-sha1 = "e7fd7b2881fa2eaa72717420894d3938177862d1"
uuid = "2def613f-5ad1-5310-b15b-b15d46f528f5"
version = "0.4.0+1"

[[Xorg_xcb_util_keysyms_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "d1151e2c45a544f32441a567d1690e701ec89b00"
uuid = "975044d2-76e6-5fbe-bf08-97ce7c6574c7"
version = "0.4.0+1"

[[Xorg_xcb_util_renderutil_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "dfd7a8f38d4613b6a575253b3174dd991ca6183e"
uuid = "0d47668e-0667-5a69-a72c-f761630bfb7e"
version = "0.3.9+1"

[[Xorg_xcb_util_wm_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "e78d10aab01a4a154142c5006ed44fd9e8e31b67"
uuid = "c22f9ab0-d5fe-5066-847c-f4bb1cd4e361"
version = "0.4.1+1"

[[Xorg_xkbcomp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxkbfile_jll"]
git-tree-sha1 = "4bcbf660f6c2e714f87e960a171b119d06ee163b"
uuid = "35661453-b289-5fab-8a00-3d9160c6a3a4"
version = "1.4.2+4"

[[Xorg_xkeyboard_config_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xkbcomp_jll"]
git-tree-sha1 = "5c8424f8a67c3f2209646d4425f3d415fee5931d"
uuid = "33bec58e-1273-512f-9401-5d533626f822"
version = "2.27.0+4"

[[Xorg_xtrans_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845"
uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
version = "1.4.0+3"

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

[[Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.0+0"

[[libass_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "acc685bcf777b2202a904cdcb49ad34c2fa1880c"
uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
version = "0.14.0+4"

[[libfdk_aac_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7a5780a0d9c6864184b3a2eeeb833a0c871f00ab"
uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
version = "0.1.6+4"

[[libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[libvorbis_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
git-tree-sha1 = "c45f4e40e7aafe9d086379e5578947ec8b95a8fb"
uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
version = "1.3.7+0"

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

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

[[x264_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d713c1ce4deac133e3334ee12f4adff07f81778f"
uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
version = "2020.7.14+2"

[[x265_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "487da2f8f2f0c8ee0e83f39d13037d6bbf0a45ab"
uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
version = "3.0.0+3"

[[xkbcommon_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"]
git-tree-sha1 = "ece2350174195bb31de1a63bea3a41ae1aa593b6"
uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd"
version = "0.9.1+5"
"""

# ╔═╡ Cell order:
# ╠═12032637-059c-4227-840a-55c32a5c402c
# ╠═f502a2cc-ea34-11eb-0eb3-93d1cb0a9ad2
# ╠═e78d0da4-69b2-450a-b9c3-a7087a698618
# ╠═65c4e11f-7c12-46ef-8c9c-e97d139a80d5
# ╠═1bd45843-5bdb-4f19-9f04-23d3917650db
# ╠═47bc3ba1-7d7c-435d-8066-87fb07baa293
# ╠═a7f2bdba-91ea-4225-8a82-a0d7bc2ed1d2
# ╠═0ad0cad5-1a38-485a-89af-3b9d4d8b1b08
# ╠═7c75d0ce-58fc-4849-90db-8a66cefc2334
# ╠═8f84f330-3158-4420-9e48-32133dfad8cc
# ╠═d7ff514d-cdce-432a-b046-92ccc3205cfe
# ╠═f3f76fab-a1f4-445f-ab84-4dae6630c2c2
# ╠═1ade71f3-2a7c-4a08-94de-4ab00c6eaea6
# ╠═a5ec6c2c-5c68-42c8-8549-12c380adb0fb
# ╠═f75da0f6-9a4f-461b-8503-4c36d89ac268
# ╠═e5f7be60-31a6-42db-8149-bf9a153d41c2
# ╠═1f058d22-aa90-4cb2-9719-3cdfbb10bb9e
# ╠═6f2b3f2a-a580-4005-93e4-f39f90305ae6
# ╠═cdde1708-43d3-4903-826c-1b9ef1af24f7
# ╟─a9b12d92-2cbd-4d8b-a7da-38b03a8bd242
# ╠═43666ed8-122f-4a40-b413-27a386e3913a
# ╠═cf6fdd02-6a0a-48ce-8c91-a6a36bdd71d5
# ╠═368ff9d5-f29d-4975-b119-bcb4c5b86150
# ╠═77b8c0ab-6687-4238-8c3c-68d503eb16fd
# ╠═9ac69ef5-1c27-4eea-b58f-83e7172b708e
# ╠═5713564a-fc5c-4e6d-a4cf-c137dcf6cffc
# ╠═b6586485-cd85-45c4-9a55-3e6b7050006d
# ╠═1ad022ea-7eed-4f10-b18a-172692e39921
# ╠═c6e8a34b-2904-45cd-acf9-e73cab338e3a
# ╠═5d526be1-3d3e-4d30-8e7d-ee8a8743cd3a
# ╠═fdbd5aee-188e-44a7-bedb-375e8dbafb83
# ╠═b3f318d5-83ec-44bf-907f-58706b386018
# ╠═740e3274-a134-4c17-9a0c-1256fb16464f
# ╠═8e7af716-ce34-46b4-97ad-6d9d19236ae7
# ╠═72b56571-dcae-4998-bad9-b22bf64f42c9
# ╠═7cc40b50-80cc-496f-96ac-4c6fb502e638
# ╟─3fd9fec9-80d2-43ee-98c5-1392a9f77c65
# ╠═4d76574b-df99-4e37-ba4b-d8a39f60efe1
# ╠═a706728a-b042-46bf-b6a6-eb865427ab71
# ╠═74d23fa6-71c0-4200-a8a2-6eb5b693b33b
# ╠═df4639f4-e766-40c2-b8ea-d878032ff44d
# ╠═8556c18b-169a-4b28-bd4c-93baa60d80c3
# ╠═fb525a7e-3b0b-4886-bffe-ee28110ec8f9
# ╠═26db9d45-8d83-4d29-b89a-6ae0c9dcf4df
# ╠═d2634e23-e155-4614-aaa9-55b0bfe90c89
# ╠═e78cb3a9-7319-459e-b772-401f96f89cfb
# ╠═6dca6932-d9fc-435a-972b-7c8caa797512
# ╠═1aab13b2-3e7b-4fc3-9287-1676318501bd
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl

max_stars_repo_name: bit-player/Covidcurve

max_issues_repo_name: bit-player/Covidcurve

max_forks_repo_name: bit-player/Covidcurve


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

using Markdown
using InteractiveUtils

# ╔═╡ 1a6379ce-0ffc-11eb-0fce-853d047f252a
begin
	using PyPlot
	using PyCall #I'll gonna use SpecialFunctions.jl and HypergeometricFunctions.jl
	using LaTeXStrings
end

# ╔═╡ cd516e2c-0ffe-11eb-2415-832b4823d2ca
md"### Legacy notebook using PyCall and mpmath library

I will keep this record for future reference"

# ╔═╡ 3e643f7a-0ffc-11eb-021f-636a88842643
mp = pyimport("mpmath")

# ╔═╡ 5ccf83de-0ffc-11eb-2b02-8f1ff10f54af
begin
	rc("text", usetex=true)
	ion()
end

# ╔═╡ 61e49ab2-0ffc-11eb-157c-db7651aee125
begin
	N = 32
	dim = N+1
	const p = 1/2.0
end;

# ╔═╡ 67fd8670-0ffc-11eb-2b6e-ffefab121f31
function nbinom(n,k); # Binomial coefficient in terms of the Gamma function
    return 1/((n+1)*convert(Float64,mp.beta(k+1, n-k+1)))
end;

# ╔═╡ 6ef5f2ac-0ffc-11eb-2dcc-39e1895c1b70
function krav(i,j) # Kravchuk symmetric polynomial in terms of 2F1
    return convert(Float64,mp.hyp2f1(-j,-i,-N,1/p,zeroprec = 250))
end;

# ╔═╡ 74c15e90-0ffc-11eb-168c-afd38ba93c9b
function omega(j)
    return nbinom(N,j)*(p^j)*(1-p)^(N-j)
end;

# ╔═╡ 7a746828-0ffc-11eb-2c9e-87213ad35743
function h(i)
    return (((1-p)/p)^i)/(nbinom(N,i))
end;

# ╔═╡ 80054320-0ffc-11eb-1837-b7f067634a26
function K(i,j)
    return sqrt(omega(j)/h(i))*krav(i,j)
end;

# ╔═╡ 8804e044-0ffc-11eb-3cf9-efd4e3a8585d
begin
	KT = zeros(Float64,(dim,dim))
	@timed for l in 0:N
    	for m in 0:N
        	KT[l+1,m+1] = K(l,m)
    	end
	end
end

# ╔═╡ 8ed657c0-0ffc-11eb-1974-3154511ac71a
begin
	fig2 = figure(figsize=(5,5))
	ax2 = gca()
	ax2.imshow(KT, cmap="gray")
	title(L"\textbf{Kravchuk Matrix}\; K(i,j)\;\textbf{(The trilobite)}")
	ax2.set_xlabel(L"\textbf{Pseudo Energy Level:}\; 0\leq i\leq N")
	ax2.set_ylabel(L"\textbf{Position Index:}\; 0\leq j\leq N")
	ax2.tick_params(direction="out",length=5,width=2,labelsize=10)
	ax2.set_xticks(0:8:N, minor = false)
	ax2.set_yticks(0:8:N, minor = false)

	axins1 = ax2.inset_axes([0, 0.8, 0.2, 0.2])
	axins1.plot(KT[1,:],"k:",linewidth=1)
	axins1.set_xticklabels("")
	axins1.set_yticklabels("")
	axins1.text(7.5,0.01,L"K(0,j)",fontsize=9)

	axins2 = ax2.inset_axes([0.8, 0, 0.2, 0.2])
	axins2.plot(KT[N+1,:],"k:",linewidth=1)
	axins2.set_xticklabels("")
	axins2.set_yticklabels("")
	axins2.text(15,-0.35,L"K(N,j)",fontsize=9)

	tight_layout()
end

# ╔═╡ 36310bca-0ffd-11eb-1acf-dfc63763f5df
md"Using PyCall with mpmath library, we obtain, for N=32, a Kravchuk matrix of NxN in 0.673462 s. With the native libraries, Special and Hypergeometric, this same matrix is obtained in 0.0028436 s (see FrFK.jl). he native evaluation is around 237 times faster. 😲"

# ╔═╡ Cell order:
# ╠═cd516e2c-0ffe-11eb-2415-832b4823d2ca
# ╠═1a6379ce-0ffc-11eb-0fce-853d047f252a
# ╟─3e643f7a-0ffc-11eb-021f-636a88842643
# ╠═5ccf83de-0ffc-11eb-2b02-8f1ff10f54af
# ╠═61e49ab2-0ffc-11eb-157c-db7651aee125
# ╠═67fd8670-0ffc-11eb-2b6e-ffefab121f31
# ╠═6ef5f2ac-0ffc-11eb-2dcc-39e1895c1b70
# ╠═74c15e90-0ffc-11eb-168c-afd38ba93c9b
# ╠═7a746828-0ffc-11eb-2c9e-87213ad35743
# ╠═80054320-0ffc-11eb-1837-b7f067634a26
# ╠═8804e044-0ffc-11eb-3cf9-efd4e3a8585d
# ╠═8ed657c0-0ffc-11eb-1974-3154511ac71a
# ╟─36310bca-0ffd-11eb-1acf-dfc63763f5df


text:
module VirtualBrain

include("datasets.jl")
include("mpr.jl")
include("kuramoto.jl")
include("epileptor.jl")

# TODO surface & field modules
# TODO clean up & unify approaches
# TODO include gradient calc examples

end


text:
"""
    log_pseudolikelihood(rbm, v; exact = false)

Log-pseudolikelihood of `v`. If `exact` is `true`, the exact pseudolikelihood is returned.
But this is slow if `v` consists of many samples. Therefore by default `exact` is `false`,
in which case the result is a stochastic approximation, where a random site is selected
for each sample, and its conditional probability is calculated. In average the results
with `exact = false` coincide with the deterministic result, and the estimate is more
precise as the number of samples increases.
"""
function log_pseudolikelihood(rbm::RBM, v::AbstractArray; exact::Bool=false)
    if exact
        return log_pseudolikelihood_exact(rbm, v)
    else
        return log_pseudolikelihood_stoch(rbm, v)
    end
end

"""
    log_pseudolikelihood_stoch(rbm, v)

Log-pseudolikelihood of `v`. This function computes an stochastic approximation, by doing
a trace over random sites for each sample. For large number of samples, this is in average
close to the exact value of the pseudolikelihood.
"""
function log_pseudolikelihood_stoch(rbm::RBM, v::AbstractArray)
    @assert size(visible(rbm)) == size(v)[1:ndims(visible(rbm))]
    sites = [
        rand(CartesianIndices(sitesize(visible(rbm))))
        for _ in CartesianIndices(batch_size(visible(rbm), v))
    ]
    return log_pseudolikelihood_sites(rbm, v, sites)
end

"""
    log_pseudolikelihood_sites(rbm, v, sites)

Log-pseudolikelihood of a site conditioned on the other sites, where `sites`
is an array of site indices (CartesianIndex), one for each sample.
Returns an array of log-pseudolikelihood values, for each sample.
"""
function log_pseudolikelihood_sites(
    rbm::RBM, v::AbstractArray, sites::AbstractArray{<:CartesianIndex}
)
    @assert size(visible(rbm)) == size(v)[1:ndims(visible(rbm))]
    @assert size(sites) == batch_size(visible(rbm), v)
    ΔE = substitution_matrix_sites(rbm, v, sites)
    @assert size(ΔE) == (colors(visible(rbm)), batch_size(visible(rbm), v)...)
    lPL = -logsumexp(-ΔE; dims=1)
    @assert size(lPL) == (1, batch_size(visible(rbm), v)...)
    return reshape(lPL, batch_size(visible(rbm), v))
end

"""
    log_pseudolikelihood_exact(rbm, v)

Log-pseudolikelihood of `v`. This function computes the exact pseudolikelihood, doing
traces over all sites. Note that this can be slow for large number of samples.
"""
function log_pseudolikelihood_exact(rbm::RBM, v::AbstractArray)
    @assert size(visible(rbm)) == size(v)[1:ndims(visible(rbm))]
    ΔE = substitution_matrix_exhaustive(rbm, v)
    @assert size(ΔE) == (
        colors(visible(rbm)), sitesize(visible(rbm))..., batch_size(visible(rbm), v)...
    )
    lPLsites = -logsumexp(-ΔE; dims=1)
    @assert size(lPLsites) == (1, sitesize(visible(rbm))..., batch_size(visible(rbm), v)...)
    lPL = mean(lPLsites; dims=2:(sitedims(visible(rbm)) + 1))
    return reshape(lPL, batch_size(visible(rbm), v))
end

"""
    substitution_matrix_sites(rbm, v, sites)

Returns an q x B matrix of free energies `F`, where `q` is the number of possible values
of each site, and `B` the number of data points. The entry `F[x,b]` equals the free energy
cost of flipping `site[b]` of `v[b]` to `x`, that is (schemetically):

    F[x, b] = free_energy(rbm, v_) - free_energy(rbm, v)

where `v = v[b]`, and `v_` is the same as `v` in all sites except `site[b]`,
where `v_` has the value `x`.
"""
function substitution_matrix_sites end

function substitution_matrix_sites(
    rbm::RBM{<:Binary}, v::AbstractArray, sites::AbstractArray{<:CartesianIndex}
)
    @assert size(visible(rbm)) == size(v)[1:ndims(visible(rbm))]
    @assert size(sites) == batch_size(visible(rbm), v)
    E_ = zeros(2, batch_size(visible(rbm), v)...)
    for (k, x) in enumerate((false, true))
        v_ = copy(v)
        for (b, i) in pairs(sites)
            v_[i, b] = x
        end
        E_[k, ..] .= free_energy(rbm, v_)
    end
    E = [E_[(v[i, b] > 0) + 1, b] for (b, i) in pairs(sites)]
    return E_ .- reshape(E, 1, batch_size(visible(rbm), v)...)
end

function substitution_matrix_sites(
    rbm::RBM{<:Spin}, v::AbstractArray, sites::AbstractArray{<:CartesianIndex}
)
    @assert size(visible(rbm)) == size(v)[1:ndims(visible(rbm))]
    @assert size(sites) == batch_size(visible(rbm), v)
    E_ = zeros(2, batch_size(visible(rbm), v)...)
    for (k, x) in enumerate((-1, 1))
        v_ = copy(v)
        for (b, i) in pairs(sites)
            v_[i, b] = x
        end
        E_[k, ..] .= free_energy(rbm, v_)
    end
    E = [E_[(v[i, b] > 0) + 1, b] for (b, i) in pairs(sites)]
    return E_ .- reshape(E, 1, batch_size(visible(rbm), v)...)
end

function substitution_matrix_sites(
    rbm::RBM{<:Potts}, v::AbstractArray, sites::AbstractArray{<:CartesianIndex}
)
    @assert size(visible(rbm)) == size(v)[1:ndims(visible(rbm))]
    @assert size(sites) == batch_size(visible(rbm), v)
    E_ = zeros(size(rbm.visible, 1), batch_size(visible(rbm), v)...)
    for k in 1:size(rbm.visible, 1)
        v_ = copy(v)
        for (b, i) in pairs(sites)
            v_[:, i, b] .= false
            v_[k, i, b] = true
        end
        E_[k, ..] .= free_energy(rbm, v_)
    end
    c = onehot_decode(v)
    E = [E_[c[i, b], b] for (b, i) in pairs(sites)]
    return E_ .- reshape(E, 1, batch_size(visible(rbm), v)...)
end

"""
    substitution_matrix_exhaustive(rbm, v)

Returns an q x N x B tensor of free energies `F`, where `q` is the number of possible
values of each site, `B` the number of data points, and `N` the sequence length:

````
q, N, B = size(v)
```

Thus `F` and `v` have the same size.
The entry `F[x,i,b]` gives the free energy cost of flipping site `i` to `x`
of `v[b]` from its original value to `x`, that is:

    F[x,i,b] = free_energy(rbm, v_) - free_energy(rbm, v[b])

where `v_` is the same as `v[b]` in all sites but `i`, where `v_` has the value `x`.

Note that `i` can be a set of indices.
"""
function substitution_matrix_exhaustive end

function substitution_matrix_exhaustive(rbm::RBM{<:Binary}, v::AbstractArray)
    @assert size(visible(rbm)) == size(v)[1:ndims(visible(rbm))]
    batch_indices = CartesianIndices(batch_size(visible(rbm), v))
    E_ = zeros(2, size(v)...)
    for i in CartesianIndices(size(visible(rbm)))
        v_ = copy(v)
        for (k, x) in enumerate((false, true))
            v_[i, batch_indices] .= x
            E_[k, i, batch_indices] .= free_energy(rbm, v_)
        end
    end
    E = [E_[(v[k] > 0) + 1, k] for k in CartesianIndices(v)]
    return E_ .- reshape(E, 1, size(v)...)
end

function substitution_matrix_exhaustive(rbm::RBM{<:Spin}, v::AbstractArray)
    @assert size(visible(rbm)) == size(v)[1:ndims(visible(rbm))]
    batch_indices = CartesianIndices(batch_size(visible(rbm), v))
    E_ = zeros(2, size(v)...)
    for i in CartesianIndices(sitesize(visible(rbm)))
        v_ = copy(v)
        for (k, x) in enumerate((Int8(-1), Int8(1)))
            v_[i, batch_indices] .= x
            E_[k, i, batch_indices] .= free_energy(rbm, v_)
        end
    end
    E = [E_[(v[k] > 0) + 1, k] for k in CartesianIndices(v)]
    return E_ .- reshape(E, 1, size(v)...)
end

function substitution_matrix_exhaustive(rbm::RBM{<:Potts}, v::AbstractArray)
    @assert size(visible(rbm)) == size(v)[1:ndims(visible(rbm))]
    batch_indices = CartesianIndices(batch_size(visible(rbm), v))
    E_ = zeros(size(v))
    for i in CartesianIndices(sitesize(visible(rbm)))
        v_ = copy(v)
        for x in 1:size(rbm.visible, 1)
            v_[:, i, batch_indices] .= false
            v_[x, i, batch_indices] .= true
            E_[x, i, batch_indices] .= free_energy(rbm, v_)
        end
    end
    c = onehot_decode(v)
    E = [E_[c[k], k] for k in CartesianIndices(c)]
    return E_ .- reshape(E, 1, size(E)...)
end

#= ***
For Binary and Spin layers, a specialized log_pseudolikelihood_sites is a bit faster.
*** =#

function log_pseudolikelihood_sites(
    rbm::RBM{<:Binary}, v::AbstractArray, sites::AbstractArray{<:CartesianIndex}
)
    @assert size(visible(rbm)) == size(v)[1:ndims(visible(rbm))]
    @assert size(sites) == batch_size(visible(rbm), v)
    v_ = copy(v)
    for (b, i) in pairs(sites)
        v_[i, b] = 1 - v_[i, b]
    end
    F = free_energy(rbm, v)
    F_ = free_energy(rbm, v_)
    return -log1pexp.(F - F_)
end

function log_pseudolikelihood_sites(
    rbm::RBM{<:Spin}, v::AbstractArray, sites::AbstractVector{<:CartesianIndex}
)
    @assert size(visible(rbm)) == size(v)[1:ndims(visible(rbm))]
    @assert size(sites) == batch_size(visible(rbm), v)
    v_ = copy(v)
    for (b, i) in pairs(sites)
        v_[i, b] = -v_[i, b]
    end
    F = free_energy(rbm, v)
    F_ = free_energy(rbm, v_)
    return -log1pexp.(F - F_)
end


text:
using Test
using ManifoldsBase

struct NonManifold <: AbstractManifold{ManifoldsBase.ℝ} end

@testset "NotImplemented Errors" begin
    M = NonManifold()
    p = [1.0]
    q = similar(p)
    X = [2.0]
    Y = similar(X)
    for B in [
        VeeOrthogonalBasis(),
        DefaultBasis(),
        DefaultOrthogonalBasis(),
        DefaultOrthonormalBasis(),
        DiagonalizingOrthonormalBasis(X),
        CachedBasis(DefaultBasis(), X),
    ]
        if !(B isa CachedBasis)
            @test_throws MethodError get_basis(M, p, B)
            @test_throws MethodError get_vector(M, p, X, B)
            @test_throws MethodError get_vector!(M, Y, p, X, B)
        else
            @test get_basis(M, p, B) == B
            @test get_vector(M, p, X, B) == 4 # since we have 1 vector
            @test get_vector!(M, Y, p, X, B) == [4] # since Y is a vector
        end
        @test_throws MethodError get_coordinates(M, p, X, B)
        @test_throws MethodError get_coordinates!(M, Y, p, X, B)
    end
    @test_throws MethodError inverse_retract(M, p, q)
    @test_throws MethodError inverse_retract!(M, Y, p, q)
    @test_throws ErrorException PadeInverseRetraction(0)
    for IR in [
        CayleyInverseRetraction(),
        EmbeddedInverseRetraction(ProjectionInverseRetraction()),
        LogarithmicInverseRetraction(),
        PadeInverseRetraction(2),
        PolarInverseRetraction(),
        ProjectionInverseRetraction(),
        QRInverseRetraction(),
        NLSolveInverseRetraction(ExponentialRetraction()),
        SoftmaxInverseRetraction(),
    ]
        @test_throws MethodError inverse_retract(M, p, q, IR)
        @test_throws MethodError inverse_retract!(M, Y, p, q, IR)
    end
    @test_throws MethodError retract(M, p, X)
    @test_throws MethodError retract!(M, q, p, X)
    @test_throws ErrorException PadeRetraction(0)
    for R in [
        CayleyRetraction(),
        EmbeddedRetraction(ProjectionRetraction()),
        ExponentialRetraction(),
        ODEExponentialRetraction(ProjectionRetraction(), DefaultBasis()),
        ODEExponentialRetraction(ProjectionRetraction()),
        PadeRetraction(2),
        PolarRetraction(),
        ProjectionRetraction(),
        QRRetraction(),
        SoftmaxRetraction(),
        CayleyRetraction(),
        PadeRetraction(2),
    ]
        @test_throws MethodError retract(M, p, X, R)
        @test_throws MethodError retract!(M, q, p, X, R)
    end
    for VT in [
        DifferentiatedRetractionVectorTransport(ExponentialRetraction()),
        ProjectionTransport(),
        ParallelTransport(),
    ]
        @test_throws MethodError vector_transport_along(M, p, X, :curve, VT)
        @test_throws MethodError vector_transport_along!(M, Y, p, X, :curve, VT)
        @test_throws MethodError vector_transport_along!(M, Y, p, X, [p], VT)
        @test_throws MethodError vector_transport_direction(M, p, X, X, VT)
        @test_throws MethodError vector_transport_direction!(M, Y, p, X, X, VT)
        @test_throws MethodError vector_transport_to(M, p, X, q, VT)
        @test_throws MethodError vector_transport_to!(M, Y, p, X, q, VT)
    end
end

@testset "Default Fallbacks and Error Messages" begin
    M = ManifoldsBase.DefaultManifold(3)
    p = [1.0, 0.0, 0.0]
    @test number_of_coordinates(M, ManifoldsBase.ℝ) == 3
    B = get_basis(M, p, DefaultBasis())
    @test_throws DomainError ODEExponentialRetraction(ProjectionRetraction(), B)
    @test_throws DomainError ODEExponentialRetraction(
        ExponentialRetraction(),
        DefaultBasis(),
    )
    @test_throws DomainError ODEExponentialRetraction(ExponentialRetraction(), B)
    @test_throws ErrorException PadeRetraction(0)
end


text:
# Implementation empirical bond length distributions, supporting pdf() functions.

using Combinatorics
using CSV
using DataFrames
using LinearAlgebra

using SmuUtilities

export BondLengthDistribution, AllBondLengthDistributions
export from_file!, from_arrays!, pdf
export add_from_files!

"""A single bond length distribution, holding a relationship between
  distance and value.
  It is assumed that the X values are uniformly distributed.
"""
mutable struct BondLengthDistribution
  bucket_size::Float32  # Distance between X values.
  min_distance::Float32
  max_distance::Float32
  pdf::Vector{Float32}  # Normalized prevalence values.
end

BondLengthDistribution() = BondLengthDistribution(0.0, 0.0, 0.0, [])

function from_file!(fname::AbstractString, distribution::BondLengthDistribution)::Bool
  df = CSV.File(fname) |> DataFrame
  return from_arrays!(df[!,1], df[!,2], distribution)
end

function from_arrays!(lengths::Vector, counts::Vector, distribution::BondLengthDistribution)::Bool
  min_diff, max_diff = extrema(diff(lengths))
  if abs(max_diff - min_diff) > 1.0e-05
    @warn("non uniform buckets $(min_diff) $(max_diff)")
    return false
  end
  distribution.bucket_size = (max_diff + min_diff) / 2.0  # For accuracy.
  distribution.min_distance = lengths[1]
  distribution.max_distance =  lengths[end] + distribution.bucket_size
  # During testing/development, fill in any missing values with 1
  distribution.pdf = Float32.(map((x)-> x==0 ? 1 : x, counts))
  normalize!(distribution.pdf, 1)
  return true
end

function pdf(distribution::BondLengthDistribution, distance::Real)::Float32
  distance <= distribution.min_distance && return 0.0
  distance >= distribution.max_distance && return 0.0

  idx = round(Int32, (distance - distribution.min_distance) / distribution.bucket_size) + 1
  idx > length(distribution.pdf) && return 0.0  # will happen due to fp inaccuracy
  return distribution.pdf[idx]
end

"""The AllBondLengthDistributions struct keeps a map from this type
  to the corresponding BondLengthDistribution.
"""
struct AtypesBtype
  type1::Int
  type2::Int
  btype::Int
end

# functions to facilitate hashing
function Base.isequal(t1::AtypesBtype, t2::AtypesBtype)::Bool
  t1.type1 == t1.type1 && t1.type2 == t2.type2 && t1.btype == t2.btype
end

function Base.hash(atypes::AtypesBtype)::UInt64
  return 200 * atypes.type1 + 100 * atypes.type2 + atypes.btype
end

"""The main interface to bond length distributions.
  Holds a mapping from atom and bond types to a distribution, which si then
  used by pdf().
"""
mutable struct AllBondLengthDistributions
  # Whether or not the non-bonded types are read during construction
  # from files.
  include_non_bonded::Bool
  # For each atom types/bond type combination, a distribution.
  distribution::Dict{AtypesBtype, BondLengthDistribution}
end
AllBondLengthDistributions() = AllBondLengthDistributions(true, Dict{AtypesBtype, BondLengthDistribution}())

"""Add a bond length distribution to `distributions` by reading data from `fname`.
Args:
Returns:
"""
function add_file!(fname::String, distributions::AllBondLengthDistributions, key::AtypesBtype)::Bool
  dist = BondLengthDistribution()
  if ! from_file!(fname, dist)
    @warn("cannot read distribution from $(fname)")
    return false
  end

  @debug("add_file got $key")
  distributions.distribution[key] = dist

  return true
end

"""Given a file name `stem` read all files that would be present for bond length
distributions and construct `distributions`.

Note the very complex Iterators.product was inspired by
https://discourse.julialang.org/t/cleanest-way-to-generate-all-combinations-of-n-arrays/20127/8

No idea if this would be efficient for larger arrays, here it does not matter, and avoids a
doubly nested loop.

Args:
  stem: file name stem for the bond length distributions files. stem.6.1.6 would be for
    carbon-carbon single bonds
  distributions: The data read will be stored here.
Returns:
  true if successful
"""
function add_from_files!(stem::String, distributions::AllBondLengthDistributions)::Bool
  bstart  = distributions.include_non_bonded ? 0 : 1

  bond_types = [0, 1, 2, 3]

  atomic_numbers = [1, 6, 7, 8, 9]
  for (t1, t2) in filter(x->x[1] <= x[2], collect(Iterators.product(ntuple(i->atomic_numbers, 2)...))[:])
    for btype in bstart:3
      fname = "$stem.$(t1).$(btype).$(t2)"
      @info("Processing $fname")
      if ! isfile(fname)
        @warn("skipping non existent file $(fname)")
        continue
      end
      key = AtypesBtype(t1, t2, btype)
#     @debug("Adding to key $key")
      if ! add_file!(fname, distributions, key)
        @error("error adding $(fname)")
        return false
      end
    end
  end
  true
end

function pdf(distributions::AllBondLengthDistributions, type1::Integer, type2::Integer, btype::Integer, distance::Real)::Float32
  key = AtypesBtype(minmax(smu_atype_to_atomic_number(type1), smu_atype_to_atomic_number(type2))..., btype)
  value = get(distributions.distribution, key, nothing)
# @debug("value for $type1 $type2 $btype $value")
  value == nothing && return 0.0

  return pdf(value, distance)
end

max_stars_repo_name: IanAWatson/smu

max_issues_repo_name: IanAWatson/smu

max_forks_repo_name: IanAWatson/smu


text:
"""
ptcsolsc(f, x0, fp=difffp; rtol=1.e-6, atol=1.e-12, maxit=100,
        delta0=1.e-6, dx=1.e-7, pdata=nothing, printerr = true, keepsolhist=true)

C. T. Kelley, 2021

Scalar pseudo-transient continuation solver. PTC is designed to find
stable steady state solutions of

dx/dt = - f(x)

The scalar code is a simple wrapper around a call to ptcsol.jl, the
PTC solver for systems.

--> PTC is ABSOLUTELY NOT a general purpose nonlinear solver.

Input:\n
f: function\n
x: initial iterate/data\n
fp: derivative. If your derivative function is fp, you give me
its name. For example fp=foobar tells me that foobar is your
function for the derivative. The default is a forward difference
Jacobian that I provide.\n

Keyword Arguments:\n
rtol, atol: real and absolute error tolerances\n

maxit: upper bound on number of nonlinear iterations. This is
coupled to delta0. If your choice of delta0 is too small (conservative)
then you'll need many iterations to converge and will need a larger
value of maxit.

delta0: initial pseudo time step. The default value of 1.e-3 is a bit
conservative and is one option you really should play with. Look at the example
where I set it to 1.0!\n

dx: default = 1.e-7\n
difference increment in finite-difference derivatives
      h=dx*norm(x)+1.e-6

pdata:\n
precomputed data for the function/derivative.
Things will go better if you use this rather than hide the data
in global variables within the module for your function/derivative
If you use this option your function and derivative must take pdata
as a second argument. eg f(x,pdata) and fp(x,pdata)

printerr: default = true\n
I print a helpful message when the solver fails. To supress that
message set printerr to false.

keepsolhist: if true you get the history of the iteration in the output
tuple. This is on by default for scalar equations and off for systems.
Only turn it on if you have use for the data, which can get REALLY LARGE.

Output: A tuple (solution, functionval, history, idid, errcode, solhist) where
history is the array of absolute function values |f(x)|
of residual norms and time steps. Unless something has gone badly wrong,
delta approx |f(x_0)|/|f(x)|.

idid=true if the iteration succeeded and false if not.

errcode = 0 if if the iteration succeeded
        = -1 if the initial iterate satisifies the termination criteria
        = 10 if no convergence after maxit iterations

solhist=entire history of the iteration if keepsolhist=true\n
ptcsolsc builds solhist with a function from the Tools directory. For
systems, solhist is an N x K array where N is the length of x and K
is the number of iteration + 1. So, for scalar equations (N=1), solhist
is a row vector. Hence I use [ptcout.solhist' ptcout.history] in the
example below.

If the iteration fails it's time to play with the tolerances, delta0, and maxit.
You are certain to fail if there is no stable solution to the equation.

### Examples for ptcsolsc

```jldoctest
julia> ptcout=ptcsolsc(sptest,.2;delta0=2.0,rtol=1.e-3,atol=1.e-3);

julia> [ptcout.solhist' ptcout.history]
7×2 Array{Float64,2}:
 2.00000e-01  9.20000e-02
 9.66666e-01  4.19962e-01
 8.75086e-01  2.32577e-01
 7.99114e-01  1.10743e-01
 7.44225e-01  4.00926e-02
 7.15163e-01  8.19395e-03
 7.07568e-01  4.61523e-04
```

"""
function ptcsolsc(
    f,
    x0,
    fp = difffp;
    rtol = 1.e-6,
    atol = 1.e-12,
    maxit = 100,
    delta0 = 1.e-3,
    dx = 1.e-7,
    pdata = nothing,
    printerr = true,
    keepsolhist = true,
)
    #
    # The scalar code is a simple wrapper for the real code (ptcsol). The
    # wrapper puts placeholders for the memory allocations and the precomputed
    # data.
    #
    fp0 = copy(x0)
    fpp0 = copy(x0)
    itout = ptcsol(
        f,
        x0,
        fp0,
        fpp0,
        fp;
        rtol = rtol,
        atol = atol,
        maxit = maxit,
        delta0 = delta0,
        dx = dx,
        pdata = pdata,
        printerr = printerr,
        keepsolhist = keepsolhist,
    )
    #             printerr=printerr,keepsolhist=keepsolhist)
    return itout
end