text:
function minimize(cur_reactor, cur_field, init_reactor=nothing)
  cur_count = 0

  if is_present(init_reactor) && cur_reactor.constraint != init_reactor.constraint
    tmp_value = getfield(init_reactor, cur_field)

    tmp_reactor = deepcopy(cur_reactor)

    if tmp_reactor.is_consistent
      tmp_reactor = hone(tmp_reactor, cur_field, tmp_value)
    else
      tmp_reactors = match(tmp_reactor, cur_field, tmp_value, false)
      tmp_reactor = isempty(tmp_reactors) ? nothing : first(tmp_reactors)
    end

    cur_count += 1

    is_nothing(tmp_reactor) && return
    tmp_reactor.is_valid || return
    tmp_reactor.is_good || return
  end

  cur_T_bar_list = bisect_reorder!(
    collect(linspace(min_T_bar, max_T_bar, 31))
  )

  for cur_T_bar in cur_T_bar_list
    is_present(init_reactor) && break

    tmp_reactor = deepcopy(cur_reactor)
    tmp_reactor.T_bar = cur_T_bar

    solve!(tmp_reactor)
    cur_count += 1

    tmp_reactor.is_valid || continue
    tmp_reactor.is_good || continue

    init_reactor = tmp_reactor
  end

  is_nothing(init_reactor) && return

  min_reactors = [init_reactor]

  good_value = getfield(init_reactor, cur_field)
  bad_value = NaN

  while isnan(bad_value)
    tmp_value = good_value / 2

    tmp_reactor = deepcopy(cur_reactor)

    if tmp_reactor.is_consistent
      tmp_reactor = hone(tmp_reactor, cur_field, tmp_value)
      tmp_reactors = [tmp_reactor]
    else
      tmp_reactors = match(tmp_reactor, cur_field, tmp_value, false)
      tmp_reactor = isempty(tmp_reactors) ? nothing : first(tmp_reactors)
    end

    cur_count += 1

    if !is_nothing(tmp_reactor) && tmp_reactor.is_valid && tmp_reactor.is_good
      min_reactors = tmp_reactors
      good_value = tmp_value
    else
      bad_value = tmp_value
    end
  end

  while ( abs( good_value - bad_value ) / good_value ) > 1e-8
    tmp_value = middle(good_value, bad_value)

    tmp_reactor = deepcopy(cur_reactor)

    if tmp_reactor.is_consistent
      tmp_reactor = hone(tmp_reactor, cur_field, tmp_value)
      tmp_reactors = [tmp_reactor]
    else
      tmp_reactors = match(tmp_reactor, cur_field, tmp_value, false)
      tmp_reactor = isempty(tmp_reactors) ? nothing : first(tmp_reactors)
    end

    cur_count += 1

    if !is_nothing(tmp_reactor) && tmp_reactor.is_valid && tmp_reactor.is_good
      min_reactors = tmp_reactors
      good_value = tmp_value
    else
      bad_value = tmp_value
    end
  end

  ( length(min_reactors) > 1 ) && custom_log("Minimize Bug...")
  ( length(min_reactors) > 1 ) && custom_log(min_reactors)

  min_reactors[1]
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 81af4446-7384-11eb-2b09-4916d5e47ee0
begin

	using Dates 		# For timestamps
	using DotEnv 		# For loading the OpenAI API key
	using PyCall 		# For interacting with OpenAI's Python API
	using PlutoUI 		# For embedding the tweet above, and other fun visuals

end

# ╔═╡ 95ac84ac-7381-11eb-3087-6bdfce67861a
md"""
# Meet Lord Robot Andy!

🤖 💻 ✍ 📔 🎭 🎫 😺

"""

# ╔═╡ 2e148102-7384-11eb-3683-63f615cfeaeb
html"""

1. Did you see Cats in theaters? Have you been craving more content from Andrew Lloyd Webber? No? Just @CaseyReiland and @_natmiller? Luckily for them, more content is just a few lines of code away! Introducing: Lord (Robot) Andy! https://t.co/aMd6g6IXCA
— Joe(y) Carpinelli (@code_lowered) February 8, 2021

"""

# ╔═╡ 78d4a6f6-7382-11eb-2568-ed3cae869fed
md"""
## Introduction
_An FAQ of sorts to introduce this project, and all it stands for._

#### Why, God why?
After watching Cats several times, two wonderful people asked for _even more_ Andrew Lloyd Webber content. Who are we to deny them?

#### How does Robot Andy work?

A stiff cup of coffee and some soft lofi music.

#### What software tools are used?
There are lots of free AI text generators out today, including Tensorflow and [textgenrnn](https://github.com/minimaxir/textgenrnn), and [OpenAI's](https://openai.com) GPT series. Some are used here — these tools are identified as they are used. See the end of this document for dependency setup and imports.

"""

# ╔═╡ 5a4acfc0-7390-11eb-0ebb-7b3fdde32b18
md"""
## Love Never (Ever) Dies
_A script written by Lord Andy themself._
"""

# ╔═╡ 7c649da6-7384-11eb-1c64-0dc21d83e13f
md"""
## Project Dependencies
_Let's load the Julia packages and other software tools that we need._
"""

# ╔═╡ 25f8a6b4-7385-11eb-0433-913d70041f71
md"""
#### Julia Packages
"""

# ╔═╡ 2d3c7bf8-7385-11eb-3eac-bb9a2642a75d
md"""
#### OpenAI API Setup
"""

# ╔═╡ e8eeb524-7384-11eb-06ba-2f24907cfc8d
begin

	# Top level of Git repository
	TOPLEVEL = string(read(`git rev-parse --show-toplevel`, String)[1:end-1])

	# Collect variables for .env file (not included in repo)
	DotEnv.config(path = joinpath(TOPLEVEL,".env"))
	ORG = ENV["ORG_ID"]
	KEY = ENV["API_KEY"]

	# Include OpenAI's Python API
	openai = pyimport("openai")
	openai.api_key = KEY

	# Create a logs directory (if needed)
	LOGSDIR = joinpath(TOPLEVEL, "logs")
	isdir(LOGSDIR) || mkdir(LOGSDIR)

end;

# ╔═╡ 0447062e-7386-11eb-0e1e-89c469497206
md"""
#### Helper Functions
"""

# ╔═╡ ec24f400-738a-11eb-1e56-ad0bbea7accf
md"""
###### Get a timestamp (yyyy-mm-dd-Thh-mm-ss-ms)
"""

# ╔═╡ 066fae54-738b-11eb-14b2-034e13166d6c
function timestamp(dt=now())

	return string(dt) |>
			str->replace(str, ":"=>"-") |>
			str->replace(str, "."=>"-") |>
			str->replace(str, "T"=>"-")

end;

# ╔═╡ 1aeba388-7386-11eb-1232-550d54e919ac
md"""
###### Give a Prompt, Get a Response
"""

# ╔═╡ 3153cec4-7391-11eb-03e4-a7c12786b9e6
md"""
###### Word Count
"""

# ╔═╡ 34aceb82-7391-11eb-20eb-63136247f19e
function wordcount(str)
	return length(split(str, " "))
end;

# ╔═╡ 819f260c-7388-11eb-3de0-a1a8694f7932
function GPT3(prompt::String;
			  engine 		=	"davinci",
			  max_tokens 	= 	wordcount(prompt),
			  log 			= 	true,
	 		  filename 		=   joinpath(LOGSDIR, "GPT-3" * timestamp() * ".txt"),
		      temperature   =   0.9)

	response = openai.Completion.create(
		engine 		=	engine,
		prompt 		=	prompt,
		max_tokens 	= 	max_tokens,
		n 			= 	1,
		temperature =   0.9)

	response_text = string(first(response["choices"])["text"])
	if log
		open(filename, "w") do io
			write(io,
				"""
				This file was autogenerated with GPT-3! All content
				is checked for obscenities, and no prompts should
				produce any inappropriate content. Still, if you
				find ANY offensive content in this file, please
				contact the project's maintainer (@cadojo on GitHub).

				This file has the following format: <prompt>... <response>

				---

				$prompt...$response_text

				"""
			)
	   	end
	end

	return response_text

end;

# ╔═╡ 797c2706-7390-11eb-3c3d-d5ac477ff030
LNED = let

	LND = readlines(joinpath(TOPLEVEL, "scripts", "love-never-dies.txt"))
	LNED = Vector{String}()

	# Add opening line...
	push!(LNED, GPT3(first(LND)))

	# Find the ~next~ line which completes a sentence
	prompt_end = 2
	for line_id ∈ 2:length(LND)
		if any(map(ch->contains(LND[line_id], ch), (".","!","?")))
			prompt_end = line_id
			break
		end
	end

	# Collect the prompt
	prompt = reduce((x,y)->string(x,"\n",y), collect(LND[2:prompt_end]))

	# Generate a script!
	push!(LNED, prompt * GPT3(prompt * "\n\n"; max_tokens = 1800))

	LNED

end;

# ╔═╡ 8429e00e-7392-11eb-147c-879efb564798
with_terminal(print, string(LNED...))

# ╔═╡ Cell order:
# ╟─95ac84ac-7381-11eb-3087-6bdfce67861a
# ╟─2e148102-7384-11eb-3683-63f615cfeaeb
# ╟─78d4a6f6-7382-11eb-2568-ed3cae869fed
# ╟─5a4acfc0-7390-11eb-0ebb-7b3fdde32b18
# ╠═797c2706-7390-11eb-3c3d-d5ac477ff030
# ╠═8429e00e-7392-11eb-147c-879efb564798
# ╟─7c649da6-7384-11eb-1c64-0dc21d83e13f
# ╟─25f8a6b4-7385-11eb-0433-913d70041f71
# ╠═81af4446-7384-11eb-2b09-4916d5e47ee0
# ╟─2d3c7bf8-7385-11eb-3eac-bb9a2642a75d
# ╠═e8eeb524-7384-11eb-06ba-2f24907cfc8d
# ╟─0447062e-7386-11eb-0e1e-89c469497206
# ╟─ec24f400-738a-11eb-1e56-ad0bbea7accf
# ╠═066fae54-738b-11eb-14b2-034e13166d6c
# ╟─1aeba388-7386-11eb-1232-550d54e919ac
# ╠═819f260c-7388-11eb-3de0-a1a8694f7932
# ╟─3153cec4-7391-11eb-03e4-a7c12786b9e6
# ╠═34aceb82-7391-11eb-20eb-63136247f19e

ext: jl

max_stars_repo_path: src/Audacious Screenplay.jl

max_stars_repo_name: cadojo/Lord-Robot-Andy

max_issues_repo_path: src/Audacious Screenplay.jl

max_issues_repo_name: cadojo/Lord-Robot-Andy

max_forks_repo_path: src/Audacious Screenplay.jl

max_forks_repo_name: cadojo/Lord-Robot-Andy


text:
export TraceContextTextMapPropagator

"""
    TraceContextTextMapPropagator()

This propagator follows the [W3C format](https://www.w3.org/TR/trace-context/#trace-context-http-headers-format)
"""
struct TraceContextTextMapPropagator <: AbstractPropagator end

function inject!(
    carrier::Union{
        AbstractVector{<:Pair{<:AbstractString,<:AbstractString}},
        AbstractDict{<:AbstractString,<:AbstractString},
    },
    propagator::TraceContextTextMapPropagator,
    ctx::Context = current_context(),
)
    sc = span_context(ctx)
    if !isnothing(sc)
        s_trace_parent = "00-$(string(sc.trace_id, pad=32))-$(string(sc.span_id,pad=16))-$(sc.trace_flag.sampled ? "01" : "00")"
        push!(carrier, "traceparent" => s_trace_parent)
        s_trace_state = string(sc.trace_state)
        if !isempty(s_trace_state)
            push!(carrier, "tracestate" => s_trace_state)
        end
    end
    carrier
end

# fallback
function inject!(
    carrier::T,
    ::TraceContextTextMapPropagator,
    ctx::Context = current_context(),
) where {T}
    @warn "unknown carrier type $T"
    carrier
end

TRACEPARENT_HEADER_FORMAT =
    r"^[ \t]*(?P<version>[0-9a-f]{2})-(?P<trace_id>[0-9a-f]{32})-(?P<span_id>[0-9a-f]{16})-(?P<trace_flag>[0-9a-f]{2})(?P<rest>-.*)?[ \t]*$"

function extract(
    carrier::Union{
        AbstractVector{<:Pair{<:AbstractString,<:AbstractString}},
        AbstractDict{<:AbstractString,<:AbstractString},
    },
    propagator::TraceContextTextMapPropagator,
    ctx::Context = current_context(),
)
    trace_id, span_id, trace_flag, trace_state = nothing, nothing, nothing, TraceState()
    for (k, v) in carrier
        if k == "traceparent"
            m = match(TRACEPARENT_HEADER_FORMAT, v)
            if !isnothing(m)
                if m["version"] == "00" && isnothing(m["rest"])
                    trace_id = parse(TraceIdType, m["trace_id"])
                    span_id = parse(SpanIdType, m["span_id"])
                    trace_flag = parse(TraceFlag, m["trace_flag"])
                end
            end
        elseif k == "tracestate"
            trace_state = parse(TraceState, v)
        end
    end
    if isnothing(trace_id) || isnothing(span_id) || isnothing(trace_flag)
        ctx
    else
        merge(
            ctx,
            Context((;
                SPAN_KEY_IN_CONTEXT => NonRecordingSpan(
                    "",
                    SpanContext(span_id, trace_id, false, trace_flag, trace_state),
                    nothing,
                )
            )),
        )
    end
end

# fallback
function extract(
    carrier::T,
    ::TraceContextTextMapPropagator,
    ctx::Context = current_context(),
) where {T}
    @warn "unknown carrier type $T"
    ctx
end


text:
# 本代码主要是针对此项目制定一个通用命名标准

# 桨叶的物理描述
# 本项目中旋翼转向为逆时针旋转，也就是俗称右旋旋翼
# 旋翼一般都是沿径向扭转的，且在实际分析中多为负扭转，Δθ=θtw*r，常规直升机一般为负扭转
R = 10  #the rotor radius; the length of the blade, measured from the rotation to tip
Omega = 66  #the rotor ratational speed or angular velocity
rho = 1.25  #air density
psi = 0 #azimuth angle of the blade; for constant rotational speed, ψ = Ω*t
r = 0.1 #radial location on the blade, measured from the center of ratation(r = 0)
ch = 0.1 #blade chord
NR = 2 #number of blades
Nbe = 10 #number of blade elements
m = 0.3 #blade mass per unit length
Ib = 10 #∫₀ᴿ(r^2*m); moment of inertia of the blade about the center of rotation
A = π*R^2 #rotor disk area
sigma = Nch/(π*R) #rotor solidity
gama = rho*cla*ch*R^4/Ib #blade lock number (洛克数代表了气动力和惯性力的比值)

# 桨叶气动参数
cla = 0.6 #blade section two-dimensional lift curve slope
alpha = 0.5 #blade section angle-of-attack
Ma = 0.3 #blade section Mach number

# 桨叶运动参数
# 在本项目中，桨叶运动被认为是绕桨毂的刚体运动
beta = 0.5 #blade flap angle
zeta = 0.5 #blade lag angle
theta = 0.5 #blade pitch angle
thetatw = -5 #blade twist
beta = beta0+betalon*cos(psi)+betalat*sin(psi)
theta = theta0+thetalat*cos(psi)+thetalon*cos(psi)+thetatw*r/R

# 旋翼迎角和速度
# 下标s表示旋翼固定坐标系
alphai = 0.5 #rotor disk plane incidence angle(低头为正)
vair = 20 #rotor or helicopter velocity with respect to the air
vind_s = 10 #rotor induced velocity, normal to the disk plane
mu = vair*cos(alphai)/(Omega*R) #rotor advance ratio
muz = vair*sin(alphai)/(Omega*R) #normal velocity ratio
lambda = (vair*sin(alphai)+vind_s)/(Omega*R) #rotor inflow ratio
lambdai = vind_s/(Omega*R) #induced inflow ratio

# 旋翼力和力矩
T = 100 #rotor thrust
H = 10 #rotor drag force in the disk plane
Y = 1 #rotor side force in the disk plane
Q = 5 #rotor shaft torque
P = 10 #rotor shaft power(P=Ω*Q)
CT = T/(rho*A*(Omega*R)^2)
CH = H/(rho*A*(Omega*R)^2)
CY = Y/(rho*A*(Omega*R)^2)
CQ = Q/(rho*A*(Omega*R)^2*R)
CP = P/(rho*A*(Omega*R)^3) #CP=CQ

#  旋翼桨盘平面
# tpp tip-path plane
# nfp no-feathering plane
# hp hub plane
# cp control plane
#
#  各坐标系下标
# s system(固定坐标系)
# r rotation(旋转坐标系)
# β beta(挥舞坐标系)
# be blade element(叶素当地坐标系)

# 其他下标说明
# all means total
# x means x coordination
# y means y coordination
# z means z coordination
# aero表示该项是由气动导致的
# ui表示该项是均匀入流
# ind代表诱导速度
# hover代表悬停状态
# tl代表tip loss

max_stars_repo_name: dingzhw/RotorTrim

max_issues_repo_name: helicopterliu/RotorTrim

max_forks_repo_name: helicopterliu/RotorTrim


text:
export LogZero
export LogSum2Exp
export LogSum3Exp
export LogSumExp

LogZero(T::DataType) = - floatmax(T)

"""
    LogSum2Exp(a::Real, b::Real) -> max(a,b) + log(1.0 + exp(-abs(a-b)))
```julia
julia> LogSum2Exp(Float32(1.2),Float64(3.3))
3.4155195283818967
```
"""
function LogSum2Exp(a::Real, b::Real)
    isinf(a) && return b
    isinf(b) && return a
    if a < b
        a, b = b, a
    end
    return (a + log(1.0 + exp(b-a)))
end

function LogSum3Exp(a::Real, b::Real, c::Real)
    return LogSum2Exp(LogSum2Exp(a,b),c)
end

function LogSumExp(a)
    tmp = LogZero(eltype(a))
    for i = 1:length(a)
        tmp = LogSum2Exp(tmp, a[i])
    end
    return tmp
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
# Copyright (c) 2015: AmplNLWriter.jl contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

import AmplNLWriter
import MINLPTests
using Test

const TERMINATION_TARGET = Dict(
    MINLPTests.FEASIBLE_PROBLEM => AmplNLWriter.MOI.LOCALLY_SOLVED,
    MINLPTests.INFEASIBLE_PROBLEM => AmplNLWriter.MOI.INFEASIBLE,
)

const PRIMAL_TARGET = Dict(
    MINLPTests.FEASIBLE_PROBLEM => AmplNLWriter.MOI.FEASIBLE_POINT,
    MINLPTests.INFEASIBLE_PROBLEM => AmplNLWriter.MOI.NO_SOLUTION,
)

# Common reasons for exclusion:
# nlp/005_011     : Uses the function `\`
# nlp/006_010     : Uses a user-defined function
# nlp/007_010     : Ipopt returns an infeasible point, not NO_SOLUTION.
# nlp/008_010     : Couenne fails to converge
# nlp/008_011     : Couenne fails to converge
# nlp-cvx/109_010 : Ipopt fails to converge
# nlp-cvx/206_010 : Couenne can't evaluate pow
# nlp-mi/001_010  : Couenne fails to converge

const CONFIG = Dict{String,Any}()

import Bonmin_jll
CONFIG["Bonmin"] = Dict(
    "amplexe" => Bonmin_jll.amplexe,
    "options" => String["bonmin.nlp_log_level=0"],
    "tol" => 1e-5,
    "dual_tol" => NaN,
    "nlp_exclude" => ["005_011", "006_010"],
    "nlpcvx_exclude" => ["109_010"],
    # 004_010 and 004_011 are tolerance failures on Bonmin
    "nlpmi_exclude" => ["004_010", "004_011", "005_011", "006_010"],
    "infeasible_point" => AmplNLWriter.MOI.NO_SOLUTION,
)

import Couenne_jll
CONFIG["Couenne"] = Dict(
    "amplexe" => Couenne_jll.amplexe,
    "options" => String[],
    "tol" => 1e-2,
    "dual_tol" => NaN,
    "nlp_exclude" => ["005_011", "006_010", "008_010", "008_011"],
    "nlpcvx_exclude" => ["109_010", "206_010"],
    "nlpmi_exclude" => ["001_010", "005_011", "006_010"],
    "infeasible_point" => AmplNLWriter.MOI.NO_SOLUTION,
)

import Ipopt_jll
CONFIG["Ipopt"] = Dict(
    "amplexe" => Ipopt_jll.amplexe,
    "options" => String["print_level=0"],
    "tol" => 1e-5,
    "dual_tol" => 1e-5,
    "nlp_exclude" => ["005_011", "006_010", "007_010"],
    "nlpcvx_exclude" => ["109_010"],
    "nlpmi_exclude" => ["005_011", "006_010"],
    "infeasible_point" => AmplNLWriter.MOI.NO_SOLUTION,
)

# SHOT fails too many tests to recommend using it.
# e.g., https://github.com/coin-or/SHOT/issues/134
# Even problems such as `@variable(model, x); @objective(model, Min, (x-1)^2)`
#
# import SHOT_jll
# CONFIG["SHOT"] = Dict(
#     "amplexe" => SHOT_jll.amplexe,
#     "options" => String[
#         "Output.Console.LogLevel=6",
#         "Output.File.LogLevel=6",
#         "Termination.ObjectiveGap.Absolute=1e-6",
#         "Termination.ObjectiveGap.Relative=1e-6",
#     ],
#     "tol" => 1e-2,
#     "dual_tol" => NaN,
#     "nlp_exclude" => [
#         "005_011",  # `\` function
#         "006_010",  # User-defined function
#     ],
#     "nlpcvx_exclude" => [
#         "501_011",  # `\` function
#     ],
#     "nlpmi_exclude" => [
#         "005_011",  # `\` function
#         "006_010",  # User-defined function
#     ],
#     "infeasible_point" => AmplNLWriter.MOI.UNKNOWN_RESULT_STATUS,
# )

@testset "$(name)" for name in ["Ipopt", "Bonmin", "Couenne"]
    config = CONFIG[name]
    OPTIMIZER =
        () -> AmplNLWriter.Optimizer(config["amplexe"], config["options"])
    PRIMAL_TARGET[MINLPTests.INFEASIBLE_PROBLEM] = config["infeasible_point"]
    @testset "NLP" begin
        MINLPTests.test_nlp(
            OPTIMIZER,
            exclude = config["nlp_exclude"],
            termination_target = TERMINATION_TARGET,
            primal_target = PRIMAL_TARGET,
            objective_tol = config["tol"],
            primal_tol = config["tol"],
            dual_tol = config["dual_tol"],
        )
    end
    @testset "NLP-CVX" begin
        MINLPTests.test_nlp_cvx(
            OPTIMIZER,
            exclude = config["nlpcvx_exclude"],
            termination_target = TERMINATION_TARGET,
            primal_target = PRIMAL_TARGET,
            objective_tol = config["tol"],
            primal_tol = config["tol"],
            dual_tol = config["dual_tol"],
        )
    end
    if name != "Ipopt"
        @testset "NLP-MI" begin
            MINLPTests.test_nlp_mi(
                OPTIMIZER,
                exclude = config["nlpmi_exclude"],
                termination_target = TERMINATION_TARGET,
                primal_target = PRIMAL_TARGET,
                objective_tol = config["tol"],
                primal_tol = config["tol"],
                dual_tol = config["dual_tol"],
            )
        end
    end
end


text:
function test_membership(s::MinHash.HashSet, positives::Set, superset)
    mistakes = false
    for i in superset
        mistakes |= ((i in positives) ⊻ (i in s))
    end
    @test !mistakes
end

@testset "Instantiation" begin
    function test_hashset_length(s::MinHash.HashSet)
        @test s.mask + 1 == length(s.data)
        @test count_ones(length(s.data)) == 1
        return length(s.data)
    end

    for len in [1, 10, 16, 17, 100]
        @test test_hashset_length(MinHash.HashSet(len)) ≥ len
    end
end

@testset "Basic functionality" begin
    for len in [10, 100, 500, 5000]
        rnge = UInt(1):UInt(len)
        rf = Set(rand(rnge, div(len, 3)))
        s = MinHash.HashSet(len)
        for i in rf
            MinHash.unsafe_push!(s, i)
        end
        test_membership(s, rf, rnge)
        @test length(s) == length(rf)
    end
end

@testset "Repopulation" begin
    rnge = UInt(1):UInt(1000)
    integers = collect(Set(rand(rnge, 400)))
    keep = Set(integers[1:100])
    heap = heapify!(collect(keep), Base.Order.Reverse)
    s = MinHash.HashSet(1000)
    for i in integers
        MinHash.unsafe_push!(s, i)
    end
    MinHash.repopulate!(s, heap)
    test_membership(s, keep, integers)
    @test length(s) == length(keep)

    # Here, test the second if statement in push!
    s = MinHash.HashSet(1000)
    for i in UInt(1):UInt(500)
        MinHash.unsafe_push!(s, i)
    end
    heap = UInt[]
    for i in UInt(501):UInt(1000)
        heappush!(heap, i, Base.Order.Reverse)
        push!(s, i, heap)
    end
end

@testset "Pushing" begin
    s = MinHash.HashSet(50)
    heap = heapify!(rand(UInt, 25), Base.Order.Reverse)
    for i in 1:1000
        largest = pop!(heap)
        smaller = largest - 10
        heappush!(heap, smaller, Base.Order.Reverse)
        push!(s, smaller, heap)
    end
    @test all(i in s for i in heap)
end


text:
module PeaceCypher

using Random
using CryptoGroups
using CryptoSignatures
using DiffieHellman
using SecureIO
using Nettle
using Pkg.TOML

### I could promote Notary and CypherSuite as abstract types.

struct ID#{T}
    id::BigInt # ::T
end

function Base.string(id::ID; length=nothing, kwargs...)
    str = string(id.id; kwargs...)
    if length!=nothing
        len = Base.length(str)
        @assert len<=length
        pre = "0"^(length-len)
    else
        pre = ""
    end
    return "$pre$str"
end

ID(str::AbstractString; kwargs...) = ID(parse(BigInt,str; kwargs...))

Base.Dict(id::ID) = Dict("id"=>string(id,base=16))
ID(dict::Dict) = ID(dict["id"],base=16)

Base.Vector{UInt8}(id::ID; kwargs...) = Vector{UInt8}(string(id; kwargs...))
ID(bytes::Vector{UInt8}; kwargs...) = ID(String(copy(bytes)); kwargs...)

Base.:(==)(a::ID,b::ID) = a.id==b.id
Base.hash(a::ID,h::UInt) = Base.hash(a.id, Base.hash(:ID,h))
Base.in(a::ID,b::ID) = a==b

abstract type Notary end
abstract type CypherSuite end
abstract type Layer end

struct CryptoNotary <: Notary
    G::AbstractGroup
end

const Key = CryptoSignatures.Signer

struct Signer
    key::Key
    notary::Notary
end

# This should allow to improve notation in the future if necessary
struct Signature
    sig::Dict
    notary::Notary
end

Base.Dict(s::Signature) = s.sig

function binary(s::Signature)
    io = IOBuffer()
    TOML.print(io,Dict(s))
    return take!(io)
end

function Signature(binary::Vector{UInt8}, notary::Notary)
    dict = TOML.parse(String(copy(binary)))
    return Signature(dict, notary)
end

### In future Notary shall accept configuration intialize it from that
Notary() = CryptoNotary(CryptoGroups.Scep256k1Group())

### Some arbitrary defaults
hash(x::AbstractString, crypto::CryptoNotary) = BigInt(Base.hash(x))
hash(x, crypto::Notary) = hash(string(x), crypto)

newsigner(crypto::CryptoNotary) = Signer(Key(crypto.G),crypto)

function sign(value, notary::CryptoNotary, signer::Key)
    signature = DSASignature(hash(value, notary), signer)
    signaturedict = Dict(signature)
    return Signature(signaturedict, notary)
end

### Perhaps I could use certify
sign(value, signer::Signer) = sign(value, signer.notary, signer.key)

# In a real life we are willing to sacrifice some conviniance for security in fear of bugs. That could be achevead by typing the value and adding validation checks with methods for CryptoNotary. For example the program asks to sign DHParameter(a, G^a) where the sign method first would validate G^a=G^a. Similalry a validation check could be done on the Braid, but that must be done in the context of braid method.

function verify(value, signature::Dict, crypto::CryptoNotary)
    signature = DSASignature{BigInt}(signature)
    return CryptoSignatures.verify(signature,crypto.G) && signature.hash==hash(value, crypto)
end

verify(value, s::Signature) = verify(value, s.sig, s.notary)

id(s::Dict, crypto::CryptoNotary) = ID(hash(string(parse(BigInt,s["pub"],base=16)), crypto))
id(s::Signature) = id(s.sig, s.notary)
id(s::Signer) = ID(hash("$(s.key.pubkey)", s.notary))

struct DiffieHellmanMerkle <: CypherSuite
    notary::Notary
end

### Again in future the correct cypher suite would be choosen uppon the passed configuration in Dict form which one could easally parse from a TOML file.
CypherSuite(notary::Notary) = DiffieHellmanMerkle(notary)

function fold(value::BigInt, signature::Dict)
    dict = Dict("value"=>string(value, base=16),"signature"=>signature)
    io = IOBuffer()
    TOML.print(io,dict)
    return take!(io)
end

function unfold(envelope::Vector{UInt8})
    dict = TOML.parse(String(copy(envelope)))
    value = parse(BigInt,dict["value"], base=16)
    signature = dict["signature"]
    return value, signature
end

function seal(value::BigInt, signer::Signer)
    signature = sign(value, signer)
    return fold(value, Dict(signature))
end

function unseal(envelope::Vector{UInt8}, crypto::Notary)
    value, sdict = unfold(envelope)
    signature = Signature(sdict, crypto)

    if verify(value, signature) ### The id is of ID type thus there shall be no problem
        return value, id(signature)
    else
        return value, nothing
    end
end

function seal(value::BigInt)
    str = string(value, base=16)
    return Vector{UInt8}(str)
end

function unseal(envelope::Vector{UInt8})
    value = parse(BigInt,String(copy(envelope)),base=16)
    return value, nothing
end

function rngint(len::Integer)
    max_n = ( BigInt(1) << len ) - 1
    if len > 2
        min_n = BigInt(1) << (len - 1)
        return rand(min_n:max_n)
    end
    return rand(1:max_n)
end

# Where did I need this hack?
import Base.in
in(x::Nothing,y::Nothing) = true

function secure(socket::IO, crypto::DiffieHellmanMerkle, signer::Signer, id)
    G = crypto.notary.G

    dh = DH(value->seal(value,signer),envelope->unseal(envelope,crypto.notary),G,(x,y,z)->seal(hash("$x $y $z",crypto.notary)),() -> rngint(100))

    key, idr = diffiehellman(socket, dh)
    @assert idr in id "$idr not in $id"
    securesocket = SecureSocket(socket, key)
    return securesocket
end

function secure(socket::IO, crypto::DiffieHellmanMerkle, signer::Signer)
    G = crypto.notary.G

    dh = DH(value->seal(value,signer),envelope->unseal(envelope),G,(x,y,z)->seal(hash("$x $y $z",crypto.notary)),() -> rngint(100))

    key, idr = diffiehellman(socket, dh)
    securesocket = SecureSocket(socket, key)
    return securesocket
end

function secure(socket::IO, crypto::DiffieHellmanMerkle, id)
    G = crypto.notary.G

    dh = DH(value->seal(value),envelope->unseal(envelope,crypto.notary),G,(x,y,z)->seal(hash("$x $y $z",crypto.notary)),() -> rngint(100))

    key, idr = diffiehellman(socket, dh)
    @assert idr in id "$idr not in $id"
    securesocket = SecureSocket(socket, key)
    return securesocket
end

abstract type KeyParameter end

mutable struct DHM <: KeyParameter
    a
    # tbegin
    # tend
    A
end

mutable struct Secret
    cs::CypherSuite
    par::KeyParameter
    signature::Signature
end

### This is a DIffie-Hellman-Merkel authentificated secret

# This method shall generate a secret for the cypher suite which, for example, could be done on aproval basis occasionly.
function Secret(cs::CypherSuite, signer::Signer)
    return nothing
end

# This method updates the layer with
function update!(secret::Secret, signer::Signer) ### Could use parametric types to be more precise
    return nothing
end

function update!(secret::Secret, par::KeyParameter, signature::Signature)
    secret.par = par
    secret.signature = signature
end

### It is undesirable to expose the signer with in the protocol as that could in principle could be used to perform a timing attack. Instead we could generate a secret with authetification like
# secret = generatesecret(cs::CypherSuite)

abstract type SymmetricLayer <: Layer end
abstract type AsymmetricLayer <: Layer end

struct SecureLayerSymmetric <: SymmetricLayer
    crypto::CypherSuite
    signer::Union{Signer,Secret}
    id
end

struct SecureLayerMaster <: AsymmetricLayer
    crypto::CypherSuite
    signer::Union{Signer,Secret}
end

struct SecureLayerSlave <: AsymmetricLayer
    crypto::CypherSuite
    id
end

Layer(crypto::CypherSuite, signer::Signer, id) = SecureLayerSymmetric(crypto, signer, id)
Layer(crypto::CypherSuite, signer::Signer) = SecureLayerMaster(crypto, signer)
Layer(crypto::CypherSuite, id) = SecureLayerSlave(crypto, id)

secure(socket::IO, sc::SecureLayerSymmetric) = secure(socket, sc.crypto, sc.signer, sc.id)
secure(socket::IO, sc::SecureLayerMaster) = secure(socket, sc.crypto, sc.signer)
secure(socket::IO, sc::SecureLayerSlave) = secure(socket, sc.crypto, sc.id)

export Notary, hash, verify, newsigner, id, ID
export Signer, sign
export Signature, binary, verify, id # Dict
export CypherSuite, Layer, secure

end # module


text:
"""
Kernel
"""
abstract type KernelBase end

struct LogisticKernel <: KernelBase
    o::PyObject
end

function LogisticKernel(lambda::Real)
    return LogisticKernel(sparse_ir.LogisticKernel(Float64(lambda)))
end

struct RegularizedBoseKernel <: KernelBase
    o::PyObject
end

function RegularizedBoseKernel(lambda::Real)
    return RegularizedBoseKernel(sparse_ir.RegularizedBoseKernel(Float64(lambda)))
end

(k::KernelBase)(x::Float64, y::Float64) = k.o(x, y)


text:
module TestCTree

using Copulas
using Base.Test

##################
## constructors ##
##################

## CTreePaths
##------------

paths = Array{Int, 1}[[1, 4], [1, 2, 3], [5, 6], [5, 7, 8]]

tP = Copulas.CTreePaths(9, paths)
tP1 = Copulas.CTreePaths(9, [1, 4], [1, 2, 3], [5, 6], [5, 7, 8])

@test tP == tP1
@test Copulas.dg(tP)[3] == [5, 6]

## sorting paths
paths = Array{Int, 1}[[1, 4], [1, 2, 3], [5, 6], [5, 7, 8]]
sortedPaths = Array{Int, 1}[[1, 2, 3], [1, 4], [5, 6], [5, 7, 8]]
@test sortedPaths == Copulas.sortPaths(paths)

## CTreeParRef
##------------

## test conversion
tPar = convert(Copulas.CTreeParRef, tP)
parNot = [9, 1, 2, 1, 9, 5, 5, 7, 0]
@test Copulas.dg(tPar) == parNot

tPar2 = Copulas.CTreeParRef(parNot)
@test tPar == tPar2

## additional conversion tests
##----------------------------

tPar = Copulas.CTreeParRef([9, 1, 2, 1, 9, -1, 5, 7, 0])
tPaths = Copulas.CTreePaths(9, [1, 4], [1, 2, 3], [5, 7, 8])
@test !(tPar == tPaths)
tPar2 = convert(Copulas.CTreeParRef, tPaths)
@test tPar == tPar2

@test tPar == convert(Copulas.CTreeParRef, tPar)
@test tPaths == convert(Copulas.CTreePaths, tPaths)

#########################
## dimension functions ##
#########################

## CTreePaths
##------------

paths = Array{Int, 1}[[1, 4], [1, 2, 3], [5, 6], [5, 7, 8]]
tr1 = Copulas.CTreePaths(9, paths)
tr2 = convert(Copulas.CTreeParRef, tr1)

## conversion to matrix
pathsArr = [1 1 5 5; 4 2 6 7; 0 3 0 8]
@test Copulas._cTreePaths2arr(paths) == pathsArr
@test Copulas._arr2CTreePaths(pathsArr) == paths

allTrees = [tr1, tr2]

for tr in allTrees
   @test Copulas.maxDepth(tr) == 3
   @test Copulas.width(tr) == 4
   @test Copulas.getDepths(tr) == [3, 2, 2, 3]
end

#############
## display ##
#############

## CTreePaths
##------------

paths = Array{Int, 1}[[1, 4], [1, 2, 3], [5, 6], [5, 7, 8]]
tP = Copulas.CTreePaths(9, paths)

display(tP)
display([tP, tP])

## CTreeParRef
##------------

tPar = convert(Copulas.CTreeParRef, tP)
display(tPar)
display([tPar, tPar])

##############
## getindex ##
##############

paths = Array{Int, 1}[[1, 4], [1, 2, 3], [5, 6], [5, 7, 8]]
tr = Copulas.CTreePaths(9, paths)

@test tr[3, 2] == 6
@test tr[3] == [5, 6]
@test tr[[3, 4], [2, 3]] == [6, 8]

################
## tree nodes ##
################

paths = Array{Int, 1}[[1, 4], [1, 2, 3], [5, 6], [5, 7, 8]]
tr1 = Copulas.CTreePaths(9, paths)
tr2 = convert(Copulas.CTreeParRef, tr1)

allTrees = [tr1, tr2]

for tr in allTrees
    @test Copulas.allNodes(tr) == [1, 2, 3, 4, 5, 6, 7, 8, 9]
    @test Copulas.allPathNodes(tr) == [1, 2, 3, 4, 5, 6, 7, 8]
end

@test 9 == Copulas.getRoot(tr1)
@test 9 == Copulas.getRoot(tr2)

#######################
## conditioning sets ##
#######################

## condSetChk
##-----------

tr = Copulas.CTreePaths(3, Array{Int, 1}[[4, 5, 6],[2, 1]])
Copulas.viz(tr)

@test Copulas.condSetChk([4, 5], tr)
@test Copulas.condSetChk([1, 2], tr)
@test Copulas.condSetChk([2, 1], tr)

tr = Copulas.CTreePaths(3, Array{Int, 1}[[4, 5, 6],[2, 1]])
tr = convert(Copulas.CTreeParRef, tr)
@test !Copulas.condSetChk([2, 4], tr)

## findAndSortCondSet
##-------------------

tr = Copulas.CTreePaths(3, Array{Int, 1}[[4, 5, 6],[2, 1]])
tr = convert(Copulas.CTreeParRef, tr)

@test [2; 1] == Copulas.findAndSortCondSet(tr, [2, 1])
@test [4; 5] == Copulas.findAndSortCondSet(tr, [5, 4])
@test_throws Exception Copulas.findAndSortCondSet(tr, [4, 6])

## getPathToRoot
##--------------

tr = Copulas.CTreePaths(3, Array{Int, 1}[[4, 5, 6],[2, 1]])
path = Copulas.getPathToRoot(tr, 6)
@test path == [4, 5]

tr = Copulas.CTreePaths(3, Array{Int, 1}[[4, 6],[2, 1], [2, 5]])
path = Copulas.getPathToRoot(tr, 2)
@test path == []

tr = Copulas.CTreePaths(3, Array{Int, 1}[[2, 1, 4],[2, 1, 5], [2, 6]])
path = Copulas.getPathToRoot(tr, 1)
@test path == [2]

path = Copulas.getPathToRoot(tr, 5)
@test path == [2, 1]

#############
## attach! ##
#############

## attach to end of path
tr = Copulas.CTreePaths(1, [2 3], [4])
Copulas.attach!(tr, 5, [2,3])

expOut = Copulas.CTreePaths(1, [2 3 5], [4])
@test tr == expOut

## attach with empty condset
tr = Copulas.CTreePaths(1, [2 3], [4])
Copulas.attach!(tr, 5, Array(Int, 0))

expOut = Copulas.CTreePaths(1, [2 3], [4], [5])
@test tr == expOut

## attach with non-matching condset
tr = Copulas.CTreePaths(1, [2 3], [4])
Copulas.attach!(tr, 6, [7, 8])

expOut = Copulas.CTreePaths(1, [2 3], [4])
@test tr == expOut

## attach to middle of path
tr = Copulas.CTreePaths(1, [2 3], [4])
Copulas.attach!(tr, 5, [2]) # does not sort

expOut = Copulas.CTreePaths(1, [2 3], [4], [2 5])
@test tr == expOut

end


text:
# ---------------
# Global Commands
# ---------------

handle("julia.set-global-client") do req, data
  global global_client = req[1]
  Base.notify(connect_notify)
  global_client
end

handle("cwd") do _, path
  path == nothing && return
  path = isfile(path) ? dirname(path) : path
  cd(path)
end

function command(cmd, data = Dict())
  data[:cmd] = cmd
  raise(global_client, "editor.eval.julia.command", data)
end

function popup(header, body="", buttons = [@d(:label => "Ok")])
  command("popup",
          @d(:header => header,
             :body => body,
             :buttons => buttons))
end

function ltprint(message; error = false)
  command("print",
          @d(:value => message,
             :error => error))
end

function notify_done(msg = nothing)
  command("done", @d(:msg => msg))
end

function notify(message; class = "")
  command("notify",
          @d(:msg => message,
             :class => class))
end

function notify_error(message)
  notify(message, class = "error")
end

function console(value::AbstractString; html = false)
  command("console",
          @d("value" => value,
             "html"  => html))
end


text:
printstyled("  uniqueness of memory allocation\n", color=:light_green)
nx = 1024

S = SeisData(2)
for i=1:S.n
  # Check arrays of same type
  @test pointer_from_objref(S.loc[i]) != pointer_from_objref(S.x[i])
  @test pointer_from_objref(S.gain) != pointer_from_objref(S.fs)
end

for i in datafields
  v = getfield(S,i)
  for j in datafields
    if i != j
      @test pointer_from_objref(getfield(S,i)) != pointer_from_objref(getfield(S,j))
    end
  end
end

# Manually set all fields in the stupidest way possible
S.id[1] *= "...YYY"
S.name[1] *= "New Station"
S.loc[1] = UTMLoc()
S.fs[1] += 100.0
S.gain[1] *= 2.0
S.resp[1] = GenResp(resp=[0.0+0.0*im 1.0+1.0*im; 0.0+0.0*im 1.0-1.0*im])
S.units[1] *= "Unknown"
S.src[1] *= "www"
append!(S.notes[1], ["sadfasgasfg","kn4ntl42ntlk4n"])
S.misc[1]["P"] = 2.0
S.t[1] = [S.t[1]; [1 round(Int, time()*1.0e6); nx 0]]
append!(S.x[1], rand(nx))

for f in datafields
  v = getfield(S, f)
  @test v[1] != v[2]
end

# Now set all fields like a smart person
S = SeisData(2)
S.id[1] = "...YYY"
S.name[1] = "New Station"
S.loc[1] = GeoLoc()
S.fs[1] = 100.0
S.gain[1] = 2.0
S.resp[1] = PZResp([0.0+0.0*im 1.0+0.767*im; 0.0+0.0*im 1.0-0.767*im])
S.units[1] = "Unknown"
S.src[1] = "www"
S.notes[1] = ["sadfasgasfg","kn4ntl42ntlk4n"]
S.misc[1] = Dict{String,Any}("P" => 2.0)
S.t[1] = [1 round(Int, time()*1.0e6); nx 0]
S.x[1] = rand(nx)

S.id[2] = "...zzz"
S.name[2] = "Old Station"
S.loc[2] = GeoLoc()
S.fs[2] = 50.0
S.gain[2] = 22.0
S.resp[2] = PZResp64(z = [0.0+0.0*im, 0.0+0.0*im], p = [1.0+1.0*im, 1.0-1.0*im])
S.units[2] = "ms/2"
S.src[2] = "file"
S.notes[2] = ["0913840183","klnelgng"]
S.misc[2] = Dict{String,Any}("S" => 6.5)
S.t[2] = [1 round(Int, time()*1.0e6); nx 0]
S.x[2] = rand(nx)

for i in datafields
  v = getfield(S,i)
  if !isimmutable(v[1])
    @test pointer_from_objref(v[1]) != pointer_from_objref(v[2])
  end
end

# Unambiguous numbered initialization
S1 = SeisData(SeisChannel(), SeisChannel())
S2 = SeisData(2)
S3 = SeisData(S2)
S4 = SeisData(SeisData(1), SeisChannel())
S5 = SeisData(SeisChannel(), SeisData(1))

@test (S1 == S2 == S3 == S4 == S5)


text:
@testset "Checking outputs on several examples" begin
    test_cases = []

    #---------------------------------------

    ode = @ODEmodel(
        x1'(t) = a * x1(t) - b * x1(t) * x2(t),
        x2'(t) = -c * x2(t) + d * x1(t) * x2(t),
        y1(t) = x1(t)
    )
    correct = Dict(
        "globally" => Set([a, c, d, x1]),
        "locally_not_globally" => Set(),
        "nonidentifiable" => Set([b, x2])
    )
    push!(
        test_cases,
        (ode, correct)
    )

    #---------------------------------------

    ode = @ODEmodel(
        x1'(t) = 0,
        y1(t) = x1(t),
        y2(t) = a * x1(t) + a^2
    )
    correct = Dict(
        "globally" => Set([x1]),
        "locally_not_globally" => Set([a]),
        "nonidentifiable" => Set()
    )
    push!(
        test_cases,
        (ode, correct)
    )

    #---------------------------------------

    ode = @ODEmodel(
        S'(t) = -b * S(t) * In(t) / N(t),
        E'(t) = b * S(t) * In(t) / N(t) - nu * E(t),
        In'(t) = nu * E(t) - a * In(t),
        N'(t) = 0,
        y1(t) = In(t),
        y2(t) = N(t)
    )
    correct = Dict(
        "globally" => Set([b, In, N]),
        "locally_not_globally" => Set([nu, S, a, E]),
        "nonidentifiable" => Set()
    )
    push!(
        test_cases,
        (ode, correct)
    )

    #---------------------------------------

    ode = @ODEmodel(
        xA'(t) = -k1 * xA(t),
        xB'(t) = k1 * xA(t) - k2 * xB(t),
        xC'(t) = k2 * xB(t),
        eA'(t) = 0,
        eC'(t) = 0,
        y1(t) = xC(t),
        y2(t) = eA(t) * xA(t) + eB * xB(t) + eC(t) * xC(t),
        y3(t) = eA(t),
        y4(t) = eC(t)
    )
    correct = Dict(
        "globally" => Set([xC, eC, eA]),
        "locally_not_globally" => Set([xA, k1, xB, eB, k2]),
        "nonidentifiable" => Set()
    )
    push!(
        test_cases,
        (ode, correct)
    )

    #---------------------------------------

    for case in test_cases
        ode = case[1]
        result = identifiability_ode(ode, get_parameters(ode); nthrds = 1)
        @test case[2] == result
    end
end

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
using Arblib, ArbExtras, Random, SpecialFunctions, Test

include("demo_problems.jl")

@testset "ArbExtras" begin
    include("utilities.jl")

    include("isolate_roots.jl")
    include("refine_root.jl")

    include("extrema_polynomial.jl")
    include("extrema_series.jl")
    include("extrema_enclosure.jl")

    include("bounded_by.jl")

    include("special_functions.jl")

    include("integrate.jl")
end


text:
function validate_file(path, filetype = :text)
    isfile(path) || throw(ArgumentError("$(path) is not a valid file"))
    check_ext(path, filetype)
end

function check_ext(path, filetype = :text)
    _, ext = splitext(path)
    if filetype == :text
        return true
    elseif filetype == :markdown
        return lowercase(ext) in markdown_exts || throw(ArgumentError("$(path) is not a valid markdown file"))
    else
        throw(ArgumentError("Unrecognized filetype $(filetype)"))
    end
end

### common utils for DemoPage and DemoSection

function validate_order(order::AbstractArray, x::Union{DemoPage, DemoSection})
    default_order = get_default_order(x)
    if intersect(order, default_order) == union(order, default_order)
        return true
    else
        config_filepath = joinpath(x.root, config_filename)

        entries = join(string.(setdiff(order, default_order)), "\n")
        isempty(entries) || @warn("The following entries in $(config_filepath) are not used anymore:\n$(entries)")

        entries = join(string.(setdiff(default_order, order)), "\n")
        isempty(entries) || @warn("The following entries in $(config_filepath) are missing:\n$(entries)")

        throw(ArgumentError("incorrect order in $(config_filepath), please check the previous warning message."))
    end
end

"""
    walkpage([f=identity], page; flatten=true, max_depth=Inf)

Walk through the page structure and apply function f to each of the item.

By default, the page structure is not preserved in the result. To preserve that, you could pass
`flatten=false` to the function.

# Examples

This is particulaly useful to generate information for the page structure. For example:

```julia
julia> page = DemoPage("docs/quickstart/");

julia> walkpage(page; flatten=true) do item
    basename(item.path)
end

"Quick Start" => [
    "simple_markdown_demo.md",
    "configure_card.md",
    "configure_sec_and_page.md",
    "hide_your_card_from_index.md",
    "hidden_card.jl",
    "1.julia_demo.jl",
    "2.cover_on_the_fly.jl"
]
```

If `flatten=false`, then it gives you something like this:

```
"Quick Start" => [
    "Usage example" => [
        "Basics" => [
            "simple_markdown_demo.md",
            "configure_card.md",
            "configure_sec_and_page.md",
            "hide_your_card_from_index.md",
            "hidden_card.jl"
        ],
        "Julia demos" => [
            "1.julia_demo.jl",
            "2.cover_on_the_fly.jl"
        ]
    ]
]
```
"""
walkpage(page::DemoPage; kwargs...) = walkpage(identity, page; kwargs...)
function walkpage(f, page::DemoPage; flatten=true, kwargs...)
    nodes = _walkpage.(f, page.sections, 1; flatten=flatten, kwargs...)
    if flatten
        nodes = reduce(vcat, nodes)
        # compensate back to a vector when the first level section has only one section and one card
        nodes isa AbstractArray || (nodes = [nodes])
    end
    page.title => nodes
end

function _walkpage(f, sec::DemoSection, depth; max_depth=Inf, flatten=true)
    depth+1 > max_depth && return flatten ? f(sec) : [f(sec), ]
    if !isempty(sec.cards)
        return flatten ? mapreduce(f, vcat, sec.cards) : sec.title => f.(sec.cards)
    else
        map_fun(section) =  _walkpage(f, section, depth+1; max_depth=max_depth, flatten=flatten)
        return flatten ? mapreduce(map_fun, vcat, sec.subsections) : sec.title => map_fun.(sec.subsections)
    end
end

### regexes and configuration parsers

"""regex that capture the URL of Edit On GitHub button from Documenter generated html files"""
const regex_edit_on_github = r"\<a class=\"docs-edit-link\" href=\"(.*)\"\s*title=..*Edit on GitHub\</span\>\</a\>"

raw"""
    get_regex(card::AbstractDemoCard, regex_type)

Return regex used to parse `card`.

`regex_type` and capture fields:

* `:image`: $1=>title, $2=>path
* `:title`: $1=>title, (Optionally, $2=>id)
* `:content`: $1=>content
"""
function get_regex(::Val{:Markdown}, regex_type)
    if regex_type == :image
        # Example: ![title](path)
        return r"^\s*!\[(?<title>[^\]]*)\]\((?<path>[^\s]*)\)"
    elseif regex_type == :title
        # Example: # [title](@id id)
        regex_title = r"\s*#+\s*\[(?<title>[^\]]+)\]\(\@id\s+(?<id>[^\s\)]+)\)"
        # Example: # title
        regex_simple_title = r"\s*#+\s*(?<title>[^#]+)"

        # Note: return the complete one first
        return (regex_title, regex_simple_title)
    elseif regex_type == :content
        # lines that are not title, image, link, list
        # FIXME: list is also captured by this regex
        return r"^\s*(?<content>[^#\-*!].*)"
    else
        error("Unrecognized regex type: $(regex_type)")
    end
end

function get_regex(::Val{:Julia}, regex_type)
    if regex_type == :image
        # Example: #md ![title](path)
        return r"^#!?\w*\s+!\[(?<title>[^\]]*)\]\((?<path>[^\s]*)\)"
    elseif regex_type == :title
        # Example: # [title](@id id)
        regex_title = r"\s*#!?\w*\s+#+\s*\[(?<title>[^\]]+)\]\(\@id\s+(?<id>[^\s\)]+)\)"
        # Example: # title
        regex_simple_title = r"\s*#!?\w*\s+#+\s*(?<title>[^#]+)"

        # Note: return the complete one first
        return (regex_title, regex_simple_title)
    elseif regex_type == :content
        # lines that are not title, image, link, list
        # FIXME: list is also captured by this regex
        return r"^\s*#!?\w*\s+(?<content>[^#\-\*!]+.*)"
    else
        error("Unrecognized regex type: $(regex_type)")
    end
end

# YAML frontmatter
# 1. markdown: ---
# 2. julia: # ---
const regex_yaml = r"^#?\s*---"

# markdown URL: [text](url)
const regex_md_url = r"\[(?<text>[^\]]*)\]\((?<url>[^\)]*)\)"

"""
    parse([T::Val], card::AbstractDemoCard)
    parse(T::Val, contents)
    parse(T::Val, path)

Parse the content of `card` and return the configuration.

Currently supported items are: `title`, `id`, `cover`, `description`.

!!! note

    Users of this function need to use `haskey` to check if keys are existed.
    They also need to validate the values.
"""
function parse(T::Val, card::AbstractDemoCard)
    header, frontmatter, body = split_frontmatter(readlines(card.path))
    config = parse(T, body)
    # frontmatter has higher priority
    if !isempty(frontmatter)
        yaml_config = try
            YAML.load(join(frontmatter, "\n"))
        catch err
            @warn "failed to parse YAML frontmatter, please double check the format"
            println(stderr, frontmatter)
            rethrow(err)
        end
        merge!(config, yaml_config)
    end

    return config
end
parse(card::JuliaDemoCard) = parse(Val(:Julia), card)
parse(card::MarkdownDemoCard) = parse(Val(:Markdown), card)

function parse(T::Val, contents::String)
    # if we just check isfile then it's likely to complain that filename is too long
    # compat: for Julia > 1.1 `endswith` supports regex as well
    if match(r"\.(\w+)$", contents) isa RegexMatch && isfile(contents)
        contents = readlines(contents)
    else
        contents = split(contents, "\n")
    end
    parse(T, contents)
end
function parse(T::Val, contents::AbstractArray{<:AbstractString})::Dict
    config = Dict()
    # only parse relevant information from the main contents body
    _, _, contents = split_frontmatter(contents)

    # it's too complicated to regex title correctly from body contents, so a simple
    # strategy is to limit possible titles to lines before the contents
    body_matches = map(contents) do line
        re = get_regex(T, :content)
        m = match(re, line)
    end
    body_lines = findall(map(x->x isa RegexMatch, body_matches))
    start_of_body = isempty(body_lines) ? 1 : minimum(body_lines)

    # The first title line in markdown format is parse out as the title
    title_matches = map(contents[1:start_of_body-1]) do line
        # try to match the most complete pattern first
        for re in get_regex(T, :title)
            m = match(re, line)
            m isa RegexMatch && return m
        end
        return nothing
    end
    title_lines = findall(map(x->x isa RegexMatch, title_matches))

    if !isempty(title_lines)
        m = title_matches[title_lines[1]]
        config["title"] = m["title"]
        if length(m.captures) == 2
            config["id"] = m["id"]
        end
    end

    # The first valid image link is parsed out as card cover
    image_matches = map(contents) do line
        m = match(get_regex(T, :image), line)
        m isa RegexMatch && return m
        return nothing
    end
    image_lines = findall(map(x->x isa RegexMatch, image_matches))
    if !isempty(image_lines)
        config["cover"] = image_matches[image_lines[1]]["path"]
    end

    description = parse_description(contents, get_regex(T, :content))
    if !isnothing(description)
        config["description"] = description
    end

    return config
end

function parse_description(contents::AbstractArray{<:AbstractString}, regex)
    # description as the first paragraph that is not a title, image, list or codes
    content_lines = map(x->match(regex, x) isa RegexMatch, contents)
    code_lines = map(contents) do line
        match(r"#?!?\w*```", lstrip(line)) isa RegexMatch
    end |> findall
    for (i,j) in zip(code_lines[1:2:end], code_lines[2:2:end])
        # mark code lines as non-content lines
        content_lines[i:j] .= 0
    end
    m = findall(content_lines)
    isempty(m) && return nothing

    paragraph_line = findall(x->x!=1, diff(m))
    offset = isempty(paragraph_line) ? length(m) : paragraph_line[1]
    description = map(contents[m[1:offset]]) do line
        m = match(regex, line)
        m["content"]
    end
    description = join(description, " ")
    return replace(description, regex_md_url => s"\g<text>")
end

"""
    split_frontmatter(contents) -> frontmatter, body

splits the YAML frontmatter out from markdown and julia source code. Leading `# ` will be
stripped for julia codes.

`contents` can be `String` or vector of `String`. Outputs have the same type of `contents`.
"""
function split_frontmatter(contents::String)
    header, frontmatter, body = split_frontmatter(split(contents, "\n"))
    return join(header, "\n"), join(frontmatter, "\n"), join(body, "\n")
end
function split_frontmatter(contents::AbstractArray{<:AbstractString})
    offsets = map(contents) do line
        m = match(regex_yaml, line)
        m isa RegexMatch
    end
    offsets = findall(offsets)

    length(offsets) < 2 && return "", "", contents
    if offsets[1] != 1 && !all(x->isempty(strip(x)) || startswith(strip(x), "#"), contents[1:offsets[1]-1])
        # For julia demo, comments and empty lines are allowed before frontmatter
        # For markdown demo, only empty lines are allowed before frontmatter
        return "", "", contents
    end

    # infer how many spaces we need to trim from the start line. E.g.,  "# ---" => 1
    start_line = contents[offsets[1]]
    m = match(r"^#?(\s*)", start_line)
    # make sure we strip the same amount of whitespaces -- preserve indentation
    indent_spaces = isnothing(m) ? "" : m.captures[1]
    frontmatter = map(contents[offsets[1]: offsets[2]]) do line
        m = match(Regex("^#?$(indent_spaces)(.*)"), line)
        if isnothing(m)
            @warn "probably incorrect YAML syntax or indentation" line
            # we don't know much about what `line` is, so do nothing here and let YAML complain about it
            line
        else
            m.captures[1]
        end
    end

    return contents[1:offsets[1]-1], frontmatter, contents[offsets[2]+1:end]
end

function get_default_title(x::Union{AbstractDemoCard, DemoSection, DemoPage})
    name_without_ext = splitext(basename(x))[1]
    strip(replace(uppercasefirst(name_without_ext), r"[_-]" => " "))
end

is_remote_url(path) = !isnothing(match(r"^https?://", path))

"""
    input_bool(prompt)::Bool

Print `prompt` message and trigger user input for confirmation.
"""
function input_bool(prompt)
    while true
        println(prompt, " [y/n]")
        response = readline()
        length(response) == 0 && continue
        reply = lowercase(first(strip(response)))
        if reply == 'y'
            return true
        elseif reply =='n'
            return false
        end
        # Otherwise loop and repeat the prompt
    end
end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function test_intersections(fname)
    Drawing(2000, 2000, fname)
    origin()
    sethue("black")
    setline(2)
    setopacity(0.95)

    tiles = Tiler(2000, 2000, 5, 5)

    for (pos, n) in tiles
        gsave()
        translate(pos)
        box(O, tiles.tilewidth, tiles.tileheight, :clip)
        p1 = Point(rand(-tiles.tilewidth/2:tiles.tilewidth/2), rand(-tiles.tileheight/2:tiles.tilewidth/2))
        p2 = Point(rand(-tiles.tilewidth/2:tiles.tilewidth/2), rand(-tiles.tileheight/2:tiles.tilewidth/2))
        p3 = Point(rand(-tiles.tilewidth/2:tiles.tilewidth/2), rand(-tiles.tileheight/2:tiles.tilewidth/2))
        map(pt -> circle(pt, 2, :fill), [p1, p2, p3])

        gsave()
        randomhue()
        rad = rand(20:tiles.tilewidth/3)
        circle(p2, rad, :fill)
        sethue("black")
        line(p1, p3, :stroke)

        sethue("green")
        setdash("dot")
        setline(0.5)
        line(between(p1, p3, -5), between(p1, p3, 5), :stroke)
        grestore()

        sethue("red")
        n, int1, int2 = intersectionlinecircle(p1, p3, p2, rad)
        if n == 1
            circle(int1, 10, :fill)
        elseif n == 2
            circle(int1, 10, :fill)
            circle(int2, 10, :fill)
        end
        sethue("black")
        text("$n intersection(s)")
        box(O, tiles.tilewidth, tiles.tileheight, :stroke)
        grestore()
    end
    l1 = Luxor.Point(-100.0, -100.0)
    l2 = Luxor.Point(300.0, 200.0)
    cpoint = Point(0, 0)
    @test intersectionlinecircle(l1, l2, cpoint, 15)[1] == 0
    @test intersectionlinecircle(l1, l2, cpoint, 20)[1] == 1
    @test intersectionlinecircle(l1, l2, cpoint, 25)[1] == 2
    @test finish() == true
end

test_intersections("intersection_line_circle_test.pdf")
println("...finished test: output in $(fname)")


text:
using Pkg
pwd()
cd("../ChemometricsTools.jl/")
Pkg.activate(".")

using ChemometricsTools

#View the data in the package space
ChemometricsToolsDatasets()
#Load the data from package-space
DF = ChemometricsToolsDataset("tecator.csv");

Xindices = findall(occursin.( "absorbance", string.(names(DF) ) ) );
Tecator = convert(Matrix,  DF[:,Xindices]);

#Let's check the Tecator data for duplicate entries ;)
length(AssessHealth(Tecator)["Duplicates"])
#Eeek plenty of duplicates, oh well this is a shoot out not a real analysis...

using Plots
plot(Tecator', legend = false, ylabel = "Absorbance", xlabel = "bins")
using Statistics
BC = boxcar(Tecator[1,:]; windowsize = 5, fn = mean);
AS = ALSSmoother(Tecator[1,:]; lambda = 100, p = 0.1, maxiters = 10);
PS = PerfectSmoother(Tecator[1,:]; lambda = 100);
SG = SavitzkyGolay(convert.(Float64, Tecator[1,:])', 0, 1, 3);

plot( [ Tecator[1,:], BC', AS', PS', SG' ], xlabel = "bins", ylabel = "Smoothed Signal",
        label = [ "Original", "BC", "AS", "PS", "SG" ], ylim = (2.5,3.75) )

FD = FirstDerivative(Tecator[1,:]);
SD = SecondDerivative(Tecator[1,:]);
HD = FractionalDerivative(convert.(Float64, Tecator[1,:]); Order = 1.5);

plot( [ FD', SD', HD' ], xlabel = "bins", ylabel = "Differentiated Signal",
        label = [ "First Derivative", "Second Derivative", "1(1/2) Derivative" ],
        ylim = (-0.1,0.1) )


text:
# Replicate the Biogeme Swissmetro example
# https://biogeme.epfl.ch/examples/swissmetro/01logit.html

@testset "Biogeme swissmetro" begin
    using CSV
    using DataFrames
    using DiscreteChoiceModels
    using Test
    using StatsBase

    data = CSV.File(joinpath(dirname(Base.source_path()), "../data/biogeme_swissmetro.dat"), delim='\t') |> DataFrame
    data = data[in.(data.PURPOSE, [Set([1, 3])]) .& (data.CHOICE .!= 0), :]

    @test nrow(data) == 6768

    data.avtr = (data.TRAIN_AV .== 1) .& (data.SP .!= 0)
    data.avsm = data.SM_AV .== 1
    data.avcar = (data.CAR_AV .== 1) .& (data.SP .!= 0)

    model = multinomial_logit(
        @utility(begin
            1 ~ αtrain + βtravel_time * TRAIN_TT / 100 + βcost * (TRAIN_CO * (GA == 0)) / 100
            2 ~ αswissmetro + βtravel_time * SM_TT / 100 + βcost * SM_CO * (GA == 0) / 100
            3 ~ αcar + βtravel_time * CAR_TT / 100 + βcost * CAR_CO / 100

            αswissmetro = 0, fixed  # fix swissmetro ASC to zero
        end),
        :CHOICE,
        data,
        availability=[
            1 => :avtr,
            2 => :avsm,
            3 => :avcar,
        ]
    )

    coefs = Dict(zip(coefnames(model), round.(coef(model), sigdigits=3)))
    ses = Dict(zip(coefnames(model), round.(stderror(model), sigdigits=3)))

    @test coefs[:αcar] ≈ -0.155 # biogeme says 0.155 which would round to 0.16, but it's actually 0.1546
    @test coefs[:αswissmetro] == 0  # fixed
    @test coefs[:αtrain] ≈ -0.701
    @test coefs[:βcost] ≈ -1.08
    @test coefs[:βtravel_time] ≈ -1.28

    @test ses[:αcar] ≈ 0.0432
    @test isnan(ses[:αswissmetro])  # fixed
    @test ses[:αtrain] ≈ 0.0549
    @test ses[:βcost] ≈ 0.0518
    @test ses[:βtravel_time] ≈ 0.0569

    idx = Dict(zip(coefnames(model), 1:length(coefs)))
    vcovmat = round.(vcov(model), sigdigits=3)
    @test vcovmat[idx[:αtrain], idx[:αcar]] ≈ 0.00138
    @test vcovmat[idx[:βcost], idx[:αcar]] ≈ 0.000485
    @test vcovmat[idx[:βcost], idx[:αtrain]] ≈ 8.22e-6
    @test vcovmat[idx[:βtravel_time], idx[:αcar]] ≈ -0.00144
    @test vcovmat[idx[:βtravel_time], idx[:αtrain]] ≈ -0.00225
    @test vcovmat[idx[:βtravel_time], idx[:βcost]] ≈ 0.000550

    @test all(isnan.(vcovmat[idx[:αswissmetro], :]))
    @test all(isnan.(vcovmat[:, idx[:αswissmetro]]))

    @test round(loglikelihood(model), digits=3) ≈ -5331.252
    @test round(model.init_ll, digits=3) ≈ -6964.663
end


text:
type BernoulliArm <: BanditArm
  p::Float64
end

function draw(arm::BernoulliArm)
  if rand() > arm.p
    0
  else
    1
  end
end

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


text:
using TensorKit
using Zygote
using LinearAlgebra
using Optim
using JLD2
using Random

using Revise
using statmech_tm_solver

# MPO for the triangular AF Ising
T = mpo_triangular_AF_ising()
T_adj = reshape(T.data, (2, 2, 2, 2))
T_adj = permutedims(conj.(T_adj), (1, 3, 2, 4))
T_adj = TensorMap(T_adj, ℂ^2*ℂ^2, ℂ^2*ℂ^2)

function f(arr::Array{ComplexF64, 3})
    phi = arr_to_TensorMap(arr)
    #phi = act(Tn, phi)
    return nonherm_cost_func(T, phi)
end

function fT(arr::Array{ComplexF64, 3})
    phi = arr_to_TensorMap(arr)
    return nonherm_cost_func(T_adj, phi)
end

function g!(grad::Array{ComplexF64, 3}, psidata::Array{ComplexF64, 3})
    grad .= gradient(f, psidata)[1]
end

function gT!(grad::Array{ComplexF64, 3}, psidata::Array{ComplexF64, 3})
    grad .= gradient(fT, psidata)[1]
end

psidata = rand(MersenneTwister(1), ComplexF64, (2, 2, 2))

#io = open("result_variational.txt", "w")
#close(io)
for chi in [2, 4, 8, 16, 32, 64, 128, 256]
    global psidata
    io = open("result_variational.txt", "a+")

    println(size(psidata))

    res_f = optimize(f, g!, psidata, LBFGS(), Optim.Options(show_trace=true, iterations=200))
    res_fT = optimize(fT, gT!, psidata, LBFGS(), Optim.Options(show_trace=true, iterations=200))

    psidata_final = Optim.minimizer(res_f)
    psi = arr_to_TensorMap(psidata_final)

    psiLdata_final = Optim.minimizer(res_fT)
    psiL = arr_to_TensorMap(psiLdata_final)

    cost_func_final = nonherm_cost_func(T, psi)
    F_final_1 = free_energy(T, psi)
    F_final = ovlp(psiL, act(T, psi)) / ovlp(psiL, psi) |> norm |> log
    println(io, chi, ' ', F_final, ' ', F_final_1, ' ', cost_func_final )
    quicksave("ckpt_variational_chi$(chi)", psi)
    quicksave("ckpt_variational_left_chi$(chi)", psiL)

    psidata = toarray(act(T, psi))

    close(io)
end


text:
# Tests for Sinusoidal arms

using MAB
import PyPlot

# Create different Sinusoidal arms
a1  = Arms.Sinusoidal( 10 )
a2  = Arms.Sinusoidal( 8 )
a3  = Arms.Sinusoidal( 12, π )
a4  = Arms.Sinusoidal( 14, π/3 )

# Plot the rewards
noOfTimeSteps   = 40
fig     = PyPlot.figure()
PyPlot.plot( 1:noOfTimeSteps, [ Arms.pull!(a1) for i=1:noOfTimeSteps ], label = "a1" )
PyPlot.plot( 1:noOfTimeSteps, [ Arms.pull!(a2) for i=1:noOfTimeSteps ], label = "a2" )
PyPlot.plot( 1:noOfTimeSteps, [ Arms.pull!(a3) for i=1:noOfTimeSteps ], label = "a3" )
PyPlot.plot( 1:noOfTimeSteps, [ Arms.pull!(a4) for i=1:noOfTimeSteps ], label = "a4" )
ax = PyPlot.gca()
PyPlot.grid()
ax[:set_ylim]( [-1.20,+1.20] )
PyPlot.legend()


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iot_data_plane
using AWS.Compat
using AWS.UUIDs

"""
    delete_thing_shadow(thing_name)
    delete_thing_shadow(thing_name, params::Dict{String,<:Any})

Deletes the shadow for the specified thing. Requires permission to access the
DeleteThingShadow action. For more information, see DeleteThingShadow in the IoT Developer
Guide.

# Arguments
- `thing_name`: The name of the thing.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"name"`: The name of the shadow.
"""
function delete_thing_shadow(thingName; aws_config::AbstractAWSConfig=global_aws_config())
    return iot_data_plane(
        "DELETE",
        "/things/$(thingName)/shadow";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_thing_shadow(
    thingName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iot_data_plane(
        "DELETE",
        "/things/$(thingName)/shadow",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_retained_message(topic)
    get_retained_message(topic, params::Dict{String,<:Any})

Gets the details of a single retained message for the specified topic. This action returns
the message payload of the retained message, which can incur messaging costs. To list only
the topic names of the retained messages, call ListRetainedMessages. Requires permission to
access the GetRetainedMessage action. For more information about messaging costs, see IoT
Core pricing - Messaging.

# Arguments
- `topic`: The topic name of the retained message to retrieve.

"""
function get_retained_message(topic; aws_config::AbstractAWSConfig=global_aws_config())
    return iot_data_plane(
        "GET",
        "/retainedMessage/$(topic)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_retained_message(
    topic, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iot_data_plane(
        "GET",
        "/retainedMessage/$(topic)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_thing_shadow(thing_name)
    get_thing_shadow(thing_name, params::Dict{String,<:Any})

Gets the shadow for the specified thing. Requires permission to access the GetThingShadow
action. For more information, see GetThingShadow in the IoT Developer Guide.

# Arguments
- `thing_name`: The name of the thing.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"name"`: The name of the shadow.
"""
function get_thing_shadow(thingName; aws_config::AbstractAWSConfig=global_aws_config())
    return iot_data_plane(
        "GET",
        "/things/$(thingName)/shadow";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_thing_shadow(
    thingName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iot_data_plane(
        "GET",
        "/things/$(thingName)/shadow",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_named_shadows_for_thing(thing_name)
    list_named_shadows_for_thing(thing_name, params::Dict{String,<:Any})

Lists the shadows for the specified thing. Requires permission to access the
ListNamedShadowsForThing action.

# Arguments
- `thing_name`: The name of the thing.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"nextToken"`: The token to retrieve the next set of results.
- `"pageSize"`: The result page size.
"""
function list_named_shadows_for_thing(
    thingName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iot_data_plane(
        "GET",
        "/api/things/shadow/ListNamedShadowsForThing/$(thingName)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_named_shadows_for_thing(
    thingName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iot_data_plane(
        "GET",
        "/api/things/shadow/ListNamedShadowsForThing/$(thingName)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_retained_messages()
    list_retained_messages(params::Dict{String,<:Any})

Lists summary information about the retained messages stored for the account. This action
returns only the topic names of the retained messages. It doesn't return any message
payloads. Although this action doesn't return a message payload, it can still incur
messaging costs. To get the message payload of a retained message, call GetRetainedMessage
with the topic name of the retained message. Requires permission to access the
ListRetainedMessages action. For more information about messaging costs, see IoT Core
pricing - Messaging.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return at one time.
- `"nextToken"`: To retrieve the next set of results, the nextToken value from a previous
  response; otherwise null to receive the first set of results.
"""
function list_retained_messages(; aws_config::AbstractAWSConfig=global_aws_config())
    return iot_data_plane(
        "GET", "/retainedMessage"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_retained_messages(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iot_data_plane(
        "GET",
        "/retainedMessage",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    publish(topic)
    publish(topic, params::Dict{String,<:Any})

Publishes an MQTT message. Requires permission to access the Publish action. For more
information about MQTT messages, see MQTT Protocol in the IoT Developer Guide. For more
information about messaging costs, see IoT Core pricing - Messaging.

# Arguments
- `topic`: The name of the MQTT topic.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"payload"`: The message body. MQTT accepts text, binary, and empty (null) message
  payloads. Publishing an empty (null) payload with retain = true deletes the retained
  message identified by topic from IoT Core.
- `"qos"`: The Quality of Service (QoS) level.
- `"retain"`: A Boolean value that determines whether to set the RETAIN flag when the
  message is published. Setting the RETAIN flag causes the message to be retained and sent to
  new subscribers to the topic. Valid values: true | false  Default value: false
"""
function publish(topic; aws_config::AbstractAWSConfig=global_aws_config())
    return iot_data_plane(
        "POST", "/topics/$(topic)"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function publish(
    topic, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iot_data_plane(
        "POST",
        "/topics/$(topic)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_thing_shadow(payload, thing_name)
    update_thing_shadow(payload, thing_name, params::Dict{String,<:Any})

Updates the shadow for the specified thing. Requires permission to access the
UpdateThingShadow action. For more information, see UpdateThingShadow in the IoT Developer
Guide.

# Arguments
- `payload`: The state information, in JSON format.
- `thing_name`: The name of the thing.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"name"`: The name of the shadow.
"""
function update_thing_shadow(
    payload, thingName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iot_data_plane(
        "POST",
        "/things/$(thingName)/shadow",
        Dict{String,Any}("payload" => payload);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_thing_shadow(
    payload,
    thingName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iot_data_plane(
        "POST",
        "/things/$(thingName)/shadow",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("payload" => payload), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
module ReactiveMPVariableTest

using Test
using ReactiveMP
using Rocket

@testset "Variable" begin

    @testset "setmarginal!" begin

        begin
            r = randomvar(:r)
            setmarginal!(r, NormalMeanVariance(-2.0, 3.0))
            subscribe!(getmarginal(r, IncludeAll()), (marginal) -> begin
                @test typeof(marginal) <: Marginal{ <: NormalMeanVariance }
                @test mean(marginal) === -2.0
                @test var(marginal) === 3.0
            end)
        end

        begin
            r = randomvar(:r)
            setmarginal!(r, NormalMeanPrecision(-2.0, 3.0))
            subscribe!(getmarginal(r, IncludeAll()), (marginal) -> begin
                @test typeof(marginal) <: Marginal{ <: NormalMeanPrecision }
                @test mean(marginal) === -2.0
                @test precision(marginal) === 3.0
            end)
        end

        begin
            rs = randomvar(:r, 2)
            setmarginals!(rs, NormalMeanVariance(-2.0, 3.0))
            subscribe!(getmarginals(rs, IncludeAll()), (marginals) -> begin
                foreach(marginals) do marginal
                    @test typeof(marginal) <: Marginal{ <: NormalMeanVariance }
                    @test mean(marginal) === -2.0
                    @test var(marginal) === 3.0
                end
            end)
        end

        begin
            rs = randomvar(:r, 2)
            setmarginals!(rs, NormalMeanPrecision(-2.0, 3.0))
            subscribe!(getmarginals(rs, IncludeAll()), (marginals) -> begin
                foreach(marginals) do marginal
                    @test typeof(marginal) <: Marginal{ <: NormalMeanPrecision }
                    @test mean(marginal) === -2.0
                    @test precision(marginal) === 3.0
                end
            end)
        end

        # getmarginals is broken for matrices
        # begin
        #     rs = randomvar(:r, 2, 2)
        #     setmarginals!(rs, NormalMeanVariance(-2.0, 3.0))
        #     subscribe!(getmarginals(rs, IncludeAll()), (marginals) -> begin
        #         foreach(marginals) do marginal
        #             @test typeof(marginal) <: Marginal{ <: NormalMeanPrecision }
        #             @test mean(marginal) === -2.0
        #             @test variance(marginal) === 3.0
        #         end
        #     end)
        # end

        # begin
        #     rs = randomvar(:r, 2, 2)
        #     setmarginals!(rs, NormalMeanPrecision(-2.0, 3.0))
        #     subscribe!(getmarginals(rs, IncludeAll()), (marginals) -> begin
        #         foreach(marginals) do marginal
        #             @test typeof(marginal) <: Marginal{ <: NormalMeanPrecision }
        #             @test mean(marginal) === -2.0
        #             @test precision(marginal) === 3.0
        #         end
        #     end)
        # end

    end

end

end


text:
# Function to relabel a lattice to bipartite labeling
function relabelSitesBipartite!(
            lattice :: L,
            label_A :: LS,
            label_B :: LS,
            seed_site_A :: Int64 = 1
        ) where {LS,D,S<:AbstractSite{LS,D},B,U,L<:AbstractLattice{S,B,U}}

    #-------------
    # PREPARATION
    #-------------

    # indicate that all labels are not set yet
    relabeled  = Bool[false for i in 1:numSites(lattice)]
    # save the sublattice as an integer
    sublattice = Int64[-1 for i in 1:numSites(lattice)]

    # get the ordered connection list of the lattice
    bonds_organized = organizedBondsFrom(lattice)

    # set the first label to be sublattice A
    label!(site(lattice, 1), deepcopy(label_A))
    # set the sublattice to be 1
    sublattice[1] = 1
    # indicate that this site has indeed been relabeled
    relabeled[1] = true

    #------------
    # RELABELING
    #------------

    # iterate while there are unlabeled labels left
    while false in relabeled
        # look through all sites
        for i in 1:numSites(lattice)
            # continue if label already set
            if relabeled[i]
                continue
            end
            # check the neighbors for their sublattice
            neighbor_sublattice = -1
            for b in bonds_organized[i]
                if relabeled[to(b)]
                    # check if the sublattice has to be set
                    if neighbor_sublattice < 0
                        neighbor_sublattice = sublattice[to(b)]
                    elseif neighbor_sublattice != sublattice[to(b)]
                        # ERROR - lattice is not bipartite
                        error("The given lattice cannot be labeled to bipartite")
                    end
                end
            end
            # check if neighbors sublattice could be identified
            if neighbor_sublattice > 0
                # set the label depending on the neighbor sublattice
                if neighbor_sublattice == 1
                    # set the sublattice to 2
                    sublattice[i] = 2
                    # set the label B
                    label!(site(lattice, i), deepcopy(label_B))
                    # indicate that the label has been found
                    relabeled[i] = true
                    # set all neighbor sublattices
                    for b in bonds_organized[i]
                        # set the sublattice to 1
                        sublattice[to(b)] = 1
                        # set the label A
                        label!(site(lattice, to(b)), deepcopy(label_A))
                        # indicate that the label has been found
                        relabeled[to(b)] = true
                    end
                elseif neighbor_sublattice == 2
                    # set the sublattice to 1
                    sublattice[i] = 1
                    # set the label A
                    label!(site(lattice, i), deepcopy(label_A))
                    # indicate that the label has been found
                    relabeled[i] = true
                    # set all neighbor sublattices
                    for b in bonds_organized[i]
                        # set the sublattice to 2
                        sublattice[to(b)] = 2
                        # set the label B
                        label!(site(lattice, to(b)), deepcopy(label_B))
                        # indicate that the label has been found
                        relabeled[to(b)] = true
                    end
                else
                    # give an error
                    error("somewhere the neighbor sublattice took the value "*string(neighbor_sublattice))
                end
            end
        end
    end
end

# wrapper using the default labels A and B
function relabelSitesBipartite!(
        lattice     :: L,
        seed_site_A :: Int64 = 1
    ) where {LS,D,S<:AbstractSite{LS,D},B,U,L<:AbstractLattice{S,B,U}}

    # call the general function
    return relabelSitesBipartite!(lattice, getDefaultLabelA(LS), getDefaultLabelB(LS), seed_site_A)
end

# Function to relabel a unitcell to bipartite labeling
function relabelSitesBipartite!(
            uc      :: U,
            label_A :: LS,
            label_B :: LS,
            seed_site_A :: Int64 = 1
        ) where {LS,D,S<:AbstractSite{LS,D},B,U<:AbstractUnitcell{S,B}}

    #-------------
    # PREPARATION
    #-------------

    # indicate that all labels are not set yet
    relabeled  = Bool[false for i in 1:numSites(uc)]
    # save the sublattice as an integer
    sublattice = Int64[-1 for i in 1:numSites(uc)]

    # get the ordered connection list of the lattice
    bonds_organized = organizedBondsFrom(uc)

    # set the first label to be sublattice A
    label!(site(uc, 1), deepcopy(label_A))
    # set the sublattice to be 1
    sublattice[1] = 1
    # indicate that this site has indeed been relabeled
    relabeled[1] = true

    #------------
    # RELABELING
    #------------

    # iterate while there are unlabeled labels left
    while false in relabeled
        # look through all sites
        for i in 1:numSites(uc)
            # continue if label already set
            if relabeled[i]
                continue
            end
            # check the neighbors for their sublattice
            neighbor_sublattice = -1
            for b in bonds_organized[i]
                if relabeled[to(b)]
                    # check if the sublattice has to be set
                    if neighbor_sublattice < 0
                        neighbor_sublattice = sublattice[to(b)]
                    elseif neighbor_sublattice != sublattice[to(b)]
                        # ERROR - unitcell is not bipartite
                        error("The given unitcell cannot be labeled to bipartite")
                    end
                end
            end
            # check if neighbors sublattice could be identified
            if neighbor_sublattice > 0
                # set the label depending on the neighbor sublattice
                if neighbor_sublattice == 1
                    # set the sublattice to 2
                    sublattice[i] = 2
                    # set the label B
                    label!(site(uc, i), deepcopy(label_B))
                    # indicate that the label has been found
                    relabeled[i] = true
                    # set all neighbor sublattices
                    for b in bonds_organized[i]
                        # set the sublattice to 1
                        sublattice[to(b)] = 1
                        # set the label A
                        label!(site(uc, to(b)), deepcopy(label_A))
                        # indicate that the label has been found
                        relabeled[to(b)] = true
                    end
                elseif neighbor_sublattice == 2
                    # set the sublattice to 1
                    sublattice[i] = 1
                    # set the label A
                    label!(site(uc, i), deepcopy(label_A))
                    # indicate that the label has been found
                    relabeled[i] = true
                    # set all neighbor sublattices
                    for b in bonds_organized[i]
                        # set the sublattice to 2
                        sublattice[to(b)] = 2
                        # set the label B
                        label!(site(uc, to(b)), deepcopy(label_B))
                        # indicate that the label has been found
                        relabeled[to(b)] = true
                    end
                else
                    # give an error
                    error("somewhere the neighbor sublattice took the value "*string(neighbor_sublattice))
                end
            end
        end
    end
end

# wrapper using the default labels A and B
function relabelSitesBipartite!(
            uc          :: U,
            seed_site_A :: Int64 = 1
        ) where {LS,D,S<:AbstractSite{LS,D},B,U<:AbstractUnitcell{S,B}}

    # call the general function
    return relabelSitesBipartite!(uc, getDefaultLabelA(LS), getDefaultLabelB(LS), seed_site_A)
end

# export the function
export relabelSitesBipartite!


text:
using StatsBase: sample
using DataFrames: nrow

function select_row_sampling_strategy(subsample)
    if 0 < subsample < 1
        row_sampling_bytree_strategy =
            function(df, args...; kwargs...)
                rows = sample(1:nrow(df), round(Int, nrow(df)*subsample); replace = false)
            end
    elseif subsample == 1
        row_sampling_bytree_strategy =
            function(df, args...; kwargs...)
                df
            end
    else
        error("`subsample` must be within [0, 1)")
    end

    row_sampling_bytree_strategy
end


text:
struct Data
    turing_data
    cases
    deaths
    hospit
    num_observations
    num_total
    dates
    pop
    predictors
    seroprev
    regions
    lockdown
    observations_end
    cases_start
    startdates
end

# ============================================================================
#                      OBSERVABLES
# ============================================================================
# data cleaning

function replace_missing!(df, v=0)
    @assert !any( ismissing.( df.date ) )
    for c in eachcol(df)
        replace!(c, missing=>v)
    end
    disallowmissing!(df)
end

function insert_missing_dates!(df)
    ts = df.date[1] : Day(1) : df.date[end]
    missing_date = .! in.(ts, Ref(df.date))
    if any(missing_date)
        for t in ts[missing_date]
            insertrow!(df, t)
        end
    end
end

function insertrow!(df, d::Date, v=0)
    i = findfirst(==(d-Day(1)), df.date)
    v = [eltype(c) == Date ? d : v for c in eachcol(df)]
    insert!.(eachcol(df), i+1, v)
    nothing
end

# ============================================================================
# determine start and end date

function startdate(df; epidemicstart = 30, epidemic_threshold = 30, regions = names(df, Not(:date)))
    start_idx = zeros(Int64, length(regions))
    for (i, r) in enumerate(regions)
        cumulative = cumsum(df[:,r])
        start_idx[i] = findfirst(>=(epidemic_threshold), cumulative) #
    end
    df.date[start_idx] .- Day(epidemicstart)
end

function enddate(df, args...)
    isnothing(df) && enddate(args...)
    min(df.date[end], enddate(args...))
end

function enddate(df1, df2)
    t1 = isnothing(df1) ? Date(10^5) : df1.date[end]
    t2 = isnothing(df2) ? Date(10^5) : df2.date[end]
    min(t1, t2)
end
# ============================================================================
# make observables consistent
# limit observations and add zero padding if necessary

function consistent!(data; epidemicstart = 30, epidemic_threshold = 30, startdates = nothing)
    @unpack cases, hospit, deaths, regions = data
    if isnothing(startdates)
        startdates = startdate(hospit; epidemicstart, epidemic_threshold, regions)
    elseif startdates isa Date
        startdates = fill(startdates, length(regions))
        epidemicstart = 1 # start epidemic with given startdate
    else
        @assert startdates isa Vector{Date}
        epidemicstart = 1 # start epidemic with given startdates
    end
    data["epidemicstart"] = epidemicstart
    data["startdates"] = startdates
    s      = minimum(startdates)
    e      = enddate(deaths, cases, hospit)
    data["deaths"] = limit(deaths, s, e)
    data["hospit"] = limit(hospit, s, e)
    data["cases"]  = limit(cases, s, e)
    data["dates"]  = s:Day(1):e
    # @assert length(data["dates"]) == size(data["cases"],1) == size(data["hospit"],1) == size(data["deaths"],1)
    return nothing
end

# function consistent!(data; epidemicstart = Date("2020-06-01"))
#     @unpack cases, hospit, deaths, regions = data
#     data["epidemicstart"] = epidemicstart
#     s      = epidemicstart
#     e      = enddate(deaths, hospit, cases)
#     data["deaths"] = limit(deaths, s, e)
#     data["hospit"] = limit(hospit, s, e)
#     data["cases"]  = limit(cases, s, e)
#     data["dates"]  = s:Day(1):e
#     # @assert length(data["dates"]) == size(data["cases"],1) == size(data["hospit"],1) == size(data["deaths"],1)
#     return nothing
# end

function limit(df, s::Date, e::Date, pad=true)
    isnothing(df) && return
    pad && df.date[1] > s && ( df = padzeros(df, s) )
    is = findfirst(==(s), df.date)
    ie = findlast(==(e), df.date)
    df[is:ie, :]
end

function padzeros(df::DataFrame, s::Date)
    npad = (df.date[1] - s).value
    dfpad = DataFrame( :date => s:Day(1):df.date[end] )
    for n in names(df)
        n == "date" && continue
        dfpad[!,n] = padzeros(df[!,n], npad)
    end
    dfpad
end

function padzeros(xs::AbstractVector, padfirst=0, total=length(xs)+padfirst)
    n  = length(xs)
    ys = zeros(eltype(xs), total)
    ys[padfirst+1:n+padfirst] = xs
    ys
end

# ============================================================================
# transform observables to Turing format

function turingformat!(data, observationsend = nothing, predictive=false)
    @unpack dates, cases, hospit, deaths = data
    e = data["observationsend"] = isnothing(observationsend) ? last(dates) : Date(observationsend)
    s = data["startdates"]
    data["num_tot"]            = totaldays(dates, s)
    ts = data["dates_turing"]  = regional_timeseries(dates, s, e)
    cs = df2vec(cases, s, e)
    hs = df2vec(hospit, s, e)
    ds = df2vec(deaths, s, e)
    for (t,c,h,d) in zip(ts,cs,hs,ds)
        @assert length(t) == length(c) == length(h) == length(d)
    end

    data["cases_turing"]  = !predictive ? cs : [similar(c, Missing) for c in cs]
    data["hospit_turing"] = !predictive ? hs : [similar(h, Missing) for h in hs]
    data["deaths_turing"] = !predictive ? ds : [similar(d, Missing) for d in ds]
    data["num_observations"]  = length.(data["dates_turing"])
    return nothing
end

function df2vec(df, startdates, observationsend)
    vs = Vector{Vector{Int64}}()
    isnothing(df) && return vs
    regions = names(df, Not(:date))
    for (s,r) in zip(startdates, regions)
        is = findfirst(==(s), df.date)
        ie = findfirst(==(observationsend), df.date)
        v = vec(df[is:ie, r])
        push!(vs, v)
    end
    vs
end

function totaldays(dates, startdates)
    observationsend = last(dates)
    regional_dates = regional_timeseries(dates, startdates, observationsend)
    return length.(regional_dates)
end

function regional_timeseries(dates, startdates, observationsend)
    ts = Vector{Vector{Dates.Date}}()
    for s in startdates
        is = findfirst(==(s), dates)
        ie = findfirst(==(observationsend), dates)
        t = dates[is:ie]
        push!(ts, t)
    end
    ts
end

# ============================================================================
#                        COVARIATES
# ============================================================================
# main function
function covariates!(data, predictors; covariates_kwargs...)
    @unpack dates_turing, num_regions = data
    if isnothing(predictors)
        data["covariates"]       = [Array{Float64,2}(undef, 1,1) for i in 1:num_regions]
        data["num_covariates"]   = 0
        data["covariates_starts"] = fill(0, num_regions)
    else
        regional_covariates = Vector{Array{Float64,2}}()
        covariates_start_indices = Vector{Int64}()
        covariates_end_indices = Vector{Int64}()

        for i in 1:num_regions
            covariates = readcovariates(
                ; fname_survey = covariates_kwargs[:fname][i]
                , fname_mobility = covariates_kwargs[:mobility][i]
                , shift = covariates_kwargs[:shift]
                , datecol = covariates_kwargs[:datecol]
                #, aggregation = covariates_kwargs[:aggregation]
            )

            limit!(covariates, data, covariates_kwargs[:conditions])

            covariates_array = select_predictors(
                              covariates, predictors
                            ; standardize = covariates_kwargs[:standardize]
                            , normalize = covariates_kwargs[:normalize]
            ) #num_obs+shift_covariates, num_tot

            covariates_start_index = findfirst(==(covariates.date[1]), dates_turing[i])
            push!(regional_covariates, covariates_array)
            push!(covariates_start_indices, covariates_start_index)
        end
        data["covariates"] = regional_covariates
        @assert length(unique([size(regional_covariates[i],2) for i in num_regions])) == 1
        data["num_covariates"]   = size(regional_covariates[1],2)
        data["covariates_starts"] = covariates_start_indices
    end
    return nothing
end

# ============================================================================
# read covariates

function readcovariates(
        ; fname_survey=nothing
        , fname_mobility=nothing
        , shift=0
        , datecol=:date
        # , aggregation=1
        , kwargs...
)

    @_ fname_survey |>
        CSV.read(__, DataFrame) |>
        DataFrames.rename!(__, datecol=>"date") |>
        transform(__, :date =>( x -> x + Day(shift) )=> :date) |>
        leftjoin(__,
            CSV.read(fname_mobility, DataFrame); on=:date) |>
        sort(__, :date)
end

function limit!(covariates, data, conditions)
    @unpack observationsend = data
    filter!(:date => x->x<=Date(observationsend), covariates)
    if conditions isa Pair
        filter!(conditions, covariates)
    elseif conditions isa Vector && !isempty(conditions)
        filter!.(conditions, Ref(covariates))
    end
    return nothing
end

function consistencycheck!(covariates, data)
    @assert issorted(df.date)
    @assert !any( ismissing.(Array(df)) )
    disallowmissing!(df)
end

# ============================================================================
# preprocess covariates

function surveynorm!(xs; standardize=false)
    i = findfirst(x->!ismissing(x), xs)
    xs[i] == 0. && return nothing
    xs ./= abs(xs[i])
    xs .-= xs[i]
    standardize && (xs ./= std(xs))
    return nothing
end

function select_predictors(covariates, predictors, num_observations = nothing, num_tot = nothing; standardize=false, normalize = true)
    predictors = String.(predictors)

    cov_names = filter(x->x!=="date", names(covariates))
    for k in cov_names
        covariates[!,k] = convert.(Float64, covariates[!,k])
        normalize && surveynorm!(covariates[!,k]; standardize)
        # if isnothing(num_observations) || isnothing(num_tot)
        #     surveynorm!(covariates[!,k]; standartize)
        # else
        #     covariates[!,k] = preproc(covariates[:,k], num_observations, num_tot; normalize=surveynorm!)
        # end
    end

    miss = .!( in.(predictors, Ref(cov_names) ) )
    any(miss) && (@error "predictor $(predictors[miss]) not found. Select either nonresidential_mobility, residential_mobility or from survey: $(surveyvars)")
    covdata = Array( covariates[!,predictors] )
    return convert(Array{Float64,2}, covdata)
end

# ============================================================================

function serialinterval(nbins = 30)
    # Imperial Report
    r = vcat([0], 1.5:1:(nbins + .5))
    # diff(cdf.(GammaMeanStd(6.5, .62), r))
    # return p / sum(p)
    ## Brauer et al. Science (2020)
    p = diff(cdf.(GammaMeanStd(5.06, 2.11), r))
    return p# / sum(p)
end

# ============================================================================
#               RANDOM-WALK
# ============================================================================

function randomwalk!(data, stepsize, covariates_kwargs)
    @unpack num_observations, lockdown_indices, num_covariates = data
    ns = Vector{Int64}()
    is = Vector{Vector{Int64}}()
    for (i, (no, l)) in enumerate(zip(num_observations, lockdown_indices))
        n = if covariates_kwargs[:semiparametric] && data["num_covariates"] > 0
            data["covariates_starts"][i] - l - 1
        else
            no - l
        end
        i = stepindex(n, stepsize)
        push!(ns, length(unique(i)))
        push!(is, i)
    end
    data["num_rt_steps"] = ns
    data["rt_step_index"] = is
end

function stepindex(n, stepsize)
    index = Vector{Int64}(undef, n)
    step = 1
    for i in 1:n
        index[i] = step
        i % stepsize == 0 && ( step += 1 )
    end
    index
end

# ==============================================================================
# observation model

struct ObservParams{SA,SO,N,DD,DL,DP,P,H,W}
    starts::SA
    stops::SO
    num_regions::N
    delay_dist::DD
    delay_length::DL
    delay_dispersion::DP
    populations::P
    holidays::H
    weekdays::W
end

function ObservParams(data, o::ObservationInit)
    @unpack dates_turing, dates, num_regions = data
    ds  = [collect(first(dt):Day(1):last(dates)) for dt in dates_turing]

    params = (
        starts           = startindices(data, o),
        stops            = stopindices(data, o),
        num_regions      = num_regions,
        delay_dist       = o.delay_dist,
        delay_length     = o.delay_length,
        delay_dispersion = o.delay_dispersion,
        # observations     = observations(data, o),
        populations      = o.populations,
        holidays         = holiday.(ds),
        weekdays         = [dayofweek.(d) for d in ds]
    )
    return ObservParams( params... )
end

function startindices(data, o)
    @unpack dates_turing, num_regions = data
    observation_startdates = if isempty(o.obs_start)
        data["startdates"] .+ Day(data["epidemicstart"])
    else
        fill(Date(o.obs_start), num_regions)
    end

    is = Vector{Int64}()
    for m in 1:num_regions
        d = observation_startdates[m]
        i = findfirst(==(d), dates_turing[m])
        push!(is,i)
    end
    return is
end

function stopindices(data, o)
    @unpack dates_turing, num_regions = data
    if isempty(o.obs_stop)
        return data["num_observations"]
    else
        is = Vector{Int64}()
        d = Date(o.obs_stop)
        for m in 1:num_regions
            i = findfirst(==(d), dates_turing[m])
            push!(is,i)
        end
        return is
    end
end

# ============================================================================

function holiday(dates)
    specialdays = [Date("2020-12-24"), Date("2020-12-25"), Date("2020-12-31"), Date("2021-01-01")]
    holidays = dates .∈ Ref(specialdays)
    convert(Vector{Int64}, holidays)
end

## ============================================================================
#               depricated
# =============================================================================

# function nconvolve(f1, f2, ts)
#     nbins = length(ts)
#     f(τ, t) = f1(τ) * f2(t - τ)
#     f3 = zeros(nbins)
#     for (i, t) in enumerate(ts)
#         integral, err = quadgk(x->f(x, t), 0, t, rtol=1e-8)
#         f3[i] = integral
#     end
#     return f3
# end
#
# function inf2hosp(nbins=50, ihr = 1)
#     # SSI: ihr ∈ 0.05%-0.5%
#     # derived from Verity et al. See parameters.jl for details
#     # incubatine time ~ gamma(5.1, .86) see Flaxman et al
#     # onset to hospitalization see Zhang et al. Lancet (2020):
#     # "Evolving epidemiology and transmission dynamics of coronavirus disease 2019 outside Hubei province, China: a descriptive and modelling study"
#     ts = 1.5:1:(nbins + .5)
#     i2c(x) = pdf(GammaMeanCv(5.1, .86), x)
#     c2h(x) = pdf(GammaMeanCv(2.6, .92), x) * ihr
#     p =  nconvolve(i2c, c2h, ts)
#     return p / sum(p)
# end
#
# function inf2death(nbins=60, ihr = 1)
#     # ifr = 89 / 100_000 for age < 70 see Erikstrup et al. in Clin. Inf. Dis.
#
#     ## report34: onset to death = 19.8, Cv = .85
#     ts = 1.5:1:(nbins + .5)
#     i2c(x) = pdf(GammaMeanCv(5.1, .86), x)
#     c2d(x) = pdf(GammaMeanCv(19.8, .85), x) * ihr
#     return nconvolve(i2c, c2d, ts)
#
#     # Brauer et al.
#     # r = vcat([0], 1.5:1:(nbins + .5))
#     # p = diff(cdf.(NegativeBinomial2(21.82, 14.26), r))
#     # return p / sum(p)
#
#     ## S. Abbott
#     # i2o(x) = pdf(LogNormal(1.621, 0.418),x)
#     # o2d(x) = pdf(LogNormalMeanStd(13.1, 11.7),x)
#     # return nconvolve(i2o, o2d, ts)
#
#     ## Nouvelette
#     # GammaMeanStd(18.8, 8.48)
# end
#
# function inf2case(nbins = 30)
#     # Brauner et al.
#     r = vcat([0], 1.5:1:(nbins + .5))
#     p = diff(cdf.(NegativeBinomial2(10.92, 5.41), r))
#     return p / sum(p)
#     ## Imperial Report
#     # r = vcat([0], 1.5:1:(nbins + .5))
#     # return diff(cdf.(gamma(5.1, .86), r))
#     ## Brauner et al.
#     # return diff(cdf.(NegativeBinomial2(10.92, 5.41), r))
#     ## S. Abbott
#     # i2o(x) = pdf(LogNormal(1.621, 0.418),x)
#     # o2c(x) = pdf(LogNormalMeanStd(6.5, 17),x)
#     # return nconvolve(i2o, o2c, ts)
# end
#
# """
# seropos()[t] is the probablity of being seropositive t days after infection.
# For details see imperial report 34
# """
# function seropos(nbins=200)
#     # λ = 13.3
#     # μ = 300 #143.7 - λ
#     # κ = 3.67
#     # f(x, t, κ, μ, λ) = 1/λ*exp(-x/λ-((t-x)/μ)^κ)
#     # c2p(t) = quadgk(x->f(x, t, κ, μ, λ), 0, t, rtol=1e-8)[1]
#     # i2c(t) = pdf(GammaMeanCv(5.1, .86), t)
#     # ts = 1.5:1:(nbins + .5)
#     # p = nconvolve(i2c, c2p, ts)
#     # return p / sum(p)
#
#     i2s = zeros(nbins)
#     delay = 14
#     i2s[delay:end] .= 1.
#     return i2s
# end

# function padlast(xs, l)
#     n  = length(xs)
#     i  = findlast(x -> typeof(x) <: Real, xs) #avoid catching a missing value
#     ys = Vector{Union{Float64, Missing}}(undef, l)
#     ys[1:n] = xs
#     ys[n+1:end] .= xs[i]
#     ys
# end
#
# function preproc(xs::AbstractVector, num_observations, num_tot; normalize=nothing)
#     mask = 1:num_observations
#     xs = padlast(xs[mask], num_tot)
#     !isnothing(normalize) && normalize(xs)
#     xs
# end

# function cases_start_date!(data, cases_start)
#     @unpack dates_turing, num_observations = data
#     if isnothing(cases_start)
#         data["num_case_obs"] = zeros(Int64, length(regions))
#         data["cases_start_idx"] = nothing
#     else
#         is = Vector{Int64}()
#         ns = Vector{Int64}()
#         for (d, no) in zip(dates_turing, num_observations)
#             i  = findfirst(==(Date.(cases_start)), d)
#             n = no - i + 1
#             push!(is, i)
#             push!(ns, n)
#         end
#         data["cases_start_idx"] = is
#         data["num_case_obs"] = ns
#     end
#     return nothing
# end
#
#
# function rw_step_idx(num_observations, start_idx, iar_step)
#     # iar_idx = zeros(Int64, num_tot)
#     # iar_idx[1:start_idx-1] .= 1
#     idx = 2
#     for (j,i) in enumerate(start_idx:num_observations)
#          iar_idx[i] = idx
#          j % iar_step == 0 && (idx += 1)
#     end
#     if num_observations < num_tot
#         iar_idx[num_observations+1:num_tot] .= iar_idx[num_observations]
#     end
#     iar_idx
# end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
push!(ARGS, "../../input_files/dynamic/basin_refinement/16GPa/50n.dat")
include("../../Basin.jl")