text:
"""
    innereltype(x)

Recursively determine the 'innermost' type in by the collection `x` (which may be, for example,
a collection of a collection).
"""
function innereltype(x)
  T = eltype(x)
  T <: AbstractArray ? innereltype(T) : return T
end

"""
    cxtype(T)

Returns `T` when `T` is `Complex`, or `Complex{T}` when `T` is `Real`.
"""
cxtype(::Type{T}) where T<:Number = T
cxtype(::Type{T}) where T<:Real = Complex{T}

"""
    fltype(T)

Returns `T` when `T<:AbstractFloat` or `Tf` when `T<:Complex{Tf}`.
"""
fltype(::Type{T})          where T<:AbstractFloat = T
fltype(::Type{Complex{T}}) where T<:AbstractFloat = T
fltype(T::Tuple) = fltype(T[1])

cxeltype(x) = cxtype(innereltype(x))
fleltype(x) = fltype(innereltype(x))

"""
    superzeros(T, A)

Returns an array like `A`, but full of zeros. If `innereltype(A)` can be promoted to `T`, then
the innermost elements of the array will have type `T`.
"""
superzeros(T, A::AbstractArray) = T(0)*A
superzeros(A::AbstractArray) = superzeros(innereltype(A), A)
superzeros(T, dims::Tuple) = eltype(dims) <: Tuple ? [ superzeros(T, d) for d in dims ] : zeros(T, dims)
superzeros(dims::Tuple) = superzeros(Float64, dims) # default
superzeros(T::Tuple, dims::Tuple) = [ superzeros(T[i], dims[i]) for i=1:length(dims) ]

"""
    @superzeros T a b c d...
    @superzeros T dims b c d...

Generate arrays `b, c, d...` with the super-dimensions of `a` and innereltype `T`.
"""
macro superzeros(T, ad, vars...)
  expr = Expr(:block)
  append!(expr.args, [:( $(esc(var)) = superzeros($(esc(T)), $(esc(ad))); ) for var in vars])
  expr
end
supersize(a) = Tuple([size(ai) for ai in a])
supersize(a::Array{T}) where T<:AbstractArray = Tuple([size(ai) for ai in a])
supersize(a::Array{T}) where T<:Number = size(a)

macro createarrays(T, dims, vars...)
  expr = Expr(:block)
  append!(expr.args, [:($(esc(var)) = zeros($(esc(T)), $(esc(dims))); ) for var in vars])
  expr
end

"""
    @zeros T dims a b c...

Create arrays of all zeros with element type `T`, size `dims`, and global names
`a`, `b`, `c` (for example). An arbitrary number of arrays may be created.
"""
macro zeros(T, dims, vars...)
  expr = Expr(:block)
  append!(expr.args, [:($(esc(var)) = zeros($(esc(T)), $(esc(dims))); ) for var in vars])
  expr
end

Base.zeros(::CPU, T, dims) = zeros(T, dims)

"""
    devzeros(dev, T, dims)

Returns an array like `A` of type `T`, but full of zeros.
"""
devzeros(dev, T, dims) = zeros(dev, T, dims)

"""
    @devzeros dev T dims a b c...

Create arrays of all zeros with element type `T`, size `dims`, and global names
`a`, `b`, `c` (for example) on device `dev`.
"""
macro devzeros(dev, T, dims, vars...)
  expr = Expr(:block)
  append!(expr.args, [:($(esc(var)) = zeros($(esc(dev))(), $(esc(T)), $(esc(dims))); ) for var in vars])
  expr
end

"""
    varsexpression(name, fieldspecs; parent=:AbstractVars, typeparams=nothing)

    varsexpression(name, fieldspecs...; parent=:AbstractVars, typeparams=nothing)

    varsexpression(name, physicalfields, fourierfields; physicaltype=:Tp, fouriertype=:Tf,
                   parent=:AbstractVars, typeparams=[physicaltype, fouriertype])

Returns an expression that defines an `AbstractVars` type.
"""
function varsexpression(name, fieldspecs; parent=:AbstractVars,
                        typeparams::Union{Nothing,Symbol,Array{Symbol,1}}=nothing)
  if typeparams == nothing
    signature = name
  else
    try
      signature = Expr(:curly, name, typeparams...)
    catch
      signature = Expr(:curly, name, typeparams) # only one typeparam given?
    end
  end

  fieldexprs = [ :( $(spec[1])::$(spec[2]); ) for spec in fieldspecs ]

  :(struct $signature <: $parent; $(fieldexprs...); end)
end

function varsexpression(name, physicalfields, fourierfields; parent=:AbstractVars,  physicaltype=:Tp, fouriertype=:Tf,
                        typeparams=[physicaltype, fouriertype])
  physicalfieldspecs = getfieldspecs(physicalfields, physicaltype)
   fourierfieldspecs = getfieldspecs(fourierfields, fouriertype)
  varsexpression(name, cat(physicalfieldspecs, fourierfieldspecs; dims=1); parent=parent, typeparams=typeparams)
end

"""
    getfieldspecs(fieldnames, fieldtype)

Returns an array of (fieldname[i], fieldtype) tuples.
"""
getfieldspecs(fieldnames::AbstractArray, fieldtype) = [ (name, fieldtype) for name in fieldnames ]

"""
    parsevalsum2(uh, g)

Returns `∫|u|² = Σ|uh|²` on the grid `g`, where `uh` is the Fourier transform of `u`.
"""
function parsevalsum2(uh, g::TwoDGrid)
  if size(uh, 1) == g.nkr # uh is in conjugate symmetric form
    U = @views sum(abs2, uh[1, :])           # k=0 modes
    U += @views 2*sum(abs2, uh[2:end, :])    # sum k>0 modes twice
  else # count every mode once
    U = sum(abs2, uh)
  end
  norm = g.Lx*g.Ly/(g.nx^2*g.ny^2)    # normalization for dft
  norm*U
end

function parsevalsum2(uh, g::OneDGrid)
  if size(uh, 1) == g.nkr # uh is conjugate symmetric
    U = sum(abs2, uh[1])                  # k=0 modes
    U += @views 2*sum(abs2, uh[2:end])    # sum k>0 modes twice
  else # count every mode once
    U = sum(abs2, uh)
  end
  norm = g.Lx/g.nx^2 # normalization for dft
  norm*U
end

"""
    parsevalsum(uh, g)

Returns `real(Σ uh)` on the grid `g`.
"""
function parsevalsum(uh, g::TwoDGrid)
  if size(uh, 1) == g.nkr       # uh is conjugate symmetric
    U = sum(uh[1, :])           # k=0 modes
    U += 2*sum(uh[2:end, :])    # sum k>0 modes twice
  else # count every mode once
    U = sum(uh)
  end

  norm = g.Lx*g.Ly/(g.nx^2*g.ny^2) # weird normalization for dft
  norm*real(U)
end

"""
    jacobianh(a, b, g)

Returns the transform of the Jacobian of two fields a, b on the grid g.
"""
function jacobianh(a, b, g::TwoDGrid)
  if eltype(a) <: Real
    bh = rfft(b)
    bx = irfft(im*g.kr.*bh, g.nx)
    by = irfft(im*g.l.*bh, g.nx)
    return im*g.kr.*rfft(a.*by)-im*g.l.*rfft(a.*bx)
  else
    # J(a, b) = dx(a b_y) - dy(a b_x)
    bh = fft(b)
    bx = ifft(im*g.k.*bh)
    by = ifft(im*g.l.*bh)
    return im*g.k.*fft(a.*by).-im*g.l.*fft(a.*bx)
  end
end

"""
    jacobian(a, b, g)

Returns the Jacobian of a and b.
"""
function jacobian(a, b, g::TwoDGrid)
  if eltype(a) <: Real
   return irfft(jacobianh(a, b, g), g.nx)
  else
   return ifft(jacobianh(a, b, g))
  end
end

"""
    radialspectrum(ah, g; n=nothing, m=nothing, refinement=2)

Returns `aρ = ∫ ah(ρ,θ) ρ dρ dθ`, the radial spectrum of `ah` known on the
Cartesian wavenumber grid (k,l).

`aρ` is found by intepolating `ah` onto a polar wavenumber grid (ρ,θ), and
then integrating over `θ` to find `aρ`. The default resolution (n,m) for the
polar wave number grid is `n=refinement*maximum(nk, nl),
m=refinement*maximum(nk, nl)`, where `refinement=2` by default. If
`ah` is in conjugate symmetric form only the upper half plane in `θ` is
represented on the polar grid.
"""
function radialspectrum(ah, g::TwoDGrid; n=nothing, m=nothing, refinement=2)

  n = n == nothing ? refinement*maximum([g.nk, g.nl]) : n
  m = m == nothing ? refinement*maximum([g.nk, g.nl]) : m

  # Calcualte shifted k and l
  lshift = range(-g.nl/2+1, stop=g.nl/2, length=g.nl)*2π/g.Ly

  if size(ah)[1] == g.nkr # conjugate symmetric form
    m = Int(m/2)                         # => half resolution in θ
    θ = range(-π/2, stop=π/2, length=m)  # θ-grid from k=0 to max(kr)
    ahshift = fftshift(ah, 2)            # shifted ah
    kshift = range(0, stop=g.nkr-1, length=g.nkr)*2π/g.Lx
  else # ordinary form
    θ = range(0, stop=2π, length=m)      # θ grid
    ahshift = fftshift(ah, [1, 2])       # shifted ah
    kshift = range(-g.nk/2+1, stop=g.nk/2, length=g.nk)*2π/g.Lx
  end

  # Interpolator for ah
  itp = scale(interpolate(ahshift, BSpline(Linear())), kshift, lshift)

  # Get radial wavenumber vector
  ρmax = minimum([(g.nk/2-1)*2π/g.Lx, (g.nl/2-1)*2π/g.Ly])
  ρ = range(0, stop=ρmax, length=n)

  # Interpolate ah onto fine grid in (ρ,θ).
  ahρθ = zeros(eltype(ahshift), (n, m))

  for i=2:n, j=1:m # ignore zeroth mode
    kk = ρ[i]*cos(θ[j])
    ll = ρ[i]*sin(θ[j])
    ahρθ[i, j] = itp(kk, ll)
  end

  # ahρ = ρ ∫ ah(ρ,θ) dθ  =>  Ah = ∫ ahρ dρ = ∫∫ ah dk dl
  dθ = θ[2]-θ[1]
  if size(ah)[1] == g.nkr
    ahρ = 2ρ.*sum(ahρθ, dims=2)*dθ # multiply by 2 for conjugate symmetry
  else
    ahρ = ρ.*sum(ahρθ, dims=2)*dθ
  end

  ahρ[1] = ah[1, 1] # zeroth mode

  ρ, ahρ
end

# Moments and cumulants
"Compute the average of `c` on the grid `g`."
domainaverage(c, g) = g.dx*g.dy*sum(c)/(g.Lx*g.Ly)

"Compute the `n`th x-moment of `c` on the grid `g`."
xmoment(c, g, n=1) = sum(g.X.^n.*c)/sum(c)

"Compute the `n`th y-moment of `c` on the grid `g`."
ymoment(c, g, n=1) = sum(g.Y.^n.*c)/sum(c)

ArrayType(::CPU) = Array
ArrayType(::CPU, T, dim) = Array{T, dim}


text:
mutable struct DataPoint{T,E}
    @atomic value::T
    start_time_unix_nano::UInt
    time_unix_nano::UInt
    exemplar_reservoir::E
end

function DataPoint{T}(exemplar_reservoir = nothing) where {T}
    t = UInt(time() * 10^9)
    DataPoint(zero(T), t, t, exemplar_reservoir)
end

function _add_to_datapoint!(p::DataPoint, v, t)
    @atomic p.value += v
    p.time_unix_nano = t
end

function _update_to_datapoint!(p::DataPoint, v, t)
    @atomic p.value = v
    p.time_unix_nano = t
end

ext: jl


text:
using Debugger: DebuggerState

copy_stack(::Nothing) = nothing

copy_stack(fd::FrameData) =
  FrameData(copy(fd.locals), copy(fd.ssavalues), copy(fd.sparams),
  fd.exception_frames, fd.last_exception, fd.last_reference,
  fd.callargs)

copy_stack(st::Frame) =
  Frame(st.framecode, copy_stack(st.framedata),
  st.pc, copy_stack(st.caller), st.callee)

copy_stack(st::DebuggerState) =
  DebuggerState(copy_stack(st.frame), st.level, st.broke_on_error,
  st.watch_list, st.lowered_status, st.mode,
  st.repl, st.terminal, st.main_mode,
  st.julia_prompt, st.standard_keymap, st.overall_result)

function reset_(f)
  try f()
  catch e
    e isa InterpreterError || rethrow()
    p = e.err
    p isa Tuple{Any,Continuation} || rethrow()
    return p[1](p[2])
  end
end

struct Continuation
  st::DebuggerState
end

Base.show(io::IO, ::Continuation) = print(io, "Continuation()")

function (c::Continuation)(x = nothing)
  st = copy_stack(c.st)
  provide_result!(st, x)
  inc_pc!(st)
  reset_(() -> runall(Continuations(), st))
end

struct Continuations end

function reset(f)
  reset_(() -> overdub(Continuations(), f))
end

@primitive Continuations reset(f) = overdub(Continuations(), f)

shift(f) = error("`shift` only works inside `reset`")

isprimitive(::Continuations, ::typeof(shift), f) = true

primitive_(::Continuations, state, ::typeof(shift), f) =
  throw((f, Continuation(state)))

macro reset(ex)
  :(reset(() -> $(esc(ex))))
end


text:
@testset "test_matrix_comps" begin
A = [-0.21 0.2; 0.2 -0.21]
B = 0.01*[1 0; 0 1]
C = [1 0; 0 1]
D = 0
sys = ss(A,B,C,D)
sysr, G = balreal(sys)

@test gram(sysr, :c) ≈ G
@test gram(sysr, :o) ≈ G
@test sort(pole(sysr)) ≈ sort(pole(sys))

sysb,T = ControlSystems.balance_statespace(sys)
Ab,Bb,Cb,T = ControlSystems.balance_statespace(A,B,C)

@test Ab*T ≈ T*A
@test Bb ≈ T*B
@test Cb*T ≈ C

@test sysb.A ≈ Ab

@test ControlSystems.balance_transform(A,B,C) ≈ ControlSystems.balance_transform(sys)

W = [1 0; 0 1]
@test covar(sys, W) ≈ [0.002560975609756 0.002439024390244; 0.002439024390244 0.002560975609756]
D2 = [1 0; 0 1]
@test covar(ss(A,B,C,D2, 1), W) ≈ [1.000110108378310 -0.000010098377310; -0.000010098377310 1.000110108378310]
# Direct term means infinite covariance
@test covar(ss(A,B,C,D2), W) ≈ [Inf Inf; Inf Inf]

# No noise on second output should give finite variance
@test covar(ss(A,B,C,[1 0; 0 0]), W) ≈ [Inf Inf; Inf 0.002560975609756]

# Unstable system has inf covar
@test covar(ss([1 0; 0 1],B,C,0), W) == [Inf Inf; Inf Inf]

# Discrete system can have direct term
@test covar(ss(A,B,C,D2,0.1),W) ≈ [1.00011010837831 -1.0098377309782909e-5; -1.0098377309782909e-5 1.00011010837831]

# TODO test in Julia 0.7 to see if supported
# # Test special matrices
As = sparse(A)
Bs = sparse(B)
Cs = sparse(C)
Asb,Bsb,Csb,Ts = ControlSystems.balance_statespace(As,Bs,Cs) #Error no LAPACK function
#
# @test Abs*Ts ≈ Ts*As
# @test Bbs ≈ Ts*Bs
# @test Cbs*Ts ≈ Cs

# Test special values
Ar = rationalize.(A)
Br = rationalize.(B)
Cr = rationalize.(Float64.(C))    # When did ever rationalize work on Int?
Arb,Brb,Crb,Tr = ControlSystems.balance_statespace(Ar,Br,Cr)

@test Arb*Tr ≈ Tr*Ar
@test Brb ≈ Tr*Br
@test Crb*Tr ≈ Cr

Tr = randn(2,2)
syst = similarity_transform(sys, Tr)
@test sys.A*Tr ≈ Tr*syst.A
@test sys.B ≈ Tr*syst.B
@test sys.C*Tr ≈ syst.C

sys = ss([1 0.1; 0 1], ones(2), [1. 0], 0)
sysi = ControlSystems.innovation_form(sys, I, I)
@test sysi.A ≈ sysi.A
@test sysi.B ≈ [4.415675759647131
 48.334204475215365]

 sysi = innovation_form(sys)
 @test sysi.B ≈ [4.415675759647131
 48.334204475215365]

 sysi = innovation_form(sys, R2 = 2I)
 @test sysi.B ≈ [4.225661436353894
 44.52445850991302]

 sysi = innovation_form(sys, R1 = 2I)
 @test sysi.B ≈ [4.734159731874057
 54.719744515739514]

 sysi = ControlSystems.innovation_form(sys, R1=2I, R2=2I)
 @test sysi.B ≈ [4.415675759647131
 48.334204475215365]

# Test with noise filters
sysw = ss([0.5 0.1; 0 0.5], [0,1], eye_(2), 0, 1)
sysi = ControlSystems.innovation_form(sys, sysw=sysw)
@test sysi.A ≈ sys.A
@test sysi.B ≈ [4.01361818808572
 40.26132476965486]

end


text:
using Documenter, ClickTable

makedocs(;
    modules=[ClickTable],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/gabomgp4@gmail.com/ClickTable.jl/blob/{commit}{path}#L{line}",
    sitename="ClickTable.jl",
    authors="Gabriel Gómez",
    assets=String[],
)

deploydocs(;
    repo="github.com/gabomgp4@gmail.com/ClickTable.jl",
)


text:
module FTFont

using FreeType
using GRU.Math2D

mutable struct Glyph
	char::Char
	box::Rect{Int}
	origin::Vec2{Float32}
	advance::Vec2{Float32}
end

mutable struct Font
	family::AbstractString
	style::AbstractString
	size::Vec2{Float32}
	lineDistance::Float32
	ascent::Float32
	descent::Float32
	glyphs::Dict{Char, Glyph}
	kerning::Dict{Tuple{Char, Char}, Vec2{Float32}}
	bitmap::Array{UInt8, 2}
	fallbackGlyph::Glyph

	Font(family, style, sizeX, sizeY, lineDistance, ascent, descent, bmpWidth, bmpHeight) =
		new(family, style, Vec2{Float32}(sizeX, sizeY), lineDistance, ascent, descent, Dict{Char, Glyph}(), Dict{Tuple{Char, Char}, Vec2{Float32}}(), zeros(UInt8, bmpWidth, bmpHeight))
end

Font(family, style, sizeX, sizeY, lineDistance, ascent, descent, maxCharWidth, maxCharHeight, charCount) =
	Font(family, style, sizeX, sizeY, lineDistance, ascent, descent, font_bitmap_size(maxCharWidth + 1, maxCharHeight + 1, charCount)...)

function font_bitmap_size(charWidth, charHeight, charCount)
	totalPixels = charWidth * charHeight * charCount
	width = ceil(sqrt(totalPixels) / charWidth) * charWidth
	height = ceil((totalPixels / width) / charHeight) * charHeight

	@assert (width / charWidth) * (height / charHeight) >= charCount

	return Int(width), Int(height)
end

mutable struct GlyphPosition
	x::Int
	y::Int
	yMax::Int

	GlyphPosition() = new(1, 1, 1)
end

function addglyph(font::Font, char::Char, bmp::Array{UInt8, 2}, origin::Vec2{Float32}, advance::Vec2{Float32}, pos::GlyphPosition)
	@assert size(bmp, 1) < size(font.bitmap, 1)
	@assert size(bmp, 2) < size(font.bitmap, 2)
	if pos.x + size(bmp, 1) > size(font.bitmap, 1)
		pos.x = 1
		pos.y = pos.yMax + 2
	end
	@assert pos.x <= size(font.bitmap, 1)
	@assert pos.y <= size(font.bitmap, 2)

	xmax = pos.x + size(bmp, 1) - 1
	ymax = pos.y + size(bmp, 2) - 1

	glyph = Glyph(char, rect(Int, pos.x, pos.y, xmax, ymax), origin, advance)
	font.glyphs[char] = glyph
	font.bitmap[pos.x:xmax, pos.y:ymax] = bmp
	# the first added glyph is used as a fallback when a glyph is missing while drawing
	if !isdefined(font, :fallbackGlyph)
		font.fallbackGlyph = glyph
	end

	pos.x = xmax + 2
	pos.yMax = max(pos.yMax, ymax)
end

function addkerning(font::Font, c1::Char, c2::Char, distance::Vec2{Float32})
	@assert haskey(font.glyphs, c1) && haskey(font.glyphs, c2)
	@assert !haskey(font.kerning, (c1, c2))
	font.kerning[(c1, c2)] = distance
end

function adding_done(font::Font)
	# maybe calculate the real ascent and descent from the glyph data as well, since the FreeType data can be off
	corner = one(Vec2{Int})
	for g in values(font.glyphs)
		corner = max(corner, g.box.max)
	end
	corner += one(Vec2{Int})
	if corner.x < size(font.bitmap, 1) || corner.y < size(font.bitmap, 2)
		font.bitmap = font.bitmap[1:corner.x, 1:corner.y]
	end
end

mutable struct TextCursor
	pos::Vec2{Float32}
	lastChar::Char

	TextCursor(x, y) = new(Vec2{Float32}(x, y), 0)
	TextCursor(cursor::TextCursor) = new(cursor.pos, cursor.lastChar)
end

function setpos(cursor::TextCursor, x, y)
	if x != cursor.pos.x || y != cursor.pos.y
		# we reset the last drawn character when we move the cursor so kerning only takes place when we draw continuous text
		cursor.pos = Vec2{Float32}(x, y)
		cursor.lastChar = Char(0)
	end
end

function drawtext(drawRect::Function, font::Font, cursor::TextCursor, s::AbstractString)
	for c in s
		glyph = get(font.glyphs, c, font.fallbackGlyph)
		if cursor.lastChar != Char(0)
			cursor.pos += get(font.kerning, (cursor.lastChar, c), zero(Vec2{Float32}))
		end
		drawRect(cursor.pos - glyph.origin, font.bitmap, glyph.box)
		cursor.pos += glyph.advance
		cursor.lastChar = c
	end
end

function textbox(font::Font, cursor::TextCursor, s::AbstractString)
	tempCursor = TextCursor(cursor)
	boxMax = rect(Float32)
	drawtext(font, tempCursor, s) do pos, bmp, box
		boxMax.min = min(boxMax.min, pos)
		boxMax.max = max(boxMax.max, pos + size(box))
	end
	return boxMax
end

fontname(font::Font) = "$(font.family)_$(font.style)($(font.size.x)x$(font.size.y))"

const ftLib = (FT_Library)[C_NULL]

function init()
	global ftLib
	@assert ftLib[1] == C_NULL
	err = FT_Init_FreeType(ftLib)
	return err == 0
end

function done()
	global ftLib
	@assert ftLib[1] != C_NULL
	err = FT_Done_FreeType(ftLib[1])
	ftLib[1] = C_NULL
	return err == 0
end

function max_glyph_size(face::FT_Face, faceRec::FreeType.FT_FaceRec, chars)
	width, height = 0, 0
	for c in chars
		err = FT_Load_Char(face, c, FT_LOAD_RENDER)
		@assert err == 0
		glyphRec = unsafe_load(faceRec.glyph)
		width = max(width, glyphRec.bitmap.width)
		height = max(height, glyphRec.bitmap.rows)
	end
	return width, height
end

function glyph_bitmap(bmpRec::FreeType.FT_Bitmap)
	@assert bmpRec.pixel_mode == FreeType.FT_PIXEL_MODE_GRAY
	bmp = Array{UInt8}(undef, bmpRec.width, bmpRec.rows)
	row = bmpRec.buffer
	if bmpRec.pitch < 0
		row -= bmpRec.pitch * (rbmpRec.rows - 1)
	end

	for r = 1:bmpRec.rows
		srcArray = unsafe_wrap(Array{UInt8}, row, bmpRec.width) #pointer_to_array(row, bmpRec.width)
		bmp[:, r] = srcArray
		row += bmpRec.pitch
	end
	return bmp
end

function get_kerning(face::FT_Face, c1::Char, c2::Char, divisor::Float32)
	i1 = FT_Get_Char_Index(face, c1)
	i2 = FT_Get_Char_Index(face, c2)
	kernVec = Array{FreeType.FT_Vector}(undef, 1)
	err = FT_Get_Kerning(face, i1, i2, FreeType.FT_KERNING_DEFAULT, kernVec)
	if err != 0
		return zero(Vec2{Float})
	end
	return Vec2{Float32}(kernVec[1].x / divisor, kernVec[1].y / divisor)
end

function loadfont(faceName::AbstractString; sizeXY::Tuple{Real, Real} = (32, 32), faceIndex::Real = 0, chars = '\u0000':'\u00ff')
	face = (FT_Face)[C_NULL]
	err = FT_New_Face(ftLib[1], faceName, Int32(faceIndex), face)
	if err != 0
		@info("Couldn't load font $faceName with error $err")
		return nothing
	end

	err = FT_Set_Pixel_Sizes(face[1], UInt32(sizeXY[1]), UInt32(sizeXY[2]))
	font = nothing
	if err != 0
		@info("Couldn't set the pixel size for font $faceName with error $err")
	else
		faceRec = unsafe_load(face[1])

		maxCharWidth, maxCharHeight = max_glyph_size(face[1], faceRec, chars)

		emScale = Float32(sizeXY[2]) / faceRec.units_per_EM
		lineDist = round(faceRec.height * emScale)
		ascent = round(faceRec.ascender * emScale)
		descent = round(faceRec.descender * emScale)

		font = Font(unsafe_string(faceRec.family_name),
					unsafe_string(faceRec.style_name),
					sizeXY[1], sizeXY[2],
					lineDist, ascent, descent,
					maxCharWidth, maxCharHeight, length(chars))

		# load glyphs
		charPos = GlyphPosition()
		for c in chars
			err = FT_Load_Char(face[1], c, FT_LOAD_RENDER)
			@assert err == 0
			glyphRec = unsafe_load(faceRec.glyph)
			@assert glyphRec.format == FreeType.FT_GLYPH_FORMAT_BITMAP
			glyphBmp = glyph_bitmap(glyphRec.bitmap)

			addglyph(font, c, glyphBmp,
					 Vec2{Float32}(-glyphRec.bitmap_left, glyphRec.bitmap_top),
					 Vec2{Float32}(glyphRec.advance.x / 64f0, glyphRec.advance.y / 64f0),
					 charPos)
		end

		# query kerning info
		if faceRec.face_flags & FreeType.FT_FACE_FLAG_KERNING != 0
			kernDivisor = (faceRec.face_flags & FreeType.FT_FACE_FLAG_SCALABLE != 0) ? 64f0 : 1f0
			for c1 in chars, c2 in chars
				kerning = get_kerning(face[1], c1, c2, kernDivisor)
				if kerning != zero(Vec2{Float32})
					addkerning(font, c1, c2, kerning)
				end
			end
		end

		adding_done(font)
	end

	err = FT_Done_Face(face[1])
	@assert err == 0
	return font
end

end


text:
using Documenter
using JLD2

makedocs(;
    modules = [JLD2],
    sitename = "Julia Data Format",
    build = "build",
    clean = true,
    doctest = true,
    #root="./",
    #src="src",
    repo = "https://github.com/JuliaIO/JLD2.jl",
    pages = [
        "Basics" => "index.md",
	"Custom Serialization" => "customserialization.md",
    "Compression" => "compression.md",
	"Internals & Design" => "internals.md",
	"HDF5 Compatibility" => "hdf5compat.md",
    "Advanced Usage" => "advanced.md",
    "Legacy" => "legacy.md",
    	],
    )

deploydocs(
    repo = "github.com/JuliaIO/JLD2.jl.git",
    target = "build",
    push_preview = true,
)


text:
import Dates, LinearAlgebra
Base.@deprecate_binding invokelatest Base.invokelatest false
Base.@deprecate_binding MathConstants Base.MathConstants false
Base.@deprecate_binding Fix2 Base.Fix2 false
Base.@deprecate_binding Sys Base.Sys false
Base.@deprecate_binding AbstractDateTime Dates.AbstractDateTime false
Base.@deprecate_binding tr LinearAlgebra.tr false
Base.@deprecate_binding IteratorSize Base.IteratorSize false
Base.@deprecate_binding IteratorEltype Base.IteratorEltype false
Base.@deprecate_binding opnorm LinearAlgebra.opnorm false
Base.@deprecate_binding norm LinearAlgebra.norm false
Base.@deprecate_binding dot LinearAlgebra.dot false
Base.@deprecate_binding (⋅) LinearAlgebra.dot false
Base.@deprecate_binding notnothing Base.notnothing false
Base.@deprecate_binding qr LinearAlgebra.qr false
Base.@deprecate_binding rmul! LinearAlgebra.rmul! false

for stdlib in [:Base64, :Dates, :DelimitedFiles, :Distributed, :InteractiveUtils, :Libdl,
    :LibGit2, :LinearAlgebra, :Markdown, :Mmap, :Pkg, :Printf, :Random, :REPL,
    :Serialization, :SharedArrays, :Sockets, :SparseArrays, :Statistics, :Test, :Unicode,
    :UUIDs]
    @eval begin
        import $stdlib
        Base.@deprecate_binding $stdlib $stdlib false
    end
end

Base.@deprecate_binding macros_have_sourceloc true false
Base.@deprecate enable_debug(x::Bool) x false

module TypeUtils
    Base.@deprecate_binding isabstract isabstracttype
    Base.@deprecate_binding typename Base.typename
    Base.@pure function _parameter_upper_bound(t::UnionAll, idx)
        return Base.rewrap_unionall((Base.unwrap_unionall(t)::DataType).parameters[idx], t)
    end
    Base.@deprecate_binding(parameter_upper_bound, _parameter_upper_bound, true,
        ", rewrite your code to use static parameters in dispatch or use `Base.rewrap_unionall((Base.unwrap_unionall(t)::DataType).parameters[idx], t)`.")
end # module TypeUtils
Base.@deprecate_binding TypeUtils TypeUtils false ", call the respective Base functions directly"


text:
# possible chart types
@enum tb_chart_type tb_multiline=1 tb_margin=2

"""
    log_custom_scalar(logger, layout::AbstractDict; step = step(logger))

Groups multiple scalars in the same plot to be visualized by the CUSTOM_SCALARS
plugin. Note that this function sets the metadata: the actual values must be
logged separately with `log_value` and referenced with the correct tag.

The `layout` argument is structured as follows:

    layout = Dict(category => Dict(name => (chart_type, [tag1, tag2, ...])))

where `category` is the main tag for the plot, `name` is the plot's name,
`chart_type` is one between `tb_multiline` and `tb_margin` and the array of tags
contains the actual references to the logged scalars.
"""
function log_custom_scalar(logger::TBLogger, layout::AbstractDict; step = nothing)
    summ = SummaryCollection(custom_scalar_summary(layout))
    write_event(logger.file, make_event(logger, summ, step = step))
end

function chart(name::String, metadata::Tuple{tb_chart_type, AbstractArray})
    chart_type, tags = metadata

    if chart_type == tb_multiline
        content = MultilineChartContent(tag = tags)
        return Chart(title = name, multiline = content)
    elseif chart_type == tb_margin
        @assert length(tags) == 3
        args = Dict(k => v for (k, v) in zip([:value, :lower, :upper], tags))
        content = MarginChartContent(
                    series = [MarginChartContent_Series(; args...)])
        return Chart(title = name, margin = content)
    else
        @error "The chart type must be `tb_multiline` or `tb_margin`"
    end
end

function charts(dict::AbstractDict)
    [chart(name, meta) for (name, meta) in zip(keys(dict), values(dict))]
end

function custom_scalar_summary(layout)
    cat_spec = zip(keys(layout), values(layout))
    categories = [Category(title = k, chart = charts(c)) for (k, c) in cat_spec]

    layout = Layout(category = categories)
    plugin_data = SummaryMetadata_PluginData(plugin_name = "custom_scalars")
    smd = SummaryMetadata(plugin_data = plugin_data)
    cs_tensor = TensorProto(dtype = _DataType.DT_STRING,
                            string_val = [serialize_proto(layout)],
                            tensor_shape = TensorShapeProto())

    Summary_Value(tag = "custom_scalars__config__",
                  tensor = cs_tensor,
                  metadata = smd)
end


text:
using Turing, ContinuousBenchmarks

include(joinpath(ContinuousBenchmarks.DATA_DIR, "toy-data", "gauss.data.jl"))
include(joinpath(ContinuousBenchmarks.MODELS_DIR, "toy-models", "gauss.model.jl"))

@tbenchmark(PG(20, 2000), gaussmodel, gaussdata)
#bench_res = @tbenchmark(PG(20, 2000), gaussmodel, gaussdata)
#logd = build_logd("Gaussian Model", bench_res...)
#print_log(logd)

@tbenchmark(HMC(20, 0.25, 5), gaussmodel, gaussdata)
#bench_res = @tbenchmark(HMC(2000, 0.25, 5), gaussmodel, gaussdata)
#logd = build_logd("Gaussian Model", bench_res...)
#print_log(logd)

@tbenchmark(Gibbs(20, HMC(10, 0.25, 5, :mu), PG(20, 10, :lam)), gaussmodel, gaussdata)
bench_res = @tbenchmark(Gibbs(200, HMC(10, 0.25, 5, :mu), PG(20, 10, :lam)), gaussmodel, gaussdata)
logd = build_logd("Gaussian Model", bench_res...)
print_log(logd)
send_log(logd)


text:
#This function uses linked lists relationed to the initial_point to reduces the memory usage and speed up the LJ potential energy calculations

function linkedlists(p, Data)
  Nt = trunc(Int64, Data.side / Data.cutoff)
  first_atom = zeros(Int64,Nt, Nt)
  next_atom = zeros(Int64, Data.N)

  for iat in 1:Data.N
    icell = trunc(Int64, p[iat][1]/Data.cutoff) + 1
    jcell = trunc(Int64, p[iat][2]/Data.cutoff) + 1
    next_atom[iat] = first_atom[icell,jcell]
    first_atom[icell,jcell] = iat
  end
 return Nt, first_atom,next_atom
end


text:
using Printf, LinearAlgebra, Statistics, Distributed, PositiveFactorizations,
      Distances, NearestNeighbors, DelimitedFiles, StatsBase

using .GP

abstract type Options end

mutable struct OptionsStat <: Options
    nmin                :: Int
    nmax                :: Int
    xbounds             :: Array{Float64, 2}
    fbounds             :: Array{Float64, 2}
    xall                :: Array{Float64, 2}
    λ²                  :: Array{Float64}
    δ                   :: Float64
    demean              :: Bool
    sdev_prop           :: Array{Float64, 1}
    sdev_pos            :: Array{Float64, 1}
    sdev_dc             :: Array{Float64, 1}
    pnorm               :: Float64
    stat_window         :: Int
    dispstatstoscreen   :: Bool
    report_freq         :: Int
    save_freq           :: Int
    fdataname           :: String
    history_mode        :: String
    costs_filename      :: String
    fstar_filename      :: String
    x_ftrain_filename   :: String
    debug               :: Bool
    quasimultid         :: Bool
    influenceradius     :: Array{Float64, 1}
    K                   :: GP.Kernel
    kdtree              :: KDTree
    balltree            :: BallTree
    timesλ              :: Float64
    needλ²fromlog       :: Bool
    updatenonstat       :: Bool
    updatenuisances     :: Bool
    peskycholesky       :: Bool
    dcvalue             :: Array{Float64}
    sampledc            :: Bool
end

function OptionsStat(;
        nmin               = 1,
        nmax               = 50,
        xbounds            = [0 1.794; 0 1.2],
        fbounds            = [-1 1],
        xall               = nothing,
        λ                 = [0.3],
        δ                 = 0.1,
        demean             = true,
        sdev_prop          = [0.01],
        sdev_pos           = [0.05;0.05],
        sdev_dc            = nothing,
        pnorm              = 2,
        stat_window        = 100,
        dispstatstoscreen  = true,
        report_freq        = 10,
        save_freq          = 100,
        history_mode       = "w",
        fdataname          = "",
        debug              = false,
        quasimultid        = "",
        influenceradius    = [-9.9],
        K                  = SqEuclidean(),
        timesλ             = 3.7,
        needλ²fromlog      = false,
        updatenonstat      = false,
        updatenuisances    = false,
        peskycholesky      = true,
        dcvalue            = nothing,
        sampledc           = false
        )

        @assert xall != nothing
        @assert all(diff(xbounds, dims=2) .> 0)
        @assert all(diff(fbounds, dims=2) .> 0)
        @assert ndims(sdev_prop) == 1
        @assert length(sdev_prop) == size(fbounds, 1)
        @assert ndims(sdev_pos) == 1
        @assert length(sdev_pos) == size(xbounds, 1)
        @assert length(λ) == size(xbounds, 1)
        if sdev_dc == nothing
            sdev_dc = copy(sdev_prop)
        else
            @assert length(sdev_dc) == size(fbounds, 1)
        end
        @assert ndims(λ) == 1
        @assert quasimultid != "" "specify true or false explicitly"
        if quasimultid
            @assert influenceradius[1] > 0.0
            @assert length(influenceradius) == size(xall, 1) - 1
        end
        @assert typeof(K) <: GP.Kernel
        if dcvalue == nothing
            dcvalue = mean(fbounds, dims=2)
        else
            @assert all((dcvalue .- fbounds[:,1]).>0)
            @assert all((dcvalue .- fbounds[:,2]).<0)
        end
        if needλ²fromlog && updatenonstat
            @assert sampledc == false
            @assert demean == false
        end
        if sampledc == true
            @assert demean == false
        end
        costs_filename = "misfits_"*fdataname*".bin"
        fstar_filename = "models_"*fdataname*".bin"
        x_ftrain_filename = "points_"*fdataname*".bin"
        kdtree = KDTree(xall)
        balltree = BallTree(xall./λ)
        OptionsStat(nmin, nmax, xbounds, fbounds, xall, λ.^2 , δ, demean, sdev_prop, sdev_pos, sdev_dc, pnorm,
                stat_window, dispstatstoscreen, report_freq, save_freq,
                fdataname, history_mode, costs_filename, fstar_filename, x_ftrain_filename,
                debug, quasimultid, influenceradius, K, kdtree, balltree, timesλ, needλ²fromlog,
                updatenonstat, updatenuisances, peskycholesky, dcvalue, sampledc)
end

mutable struct OptionsNonstat <: Options
    nmin                :: Int
    nmax                :: Int
    xbounds             :: Array{Float64, 2}
    fbounds             :: Array{Float64, 2}
    xall                :: Array{Float64, 2}
    δ                   :: Float64
    demean              :: Bool
    sdev_prop           :: Array{Float64, 1}
    sdev_pos            :: Array{Float64, 1}
    sdev_dc             :: Array{Float64, 1}
    pnorm               :: Float64
    stat_window         :: Int
    dispstatstoscreen   :: Bool
    report_freq         :: Int
    save_freq           :: Int
    fdataname           :: String
    history_mode        :: String
    costs_filename      :: String
    fstar_filename      :: String
    x_ftrain_filename   :: String
    debug               :: Bool
    quasimultid         :: Bool
    influenceradius     :: Array{Float64, 1}
    K                   :: GP.Kernel
    kdtree              :: KDTree
    needλ²fromlog       :: Bool
    updatenonstat       :: Bool
    updatenuisances     :: Bool
    peskycholesky       :: Bool
    dcvalue             :: Array{Float64}
    sampledc            :: Bool
end

function OptionsNonstat(opt::OptionsStat;
        nmin               = 1,
        nmax               = 50,
        fbounds            = [-1 1],
        δ                  = 0.1,
        demean             = true,
        sdev_prop          = [0.01],
        sdev_pos           = [0.05;0.05],
        sdev_dc            = nothing,
        pnorm              = 2,
        influenceradius    = [-9.9],
        K                  = SqEuclidean(),
        dcvalue            = nothing,
        sampledc           = false
        )

        @assert all(diff(fbounds, dims=2) .> 0)
        @assert ndims(sdev_prop) == 1
        @assert length(sdev_prop) == size(fbounds, 1)
        @assert ndims(sdev_pos) == 1
        @assert length(sdev_pos) == size(opt.xbounds, 1)
        if sdev_dc == nothing
            sdev_dc = copy(sdev_prop)
        else
            @assert length(sdev_dc) == size(fbounds, 1)
        end
        @assert typeof(K) <: GP.Kernel
        if dcvalue == nothing
            dcvalue = mean(fbounds, dims=2)
        else
            @assert all((dcvalue .- fbounds[:,1]).>0)
            @assert all((dcvalue .- fbounds[:,2]).<0)
        end
        if sampledc == true
            @assert demean == false
        end
        costs_filename = "misfits_ns_"*opt.fdataname*".bin"
        fstar_filename = "models_ns_"*opt.fdataname*".bin"
        x_ftrain_filename = "points_ns_"*opt.fdataname*".bin"
        OptionsNonstat(nmin, nmax, opt.xbounds, fbounds, opt.xall, δ, demean, sdev_prop, sdev_pos, sdev_dc, pnorm,
                opt.stat_window, opt.dispstatstoscreen, opt.report_freq, opt.save_freq,
                opt.fdataname, opt.history_mode, opt.costs_filename, opt.fstar_filename, opt.x_ftrain_filename,
                opt.debug, opt.quasimultid, influenceradius, K, opt.kdtree, opt.needλ²fromlog, opt.updatenonstat,
                opt.updatenuisances, opt.peskycholesky, dcvalue, sampledc)
end

mutable struct OptionsNuisance
    # not subtyping Options because methods that expect Options are expecting
    # options for a Gaussian process model, not fixed-dimension MCMC
    sdev                   :: Array{Float64,1}
    bounds                 :: Array{Float64,2}
    nnu                    :: Int64

    updatenuisances        :: Bool
    updatenonstat          :: Bool
    debug                  :: Bool

    stat_window            :: Int
    dispstatstoscreen      :: Bool

    report_freq            :: Int
    save_freq              :: Int
    history_mode           :: String
    costs_filename         :: String
    vals_filename          :: String
    fdataname              :: String

    W                      :: Array{Float64, 2}
    idxnotzero             :: Array{Int, 1}
    rotatebounds           :: Array{Float64, 2}
    Xbar                   :: Array{Float64, 1}
end

function OptionsNuisance(opt::OptionsStat;
        C = nothing, nsamples=100_000, updatenonstat=false,
        sdev = [0.1, 0.1],
        bounds = [0 1.; 0 1.],
        updatenuisances = false)

    @assert all(sdev .>= 0.)
    @assert all(diff(bounds, dims=2) .>= 0)
    @assert length(sdev) == size(bounds, 1)
    @assert all(sdev .< 1.) # is a fraction of the bounds
    nnu = length(sdev)

    # make an empty nuisance options struct
    optn = OptionsNuisance(sdev, bounds, nnu, updatenuisances,
     updatenonstat, opt.debug, opt.stat_window, opt.dispstatstoscreen, opt.report_freq, opt.save_freq, opt.history_mode, "misfits_nuisance_"*opt.fdataname*".bin",
     "values_nuisance_"*opt.fdataname*".bin", opt.fdataname, [0 0.], [0], [0. 0.], [0.])

    # fill remaining fields with sensible default values
    idxnotzero = zeros(Int,0)
    for i = 1:optn.nnu
        numin, numax = extrema(optn.bounds[i,:])
        if abs(numin-numax)>1e-12
             push!(idxnotzero, i)
        end
    end
    optn.idxnotzero = findidxnotzero(optn.nnu, optn.bounds)
    nnonzero = length(idxnotzero)

    # Identity rotation matrix if norough covariance estimate provided
    W = Matrix(1.0I, nnonzero, nnonzero)
    if C != nothing
        W = eigen(C).vectors
    end
    optn.W = W

    # Generate random samples from uniform nuisance priors
    X = zeros(nsamples, nnonzero)
    for (i,idx) in enumerate(idxnotzero)
        X[:,i] = optn.bounds[idx, 1] .+
                   diff(optn.bounds[idx,:])[:].*rand(nsamples)
    end

    # Demean them and transform to principal directions
    Xbar = mean(optn.bounds[idxnotzero,:], dims=2)
    optn.Xbar = vec(Xbar)
    X = X .- Xbar'
    T = X*W

    # Find the bounds in rotated space
    rotatebounds = zeros(nnonzero, 2)
    for i in 1:nnonzero
        rotatebounds[i,:] .= extrema(T[:,i])
    end
    optn.rotatebounds = rotatebounds

    # scale the bounds to the sdev fraction
    optn.sdev[idxnotzero] = optn.sdev[idxnotzero].*
                            diff(rotatebounds, dims=2)[:]
    optn
end

function findidxnotzero(nnu, bounds)
    idxnotzero = zeros(Int,0)
    for i = 1:nnu
        numin, numax = extrema(bounds[i,:])
        if abs(numin-numax)>1e-12
            push!(idxnotzero, i)
        end
    end
    idxnotzero
end

abstract type Model
    #"Model" means a Gaussian-process parametrised function
end

mutable struct ModelStat <:Model
    fstar         :: Array{Float64}
    xtrain        :: Array{Float64}
    ftrain        :: Array{Float64}
    K_y           :: Array{Float64, 2}
    Kstar         :: Array{Float64, 2}
    n             :: Int
    ftrain_old    :: Array{Float64}
    xtrain_old    :: Array{Float64}
    iremember     :: Int # stores old changed point index to recover state
    xtrain_focus  :: Array{Float64} # for quasimultid to keep track of
    dcvalue       :: Array{Float64}
    dcvalue_old   :: Array{Float64}
end

mutable struct ModelNonstat <:Model
    fstar         :: Array{Float64}
    xtrain        :: Array{Float64}
    ftrain        :: Array{Float64}
    K_y           :: Array{Float64, 2}
    Kstar         :: Array{Float64, 2}
    n             :: Int
    ftrain_old    :: Array{Float64}
    xtrain_old    :: Array{Float64}
    iremember     :: Int # stores old changed point index to recover state
    xtrain_focus  :: Array{Float64} # for quasimultid to keep track of
    K_y_old       :: Array{Float64, 2}
    Kstar_old     :: Array{Float64, 2}
    dcvalue       :: Array{Float64}
    dcvalue_old   :: Array{Float64}
end

mutable struct ModelNuisance
    # "nuisances" are a fixed-dimensional part of the model, with the only
    # allowed move being a value change. Contrasting with the GP parametrisation
    # of the main part of the model which allows birth and death of GP nodes
    nuisance      :: Array{Float64}
    nu_old        :: Array{Float64} #prev val of nidx_mem
end

mutable struct Stats
    move_tries::Array{Int, 1}
    accepted_moves::Array{Int, 1}
    accept_rate::Array{Float64, 1}
end

function Stats(;nmoves=5)
    Stats(zeros(Int, nmoves), zeros(Int, nmoves), zeros(Float64, nmoves))
end

mutable struct Writepointers
    fp_costs      :: IOStream
    fp_fstar      :: IOStream
    fp_x_ftrain   :: IOStream
end

mutable struct Writepointers_nuisance
    fp_costs      :: IOStream
    fp_vals       :: IOStream
end

# Stationary GP functions, i.e., for λ
function init(opt::OptionsStat)
    n, xtrain, ftrain, dcvalue = initvalues(opt)
    K_y = zeros(opt.nmax, opt.nmax)
    map!(x->GP.κ(opt.K, x),K_y,pairwise(WeightedEuclidean(1 ./opt.λ² ), xtrain, dims=2))
    K_y[diagind(K_y)] .+= opt.δ^2
    Kstar = zeros(Float64, size(opt.xall,2), opt.nmax)
    xtest = opt.xall
    map!(x->GP.κ(opt.K, x),Kstar,pairwise(WeightedEuclidean(1 ./opt.λ² ), xtest, xtrain, dims=2))
    fstar = zeros(size(opt.fbounds, 1), size(opt.xall, 2))
    calcfstar!(fstar, ftrain, opt, K_y, Kstar, n, dcvalue)
    return ModelStat(fstar, xtrain, ftrain, K_y, Kstar, n,
                 [0.0], zeros(Float64, size(opt.xbounds, 1)), 0,
                 zeros(Float64, size(opt.xbounds, 1)), dcvalue, copy(dcvalue))
end

function init(opt::OptionsNuisance)
    nuisance = zeros(0)
    if opt.nnu > 0 # not a dummy option
        if opt.history_mode == "w" # fresh start
            nuisance = opt.bounds[:,1] + diff(opt.bounds, dims = 2)[:].*rand(opt.nnu)
        else # is a restart
            # TODO hacky I don't like this
            @info "reading $(opt.vals_filename)"
            nuisance = vec(readdlm(opt.vals_filename)[end,2:1+opt.nnu])
        end
    end
    return ModelNuisance(nuisance, copy(nuisance))
end

function getmean(opt, n, ftrain, dcvalue)
    if opt.demean && n>1
        mf = mean(ftrain[:,1:n], dims=2)
    else
        if opt.sampledc
            mf = dcvalue # from what is stored in model and passed to it
        else
            mf = opt.dcvalue # supplied by user or mid point of prior range by default
        end
    end
    mf
end

function commoncalc(ftrain::Array{Float64,2}, opt::Options, K_y::Array{Float64,2},
    Kstar::Array{Float64, 2}, n::Int, dcvalue::Array{Float64})
    mf = getmean(opt, n, ftrain, dcvalue)
    rhs = ftrain[:,1:n] .- mf
    ky = @view K_y[1:n,1:n]
    if opt.peskycholesky
        U = cholesky(Positive, ky, Val{false}).U
    else
        U = cholesky(ky).U
    end
    ks = @view Kstar[:,1:n]
    rhs, mf, U, ks
end

function calcfstar!(fstar::Array{Float64,2}, ftrain::Array{Float64,2},
                    opt::OptionsStat, K_y::Array{Float64,2},
                    Kstar::Array{Float64, 2}, n::Int, dcvalue::Array{Float64})
    rhs, mf, U, ks = commoncalc(ftrain, opt, K_y, Kstar, n, dcvalue)
    if opt.needλ²fromlog
        copy!(fstar, 10 .^(2(mf' .+ ks*(U\(U'\rhs'))))' )
    else
        copy!(fstar, (mf' .+ ks*(U\(U'\rhs')))' )
    end
    nothing
end

function initvalues(opt::Options)
    xtrain = zeros(Float64, size(opt.xbounds,1), opt.nmax)
    ftrain = zeros(Float64, size(opt.fbounds,1), opt.nmax)
    if opt.history_mode == "w" # new start
        n = opt.nmin
        xtrain[:,1:n] = opt.xbounds[:,1] .+ diff(opt.xbounds, dims=2).*rand(size(opt.xbounds, 1), n)
        ftrain[:,1:n] = opt.fbounds[:,1] .+ diff(opt.fbounds, dims=2).*rand(size(opt.fbounds, 1), n)
        dcvalue       = opt.fbounds[:,1] .+ diff(opt.fbounds, dims=2).*rand(size(opt.fbounds, 1), 1)
    else # restart
        @info "opening $(opt.x_ftrain_filename)"
        n = history(opt, stat=:nodes)[end]
        xft = history(opt, stat=:x_ftrain)[end]
        dcvalue = history(opt, stat=:dcvalue)[end,:]
        xtrain[:,1:n] = xft[1:size(opt.xall, 1), 1:n]
        ftrain[:,1:n] = xft[size(opt.xall, 1)+1:end, 1:n]
    end
    n, xtrain, ftrain, dcvalue
end

function updatenskernels!(opt::OptionsStat, m::ModelStat, ipoint::Union{Int, Array{Int, 1}},
                          optns::OptionsNonstat, mns::ModelNonstat; doall=false, isposchange=false)
    xtrain, λ², xtest = m.xtrain, m.fstar, opt.xall
    copyto!(mns.K_y_old, CartesianIndices((mns.n, mns.n)), mns.K_y, CartesianIndices((mns.n, mns.n)))
    copyto!(mns.Kstar_old, CartesianIndices((size(mns.Kstar_old, 1), mns.n)),
            mns.Kstar, CartesianIndices((size(mns.Kstar_old, 1), mns.n)))
    xt = xtrain[:,ipoint]
    isposchange && (xt = hcat(xt, m.xtrain_old))

    if doall
        kstarchangeidx = 1:size(mns.Kstar, 1)
        kychangeidx = 1:mns.n
    else
        kstarchangeidx = unique(vcat(inrange(opt.balltree, xt./sqrt.(opt.λ²), opt.timesλ)...))
        balltree = BallTree(mns.xtrain[:,1:mns.n]./sqrt.(opt.λ²))
        kychangeidx = inrange(balltree, xt./sqrt.(opt.λ²), opt.timesλ)
    end
    idxs = gettrainidx(opt.kdtree, mns.xtrain, mns.n)
    # changes where test points are in influence radius of lscale change
    ks = view(mns.Kstar, kstarchangeidx, 1:mns.n)
    λ²train = @view λ²[:,idxs]
    λ²test = @view λ²[:,kstarchangeidx]
    xtr = @view mns.xtrain[:,1:mns.n]
    xte = @view xtest[:,kstarchangeidx]
    GP.pairwise(ks, optns.K, xtr, xte, λ²train, λ²test)
    if length(kychangeidx) > 0
        isposchange && (kychangeidx = reduce(vcat, kychangeidx))
        # set 1 of changes where training points are in influence radius of lscale change
        ks2 = view(mns.Kstar, :, kychangeidx)
        xtr2 = @view mns.xtrain[:,kychangeidx]
        λ²train2 = @view λ²[:,idxs][:,kychangeidx]
        GP.pairwise(ks2, optns.K, xtr2, xtest, λ²train2, λ²)
        # set 2 of changes where training points are in influence radius of lscale change
        ky = view(mns.K_y, 1:mns.n , kychangeidx)
        GP.pairwise(ky, optns.K, xtr2, xtr, λ²train2, λ²train)
        mns.K_y[kychangeidx,1:mns.n] .= mns.K_y[1:mns.n,kychangeidx]'
        # nugget add
        ky = view(mns.K_y, 1:mns.n , 1:mns.n)
        ky[diagind(ky)] .= 1. + optns.δ^2
    end
    sync_model!(mns, optns)
    nothing
end

function birth!(m::ModelStat, opt::OptionsStat)
    # stationary only
    xtrain, ftrain, K_y,  Kstar, n = m.xtrain, m.ftrain, m.K_y,  m.Kstar, m.n
    xtrain[:,n+1] = opt.xbounds[:,1] + diff(opt.xbounds, dims=2).*rand(size(opt.xbounds, 1))
    copy!(m.xtrain_focus, xtrain[:,n+1])
    ftrain[:,n+1] = opt.fbounds[:,1] + diff(opt.fbounds, dims=2).*rand(size(opt.fbounds, 1))
    xtest = opt.xall
    Kstarv = @view Kstar[:,n+1]
    map!(x->GP.κ(opt.K, x),Kstarv,colwise(WeightedEuclidean(1 ./opt.λ² ), xtrain[:,n+1], xtest))
    K_yv = @view K_y[n+1,1:n+1]
    map!(x->GP.κ(opt.K, x),K_yv,colwise(WeightedEuclidean(1 ./opt.λ² ), xtrain[:,n+1], xtrain[:,1:n+1]))
    K_y[1:n+1,n+1] = K_y[n+1,1:n+1]
    K_y[n+1,n+1] = K_y[n+1,n+1] + opt.δ^2
    calcfstar!(m.fstar, m.ftrain, opt, K_y, Kstar, n+1, m.dcvalue)
    m.n = n+1
    nothing
end

function birth!(m::ModelStat, opt::OptionsStat, mns::ModelNonstat, optns::OptionsNonstat; doall=false)
    # stationary and nonstat together but for stat model update
    birth!(m, opt)
    # updating the nonstationary kernels now
    updatenskernels!(opt, m, m.n, optns, mns, doall=doall)
    nothing
end

function undo_birth!(m::ModelStat)
    m.n = m.n - 1
    nothing
end

function undo_birth!(m::ModelStat, mns::ModelNonstat)
    # stationary and nonstat together but for stat model update
    undo_birth!(m)
    # updating the nonstationary kernels now
    copyto!(mns.K_y, CartesianIndices((mns.n, mns.n)), mns.K_y_old, CartesianIndices((mns.n, mns.n)))
    copyto!(mns.Kstar, CartesianIndices((size(mns.Kstar_old, 1), mns.n)),
                mns.Kstar_old, CartesianIndices((size(mns.Kstar_old, 1), mns.n)))
    nothing
end

function death!(m::Model, opt::Options)
    xtrain, ftrain, K_y,  Kstar, n = m.xtrain, m.ftrain, m.K_y,  m.Kstar, m.n
    ipoint = rand(1:n)
    m.iremember = ipoint
    copy!(m.xtrain_focus, xtrain[:,ipoint])
    xtrain[:,ipoint], xtrain[:,n] = xtrain[:,n], xtrain[:,ipoint]
    ftrain[:,ipoint], ftrain[:,n] = ftrain[:,n], ftrain[:,ipoint]
    Kstar[:,ipoint], Kstar[:,n] = Kstar[:,n], Kstar[:,ipoint]
    K_y[ipoint,1:n], K_y[n,1:n] = K_y[n,1:n], K_y[ipoint,1:n]
    K_y[1:n-1,ipoint] = K_y[ipoint,1:n-1]
    K_y[ipoint,ipoint] = 1.0 + opt.δ^2
    calcfstar!(m.fstar, m.ftrain, opt, K_y, Kstar, n-1, m.dcvalue)
    m.n = n-1
    nothing
end

function death!(m::ModelStat, opt::OptionsStat,
                mns::ModelNonstat, optns::OptionsNonstat; doall=false)
    # stationary and nonstat together but for stat model update
    death!(m, opt)
    # updating the nonstationary kernels now
    updatenskernels!(opt, m, m.n+1, optns, mns, doall=doall)
    nothing
end

function undo_death!(m::Model, opt::Options)
    m.n = m.n + 1
    ftrain, xtrain, Kstar, K_y, n, ipoint = m.ftrain, m.xtrain, m.Kstar, m.K_y, m.n, m.iremember
    xtrain[:,ipoint], xtrain[:,n] = xtrain[:,n], xtrain[:,ipoint]
    ftrain[:,ipoint], ftrain[:,n] = ftrain[:,n], ftrain[:,ipoint]
    Kstar[:,ipoint], Kstar[:,n] = Kstar[:,n], Kstar[:,ipoint]
    K_y[ipoint,1:n], K_y[n,1:n] = K_y[n,1:n], K_y[ipoint,1:n]
    K_y[1:n,ipoint] = K_y[ipoint,1:n]
    K_y[ipoint,ipoint] = 1.0 + opt.δ^2
    K_y[n,n] = 1.0 + opt.δ^2
    nothing
end

function undo_death!(m::ModelStat, opt::OptionsStat, mns::ModelNonstat)
    undo_death!(m, opt)
    # updating the nonstationary kernels now
    copyto!(mns.K_y, CartesianIndices((mns.n, mns.n)), mns.K_y_old, CartesianIndices((mns.n, mns.n)))
    copyto!(mns.Kstar, CartesianIndices((size(mns.Kstar_old, 1), mns.n)),
            mns.Kstar_old, CartesianIndices((size(mns.Kstar_old, 1), mns.n)))
    nothing
end

function property_change!(m::Model, opt::Options)
    ftrain, K_y, Kstar, n = m.ftrain, m.K_y, m. Kstar, m.n
    ipoint = 1 + floor(Int, rand()*n)
    m.iremember = ipoint
    copy!(m.ftrain_old, ftrain[:,ipoint])
    ftrain[:,ipoint] = ftrain[:,ipoint] + opt.sdev_prop.*randn(size(opt.fbounds, 1))
    for i in eachindex(ftrain[:,ipoint])
        while (ftrain[i,ipoint]<opt.fbounds[i,1]) || (ftrain[i,ipoint]>opt.fbounds[i,2])
                (ftrain[i,ipoint]<opt.fbounds[i,1]) && (ftrain[i,ipoint] = 2*opt.fbounds[i,1] - ftrain[i,ipoint])
                (ftrain[i,ipoint]>opt.fbounds[i,2]) && (ftrain[i,ipoint] = 2*opt.fbounds[i,2] - ftrain[i,ipoint])
        end
    end
    calcfstar!(m.fstar, m.ftrain, opt, K_y, Kstar, n, m.dcvalue)
    nothing
end

function property_change!(m::ModelStat, opt::OptionsStat,
                          mns::ModelNonstat, optns::OptionsNonstat; doall=false)
    property_change!(m, opt)
    # updating the nonstationary kernels now
    updatenskernels!(opt, m, m.iremember, optns, mns, doall=doall)
    nothing
end

function undo_property_change!(m::Model)
    ipoint, ftrain = m.iremember, m.ftrain
    ftrain[:,ipoint] = m.ftrain_old
    nothing
end

function undo_property_change!(m::ModelStat, mns::ModelNonstat)
    undo_property_change!(m)
    # updating the nonstationary kernels now
    copyto!(mns.K_y, CartesianIndices((mns.n, mns.n)), mns.K_y_old, CartesianIndices((mns.n, mns.n)))
    copyto!(mns.Kstar, CartesianIndices((size(mns.Kstar_old, 1), mns.n)),
            mns.Kstar_old, CartesianIndices((size(mns.Kstar_old, 1), mns.n)))
    nothing
end

function position_change!(m::ModelStat, opt::OptionsStat)
    xtrain, ftrain, K_y, Kstar, n = m.xtrain, m.ftrain, m.K_y, m.Kstar, m.n
    ipoint = 1 + floor(Int, rand()*n)
    m.iremember = ipoint
    m.xtrain_old[:] = xtrain[:,ipoint]
    xtrain[:,ipoint] = xtrain[:,ipoint] + opt.sdev_pos.*randn(size(opt.xbounds, 1))
    copy!(m.xtrain_focus, xtrain[:,ipoint])
    for i in eachindex(xtrain[:,ipoint])
        while (xtrain[i,ipoint]<opt.xbounds[i,1]) || (xtrain[i,ipoint]>opt.xbounds[i,2])
                (xtrain[i,ipoint]<opt.xbounds[i,1]) && (xtrain[i,ipoint] = 2*opt.xbounds[i,1] - xtrain[i,ipoint])
                (xtrain[i,ipoint]>opt.xbounds[i,2]) && (xtrain[i,ipoint] = 2*opt.xbounds[i,2] - xtrain[i,ipoint])
        end
    end
    xtest = opt.xall
    Kstarv = @view Kstar[:,ipoint]
    map!(x->GP.κ(opt.K, x),Kstarv,colwise(WeightedEuclidean(1 ./opt.λ² ), xtrain[:,ipoint], xtest))
    K_yv = @view K_y[ipoint,1:n]
    map!(x->GP.κ(opt.K, x),K_yv,colwise(WeightedEuclidean(1 ./opt.λ² ), xtrain[:,ipoint], xtrain[:,1:n]))
    K_y[1:n,ipoint] = K_y[ipoint,1:n]
    K_y[ipoint,ipoint] = K_y[ipoint,ipoint] + opt.δ^2
    calcfstar!(m.fstar, m.ftrain, opt, K_y, Kstar, n, m.dcvalue)
    nothing
end

function position_change!(m::ModelStat, opt::OptionsStat,
                          mns::ModelNonstat, optns::OptionsNonstat; doall=false)
    position_change!(m, opt)
    # updating the nonstationary kernels now
    updatenskernels!(opt, m, m.iremember, optns, mns, doall=doall, isposchange=true)
    nothing
end

function undo_position_change!(m::ModelStat, opt::OptionsStat)
    xtrain, K_y, Kstar, n = m.xtrain, m.K_y, m.Kstar, m.n
    ipoint = m.iremember
    xtrain[:,ipoint] = m.xtrain_old
    xtest = opt.xall
    Kstarv = @view Kstar[:,ipoint]
    map!(x->GP.κ(opt.K, x),Kstarv,colwise(WeightedEuclidean(1 ./opt.λ² ), xtrain[:,ipoint], xtest))
    K_yv = @view K_y[ipoint,1:n]
    map!(x->GP.κ(opt.K, x),K_yv,colwise(WeightedEuclidean(1 ./opt.λ² ), xtrain[:,ipoint], xtrain[:,1:n]))
    K_y[1:n,ipoint] = K_y[ipoint,1:n]
    K_y[ipoint,ipoint] = K_y[ipoint,ipoint] + opt.δ^2
    nothing
end

function undo_position_change!(m::ModelStat, opt::OptionsStat, mns::ModelNonstat)
    undo_position_change!(m, opt)
    # updating the nonstationary kernels now
    copyto!(mns.K_y, CartesianIndices((mns.n, mns.n)), mns.K_y_old, CartesianIndices((mns.n, mns.n)))
    copyto!(mns.Kstar, CartesianIndices((size(mns.Kstar_old, 1), mns.n)),
            mns.Kstar_old, CartesianIndices((size(mns.Kstar_old, 1), mns.n)))
    nothing
end

function dc_change!(m::Model, opt::Options)
    ftrain, K_y, Kstar, n, dcvalue = m.ftrain, m.K_y, m. Kstar, m.n, m.dcvalue
    copy!(m.dcvalue_old, dcvalue)
    dcvalue .= m.dcvalue_old .+ opt.sdev_dc.*randn(size(opt.fbounds, 1))
    for i in eachindex(dcvalue)
        while (dcvalue[i]<opt.fbounds[i,1]) || (dcvalue[i]>opt.fbounds[i,2])
                (dcvalue[i]<opt.fbounds[i,1]) && (dcvalue[i] = 2*opt.fbounds[i,1] - dcvalue[i])
                (dcvalue[i]>opt.fbounds[i,2]) && (dcvalue[i] = 2*opt.fbounds[i,2] - dcvalue[i])
        end
    end
    calcfstar!(m.fstar, m.ftrain, opt, K_y, Kstar, n, m.dcvalue)
end

function undo_dc_change!(m::Model)
    m.dcvalue .= m.dcvalue_old
    nothing
end

# Non stationary GP functions, i.e., for mns.fstar
function gettrainidx(kdtree::KDTree, xtrain::Array{Float64, 2}, n::Int)
    idxs,  = knn(kdtree, xtrain[:,1:n], 1)
    reduce(vcat, idxs)
end

function init(opt::OptionsNonstat, m::ModelStat)
    donotinit = !opt.needλ²fromlog && !opt.updatenonstat
    if !donotinit
        λ² = m.fstar
        n, xtrain, ftrain, dcvalue = initvalues(opt)
        K_y = zeros(opt.nmax, opt.nmax)
        idxs = gettrainidx(opt.kdtree, xtrain, n)
        ky = view(K_y, 1:n, 1:n)
        GP.pairwise(ky, opt.K, xtrain[:,1:n], xtrain[:,1:n], λ²[:,idxs], λ²[:,idxs])
        K_y[diagind(K_y)] .+= opt.δ^2
        Kstar = zeros(Float64, size(opt.xall,2), opt.nmax)
        xtest = opt.xall
        ks = view(Kstar, :, 1:n)
        GP.pairwise(ks, opt.K, xtrain[:,1:n], xtest, λ²[:,idxs], λ²)
        fstar = zeros(size(opt.xall, 2), size(opt.fbounds, 1))
        calcfstar!(fstar, ftrain, opt, K_y, Kstar, n, dcvalue)
        return ModelNonstat(fstar, xtrain, ftrain, K_y, Kstar, n,
                 [0.0], zeros(Float64, size(opt.xbounds, 1)), 0, zeros(Float64, size(opt.xbounds, 1)),
                 copy(K_y), copy(Kstar), dcvalue, copy(dcvalue))
    else
        dummy2d = [0.0 0.0]
        return ModelNonstat(dummy2d, dummy2d, dummy2d, dummy2d, dummy2d, 0,
                 dummy2d, dummy2d, 0, dummy2d,
                 dummy2d, dummy2d, dummy2d, dummy2d)
    end
end

function calcfstar!(fstar::Array{Float64,2}, ftrain::Array{Float64,2},
                    opt::OptionsNonstat, K_y::Array{Float64,2},
                    Kstar::Array{Float64, 2}, n::Int, dcvalue::Array{Float64})
    rhs, mf, U, ks = commoncalc(ftrain, opt, K_y, Kstar, n, dcvalue)
    copy!(fstar, mf' .+ ks*(U\(U'\rhs')) )
    nothing
end

function birth!(m::ModelNonstat, opt::OptionsNonstat, l::ModelStat)
    λ² = l.fstar
    xtrain, ftrain, K_y,  Kstar, n = m.xtrain, m.ftrain, m.K_y,  m.Kstar, m.n
    xtrain[:,n+1] = opt.xbounds[:,1] + diff(opt.xbounds, dims=2).*rand(size(opt.xbounds, 1))
    copy!(m.xtrain_focus, xtrain[:,n+1])
    ftrain[:,n+1] = opt.fbounds[:,1] + diff(opt.fbounds, dims=2).*rand(size(opt.fbounds, 1))
    xtest = opt.xall
    Kstarv = @view Kstar[:,n+1]
    idxs = gettrainidx(opt.kdtree, xtrain, n+1)
    K_yv = @view K_y[n+1,1:n+1]
    @views begin
    GP.colwise!(Kstarv, opt.K, xtrain[:,n+1], xtest, λ²[:,idxs[end]], λ²)
    GP.colwise!(K_yv, opt.K, xtrain[:,n+1], xtrain[:,1:n+1], λ²[:,idxs[end]], λ²[:,idxs])
    end
    K_y[1:n+1,n+1] = K_y[n+1,1:n+1]
    K_y[n+1,n+1] = K_y[n+1,n+1] + opt.δ^2
    calcfstar!(m.fstar, m.ftrain, opt, K_y, Kstar, n+1, m.dcvalue)
    m.n = n+1
    nothing
end

function undo_birth!(m::ModelNonstat)
    m.n = m.n - 1
    nothing
end

function position_change!(m::ModelNonstat, opt::OptionsNonstat, l::ModelStat)
    λ² = l.fstar
    xtrain, ftrain, K_y, Kstar, n = m.xtrain, m.ftrain, m.K_y, m.Kstar, m.n
    ipoint = 1 + floor(Int, rand()*n)
    m.iremember = ipoint
    m.xtrain_old[:] = xtrain[:,ipoint]
    xtrain[:,ipoint] = xtrain[:,ipoint] + opt.sdev_pos.*randn(size(opt.xbounds, 1))
    copy!(m.xtrain_focus, xtrain[:,ipoint])
    for i in eachindex(xtrain[:,ipoint])
        while (xtrain[i,ipoint]<opt.xbounds[i,1]) || (xtrain[i,ipoint]>opt.xbounds[i,2])
                (xtrain[i,ipoint]<opt.xbounds[i,1]) && (xtrain[i,ipoint] = 2*opt.xbounds[i,1] - xtrain[i,ipoint])
                (xtrain[i,ipoint]>opt.xbounds[i,2]) && (xtrain[i,ipoint] = 2*opt.xbounds[i,2] - xtrain[i,ipoint])
        end
    end
    xtest = opt.xall
    Kstarv = @view Kstar[:,ipoint]
    idxs = gettrainidx(opt.kdtree, xtrain, n)
    K_yv = @view K_y[ipoint,1:n]
    @views begin
    GP.colwise!(Kstarv, opt.K, xtrain[:,ipoint], xtest, λ²[:,idxs[ipoint]], λ²)
    GP.colwise!(K_yv, opt.K, xtrain[:,ipoint], xtrain[:,1:n], λ²[:,idxs[ipoint]], λ²[:,idxs])
    end
    K_y[1:n,ipoint] = K_y[ipoint,1:n]
    K_y[ipoint,ipoint] = K_y[ipoint,ipoint] + opt.δ^2
    calcfstar!(m.fstar, m.ftrain, opt, K_y, Kstar, n, m.dcvalue)
    nothing
end

function undo_position_change!(m::ModelNonstat, opt::OptionsNonstat, l::ModelStat)
    λ² = l.fstar
    xtrain, K_y, Kstar, n = m.xtrain, m.K_y, m.Kstar, m.n
    ipoint = m.iremember
    xtrain[:,ipoint] = m.xtrain_old
    xtest = opt.xall
    Kstarv = @view Kstar[:,ipoint]
    idxs = gettrainidx(opt.kdtree, xtrain, n)
    K_yv = @view K_y[ipoint,1:n]
    @views begin
    GP.colwise!(Kstarv, opt.K, xtrain[:,ipoint], xtest, λ²[:,idxs[ipoint]], λ²)
    GP.colwise!(K_yv, opt.K, xtrain[:,ipoint], xtrain[:,1:n], λ²[:,idxs[ipoint]], λ²[:,idxs])
    end
    K_y[1:n,ipoint] = K_y[ipoint,1:n]
    K_y[ipoint,ipoint] = K_y[ipoint,ipoint] + opt.δ^2
    nothing
end

function sync_model!(m::Model, opt::Options)
    ftrain, K_y, Kstar, n = m.ftrain, m.K_y, m.Kstar, m.n
    calcfstar!(m.fstar, m.ftrain, opt, K_y, Kstar, n, m.dcvalue)
    nothing
end

function testupdate(optns::OptionsNonstat, l::ModelStat, mns::ModelNonstat)
    λ² = l.fstar
    mf = getmean(optns, mns.n, mns.ftrain, mns.dcvalue)
    idxs = gettrainidx(optns.kdtree, mns.xtrain, mns.n)
    ftest, = GP.GPfit(optns.K, mns.ftrain[:,1:mns.n], mns.xtrain[:,1:mns.n],
            optns.xall, λ², λ²[:,idxs], optns.δ, p=2, demean=optns.demean, nogetvars=true,
            my=mf)
    return ftest
end

function testupdate(opt::OptionsStat, m::ModelStat)
    mf = getmean(opt, m.n, m.ftrain, m.dcvalue)
    ftest, = GP.GPfit(opt.K, m.ftrain[:,1:m.n], m.xtrain[:,1:m.n],
            opt.xall, opt.λ², opt.δ, nogetvars=true, demean=opt.demean, p=2,
            my=mf)
    return ftest
end

function do_move!(mns::ModelNonstat, m::ModelStat, optns::OptionsNonstat, statns::Stats)
    # purely nonstationary GP moves
    cumprob = [0.20, 0.40, 0.60, 0.8] # if sampledc and 5 total noves
    if !optns.sampledc
        cumprob +=  [0.05, 0.10, 0.15, 0.2] # add these if no sampledc, 4 total moves
    end
    unifrand = rand()
    movetype, priorviolate = 0, false
    if unifrand<cumprob[1]
        if mns.n<optns.nmax
            birth!(mns, optns, m)
        else
            priorviolate = true
        end
        movetype = 1
    elseif unifrand<cumprob[2]
        if mns.n>optns.nmin
            death!(mns, optns)
        else
            priorviolate = true
        end
        movetype = 2
    elseif unifrand<cumprob[3]
        position_change!(mns, optns, m)
        movetype = 3
    elseif unifrand<cumprob[4]
        property_change!(mns, optns)
        movetype = 4
    else
        dc_change!(mns, optns)
        movetype = 5
    end
    statns.move_tries[movetype] += 1
    return movetype, priorviolate
end

function undo_move!(movetype::Int, mns::ModelNonstat, optns::OptionsNonstat, m::ModelStat)
    # purely nonstationary GP moves
    if movetype == 1
        undo_birth!(mns)
    elseif movetype == 2
        undo_death!(mns, optns)
    elseif movetype == 3
        undo_position_change!(mns, optns, m)
    elseif movetype == 4
        undo_property_change!(mns)
    else
        undo_dc_change!(mns)
    end
    sync_model!(mns, optns)
    nothing
end

function do_move!(m::ModelStat, opt::OptionsStat, stat::Stats,
                  mns::ModelNonstat, optns::OptionsNonstat)
    # Stationary GP changes which update nonstationary GP
    cumprob = [0.20, 0.40, 0.60, 0.8] # if sampledc and 5 total noves
    if !opt.sampledc
        cumprob +=  [0.05, 0.10, 0.15, 0.2] # add these if no sampledc, 4 total moves
    end
    unifrand = rand()
    movetype, priorviolate = 0, false
    if unifrand<cumprob[1]
        if m.n<opt.nmax
            birth!(m, opt, mns, optns)
        else
            priorviolate = true
        end
        movetype = 1
    elseif unifrand<cumprob[2]
        if m.n>opt.nmin
            death!(m, opt, mns, optns)
        else
            priorviolate = true
        end
        movetype = 2
    elseif unifrand<cumprob[3]
        position_change!(m, opt, mns, optns)
        movetype = 3
    elseif unifrand<cumprob[4]
        property_change!(m, opt, mns, optns)
        movetype = 4
    else
        dc_change!(m, opt)
        movetype = 5
    end
    stat.move_tries[movetype] += 1
    return movetype, priorviolate
end

function undo_move!(movetype::Int, m::ModelStat, opt::OptionsStat,
                    mns::ModelNonstat, optns::OptionsNonstat)
    # Stationary GP changes which update nonstationary GP
    if movetype == 1
        undo_birth!(m, mns)
    elseif movetype == 2
        undo_death!(m, opt, mns)
    elseif movetype == 3
        undo_position_change!(m, opt, mns)
    elseif movetype == 4
        undo_property_change!(m, mns)
    else
        undo_dc_change!(m)
    end
    sync_model!(m, opt)
    opt.updatenonstat && sync_model!(mns, optns)
    nothing
end

function do_move!(m::ModelStat, opt::OptionsStat, stat::Stats)
    # purely stationary GP moves
    cumprob = [0.20, 0.40, 0.60, 0.8] # if sampledc and 5 total noves
    if !opt.sampledc
        cumprob +=  [0.05, 0.10, 0.15, 0.2] # add these if no sampledc, 4 total moves
    end
    unifrand = rand()
    movetype, priorviolate = 0, false
    if unifrand<cumprob[1]
        if m.n<opt.nmax
            birth!(m, opt)
        else
            priorviolate = true
        end
        movetype = 1
    elseif unifrand<cumprob[2]
        if m.n>opt.nmin
            death!(m, opt)
        else
            priorviolate = true
        end
        movetype = 2
    elseif unifrand<cumprob[3]
        position_change!(m, opt)
        movetype = 3
    elseif unifrand<cumprob[4]
        property_change!(m, opt)
        movetype = 4
    else
        dc_change!(m, opt)
        movetype = 5
    end
    stat.move_tries[movetype] += 1
    return movetype, priorviolate
end

function undo_move!(movetype::Int, m::ModelStat, opt::OptionsStat)
    # purely stationary GP moves
    if movetype == 1
        undo_birth!(m)
    elseif movetype == 2
        undo_death!(m, opt)
    elseif movetype == 3
        undo_position_change!(m, opt)
    elseif movetype == 4
        undo_property_change!(m)
    else
        undo_dc_change!(m)
    end
    sync_model!(m, opt)
    opt.updatenonstat && sync_model!(mns, optns)
    nothing
end

function do_oneaxis_move!(mn::ModelNuisance, optn::OptionsNuisance, statn::Stats)
    nidx = rand(optn.idxnotzero)
    new_nval = mn.nuisance[nidx] + optn.sdev[nidx]*randn()
    while (new_nval<optn.bounds[nidx,1]) || (new_nval>optn.bounds[nidx,2])
        (new_nval<optn.bounds[nidx,1]) && (new_nval = 2*optn.bounds[nidx,1] - new_nval)
        (new_nval>optn.bounds[nidx,2]) && (new_nval = 2*optn.bounds[nidx,2] - new_nval)
    end

    mn.nu_old[:] = mn.nuisance
    mn.nuisance[nidx] = new_nval
    statn.move_tries[nidx] += 1
    return nidx
end

function do_move!(mn::ModelNuisance, optn::OptionsNuisance, statn::Stats)
    idxnotzero = optn.idxnotzero
    boundsmean = optn.Xbar
    W = optn.W

    newoptn = deepcopy(optn)
    newoptn.bounds[idxnotzero,:] = optn.rotatebounds
    newmn = deepcopy(mn)
    newnuisance = (mn.nuisance[idxnotzero] - boundsmean)'*W
    newmn.nuisance[idxnotzero] = newnuisance'
    nidx = do_oneaxis_move!(newmn, newoptn, statn)
    newnuisance = W*newmn.nuisance[idxnotzero] + boundsmean
    mn.nu_old[:] = mn.nuisance
    mn.nuisance[idxnotzero] = newnuisance
    nidx
end

function undo_move!(mn::ModelNuisance)
    mn.nuisance[:] = mn.nu_old
end

function get_acceptance_stats!(isample::Int, opt::Options, stat::Stats)
    if mod(isample-1, opt.stat_window) == 0
        stat.accept_rate[:] = 100. *stat.accepted_moves./stat.move_tries
        if opt.dispstatstoscreen
            msg = @sprintf("ARs Birth %5.2f Death %5.2f Position %5.2f Property %5.2f DC %5.2f",
                            stat.accept_rate[1],
                            stat.accept_rate[2],
                            stat.accept_rate[3],
                            stat.accept_rate[4],
                            stat.accept_rate[5])
            @info typeof(opt), msg
        end
        fill!(stat.move_tries, 0)
        fill!(stat.accepted_moves, 0)
    end
end

function get_acceptance_stats!(isample::Int, opt::OptionsNuisance, stat::Stats)
    if mod(isample-1, opt.stat_window) == 0
        stat.accept_rate[:] = 100. *stat.accepted_moves./stat.move_tries
        ar = stat.accept_rate
        if opt.dispstatstoscreen
            @info typeof(opt), round.(Int, ar[.!isnan.(ar)])
        end
        fill!(stat.move_tries, 0)
        fill!(stat.accepted_moves, 0)
    end
end

# history methods
function open_history(opt::Options)
    if isfile(opt.costs_filename)
        @assert (opt.history_mode=="a") "$(opt.costs_filename) exists"
    end
    if opt.report_freq > 0
        @info("running transD_sampler...")
    end
    fp_costs = nothing
    if length(opt.costs_filename) > 0
        fp_costs = open(opt.costs_filename, opt.history_mode)
    end
    fp_models = nothing
    if length(opt.fstar_filename) > 0
        fp_fstar = open(opt.fstar_filename, opt.history_mode)
    end
    fp_x_ftrain = nothing
    if length(opt.x_ftrain_filename) > 0
        fp_x_ftrain = open(opt.x_ftrain_filename, opt.history_mode)
    end
    return Writepointers(fp_costs, fp_fstar, fp_x_ftrain)
end

#history for nuisance models (which have different internal structure)
function open_history(optn::OptionsNuisance)
    if isfile(optn.costs_filename)
        @assert (optn.history_mode == "a") "$(optn.costs_filename) exists"
    end
    if optn.report_freq > 0
        @info("running nuisance sampler...")
    end
    fp_costs = nothing
    if length(optn.costs_filename) > 0
        fp_costs = open(optn.costs_filename, optn.history_mode)
    end
    fp_vals = nothing
    if length(optn.vals_filename) > 0
        fp_vals = open(optn.vals_filename, optn.history_mode)
    end
    return Writepointers_nuisance(fp_costs, fp_vals)
end

function close_history(wp::Writepointers)
    if wp.fp_costs != nothing
        close(wp.fp_costs)
    end
    if wp.fp_fstar != nothing
        close(wp.fp_fstar)
    end
    if wp.fp_x_ftrain != nothing
        close(wp.fp_x_ftrain)
    end
    @info "closed files"
end

function setrestartflag(opt)
    opt.history_mode = "a"
end

function close_history(wpn::Writepointers_nuisance)
    if wpn.fp_costs != nothing
        close(wpn.fp_costs)
    end
    if wpn.fp_vals != nothing
        close(wpn.fp_vals)
    end
    @info "closed files"
end

function clear_history(opt::Options)
    if length(opt.costs_filename) != 0 && isfile(opt.costs_filename) == true
        rm(opt.costs_filename)
    end
    if length(opt.fstar_filename) != 0 && isfile(opt.fstar_filename) == true
        rm(opt.models_decompr_filename)
    end
    if length(opt.x_ftrain_filename) != 0 && isfile(opt.x_ftrain_filename) == true
        rm(opt.x_ftrain_filename)
    end
end

function mode_history(opt::Options, mode::String)
    @assert mode == "w" || mode == "a"
    opt.history_mode = mode
end

function write_history(isample::Int, opt::Options, m::Model, misfit::Float64,
                        stat::Stats, wp::Writepointers, T::Float64, writemodel::Bool)
    write_history(opt, m.fstar, [m.xtrain; m.ftrain], m.dcvalue, misfit, stat.accept_rate[1],
                        stat.accept_rate[2], stat.accept_rate[3], stat.accept_rate[4], stat.accept_rate[5], m.n,
                       isample, wp.fp_costs, wp.fp_fstar, wp.fp_x_ftrain, T, writemodel)
end

function write_history(isample::Int, optn::OptionsNuisance, mn::ModelNuisance, misfit::Float64,
                    statn::Stats, wpn::Writepointers_nuisance, T::Float64, writemodel::Bool)

    write_history(optn, mn.nuisance, misfit, statn.accept_rate, isample, wpn.fp_costs,
                wpn.fp_vals, T, writemodel)
end

function write_history(opt::Options, fstar::AbstractArray, x_ftrain::AbstractArray, dcvalue::AbstractArray, U::Float64, acceptanceRateBirth::Float64,
                    acceptanceRateDeath::Float64, acceptanceRatePosition::Float64, acceptanceRateProperty::Float64, ARdc, nodes::Int,
                    iter::Int, fp_costs::Union{IOStream, Nothing}, fp_fstar::Union{IOStream, Nothing},
                    fp_x_ftrain::Union{IOStream, Nothing}, T::Float64, writemodel::Bool)
    if (mod(iter-1, opt.save_freq) == 0 || iter == 1)
        if fp_costs != nothing
            msg = @sprintf("%d %e %e %e %e %e %d %e %e", iter, acceptanceRateBirth, acceptanceRateDeath,
                                        acceptanceRatePosition, acceptanceRateProperty, ARdc, nodes, U, T)
            for dc in dcvalue # saves dcvalue vector after last cost in msg above on same line
                msg *= @sprintf(" %e", dc)
            end
            msg *= "\n"
            write(fp_costs, msg)
            flush(fp_costs)
        end
        if fp_x_ftrain != nothing
            write(fp_x_ftrain, convert(Array{eltype(Float64)},x_ftrain))
            flush(fp_x_ftrain)
        end
        if writemodel
            if fp_fstar != nothing
                write(fp_fstar, convert(Array{eltype(Float64)},fstar))
                flush(fp_fstar)
            end
        end
    end
end

function write_history(optn::OptionsNuisance, nvals::Array{Float64,1}, misfit::Float64,
                    acceptanceRate::Array{Float64, 1}, iter::Int, fp_costs::Union{IOStream, Nothing},
                    fp_vals::Union{IOStream, Nothing}, T::Float64, writemodel::Bool)
    if (mod(iter-1, optn.save_freq) == 0 || iter == 1)
        ars = acceptanceRate # TODO hacky for now, would like all acceptance rates
        if fp_costs != nothing
            msg = @sprintf("%d %e %e", iter, misfit, T)
            for ar in ars
                msg *= @sprintf(" %e", ar)
            end
            msg *= "\n"
            write(fp_costs, msg)
            flush(fp_costs)
        end
        if fp_vals != nothing
            msg = @sprintf("%d", iter)
            for nval = nvals
                msg *= @sprintf(" %e", nval)
            end
            msg *= "\n"
            write(fp_vals, msg)
            flush(fp_vals)
        end
    end
end

function history(optn::OptionsNuisance; stat=:misfit)
    idxlast = 0
    for (statname, idx) in ((:iter,   1),
                            (:misfit, 2),
                            (:T,      3))
        if stat == statname
            data = readdlm(optn.costs_filename, ' ', Float64)[:,idx]
        end
        idxlast = idx
    end
    if stat == :iter
        data = Int.(data)
    end
    if stat == :acceptanceRate
        data = readdlm(optn.costs_filename, ' ', Float64)[:,idxlast+1:end]
    end
    return data
end

function history(opt::Options; stat=:U)
    idxlast = 0
    for (statname, el, idx) in ((:iter,                   Int,      1),
                                (:acceptanceRateBirth,    Float64,  2),
                                (:acceptanceRateDeath,    Float64,  3),
                                (:acceptanceRatePosition, Float64,  4),
                                (:acceptanceRateProperty, Float64,  5),
                                (:acceptanceRateDC,       Float64,  6),
                                (:nodes,                  Int,      7),
                                (:U,                      Float64,  8),
                                (:T,                      Float64,  9))

        if stat == statname
            if length(opt.costs_filename) == 0
                @warn("history, requested $(statname), but you haven't stored this information.")
                return []
            end
            fp_costs = open(opt.costs_filename)
            mark(fp_costs)
            X = Array{el}(undef, countlines(fp_costs))
            reset(fp_costs)
            i = 1
            for line in readlines(fp_costs)
                X[i] = el == Float64 ? parse(Float64,split(line)[idx]) : parse(Int,split(line)[idx])
                i += 1
            end
            return X
        end
        idxlast = idx # here only if noe of the statnames in statname are met
    end
    if stat == :dcvalue # assumes dcvalue is stored last after all the regular costs
        return readdlm(opt.costs_filename)[:,idxlast+1:end]
    end
    if stat == :fstar
        if length(opt.fstar_filename) == 0
            @warn("history, requested fstar, but you haven't stored this information.")
            return []
        end
        iters, rem = divrem(filesize(opt.fstar_filename), size(opt.xall,2) * size(opt.fbounds, 1) * sizeof(Float64))
        @assert rem == 0
        fp_models = open(opt.fstar_filename)
        fstar = Array{Array{Float64}}(undef, iters)
        for i = 1:iters
            if typeof(opt) == OptionsStat
                fstar[i] = zeros(Float64, (size(opt.fbounds, 1), size(opt.xall,2)))
            else
                fstar[i] = zeros(Float64, (size(opt.xall,2), size(opt.fbounds, 1)))
            end
            read!(fp_models, fstar[i])
        end
        return fstar
    end
    if stat == :x_ftrain
        if length(opt.x_ftrain_filename) == 0
            @warn("history, requested x_ftrain, but you haven't stored this information.")
            return []
        end
        iters, rem = divrem(filesize(opt.x_ftrain_filename), opt.nmax * sizeof(Float64) * (size(opt.fbounds, 1) +
                                                                                size(opt.xbounds, 1)))
        @assert rem == 0
        fp_models = open(opt.x_ftrain_filename)
        x_ftrain = Array{Array{Float64,2},1}(undef, iters)
        for i = 1:iters
            x_ftrain[i] = zeros(Float64, size(opt.fbounds, 1) + size(opt.xbounds, 1), opt.nmax)
            read!(fp_models, x_ftrain[i])
        end
        return x_ftrain
    end
    @warn("history, requested stat: $(stat) is not recognized.")
    return []
end


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using CovidSurvey
using CSV
using DataFrames
using Dates
# using ArviZ
using PrettyTables
import PyPlot as plt
# using RCall
using Turing
using StatsPlots
using Plots.PlotMeasures
using Random
using StatsBase
using BSON
using ColorSchemes
using Measures
colors = ColorSchemes.tableau_10.colors

using Underscores
import Downloads
using RData

gr()
pgfplotsx()
## ==========================================================================
#          Introduction

fdir, ps, ignores = Regional.parse_fname(fname)

pp = let
	fname = "/home/ankoh/dev/CovidSurvey/reports/regional_prediction_end_of_peak/2021-09-07/CHAIN_chains=10_model=hospit_observ=2021-03-01_regional=true_rwstep=1_semipar=true_std=false_steps=500_warmup=2000.jls"
	chain = read(fname, Chains)

  data_params = (
        observationsend  = ps.observ
      , predictors        = ps.preds|> CovidSurvey.parse_predictors
      # , cases_start       = ps.cases
      # , hospitmodel       = Regional.HospitInit1(obs_stop="2020-07-01")
      , casemodel         = Regional.CaseInit(obs_start="2020-06-01")
      , seromodel         = Regional.SeroInit(delay=0, std=1.)
      , rwstep            = ps.rwstep
      , epidemicstart    = 20
      , numimpute        = 6
      , include_early_dynamic = false#ps.include
      , link              = KLogistic(4.5)
      , invlink           = KLogit(4.5)
      , predictive        = false
      , covariates_kwargs = Dict(
        :semiparametric => ps.semipar,
        :fname          => [projectdir("data/contacts/dk/", "mean_contact_rates_region=$(r).csv" ) for r in 1:Regional.nregions],#[projectdir("data", "mean_contact_rates_region=$r.csv" ) for r in 1:Regional.nregions],
        :shift          => -1,
        :startdate      => "2020-11-10", # >= 1000 cases / day
        :datecol        => "dates",
        :aggregation    => ps.rwstep,
        :standartize    => ps.std,
        :mobility       => [projectdir("data/mobility/mobility_region=$(Regional.regions[i]).csv") for i in 1:Regional.nregions]
        )
      )
  data = Regional.load_data(; data_params... )

	model = Regional.name2model[ps.model]
	Regional.PostProcessing2(fdir, ps, ignores, chain, data, model)
end

function Plots.plot!(p::Plots.Plot, r::Regional.ObservationPlottingRecipe1; plot_only_fit=false)
    o  = r.observed
    e  = r.expected
    pr = r.predicted
    ed = r.enddate
    sd = r.startdate

    if !plot_only_fit
        !isnothing(ed) && vline!(p, [ed], lab="end observations", lw=2, lc=:black, hover="$ed")
        !isnothing(ed) && vline!(p, [sd], lab="start observations", lw=2, lc=:black, hover="$sd", ls=:dash)
    end
    plot!(p, o.dates, o.values, α=0.5, lc=:match, lab="observed $(r.label)", c=:midnightblue, lw=4, ylab="cases")
    # plot_confidence_timeseries!(p, pr.dates, pr.values; label = "predicted $(r.label)", c=:midnightblue) #Dict(hover=>strdates)
    plot_confidence_timeseries!(p, e.dates, e.values; c=colors[2], label = "expected $(r.label)") #Dict(hover=>strdates)
    if plot_only_fit
        xlo = Dates.value.(sd)
        xup = Dates.value.(ed)
        yup = maximum( o.values[o.dates .< ed] ) * 1.1
        xlims!(p, xlo, xup)
        ylims!(p, 0, yup)
    end
    return p
end

p = let region = "capital"
    gp = Regional.generate_posterior(pp)
    @unpack fdir, ps, ignores, chain, data = pp

    p1 = let
        recipe = Regional.RtPlottingRecipe1(data, gp, region, "Rt")
        p = plot(title = "Reproduction Number", bottom_margin = 20mm)
        plot!(p, recipe; plot_only_fit=true)
    end

    p2 = let
        r = recipe = Regional.LatentInfectionsPlottingRecipe3(pp.data, gp, region, "infections")
        p = plot(title = "Latent Infections", bottom_margin = 20mm)
        plot!(p, recipe; plot_only_fit=true)
    end

    p3 = let
        recipe = Regional.ObservationPlottingRecipe1(pp.data, gp, region, "hospitalizations")
        p = plot(title = "Observed Hospitalizations")
        plot!(p, recipe; plot_only_fit=true)
    end

    plot(p1, p2, p3, layout=(3,1), size=(800, 3*250), sharex=true, link=:x, legend = nothing, dpi=900)
end

savefig(p, projectdir("figures/SSI-talk/intro.png"))

## ==========================================================================
#          Qualitative Analysis
# ==========================================================================
function StatsBase.zscore(df::DataFrame)
    y = similar(df)
    for n in names(df)
        y[!,n] = zscore(df[:,n])
    end
    return y
end

function conform(args...; normalize = false, stopdate=nothing, startdate=nothing)
    s = isnothing(startdate) ? maximum( first.( select.(args, Ref(:date) ) ) ).date : Date(startdate)
    e = isnothing(stopdate) ? minimum( last.( select.(args, Ref(:date) ) ) ).date : Date(stopdate)
    x = s:Day(1):e
    y = filter.(Ref(:date => >=(s)), args)
    y = filter.(Ref(:date => <=(e)), y)
    normalize && ( y = zscore.( select.( y, Ref(Not(:date)) ) ) )
    (x, y...)
end

function load_regional_data(i; kwargs...)

    observables = Dict()
    Regional.observables!(observables)
    hospit = observables["hospit"]
    cases = observables["cases"]

    # fname = projectdir("reports/2021-08-01/Rt-$(uppercase(Regional.regions[i]))_chains=10_epidemicstart=10_folder=2021-08-01_model=cases_numimpute=6_rwstep=1_steps=100_warmup=500.csv")
    fname = projectdir("/home/ankoh/dev/CovidSurvey/reports/regional_prediction_end_of_peak/2021-09-07/Rt-$(uppercase(Regional.regions[i]))_chains=10_model=hospit_observ=2021-03-01_regional=true_rwstep=1_semipar=true_std=false_steps=500_warmup=2000.csv")
    rt = CSV.File(fname) |> DataFrame
    rename!(rt, "dates"=>"date")

    # fname = projectdir("data/raw_contact_rates_region=$(i)_aggregation=7.csv")
    fname = projectdir("data/contacts/dk/averaged_contact_rates_region=$i.csv")
    surv = CSV.File(fname) |> DataFrame
    "dates" in names(surv) && rename!(surv, "dates"=>"date")

    fname = projectdir("data/mobility/mobility_region=$(Regional.regions[i]).csv")
    mobil = CSV.File(fname) |> DataFrame

    return conform(cases, hospit, rt, surv, mobil; kwargs...)
end

function plot_comparison(date, cases, hospit, rt, surv, mobil; column=:country, plotactive = false)
    plots  = Vector{Plots.Plot}()
    ds = Date("2020-11-10")
    de = Date("2021-02-01")
    p = plot(; xaxis = true, title="observed hospitalizations", xticks=:native, yticks=:native, bottom_margin = 10mm)
    plotactive && vspan!([ds, de], label="active epidemic", α=0.2, fc=:midnightblue)
    #plot!(date, zscore(cases[:,column]), lab="cases", lw=2)
    plot!(date, hospit[:,column], lab="hospitalizations", lw=2, c=colors[1], legend=false)
    push!(plots, p)

    p = plot(; xaxis = true, title="reproduction number", xticks=:native, bottom_margin = 10mm,xformatter=_->"", legend=false)
    plot!(rt.date, rt.mq, ribbon=(rt.mq - rt.lq, rt.uq - rt.mq), c=colors[2], yticks=:native, lab="reproduction number", lw=2)
    plotactive && vspan!([ds, de], label="active epidemic", α=0.2, fc=:midnightblue)
    push!(plots, p)

    # p = plot(; xaxis = true, title="contact survey", xticks=:native, yticks=:native, bottom_margin = 0mm, legend = false)
    # for (i, c) in enumerate(["family", "colleagues", "friends"])#names(surv, Not(:date))
    #     plot!(surv.date, surv[:,c], lab=c, lw=2, c=colors[i+3])
    # end
    # plot!(surv.date, surv[:,"strangers"], lab="strangers", lw=2, c=colors[3])
    # plotactive && vspan!([ds, de], label="active epidemic", α=0.2, fc=:midnightblue)
    # push!(plots, p)

    p = plot(; xaxis = true, title="mobility", xticks=:native, yticks=:native, legend = :bottomleft)
    plot!(mobil.date, mobil.apple, lab="apple", lw=2, c=colors[5])
    plot!(mobil.date, mobil.telco, lab="telco", lw=2, c=colors[4])
    plot!(mobil.date, mobil.google, lab="google", lw=3, c=colors[3])
    mobil.telco
    plotactive && vspan!([ds, de], label="active epidemic", α=0.2, fc=:midnightblue)
    push!(plots, p)

    return plot(plots..., layout=(length(plots),1), size=(1200, length(plots)*250), sharex=true, link=:x, dpi=900)
end

i = 1
r = Regional.regions[i]
data  = load_regional_data(i; stopdate="2021-02-01")
p     = plot_comparison( data...; column=Regional.regions[i], plotactive=true)
savefig(p, projectdir("figures/SSI-talk/regionaldata_#4.png"))

## ==========================================================================
#         International Data Qualitative Analysis
# ==========================================================================

function parsetime(sraw)
    svec = split(sraw, "-")
    if length(svec) == 3
        d,m,y = svec
        sclean = join([d,m,"20$y"], "-")
        return Date(sclean, "d-u-y")
    else
        d = sraw[1:2]
        m = sraw[3:5]
        y = sraw[6:9]
        sclean = join([d,m,y], "-")
        return Date(sclean, "d-u-y")
    end
end

function readcontacts( fname )
    @_ fname |>
        CSV.read(__, DataFrame) |>
        DataFrames.select(__,
            :dataRound_glb => :batch,
            :edate => ByRow(x->parsetime(x)) => :date,
            "country_string" => :region,
            "qNew_3_2_resp" => :family,
            "qNew_3_3_resp" => :colleagues,
            "qNew_3_4_resp" => :friends,
            "qNew_3_5_resp" => :strangers,
        ) |>
    sort(__, :date)
end

function loadrt(fname)
    data = RData.load(fname)
    dates      = data["time"]
    values     = @_ data["draws"] |>
                    transpose |>
                    Array |>
                    __[:, 2001 : end]
    recipe = Regional.RtPlottingRecipe1( nothing, nothing, nothing, (; dates, values), "Rt")
    return (; recipe, dates)
end

let
    ctype = "strangers"
    regions = ["Hungary", "UK", "France", "Sweden", "Denmark", "Germany", "Italy", "USA"]
    for region in regions
        # plot raw data

        # plot case counts
        epidemic = CSV.read(projectdir("data/epidata/epidata_region=$(region).csv"), DataFrame)
        is = findfirst(>(22), epidemic.cases_per_100000)
        ds = epidemic.date[is]
        de1 = epidemic.date[is:end][findfirst(<(10), epidemic.cases_per_100000[is:end])]
        de2 = epidemic.date[findfirst(>(0.05), epidemic.fraction_vaccinated)]
        de  = min(de1, de2)
        p1 = plot(title="observed cases")
        vspan!([ds, de], label="active epidemic", α=0.2, fc=:midnightblue)
        plot!(epidemic.date, epidemic.deaths, fc=colors[1], lc=colors[1], label="smoothed cases / 100_000", lw=4)#, hover=epidemic.date)
        #plot!(twinx(), epidemic.date, epidemic.cases, fc=colors[1], lc=colors[1], label="cases", lw=0.8, α=0.5)

		fname = projectdir("reports/stanfit/2021-10-05/sampled_rt_region=$(region).rds")
		rt = loadrt(fname)
		p2 = plot(title="Reproduction Number")
		plot!(p2, rt.recipe)
        vspan!([ds, de], label="active epidemic", α=0.2, fc=:midnightblue)

		# plot predictors
		smoothed = CSV.read(projectdir("data/contacts/int/smoothed_contact_rates_region=$(region).csv"), DataFrame)
        # suffix = region != "Denmark" ? "" : "_low_resolution"
		mobility = CSV.read(projectdir("data/mobility/mobility_country=$(region).csv"), DataFrame)
        p3 = plot(title="predictors")
        vspan!([ds, de], label="active epidemic", α=0.2, fc=:midnightblue)
        plot!(p3, smoothed.date, zscore(smoothed[:,ctype]); label="contact to $ctype", hover=smoothed.date, lw=4)
        plot!(p3, mobility.date, zscore(mobility.apple); label="apple", lw=1, hover=mobility.date)
        plot!(p3, mobility.date, zscore(mobility.google); label="google", lw=1, hover=mobility.date)

		p4 = plot(title="fraction of vaccinated")
		vspan!([ds, de], label="active epidemic", α=0.2, fc=:midnightblue)
		plot!(epidemic.date, epidemic.fraction_vaccinated * 100, fc=colors[1], lc=colors[1], label="vaccinated", lw=4)#, hover=epidemic.date)

        p = plot(p1,p2,p3, p4, layout=(4,1), size=(1200, 4*250), sharex=true, link=:x, xlim=(rt.dates[1],Date("2021-06-01")), legend=false, dpi=900)
        # firefox(p, fname="$region.html")
		savefig(p, projectdir("figures/SSI-talk/international_data_region=$region.png"))
    end
end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
module BayesianLinearRegression

# using Statistics
using LinearAlgebra
using Measurements
using Printf
using LazyArrays

include("callbacks.jl")

export BayesianLinReg

mutable struct BayesianLinReg{T}
    XtX :: SubArray{Float64,2,Array{Float64,2},Tuple{Array{Int64,1},Array{Int64,1}},false}
    Xty :: SubArray{Float64,1,Array{Float64,1},Tuple{Array{Int64,1}},false}
    yty :: Float64
    N   :: Int

    Hinv:: SubArray{Float64,2,Array{Float64,2},Tuple{Array{Int64,1},Array{Int64,1}},false}

    eig :: LinearAlgebra.Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}

    priorPrecision :: T
    updatePrior :: Bool

    noisePrecision :: T
    updateNoise :: Bool

    weights

    iterationCount :: Int
    done :: Bool

    uncertaintyBasis :: Vector{Measurement{Float64}}

    active :: Vector{Int}
end

function symmetric!(S)
    S .+= S'
    S ./= 2
end

###########################################################################################
# Hessian

export hessian

function hessian(eig, α, β)
    Q = eig.vectors
    D = Diagonal(α .+ β .* eig.values)

    H = Q * D * Q'
    symmetric!(H)

    return H
end

function hessian(m::BayesianLinReg)
    α = m.priorPrecision
    β = m.noisePrecision

    return hessian(m.eig, α, β)
end

###########################################################################################
# Inverse Hessian

function hessianinv!(m::BayesianLinReg)
    Λ = m.eig.values
    Q = m.eig.vectors

    α = m.priorPrecision
    β = m.noisePrecision

    D = Diagonal(inv.(α .+ β .* Λ))

    Hinv = m.Hinv
    Hinv .= @~ Q * D * Q'
    symmetric!(Hinv)

    return Hinv
end

###########################################################################################
# Log-determinant of Hessian

function logdetH(m::BayesianLinReg)
    Λ = m.eig.values
    α = m.priorPrecision
    β = m.noisePrecision
    return logdetH(α, β, Λ)
end

function logdetH(α, β, Λ)
    sum(log, α .+ β .* Λ)
end

function updateWeights!(m::BayesianLinReg)
    β = m.noisePrecision
    Hinv = hessianinv!(m)

    m.weights .= @~ β .* (Hinv * m.Xty)
    return m.weights
end

function BayesianLinReg(
          X::Matrix{T}
        , y::Vector{T}
        ; updatePrior=true
        , updateNoise=true
        ) where {T}

    (N, p) = size(X)

    XtX = X' * X
    Xty = X' * y
    yty = normSquared(y)

    symmetric!(XtX)

    return BayesianLinReg(XtX, Xty, yty, N)
end

function BayesianLinReg(
      XtX
    , Xty
    , yty
    , N
    ; updatePrior=true
    , updateNoise=true)
    eig = eigen(XtX)

    p = length(Xty)
    ps = [1:p;]

    XtX = view(XtX, ps, ps)
    Xty = view(Xty, ps)

    Λ  = eig.values
    Λ .= max.(Λ, 0.0)
    Q = eig.vectors

    α = 1.0
    β = 1.0

    ps = [1:p;]
    D = Diagonal(inv.(α .+ β .* Λ))

    Hinv = Q * D * Q'
    symmetric!(Hinv)
    Hinv = view(Hinv, ps, ps)

    weights = view(β .* (Hinv * Xty), ps)

    BayesianLinReg(
        XtX
      , Xty
      , yty
      , N
      , Hinv
      , eig
      , α
      , updatePrior
      , β
      , updateNoise
      , weights
      , 0
      , false
      , zeros(p) .± 1.0
      , [1:p;]
  )

end

"""
Sum of squared residuals
"""
function ssr(m)
    XtX = m.XtX
    Xty = m.Xty
    yty = m.yty
    w = m.weights
    return yty - 2 * mydot(w, Xty) + w' * XtX * w
end

function Base.iterate(m::BayesianLinReg{T}, iteration=1) where {T}
    m.done && return nothing

    N = m.N
    ps = m.active
    α = m.priorPrecision
    β = m.noisePrecision

    gamma(m) = effectiveNumParameters(m)

    XtX = m.XtX
    Xty = m.Xty
    yty = m.yty

    w = m.weights

    wᵀw = normSquared(w)
    rᵀr = ssr(m)
    if m.updatePrior
        m.priorPrecision = gamma(m) / wᵀw
    end

    if m.updateNoise
        m.noisePrecision = (N - gamma(m)) / rᵀr
    end

    updateWeights!(m)

    m.iterationCount += 1
    return (m, iteration + 1)
end

export fit!

function fit!(m::BayesianLinReg; kwargs...)
    m.done = false
    callback = get(kwargs, :callback, stopAtIteration(10))

    if m.updatePrior
        m.priorPrecision = 1.0
    end

    if m.updateNoise
        m.noisePrecision = 1.0
    end

    try
        for iter in m
            callback(iter)
        end
    catch e
        if e isa InterruptException
            @warn "Computation interrupted"
            return m
        else
            rethrow()
        end
    end
    return m
end

export logEvidence

function logEvidence(m::BayesianLinReg{T}) where {T}
    N = m.N
    α = m.priorPrecision
    β = m.noisePrecision
    rtr = ssr(m)
    return _logEv(N, rtr, m.active, α, β, m.eig.values, m.weights)
end

const log2π = log(2π)

function _logEv(N, rtr, active, α, β, Λ, w)
    p = length(active)

    logEv = 0.5 *
        ( p * log(α)
        + N * log(β)
        - (β * rtr + α * normSquared(w))
        - logdetH(α, β, Λ)
        - N * log2π
        )
    return logEv
end

export effectiveNumParameters

function effectiveNumParameters(m::BayesianLinReg)
    α_over_β = m.priorPrecision / m.noisePrecision

    Λ = m.eig.values

    return sum(λ -> λ / (α_over_β + λ), Λ)
end

export posteriorPrecision

function posteriorPrecision(m::BayesianLinReg)
    return hessian(m)
end

export posteriorVariance

posteriorVariance(m::BayesianLinReg) = hessianinv!(m)

export posteriorWeights

function posteriorWeights(m)
    p = length(m.active)
    ϕ = posteriorPrecision(m)
    U = cholesky!(ϕ).U

    w = m.weights
    ε = inv(U) * view(m.uncertaintyBasis, m.active)
    return w + ε
end

# TODO: Why is this slower than `posteriorWeights`?
function postWeights(m)
    α = m.priorPrecision
    β = m.noisePrecision
    Λ = m.Λ
    Vt = m.Vt
    V = Vt'

    w = view(m.weights, m.active)

    # S = V * diagm(sqrt.(inv.(α .+ β .* Λ))) * Vt
    # ε = S * view(m.uncertaintyBasis, m.active)

    ε = view(m.uncertaintyBasis, m.active)
    ε .= @~ Vt * ε
    ε ./= sqrt.(α .+ β .* Λ)
    ε .= @~ V * ε

    return w + ε
end

export predict

function predict(m::BayesianLinReg, X; uncertainty=true, noise=false)
    noise &= uncertainty
    # dispatch to avoid type instability
    return _predict(m, X, Val(uncertainty), Val(noise))
end

function _predict(m, X, ::Val{true}, ::Val{true})
    yhat = _predict(m, X, Val(true), Val(false))
    n = length(yhat)
    noise = zeros(n) .± noiseScale(m)
    yhat .+= noise
    return yhat
end

_predict(m, X, ::Val{true}, ::Val{false}) = X * posteriorWeights(m)

_predict(m, X, ::Val{false}, ::Val{true}) = @error "Noise requires uncertainty"

_predict(m, X, ::Val{false}, ::Val{false}) = X * m.weights

export priorPrecision
priorPrecision(m::BayesianLinReg) = m.priorPrecision

export priorVariance
priorVariance(m::BayesianLinReg) = 1/m.priorPrecision

export priorScale
priorScale(m::BayesianLinReg) = sqrt(priorVariance(m))

export noisePrecision
noisePrecision(m::BayesianLinReg) = m.noisePrecision

export noiseVariance
noiseVariance(m::BayesianLinReg) = 1/m.noisePrecision

export noiseScale
noiseScale(m::BayesianLinReg) = sqrt(noiseVariance(m))

function Base.show(io::IO, m::BayesianLinReg{T}) where {T}
    @printf io "BayesianLinReg model\n"
    @printf io "\n"
    @printf io "Log evidence: %3.2f\n" logEvidence(m)
    @printf io "Prior scale: %5.2f\n" priorScale(m)
    @printf io "Noise scale: %5.2f\n" noiseScale(m)
    @printf io "\n"
    @printf io "Coefficients:\n"
    weights = posteriorWeights(m)
    for (j,w) in zip(m.active, weights)
        @printf io "%3d: %5.2f ± %4.2f\n" j w.val w.err
    end
end

###########################################################################################
# Helpers

function normSquared(v::AbstractVector{T}) where {T}
    s = zero(T)
    @inbounds @simd for i ∈ eachindex(v)
        s += v[i]*v[i]
    end
    return s
end

mydot(x,y) = mydot(promote(x,y...)...)

function mydot(x::AbstractVector{T},y::AbstractVector{T}) where {T}
    s = zero(T)
    @inbounds @simd for i ∈ eachindex(x)
        s += x[i]*y[i]
    end
    return s
end

export update!
function update!(m)
    ps = m.active
    XtX = m.XtX = view(m.XtX.parent, ps, ps)
    Xty = m.Xty = view(m.Xty.parent, ps)
    weights = m.weights = view(m.weights.parent, ps)
    Hinv = m.Hinv = view(m.Hinv.parent, ps, ps)

    m.eig = eigen(collect(m.XtX))
    m.eig.values .= max.(m.eig.values, 0.0)

    α = m.priorPrecision = 1.0
    β = m.noisePrecision = 1.0

    Q = m.eig.vectors
    Λ = m.eig.values
    D = Diagonal(inv.(α .+ β .* Λ))
    Hinv .= @~ Q * D * Q'

    m.weights .= @~ β .* (Hinv * Xty)
    return m
end

end # module