text:
# Julia wrapper for header: /usr/local/include/libavutil/rational.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_cmp_q,
    av_q2d,
    av_reduce,
    av_mul_q,
    av_div_q,
    av_add_q,
    av_sub_q,
    av_inv_q,
    av_d2q,
    av_nearer_q,
    av_find_nearest_q_idx

function av_cmp_q(a::AVRational,b::AVRational)
    ccall((:av_cmp_q,libavutil),Cint,(AVRational,AVRational),a,b)
end

function av_q2d(a::AVRational)
    ccall((:av_q2d,libavutil),Cdouble,(AVRational,),a)
end

function av_reduce(dst_num,dst_den,num::Int64,den::Int64,max::Int64)
    ccall((:av_reduce,libavutil),Cint,(Ptr{Cint},Ptr{Cint},Int64,Int64,Int64),dst_num,dst_den,num,den,max)
end

function av_mul_q(b::AVRational,c::AVRational)
    ccall((:av_mul_q,libavutil),AVRational,(AVRational,AVRational),b,c)
end

function av_div_q(b::AVRational,c::AVRational)
    ccall((:av_div_q,libavutil),AVRational,(AVRational,AVRational),b,c)
end

function av_add_q(b::AVRational,c::AVRational)
    ccall((:av_add_q,libavutil),AVRational,(AVRational,AVRational),b,c)
end

function av_sub_q(b::AVRational,c::AVRational)
    ccall((:av_sub_q,libavutil),AVRational,(AVRational,AVRational),b,c)
end

function av_inv_q(q::AVRational)
    ccall((:av_inv_q,libavutil),AVRational,(AVRational,),q)
end

function av_d2q(d::Cdouble,max::Integer)
    ccall((:av_d2q,libavutil),AVRational,(Cdouble,Cint),d,max)
end

function av_nearer_q(q::AVRational,q1::AVRational,q2::AVRational)
    ccall((:av_nearer_q,libavutil),Cint,(AVRational,AVRational,AVRational),q,q1,q2)
end

function av_find_nearest_q_idx(q::AVRational,q_list)
    ccall((:av_find_nearest_q_idx,libavutil),Cint,(AVRational,Ptr{AVRational}),q,q_list)
end


text:
function visualize_def(grid, t, kw_args=Dict())
    println("def: ", myid())
end
function teeest(_)
	data = rand(10^6)
	for i=1:10
		sin(data+data)
	end
	visualize_def(rand(Float32, 10,10), 1)
end
println("start: ", myid())
t = Timer(teeest)
start_timer(t, 0, 1)
for i=1:10
	println("for: ", myid())
	sleep(1)
end
stop_timer(t)

ext: jl


text:
"""
    fetch_series(code::AbstractString, vintage::AbstractString)

Fetch series with series_id `code`.

Vintage can be any of `current` (alias `latest`, default), `previous` or a date in YYYY-MM-DD
form.

Returns a dataframe with four colums: `dates`, `vintage`, `series_id`, `values`.
"""
function fetch_series(code::AbstractString, vintage::AbstractString = "latest")
    fetch_series([code => vintage])
end

"""
    fetch_series(codes::AbstractArray{<:AbstractString, 1}, vintage::AbstractString)

Fetch a list of series with codes given by the array `codes`. For all series, fetch the
vintage given by `vintage`.

Vintage can be any of `current` (alias `latest`, default), `previous` or a date in YYYY-MM-DD
form.

Returns a dataframe with four colums: `dates`, `as_at_date`, `series_id`, `values`.
"""
function fetch_series(codes::AbstractArray{<:AbstractString, 1}, vintage::AbstractString = "latest")
    series = map(c -> c => vintage, codes)

    fetch_series(series)
end

"""
    fetch_series(series::AbstractArray{T, 1}) where T <: Pair{<:AbstractString, <:AbstractString}

Fetch a list of series with specified series codes and vintages.

Each element of the list is a pair `series_code` => `vintage`. Vintage must be explicitly specified.

Vintage can be any of `current` (alias `latest`), `previous` or a date in YYYY-MM-DD
form.

Returns a dataframe with four colums: `dates`, `as_at_date`, `series_id`, `values`.

# Example
```
fetch_series(["CPI" => "latest", "CPI_SYD" => "previous"])
````
"""
function fetch_series(series::AbstractArray{T, 1}) where T <: Pair{<:AbstractString, <:AbstractString}
    map(p -> begin
        check_valid_code(p[1])
        check_valid_vintage(p[2])
    end, series)

    series_req = join(
        map(p -> """{ "id" : "$(p[1])", "vintage" : "$(p[2])" }""", series),
        ", "
    )

    res = HTTP.request(
        "POST",
        "https://aeconomica.io/api/v1/fetchseries",
        [("Content-Type", "application/json")],
        """{
        "series": [
            $series_req
        ],
        "apikey" : "$(apikey())"}""",
        status_exception = false)
    if res.status == 200
        response = JSON.parse(String(res.body), null = missing)
        df = reduce(vcat, map(x -> DataFrames.DataFrame(x), response))
        df.dates = Dates.Date.(df.dates)
        df.vintage = Dates.Date.(df.vintage)
        df.values = if any(ismissing.(df.values))
            Vector{Union{Float64, Missing}}(df.values)
        else
            Vector{Float64}(df.values)
        end

        return df
    else
        error = JSON.parse(String(res.body))["error"]
        if error[1:21] == "500 Internal Error - "
            error = error[22:end]
            throw(ErrorException(error))
        else
            throw(ErrorException(error))
        end
    end
end


text:
# Automatically generated using Clang.jl

const CUDD_TRUE = 1
const CUDD_FALSE = 0
const CUDD_OUT_OF_MEM = -1
const CUDD_UNIQUE_SLOTS = 256
const CUDD_CACHE_SLOTS = 262144
const CUDD_RESIDUE_DEFAULT = 0
const CUDD_RESIDUE_MSB = 1
const CUDD_RESIDUE_TC = 2

# Skipping MacroDefinition: Cudd_Not ( node ) ( ( DdNode * ) ( ( uintptr_t ) ( node ) ^ ( uintptr_t ) 01 ) )
# Skipping MacroDefinition: Cudd_NotCond ( node , c ) ( ( DdNode * ) ( ( uintptr_t ) ( node ) ^ ( uintptr_t ) ( c ) ) )
# Skipping MacroDefinition: Cudd_Regular ( node ) ( ( DdNode * ) ( ( uintptr_t ) ( node ) & ~ ( uintptr_t ) 01 ) )
# Skipping MacroDefinition: Cudd_Complement ( node ) ( ( DdNode * ) ( ( uintptr_t ) ( node ) | ( uintptr_t ) 01 ) )
# Skipping MacroDefinition: Cudd_IsComplement ( node ) ( ( int ) ( ( uintptr_t ) ( node ) & ( uintptr_t ) 01 ) )
# Skipping MacroDefinition: Cudd_ReadIndex ( dd , index ) ( Cudd_ReadPerm ( dd , index ) )
# Skipping MacroDefinition: Cudd_ForeachCube ( manager , f , gen , cube , value ) for ( ( gen ) = Cudd_FirstCube ( manager , f , & cube , & value ) ; Cudd_IsGenEmpty ( gen ) ? Cudd_GenFree ( gen ) : CUDD_TRUE ; ( void ) Cudd_NextCube ( gen , & cube , & value ) )
# Skipping MacroDefinition: Cudd_ForeachPrime ( manager , l , u , gen , cube ) for ( ( gen ) = Cudd_FirstPrime ( manager , l , u , & cube ) ; Cudd_IsGenEmpty ( gen ) ? Cudd_GenFree ( gen ) : CUDD_TRUE ; ( void ) Cudd_NextPrime ( gen , & cube ) )
# Skipping MacroDefinition: Cudd_ForeachNode ( manager , f , gen , node ) for ( ( gen ) = Cudd_FirstNode ( manager , f , & node ) ; Cudd_IsGenEmpty ( gen ) ? Cudd_GenFree ( gen ) : CUDD_TRUE ; ( void ) Cudd_NextNode ( gen , & node ) )
# Skipping MacroDefinition: Cudd_zddForeachPath ( manager , f , gen , path ) for ( ( gen ) = Cudd_zddFirstPath ( manager , f , & path ) ; Cudd_IsGenEmpty ( gen ) ? Cudd_GenFree ( gen ) : CUDD_TRUE ; ( void ) Cudd_zddNextPath ( gen , & path ) )

@cenum Cudd_ReorderingType::UInt32 begin
    CUDD_REORDER_SAME = 0
    CUDD_REORDER_NONE = 1
    CUDD_REORDER_RANDOM = 2
    CUDD_REORDER_RANDOM_PIVOT = 3
    CUDD_REORDER_SIFT = 4
    CUDD_REORDER_SIFT_CONVERGE = 5
    CUDD_REORDER_SYMM_SIFT = 6
    CUDD_REORDER_SYMM_SIFT_CONV = 7
    CUDD_REORDER_WINDOW2 = 8
    CUDD_REORDER_WINDOW3 = 9
    CUDD_REORDER_WINDOW4 = 10
    CUDD_REORDER_WINDOW2_CONV = 11
    CUDD_REORDER_WINDOW3_CONV = 12
    CUDD_REORDER_WINDOW4_CONV = 13
    CUDD_REORDER_GROUP_SIFT = 14
    CUDD_REORDER_GROUP_SIFT_CONV = 15
    CUDD_REORDER_ANNEALING = 16
    CUDD_REORDER_GENETIC = 17
    CUDD_REORDER_LINEAR = 18
    CUDD_REORDER_LINEAR_CONVERGE = 19
    CUDD_REORDER_LAZY_SIFT = 20
    CUDD_REORDER_EXACT = 21
end

@cenum Cudd_AggregationType::UInt32 begin
    CUDD_NO_CHECK = 0
    CUDD_GROUP_CHECK = 1
    CUDD_GROUP_CHECK2 = 2
    CUDD_GROUP_CHECK3 = 3
    CUDD_GROUP_CHECK4 = 4
    CUDD_GROUP_CHECK5 = 5
    CUDD_GROUP_CHECK6 = 6
    CUDD_GROUP_CHECK7 = 7
    CUDD_GROUP_CHECK8 = 8
    CUDD_GROUP_CHECK9 = 9
end

@cenum Cudd_HookType::UInt32 begin
    CUDD_PRE_GC_HOOK = 0
    CUDD_POST_GC_HOOK = 1
    CUDD_PRE_REORDERING_HOOK = 2
    CUDD_POST_REORDERING_HOOK = 3
end

@cenum Cudd_ErrorType::UInt32 begin
    CUDD_NO_ERROR = 0
    CUDD_MEMORY_OUT = 1
    CUDD_TOO_MANY_NODES = 2
    CUDD_MAX_MEM_EXCEEDED = 3
    CUDD_TIMEOUT_EXPIRED = 4
    CUDD_TERMINATION = 5
    CUDD_INVALID_ARG = 6
    CUDD_INTERNAL_ERROR = 7
end

@cenum Cudd_LazyGroupType::UInt32 begin
    CUDD_LAZY_NONE = 0
    CUDD_LAZY_SOFT_GROUP = 1
    CUDD_LAZY_HARD_GROUP = 2
    CUDD_LAZY_UNGROUP = 3
end

@cenum Cudd_VariableType::UInt32 begin
    CUDD_VAR_PRIMARY_INPUT = 0
    CUDD_VAR_PRESENT_STATE = 1
    CUDD_VAR_NEXT_STATE = 2
end

const CUDD_VALUE_TYPE = Cdouble
const DdNode = Cvoid
const DdNodePtr = Ptr{DdNode}
const DdManager = Cvoid
const DdGen = Cvoid
const DdApaDigit = UInt32
const DdApaNumber = Ptr{DdApaDigit}
const DdConstApaNumber = Ptr{DdApaDigit}
const DdTlcInfo = Cvoid
const DD_HFP = Ptr{Cvoid}
const DD_PRFP = Ptr{Cvoid}
const DD_AOP = Ptr{Cvoid}
const DD_MAOP = Ptr{Cvoid}
const DD_CTFP = Ptr{Cvoid}
const DD_CTFP1 = Ptr{Cvoid}
const DD_OOMFP = Ptr{Cvoid}
const DD_QSFP = Ptr{Cvoid}
const DD_THFP = Ptr{Cvoid}
const DD_TOHFP = Ptr{Cvoid}


text:
using Oolong
using Test
using Base.Threads

@testset "Oolong.jl" begin
    include("core.jl")
end


text:
a = DataVector[1.0, 2.0, 3.0]
v_a = [1.0, 2.0, 3.0]

b = dataeye(3, 3)
m_b = eye(3, 3)

#
# Transposes
#

@assert all(a' .== v_a')
@assert all(a'' .== v_a'') #'
@assert all(b' .== m_b')
@assert all(b'' .== m_b'') #'

#
# DataVector * DataMatrix
#

@assert all(a * b[1, :] .== v_a * m_b[1, :])

#
# DataMatrix * DataVector
#

@assert all(b * a .== m_b * v_a)
@assert all(vector(b * a) .== m_b * v_a)

#
# DataMatrix * DataMatrix
#

@assert all(b * b .== m_b * m_b)

#
# DataVector * DataMatrix w/ NA's
#

b[1, 1] = NA
res = a * b[1, :]
@assert all(isna(res[:, 1]))
@assert all(!isna(res[:, 2]))
@assert all(!isna(res[:, 3]))
res = a * b[2, :]
@assert all(!isna(res))

#
# DataMatrix w NA's * DataVector
#

res = b * a
@assert isna(res[1])
@assert !isna(res[2])
@assert !isna(res[3])

#
# DataMatrix * DataMatrix
#

res = b * b
# 3x3 Float64 DataMatrix:
#  NA   NA   NA
#  NA  1.0  0.0
#  NA  0.0  1.0
@assert isna(res[1, 1])
@assert isna(res[1, 2])
@assert isna(res[1, 3])
@assert isna(res[2, 1])
@assert !isna(res[2, 2])
@assert !isna(res[2, 3])
@assert isna(res[3, 1])
@assert !isna(res[3, 2])
@assert !isna(res[3, 3])

res = b * dataeye(3)
# 3x3 Float64 DataMatrix:
#   NA   NA   NA
#  0.0  1.0  0.0
#  0.0  0.0  1.0
@assert isna(res[1, 1])
@assert isna(res[1, 2])
@assert isna(res[1, 3])
@assert !isna(res[2, 1])
@assert !isna(res[2, 2])
@assert !isna(res[2, 3])
@assert !isna(res[3, 1])
@assert !isna(res[3, 2])
@assert !isna(res[3, 3])

res = dataeye(3) * b
# julia> dataeye(3) * b
# 3x3 Float64 DataMatrix:
#  NA  0.0  0.0
#  NA  1.0  0.0
#  NA  0.0  1.0
@assert isna(res[1, 1])
@assert !isna(res[1, 2])
@assert !isna(res[1, 3])
@assert isna(res[2, 1])
@assert !isna(res[2, 2])
@assert !isna(res[2, 3])
@assert isna(res[3, 1])
@assert !isna(res[3, 2])
@assert !isna(res[3, 3])

# Test row operations
dm = dataeye(6, 2)
rowmeans(dm)

# Test column operations
dm = dataeye(6, 2)
colmeans(dm)

# Test linear algebra
du, dd, dv = svd(dataeye(3, 3))
u, d, v = svd(eye(3, 3))
@assert all(du .== u)
@assert all(dd .== d)
@assert all(dv .== v)

# Test elementary functions
dm = -dataeye(5, 5)
@assert all(abs(dm) .== eye(5, 5))


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*-----------------------------------------------------------------------
#* Author: Tiago Os򱨯rio
#* $Id$
#*---------------------------------------------------------------------
type ISE
	ISE()=begin
		[
			:(output=0*"s"),
		],
		[
			"",
		],
		new(
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"input signal"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"output signal",
				:Unit=>"s"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"setpoint"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"error^2"
			)),
			[
				:(sqError = (setPoint-input)^2),
				:(diff(output)=sqError),
			],
			[
				"Square error definition","Calculate output",
			],
			[:input,:output,:setPoint,:sqError,]
		)
	end
	input::DanaReal
	output::DanaReal
	setPoint::DanaReal
	sqError::DanaReal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	initials::Array{Expr,1}
	initialNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export ISE
function setEquationFlow(in::ISE)
	addEquation(1)
	addEquation(2)
end
function initial(in::ISE)
	addEquation(1)
end
function atributes(in::ISE,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Icon]="icon/ISE"
	fields[:Brief]="Model ISE."
	fields[:Info]="== Inputs ==
* One input signal.
* One setpoint signal.

== Outputs ==
* One output signal.
"
	drive!(fields,_)
	return fields
end
ISE(_::Dict{Symbol,Any})=begin
	newModel=ISE()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using RegistryCI
using Documenter

DocMeta.setdocmeta!(RegistryCI, :DocTestSetup, :(using RegistryCI); recursive=true)

makedocs(;
    modules=[RegistryCI],
    authors="Dilum Aluthge <dilum@aluthge.com>, Fredrik Ekre <ekrefredrik@gmail.com>, contributors",
    repo="https://github.com/JuliaRegistries/RegistryCI.jl/blob/{commit}{path}#{line}",
    sitename="RegistryCI.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaRegistries.github.io/RegistryCI.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Automatic merging guidelines" => "guidelines.md",
        "Regexes" => "regexes.md",
        "Using RegistryCI on your own package registry" => "private-registries.md",
        "Public API" => "public.md",
        "Internals (Private)" => "internals.md",
    ],
    strict=true,
)

deploydocs(;
    repo="github.com/JuliaRegistries/RegistryCI.jl",
)


text:
struct Groups{NG} <: AbstractVector{Int8}
    groups::Vector{Int8}
    Groups{NG}(::UndefInitializer, N) where NG = new(Vector{Int8}(undef, N))
end
Base.length(g::Groups{NG}) where NG = length(g.groups)
Base.size(g::Groups{NG}) where NG = size(g.groups)
@inline Base.getindex(g::Groups, i) = @inbounds g.groups[i]
@inline Base.setindex!(g::Groups, v, i) = g.groups[i] = v
@inline Base.IndexStyle(::Groups) = IndexLinear()
@inline Base.pointer(g::Groups) = pointer(g.groups)
@inline Base.resize!(g::Groups, i::Int) = resize!(g.groups, i)

# Should these functions be located here, or in rng.jl?

@generated function Random.rand!(rng::AbstractRNG, group::Groups{NG}, unifs::Vector{T},
                                                    probabilities::AbstractMatrix{T}) where {NG,T}
    W = LoopVectorization.pick_vector_width(T)
    Vi = Vec{W,Int8}
    Vf = Vec{W,T}
    p_NG = Symbol(:p_, NG)
    vp_NG = Symbol(:vp_, NG)
    vg_NG = Symbol(:vg_, NG)
    quote
        N = length(group)
        rand!(rng, unifs)
        @nexprs $NG g -> vg_g = vbroadcast($Vi, g)
        pg = pointer(group)
        for n ∈ 1:($W):N+$(1-W)
            vp_1 = vload($Vf, probabilities, n)
            @nexprs $(NG-1) g -> vp_{g+1} = SIMDPirates.vadd(vp_g, vload($Vf, probabilities, (n,g+1)))
            vu = SIMDPirates.vmul(vload($Vf, unifs, n), $vp_NG)
            vg = $vg_NG
            @nexprs $(NG-1) g -> vg = vifelse( SIMDPirates.vless(vu, vp_{$NG-g}), vg_{$NG-g}, vg )
            vstore!(pg - 1 + n, vg)
        end
        for n ∈ N+1-(N % $W):N
            p_1 = probabilities[n]
            @nexprs $(NG-1) g -> p_{g+1} = p_g + probabilities[n,g+1]
            ptup = @ntuple $NG p
            u = unifs[n] * $p_NG
            g = Int8(1)
            while u > ptup[g]
                g += Int8(1)
            end
            group[n] = g
        end
    end
end

@generated function Random.rand!(rng::AbstractRNG, group::Groups{NG}, unifs::AbstractVector{T},
                                                    probabilities::AbstractVector{T}) where {NG,T}
    W = LoopVectorization.pick_vector_width(T)
    Vi = Vec{W,Int8}
    Vf = Vec{W,T}
    p_NG = Symbol(:p_, NG)
    vp_NG = Symbol(:vp_, NG)
    vg_NG = Symbol(:vg_, NG)
    quote
        N = length(group)
        rand!(rng, unifs)
        @nexprs $NG g -> vg_g = vbroadcast($Vi, g)

        p_1 = probabilities[1]
        @nexprs $(NG-1) g -> p_{g+1} = p_g + probabilities[g+1]
        @nexprs $NG g -> vp_g = vbroadcast($Vf, p_g)
        pg = pointer(group)
        for n ∈ 0:($W):N-$W
            vu = SIMDPirates.vmul(vload($Vf, unifs, n), $vp_NG)
            vg = $vg_NG
            @nexprs $(NG-1) g -> vg = vifelse( SIMDPirates.vless(vu, vp_{$NG-g}), vg_{$NG-g}, vg )
            vstore!(pg + n, vg)
        end
        r = N % $W
        if r > 0
            ptup = @ntuple $NG p
            for n ∈ N+1-r:N
                u = unifs[n] * $p_NG
                g = Int8(1)
                while u > ptup[g]
                    g += Int8(1)
                end
                group[n] = g
            end
        end
    end
end

loop_max_expr(M, W, N) = :($N - $((M-1)*W+1))
function remloop_min_expr(M, W, N)
    NpW = gensym(:NpW)
    :($NpW = $N + $(W - 1); $NpW - ($NpW & $(W * M - 1)))
end
@generated function Random.rand!(rng::VectorizedRNG.PCG{M}, group::Groups{NG}, probabilities::AbstractVector{T}) where {M,NG,T}
    W = LoopVectorization.pick_vector_width(T)
    Vi = Vec{W,Int8}
    Vf = Vec{W,T}
    p_NG = Symbol(:p_, NG)
    vp_NG = Symbol(:vp_, NG)
    vg_NG = Symbol(:vg_, NG)
    quote

        @nexprs $NG g -> vg_g = vbroadcast($Vi, g)
        N = length(group)

        p_1 = probabilities[1]
        @nexprs $(NG-1) g -> p_{g+1} = p_g + probabilities[g+1]
        @nexprs $NG g -> vp_g = vbroadcast($Vf, p_g)
        pg = pointer(group)
        # for n ∈ 0:$(M*W):N-$(M*W)
        for n ∈ 0:$(M*W):$(loop_max_expr(M, W, :N))
            vufull = rand(rng, Vec{$(M*W),$T})
            Base.Cartesian.@nexprs $M u -> begin
                vu = Base.Cartesian.@ntuple $W w -> @inbounds vufull[w + (u-1)*$W]
                vu = SIMDPirates.vmul( vu, $vp_NG)
                vg = $vg_NG
                @nexprs $(NG-1) g -> vg = vifelse( SIMDPirates.vless(vu, vp_{$NG-g}), vg_{$NG-g}, vg )
                vstore!(pg + n + (u-1)*$W, vg)
            end
        end
        # Requires extra space beyond the loop to be allocated.
        for n ∈ $(remloop_min_expr(M, W, :N)):$W:N-1
            vu = rand(rng, Vec{$W,$T})
            # @show typeof(vu)
            # @show typeof($vp_NG)
            vu = SIMDPirates.vmul(vu, $vp_NG)
            vg = $vg_NG
            @nexprs $(NG-1) g -> vg = vifelse( SIMDPirates.vless(vu, vp_{$NG-g}), vg_{$NG-g}, vg )
            vstore!(pg + n, vg)
        end
    end
end

@generated function Random.rand!(rng::VectorizedRNG.PCG{M}, group::Groups{NG}, probabilities::AbstractMatrix{T}) where {M,NG,T}
    W = LoopVectorization.pick_vector_width(T)
    Vi = Vec{W,Int8}
    Vf = Vec{W,T}
    p_NG = Symbol(:p_, NG)
    vp_NG = Symbol(:vp_, NG)
    vg_NG = Symbol(:vg_, NG)
    quote
        @nexprs $NG g -> vg_g = vbroadcast($Vi, g)
        N = length(group)

        pg = pointer(group)
        for n ∈ 0:$(M*W):$(loop_max_expr(M, W, :N))
            vufull = rand(rng, Vec{$(M*W),$T})
            Base.Cartesian.@nexprs $M u -> begin
                vp_1 = vload($Vf, probabilities, n + (u-1)*$W + 1)
                @nexprs $(NG-1) g -> vp_{g+1} = SIMDPirates.vadd(vp_g, vload($Vf, probabilities, (n+(u-1)*$W,g+1)))
                vu = SIMDPirates.vmul((Base.Cartesian.@ntuple $W w -> @inbounds vufull[w + (u-1)*$W]), $vp_NG)
                vg = $vg_NG
                @nexprs $(NG-1) g -> vg = vifelse( SIMDPirates.vless(vu, vp_{$NG-g}), vg_{$NG-g}, vg )
                vstore!(pg + n + (u-1)*$W, vg)
            end
        end
        # back tracking loop. This approach imposes a minimum sample size limit.
        # for n ∈ (N - (N & $W) - ((N & $(W-1)) > 0) * $W):$W:N-$W
        # because we preallocate extra space, we can instead extend past the end without
        # segfaulting. This however forces us to intialize probs data as zeros, to avoid subnormals.
        for n ∈ $(remloop_min_expr(M, W, :N)):$W:N-1
            vp_1 = vload($Vf, probabilities, n)
            @nexprs $(NG-1) g -> vp_{g+1} = SIMDPirates.vadd(vp_g, vload($Vf, probabilities, (n,g+1)))
            vu = SIMDPirates.vmul(rand(rng, Vec{$W,$T}), $vp_NG)
            vg = $vg_NG
            @nexprs $(NG-1) g -> vg = vifelse( SIMDPirates.vless(vu, vp_{$NG-g}), vg_{$NG-g}, vg )
            vstore!(pg + n, vg)
        end
    end
end

Random.rand!(g::Groups, α::AbstractArray) = Random.rand!(Random.GLOBAL_RNG, g, similar(α, size(α,1)), α)
function Random.rand!(g::Groups, unifs::AbstractVector, α::AbstractArray)
        Random.rand!(Random.GLOBAL_RNG, g, unifs, α)
end
Random.rand!(g::Groups) = Random.rand!(Random.GLOBAL_RNG, g)
function Random.rand!(rng::AbstractRNG, g::Groups{NG}) where NG
    @inbounds for n ∈ eachindex(g.groups)
        g[n] = rand(rng, Int8(1):Int8(NG))
    end
end
function Random.rand(::Type{Groups{NG}}, α::AbstractMatrix) where NG
    Random.rand!(Random.GLOBAL_RNG, Groups{NG}(undef,size(α,1)), α)
end
function Random.rand(rng::AbstractRNG, ::Type{Groups{NG}}, α::AbstractMatrix) where NG
    Random.rand!(rng, Groups{NG}(undef,size(α,1)), α)
end
function Random.rand(::Type{Groups{NG}}, N::Integer) where NG
    g = Groups{NG}(undef, N)
    rand!(g)
    g
end


text:
function _printconverged(io::IO, converged, niters)
  symb = converged ? "✓ " : "✗ "
  col  = converged ? :green : :red
  str  = converged ? "converged" : "could not converge"
  print(io, "  ")
  print_with_color(col, io,  symb)
  print(io, " ", str, " after ")
  print_with_color(:white, io, string(niters))
  println(io, " iterations")
end

function _printvariable(io, i, symbol, coef; colored = true, newline = true)
  len = i - length(string(symbol))
  print(io, repeat(" ", len))
  if colored
    print_with_color(:cyan, io, string(symbol))
  else
    print(io,  string(symbol))
  end
  print(io, " = ")
  if typeof(coef) <: AbstractArray
    if size(coef,1) <= 10 && size(coef,2) == 1
      if typeof(coef[1]) <: Float64
        print_with_color(:white, io, string(map(x -> round(x, 4), coef)))
      else
        print_with_color(:white, io, string(coef))
      end
    else
      print_with_color(:white, io, "$(size(coef,1))x$(size(coef,2)), $(typeof(coef))")
    end
  else
    print_with_color(:white, io, string(coef))
  end
  if newline; print(io, "\n") end
end


text:
using Test

include("const.jl")
include("downloads.jl")


text:
using CurveProximityQueries
using ConvexBodyProximityQueries
using IntervalArithmetic
using Plots
using StaticArrays
import Random: seed!

pgfplots()
seed!(22)

curve = Bernstein([[-0.6, 0.0], [-0.6, 0.0], [0.1, -0.4], [0.4, 1.6], [-0.2, -1.6], [-0.1, 0.4], [0.8, 0.0], [0.8, 0.0]])

poly = randpoly([0.3, 0.7]; n=20, scale=0.5)
pts = closest_points(curve(0.5), poly, @SVector(ones(2)))

plot(curve, 0.0..0.2,
     linewidth=2.0,
     linestyle=:dash,
     linecolor=Gray(0.5),
    )
plot!(curve, 0.2..0.8,
      linewidth=2.0,
      linecolor=:black,
     )
plot!(curve, 0.8..1.0,
      linewidth=2.0,
      linestyle=:dash,
      linecolor=Gray(0.5),
     )
scatter!(curve(0.2), markersize=4, markercolor=:black)
scatter!(curve(0.8), markersize=4, markercolor=:black)

plot!(poly,
      fillalpha=1.0,
      fillcolor=Gray(0.7),
     )
plot!(pts,
      linestyle=:dash,
      markersize=8.0,
      linewidth=2.0,
     )
savefig("out/fig5a.pdf")


text:
export amoeba_carcase_domain_heuristic, archimedean_tropical_polynomial, archimedean_tropical_curve, archimedean_neighbourhood, archimedean_tropical_hypersurface

"""
    archimedean_tropical_polynomial(f)

Compute the Archimedean tropical polynomial of `f`.
"""
function archimedean_tropical_polynomial(f::MP.AbstractPolynomial{<:Union{Real, Complex}})
    mapcoefficients(c -> Tropical(logabs.(c)), f)
end

"""
    archimedean_tropical_curve(f)

Compute the curve associated to the Archimedean tropical polynomial of `f`.
"""
function archimedean_tropical_curve(f::MP.AbstractPolynomial{<:Union{Real, Complex}})
    MP.nvariables(f) == 2 || throw(error("Only bivariate polynomials are currently supported"))
    TropicalCurve(archimedean_tropical_polynomial(f))
end
archimedean_tropical_curve(f::AbstractTropicalPolynomial) = TropicalCurve(f)

"""
    archimedean_tropical_curve(f)

Compute the hypersurface associated to the Archimedean tropical polynomial of `f`.
"""
function archimedean_tropical_hypersurface(f::MP.AbstractPolynomial{<:Union{Real, Complex}})
    MP.nvariables(f) == 2 || throw(error("Only bivariate polynomials are currently supported"))
    archimedean_tropical_curve(f)
end

"""
    amoeba_carcase_domain_heuristic(f; factor=1.5, aspect_ratio=:default)

Compute the boundary `(xmin, xmax, ymin, ymax)` of domain Ω such that
Ω ∩ Amoeba(f) is probably the carcase of Amoeba(f).
"""
function amoeba_carcase_domain_heuristic(f::MP.AbstractPolynomial; kwargs...)
    amoeba_carcase_domain_heuristic(archimedean_tropical_curve(f); kwargs...)
end
function amoeba_carcase_domain_heuristic(curve::TropicalCurve; factor=1.5, aspect_ratio=:default)
    t = MP.nterms(curve.polynomial)
    r = factor * log(t - 1)

    xmin, xmax = Inf, -Inf
    ymin, ymax = Inf, -Inf
    for node in curve.vertices
        x, y = node
        xmin = min(xmin, x - r)
        xmax = max(xmax, x + r)
        ymin = min(ymin, y - r)
        ymax = max(ymax, y + r)
    end

    if aspect_ratio == :equal
        return equal_aspect_ratio(xmin, xmax, ymin, ymax)
    end

    xmin, xmax, ymin, ymax
end

"""
    archimedean_neighbourhood(archtrop, grid)

Compute the neighbourhood of `archtrop` in which the amoeba of `archtrop.polynomial` is
guaranteed to lie. See Avendaño, Kogan, Nisse, Rojas [^1] for details.

[^1]: Avendano, Martin, et al. "Metric Estimates and membership complexity
for archimedean amoebae and tropical hypersurfaces." arXiv preprint arXiv:1307.3681 (2013).
"""
function archimedean_neighbourhood(curve::TropicalCurve, grid)
    t = MP.nterms(curve.polynomial)
    r = log(t - 1)
    nb = archimedean_neighbourhood_rectangles(curve, r, limits(grid)...)
    M = falses(size(grid)...)
    for (px, py) in nb
        fillrect!(M, xgridpoint.(grid, px), ygridpoint.(grid, py))
    end
    Bitmap2D(M', grid)
end
function archimedean_neighbourhood(f::MP.AbstractPolynomial, grid)
    archimedean_neighbourhood(archimedean_tropical_curve(f), grid)
end

"Compute the neighbourhoods with radius `r` around the curve"
function archimedean_neighbourhood_rectangles(curve::TropicalCurve, r, xmin, xmax, ymin, ymax)
    # This doesn't really work if we the ray starts from outside the visible area
    nbs = Vector{NTuple{2, SVector{4,Float64}}}()
    for (raynode, raydirection) in curve.halfrays
        a = curve.vertices[raynode]
        b = a + maxlength(a, xmin, xmax, ymin, ymax) * raydirection
        push!(nbs, _rectline(a, b, r))
    end

    for segment in curve.segments
        a, b = curve.vertices[segment]
        # in degenerate cases it can happen that vertices coincide
        if norm(a-b) > 1e-8
            push!(nbs, _rectline(a, b, r))
        end
    end
    nbs
end
@inline function maxlength(a::SVector{2}, xmin, xmax, ymin, ymax)
    max(
        norm(a - SVector(xmin, ymin)),
        norm(a - SVector(xmin, ymax)),
        norm(a - SVector(xmax, ymin)),
        norm(a - SVector(xmax, ymax))
        )
end

"A rectangle centered around the line from `a` to `b` with height `2r`"
function _rectline(a, b, r)
    v1, v2 = normalize(a - b)
    d1 = -r*v2
    d2 = r*v1
    (SVector(a[1] + d1, b[1] + d1, b[1] - d1, a[1] - d1),
     SVector(a[2] + d2, b[2] + d2, b[2] - d2, a[2] - d2))
end

# function Base.show(io::IO, mime::MIME"text/html", A::ArchTrop)
#     if inIJulia() && plotsdefined()
#         show(io, mime, Main.Plots.plot(curve(A)))
#     else
#         show(io, A)
#     end
# end
# function Base.show(io::IO, A::ArchTrop)
#     #println(typeof(A), ":")
#     println(io, "ArchTrop of $(A.polynomial) is $(A.curve.polynomial)")
# end


text:
using NLsolve

include("functions_endoL.jl")

hh = Household();
r_guess = 0.005;
A_guess = 0.45726596004725534;
cpol,apol,lpol = iterate_egm(hh,A_guess;r=r_guess)

plot(hh.Amat[:,1], hh.Amat[:,1], label = "", color = :black, linestyle = :dash)
plot!(hh.Amat[:,1], apol[:,1], label = "ϵ = $(round(hh.Ymat[1,1], digits = 3))")
plot!(hh.Amat[:,1], apol[:,5], label = "ϵ = $(round(hh.Ymat[1,5], digits = 3))")

plot(hh.Amat[:,1], lpol[:,1], label = "ϵ = $(round(hh.Ymat[1,1], digits = 3))")
plot!(hh.Amat[:,1], lpol[:,2], label = "ϵ = $(round(hh.Ymat[1,2], digits = 3))")

# Calibrating A such that Y = 1, A = 0.45726596004725534
"""
dist_Y = 10;
while dist_Y > 10^-5
    r,Y = market_clearing(hh,A_guess,r_guess);
    dist_Y = abs(Y-1);
    A_guess = 0.5*A_guess+0.5*A_guess/Y;
    println("Y is $Y, new A is $A_guess")
end
"""

#market_clearing(hh,A_guess,r_guess)
#plot_market_clearing(hh,A_guess)

function f!(F,x)
    r = x[1]
    A = x[2]
    r_eq,K,N = market_clearing(hh,A,r)
    F[1] = r_eq-r
    F[2] = 1-A*(K^0.3)*(N^(1-0.3))
end

nlsolve(f!,[0.015,0.4])


text:
using Base.Test
using Compat  # for breaking changes in julia 0.7
using SoilTracers

anyerrors = false

tests = ["consts.jl",
         "soilgrid/soilgrid.jl",
         "physchem/thermodyn.jl",
         "physchem/diffus.jl",
         "physchem/solub.jl",
         "physchem/water.jl",
         "physchem/rxn.jl",
         "physchem/physchem.jl"]

println("Running tests:")

for t in tests
    try
        include(t)
        println("\t\033[1m\033[32mPASSED\033[0m: $(t)")
    catch e
        global anyerrors = true
        println("\t\033[1m\033[31mFAILED\033[0m: $(t)")
        showerror(stdout, e, backtrace())
    end
end

if anyerrors
    throw(error("Tests failed"))
end


text:
#directory info

SRCDIR = dirname(pathof(TightlyBound))
EXAMPLESDIR = joinpath(dirname(pathof(TightlyBound)), "..", "examples")
TESTDIR = joinpath(dirname(pathof(TightlyBound)), "..", "test")
TEMPLATEDIR = joinpath(dirname(pathof(TightlyBound)), "..", "template_inputs")
STRUCTDIR = joinpath(dirname(pathof(TightlyBound)), "..", "reference_structures")

DATSDIR1 = joinpath(dirname(pathof(TightlyBound)), "..", "dats", "pbesol", "v1.1")

#DATSDIR1 = joinpath(dirname(pathof(TightlyBound)), "..", "dats", "pbesol", "v0.9")

#DATSDIR1 = joinpath(dirname(pathof(TightlyBound)), "..", "dats", "pbesol", "v0.4")

DATSDIR2 = joinpath(dirname(pathof(TightlyBound)), "..", "dats", "pbesol", "v1.0")

DOCSDIR = joinpath(dirname(pathof(TightlyBound)), "..", "docs")

#DATSDIR1 = "/home/kfg/codes/TB_fit/fit_dimer/datab_2/"
#DATSDIR2 = ""

#DATSDIR2 = "/home/kfg/codes/TB_fit/binary_v10/datab/"


text: include("projectednewton.jl")

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
"""
Returns a vector of default scenarios for output.
"""
function get_scenarios()::Vector{NetworkParameters}

    scenario1 = NetworkParameters(  L=3,
                                    gamma_scv = 3.0,
                                    λ = NaN,
                                    η = 4.0,
                                    μ_vector = ones(3),
                                    P = [0 1.0 0;
                                        0 0 1.0;
                                        0 0 0],
                                    Q = zeros(3,3),
                                    p_e = [1.0, 0, 0],
                                    K = fill(5,3))

    scenario2 = NetworkParameters(  L=3,
                                    gamma_scv = 3.0,
                                    λ = NaN,
                                    η = 4.0,
                                    μ_vector = ones(3),
                                    P = [0 1.0 0;
                                        0 0 1.0;
                                        0.5 0 0],
                                    Q = zeros(3,3),
                                    p_e = [1.0, 0, 0],
                                    K = fill(5,3))

    scenario3 = NetworkParameters(  L=3,
                                    gamma_scv = 3.0,
                                    λ = NaN,
                                    η = 4.0,
                                    μ_vector = ones(3),
                                    P = [0 1.0 0;
                                        0 0 1.0;
                                        0.5 0 0],
                                    Q = [0 0.5 0;
                                        0 0 0.5;
                                        0.5 0 0],
                                    p_e = [1.0, 0, 0],
                                    K = fill(5,3))

    scenario4 = NetworkParameters(  L=5,
                                    gamma_scv = 3.0,
                                    λ = NaN,
                                    η = 4.0,
                                    μ_vector = collect(5:-1:1),
                                    P = [0   0.5 0.5 0   0;
                                        0   0   0   1   0;
                                        0   0   0   0   1;
                                        0.5 0   0   0   0;
                                        0.2 0.2 0.2 0.2 0.2],
                                    Q = [0 0 0 0 0;
                                        1 0 0 0 0;
                                        1 0 0 0 0;
                                        1 0 0 0 0;
                                        1 0 0 0 0],
                                    p_e = [0.2, 0.2, 0, 0, 0.6],
                                    K = [-1, -1, 10, 10, 10])

    scenario5 = NetworkParameters(  L=20,
                                    gamma_scv = 3.0,
                                    λ = NaN,
                                    η = 4.0,
                                    μ_vector = ones(20),
                                    P = zeros(20,20),
                                    Q = diagm(3=>0.8*ones(17), -17=>ones(3)), # changed from Yoni's advice
                                    p_e = vcat(1, zeros(19)),
                                    K = fill(5,20))

    return [scenario1, scenario2, scenario3, scenario4, scenario5]
end

"""
Runs simulations with varied arrival rate where the state keeps track of jobs locations.
Then plots the mean number of jobs in the system and Proportion of jobs orbiting between queues wrt λ.
Also plots the emprical distribution of the sojurn time of jobs in the system varied by λ.
"""
function plot_simulation_summary(scenario::NetworkParameters;
                                                max_time::Real = 10.0^7,
                                                scenario_label::String = "",
                                                lambda_range = 1.0:5.0,
                                                save_folder::Union{Symbol, String} = :none)

    mean_jobs, proportions, all_durations = Float64[], Float64[], Vector{Float64}[]

    for λ in lambda_range
        prev_time, prev_total, prev_transit = 0.0, 0.0, 0.0
        total_integral, transit_integral = 0.0, 0.0
        durations = Float64[]

        function record(time::Float64, state::TrackedState)
            time_diff = time - prev_time
            total_integral += prev_total * time_diff
            transit_integral += prev_transit * time_diff

            prev_time = time
            prev_total = total_count(state)
            prev_transit = transit_count(state)

            while length(state.left_system) > 0
                job = pop!(state.left_system)
                push!(durations, duration(job))
            end
            return nothing
        end

        Random.seed!(0)
        simulate(NetworkParameters(scenario, λ),
                job_tracking = :times,
                max_time = max_time,
                call_back = record)

        push!(mean_jobs, total_integral / max_time)
        push!(proportions, (transit_integral / total_integral))
        push!(all_durations, durations)
    end

    p1 = plot(lambda_range, mean_jobs,
                    title = "Mean number of jobs in system $scenario_label",
                    xlabel = "λ (mean arrival rate)",
                    ylabel = "Mean num jobs",
                    label = false)
    if save_folder != :none
        savefig(p1, save_folder * "/mean $(scenario_label).png")
    end
    display(p1)

    p2= plot(lambda_range, proportions,
                    title = "Proportion of jobs in orbit $scenario_label",
                    xlabel = "λ (mean arrival rate)",
                    ylabel = "Proportion",
                    label = false)
    if save_folder != :none
        savefig(p2, save_folder * "/proportion $(scenario_label).png")
    end
    display(p2)

    p3 = plot(title = "Empirical CDF of sojurn time $scenario_label", xlabel = "Duration", ylabel = "Probability")
    times = 0:0.01:max_time
    for i in LinRange(1, length(lambda_range), 5)
        index = convert(Int, floor(i))
        n = length(all_durations[index])
        #plot!(sort(all_durations[index]), (1:n)./n, label = "λ = $(lambda_range[index])")
        plot!(sort(all_durations[index])[convert.(Int, floor.(LinRange(1,n, 100)))], 1:100, label = "λ = $(round(lambda_range[index], digits = 2))")
    end
    if save_folder != :none
        savefig(p3, save_folder * "/ecdf $(scenario_label).png")
    end
    display(p3)

end

"""
Plots the simulation summarys (see plot_simulation_summary()) of all default scenarios run over a long simulation.
"""
function run_default_sims(;lambda_range = 1.0:5.0, max_time = 10.0^4)

    for (i, scenario) in enumerate(get_scenarios()[1:4])
        plot_simulation_summary(scenario,
                                     max_time = max_time,
                                     lambda_range = lambda_range,
                                     scenario_label = "scenario $i")
    end
end

"""
Runs a short simulation of the given scenario printing the full state of the system.
"""
function run_tracking_sim(scenario::NetworkParameters, λ::Real;
                            max_time::Float64 = 10.0,
                            log_times::Vector{Float64} = [max_time - 10.0^(-10)])

        Random.seed!(0)
        simulate(NetworkParameters(scenario, λ),
                job_tracking = :full,
                max_time = max_time,
                call_back = (time, state) -> nothing,
                log_times = log_times)
end

"""
Runs simulations with varied arrival rate where the state only keeps track of queue and transit totals.
Then plots the mean number of jobs in the system and Proportion of jobs orbiting between queues wrt λ.
"""
function run_default_no_tracking(;lambda_range = 1.0:5.0, max_time = 10.0^4)

    for (i, scenario) in enumerate(get_scenarios()[1:4])
        mean_jobs, proportions = Float64[], Float64[]

        for λ in lambda_range
            prev_time, prev_total, prev_transit = 0.0, 0.0, 0.0
            total_integral, transit_integral = 0.0, 0.0

            function record(time::Float64, state::NetworkState)
                time_diff = time - prev_time
                total_integral += prev_total * time_diff
                transit_integral += prev_transit * time_diff

                prev_time = time
                prev_total = total_count(state)
                prev_transit = transit_count(state)
                return nothing
            end

            Random.seed!(0)
            simulate(NetworkParameters(scenario, λ),
                    job_tracking = :none,
                    max_time = max_time,
                    call_back = record)

            push!(mean_jobs, total_integral / max_time)
            push!(proportions, (transit_integral / total_integral))
        end

        display(plot(lambda_range, mean_jobs,
                        title = "Mean number of jobs in system scenario $i",
                        xlabel = "λ",
                        ylabel = "Mean num jobs",
                        label = false))

        display(plot(lambda_range, proportions,
                        title = "Proportion of jobs in orbit scenario $i",
                        xlabel = "λ",
                        ylabel = "Proportion",
                        label = false))

    end
end


text:
import Base.==

struct Length{U, T} <: Measure
    value::T
end
Length(unit::Symbol, x::T) where T = Length{unit, T}(x)

const AbsoluteLength = Length{:mm, Float64}

Base.convert(::Type{Length{U, T1}}, x::Length{U, T2}) where {U, T1 <: Number, T2 <: Number} =
    Length{U, T1}(x.value)

==(x::Length{U}, y::Length{U}) where U = x.value == y.value
Base.isequal(x::Length{U}, y::Length{U}) where U = isequal(x.value, y.value)
Base.hash(x::Length{U}) where U = hash(x.value, hash(U))

# Operations
# ----------

Neg(x::T) where T <: Length = T(-x.value)
Div(a::Length{U}, b::Number) where U = Length(U, a.value / b)
Div(a::Length{U}, b::Length{U}) where U = a.value / b.value
Mul(a::Length{U}, b::Number) where U = Length(U, a.value * b)
Max(a::Length{U}, b::Length{U}) where U = Length(U, max(a.value, b.value))
Min(a::Length{U}, b::Length{U}) where U = Length(U, min(a.value, b.value))

Base.:+(a::Length{U}, b::Length{U}) where U = Length(U, a.value + b.value)
Base.:-(a::Length{U}, b::Length{U}) where U = Length(U, a.value - b.value)

iszero(x::Length) = x.value == 0.0

Base.abs(a::T) where {T <: Length} = T(abs(a.value))
Base.isless(a::Length{U}, b::Length{U}) where U = a.value < b.value
Base.oneunit(a::Type{Length{U,T}}) where {U,T} = Length(U, one(T))

isabsolute(::Type{Length{U,T}}) where {U,T} = !in(U, [:w,:h])

# Constants
# ---------

const mm   = AbsoluteLength(1.0)
const cm   = AbsoluteLength(10.0)
const inch = AbsoluteLength(25.4)
const pt   = inch/72.0

const w    = Length{:w, Float64}(1.0)
const h    = Length{:h, Float64}(1.0)
const d    = Length{:d, Float64}(1.0)


text:
@testset "'SampleSize' ......................... " begin

  function simonsDesign(r1, n1, r, n)
    nvec = [[n1 for x1 in 0:r1]; [n for x1 in (r1 + 1):n1]]
    cvec = [[Inf for x1 in 0:r1]; [r for x1 in (r1 + 1):n1]]
    return Design(nvec, cvec)
  end

  # Simon's designs for beta = .2, alpha = .05, p1 = p0 +0.2
  # cf. Simon, R "Optimal Two-Stage Designs for Phase II Clinical Trials",
  # Controlled Clinical Trials 10:1-10 (1989). (p. 4)
  p0 = collect(linspace(.1, .7, 7))
  sd = [
    simonsDesign( 1, 10,  5, 29), # p0 = 0.1
    simonsDesign( 3, 13, 12, 43), # p0 = 0.2
    simonsDesign( 5, 15, 18, 46), # p0 = 0.3
    simonsDesign( 7, 16, 23, 46), # p0 = 0.4
    simonsDesign( 8, 15, 26, 43), # p0 = 0.5
    simonsDesign( 7, 11, 30, 43), # p0 = 0.6
    simonsDesign( 4,  6, 22, 27)  # p0 = 0.7
  ]

  en = [15.0, 20.6, 23.6, 24.5, 23.5, 20.5, 14.8]
  for i in 1:length(p0)
    ss = SampleSize(sd[i], p0[i])
    @test abs(mean(ss) - en[i]) < .1
    # should be sufficient to ensure correctness,
    # check whether all routines can be called
    nr = rand(ss, convert(Int64, 1e4))
    @test all(minimum(ss) .<= nr .<= maximum(ss))
    pdf(ss, 25)
    quantile(ss, .5)
  end

end

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
function post(alg::LGG09{N, AM, VN, TN}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
              Δt0::TimeInterval=zeroI, kwargs...) where {N, AM, VN, TN}

    @unpack δ, approx_model, template, static, threaded, vars = alg

    # dimension check
    @assert statedim(ivp) == dim(template) "the problems' dimension $(statedim(ivp)) " *
        "doesn't match the dimension of the template directions, $(dim(template))"

    if haskey(kwargs, :NSTEPS)
        NSTEPS = kwargs[:NSTEPS]
        T = NSTEPS * δ
    else
        # get time horizon from the time span imposing that it is of the form (0, T)
        T = _get_T(tspan, check_zero=true, check_positive=true)
        NSTEPS = ceil(Int, T / δ)
    end

    # normalize system to canonical form
    ivp_norm = _normalize(ivp)

    # homogeneize the initial-value problem
    if haskey(kwargs, :homogeneize) && kwargs[:homogeneize] == true
        ivp_norm = homogeneize(ivp_norm)
    end

    # discretize system
    ivp_discr = discretize(ivp_norm, δ, approx_model)
    Φ = state_matrix(ivp_discr)
    Ω₀ = initial_state(ivp_discr)
    X = stateset(ivp_discr)

    if alg.sparse # ad-hoc conversion of Φ to sparse representation
        Φ = sparse(Φ)
    end

    cacheval = Val(alg.cache)

    # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
    got_homogeneous = !hasinput(ivp_discr)

    # NOTE: option :static currently ignored!

    # preallocate output flowpipe
    SN = SubArray{N, 1, Matrix{N}, Tuple{Base.Slice{Base.OneTo{Int}}, Int}, true}
    RT = TemplateReachSet{N, VN, TN, SN}
    F = Vector{RT}(undef, NSTEPS)

    if got_homogeneous
        ρℓ = reach_homog_LGG09!(F, template, Ω₀, Φ, NSTEPS, δ, X, Δt0, cacheval, Val(alg.threaded))
    else
        U = inputset(ivp_discr)
        @assert isa(U, LazySet)
        ρℓ = reach_inhomog_LGG09!(F, template, Ω₀, Φ, NSTEPS, δ, X, U, Δt0, cacheval, Val(alg.threaded))
    end

    return Flowpipe(F, Dict{Symbol, Any}(:sfmat => ρℓ, :alg_vars => vars))
end


text:
const _GA         = MathConstants.γ       # because I don't like unicode...
const QUADGAMMA1 =  -2.40411380631918857 # ψ^(2)(1)
const HEXGAMMA1  = -24.88626612344087823 # ψ^(4)(1)

# special methods for cosh(mu(v,x))*(x^vf) and sinh(...).
@inline coshmuxv(v, x) = (exp2(v) + exp2(-v)*x^(2*v))/2
@inline sinhmuxv(v, x) = (exp2(v) - exp2(-v)*x^(2*v))/2

const TPCOEF = (1, 0, (pi^2)/6, 0, (7*pi^4)/360, 0, (31*pi^6)/15120)
@inline tp_taylor0(v) = evalpoly(v, TPCOEF) # trig part

const G1COEF = (_GA,
                0,
                (2*_GA^3 - _GA*pi^2 - 2*QUADGAMMA1)/12,
                0,
                (12*_GA^5 - 20*(_GA^3)*pi^2 + _GA*pi^4 - 120*(_GA^2)*QUADGAMMA1 + 20*(pi^2)*QUADGAMMA1 - 12*HEXGAMMA1)/1440)
@inline g1_taylor0(v) = -evalpoly(v, G1COEF)

const G2COEF = (1.0,
                0,
                (_GA^2 - (pi^2)/6)/2,
                0,
                (60*_GA^4 - 60*(_GA*pi)^2 + pi^4 - 240*_GA*QUADGAMMA1)/1440)
@inline g2_taylor0(v) = evalpoly(v, G2COEF)

const SHCOEF = (1, 0, 1/6, 0, 1/120, 0, 1/5040, 0, 1/362880)
@inline sh_taylor0(v) = evalpoly(v, SHCOEF) # sinh part

# TODO: there is still a problem when v is not zero but z is zero. The NaN might
# be mathematically correct, and it isn't a branch that adbesselkxv hits.
@inline function f0_expansion0(v, z, modify)
  _t  = tp_taylor0(v)
  _g1 = g1_taylor0(v)
  _g2 = g2_taylor0(v)
  if !modify
    mu  = v*log(2/z)
    _cm = cosh(mu)
    _sm = sinh(mu)
    _sh = sh_taylor0(mu)*log(2/z)
  else
    _cm = coshmuxv(v, z)
    if iszero(v)
      # Because of the branching here, if v is zero, then I know that z is NOT zero.
      mu  = v*log(2/z)
      _sh = (z^v)*sh_taylor0(mu)*log(2/z)
    else
      _sh = sinhmuxv(v, z)/v
    end
  end
  _t*(_g1*_cm + _g2*_sh)
end

# EVEN Cheby coefs for computing the gamma pair thing when v is not near zero.
# (but |v| is <= 1/2).
const A2N = (1.843740587300906,
            -0.076852840844786,
             0.001271927136655,
            -0.000004971736704,
            -0.000000033126120,
             0.000000000242310,
            -0.000000000000170,
            -0.000000000000001
            )

# EVEN Cheby coefs for computing the gamma pair thing when v is not near zero.
# (but |v| is <= 1/2).
const A2Np1 = (-0.283876542276024,
                0.001706305071096,
                0.000076309597586,
               -0.000000865920800,
                0.000000001745136,
                0.000000000009161,
               -0.000000000000034)

# This gives g1 and g2 directly when v is not near zero and completely avoids
# calls to gamma functions.
@inline function temmegammas(v)
  twov = one(v)+one(v)
  _v   = twov*v
  (tv_even, tv_odd)   = (one(v), _v)
  (ser_even, ser_odd) = (A2N[1]/2, A2Np1[1]*tv_odd)
  @inbounds for n in 2:7
    # get next even value. Note that tv_even is now T_2(_v).
    tv_even   = twov*_v*tv_odd - tv_even
    # add to the even term series:
    ser_even += A2N[n]*tv_even
    # get the next odd term:
    tv_odd    = twov*_v*tv_even - tv_odd
    # add to the odd term series:
    ser_odd  += A2Np1[n]*tv_odd
  end
  # one more term for the evens:
  tv_even   = twov*_v*tv_odd - tv_even
  ser_even += A2N[8]*tv_even
  # now return:
  (ser_odd/v, ser_even)
end

function temme_pair(v, z, maxit, tol, modify=false)
  @assert abs(v) <= 1/2 "This internal routine is only for |v|<=1/2."
  # Some very low-level objects:
  onez = one(z)
  twoz = onez+onez
  zd2  = z/twoz
  _2dz = twoz/z
  # Creating the necessary things to get initial f, p, q:
  if abs(v) < 0.001
    g1 = g1_taylor0(v)
    g2 = g2_taylor0(v)
  else
    (g1, g2) = temmegammas(v)
  end
  (gp, gm) = (inv(-(g1*twoz*v - g2*twoz)/twoz), inv((g1*twoz*v + g2*twoz)/twoz))
  # p0 and q0 terms, branches for if we're modifying by (x^vf):
  if !modify
    p0 = (exp(-v*log(zd2))/twoz)*gp
    q0 = (exp(v*log(zd2))/twoz)*gm
  else
    p0 = exp2(v-one(v))*gp
    q0 = exp2(-v-one(v))*(z^(2*v))*gm
  end
  # cosh and sinh terms for f0, branches for if we're modifying by (x^vf):
  if !modify
     mu = v*log(_2dz)
    _cm = cosh(mu)
    _sm = sinh(mu)
  else
    _cm = coshmuxv(v, z)
    _sm = sinhmuxv(v, z)
  end
  # One more branch for f0, which is to check if v is near zero or z is near two.
  if iszero(z) && modify
    f0 = one(z)
  else
    if abs(v) < 0.001
      f0 = f0_expansion0(v, z, modify)
    elseif abs(z-2)<0.001
      _s = sh_taylor0(v*log(_2dz)) # manually plug in mu.
      #f0 = (v*pi/sinpi(v))*(g1*cosh(mu) + g2*log(_2dz)*_s)
      f0 = (v*pi/sinpi(v))*(g1*_cm + g2*log(_2dz)*_s)
    else
      # Temme's form is:
      #f0 = (v*pi/sinpi(v))*(g1*cosh(mu) + g2*log(_2dz)*sinh(mu)/mu)
      # But if I modify to this, I get rid of a near singularity as z->0:
      f0 = (v*pi/sinpi(v))*(g1*_cm + g2*_sm/v)
    end
  end
  (_f, _p, _q) = (f0, p0, q0)
  # a few other odds and ends for efficient looping:
  (ser_kv, ser_kvp1) = (f0, _p)
  (factk, _floatk)   = (onez, onez)
  (v2, _zd4, _z)     = (v*v, z*z/(twoz + twoz), z*z/(twoz + twoz))
  for k in 1:maxit
    _f   = (k*_f + _p + _q)/(_floatk^2 - v2)
    _p  /= (_floatk-v)
    _q  /= (_floatk+v)
    ck   =  _z/factk
    # update term for besselk(v, z).
    term_v    = ck*_f
    ser_kv   += term_v
    # update term for besselk(v+1, z).
    term_vp1  = ck*(_p - k*_f)
    ser_kvp1 += term_vp1
    if max(abs(term_v), abs(term_vp1)) < tol
      if !modify
        return (ser_kv, ser_kvp1*_2dz)
      else
        return (ser_kv, ser_kvp1*2) # note that I'm multiplying by a z, so cancel manually.
      end
    end
    _floatk += onez
    factk   *= _floatk
    _z      *= _zd4
  end
  throw(error("Term tolerance $tol not reached in $maxit iters."))
end

# NOTE: In the modified scaling of (x^v)*besselk(v,x), there is a problem for
# integer v: (x^0)*besselk(0,x) is just besselk(0,x). And that is still Inf for
# x=0. Which really breaks the whole strategy of integer derivatives for the
# rescaled Bessel here. BUT: there is a workaround! You don't actually need
# besselk(0, x) for the modified recurrence, so we just throw away that value.
function _besselk_temme(v, z, maxit, tol, mod)
  @assert v > -1/2 "This routine does not presently handle the case of v > -1/2."
  _p   = floor(v)
  (v - _p > 1/2) && (_p += one(_p))
  vf   = v - _p
  twov = one(v)+one(v)
  _v   = vf
  kvp2 = zero(v)
  (kv, kvp1) = temme_pair(_v, z, maxit, tol, mod)
  # check if any recurrence is necessary:
  v <= one(v)/twov          && return kv
  v <= (twov + one(v))/twov && return kvp1
  # if it is necessary, perform it:
  if mod
    # not necessarily the "right" way to handle this, but seems to stop the
    # propagation of NaNs in AD.
    #
    # a slightly different recurrence for (x^v)*besselk(v,x).
    if iszero(z)
      # special case for z=0:
      for _ in 1:(Int(_p)-1)
        kvp2 = twov*(_v+one(v))*kvp1
        _v  += one(v)
        kv   = kvp1
        kvp1 = kvp2
      end
    else
      z2 = z*z
      for _ in 1:(Int(_p)-1)
        kvp2 = twov*(_v+one(v))*kvp1 + z2*kv
        _v  += one(v)
        kv   = kvp1
        kvp1 = kvp2
      end
    end
  else
    for _ in 1:(Int(_p)-1)
      kvp2 = ((twov*(_v+one(v)))/z)*kvp1 + kv
      _v  += one(v)
      kv   = kvp1
      kvp1 = kvp2
    end
  end
  kvp2
end


text:
#========================================================================================#
#	Laboratory 9
#
# Working with noisy data.
#
# Author: Niall Palfreyman, 17/03/2022
#========================================================================================#
[
	Activity(
		"""
		GAs showed us how useful random processes can be for finding the solution to problems.
		Indeed, the entire field of biology is concerned with the way living organisms use
		random processes to solve problems. Notice this does not mean biological process are
		random: the breakdown of starch by ptyalin in your mouth follows a random path, yet it
		is absolutely certain that the starch will be broken down. However, we think of organisms
		as CHOOSING what to do. That is, they respond creatively to conditions by guiding their
		own random processes into EXPLORATORY processes that are neither entirely random nor
		entirely deterministic.

		To see how they do this, consider the bacterial virus lambda-phage (λ-phage). When
		λ-phage infects an E.coli bacterium, the infection can go in either of two different
		ways: lysis (rupturing the bacterium) or lysogeny (virus becomes integrated into the
		bacterium's DNA). Now comes the interesting part: Observed in a laboratory, it seems as
		if the lysis/lysogeny fate of the bacterium is purely random, driven by statistical
		variations in the concentration of certain transcription factors in the cell. This is
		very different from non-biochemical reactions, which either happen or they don't!

		Let's do some quick calculations. The E.coli cell is a cylinder with length 2μm and
		diameter 1μm. Use Julia to calculate the volume of this cell in litres:
		""",
		"Remember the formula for cylinder volume, and use the float notation 1e-6",
		x -> (1e-12 < x < 2e-12)
	),
	Activity(
		"""
		Transcription factors influnce dynamics in the cell at concentrations of around 1nM
		(i.e., nanomoles per litre). Use this value to calculate the total number of
		transcription factor molecules inside the E.coli cell:
		""",
		"Remember that you will need the Avogadro constant: 6e23 particles per mole",
		x -> (900 < x < 1000)
	),
	Activity(
		"""
		As you can see, there are only under 1000 molecules of the transcription factor inside the
		E.coli cell, so statistical variations play an important role in determining the
		decisions that the cell makes. Organisms form themselves in sizes and ways that can
		make use of random noise to explore multiple possible solutions.

		To model exploratory processes on a computer, we use PSEUDORANDOM numbers: numbers
		that are not truly random, because they are generated by a computer algorithm, but
		which are difficult to predict and uniformly distributed over some interval of
		numbers. What is the Julia function for producing random values?
		""",
		"",
		x -> x[1:4] == "rand"
	),
	Activity(
		"""
		Now we'll look at how uniformly distributed random numbers can generate patterns.
		Load the packages ElectronDisplay and CairoMakie, then create a list (Vector) of
		10 000 integers between 1 and 100. Now use the Makie function hist(data,bins=100) to
		create a histogram showing the distribution of theses numbers. Notice that the heights
		of the bars vary quite widely, but the whole histogram approximates a particular kind of
		statistical distribution. What word describes this kind of distribution?
		""",
		"Notice that although the heights vary, all bars are approximately the same height",
		x -> occursin("uniform",lowercase(x))
	),
	Activity(
		"""
		The heights of these bars vary around a mean value - what is this value?
		""",
		"",
		x -> x==100
	),
	Activity(
		"""
		Now create a list of 10 000 numbers between 1 and 100 by taking the mean of 100 integers
		between 1 and 100. The following line of code does this:

			data = [sum(rand(1:100,100))./100 for _ in 1:10000]

		Now draw a histogramm of this data: hist(data,bins=100). What kind of distribution is this?
		""",
		"",
		x -> occursin("normal",lowercase(x)) || occursin("gauss",lowercase(x))
	),
	Activity(
		"""
		So average values of a uniform distribution form a gaussian distribution. This is the
		content of the Central Limit Theorem of statistics, and is a technique that your own body
		uses to recognise accurately the signal of an antigen from amongst random noise.

		In a computer, we can use statistics to reduce noise. We shall now create a noise-free
		experiment by defining a simple linear function and a range of time and signal values:

			t = sort(rand(100));
			signal = 1.3t .- 0.3;

		Now plot this data and tell me the vertical intercept of the resulting plot:

			lines(t,signal)
		""",
		"Don't forget to load ElectronDisplay and CairoMakie first!",
		x -> x == -0.3
	),
	Activity(
		"""
		Now imagine we measured the data x and y, but didn't know the nature of the function f.
		If we know that f is a polynomial function, we can simply test the data against lots of
		possible polynomials to see which of them best fits our data. The CurveFit package does
		this job using least-squares fitting:

			using CurveFit
			fit = curve_fit(Polynomial, t, signal);

		What is the type of the curve_fit() result "fit"?
		""",
		"You may first need to open Package Manager (']') to add the package CurveFit",
		x -> x == Main.Polynomial
	),
	Activity(
		"""
		Now we can find the coefficients of our original signal by listing the components
		of fit: fit[:]. Which component of fit contains the vertical intercept of the signal?
		""",
		"Careful: Notice that Polynomials use more indices than you are used to!",
		x -> x==0
	),
	Activity(
		"""
		Right, now let's add some gaussian noise to our signal to create a noisy experiment:

			noisy = signal + 0.01randn(100);
			lines(t,noisy)

		What is your estimate of the gradient of the graph of noisy?
		""",
		"",
		x -> (1.2 < x < 1.4)
	),
	Activity(
		"""
		Now use curve_fit to tell me the approximate gradient of the noisy experiment:
		""",
		"Look back at how you performed curve-fitting earlier for the clean signal",
		x -> (1.27 < x < 1.33)
	),
	Activity(
		"""
		Our last topic in this laboratory is numerical differentiation. Something we often want
		to do is find the minimum value of a curve given just a few data points. For example,
		suppose we know that the following data has a minimum value somewhere:

			data = besselj.(6,9:.5:15);

		As a first approximation, use the findmin function to find the position of the minimum
		value. Which component of the data is the smallest?
		""",
		"",
		x -> x==6
	),
	Activity(
		"""
		But of course, our data has gaps in it between the datapoints. Is there maybe an even
		smaller value in one of these gaps? One way we can track down the minimum is to find the
		gradient of the data, then find out where this gradient is equal to zero. If t is the
		index of the datapoints in the array, we can define the gradient of the data by
		∂data≡Δ(data)/Δt. To calculate this, we subtract from each datapoint the value of the
		PREVIOUS datapoint, and then divide by the gap Δt=1 between consecutive indices of the
		data. We can do this using the circshift function:

			Δt = 1
			Δdata = data - circshift(data,1)
			∂data = Δdata / Δt

		The first components in your ∂data vector should be negative. How many negative components
		can you count, before they become positive?
		""",
		"To make the symbol Δ, type \\Delta then press <Tab>; for ∂, type \\partial, then <Tab>",
		x -> true
	),
	Activity(
		"""
		The circhshift function is an important tool, so let's take a moment to see how it works.
		circshift rotates the elements of an array circularly in the direction given by its second
		argument, for example:

			a = 1:5
			circshift(a,1)
			circhshift(a,-1)

		We can use circshift with arrays of any size. Enter the following lines and in each case,
		follow the shift in position of the number 1 from its position in the original array b.
		Notice how this shift corresponds to the second argument of circhshift:

			b = reshape(1:16,(4,4))
			circshift(b,(-1,0))
			circshift(b,(0,-2))
			circshift(b,(-1,-2))
		""",
		"",
		x -> true
	),
	Activity(
		"""
		OK, so our gradient ∂data looks good, but it has a problem. To calculate it, we subtracted
		the PREVIOUS datapoint from the current datapoint:

			∂data = (data(t)-data(t-Δt)) / Δt

		In the following graphical figure, we use this method to calculate the derivative at x=1:

			t = 0.1:.1:2.1
			lines( t, t.^2)
			lines!( t, t, color=:red)
			current_figure()

		At approximately which point between 0 and 1 is the curve parallel to the straight line?
		""",
		"It's about halfway, isn't it?",
		x -> x == 1.5
	),
	Activity(
		"""
		We just saw that subtracting the PREVIOUS datapoint from the current datapoint gives us the
		value of the gradient just BEFORE the current datapoint x. But instead, we could just as
		easily subtract the current datapoint from the NEXT datapoint:

			∂data = (data(t+Δt)-data(t)) / Δt

		The Mean Value Theorem of mathematics tells us that if we use two datapoints at
		t and t+Δt to calculate the gradient, this only gives us an accurate value for the
		gradient at ONE point somewhere between the two datapoints. To see this, enter the
		following lines of code:

			lines!( t, 3t.-2, color=:green)
			current_figure()

		Now, our new calculation gives us the value of the gradient just _______ the current
		datapoint t:
		""",
		"",
		x -> occursin("after", lowercase(x))
	),
	Activity(
		"""
		Look at the figure of your two gradient calculations. In both cases, the gradient is
		parallel to the curve at a point approximately halfway between the two datapoints. This
		idea is the basis of a much more accurate way of approximating the gradient! To calculate
		the approximate gradient AT the datapoint t, we look at the difference between the NEXT
		and the PREVIOUS datapoints:

			∂data(t) ≈ (data(t+Δt)-data(t-Δt)) / 2Δt

		You can see this idea in action by entering the following yellow line in your figure:

			lines!(x,2x,color=:yellow)
			current_figure()

		At approximately which value of t is the curve parallel to the straight line
		between 0 and 2?
		""",
		"It's about halfway, isn't it?",
		x -> x==1
	),
	Activity(
		"""
		Now let's apply our new trick to the data we created earlier. We'll use circshift to
		calculate at each datapoint the value (NEXT-PREVIOUS) / (2Δx):

			Δdata = circshift(data,-1) - circshift(data,1)
			∂data = Δdata / 2Δt

		Look at the component values of the gradient ∂data and make a guess: at what value of t do
		you estimate the minimum value of data to lie?
		""",
		"At approximately which value of t does the gradient ∂data change its sign?",
		x -> (6.4<x<6.6)
	),
	Activity(
		"""
		This idea of calculating gradients more accurately at the centre of an interval is the basis
		of the Runge-Kutta (R-K) method for solving differential equations (DEs) numerically. We
		will now implement the RK-method, and to test our implementation, we will need a simple test
		case. Consider the exponential equation with growth constant and initial value equal to 1:

			dx/dt = x;   x(0) = 1

		What would you expect to be the solution x(1) of this equation when t = 1?
		""",
		"",
		x -> abs(x-exp(1)) < 0.01
	),
	Activity(
		"""
		We can use our knowledge of the value of x(1) to test the goodness of the Euler and RK
		methods for solving a DE. First we set up Euler's method for solving the exponential
		equation over the interval from t=0 to t=1. To stress-test Euler's method, we shall only
		divide this interval into 4 time-steps:

			x = zeros(5); x[1] = 1; Δt = 0.25;
			m(x) = x
			for k in 2:5
				x[k] = x[k-1] + m(x[k-1]) * Δt
			end
			lines(0:4,x)

		You should see the standard increasing slope of the exponential function. What is the final
		value of x (after the four time-steps)?
		""",
		"You can see it by entering x[end]",
		x -> (2.4<x<2.5)
	),
	Activity(
		"""
		We would expect the last value of x to be something like e = 2.718. However, as you may
		already know, Euler's method suffers from the problem that its successive linear
		approximation steps tend to drift outwards away from the solution curve, making our graph
		inaccurate. Would this explanation lead you to expect your value of x[5] to be too low or
		to be too high?
		""",
		"The exponential function curves upwards, and Euler veers to the outside of this curve",
		x -> occursin("low",lowercase(x))
	),
	Activity(
		"""
		The RK-method reduces this problem greatly by first using Euler to approximate the gradient
		of the solution HALFWAY through the timestep, and then using this halfway value to
		approximate the gradient over a COMPLETE timestep. This usually gives excellent results:

			x = zeros(5); x[1] = 1; Δt = 0.25;
			m(x) = x
			for k in 2:5
				xHalfway = x[k-1] + m(x[k-1]) * Δt/2
				x[k] = x[k-1] + m(xHalfway) * Δt
			end
			lines(0:4,x)

		Again, you should see the increasing slope of the exponential function. What is your final
		value of x now (after the four time-steps)?
		""",
		"You can see it by entering x[end]",
		x -> true
	),
	Activity(
		"""
		As you can see, the RK result is much closer to the actual result of exp(1).

		Congratulations! You have implemented your very own DE solver!

		Before moving on to a well-earned rest, and to the next laboratory, take a couple of
		moments to change the number of time-steps in your solver and see how many you need in
		order to achieve the more accurate value of 2.718.
		""",
		"",
		x -> true
	),
]


text:
using InstantiateFromURL
# optionally add arguments to force installation: instantiate = true, precompile = true
github_project("QuantEcon/quantecon-notebooks-julia", version = "0.8.0")

using LinearAlgebra, Statistics
using QuantEcon, Plots
gr(fmt=:png);

# Exercise 1
ϕ_0, ϕ_1, ϕ_2 = 1.1, 0.8, -0.8

A = [1.0 0.0 0.0; ϕ_0 ϕ_1 ϕ_2; 0.0 1.0 0.0]
C = zeros(3,1)
G = [0 1 0]
μ_0 =  ones(3)

lss = LSS(A,C,G; mu_0 = μ_0)

x,y = simulate(lss, 50)
plot(dropdims(y, dims = 1), color = :blue)
plot!(xlabel="time", ylabel = "y_t", legend = :none)

# Exercise 2

using Random
Random.seed!(42) # For deterministic results.

ϕ1, ϕ2, ϕ3, ϕ4 = 0.5, -0.2, 0, 0.5
σ = 0.2

A = [ϕ1 ϕ2 ϕ3 ϕ4; 1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0]
C = [σ; 0.0; 0.0; 0.0]
G = [1.0 0.0 0.0 0.0]
μ_0 =  ones(4)

lss = LSS(A,C,G; mu_0 = μ_0)

x,y = simulate(lss, 200)
plot(dropdims(y, dims = 1), color = :blue)
plot!(xlabel="time", ylabel = "y_t", legend = :none)

# Exercise 3
ϕ1, ϕ2, ϕ3, ϕ4 = 0.5, -0.2, 0, 0.5
σ = 0.1

A = [ϕ1 ϕ2 ϕ3 ϕ4; 1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0]
C = [σ; 0.0; 0.0; 0.0]
G = [1.0 0.0 0.0 0.0]
μ_0 = ones(4)

I = 20
T = 50

ar = LSS(A, C, G; mu_0 = μ_0)

ensemble_mean = zeros(T)
ys = []
for ii in 1:I
        x,y = simulate(ar, T)
        y = dropdims(y, dims=1)
        push!(ys,y)
        ensemble_mean .+= y
end
ensemble_mean = ensemble_mean./I

plot(ensemble_mean, label = "y_t_bar", linewidth = 2)
plot!(ys, linewidth = 0.8, alpha = 0.2, label = "", color =:blue)

ymin, ymax = -0.5, 1.15
m = moment_sequence(ar)
pop_means = zeros(0)
for (i, t) ∈ enumerate(m)
    (μ_x, μ_y, Σ_x, Σ_y) = t
    push!(pop_means, μ_y[1])
end
plot!(pop_means, color = :green, linewidth = 2, label = "G mu_t")
plot!(ylims=(ymin, ymax), xlabel = "time", ylabel = "y_t", legendfont = font(12))

max_stars_repo_name: shanemcmiken/quantecon-notebooks-julia

max_issues_repo_name: shanemcmiken/quantecon-notebooks-julia

max_forks_repo_name: shanemcmiken/quantecon-notebooks-julia


text:
module RuntimeGeneratedFunctions

using ExprTools, Serialization, SHA

export RuntimeGeneratedFunction, @RuntimeGeneratedFunction

"""
    @RuntimeGeneratedFunction(function_expression)
    @RuntimeGeneratedFunction(context_module, function_expression, opaque_closures=true)

    RuntimeGeneratedFunction(cache_module, context_module, function_expression; opaque_closures=true)

Construct a function from `function_expression` which can be called immediately
without world age problems. Somewhat like using `eval(function_expression)` and
then calling the resulting function. The differences are:

* The result can be called immediately (immune to world age errors)
* The result is not a named generic function, and doesn't participate in
  generic function dispatch; it's more like a callable method.

You need to use `RuntimeGeneratedFunctions.init(your_module)` a single time at
the top level of `your_module` before any other uses of the macro.

If provided, `context_module` is module in which symbols within
`function_expression` will be looked up. By default this is module in which
`@RuntimeGeneratedFunction` is expanded.

`cache_module` is the module where the expression `code` will be cached. If
`RuntimeGeneratedFunction` is used during precompilation, this must be a module
which is currently being precompiled. Normally this would be set to
`@__MODULE__` using one of the macro constructors.

If `opaque_closures` is `true`, all closures in `function_expression` are
converted to
[opaque closures](https://github.com/JuliaLang/julia/pull/37849#issue-496641229).
This allows for the use of closures and generators inside the generated function,
but may not work in all cases due to slightly different semantics. This feature
requires Julia 1.7.

# Examples
```
RuntimeGeneratedFunctions.init(@__MODULE__) # Required at module top-level

function foo()
    expression = :((x,y)->x+y+1) # May be generated dynamically
    f = @RuntimeGeneratedFunction(expression)
    f(1,2) # May be called immediately
end
```
"""
struct RuntimeGeneratedFunction{argnames, cache_tag, context_tag, id} <: Function
    body::Expr
    function RuntimeGeneratedFunction(cache_tag, context_tag, ex; opaque_closures=true)
        def = splitdef(ex)
        args, body = normalize_args(def[:args]), def[:body]
        if opaque_closures && isdefined(Base, :Experimental) &&
            isdefined(Base.Experimental, Symbol("@opaque"))
            body = closures_to_opaque(body)
        end
        id = expr_to_id(body)
        cached_body = _cache_body(cache_tag, id, body)
        new{Tuple(args), cache_tag, context_tag, id}(cached_body)
    end

    # For internal use in deserialize() - doesen't check whether the body is in the cache!
    function RuntimeGeneratedFunction{argnames, cache_tag, context_tag, id}(body::Expr) where {argnames,cache_tag,context_tag,id}
        new{argnames, cache_tag, context_tag, id}(body)
    end
end

function _check_rgf_initialized(mods...)
    for mod in mods
        if !isdefined(mod, _tagname)
            error("""You must use `RuntimeGeneratedFunctions.init(@__MODULE__)` at module
                  top level before using runtime generated functions in $mod""")
        end
    end
end

function RuntimeGeneratedFunction(
        cache_module::Module, context_module::Module, code; opaque_closures=true,
    )
    _check_rgf_initialized(cache_module, context_module)
    RuntimeGeneratedFunction(
        getfield(cache_module, _tagname),
        getfield(context_module, _tagname),
        code;
        opaque_closures = opaque_closures
    )
end

macro RuntimeGeneratedFunction(code)
    quote
        RuntimeGeneratedFunction(@__MODULE__, @__MODULE__, $(esc(code)))
    end
end
macro RuntimeGeneratedFunction(context_module, code, opaque_closures=true)
    quote
        RuntimeGeneratedFunction(
            @__MODULE__, $(esc(context_module)), $(esc(code));
            opaque_closures = $(esc(opaque_closures)),
        )
    end
end

# Duplicate RuntimeGeneratedFunction docs onto @RuntimeGeneratedFunction
@eval @doc $(@doc RuntimeGeneratedFunction) var"@RuntimeGeneratedFunction"

function Base.show(io::IO, ::MIME"text/plain", f::RuntimeGeneratedFunction{argnames, cache_tag, context_tag, id}) where {argnames,cache_tag,context_tag,id}
    cache_mod = parentmodule(cache_tag)
    context_mod = parentmodule(context_tag)
    func_expr = Expr(:->, Expr(:tuple, argnames...), f.body)
    print(io, "RuntimeGeneratedFunction(#=in $cache_mod=#, #=using $context_mod=#, ", repr(func_expr), ")")
end

(f::RuntimeGeneratedFunction)(args::Vararg{Any,N}) where N = generated_callfunc(f, args...)

# We'll generate a method of this function in every module which wants to use
# @RuntimeGeneratedFunction
function generated_callfunc end

function generated_callfunc_body(argnames, cache_tag, id, __args)
    setup = (:($(argnames[i]) = @inbounds __args[$i]) for i in 1:length(argnames))
    body = _lookup_body(cache_tag, id)
    @assert body !== nothing
    quote
        $(setup...)
        $(body)
    end
end

### Body caching and lookup
#
# Looking up the body of a RuntimeGeneratedFunction based on the id is a little
# complicated because we want the `id=>body` mapping to survive precompilation.
# This means we need to store the mapping created by a module in that module
# itself.
#
# For that, we need a way to lookup the correct module from an instance of
# RuntimeGeneratedFunction. Modules can't be type parameters, but we can use
# any type which belongs to the module as a proxy "tag" for the module.
#
# (We could even abuse `typeof(__module__.eval)` for the tag, though this is a
# little non-robust to weird special cases like Main.eval being
# Base.MainInclude.eval.)

# It appears we can't use a ReentrantLock here, as contention seems to lead to
# deadlock. Perhaps because it triggers a task switch while compiling the
# @generated function.
_cache_lock = Threads.SpinLock()
_cachename = Symbol("#_RuntimeGeneratedFunctions_cache")
_tagname = Symbol("#_RGF_ModTag")

function _cache_body(cache_tag, id, body)
    lock(_cache_lock) do
        cache = getfield(parentmodule(cache_tag), _cachename)
        # Caching is tricky when `id` is the same for different AST instances:
        #
        # Tricky case #1: If a function body with the same `id` was cached
        # previously, we need to use that older instance of the body AST as the
        # canonical one rather than `body`. This ensures the lifetime of the
        # body in the cache will always cover the lifetime of the parent
        # `RuntimeGeneratedFunction`s when they share the same `id`.
        #
        # Tricky case #2: Unless we hold a separate reference to
        # `cache[id].value`, the GC can collect it (causing it to become
        # `nothing`). So root it in a local variable first.
        #
        cached_body = haskey(cache, id) ? cache[id].value : nothing
        cached_body = cached_body !== nothing ? cached_body : body
        # Use a WeakRef to allow `body` to be garbage collected. (After GC, the
        # cache will still contain an empty entry with key `id`.)
        cache[id] = WeakRef(cached_body)
        return cached_body
    end
end

function _lookup_body(cache_tag, id)
    lock(_cache_lock) do
        cache = getfield(parentmodule(cache_tag), _cachename)
        cache[id].value
    end
end

"""
    RuntimeGeneratedFunctions.init(mod)

Use this at top level to set up your module `mod` before using
`@RuntimeGeneratedFunction`.
"""
function init(mod)
    lock(_cache_lock) do
        if !isdefined(mod, _cachename)
            mod.eval(quote
                const $_cachename = Dict()
                struct $_tagname
                end

                # We create method of `generated_callfunc` in the user's module
                # so that any global symbols within the body will be looked up
                # in the user's module scope.
                #
                # This is straightforward but clunky.  A neater solution should
                # be to explicitly expand in the user's module and return a
                # CodeInfo from `generated_callfunc`, but it seems we'd need
                # `jl_expand_and_resolve` which doesn't exist until Julia 1.3
                # or so. See:
                #   https://github.com/JuliaLang/julia/pull/32902
                #   https://github.com/NHDaly/StagedFunctions.jl/blob/master/src/StagedFunctions.jl#L30
                @inline @generated function $RuntimeGeneratedFunctions.generated_callfunc(
                        f::$RuntimeGeneratedFunctions.RuntimeGeneratedFunction{argnames, cache_tag, $_tagname, id}, __args...) where {argnames, cache_tag, id}
                    $RuntimeGeneratedFunctions.generated_callfunc_body(argnames, cache_tag, id, __args)
                end
            end)
        end
    end
end

###
### Utilities
###
normalize_args(args::Vector) = map(normalize_args, args)
normalize_args(arg::Symbol) = arg
function normalize_args(arg::Expr)
    arg.head === :(::) || error("argument malformed. Got $arg")
    arg.args[1]
end

function expr_to_id(ex)
    io = IOBuffer()
    Serialization.serialize(io, ex)
    return Tuple(reinterpret(UInt32, sha1(take!(io))))
end

@nospecialize

closures_to_opaque(x, _=nothing) = x
_tconvert(T, x) = Expr(:(::), Expr(:call, GlobalRef(Base, :convert), T, x), T)
function closures_to_opaque(ex::Expr, return_type=nothing)
    head, args = ex.head, ex.args
    fdef = splitdef(ex; throw=false)
    if fdef !== nothing
        body = get(fdef, :body, nothing)
        if haskey(fdef, :rtype)
            body = _tconvert(fdef[:rtype], closures_to_opaque(body, fdef[:rtype]))
            delete!(fdef, :rtype)
        else
            body = closures_to_opaque(body)
        end
        fdef[:head] = :(->)
        fdef[:body] = body
        name = get(fdef, :name, nothing)
        name !== nothing && delete!(fdef, :name)
        _ex = Expr(:opaque_closure, combinedef(fdef))
        # TODO: emit named opaque closure for better stacktraces
        # (ref https://github.com/JuliaLang/julia/pull/40242)
        if name !== nothing
            name isa Symbol ||
                error("Unsupported function definition `$ex` in RuntimeGeneratedFunction.")
            _ex = Expr(:(=), name, _ex)
        end
        return _ex
    elseif head === :generator
        f_args = Expr(:tuple, Any[x.args[1] for x in args[2:end]]...)
        iters = Any[x.args[2] for x in args[2:end]]
        return Expr(
            :call,
            GlobalRef(Base, :Generator),
            closures_to_opaque(Expr(:(->), f_args, args[1])),
            iters...,
        )
    elseif head === :opaque_closure
        return closures_to_opaque(args[1])
    elseif head === :return && return_type !== nothing
        return Expr(:return, _tconvert(return_type, closures_to_opaque(args[1], return_type)))
    end
    return Expr(head, Any[closures_to_opaque(x, return_type) for x in args]...)
end

# We write an explicit deserialize() here to trigger caching of the body on a
# remote node when using Serialialization.jl (in Distributed.jl and elsewhere)
function Serialization.deserialize(s::AbstractSerializer,
        ::Type{RuntimeGeneratedFunction{argnames, cache_tag, context_tag, id}}) where {argnames,cache_tag,context_tag,id}
    body = deserialize(s)
    cached_body = _cache_body(cache_tag, id, body)
    RuntimeGeneratedFunction{argnames, cache_tag, context_tag, id}(cached_body)
end

@specialize

end


text:
function ComplexFormat()
    return ComplexFormat(())
end

function ComplexFormat(arg0::JString)
    return ComplexFormat((JString,), arg0)
end

function ComplexFormat(arg0::JString, arg1::NumberFormat)
    return ComplexFormat((JString, NumberFormat), arg0, arg1)
end

function ComplexFormat(arg0::JString, arg1::NumberFormat, arg2::NumberFormat)
    return ComplexFormat((JString, NumberFormat, NumberFormat), arg0, arg1, arg2)
end

function ComplexFormat(arg0::NumberFormat)
    return ComplexFormat((NumberFormat,), arg0)
end

function ComplexFormat(arg0::NumberFormat, arg1::NumberFormat)
    return ComplexFormat((NumberFormat, NumberFormat), arg0, arg1)
end

function format(obj::ComplexFormat, arg0::Double)
    return jcall(obj, "format", JString, (Double,), arg0)
end

function format(obj::ComplexFormat, arg0::JComplex)
    return jcall(obj, "format", JString, (JComplex,), arg0)
end

function format(obj::ComplexFormat, arg0::JComplex, arg1::StringBuffer, arg2::FieldPosition)
    return jcall(obj, "format", StringBuffer, (JComplex, StringBuffer, FieldPosition), arg0, arg1, arg2)
end

function format(obj::ComplexFormat, arg0::Object, arg1::StringBuffer, arg2::FieldPosition)
    return jcall(obj, "format", StringBuffer, (Object, StringBuffer, FieldPosition), arg0, arg1, arg2)
end

function get_available_locales(::Type{ComplexFormat})
    return jcall(ComplexFormat, "getAvailableLocales", Vector{Locale}, ())
end

function get_complex_format(::Type{ComplexFormat})
    return jcall(ComplexFormat, "getComplexFormat", ComplexFormat, ())
end

function get_complex_format(::Type{ComplexFormat}, arg0::JString, arg1::Locale)
    return jcall(ComplexFormat, "getComplexFormat", ComplexFormat, (JString, Locale), arg0, arg1)
end

function get_complex_format(::Type{ComplexFormat}, arg0::Locale)
    return jcall(ComplexFormat, "getComplexFormat", ComplexFormat, (Locale,), arg0)
end

function get_imaginary_character(obj::ComplexFormat)
    return jcall(obj, "getImaginaryCharacter", JString, ())
end

function get_imaginary_format(obj::ComplexFormat)
    return jcall(obj, "getImaginaryFormat", NumberFormat, ())
end

function get_real_format(obj::ComplexFormat)
    return jcall(obj, "getRealFormat", NumberFormat, ())
end

function parse(obj::ComplexFormat, arg0::JString)
    return jcall(obj, "parse", JComplex, (JString,), arg0)
end

function parse(obj::ComplexFormat, arg0::JString, arg1::ParsePosition)
    return jcall(obj, "parse", JComplex, (JString, ParsePosition), arg0, arg1)
end


text:
struct VLFResult
    # fit::Optim.MultivariateOptimizationResults
    fit::LsqFit.LsqFitResult
    param_array::AbstractArray
    data::VirusLoadData
    names::Vector
end

# H(x::AbstractFloat) = ifelse(x < 0, zero(x), ifelse(x > 0, one(x), oftype(x,0.5)))
H(x) = 0.5*(tanh(1e3*x) + 1.)

v₁(t, a₁, a₂, logVmax) = 1. + (10^logVmax - 1.0)*(tanh(6.0*(t - (a₁ + a₂)/2)/abs(a₂ - a₁)) - tanh(-3.0*(a₂ + a₁)/abs(a₂ - a₁)))/2

v₂(t, a₂, α) = 1. - H(t - a₂) + H(t - a₂)*exp(-α*(t - a₂))

v₃(t, b₁, b₂, logVmin) = 1. - (1.0 - 10^logVmin)*(tanh(6.0*(t - (b₁ + b₂)/2)/abs(b₂ - b₁)) - tanh(-3.0*(b₂ + b₁)/abs(b₂ - b₁)))/2

function assert_params(a₁, a₂, b₁, b₂, α, logVmax)
    all([a₁ < a₂, a₂ < b₁, b₁ < b₂, a₁ > 0, a₂ > 0, b₁ > 0, b₂ > 0, α > 0, logVmax > 0])
end

function LogVirusLoadFunction(t, a₁, a₂, b₁, b₂, α, logVmax)
    logVmin = -7.0
    log10(v₁(t, a₁, a₂, logVmax)*v₂(t, a₂, α)*v₃(t, b₁, b₂, logVmin))
end

function LogEmpiricalVirusLoadFunction(t::Vector, p::Vector, data::VirusLoadData)
    θ = minimum(data.v)
    logV = LogVirusLoadFunction.(t, p...)
    max.(logV, θ)
end

function rss(p, data)
    sum(abs2, LogEmpiricalVirusLoadFunction(data.t, p, data) - data.v)
end

function get_bounds(data, p0)
    a₁, a₂, b₁, b₂, α, logVmax = p0
    a₁a₂ = (a₁ + a₂)/2
    a₂b₁ = (a₂ + b₁)/2
    b₁b₂ = (b₁ + b₂)/2
    tmin = 0.0
    tmax = maximum(data.t)
    vmax = maximum(data.v)
    lb   = [tmin, a₁a₂, a₂b₁, b₁b₂, 1e-8, vmax-3.0]
    ub   = [a₁a₂, a₂b₁, b₁b₂, tmax, 1e+4, vmax+3.0]
    return lb, ub
end

function get_rand_pars(fit1, data)
    tmax = maximum(data.t)
    p = fit1.param
    s = 1.1
    a₁ = rand(truncated(Normal(p[1], s), 0, tmax))
    a₂ = rand(truncated(Normal(p[2], s), a₁, tmax))
    b₁ = rand(truncated(Normal(p[3], s), a₂, tmax))
    b₂ = rand(truncated(Normal(p[4], s), b₁, tmax))
    α = rand(truncated(Normal(p[5], 0.5), 0, Inf))
    logVmax = rand(truncated(Normal(p[6], 0.5), 0, Inf))
    return [a₁, a₂, b₁, b₂, α, logVmax]
end

function fitVLF(data::VirusLoadData, p0::Vector)
    lb, ub = get_bounds(data, p0)
    return curve_fit((t, p) -> LogEmpiricalVirusLoadFunction(t, p, data), data.t, data.v, p0, lower=lb, upper=ub)
    # cost(p) = sum(abs2, LogEmpiricalVirusLoadFunction(data.t, p, data) - data.v)
    # fit = optimize(cost, lb, ub, p0, SAMIN(), Optim.Options(iterations=10^4))
end

function fitVLF(data::VirusLoadData; c=0.15, niter = 10^3)
    # Best parameter
    ϵ = -2/length(data.t)*log(c)
    fit1 = nothing
    pb1 = Progress(niter, 0.5, "Fitting the VLF to data ")
    for iter = 1:niter
        #      a₁,     a₂,     b₁,     b₂      α,      logVmax
        par0 = vcat(sort(data.t[end]*rand(4)), rand(), maximum(data.v)*(rand()+0.5))
        fit0 = try
            VirusLoadCurve.fitVLF(data, par0)
        catch y
            nothing
        end
        if fit0 !== nothing
            if fit1 === nothing
                fit1 = fit0
            else
                if sum(abs2, fit0.resid) < sum(abs2, fit1.resid)
                    fit1 = fit0
                end
            end
        end
        next!(pb1)
    end
    # Set of parameters
    param_array = [fit1.param]
    niter = 10^6
    tol = (1.0 + ϵ)*rss(fit1.param, data)
    pb2 = Progress(niter, 0.5, "Finding possible parameter values ")
    for k = 1:niter
        par = get_rand_pars(fit1, data)
        if (rss(par, data) < tol && assert_params(par...))
            push!(param_array, par)
        end
        next!(pb2)
    end
    println("Number of possible parameters found: ", length(param_array), " (out of ", niter, " testted)")
    names = ["a₁", "a₂", "b₁", "b₂", "α", "logVmax"]
    VLFResult(fit1, param_array, data, names)
end

function param_extrema(param_array)
    if length(param_array)==0 return zip(fill(-Inf, 6), fill(Inf, 6)) end
    array = map(x -> map(y -> y[x], param_array), 1:length(param_array[1]))
    param_lower = map(minimum, array)
    param_upper = map(maximum, array)
    return zip(param_lower, param_upper)
end

@recipe function f(result::VLFResult; empirical=false, plotrange=true)
    tmin, tmax = extrema(result.data.t)
    tt = Vector(range(tmin, tmax, step=1e-3))
    func = empirical ? ((t, p) -> LogEmpiricalVirusLoadFunction(t, p, result.data)) : ((t, p) -> LogVirusLoadFunction.(t, p...))
    yy = func(tt, result.fit.param)
    vmin, vmax = extrema(result.data.v)
    xaxis := ("Time (days)", (tmin, tmax), font(14))
    yaxis := (L"\log\,V(t)", font(14))
    grid := :none

    @series begin
        x := tt
        y := yy
        linecolor := :black
        linewidth := 4
        label := ("Virus load function")
        ()
    end

    if plotrange
        yylower = yyupper = yy
        pmin = min(length(result.param_array), 1000)
        for param in result.param_array[1:pmin]
            yy0 = func(tt, param)
            yylower = min.(yylower, yy0)
            yyupper = max.(yyupper, yy0)
        end
        vmax = max(vmax, maximum(yyupper))
        ylims := (vmin-0.5, vmax+0.5)

        @series begin
            x := tt
            y := yylower
            linewidth := 0
            linecolor := 2
            linealpha := 0.4
            fillcolor := 2
            fillalpha := 0.4
            fillrange := yyupper
            label := ("Range of VLF")
            ()
        end
    end

end

function Base.summary(result::VLFResult)
    println(@sprintf "RSS = %.5e" sum(abs2, result.fit.resid))
    [println(p, " = ", v, ", CI=", ci) for (p, v, ci) in zip(result.names, result.fit.param, param_extrema(result.param_array))]
    nothing
end


text:
"""
Paper: Bayesian inference for finite mixtures of univariate and multivariate
       skew-normal and skew-t distributions, Biostatistics 2010.
skew (delta): a real number in (-1, 1)
"""
function rand_skewnormal(loc, scale, skew)
    z = rand(TruncatedNormal(0, 1, 0, Inf))
    return loc + scale * skew * z + scale * sqrt(1 - skew ^ 2) * randn()
end

function rand_skewnormal(loc, scale, skew, dims...)
    z = rand(TruncatedNormal(0, 1, 0, Inf), dims...)
    return loc .+ scale * skew * z + scale * sqrt(1 - skew ^ 2) * randn(dims...)
end


text:
## Expected periodogram ##

### Multivariate ###

"""
    EI!(store, model::TimeSeriesModel)

Compute the expected periodogram and assign to appropriate place in memory (note this computes at twice the desired resolution).
"""
function EI!(store::Union{TimeSeriesModelStorageFunction,AdditiveStorage}, model::TimeSeriesModel)
    acv!(store, model)
    _func_EI!(store)
    return nothing
end
_func_EI!(store::AdditiveStorage) = _func_EI!(store.store1)
_func_EI!(store::TimeSeriesModelStorageFunction) = _EI!(store.funcmemory, store.encodedtime.Δ)

"""
    grad_EI!(store::TimeSeriesModelStorageGradient, model::TimeSeriesModel)

Compute the gradient of the expected periodogram and assign to appropriate place in memory (note this computes at twice the desired resolution).
"""
function grad_EI!(store::TimeSeriesModelStorageGradient, model::TimeSeriesModel)
    grad_acv!(store, model)
    _EI!(store.gradmemory, store.encodedtime.Δ)
    return nothing
end
function grad_EI!(store::AdditiveStorage, model::AdditiveTimeSeriesModel)
    grad_EI!(store.store1, model.model1)
    grad_EI!(store.store2, model.model2)
    return nothing
end

"""
    hess_EI!(store::TimeSeriesModelStorageHessian, model::TimeSeriesModel)

Compute the Hessian of the expected periodogram and assign to appropriate place in memory (note this computes at twice the desired resolution).
"""
function hess_EI!(store::TimeSeriesModelStorageHessian, model::TimeSeriesModel)
    hess_acv!(store, model)
    _EI!(store.hessmemory, store.encodedtime.Δ)
    return nothing
end
function hess_EI!(store::AdditiveStorage, model::AdditiveTimeSeriesModel)
    @views hess_EI!(store.store1, model.model1)
    @views hess_EI!(store.store2, model.model2)
    return nothing
end

_EI!(store::Sdf2EIStorage, Δ) = __EI!(store.acv2EI, Δ)
_EI!(store::Acv2EIStorage, Δ) = __EI!(store, Δ)
_EI!(store::Sdf2EIStorageUni, Δ) = __EI!(store.acv2EI, Δ)
_EI!(store::Acv2EIStorageUni, Δ) = __EI!(store, Δ)

function __EI!(store, Δ)
    mult_allocated_kernel!(store.allocatedarray, store.kernel)
    store.planned_fft * store.allocatedarray # essentially fft!(acv, ndims(store.allocatedarray))
    store.allocatedarray .*= Δ / (2π)
    return nothing
end

function mult_allocated_kernel!(array::Vector{ComplexF64},kernel)
    array .*= kernel
    nothing
end
function mult_allocated_kernel!(array::Matrix{ComplexF64},kernel)
    size(array, 2) == length(kernel) || throw(DimensionMismatch("size of second dimension of array should equal length of kernel."))
    @inbounds for i ∈ 1:size(array, 2)
        @views array[:, i] .*= kernel[i]
    end
    nothing
end
function mult_allocated_kernel!(array::Array{ComplexF64,3},kernel)
    size(array, 3) == length(kernel) || throw(DimensionMismatch("size of third dimension of array should equal length of kernel."))
    @inbounds for i ∈ 1:size(array, 3)
        @views array[:, :, i] .*= kernel[i]
    end
    nothing
end

"""
    extract_EI(store::Sdf2EIStorage)
    extract_EI(store::Acv2EIStorage)
    extract_EI(store::Sdf2EIStorageUni)
    extract_EI(store::Acv2EIStorageUni)

Extract the EI from general storage.
"""
extract_EI(store) = extract_EI(extract_S(store))
extract_EI(store::Sdf2EIStorage) = extract_EI(store.acv2EI)
extract_EI(store::Acv2EIStorage) = store.hermitianarray
extract_EI(store::Sdf2EIStorageUni) = extract_EI(store.acv2EI)
extract_EI(store::Acv2EIStorageUni) = store.allocatedarray

"""
    EI(model::TimeSeriesModel, n, Δ)

Compute EI at Fourier frequencies `fftshift(fftfreq(n,2π/Δ))`.

Note internal computation provides values at twice the resolution, this function returns at the desired resolution.
"""
function EI(model::TimeSeriesModel, n, Δ)
    store = allocate_memory_EI_F(typeof(model), n , Δ)
    EI!(store, model)
    return fftshift(copy(extract_EI(store)))[1:2:end]
end


text:
include("04_one_site_operator.jl")

println("
#######################################################
# Tutorial 5                                          #
#                                                     #
# 1-site custom operator/gate definition              #
#######################################################
")

# Define your own `op`:
import ITensors: op

function op(::OpName"iX", ::SiteType"S=1/2", i::Index)
  iX = ITensor(i', dag(i))
  iX[i' => 2, i => 1] = im
  iX[i' => 1, i => 2] = im
  return iX
end

# Alternative shorthand with matrices:
function op(::OpName"iX", ::SiteType"S=1/2")
  return [
    0 im
    im 0
  ]
end

iX = ITensor(i', dag(i))
iX[i' => 2, i => 1] = im
iX[i' => 1, i => 2] = im

@show iX == op("iX", i)

@show inner(Xp', iX, Xp)
@show inner(Xp, apply(iX, Xp))


text:
@testset "Block Tridiagonal" begin
    solve = AirfoilMeshers.solve_block_tridiag!
    DIM = 5
    A = [rand(2, 2) for i in 1:DIM]
    B = [I + rand(2, 2) for i in 1:DIM]
    C = [rand(2, 2) for i in 1:DIM]
    D = [rand(2) for i in 1:DIM]
    MM = zeros(2DIM, 2DIM)
    blk(MM, i, j) = @view MM[(2i - 1):(2i), (2j - 1):(2j)]
    for i in 1:DIM
        blk(MM, i, i) .= B[i]
        if i > 1
            blk(MM, i, i - 1) .= A[i]
        end
        if i < DIM
            blk(MM, i, i + 1) .= C[i]
        end
    end
    DD = vcat(D...)
    X0 = MM \ DD
    solve(A,B,C,D)
    X1 = vcat(D...)
    @test X0 ≈ X1
end


text:
using Mmap

@testset "datafreetree" begin
    function test(data, data2, data3)
        t = DataFreeTree(KDTree, data)
        @test_throws ArgumentError injectdata(t, data2)
        @test_throws DimensionMismatch injectdata(t, data3)
        for typ in [KDTree, BallTree]
            dfilename = tempname()
            d = 2
            n = 100
            mktemp() do _, io
                data = Mmap.mmap(io, Matrix{Float32}, (d, n))
                data[:] = rand(Float32, d, n)
                t = injectdata(DataFreeTree(typ, data), data)
                tr = typ(data)
                for i = 1:n
                    @test knn(t, data[:,i], 3) == knn(tr, data[:,i], 3)
                end
                finalize(data)
            end
        end
    end
    data = rand(2,100)
    data2 = rand(2,100)
    data3 = rand(3,100)
    test(data, data2, data3)
    test(view(data, :, :), view(data2, :, :), view(data3, :, :))
end


text:
module MolSimReg

import Registrator

const _registered_modules =
   [
      "JuLIP",
      "SHIPs",
      "SKTB",
      "JuLIPMaterials",
      "NBodyIPs",
      "IPFitting",
      "MolSimReg",
      "PolyPairPots",
      "Isaac"
   ]

regpath() = joinpath(homedir(), ".julia", "registries", "MolSim")

"""
register the module `M` in the `MolSim` registry
```
register(M::Module; force = false, push = false)
```
"""
function register(M::Module; force = false, push = false)
   if !(string(M) in _registered_modules) && !force
      error("The module $(string(M)) appears not to be registered yet. Use `force = true`")
   end
   Registrator.register(M, regpath())
   if push
      push!()
   end
end

"""
push the registry to the remote git repository
"""
function push!()
   error("push! is not yet working")
   run(`cd $(regpath())`)
   run(`git push`)
   run(`cd -`)
end

end # module