text:
#= for significands in [1, 2), exponents are powers of 2
exponent_max(::Type{Float16})  =   15
exponent_max(::Type{Float32})  =  127
exponent_max(::Type{Float64})  = 1023
exponent_min(::Type{Float16})  = 1 -   15
exponent_min(::Type{Float32})  = 1 -  127
exponent_min(::Type{Float64})  = 1 - 1023
=#
# for significands in [0.5, 1)
exponent_max(::Type{Float16})  =   15+1
exponent_max(::Type{Float32})  =  127+1
exponent_max(::Type{Float64})  = 1023+1
exponent_min(::Type{Float16})  = 2 -   15
exponent_min(::Type{Float32})  = 2 -  127
exponent_min(::Type{Float64})  = 2 - 1023

for (U,F) in ((:UInt64, :Float64), (:UInt32, :Float32), (:UInt16, :Float16))
    @eval begin
        function randf(::Type{$F})
            r = $F(Inf)
            while !isfinite(r)
                r = reinterpret($F, rand($U))
            end
            return r
        end
        function randfloat(::Type{$F}, n::Int=1; emin::Int=exponent_min($F), emax::Int=exponent_max($F), signed::Bool=false)
            emin, emax = minmax(emin, emax)
            emin = max(exponent_min($F), emin)
            emax = min(exponent_max($F), emax)
            n = max(1, n)
            result = rand($F, n)
            # significands in [0.5, 1)
            result *= 0.5
            result += 0.5
            exponent = rand(emin:1:emax, n)
            for i in 1:n
                result[i] = ldexp(result[i], exponent[i])
            end
            if signed
               signs = rand(-1:2:1,n)
               for i in 1:n
                   result[i] = copysign(result[i], signs[i])
               end
            end
            return result
        end
    end
end


text:
function linear_index(ind)
    j = round(Int, floor(-0.5 + 0.5*sqrt(1 + 8*(ind - 1))) + 2)
    i = round(Int, j .* (3 - j)/2 + ind -1)
    return i, j
end

ext: jl


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

"""
Simple uninformed Monte Carlo search. Should be used as a baseline for other local solvers,
in terms of runtime performance and overall efficacy.
"""
module MonteCarloSearch

export
    MCS,
    solve

using ..Solvers
import ..Solvers.solve

using CommonRLInterface
using DataStructures: PriorityQueue, enqueue!
using ProgressMeter

"""
    MCS <: LocalSolver

Monte Carlo search solver.
"""
Base.@kwdef mutable struct MCS <: LocalSolver
    num_iterations::Int64 = 1000
    top_k::Int64          = 10
end

"""
    MCSResult <: LocalResult

Monte Carlo search result.
"""
mutable struct MCSResult <: LocalResult
    path::Vector
end

function Solvers.solve(mcs::MCS, env_fn::Function)
    mdp = env_fn()
    A = typeof(rand(actions(mdp)))
    best_paths = PriorityQueue()

    @showprogress for _ in 1:mcs.num_iterations
        reset!(mdp)
        path = A[]
        r = 0.0
        d = false
        while !d
            d = terminated(mdp)
            a = rand(actions(mdp))
            push!(path, a)
            r += act!(mdp, a)
            d && enqueue!(best_paths, MCSResult(path), r, mcs.top_k)
        end
    end

    return best_paths
end

end


text:
@userplot hornerplot

@recipe function f(j::hornerplot;
                    UnsteadyRange=false,Q=zeros(3), β=zeros(3), μ=ones(3), h=1,Rs=0.0, ct=1e-6, rw=0.583, ϕ=0.1, pr=true)

    T, P = j.args

    Tp=T[1]
    ΔT=T.-Tp
    HornerTime=(Tp.+ΔT)./ΔT
    H1hr=Tp.+1.0
    K=zeros(3)
        if UnsteadyRange!=false
       UnsteadyTime=HornerTime[(HornerTime.>=UnsteadyRange[1]) .& (HornerTime.<=UnsteadyRange[2])]
       UnsteadyPressure=P[(HornerTime.>=UnsteadyRange[1]) .& (HornerTime.<=UnsteadyRange[2])]

       m,b,r=Reg(UnsteadyTime,UnsteadyPressure,TypeReg="Logarithmic")
        UnsteadyLine=map(x->m*log10(x)+b,UnsteadyRange)
        UnsteadyExtrap=map(x->m*log10(x)+b,[UnsteadyRange[1],1])
       K[[1,2]]=map((q,be,mu)->(162.6*q*be*mu)/(abs(m)*h),Q[[1,2]],β[[1,2]],μ[[1,2]]) #liquis permeability Oil and Water
       K[3]= (162.6*(Q[3]*1000-Q[1]*Rs)*β[3]*μ[3])/(abs(m)*h)
       P1hr=m*log10(H1hr)+b
        λt=(K[1]/μ[1])+(K[2]/μ[2])+(K[3]/μ[3])
        S=1.151*(((P1hr-P[1])/abs(m))-(log10(λt/(ϕ*ct*rw^2)))+3.23)
        ΔPs=0.87*abs(m)*S
        Pi=b
        #Bg must be in bbl/scf
    end

    xlabel := "Horner Time"
    ylabel := "Pressure"
    xflip := true
    xscale := :log10
    legend := :topleft
    minorgrid := true
    gridcolor := :black
    ymirror := true

    @series  begin
        seriestype := :path
        linestyle := :dash
        linecolor := :black
        linewidth := 3
        label := "Horner Time"
        HornerTime, P
    end
    if  UnsteadyRange!=false
       @series begin
           seriestype := :path
           linestyle := :dash
           linecolor := :red
           linewidth := 1
           label := "Pi=$(round(Pi,digits=1)) psi"
           [UnsteadyRange[1],1],  UnsteadyExtrap
        end

        @series begin
            seriestype := :path
            linestyle := :solid
            linecolor := :red
            linewidth := 2
            label := "Unsteady Period"
            UnsteadyRange,  UnsteadyLine
         end
         pr==true ?  println("Unsteady State Parameters: \n Pi=$(round(Pi,digits=1)) psi  \n Ko=$(round(K[1],digits=1))md \n Kw=$(round(K[2],digits=1)) md \n Kg=$(round(K[3],digits=1)) md\n Skin=$(round(S,digits=2)) \n DeltaP Skin=$(round(ΔPs,digits=2)) psi ") : println(" ")
    end
end

@userplot derivativeplot

@recipe function f(j::derivativeplot)

        T, P = j.args

    Tp=T[1]
    Pwf=P[1]
    ΔT=T.-Tp
    ΔTe=ΔT./(1 .+(ΔT./Tp))
    ΔP=P.-Pwf
     δ=CentralDiff(ΔTe,ΔP)
    Der=((Tp.+ΔTe[2:end-1])./Tp).*δ.*ΔTe[2:end-1]

    xlabel := "Agarwal Time"
   ylabel := "Delta Pressure [psi]"
    legend := :topleft

    @series begin
        seriestype := :path
        linecolor := :black
        xscale := :log10
        yscale := :log10
        yformatter := :plain
        label := "Delta P"
        ΔTe[2:end], ΔP[2:end]
    end
      @series begin
        seriestype := :scatter
        markershape := :circle
        markersize := 2
        markercolor := :red
        xscale := :log10
        yscale := :log10
        yformatter := :plain
        label := "Derivative Function"
        ΔTe[2:end-1][Der.>0], Der[Der.>0]
    end

end


text:
"""
    retrieve(
        dset  :: AbstractString,
        dkeys :: AbstractDict,
        fnc   :: AbstractString,
        ckeys :: AbstractDict = cdskey()
    )
Retrieves dsets from the Climate Data Store, with options specified in a Julia Dictionary and saves it into a specified file.
Arguments:
    * `dset::AbstractString` : string specifies the name of the dset within the Climate Data Store that the `retrieve` function is attempting to retrieve data from
    * `dkeys::AbstractDict` : dictionary that contains the keywords that specify the properties (e.g. date, resolution, grid) of the data being retrieved
    * `fnc::AbstractString` : string that contains the path and name of the file that the data is to be saved into
    * `ckeys::AbstractDict` : dictionary that contains API Key information read from the .cdsapirc file in the home directory (optional)
"""
function retrieve(
    dset  :: AbstractString,
    dkeys :: AbstractDict,
    fnc   :: AbstractString,
    ckeys :: AbstractDict = cdskey()
)

    @info "$(now()) - CDSAPI - Welcome to the Climate Data Store"
    apikey = string("Basic ", base64encode(ckeys["key"]))

    @info "$(now()) - CDSAPI - Sending request to https://cds.climate.copernicus.eu/api/v2/resources/$(dset) ..."
    response = HTTP.request(
        "POST", ckeys["url"] * "/resources/$(dset)",
        ["Authorization" => apikey],
        body = JSON.json(dkeys),
        verbose = 0
    )
    resp_dict = JSON.parse(String(response.body))
    data = Dict("state" => "queued")

    @info "$(now()) - CDSAPI - Request is queued"
    while data["state"] == "queued"
        data = parserequest(ckeys,resp_dict,apikey)
    end

    @info "$(now()) - CDSAPI - Request is running"
    while data["state"] == "running"
        data = parserequest(ckeys,resp_dict,apikey)
    end

    if data["state"] == "completed"

        @info "$(now()) - CDSAPI - Request is completed"

        @info """$(now()) - CDSAPI - Downloading $(uppercase(dset)) data ...
          URL:         $(data["location"])
          Destination: $(fnc)
        """

        dt1 = now()
        HTTP.download(data["location"],fnc,update_period=Inf)
        dt2 = now()

        @info "$(now()) - CDSAPI - Downloaded $(@sprintf("%.1f",data["content_length"]/1e6)) MB in $(@sprintf("%.1f",Dates.value(dt2-dt1)/1000)) seconds (Rate: $(@sprintf("%.1f",data["content_length"]/1e3/Dates.value(dt2-dt1))) MB/s)"

    elseif data["state"] == "failed"

        @error "$(now()) - CDSAPI - Request failed"

    end

    return

end

"""
    cdskey() -> Dict{Any,Any}
Retrieves the CDS API ckeys from the `~/.cdsapirc` file in the home directory
"""
function cdskey()

    ckeys = Dict(); cdsapirc = joinpath(homedir(),".cdsapirc")

    @info "$(now()) - CDSAPI - Loading CDSAPI credentials from $(cdsapirc) ..."
    open(cdsapirc) do f
        for line in readlines(f)
            key,val = strip.(split(line,':',limit=2))
            ckeys[key] = val
        end
    end

    return ckeys

end

"""
    parserequest(
        ckeys :: AbstractDict,
        resp  :: AbstractDict,
        api   :: AbstractString
    ) -> Dict{Any,Any}
Get info on HTTP request, and parse the information and update the dictionary

Arguments
---------
* `ckeys` : Dictionary that contains the CDSAPI information held in `~/.cdsapirc`
* `resp` : Dictionary that contains the HTTP response
* `api`  : String that contains the API used for Climate Data Store authentication
"""
function parserequest(
    ckeys :: AbstractDict,
    resp  :: AbstractDict,
    api   :: AbstractString
)

    data = HTTP.request(
        "GET", ckeys["url"] * "/tasks/" * string(resp["request_id"]),
        ["Authorization" => api]
    )
    data = JSON.parse(String(data.body))

    return data

end


text:
"""
    get_column_id(table::Table, sym::Symbol)

Return the column id of the given symbol `sym`.
If the symbol doesn't exist throw an error.
"""
function get_column_id(table::Table, sym::Symbol)
    setup = table.setup
    ids = setup.ids
    for i in 1:length(ids)
        column_id = ids[i]
        if column_id == sym
            return i
            break
        end
    end
    error("Column symbol: $sym doesn't exist")
end

"""
    set_value!(table::Table, column::Symbol, value)

Set the current value of the given `column` to `value`
"""
function set_value!(table::Table, column::Symbol, value)
    id = get_column_id(table, column)
    table.current_values[id] = value
end

"""
    update_for_new_row(table)

Set the previous values to `current_values` by using `deepcopy` and set all `current_values` to nothing.
"""
function update_for_new_row(table)
    table.prev_values = deepcopy(table.current_values)
    table.current_values .= nothing
end

"""
    differs_enough(value, prev_value)

Return true if the `value` is significantly different from `prev_value`.
In this default case simply return `true` if they are different.
Provide your own `differs_enough` function by dispatching on both values.
"""
function differs_enough(value, prev_value)
    return value != prev_value
end

"""
    differs_enough(value::DiffX, prev_value::DiffX)

Check whether two [`DiffX`](@ref) values differ enough from each other.
"""
function differs_enough(value::DiffX, prev_value::DiffX)
    diff = abs(value.value - prev_value.value)
    return diff >= value.by
end

"""
    get_value(d::DiffX)

Taking the value of [`DiffX`](@ref) to represent it in the table.
"""
get_value(d::DiffX) = d.value
get_value(val) = val


text:
module PolyharmonicSplines

using LinearAlgebra

# Write your package code here.
include("polyharmonic.jl")
export PolyharmonicSpline

end


text:
@enum ClientCodes::UInt64 begin
    CLIENT_HELLO = 0
    CLIENT_QUERY = 1
    CLIENT_DATA = 2
    CLIENT_CANCEL = 3
    CLIENT_PING = 4
    CLIENT_TABLE_STATUS_REQ = 5
    CLIENT_KEEP_ALIVE = 6
end

@enum ServerCodes::UInt64 begin
    SERVER_HELLO = 0
    SERVER_DATA = 1
    SERVER_EXCEPTION = 2
    SERVER_PROGRESS = 3
    SERVER_PONG = 4
    SERVER_END_OF_STREAM = 5
    SERVER_PROFILE_INFO = 6
    SERVER_TOTALS = 7
    SERVER_EXTREMES = 8
    SERVER_TABLES_STATUS_REPORT = 9
    SERVER_TABLES_LOG = 10
    SERVER_TABLE_COLUMNS = 11
end

@reg_packet CLIENT_HELLO ClientHello
@reg_packet CLIENT_QUERY ClientQuery
@reg_packet CLIENT_DATA Block
@reg_packet CLIENT_PING ClientPing
#TODO CANCEL(3), TABLE_STATUS_REQ(5), KEEP_ALIVE(6)

@reg_packet SERVER_HELLO ServerInfo
@reg_packet SERVER_DATA ServerData
@reg_packet SERVER_EXCEPTION ServerException
@reg_packet SERVER_PROGRESS ServerProgress
@reg_packet SERVER_PONG ServerPong
@reg_packet SERVER_END_OF_STREAM ServerEndOfStream
@reg_packet SERVER_PROFILE_INFO ServerProfileInfo
@reg_packet SERVER_TOTALS ServerTotals
@reg_packet SERVER_EXTREMES ServerExtremes
@reg_packet SERVER_TABLE_COLUMNS ServerTableColumns

#TODO SERVER_TABLES_STATUS_REPORT(9), SERVER_TABLES_LOG(10), SERVER_TABLE_COLUMNS(11)

function read_packet(sock::ClickHouseSock, ::Type{CodeT}) where {CodeT}
    opcode = CodeT(UInt64(chread(sock, VarUInt)))
    struct_type = packet_struct(Val(opcode))
    return chread(sock, struct_type)
end

read_client_packet(sock::ClickHouseSock) = read_packet(sock, ClientCodes)

function read_server_packet(sock::ClickHouseSock)
    packet = read_packet(sock, ServerCodes)

    if typeof(packet) == ServerException
        throw(ClickHouseServerException(packet.code, packet.name, packet.message))
    end

    packet
end
function write_packet(sock::ClickHouseSock, packet::T) where {T}

    chwrite(sock, VarUInt(packet_code(T)))
    res = chwrite(sock, packet)
    flush(sock)
    return res
end


text:
include( Pkg.dir("GenomeAnnotations", "src","Config.jl") )

function setup_genome_annotations( config_ini_path = Pkg.dir("GenomeAnnotations", "conf","conf.ini"))

    local_store_path = LOCAL_STORE_PATH( config_ini_path )

    #mkdir if it doesn't exist
    make_localstore(local_store_path)
    println("Made local GenomeAnnotation store at $local_store_path")

    # create a db description file
    db_names_file = joinpath(local_store_path,"db.names")
    if isfile(db_names_file)
        println("$db_names_file already exists")
    else
        touch( db_names_file )
        println("created db description file at $db_names_file")
    end
end

function make_localstore(local_store_dir)
    if ! isdir(local_store_dir)
        mkpath(local_store_dir)
    else
       println("$local_store_dir exists, no need to create")
    end
end

function download_url ( url, dest)
    println("Running curl_cmd")
    if !isfile(dest)
        run(`curl -o $dest -L $url`)
        println("Downloaded $ensgene_url to $dest")
    else
        println("$dest exists no need to download")
    end
end

function current_annotation_dirs( )
   local_store_path = LOCAL_STORE_PATH()
   dirs = readdir( local_store_path )
   return dirs
end

function current_annotation_files()
   dirs = current_annotation_dirs()
   local_store_path = LOCAL_STORE_PATH()

   annotation_dict = Dict()
   for dir in dirs
       path = joinpath( local_store_path, dir )
       if isdir(path)
           paths = readdir( path )
           full_paths = map( x -> joinpath( path,x), paths )
           annotation_dict[dir] = full_paths
       end
   end
   return annotation_dict
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 9a0cec14-08db-11eb-3cfa-4d1c327c63f1
begin
	using Plots
	using PlutoUI
	using Statistics
end

# ╔═╡ a3b2accc-0845-11eb-229a-e97bc3943016
md"""
## From micro to macro & discrete to continuous
"""

# ╔═╡ ac4f6944-08d9-11eb-0b3c-f5e0a8b8c17e
md"""

Microscopic stochastic models, like the ones we have been studying in the homeworks, are intuitive and useful to include different effects. But often what we are actually interested in is a global picture of how an epidemic is evolving in time: how many infectious individuals there are in total at a given time, for example.

In this notebook we will see that it is often possible to summarise the **microscopic dynamics** of a **stochastic model** using **deterministic** equations for a **macroscopic** description.

Those macroscopic equations can either be in discrete time (**recurrence relations** or **difference equations**), or we can take a **continuous limit** and convert them into **ordinary differential equations**.

We're aiming for a plot like the following:
"""

# ╔═╡ f9a75ac4-08d9-11eb-3167-011eb698a32c
md"""
Let's go back to the model of recovery from an infection from the homework, with allowed transitions $I \to R$:

> we have $N$ people infected at time $0$. If each has probability $p$ to recover each day, how many are still infected at day number $t$?

Let's use a computational thinking approach: start off by coding the simulation and plotting the results:
"""

# ╔═╡ ba7ffe78-0845-11eb-2847-851a407dd2ec
bernoulli(p) = rand() < p

# ╔═╡ d088ed2e-0845-11eb-0697-310f374effbc
N = 200

# ╔═╡ e2d764d0-0845-11eb-0031-e74d2f5acaf9
function step!(infectious, p)
	for i in 1:length(infectious)

		if infectious[i] && bernoulli(p)
			infectious[i] = false
		end
	end

	return infectious
end

# ╔═╡ 9282eca0-08db-11eb-2e36-d761594b427c
T = 100

# ╔═╡ 58d8542c-08db-11eb-193a-398ce01b8635
begin
	infected = trues(N)

	results = [copy(step!(infected, 0.05)) for i in 1:T]
	pushfirst!(results, trues(N))
end

# ╔═╡ 8d6c0c06-08db-11eb-3790-c98fdc545352
@bind i Slider(1:T, show_value=true)

# ╔═╡ 7e1b61ac-08db-11eb-209e-1d6c328f5113
begin
	scatter(results[i],
		alpha=0.5, size=(300, 200), leg=false, c=Int.(results[i]))

	annotate!(N, 0.9, text("$(count(results[i]))", 10))
	annotate!(N, 0.1, text("$(N - count(results[i]))", 10))

	ylims!(-0.1, 1.1)

	xlabel!("i")
	ylabel!("X_i(t)")

end

# ╔═╡ 33f9fc36-0846-11eb-18c2-77f92fca3176
function simulate_recovery(p, T)
	infectious = trues(N)
	num_infectious = [N]

	for t in 1:T
		step!(infectious, p)
		push!(num_infectious, count(infectious))
	end

	return num_infectious
end

# ╔═╡ 39a69c2a-0846-11eb-35c1-53c68a9f71e5
p = 0.1

# ╔═╡ cb278624-08dd-11eb-3375-276bfe8d7b3a
begin
	pp = 0.05

	plot(simulate_recovery(pp, T), label="run 1", alpha=0.5, lw=2, m=:o)
	plot!(simulate_recovery(pp, T), label="run 2", alpha=0.5, lw=2, m=:o)

	xlabel!("time t")
	ylabel!("number infectious")
end

# ╔═╡ f3c85814-0846-11eb-1266-63f31f351a51
all_data = [simulate_recovery(pp, T) for i in 1:30];

# ╔═╡ 01dbe272-0847-11eb-1331-4360a575ff14
begin
	plot(all_data, alpha=0.1, leg=false, m=:o, ms=1,
		size=(500, 400), label="")
	xlabel!("time t")
	ylabel!("number infectious")
end

# ╔═╡ be8e4ac2-08dd-11eb-2f72-a9da5a750d32
plot!(mean(all_data), leg=true, label="mean",
		lw=3, c=:red, m=:o, alpha=0.5,
		size=(500, 400))

# ╔═╡ 8bc52d58-0848-11eb-3487-ef0d06061042
begin
	plot(replace.(all_data, 0.0 => NaN),
		yscale=:log10, alpha=0.3, leg=false, m=:o, ms=1,
		size=(500, 400))

	plot!(mean(all_data), yscale=:log10, lw=3, c=:red, m=:o, label="mean", alpha=0.5)

	xlabel!("time t")
	ylabel!("number infectious")
end

# ╔═╡ caa3faa2-08e5-11eb-33fe-cbbc00cfd459
md"""
## Deterministic dynamics for the mean: Intuitive derivation
"""

# ╔═╡ 2174aeba-08e6-11eb-09a9-2d6a882a2604
md"""
The mean seems to behave in a rather predictable way over time. Can we derive this?

Let $I_t$ be the number of infectious people at time $t$. This decreases because some people recover. Since people recover with probability $p$, the number of people that recover at time $t$ is, on average, $p I_t$. [Note that one time unit corresponds to one *sweep* of the simulation.]

So

$$I_{t+1} = I_t - p \, I_t$$

or

$$I_{t+1} = I_t (1 - p).$$

"""

# ╔═╡ 7e89f992-0847-11eb-3155-c5313575f767
md"""
At time $t$ there are $I_t$ infectious.
How many decay? Each decays with probability $p$, so on average $p I_t$ recover, so are removed from the number of infectious, giving the change

$$\Delta I_t = I_{t+1} - I_t = -p \, I_t$$
"""

# ╔═╡ f5756dd6-0847-11eb-0870-fd06ad10b6c7
md"""
We can rearrange and solve the recurrence relation:

$$I_{t+1} = (1 - p) I_t. $$

so

$$I_{t+1} = (1 - p) (1 - p) I_{t-1} = (1 - p)^2 I_{t-1}$$

and hence solve the recurrence relation:

$$I_t = (1-p)^t \, I_0.$$
"""

# ╔═╡ 113c31b2-08ed-11eb-35ef-6b4726128eff
md"""
Let's compare the exact and numerical results:
"""

# ╔═╡ 6a545268-0846-11eb-3861-c3d5f52c061b
exact = [N * (1-pp)^t for t in 0:T]

# ╔═╡ 4c8827b8-0847-11eb-0fd1-cfbdbdcf392e
begin
	plot(mean(all_data), m=:o, label="numerical mean",
		size=(500, 400))
	plot!(exact, lw=3, label="analytical mean")
end

# ╔═╡ 3cd1ad48-08ed-11eb-294c-f96b0e7c33bb
md"""
They agree well, as they should. The agreement is expected to be better (i.e. the fluctuations smaller) for a larger population.
"""

# ╔═╡ 57665108-08e7-11eb-2d45-311e07217c4e
md"""

## Derivation using mean of stochastic process

Instead of appealing to our intuition in the derivation, we can formulate an exact description of the stochastic process in terms of random variables: $X_t^i$ takes the value 1 if the $i$th person is infected at time $t$, and $0$ if they have recovered, as in the plot above.

Each coin flip is given by an independent Bernoulli random variable $B_t^i$, which has the value 1 (true) with probability $p$, and 0 (false) with probability $1 - p$.

The value at the next time step, $X_{t+1}^i$, is given by

$$X_{t+1}^i = \begin{cases} 0 &\text{if } X_t^i = 0 \text{ or } B_t^i = 1 \\
1 &\text{otherwise}
\end{cases}$$

We can write this as

$$X_{t+1}^i = X_t^i \, (1 - B_t^i).$$

(You should check that this is equivalent to the previous description!)

"""

# ╔═╡ 04d1ce96-0986-11eb-0c35-c73eb60d0994
md"""

$$\newcommand{mean}[1]{\left \langle #1 \right \rangle}$$

Let's denote the mean of a random variable $X$ by, denoted $\langle X \rangle$. In probability theory this is often called **expectation** (or "expected value") and written $\mathbb{E}[X]$.

Taking the expectation of both sides of the above equation, and using that the expectation of a product of independent random variables is the product of their expectations, we obtain

$$\mean{X_{t+1}^i} = \mean{X_t^i} \, \left(1 - \mean{B_t^i} \right).$$

The mean of each independent but identically distributed Bernoulli variable $B_t^i$ is $p$.

To find the mean total number of infected individuals we just sum up the boolean variables $X_t^i$:

$$I_t = \mean{\sum_{i=1}^N X_t^i}.$$
"""

# ╔═╡ 9d688084-0986-11eb-240f-23d85d9a3554
md"""
Putting all this together gives us back our previous result,

$$I_{t+1} = (1 - p) \, I_t .$$
"""

# ╔═╡ 2f980870-0848-11eb-3edb-0d4cd1ed5b3d
md"""
## Continuous time

If we look at the graph of the mean as a function of time, it seems to follow a smooth curve. Indeed it makes sense to ask not only how many people have recovered each *day*, but to aim for finer granularity.

Suppose we instead increment time in steps of $\delta t$; the above analysis was for $\delta t = 1$.

Then we will need to adjust the probability of recovery in each time step.
It turns out that to make sense in the limit $\delta t \to 0$, we need to choose the probability $p(\delta t)$ to recover in time $t$ to be proportional to $\delta t$:

$$p(\delta t) \simeq \lambda \, \delta t,$$

where $\lambda$ is the recovery **rate**. Note that a rate is a probability *per unit time*.

We get
"""

# ╔═╡ 6af30142-08b4-11eb-3759-4d2505faf5a0
md"""
$$I(t + \delta t) - I(t) \simeq -\lambda \,\delta t \, I(t)$$
"""

# ╔═╡ c6f9feb6-08f3-11eb-0930-83385ca5f032
md"""
Dividing by $\delta t$ gives

$$\frac{I(t + \delta t) - I(t)}{\delta t} \simeq -\lambda \, I(t)$$

We recognise the left-hand side as the definition of the **derivative** when $\delta t \to 0$. Taking that limit finally gives
"""

# ╔═╡ d8d8e7d8-08b4-11eb-086e-6fdb88511c6a
md"""
$$\frac{dI(t)}{dt} = -\lambda \, I(t)$$

That is, we obtain an **ordinary differential equation** that gives the solution implicitly. Solving this equation with initial condition $I(0) = I_0$ gives
"""

# ╔═╡ 780c483a-08f4-11eb-1205-0b8aaa4b1c2d
md"""
$$I(t) = I_0 \exp(-\lambda \, t).$$
"""

# ╔═╡ a13dd444-08f4-11eb-08f5-df9dd99c8ab5
md"""
Alternatively, we can derive this by recognising the exponential in the limit $\delta t \to 0$ of the following expression, which is basically the expression for compounding interest:
"""

# ╔═╡ cb99fe22-0848-11eb-1f61-5953be879f92
md"""
$$I_{t} = (1 - \lambda \, \delta t)^{(t / \delta t)} I_0$$
"""

# ╔═╡ d74bace6-08f4-11eb-2a6b-891e52952f57
md"""
## SIR model
"""

# ╔═╡ dbdf2812-08f4-11eb-25e7-811522b24627
md"""
Now let's extend the procedure to the full SIR model, $S \to I \to R$. Since we already know how to deal with recovery, consider just the SI model, where susceptible agents are infected via contact, with probability
"""

# ╔═╡ 238f0716-0903-11eb-1595-df71600f5de7
md"""
Let's denote by $S_t$ and $I_t$ be the number of susceptible and infectious people at time $t$, respectively, and by $N$ the total number of people.

On average, in each sweep each infectious individual has the chance to interact with one other individual. That individual is chosen uniformly at random from the total population of size $N$. But a new infection occurs only if that chosen individual is susceptible, which happens with probability $S_t / N$, and then if the infection is successful, with probability $b$, say.

Hence the change in the number of infectious people after that step is.

The decrease in $S_t$ is also given by $\Delta I_t$.
"""

# ╔═╡ 8e771c8a-0903-11eb-1e34-39de4f45412b
md"""
$$\Delta I_t = I_{t+1} - I_t = b \, I_t \, \left(\frac{S_t}{N} \right)$$
"""

# ╔═╡ e83fc5b8-0904-11eb-096b-8da3a1acba12
md"""
It is useful to normalize by $N$, so we define

$$s_t := \frac{S_t}{N}; \quad i_t := \frac{I_t}{N}; \quad r_t := \frac{R_t}{N}$$
"""

# ╔═╡ d1fbea7a-0904-11eb-377d-690d7a16aa7b
md"""
Including recovery with probability $c$ we obtain the **discrete-time SIR model**:
"""

# ╔═╡ dba896a4-0904-11eb-3c47-cbbf6c01e830
md"""
$$\begin{align}
s_{t+1} &= s_t - b \, s_t \, i_t \\
i_{t+1} &= i_t + b \, s_t \, i_t - c \, i_t\\
r_{t+1} &= r_t + c \, i_t
\end{align}$$
"""

# ╔═╡ 267cd19e-090d-11eb-0676-0f88b57da937
md"""
Again we can obtain this from the stochastic process by taking expectations (exercise!). [Hint: Ignore recovery to start with and take variables $Y_t^i$ that are $0$ if the person is susceptible and 1 if it is infected.]
"""

# ╔═╡ 4e3c7e62-090d-11eb-3d16-e921405a6b16
md"""
And again we can allow the processes to occur in steps of length $\delta t$ and take the limit $\delta t \to 0$. With rates $\beta$ and $\gamma$ we obtain the standard (continuous-time) **SIR model**:
"""

# ╔═╡ 72061c66-090d-11eb-14c0-df619958e2b6
md"""
$$\begin{align}
\textstyle \frac{ds(t)}{dt} &= -\beta \, s(t) \, i(t) \\
\textstyle \frac{di(t)}{dt} &= +\beta \, s(t) \, i(t) &- \gamma \, i(t)\\
\textstyle \frac{dr(t)}{dt} &= &+ \gamma \, i(t)
\end{align}$$
"""

# ╔═╡ c07367be-0987-11eb-0680-0bebd894e1be
md"""
We can think of this as a model of a chemical reaction with species S, I and R. The term $s(t) i(t)$ is known as the [**mass action**](https://en.wikipedia.org/wiki/Law_of_mass_action) form of interaction.

Note that no analytical solutions of these (simple) nonlinear ODEs are known as a function of time! (However, [parametric solutions are known](https://arxiv.org/abs/1403.2160).)
"""

# ╔═╡ f8a28ba0-0915-11eb-12d1-336f291e1d84
md"""
Below is a simulation of the discrete-time model. Note that the simplest numerical method to solve (approximately) the system of ODEs, the **Euler method**, basically reduces to solving the discrete-time model!  A whole suite of more advanced ODE solvers is provided in the [Julia `DiffEq` ecosystem](https://diffeq.sciml.ai/dev/).
"""

# ╔═╡ d994e972-090d-11eb-1b77-6d5ddb5daeab
begin
	NN = 100

	SS = NN - 1
	II = 1
	RR = 0
end

# ╔═╡ 050bffbc-0915-11eb-2925-ad11b3f67030
ss, ii, rr = SS/NN, II/NN, RR/NN

# ╔═╡ 1d0baf98-0915-11eb-2f1e-8176d14c06ad
p_infection, p_recovery = 0.1, 0.01

# ╔═╡ 28e1ec24-0915-11eb-228c-4daf9abe189b
TT = 1000

# ╔═╡ 349eb1b6-0915-11eb-36e3-1b9459c38a95
function discrete_SIR(s0, i0, r0, T=1000)

	s, i, r = s0, i0, r0

	results = [(s=s, i=i, r=r)]

	for t in 1:T

		Δi = p_infection * s * i
		Δr = p_recovery * i

		s_new = s - Δi
		i_new = i + Δi - Δr
		r_new = r      + Δr

		push!(results, (s=s_new, i=i_new, r=r_new))

		s, i, r = s_new, i_new, r_new
	end

	return results
end

# ╔═╡ 39c24ef0-0915-11eb-1a0e-c56f7dd01235
SIR = discrete_SIR(ss, ii, rr)

# ╔═╡ 442035a6-0915-11eb-21de-e11cf950f230
begin
	ts = 1:length(SIR)
	discrete_time_SIR_plot = plot(ts, [x.s for x in SIR],
		m=:o, label="S", alpha=0.2, linecolor=:blue, leg=:right, size=(400, 300))
	plot!(ts, [x.i for x in SIR], m=:o, label="I", alpha=0.2)
	plot!(ts, [x.r for x in SIR], m=:o, label="R", alpha=0.2)

	xlims!(0, 500)
end

# ╔═╡ 5f4516fe-098c-11eb-3abe-418aac994cc3
discrete_time_SIR_plot

# ╔═╡ Cell order:
# ╠═9a0cec14-08db-11eb-3cfa-4d1c327c63f1
# ╟─a3b2accc-0845-11eb-229a-e97bc3943016
# ╟─ac4f6944-08d9-11eb-0b3c-f5e0a8b8c17e
# ╠═5f4516fe-098c-11eb-3abe-418aac994cc3
# ╟─f9a75ac4-08d9-11eb-3167-011eb698a32c
# ╠═ba7ffe78-0845-11eb-2847-851a407dd2ec
# ╠═d088ed2e-0845-11eb-0697-310f374effbc
# ╠═e2d764d0-0845-11eb-0031-e74d2f5acaf9
# ╠═9282eca0-08db-11eb-2e36-d761594b427c
# ╠═58d8542c-08db-11eb-193a-398ce01b8635
# ╠═8d6c0c06-08db-11eb-3790-c98fdc545352
# ╟─7e1b61ac-08db-11eb-209e-1d6c328f5113
# ╠═33f9fc36-0846-11eb-18c2-77f92fca3176
# ╠═39a69c2a-0846-11eb-35c1-53c68a9f71e5
# ╠═cb278624-08dd-11eb-3375-276bfe8d7b3a
# ╠═f3c85814-0846-11eb-1266-63f31f351a51
# ╟─01dbe272-0847-11eb-1331-4360a575ff14
# ╟─be8e4ac2-08dd-11eb-2f72-a9da5a750d32
# ╟─8bc52d58-0848-11eb-3487-ef0d06061042
# ╟─caa3faa2-08e5-11eb-33fe-cbbc00cfd459
# ╟─2174aeba-08e6-11eb-09a9-2d6a882a2604
# ╟─7e89f992-0847-11eb-3155-c5313575f767
# ╟─f5756dd6-0847-11eb-0870-fd06ad10b6c7
# ╟─113c31b2-08ed-11eb-35ef-6b4726128eff
# ╟─6a545268-0846-11eb-3861-c3d5f52c061b
# ╟─4c8827b8-0847-11eb-0fd1-cfbdbdcf392e
# ╟─3cd1ad48-08ed-11eb-294c-f96b0e7c33bb
# ╟─57665108-08e7-11eb-2d45-311e07217c4e
# ╟─04d1ce96-0986-11eb-0c35-c73eb60d0994
# ╟─9d688084-0986-11eb-240f-23d85d9a3554
# ╟─2f980870-0848-11eb-3edb-0d4cd1ed5b3d
# ╟─6af30142-08b4-11eb-3759-4d2505faf5a0
# ╟─c6f9feb6-08f3-11eb-0930-83385ca5f032
# ╟─d8d8e7d8-08b4-11eb-086e-6fdb88511c6a
# ╟─780c483a-08f4-11eb-1205-0b8aaa4b1c2d
# ╟─a13dd444-08f4-11eb-08f5-df9dd99c8ab5
# ╟─cb99fe22-0848-11eb-1f61-5953be879f92
# ╟─d74bace6-08f4-11eb-2a6b-891e52952f57
# ╟─dbdf2812-08f4-11eb-25e7-811522b24627
# ╟─238f0716-0903-11eb-1595-df71600f5de7
# ╟─8e771c8a-0903-11eb-1e34-39de4f45412b
# ╟─e83fc5b8-0904-11eb-096b-8da3a1acba12
# ╟─d1fbea7a-0904-11eb-377d-690d7a16aa7b
# ╟─dba896a4-0904-11eb-3c47-cbbf6c01e830
# ╟─267cd19e-090d-11eb-0676-0f88b57da937
# ╟─4e3c7e62-090d-11eb-3d16-e921405a6b16
# ╟─72061c66-090d-11eb-14c0-df619958e2b6
# ╟─c07367be-0987-11eb-0680-0bebd894e1be
# ╟─f8a28ba0-0915-11eb-12d1-336f291e1d84
# ╠═442035a6-0915-11eb-21de-e11cf950f230
# ╠═d994e972-090d-11eb-1b77-6d5ddb5daeab
# ╠═050bffbc-0915-11eb-2925-ad11b3f67030
# ╠═1d0baf98-0915-11eb-2f1e-8176d14c06ad
# ╠═28e1ec24-0915-11eb-228c-4daf9abe189b
# ╠═349eb1b6-0915-11eb-36e3-1b9459c38a95
# ╠═39c24ef0-0915-11eb-1a0e-c56f7dd01235

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
struct StaticIRProjectState
    schema::Union{StaticAddressSchema,EmptyAddressSchema,AllAddressSchema}
    stmts::Vector{Any}
end

function process!(state::StaticIRProjectState, node) end

function process!(state::StaticIRProjectState, node::RandomChoiceNode)
    schema = state.schema
    @assert isa(schema, StaticAddressSchema) || isa(schema, EmptyAddressSchema) || isa(schema, AllAddressSchema)
    if isa(schema, AllAddressSchema) || (isa(schema, StaticAddressSchema) && (node.addr in keys(schema)))
        push!(state.stmts, :($weight += trace.$(get_score_fieldname(node))))
    end
end

function process!(state::StaticIRProjectState, node::GenerativeFunctionCallNode)
    schema = state.schema
    addr = QuoteNode(node.addr)
    @assert isa(schema, StaticAddressSchema) || isa(schema, EmptyAddressSchema) || isa(schema, AllAddressSchema)
    subtrace = get_subtrace_fieldname(node)
    subselection = gensym("subselection")
    if isa(schema, AllAddressSchema) || (isa(schema, StaticAddressSchema) && (node.addr in keys(schema)))
        push!(state.stmts, :($subselection = $qn_static_getindex(selection, Val($addr))))
        push!(state.stmts, :($weight += $qn_project(trace.$subtrace, $subselection)))
    else
        push!(state.stmts, :($weight += $qn_project(trace.$subtrace, $qn_empty_selection)))
    end
end

function codegen_project(trace_type::Type, selection_type::Type)
    gen_fn_type = get_gen_fn_type(trace_type)
    schema = get_address_schema(selection_type)

    # convert the selection to a static selection if it is not already one
    if !(isa(schema, StaticAddressSchema) || isa(schema, EmptyAddressSchema) || isa(schema, AllAddressSchema))
        return quote $qn_project(trace, $(QuoteNode(StaticSelection))(selection)) end
    end

    ir = get_ir(gen_fn_type)
    stmts = []

    # initialize weight
    push!(stmts, :($weight = 0.))

    # process expression nodes in topological order
    state = StaticIRProjectState(schema, stmts)
    for node in ir.nodes
        process!(state, node)
    end

    # return trace and weight
    push!(stmts, :(return $weight))

    Expr(:block, stmts...)
end

push!(generated_functions, quote

@generated function $(Expr(:(.), Gen, QuoteNode(:project)))(trace::T, selection::$(QuoteNode(Selection))) where {T <: $(QuoteNode(StaticIRTrace))}
    $(QuoteNode(codegen_project))(trace, selection)
end

function $(Expr(:(.), Gen, QuoteNode(:project)))(trace::T, selection::$(QuoteNode(EmptySelection))) where {T <: $(QuoteNode(StaticIRTrace))}
    trace.$total_noise_fieldname
end

end)

max_stars_repo_name: ali-ramadhan/Gen

max_issues_repo_name: ali-ramadhan/Gen

max_forks_repo_name: ali-ramadhan/Gen


text:
using Plots

# plot some data
plot([cumsum(rand(500) .- 0.5), cumsum(rand(500) .- 0.5)])

# save the current figure
savefig("plots.svg")
# .eps, .pdf, & .png are also supported
# we used svg here because it respects the width and height specified above

max_stars_repo_name: talgatomarov/JuliaByExample

max_issues_repo_name: talgatomarov/JuliaByExample

max_forks_repo_name: talgatomarov/JuliaByExample


text:
module LibPQ

export status, reset!, execute, prepare, async_execute, cancel,
    num_columns, num_rows, num_params, num_affected_rows

using Base: Semaphore, acquire, release
using Base.Iterators: zip, product
using Base.Threads

using Dates
using DocStringExtensions
using Decimals
using FileWatching
using Tables
using Infinity: InfExtendedTime, isposinf, ∞
using Intervals
using IterTools: imap
using LayerDicts
using Memento: Memento, getlogger, warn, info, error, debug
using OffsetArrays
using TimeZones

const Parameter = Union{String, Missing}
const LOGGER = getlogger(@__MODULE__)

function __init__()
    INTERVAL_REGEX[] = _interval_regex()
    Memento.register(LOGGER)
end

# Docstring template for types using DocStringExtensions
@template TYPES =
    """
        $(TYPEDEF)

    $(DOCSTRING)

    ## Fields:

    $(TYPEDFIELDS)
    """

include(joinpath(@__DIR__, "utils.jl"))

module libpq_c
    export Oid

    @static if VERSION < v"1.3.0"
        include(joinpath(@__DIR__, "..", "deps", "deps.jl"))

        function __init__()
            check_deps()
        end
    else
        using LibPQ_jll
    end

    include(joinpath(@__DIR__, "headers", "libpq-fe.jl"))
end

using .libpq_c

include("typemaps.jl")

const DEFAULT_CLIENT_TIME_ZONE = Ref("UTC")

"""
    const LIBPQ_TYPE_MAP::PQTypeMap

The [`PQTypeMap`](@ref) containing LibPQ-level type mappings for LibPQ.jl.
Adding type mappings to this constant will override the default type mappings for all code
using LibPQ.jl.
"""
const LIBPQ_TYPE_MAP = PQTypeMap()

"""
    const LIBPQ_CONVERSIONS::PQConversions

The [`PQConversions`](@ref) containing LibPQ-level conversion functions for LibPQ.jl.
Adding conversions to this constant will override the default conversions for all code using
LibPQ.jl.
"""
const LIBPQ_CONVERSIONS = PQConversions()

include("connections.jl")
include("results.jl")
include("statements.jl")
include("exceptions.jl")

include("parsing.jl")
include("copy.jl")
include("tables.jl")

include("asyncresults.jl")

include("deprecated.jl")

end


text:
@testset "Search space" begin
    @testset "in()" begin
        for i in 1:NumTestRepetitions
            reps = rand(1:10)
            ss1 = RectSearchSpace(reps, (0.0, 1.0))
            ind = rand_individual(ss1)
            for j in 1:reps
                @test (dimmin(ss1, j) <= ind[j] <= dimmax(ss1, j))
            end
        end
    end

    @testset "RectSearchSpace" begin
        @testset "RectSearchSpace() with given ranges per dimension" begin
            ss = RectSearchSpace([(0.0, 1.0)])
            @test dimmin(ss) == [0.0]
            @test dimmax(ss) == [1.0]
            @test dimdelta(ss) == [1.0]

            ss = RectSearchSpace([(0.0, 1.0), (0.5, 10.0)])
            @test dimmin(ss) == [0.0, 0.5]
            @test dimmax(ss) == [1.0, 10.0]
            @test dimdelta(ss) == [1.0, 9.5]
        end

        @testset "RectSearchSpace() search space with the default range" begin
            ss1 = RectSearchSpace(1)
            @test numdims(ss1) == 1
            @test dimmin(ss1) == [0.0]
            @test dimmin(ss1, 1) == 0.0
            @test dimmax(ss1) == [1.0]
            @test dimmax(ss1, 1) == 1.0
            @test dimdelta(ss1) == [1.0]
            @test dimdelta(ss1, 1) == 1.0
            @test dimrange(ss1) == [(0.0, 1.0)]
            @test dimrange(ss1, 1) == (0.0, 1.0)

            for i in 1:NumTestRepetitions
                ind = rand_individual(ss1)
                @test size(ind) == (1,)
                @test in(ind, ss1)
            end

            ss3 = RectSearchSpace(3)
            @test numdims(ss3) == 3
            @test dimmin(ss3) == [0.0, 0.0, 0.0]
            @test dimmin(ss3, 2) == 0.0
            @test dimmax(ss3) == [1.0, 1.0, 1.0]
            @test dimmax(ss3, 2) == 1.0
            @test dimdelta(ss3) == [1.0, 1.0, 1.0]
            @test dimdelta(ss3, 2) == 1.0
            @test dimrange(ss3) == [(0.0, 1.0), (0.0, 1.0), (0.0, 1.0)]
            @test dimrange(ss3, 1) == (0.0, 1.0)
            @test dimrange(ss3, 2) == (0.0, 1.0)
            @test dimrange(ss3, 3) == (0.0, 1.0)

            for i in 1:NumTestRepetitions
                ind = rand_individual(ss3)
                @test size(ind) == (3,)
                @test in(ind, ss3)
            end
        end

        @testset "ContinuousRectSearchSpace with given range" begin
            ss1 = RectSearchSpace(1, (-1.0, 1.0))
            @test_throws ArgumentError RectSearchSpace(1, (0.0, -1.0))
            @test ss1 isa ContinuousRectSearchSpace
            @test numdims(ss1) == 1
            @test dimrange(ss1) == [(-1.0, 1.0)]
            @test dimrange(ss1, 1) == (-1.0, 1.0)
            @test dimdigits(ss1) == [-1]
            @test dimdigits(ss1, 1) == -1
            #@test_throws BoundsError dimdigits(ss1, 2) # FIXME

            for i in 1:NumTestRepetitions
                dims = rand(1:100)
                a = rand()
                range = (a, a + (1-a)*rand())
                ss = RectSearchSpace(dims, range)
                @test numdims(ss) == dims
                @test all(dr -> dr == range, dimrange(ss))
            end
        end

        @testset "MixedPrecisionRectSearchSpace with given range" begin
            ss1 = RectSearchSpace(1, (-1.0, 1.0), dimdigits=2)
            @test ss1 isa MixedPrecisionRectSearchSpace
            @test numdims(ss1) == 1
            @test dimrange(ss1) == [(-1.0, 1.0)]
            @test dimmin(ss1) == fill(-1.0, numdims(ss1))
            @test dimmin(ss1, 1) == -1.0
            @test dimmax(ss1) == fill(1.0, numdims(ss1))
            @test dimmax(ss1, 1) == 1.0
            @test dimrange(ss1) == [(-1.0, 1.0)]
            @test dimdelta(ss1) == fill(2.0, numdims(ss1))
            @test dimdelta(ss1, 1) == 2.0
            @test dimdigits(ss1) == fill(2, numdims(ss1))
            @test dimdigits(ss1, 1) == 2
            @test_throws BoundsError dimdigits(ss1, 2)
        end
    end

    @testset "rand_individuals()" begin
        @testset "rand_individual($SS) is within the search space" for
                SS in (ContinuousRectSearchSpace, MixedPrecisionRectSearchSpace)
            for i in 1:NumTestRepetitions
                dims = rand(1:100)
                mm = sort!(rand(2))
                range = (mm[1], mm[2])
                ss = RectSearchSpace(dims, range, dimdigits = SS == MixedPrecisionRectSearchSpace ? rand(-1:2) : nothing)
                ind = rand_individual(ss)
                @test length(ind) == numdims(ss)
                @test in(ind, ss)
            end
        end

        @testset "rand_individuals() creates many individuals and all are within the search space" begin
            for i in 1:NumTestRepetitions
                reps = rand(1:10)
                mm = sort(rand(2,1), dims=1)
                range = (mm[1], mm[2])
                ss = RectSearchSpace(reps, range)
                numinds = rand(1:10)
                inds = rand_individuals(ss, numinds)
                @test size(inds,1) == numdims(ss)
                @test size(inds,2) == numinds
                for j in 1:numinds
                    @test in(inds[:,j], ss)
                end
            end
        end

        @testset "rand_individuals(method=:$method) correctly handles individual dimensions in $SS" for
                method in (:uniform, :latin_hypercube), SS in (ContinuousRectSearchSpace, MixedPrecisionRectSearchSpace)
            for _ in 1:NumTestRepetitions÷10
                numdimensions = rand(1:13)
                minbounds = rand(numdimensions)
                maxbounds = minbounds .+ rand(1:10, numdimensions) .* rand(numdimensions)
                digits = SS == MixedPrecisionRectSearchSpace ? rand(-1:2, numdimensions) : nothing
                if digits !== nothing
                    for d in 1:numdimensions
                        if digits[d] >= 0
                            minbounds[d] = round(minbounds[d], digits=digits[d])
                            maxbounds[d] = round(maxbounds[d], digits=digits[d])
                        end
                    end
                end
                ss = RectSearchSpace(tuple.(minbounds, maxbounds), dimdigits = digits)
                @test dimmin(ss) == minbounds
                @test dimmax(ss) == maxbounds
                @test round.(dimdelta(ss), digits=6) == round.(maxbounds .- minbounds, digits=6)

                # Now generate 100 individuals and make sure they are all within bounds
                inds = rand_individuals(ss, 100, method=method)
                @test size(inds, 2) == 100
                @inbounds for i in 1:size(inds, 2)
                    indi = view(inds, :, i)
                    @test indi == BlackBoxOptim.feasible(indi, ss)
                    @test all(minbounds .<= indi .<= maxbounds)
                end
            end
        end

        @testset "rand_individuals(method=:latin_hypercube) samples in LHS intervals" begin
            ss = RectSearchSpace([(0.0, 1.0), (2.0, 3.0), (4.0, 5.0)])

            inds = rand_individuals(ss, 2, method=:latin_hypercube)
            @test size(inds, 1) == 3
            @test size(inds, 2) == 2

            sorted = sort(inds, dims=2) # Sort per row --> in their ordered intervals
            @test (0.0 <= sorted[1,1] <= 0.5)
            @test (0.5 <= sorted[1,2] <= 1.0)

            @test (2.0 <= sorted[2,1] <= 2.5)
            @test (2.5 <= sorted[2,2] <= 3.0)

            @test (4.0 <= sorted[3,1] <= 4.5)
            @test (4.5 <= sorted[3,2] <= 5.0)
        end

        @testset "rand_individuals() unknown sampling method" begin
            @test_throws ArgumentError rand_individuals(RectSearchSpace([(0.0, 1.0), (2.0, 3.0)]), 5, method=:simple)
        end
    end

    @testset "feasible(x, ss) projects `x` to the search space `ss`" begin
        @testset "ContinuousSearchSpace" begin
            ss = RectSearchSpace([(0.0, 1.0), (2.0, 3.0), (4.0, 5.0)])

            @test BlackBoxOptim.feasible([1.1, 2.0, 4.0], ss) == [1.0, 2.0, 4.0]
            @test BlackBoxOptim.feasible([1.1, 3.0, 4.0], ss) == [1.0, 3.0, 4.0]
            @test BlackBoxOptim.feasible([1.1, 2.0, 5.0], ss) == [1.0, 2.0, 5.0]
            @test BlackBoxOptim.feasible([1.1, 3.0, 5.0], ss) == [1.0, 3.0, 5.0]

            @test BlackBoxOptim.feasible([-0.1, 2.0, 4.0], ss) == [0.0, 2.0, 4.0]
            @test BlackBoxOptim.feasible([-0.1, 3.0, 4.0], ss) == [0.0, 3.0, 4.0]
            @test BlackBoxOptim.feasible([-0.1, 2.0, 5.0], ss) == [0.0, 2.0, 5.0]
            @test BlackBoxOptim.feasible([-0.1, 3.0, 5.0], ss) == [0.0, 3.0, 5.0]

            @test BlackBoxOptim.feasible([0.0, 1.9, 4.0], ss) == [0.0, 2.0, 4.0]
            @test BlackBoxOptim.feasible([0.0, 1.9, 4.0], ss) == [0.0, 2.0, 4.0]
            @test BlackBoxOptim.feasible([1.0, 1.9, 5.0], ss) == [1.0, 2.0, 5.0]
            @test BlackBoxOptim.feasible([1.0, 1.9, 5.0], ss) == [1.0, 2.0, 5.0]

            @test BlackBoxOptim.feasible([0.0, 3.3, 4.0], ss) == [0.0, 3.0, 4.0]
            @test BlackBoxOptim.feasible([0.0, 3.2, 4.0], ss) == [0.0, 3.0, 4.0]
            @test BlackBoxOptim.feasible([1.0, 3.1, 5.0], ss) == [1.0, 3.0, 5.0]
            @test BlackBoxOptim.feasible([1.0, 3.9, 5.0], ss) == [1.0, 3.0, 5.0]

            @test BlackBoxOptim.feasible([-0.4, 3.3, 14.5], ss) == [0.0, 3.0, 5.0]
        end

        @testset "MixedPrecisionRectSearchSpace" begin
            ss = RectSearchSpace([(0.0, 1.0), (2.0, 3.0), (4.0, 5.0)], dimdigits=[-1, 2, 0])

            @test BlackBoxOptim.feasible([1.124, 2.001, 4.0], ss) == [1.0, 2.0, 4.0]
            @test BlackBoxOptim.feasible([0.124, 2.555, 3.9], ss) == [0.124, 2.56, 4.0]
            @test BlackBoxOptim.feasible([0.12345, 2.555, 3.5], ss) == [0.12345, 2.56, 4.0]
        end
    end

    @testset "dimrange()" begin
        ss = RectSearchSpace([(0.0, 1.0), (2.0, 3.0), (4.0, 5.0)])
        diams = dimdelta(ss)

        @test length(diams) == 3
        @test diams == [1.0, 1.0, 1.0]
    end

    @testset "vcat(ss1, ss2)" begin
        ss1 = RectSearchSpace([(0.0, 1.0), (2.0, 3.0), (4.0, 5.0)])
        ss2 = RectSearchSpace([(6.0, 7.0), (8.0, 9.0)])

        sscat = vcat(ss1, ss2)
        @test numdims(sscat) == 5
        @test dimmin(sscat) == [0.0, 2.0, 4.0, 6.0, 8.0]
        @test dimmax(sscat) == [1.0, 3.0, 5.0, 7.0, 9.0]
        @test dimdelta(sscat) == fill(1.0, 5)
        @test dimdigits(sscat) == fill(-1, 5)

        # concat continuous and mixed precision space
        ss3 = RectSearchSpace([(6.0, 7.0), (8.0, 9.0), (10.0, 11.0)], dimdigits=[1, 0, -1])
        sscat2 = vcat(ss1, ss3)
        @test sscat2 isa MixedPrecisionRectSearchSpace
        @test numdims(sscat2) == 6
        @test dimmin(sscat2) == [0.0, 2.0, 4.0, 6.0, 8.0, 10.0]
        @test dimmax(sscat2) == [1.0, 3.0, 5.0, 7.0, 9.0, 11.0]
        @test dimdelta(sscat2) == fill(1.0, 6)
        @test dimdigits(sscat2) == [-1, -1, -1, 1, 0, -1]
    end
end


text:
module Flo

using Revise

using GLMakie; GLMakie.activate!
using Makie.Colors
using LinearAlgebra: norm
using Statistics: quantile, mean, median
using Printf: @sprintf

include("Swarms.jl")
include("ExampleFunctions.jl")

using .Attractors
using .Swarms: Swarm, step!, RK4

include("FigureElements.jl")

const FPS = 144

"""
    findLimits(func, q)

Runs a short simulation to determine appropriate limits for displaying a particular function.
"""
function findLimits(func::Function; q = 0.001)::Tuple
    num_particles = 100

    attractor = Swarm(func, size = num_particles, positions = (rand(Float64, 3, num_particles) .-0.5) .* 2)
    history = Array{Float64}(undef, 3, num_particles, 10000)

    dims, particles, steps = map(x -> 1:size(history, x), [1,2,3])

    # Run a time-limited simulation
    for i in steps
        history[:,:,i] = attractor.positions
        step!(attractor, RK4)
    end

    high_quantile = x -> quantile(x, 1 - q)
    low_quantile = x -> quantile(x, q)

    upper_limits = mean([high_quantile(history[i,j,:]) for i in dims, j in particles], dims = 2)
    lower_limits = mean([low_quantile(history[i,j,:]) for i in dims, j in particles], dims = 2)

    limits = (Tuple([(zip(lower_limits, upper_limits)...)...]) .* 1.5)
end

"""
    dropdownMapping(func_names)

Generates a Dictionary that cleans a set of function names and creates a mapping
from styled names => function names.
"""
function dropdownMapping(func_names::Vector{Symbol})
    dropdown_style = str -> replace(str, "_" => " ") |> titlecase |> str -> replace(str, " " => "-")
    styled_names = (dropdown_style ∘ string).(func_names)
    return (Dict ∘ zip)(styled_names, func_names)
end

function julia_main()::Cint
    ode_func = Attractors.aizawa

    attractor = Observable(Swarm(ode_func, size = 2000, step_size = repeat([1e-2], 2000)))

    limits = Observable(findLimits(ode_func))

    include("src/Theme.jl")

    figure, ax, run_var = createFigure(limits = limits); display(figure)

    speed_slider, _, _, layout = labelslider!(figure.scene, "Speed", exp10.(range(-4, -1, length = 50)), startvalue = 1e-2, format = x -> @sprintf("%.1e", x))

    figure[3,2] = layout

    on(speed_slider.value) do speed
        attractor[].step_size .= speed
    end

    dropdown_dict = dropdownMapping(names(Flo.Attractors)[2:end])
    menu = createDropdown!(figure, [keys(dropdown_dict)...])

    on(menu.selection) do dropdown_item
        func = dropdown_dict[dropdown_item]
        attractor[] = eval(:(Swarm($(func))))
        ax.limits[] = findLimits(eval(:($(func))))
    end

    # createSliders!(figure, (σ = 10, ρ = 28, β = 8/3))
    solver = RK4

    points = Observable(Point3f.(Vector{Float64}.([eachcol(attractor[].positions)...])))

    scatter!(points, markersize = 0.025, markerspace = SceneSpace)

    while events(figure.scene).window_open.val
        # Main sim loop
        while run_var[] & events(figure.scene).window_open.val
            step!(attractor[], solver)
            points[] = [eachcol(attractor[].positions)...]
            notify(points)
            sleep(1/FPS)
        end
        sleep(1e-12)
    end

    return 0
end

end


text:
function dtsol(n, L::TronSparseMatrixCSC, r)
    zero = 0.0

    # Solve L'*x = r and store the result in r.

    r[n] = r[n]/L.diag_vals[n]
    @inbounds for j=n-1:-1:1
        temp = zero
        for k=L.colptr[j]:L.colptr[j+1]-1
            temp = temp + L.tril_vals[k]*r[L.rowval[k]]
        end
        r[j] = (r[j] - temp)/L.diag_vals[j]
    end

    return
end

function dtsol(n, L::TronDenseMatrix, r)
    zero = 0.0

    # Solve L'*x = r and store the result in r.

    r[n] = r[n]/L.vals[n,n]
    @inbounds for j=n-1:-1:1
        temp = zero
        for k=j+1:n
            temp = temp + L.vals[k,j]*r[k]
        end
        r[j] = (r[j] - temp)/L.vals[j,j]
    end

    return
end


text:
module LengthChannels
export LengthChannel, fullsizeof

"""
This package defines a type `LengthChannel{T} <: AbstractChannel{T}` which simply adds information about the length of the channel when it is iterated. The constructor behaves the same as the constructor for `Channel`, but takes an additional integer specifying the length. This length is not to be confused with the buffer size of the channel, referred to as `buf` in the example below. When a `LengthChannel` is iterated, it continues until it has iterated the specified number of elements, after that the channel is closed, even if there are more elements put in the channel.

Examples:

```julia
using LengthChannels, Test
len,buf = 10,4
lc = LengthChannel{Int}(len, buf) do ch
    for i = 1:100
        put!(ch, i)
    end
end

@test length(lc) == len
cc = collect(lc)
@test length(cc) == len
@test cc == 1:l
@test isopen(lc)
```

The constructor to a `LengthChannel` further takes a keyword argument `autoclose=false (default)` which determines if the channel closes automatically after having iterated for the specified length. It might be useful to keep it open if you want to iterate the specified length several times. Just make sure the channel is still being populated.

## Use as buffered iterator for machine learning
A typical use-case for a channel with length is as a buffered dataset for training machine learning models. Here, an example that reads audio data from disk and does some light pre-processing before putting it in the channel
```julia
using LengthChannels, Flux, Random
files      = readdir(path_to_datafiles)
buffersize = 10
dataset = LengthChannel{Vector{Float32}}(length(files), buffersize, spawn=true) do ch
    while true
        for file in shuffle(files)
            data = read_from_disk(file)
            put!(ch, data)
        end
    end
end
```
a batch iterator suitable for training CNN models in Flux can be obtained like so
```julia
batchsize  = 16
buffersize = 10
files      = readdir(path_to_datafiles)
dataset = LengthChannel{Array{Float32,4}}(length(files)÷batchsize, buffersize, spawn=true) do ch
    batch = Array{Float32,4}(undef,height,width,nchannels,batchsize)
    bi = 1
    while true
        for file in shuffle(files)
            data = read_from_disk(file)
            batch[:,:,:,bi] .= data
            bi += 1
            if bi > batchsize
                bi = 1
                put!(ch, copy(batch))
            end
        end
    end
end
```
where `height,width,nchannels` are integers specifying the size of your data.
You can now send the `dataset` to `Flux.train!(loss, ps, dataset, opt)`
"""
struct LengthChannel{T} <: AbstractChannel{T}
    ch::Channel{T}
    l::Int
    autoclose::Bool
end

"""
    LengthChannel{T}(l::Int, args...; autoclose=false, kwargs...)

Construct a channel with length `l`
"""
function LengthChannel{T}(l::Int, args...; autoclose=false, kwargs...) where T
    ch = Channel{T}(args...; kwargs...)
    LengthChannel{T}(ch, l, autoclose)
end

function LengthChannel(l::Int, args...; autoclose=false, kwargs...)
    ch = Channel{Any}(args...; kwargs...)
    LengthChannel{Any}(ch, l, autoclose)
end

"""
    LengthChannel{T}(f::Function, l::Int, args...; autoclose=false, kwargs...)

Construct a channel with length `l`, using function `f` to populate the channel, see [`Channel`](@ref) for help using this constructor.
"""
function LengthChannel{T}(f::Function, l::Int, args...; autoclose=false, kwargs...) where T
    ch = Channel{T}(f, args...; kwargs...)
    LengthChannel{T}(ch, l, autoclose)
end

function LengthChannel(f::Function, l::Int, args...; autoclose=false, kwargs...)
    ch = Channel{Any}(f, args...; kwargs...)
    LengthChannel{Any}(ch, l, autoclose)
end

"""
    LengthChannel(f::Function, channel::LengthChannel; kwargs...)
    LengthChannel{T}(f::Function, channel::LengthChannel; kwargs...)

Wraps a channel in a new channel, applying `f` to each element of the original channel. This can be useful to, e.g., allow the inner channel to read and pre-process data on a separate thread using `spawn=true`, while forcing the wrapper channel to operate on the main thread. This is required for `f` that are not thread safe, notably `f=cu` or `f=gpu` which puts data on a GPU.

If `T` is omitted while wrapping a channel, it is assumed that `f(eltype(dataset)) == typeof(f(::eltype(dataset)))` or in words, `f` must have a method returning the type resulting from applying `f` to an element of the wrapped channel.
"""
function LengthChannel(adaptor, channel::LengthChannel{T}; kwargs...) where T
    gT = adaptor(T)
    LengthChannel{gT}(adaptor, channel::LengthChannel; kwargs...)
end

function LengthChannel{gT}(adaptor, channel::LengthChannel; kwargs...) where gT
    LengthChannel{gT}(length(channel), channel.ch.sz_max; kwargs...) do ch
        while true
            for d in channel
                put!(ch, adaptor(d))
            end
        end
    end
end

Base.length(lc::LengthChannel) = lc.l

"""
    fullsizeof(lc::LengthChannel)

Return the sum of `sizeof` over all elements, assuming they all have the same `sizeof`.
"""
function fullsizeof(lc::LengthChannel)
    isready(lc) || error("Channel is not ready")
    sizeof(fetch(lc))*length(lc)
end

for f in (bind, close, fetch, isopen, isready, lock, popfirst!, push!, put!, take!, trylock, unlock, wait, eltype)
    f = nameof(f)
    @eval Base.$f(lc::LengthChannel, args...; kwargs...) = $(f)(lc.ch, args...; kwargs...)
end

function Base.iterate(lc::LengthChannel, state=0)
    if state == length(lc) || !isopen(lc)
        lc.autoclose && close(lc)
        return nothing
    end
     (take!(lc),state+1)
end

end


text:
export StraightLineHomotopy

"""
    StraightLineHomotopy(start, target)

Construct the homotopy `t * start + (1-t) * target`.

`start` and `target` have to match and to be one of the following
* `Vector{<:MP.AbstractPolynomial}` where `MP` is [`MultivariatePolynomials`](https://github.com/blegat/MultivariatePolynomials.jl)
* `MP.AbstractPolynomial`
* `Vector{<:FP.Polynomial}` where `FP` is [`FixedPolynomials`](https://github.com/saschatimme/FixedPolynomials.jl)

    StraightLineHomotopy{T}(start, target)

You can also force a specific coefficient type `T`.
"""
mutable struct StraightLineHomotopy{T<:Number} <: AbstractPolynomialHomotopy{T}
    start::Vector{FP.Polynomial{T}}
    target::Vector{FP.Polynomial{T}}

    function StraightLineHomotopy{T}(start::Vector{FP.Polynomial{T}}, target::Vector{FP.Polynomial{T}}) where {T<:Number}
        @assert length(start) == length(target) "Expected the same number of polynomials, but got $(length(start)) and $(length(target))"

        s_nvars = maximum(FP.nvariables.(start))
        @assert all(s_nvars .== FP.nvariables.(start)) "Not all polynomials of the start system have $(s_nvars) variables."

        t_nvars = maximum(FP.nvariables.(target))
        @assert all(t_nvars .== FP.nvariables.(target)) "Not all polynomials of the target system have $(t_nvars) variables."

        @assert s_nvars == t_nvars "Expected start and target system to have the same number of variables, but got $(s_nvars) and $(t_nvars)."
        new(start, target)
    end

    function StraightLineHomotopy{T}(start, target) where {T<:Number}
        s, t = construct(T, start, target)
        StraightLineHomotopy{T}(s, t)
    end
end

function StraightLineHomotopy(start, target)
    T, s, t = construct(start, target)
    StraightLineHomotopy{T}(s, t)
end

#
# SHOW
#
function Base.deepcopy(H::StraightLineHomotopy)
    StraightLineHomotopy(deepcopy(H.start), deepcopy(H.target))
end

#
# PROMOTION AND CONVERSION
#
Base.promote_rule(::Type{StraightLineHomotopy{T}}, ::Type{StraightLineHomotopy{S}}) where {S<:Number,T<:Number} = StraightLineHomotopy{promote_type(T,S)}
Base.promote_rule(::Type{StraightLineHomotopy}, ::Type{S}) where {S<:Number} = StraightLineHomotopy{S}
Base.promote_rule(::Type{StraightLineHomotopy{T}}, ::Type{S}) where {S<:Number,T<:Number} = StraightLineHomotopy{promote_type(T,S)}

Base.convert(::Type{StraightLineHomotopy{T}}, H::StraightLineHomotopy) where {T} = StraightLineHomotopy{T}(H.start, H.target)

#
# EVALUATION + DIFFERENTATION
#
function evaluate!(u::AbstractVector, H::StraightLineHomotopy{T}, x::Vector, t::Number) where T
    for i = 1:length(H.target)
        f = H.target[i]
        g = H.start[i]
        u[i] = (one(T) - t) * FP.evaluate(f, x) + t * FP.evaluate(g, x)
    end
    u
end
(H::StraightLineHomotopy)(x,t) = evaluate(H,x,t)

function evaluate!(u::AbstractVector{T}, H::StraightLineHomotopy, x::Vector, t::Number, cfg::PolynomialHomotopyConfig, precomputed=false) where {T<:Number}
    evaluate_start_target!(cfg, H, x, precomputed)
    u .= (one(t) - t) .* value_target(cfg) .+ t .* value_start(cfg)
end

function jacobian!(u::AbstractMatrix, H::StraightLineHomotopy{T}, x::AbstractVector, t, cfg::PolynomialHomotopyConfig, precomputed=false) where {T<:Number}
    jacobian_start_target!(cfg, H, x, precomputed)

    u .= (one(t) - t) .* jacobian_target(cfg) .+ t .* jacobian_start(cfg)
end

function jacobian!(r::JacobianDiffResult, H::StraightLineHomotopy{T}, x::AbstractVector, t, cfg::PolynomialHomotopyConfig, precomputed=false) where {T<:Number}
    evaluate_and_jacobian_start_target!(cfg, H, x)

    r.value .= (one(t) - t) .* value_target(cfg) .+ t .* value_start(cfg)
    r.jacobian .= (one(t) - t) .* jacobian_target(cfg) .+ t .* jacobian_start(cfg)
    r
end

function dt!(u, H::StraightLineHomotopy{T}, x::AbstractVector, t, cfg::PolynomialHomotopyConfig, precomputed=false) where {T<:Number}
    evaluate_start_target!(cfg, H, x, precomputed)
    u .= value_start(cfg) .- value_target(cfg)
end

function dt!(r::DtDiffResult, H::StraightLineHomotopy{T}, x::AbstractVector, t, cfg::PolynomialHomotopyConfig, precomputed=false) where {T<:Number}
    evaluate_start_target!(cfg, H, x, precomputed)
    r.value .= (one(T) - t) .* value_target(cfg) .+ t .* value_start(cfg)
    r.dt .= value_start(cfg) .- value_target(cfg)
    r
end

function weylnorm(H::StraightLineHomotopy{T})  where {T<:Number}
    f = FP.homogenize.(H.start)
    g = FP.homogenize.(H.target)
    λ_1 = FP.weyldot(f,f)
    λ_2 = FP.weyldot(f,g)
    λ_3 = FP.weyldot(g,g)

    function (t)
        sqrt(abs2(one(T) - t) * λ_1 + 2 * real((one(T) - t) * conj(t) * λ_2) + abs2(t) * λ_3)
    end
end


text:
# This file is a part of BoostFractor.jl, licensed under the MIT License (MIT).

__precompile__(true)

module BoostFractor

using FFTW

# General Interface
include("boost_fractor.jl")

# Differnt algorithms
include("boost_fractor_recurprop.jl")
include("boost_fractor_transformer.jl")

# Some convenient tools
include("beam_gaussian.jl")
include("beam_coupling.jl")

end # module


text:
function unsafe_broadcasttimes!(a::Vector{T}, b::Vector{T}) where T
    @simd for i in 1:length(a)
        Base.unsafe_setindex!(a, Base.unsafe_getindex(a, i)*Base.unsafe_getindex(b, i), i)
    end
    a
end

for op in (:(+),:(-))
    @eval begin
        function $op(B::AbstractBarycentricMatrix{T}) where T
            B2 = deepcopy(B)
            B2.F[:] = $op(B.F)
            B2
        end
        function $op(B1::EvenBarycentricMatrix{T}, B2::EvenBarycentricMatrix{T}) where T
            @assert B1.a == B2.a && B1.b == B2.b && B1.c == B2.c && B1.d == B2.d
            @assert length(B1.x) == length(B2.x)
            B3 = deepcopy(B1)
            B3.F[:] = $op(B1.F, B2.F)
            B3
        end
    end
end

# This file implements mul! overrides when A is a small matrix and u and v
# are larger vectors, but we want to operate in-place, equivalent to:
#
# u[istart:istart+size(A, 1)-1] += A*v[jstart:jstart+size(A, 2)-1]
#

# Generic Matrix

mul!(y::AbstractVecOrMat{T}, A::AbstractMatrix{T}, x::AbstractVecOrMat{T}, istart::Int, jstart::Int)  where T =
    mul!(y, A, x, istart, jstart, 1, 1)

function mul!(y::AbstractVecOrMat{T}, A::AbstractMatrix{T}, x::AbstractVecOrMat{T}, istart::Int, jstart::Int, INCX::Int, INCY::Int) where T
    m, n = size(A)
    ishift, jshift = istart-INCY, jstart-INCX
    @inbounds for j = 1:n
        xj = x[jshift+j*INCX]
        for i = 1:m
            y[ishift+i*INCY] += A[i,j]*xj
        end
    end

    y
end

mul!(y::AbstractVecOrMat{T}, At::Transpose{T,<:AbstractMatrix{T}}, x::AbstractVecOrMat{T}, istart::Int, jstart::Int) where T = mul!(y, At, x, istart, jstart, 1, 1)

function mul!(y::AbstractVecOrMat{T}, At::Transpose{T,<:AbstractMatrix{T}}, x::AbstractVecOrMat{T}, istart::Int, jstart::Int, INCX::Int, INCY::Int) where T
    A = parent(At)
    m, n = size(A)
    ishift, jshift = istart-INCY, jstart-INCX
    @inbounds for i = 1:n
        yi = zero(eltype(y))
        for j = 1:m
            yi += A[j,i]*x[jshift+j*INCX]
        end
        y[ishift+i*INCY] += yi
    end

    y
end

mul!(y::AbstractVecOrMat{T}, Ac::Adjoint{T,<:AbstractMatrix{T}}, x::AbstractVecOrMat{T}, istart::Int, jstart::Int) where T = mul!(y, Ac, x, istart, jstart, 1, 1)

function mul!(y::AbstractVecOrMat{T}, Ac::Adjoint{T,<:AbstractMatrix{T}}, x::AbstractVecOrMat{T}, istart::Int, jstart::Int, INCX::Int, INCY::Int) where T
    A = parent(Ac)
    m, n = size(A)
    ishift, jshift = istart-INCY, jstart-INCX
    @inbounds for i = 1:n
        yi = zero(eltype(y))
        for j = 1:m
            yi += conj(A[j,i])*x[jshift+j*INCX]
        end
        y[ishift+i*INCY] += yi
    end

    y
end

# LowRankMatrix

mul!(y::AbstractVecOrMat{T}, L::LowRankMatrix{T}, x::AbstractVecOrMat{T}) where T = mul!(y, L, x, 1, 1)
mul!(y::AbstractVecOrMat{T}, L::LowRankMatrix{T}, x::AbstractVecOrMat{T}, istart::Int, jstart::Int) where T = mul!(y, L, x, istart, jstart, 1, 1)

function mul!(y::AbstractVecOrMat{T}, L::LowRankMatrix{T}, x::AbstractVecOrMat{T}, istart::Int, jstart::Int, INCX::Int, INCY::Int) where T
    m, n = size(L)
    ishift, jshift = istart-INCY, jstart-INCX
    temp = L.temp

    @inbounds for k = 1:rank(L)
        temp[k] = zero(T)
        for j = 1:n
            temp[k] += L.V[j,k]*x[jshift+j*INCX]
        end
        temp[k] *= L.Σ[k,k]
    end

    @inbounds for k = 1:rank(L)
        tempk = temp[k]
        for i = 1:m
            y[ishift+i*INCY] += L.U[i,k]*tempk
        end
    end

    y
end

# BLAS'ed
# include("blas.jl")

# BarycentricMatrix

mul!(u::Vector{T}, A::AbstractBarycentricMatrix{T}, v::AbstractVector{T}) where T = mul!(u, A, v, 1, 1)

function mul!(u::Vector{T}, B::EvenBarycentricMatrix{T}, v::AbstractVector{T}, istart::Int, jstart::Int) where T
    β, W, F = B.β, B.W, B.F
    ishift, jshift, n = istart-1, jstart-1, length(β)

    if iseven(ishift+jshift)
        @inbounds for k = 1:n
            βk = zero(T)
            for j = 1:2:size(B, 2)
                βk += v[jshift+j]*F[j,k]
            end
            β[k] = βk
        end

        @inbounds for i = 1:2:size(B, 1)
            ui = zero(T)
            for k = 1:n
                ui += β[k]*W[k,i]
            end
            u[ishift+i] += ui
        end

        @inbounds for k = 1:n
            βk = zero(T)
            for j = 2:2:size(B, 2)
                βk += v[jshift+j]*F[j,k]
            end
            β[k] = βk
        end

        @inbounds for i = 2:2:size(B, 1)
            ui = zero(T)
            for k = 1:n
                ui += β[k]*W[k,i]
            end
            u[ishift+i] += ui
        end
    else
        @inbounds for k = 1:n
            βk = zero(T)
            for j = 2:2:size(B, 2)
                βk += v[jshift+j]*F[j,k]
            end
            β[k] = βk
        end

        @inbounds for i = 1:2:size(B, 1)
            ui = zero(T)
            for k = 1:n
                ui += β[k]*W[k,i]
            end
            u[ishift+i] += ui
        end

        @inbounds for k = 1:n
            βk = zero(T)
            for j = 1:2:size(B, 2)
                βk += v[jshift+j]*F[j,k]
            end
            β[k] = βk
        end

        @inbounds for i = 2:2:size(B, 1)
            ui = zero(T)
            for k = 1:n
                ui += β[k]*W[k,i]
            end
            u[ishift+i] += ui
        end
    end

    u
end

# BarycentricMatrix2D

function mul!(u::Vector{T}, B::BarycentricMatrix2D{T}, v::AbstractVector{T}, istart::Int, jstart::Int) where T
    U, F, V, temp1, temp2 = B.U, B.B.F, B.V, B.temp1, B.temp2
    ishift, jshift, r = istart-1, jstart-1, length(temp1)

    # temp1 = V'*v[jshift+j]

    for k in 1:r
        temp1k = zero(T)
        for j in 1:size(V, 1)
            temp1k += V[j,k]*v[jshift+j]
        end
        temp1[k] = temp1k
        temp2[k] = zero(T)
    end

    # temp2 = F*temp1

    for l in 1:r
        temp1l = temp1[l]
        for k in 1:r
            temp2[k] += F[k,l]*temp1l
        end
    end

    # u[ishift+i] = U*temp2

    for k in 1:r
        temp2k = temp2[k]
        for i in 1:size(U, 1)
            u[ishift+i] += U[i,k]*temp2k
        end
    end

    u
end

# C = A*Diagonal(b[jstart:jstart+size(A, 2)-1])

function scale!(C::Matrix, A::Matrix, b::Vector, jstart::Int)
    m, n = size(A)
    p, q = size(C)
    jshift = jstart-1
    @inbounds for j = 1:n
        bj = b[jshift+j]
        for i = 1:m
            C[i,j] = A[i,j]*bj
        end
    end
    C
end

function scale!(C::LowRankMatrix, A::LowRankMatrix, b::Vector, jstart::Int)
    scale!(C.V, b, A.V, jstart)
    C
end

# C = Diagonal(b[istart:istart+size(A, 1)-1])*A

function scale!(C::Matrix, b::Vector, A::Matrix, istart::Int)
    m, n = size(A)
    p, q = size(C)
    ishift = istart-1
    @inbounds for j = 1:n, i = 1:m
        C[i,j] = A[i,j]*b[ishift+i]
    end
    C
end

function scale!(C::LowRankMatrix, b::Vector, A::LowRankMatrix, istart::Int)
    scale!(C.U, b, A.U, istart)
    C
end

# Add e_j^T u[istart:istart+size(L, 1)-1] to A.

function add_col!(A::AbstractMatrix{T}, u::Vector{T}, istart::Int, j::Int) where T
    ishift, m = istart-1, size(A, 1)
    jm = (j-1)*m
    @simd for i = 1:m
        @inbounds A[i+jm] += u[i+ishift]
    end
    A
end

function add_col!(A::Matrix{T}, u::Vector{T}, istart::Int, j::Int) where T<:BlasReal
    m, n = size(A, 1), size(A, 2)
    BLAS.axpy!(m, one(T), pointer(u, istart), 1, pointer(A, (j-1)*m+1), 1)
    A
end

function add_col(A::AbstractMatrix{T}, u::Vector{T}, istart::Int, j::Int) where T
    B = deepcopy(A)
    ishift, m = istart-1, size(A, 1)
    jm = (j-1)*m
    @simd for i = 1:m
        @inbounds B[i+jm] += u[i+ishift]
    end
    B
end

function add_col(A::Matrix{T}, u::Vector{T}, istart::Int, j::Int) where T<:BlasReal
    m, n = size(A, 1), size(A, 2)
    B = Matrix{T}(undef, m, n)
    BLAS.blascopy!(m*n, A, 1, B, 1)
    BLAS.axpy!(m, one(T), pointer(u, istart), 1, pointer(B, (j-1)*m+1), 1)
    B
end

# Add e_j^T u[istart:istart+size(L, 1)-1] to UΣV^T.

function add_col(L::LowRankMatrix{T}, u::Vector{T}, istart::Int, j::Int) where T
    U, Σ, V = L.U, L.Σ, L.V
    U1 = hcat(U, u[istart:istart+size(L,1)-1])
    Σ1 = Diagonal(vcat(Σ.diag, one(T)))
    V1 = hcat(V, zeros(T, size(L, 2)))
    V1[j,rank(L)+1] = one(T)

    LowRankMatrix(U1,Σ1,V1)
end

function add_col(L::LowRankMatrix{T}, u::Vector{T}, istart::Int, j::Int) where T<:BlasReal
    U, Σ, V = L.U, L.Σ, L.V
    m, n, r, un = size(L, 1), size(L, 2), rank(L), one(T)
    U1 = Matrix{T}(undef, m, r+1)
    BLAS.blascopy!(m, pointer(u, istart), 1, U1, 1)
    BLAS.blascopy!(m*r, U, 1, pointer(U1, m+1), 1)
    diag = Vector{T}(r+1)
    BLAS.blascopy!(r, Σ.diag, 1, pointer(diag, 2), 1)
    Base.unsafe_setindex!(diag, un, 1)
    Σ1 = Diagonal(diag)
    V1 = zeros(T, n, r+1)
    BLAS.blascopy!(n*r, V, 1, pointer(V1, n+1), 1)
    Base.unsafe_setindex!(V1, un, j, 1)

    LowRankMatrix(U1,Σ1,V1)
end

function update!(::Type{T}, f::Function, A::AbstractMatrix{T}, x::Vector{T}, y::Vector{T}, ir::UnitRange{Int}, jr::UnitRange{Int}) where T
    ishift = 1-first(ir)
    jshift = 1-first(jr)

    for j in jr, i in ir
        A[i+ishift,j+jshift] = f(T, x[i], y[j])
    end

    A
end


text:
"""
```
function thin_mh_draws(m::AbstractDSGEModel, params::Matrix{Float64}; jstep::Int64 = 1)
```
Allows for a custom thinning factor (jstep) to be specified.
If not, it pulls the jstep from the model
"""
function thin_mh_draws(m::AbstractDSGEModel, params::Matrix{Float64}; jstep::Int64 = 1)
    jstep = jstep == 1 ? m.settings[:forecast_jstep].value : jstep
    n_total_draws, n_params = size(params)

    # Thin as usual if n_total_draws % jstep == 0
    # If it does not evenly divide, then start from the remainder+1-th index
    # and then take a thinned subset from there
    n_new_draws, offset = divrem(n_total_draws, jstep)
    params_thinned = Matrix{Float64}(undef, n_new_draws, n_params)
    params_offset = params[offset+1:end, :]

    for (i, j) in enumerate(1:jstep:n_total_draws)
        params_thinned[i, :] = params_offset[j, :]
    end
    return params_thinned
end

"""
```
function calculate_mode(m::Union{AbstractDSGEModel, AbstractVARModel}, data::Matrix{S},
                        starting_mode::Vector{S}, method::Symbol = :csminwel; mle::Bool = false,
                        verbose::Symbol = :none) where {S <: Real}
```
uses a numerical optimizer to calculate the mode of the posterior distribution
starting from pre-specified parameters. This function automatically saves
the mode found to an HDF5 file.
"""
function calculate_mode(m::Union{AbstractDSGEModel, AbstractVARModel}, data::Matrix{S},
                        starting_params::Vector{S}, method::Symbol = :csminwel;
                        mle::Bool = false, get_all_results::Bool = false,
                        save_results::Bool = true, check_neg_diag::Bool = true,
                        verbose::Symbol = :none) where {S <: Real}
    n_iterations = get_setting(m, :optimization_iterations)
    ftol         = get_setting(m, :optimization_ftol)
    xtol         = get_setting(m, :optimization_xtol)
    gtol         = get_setting(m, :optimization_gtol)
    step_size    = get_setting(m, :optimization_step_size)
    max_attempts = get_setting(m, :optimization_attempts)

    DSGE.update!(m, starting_params)

    modal_params  = optimize_step!(m, data, method; mle = mle, n_iterations = n_iterations,
                                   ftol = ftol, xtol = xtol, gtol = gtol, step_size = step_size,
                                   max_attempts = max_attempts, get_all_results = get_all_results,
                                   save_results = save_results, verbose = verbose)
    modal_hessian = hessian_step!(m, data; check_neg_diag = check_neg_diag, verbose = verbose)

    if get_all_results
        return modal_params..., modal_hessian
    else
        return modal_params, modal_hessian
    end
end

"""
```
function optimize_step!(m::Union{AbstractDSGEModel, AbstractVARModel}, data::Matrix{S},
                       method::Symbol = :csminwel; mle::Bool = false, n_iterations::Int = 4,
                       ftol::S = 1e-10, xtol::S = 1e-10, gtol::S = 1e-10, step_size::S = 0.01,
                       max_attempts::Int = 4, verbose::Symbol = :none) where {S <: Real}
```
maximizes the posterior or likelihood of the model object passed using
a numerical optimizer, which by default is csminwel. This function
automatically saves the mode found.
"""
function optimize_step!(m::Union{AbstractDSGEModel, AbstractVARModel}, data::Matrix{S},
                        method::Symbol = :csminwel; mle::Bool = false, n_iterations::Int = 4,
                        ftol::S = 1e-10, xtol::S = 1e-10, gtol::S = 1e-10, step_size::S = 0.01,
                        max_attempts::Int = 4, get_all_results::Bool = false,
                        save_results::Bool = true,
                        verbose::Symbol = :none) where {S <: Real}
    converged = false

    # If the algorithm stops only because we have exceeded the maximum number of
    # iterations, continue improving guess of modal parameters
    total_iterations = 0
    optimization_time = 0
    attempts = 1

    params = map(θ -> θ.value, get_parameters(m)) # define in global scope first
    out = nothing
    H   = nothing

    while !converged
        begin_time = time_ns()
        out, H = optimize!(m, data; method = method,
                           ftol = ftol, grtol = gtol, xtol = xtol,
                           iterations = n_iterations, show_trace =
                           (verbose == :none ? false : true), step_size = step_size,
                           verbose = verbose,
                           mle = mle)

        attempts += 1
        total_iterations += out.iterations
        converged = out.converged || attempts > max_attempts

        end_time = (time_ns() - begin_time)/1e9
        println(verbose, :low, @sprintf "Total iterations completed: %d\n" total_iterations)
        println(verbose, :low,
                @sprintf "Optimization time elapsed: %5.2f\n" optimization_time += end_time)

        # Write params to file after every `n_iterations` iterations
        params = map(θ -> θ.value, get_parameters(m))
        if save_results
            h5open(rawpath(m, "estimate", "paramsmode.h5"),"w") do file
                file["params"] = params
            end
        end
    end

    # write parameters to file one last time so we have the final mode
    if save_results
        h5open(rawpath(m, "estimate", lowercase(string(get_setting(m, :sampling_method))) *
                       "_" * "paramsmode.h5"),"w") do file
                           file["params"] = params
                       end
    end

    if get_all_results
        return params, out, H
    else
        return params
    end
end

"""
```
function hessian_step!(m::Union{AbstractDSGEModel, AbstractVARModel},
                       data::Matrix{S}; verbose::Symbol = :none,
                       check_neg_diag::Bool = true) where {S <: Real}
```
calculates the Hessian at the current calibration of the model object
and checks that it is positive definite.
"""
function hessian_step!(m::Union{AbstractDSGEModel, AbstractVARModel},
                       data::Matrix{S}; verbose::Symbol = :none,
                       check_neg_diag::Bool = true) where {S <: Real}

    params = map(θ -> θ.value, get_parameters(m))
    hessian, _ = hessian!(m, params, data; verbose = verbose, check_neg_diag = check_neg_diag)

    return hessian
end


text:
"""
This is a main module file. In it, you should:

- `export` the desired functions from your `include()` calls (comma-delimited)

- `include()` your `srcfile`s. There is a stubbed, vectorized `include()` here,
  but if you want to change it, just construct a manual array of source file
  names.
"""

module PkgName

export sharpe

dir = readdir(dirname(@__FILE__))
includes = filter(x -> (endswith(x, "jl") && x != "PkgName.jl"), dir)
include.(includes)

end # module PkgName

max_stars_repo_name: ryapric/boilerplate

max_issues_repo_name: ryapric/boilerplate

max_forks_repo_name: ryapric/boilerplate


text:
function MatrixMultiplyOp(in,adj;matrix=1)

	if (adj)
		out = matrix'*in
	else
		out = matrix*in
	end

	return out
end


text:
export codemirror

@api codemirror => (CodeMirror <: Widget) begin
    doc("Create a code viewer/editor")
    arg(code::AbstractString="", doc="The code to display.")
    kwarg(readonly::Bool=false, doc="If set to true, editing will be disabled.")
    kwarg(language::AbstractString="julia", doc="The language used for syntax highlighting.")
    kwarg(
        theme::AbstractString="elegant",
        doc=md"""The theme. Valid values are `"ambiance"`,`"ambiance-mobile"`,
                 `"elegant"`, `"monokai"`, `"solarized"` and `"twilight"`."""
    )
    kwarg(linenumbers::Bool=true, doc="If set to true, line numbers will be shown.")
    kwarg(tabsize::Int=4, doc="The tab size.")
end

wrapbehavior(c::CodeMirror) =
    hasstate(c, attr="immediateValue", trigger="change")

# Render to virtual DOM
render(c::CodeMirror, state) =
    Elem("code-mirror",
        attributes = @d(
            :value=>c.code,
            :mode=>c.language,
            "read-only"=>boolattr(c.readonly),
            :theme=>c.theme,
            "line-numbers"=>boolattr(c.linenumbers),
        )
    )