text:
#=
sensemap-sim.jl
Simulate coil sensitivity maps

Matlab notes:
- 2005-6-20, Jeff Fessler and Amanda Funai, University of Michigan
- 2014-08-19 JF more testing, verifying phase is correct
- 2014-09-09 modified for 3D by Mai Le
- 2016-05-03 JF fixes

2019-06-18, Jeff Fessler, University of Michigan
=#

export ir_mri_sensemap_sim

using SpecialFunctions: ellipk, ellipe

#using MIRT: ndgrid

"""
    smap = ir_mri_sensemap_sim(...)

Simulate 2D or 3D sensitivity maps for sensitivity-encoded MRI based on
grivich:00:tmf http://doi.org/10.1119/1.19461

This code makes maps for multiple coils,
but does not model coupling between coils,
so most likely it is an approximation at best.

option
- `dims::Dims`		image size; default (64, 64)
- `dx::Real`		pixel/voxel dimension; default: 3
- `dy::Real`		pixel/voxel dimension; default: `dx`
- `dz::Real`		""
- `ncoil::Int`		# of coils total; default 4
- `nring::Int`		# of rings of coils; default 1
- `rcoil::Real`		coil radius; default `dx * nx / 2 * 0.50`
- `dz_coil`			ring spacing in z; default `nz*dz/nring`
   * (3D geometry is a cylinder)
- `coil_distance::Real`		distance of coil center from isocenter
   * for central ring of coils as a multiple of `FOVx`,
   * where `FOVx=nx*dx`; default 1.2
- `orbit::Real`			default 360 [degrees]
- `orbit_start::AbstractVector{<:Real} = fill(0, nring)` [degrees]
- `scale::Symbol`
  + `:none` (default)
  +	`ssos_center` make SSoS of center = 1

out
- `smap	[dims ncoil]`	simulated sensitivity maps (complex!)

All length parameters must have same units (e.g., mm or cm)
"""
function ir_mri_sensemap_sim( ; kwargs...)
	return ir_mri_sensemap_sim(
		Vector{Tuple{Int,Int}}(undef, 0) ; # ir_ic_pair
		kwargs...,
	)[1]
end

"""
    (smap,info) = ir_mri_sensemap_sim( :all ; kwargs)
Like `ir_mri_sensemap_sim` but also returns `info` with data for all coils,
mainly for testing and plotting.
"""
function ir_mri_sensemap_sim(all::Symbol ;
	nring::Int = 1, ncoil::Int = 4, kwargs...,
)
	all == :all || throw("bad symbol argument $all")
	return ir_mri_sensemap_sim(
		vec(collect(Iterators.product(1:nring,1:ncoil))) ; # ir_ic_pair
		nring=nring, ncoil=ncoil, kwargs...,
	)
end

# handle (usual) scalar orbit_start
#=
(no: causes Method overwritten warning)
function ir_mri_sensemap_sim( ir_ic_pair::Vector{Tuple{Int,Int}} ;
	nring::Int = 1, orbit_start::Real = 0, kwargs...,
)
	return ir_mri_sensemap_sim( ir_ic_pair ; nring = nring,
		orbit_start = repeat([orbit_start], nring),
	)
end
=#

"""
    (smap,info) = ir_mri_sensemap_sim( ir_ic_pair ; kwargs)
Like `ir_mri_sensemap_sim` but also returns `info` with data for specific coils
where `ir_ic_pair::Vector{Tuple{Int,Int}}`.
(Usually used internally only.)
- `info::NamedTuple` geometry information for plots
"""
function ir_mri_sensemap_sim(
	ir_ic_pair::Vector{Tuple{Int,Int}} ; # usually length 0
	dims::Dims = (64,64), # 2D default
	dx::Real = 3,
	dy::Real = dx,
	dz::Real = 3,
	ncoil::Int = 4,
	nring::Int = 1,
	rcoil::Real = dx * dims[1] / 2 * 0.5,
	dz_coil::Real = ((length(dims) == 3) ? dims[3] : 1) * dz / nring,
	coil_distance::Real = 1.2, # multiplies fov/2
	orbit::Real = 360,
	orbit_start::AbstractVector{<:Real} = fill(0,nring),
	scale::Symbol = :none, # or :ssos_center
)

	(2 <= length(dims) <= 3) || throw("2D or 3D only")
	nx = dims[1]
	ny = dims[2]
	nz = (length(dims) == 3) ? dims[3] : 0

	coils_per_ring = round(Int, ncoil / nring)
	ncoil != nring * coils_per_ring && throw("nring must be divisor of ncoil")

	(smap, info) = ir_mri_sensemap_sim_do(
			nx, ny, nz,
			dx, dy, dz,
			ncoil, coils_per_ring, rcoil, dz_coil,
			orbit, orbit_start, coil_distance, ir_ic_pair)

	scale_center = (nz == 0) ?
		1 / sqrt(sum(abs2.(smap[round(Int,nx/2),round(Int,ny/2),:]))) :
		1 / sqrt(sum(abs2.(smap[round(Int,nx/2),round(Int,ny/2),round(Int,nz/2),:])))

	if scale === :ssos_center
		smap *= Float32(scale_center)
	elseif scale != :none
		throw("scale $scale")
	end

	return (smap, info)
end

"""
    (smap, info) = ir_mri_sensemap_sim_do()
"""
function ir_mri_sensemap_sim_do(
	nx, ny, nz,
	dx, dy, dz, ncoil, ncoilpr, rcoil, dz_coil,
	orbit, orbit_start, coil_distance, ir_ic_pair,
)

	T = Float32
	nring = Int(ncoil / ncoilpr)
	rlist = T(rcoil) * ones(T, ncoilpr, nring) # coil radii

	plist = zeros(T, ncoilpr, nring, 3) # position of coil center [x y z]
	nlist = zeros(T, ncoilpr, nring, 3) # normal vector (inward) from coil center
	olist = zeros(T, ncoilpr, nring, 3) # unit vector orthogonal to normal vector in x-y
	ulist = zeros(T, ncoilpr, nring, 3) # upward vector

	length(orbit_start) != nring && throw("bad orbit_start length")

	# cylindrical coil configuration, like abdominal coils
	alist = deg2rad.(orbit) * (0:(ncoilpr-1)) / ncoilpr # coil angles [radians]
	z_ring = ((1:nring) .- (nring+1)/2) * dz_coil
	for ir = 1:nring
		for ic = 1:ncoilpr
			phi = alist[ic] + deg2rad(orbit_start[ir])
			Rad = max(nx/2 * dx, ny/2 * dy) * coil_distance
			plist[ic,ir,:] = [Rad * [cos(phi), sin(phi)]; z_ring[ir]]
			nlist[ic,ir,:] = -[cos(phi), sin(phi), 0*z_ring[ir]] # cylinder
			olist[ic,ir,:] = [-sin(phi), cos(phi), 0]
			ulist[ic,ir,:] .= [0, 0, 1]
		end
	end

	# object coordinates
	x = T.(((1:nx) .- (nx+1)/2) * dx)
	y = T.(((1:ny) .- (ny+1)/2) * dy)
	z = (nz > 0) ? T.(((1:nz) .- (nz+1)/2) * dz) : [0]
	(xx, yy, zz) = ndgrid(x,y,z)

	data_per_coil = Array{Any}(undef, nring, ncoilpr)

	smap = zeros(Complex{T}, nx, ny, max(nz,1), ncoilpr, nring)
	for ir = 1:nring
		for ic=1:ncoilpr
			# rotate coordinates to correspond to coil orientation
			zr =	(xx .- plist[ic,ir,1]) .* nlist[ic,ir,1] +
					(yy .- plist[ic,ir,2]) .* nlist[ic,ir,2] +
					(zz .- plist[ic,ir,3]) .* nlist[ic,ir,3]
			xr = xx .* nlist[ic,ir,2] - yy .* nlist[ic,ir,1]
			yr = zz .- plist[ic,ir,3] # translate along object z axis

			# compute sensitivity vectors in coil coordinates
			tmp = ir_mri_smap1.(xr, yr, zr, rlist[ic,ir])
			sx = [p[1] for p in tmp]
			sy = [p[2] for p in tmp]
			sz = [p[3] for p in tmp]

			# coil response depends on transverse magnetization only?
			# todo: unsure if this should depend on sy and ulist in 3D
				bx = sz * nlist[ic,ir,1] + sx * olist[ic,ir,1]
				by = sz * nlist[ic,ir,2] + sx * olist[ic,ir,2]
			#	bz = sz * nlist[ic,ir,3] + sx * olist[ic,ir,3]
			smap[:,:,:,ic,ir] = complex.(bx, by)

			if (ir,ic) ∈ ir_ic_pair # save data for plotting
				data_per_coil[ir, ic] =
					(xr=xr, yr=yr, zr=zr, sx=sx, sy=sy, sz=sz, bx=bx, by=by)
			end
		end
	end

	smap *= rlist[1] / T(2π) # trick: scale so maximum is near unity
	smap = nz == 0 ?
		reshape(smap, nx, ny, ncoil) :
		reshape(smap, nx, ny, nz, ncoil)

	info = (x=x, y=y, z=z, dx=dx, dy=dy, dz=dz,
		nlist=nlist, plist=plist, rlist=rlist, olist=olist, ulist=ulist,
		nring=nring, ncoilpr=ncoilpr, rcoil=rcoil, data=data_per_coil)
	return smap, info
end

"""
    ir_mri_smap_r(r, z)
function for testing near 0
"""
function ir_mri_smap_r(r, z)
	M = 4 * r / ((1 + r)^2 + z^2) # = k^2, see ellipke
	(K,E) = (ellipk.(M), ellipe.(M))
	return 2 * z / r * ((1 + r)^2 + z^2)^(-0.5) *
		((1 + r^2 + z^2) / ((1 - r)^2 + z^2) * E - K)
end

"""
    ir_mri_smap1()
based on grivich:00:tmf

for a circular coil in "x-y plane" of radius "a"

Note that coil x-y plane is not same as object x-y plane!

Returns `(i,j,k)` components of B vector for each `(x,y,z)` location
"""
function ir_mri_smap1(x, y, z, a)
	x = x / a # normalized units
	y = y / a
	z = z / a
	r = sqrt(x^2 + y^2)
	M = 4 * r / ((1 + r)^2 + z^2) # = k^2, see ellipke
	(K,E) = (ellipk.(M), ellipe.(M))

	# the following is B_z in eqn (18) in grivich:00:tmf
	# and same as eqn [10] in wang:00:dop to within constant scale factor
	smap_z = 2 * ((1 + r)^2 + z^2)^(-0.5) *
		(K + (1 - r^2 - z^2) / ((1 - r)^2 + z^2) * E)
	smap_z /= a

	# the following is B_r in eqn (17) in grivich:00:tmf
	smap_r = 2 * z / r * ((1+r)^2 + z^2)^(-0.5) *
		((1 + r^2 + z^2) / ((1-r)^2 + z^2) * E - K)

	if abs(r) < 1e-6
		smap_r = 3 * pi * z / ((1 + z^2).^2.5) * r
	end
	smap_r /= a

	(isnan(smap_r) || isnan(smap_z)) && throw("nan")

	smap_x = smap_r * (r == 0 ? 0 : x / r)
	smap_y = smap_r * (r == 0 ? 0 : y / r)

#	phi = atan2(y, x)
#	smap_x = smap_r * cos(phi)
#	smap_y = smap_r * sin(phi)

	return Float32(smap_x), Float32(smap_y), Float32(smap_z)
end


text:
# Julia translation of
# https://stackoverflow.com/questions/16992038/inline-labels-in-matplotlib

export label_line, label_lines

using PyCall

"""
Places a text label next to the end of a line in a line plot.
"""
function label_line(line, x, label=none, align=true; xshift=0, yshift=0, extra_args...)
    ax = line["axes"]
    xdata = line[:get_xdata]()
    ydata = line[:get_ydata]()

    @assert x >= xdata[1] && x <= xdata[end] "x label location is outside data range!"

    ip = 1
    for i in 1:length(xdata)
        if x < xdata[i]
            ip = i
            break
        end
    end

    y = (ydata[ip + 1] + ydata[ip]) / 2
    xdist = xdata[ip + 1] - xdata[ip]
    if xdist != 0
        y = ydata[ip] + (ydata[ip + 1] - ydata[ip]) * (x - xdata[ip]) / xdist
    end

    if length(label) == 0
        label = line[:get_label]()
    end

    transformed_degrees = 0
    if align
        # Compute the slope
        dx = xdata[ip + 1] - xdata[ip]
        dy = ydata[ip + 1] - ydata[ip]
        degrees::Float64 = atan(dy, dx) |> rad2deg

        # Transform to screen coordinates.
        # Use reshape() here to avoid getting out a RowVector, which isn't
        # currently converted into an ndarray by PyCall.jl.
        pt = reshape([x, y], (1, 2))
        transformed_degrees = ax["transData"][:transform_angles]([degrees], pt)[1]
    end

    # Set a bunch of keyword arguments
    kwargs = Dict{Symbol, Any}(extra_args)
    if !haskey(kwargs, :color)
        kwargs[:color] = line[:get_color]()
    end

    if !haskey(kwargs, :alpha)
        kwargs[:alpha] = line[:get_alpha]()
    end

    if !haskey(kwargs, :horizontalalignment) && !haskey(kwargs, :ha)
        kwargs[:ha] = "center"
    end

    if !haskey(kwargs, :verticalalignment) && !haskey(kwargs, :va)
        kwargs[:va] = "center"
    end

    if !haskey(kwargs, :backgroundcolor)
        kwargs[:backgroundcolor] = ax[:get_facecolor]()
    end

    if !haskey(kwargs, :clip_on)
        kwargs[:clip_on] = true
    end

    if !haskey(kwargs, :zorder)
        kwargs[:zorder] = 2.5
    end

    kwargs[:rotation] = transformed_degrees

    ax.text(x + xshift, y + yshift, label; kwargs...)
end

"""
Annotates the given lines with a text label adjacent to their end points (none
if label property not set).
"""
function label_lines(lines::Vector, align=true, xvals=nothing; kwargs...)
    ax = lines[1].axes
    lab_lines = []
    labels = []

    #Take only the lines which have labels other than the default ones
    for line in lines
        label = line.get_label()
        if !contains(label, "_line")
            push!(lab_lines, line)
            push!(labels, label)
        end
    end

    if xvals == nothing
        xmin, xmax = ax.get_xlim()
        xvals = range(xmin, xmax, length=(length(lab_lines) + 2))[2:end-1]
    end

    for (line, x, label) in zip(lab_lines, xvals, labels)
        xrange = extrema(line.get_xdata())
        label_line(line, clamp(x, xrange...), label, align; kwargs...)
    end
end

ext: jl


text:
@info "IncrementalInference is adding Flux related functionality."

# the factor definitions
export FluxModelsDistribution
export MixtureFluxModels

# Required packages
using .Flux
using DataStructures: OrderedDict
using Random, Statistics

# import Base: convert
import Random: rand

const _IIFListTypes = Union{<:AbstractVector, <:Tuple, <:NTuple, <:NamedTuple}

function rand(nfb::FluxModelsDistribution,
              N::Int=1 )
  #

  # number of predictors to choose from, and choose random subset
  numModels = length(nfb.models)
  allPreds = 1:numModels |> collect # 1:Npreds |> collect
  # TODO -- compensate when there arent enough prediction models
  if numModels < N
    reps = (N ÷ numModels) + 1
    allPreds = repeat(allPreds, reps  )
    resize!(allPreds, N)
  end
  # samples for the order in which to use models, dont shuffle if N models
  # can suppress shuffle for NN training purposes
  1 < numModels && nfb.shuffle[] ? shuffle!(allPreds) : nothing

  # generate the measurements
  meas = zeros(nfb.outputDim..., N)
  for i in 1:N
    meas[:,i] = (nfb.models[allPreds[i]])(nfb.data)
  end

  return meas
end

FluxModelsDistribution( inDim::NTuple{ID,Int},
                        outDim::NTuple{OD,Int},
                        models::Vector{P},
                        data::D,
                        shuffle::Bool=true,
                        serializeHollow::Bool=false ) where {ID,OD,P,D<:AbstractArray} = FluxModelsDistribution{ID,OD,P,D}(inDim, outDim, models, data, Ref(shuffle), Ref(serializeHollow) )
#

FluxModelsDistribution( models::Vector{P},
                        inDim::NTuple{ID,Int},
                        data::D,
                        outDim::NTuple{OD,Int};
                        shuffle::Bool=true,
                        serializeHollow::Bool=false ) where {ID,OD,P,D<:AbstractArray} = FluxModelsDistribution{ID,OD,P,D}(inDim, outDim, models, data, Ref(shuffle), Ref(serializeHollow) )
#

"""
    $SIGNATURES

Helper function to construct `MixtureFluxModels` containing a `NamedTuple`, resulting in a
`::Mixture` such that `(fluxnn=FluxNNModels, c1=>MvNormal, c2=>Uniform...)` and order sensitive
`diversity=[0.7;0.2;0.1]`.  The result is the mixture heavily favors `.fluxnn` and names
`c1` and `c2` for two other components were auto generated.

Notes
- The user can specify own component names if desired (see example).
- `shuffle` is passed through to internal `FluxModelsDistribution` to command shuffling of NN models.
  - `shuffle` does not influence selection of components in the mixture.

Example:

```julia
# some made up data
data = randn(10)
# Flux models
models = [Flux.Chain(softmax, Dense(10,5,σ), Dense(5,1, tanh)) for i in 1:20]
# mixture with user defined names (optional) -- could also just pass Vector or Tuple of components
mix = MixtureFluxModels(PriorSphere1, models, (10,), data, (1,),
                        (naiveNorm=Normal(),naiveUnif=Uniform()),
                        [0.7; 0.2; 0.1],
                        shuffle=false )
#

# test by add to simple graph
fg = initfg()
addVariable!(fg, :testmix, Sphere1)
addFactor!(fg, [:testmix;], mix)

# look at proposal distribution from the only factor on :testmix
_,pts,__, = localProduct(fg, :testmix)
```

Related

Mixture, FluxModelsDistribution
"""
function MixtureFluxModels( F_::FunctorInferenceType,
                            nnModels::Vector{P},
                            inDim::NTuple{ID,Int},
                            data::D,
                            outDim::NTuple{OD,Int},
                            otherComp::_IIFListTypes,
                            diversity::Union{<:AbstractVector, <:NTuple, <:DiscreteNonParametric};
                            shuffle::Bool=true,
                            serializeHollow::Bool=false ) where {P,ID,D<:AbstractArray,OD}
  #
  # must preserve order
  allComp = OrderedDict{Symbol, Any}()

  # always add the Flux model first
  allComp[:fluxnn] = FluxModelsDistribution(nnModels,
                                            inDim,
                                            data,
                                            outDim,
                                            shuffle=shuffle,
                                            serializeHollow=serializeHollow)
  #
  isNT = otherComp isa NamedTuple
  for idx in 1:length(otherComp)
    nm = isNT ? keys(otherComp)[idx] : Symbol("c$(idx+1)")
    allComp[nm] = otherComp[idx]
  end
  # convert to named tuple
  ntup = (;allComp...)

  # construct all the internal objects
  return Mixture(F_, ntup, diversity)
end

MixtureFluxModels(::Type{F},
                  w...;
                  kw...) where F <: FunctorInferenceType = MixtureFluxModels(F(LinearAlgebra.I),w...;kw...)

#


text:
#
#--------------------------------------#
export linsolve
#--------------------------------------#
function linsolve(a,problem,solver)
    lhs,rhs = problem(a)
    return solver(lhs,rhs,false)
end
#--------------------------------------#
Zygote.@adjoint function linsolve(a,problem,solver)
    lhs,rhs = problem(a)
    u = solver(lhs,rhs,false)
    function fun(u̅)
        λ = solver(lhs,u̅,true)
        _,drda=pullback((aa)->resi(u,aa,problem),a)
        return (drda(λ)[1],nothing)
    end
    return u,fun
end
#--------------------------------------#
function resi(u,a,problem)
    lhs, rhs = problem(a)
    return rhs .- lhs(u);
end
#--------------------------------------#


text:
module BayesNegativeBinomial

using Distributions
using PolyaGammaSamplers
using LinearAlgebra
using Random

import Distributions: pdf, logpdf

include("body.jl")

end


text:
mutable struct STL{T<:TimeArray}
    ta::T         # A time array with three time series from a fitted STL model
    call::String  # method called to generate ta
end

"""
Package: Forecast

    stl(Yv, np; robust=false,
                nl=nextodd(np),
                ns=10*length(Yv)+1,
                nt=nextodd(1.5*np/(1-1.5/ns)),
                ni=robust ? 1 : 2,
                no=0,
                spm=false,
                qsmp=max(div(np,7),2))

Decompose a time series into trend, seasonal, and remainder components.

"STL has a simple design that consists of a sequence of applications of the loess smoother; the simplicity allows analysis of the properties of the procedure and allows fast computation, even for very long time series and large amounts of trend and seasonal smoothing. Other features of STL  are specification of amounts of seasonal and trend smoothing that range, in a nearly continous way, from very small amount of smoothing to a very large amount; robust estimates of the trend and seasonal components that are not distorted by aberrant behavior in the data; specification of the period of the seasonal component to any intenger multiple of the time sampling interval greater than one; and the ability to decompose time series with missing values."*

All default values are chosen following the recommendations of the original paper when those were recommended. `ns` is recommended to be chosen of the basis of knowledge of the time series and on the basis of diagnostic methods; it must nonethelessbe  always odd and at least 7. A default value is not advised on the original paper, instead the same default value used in the stl implementation in R in usere here.

for `no` the authors advise 5 ("safe value") or 10 ("near certainty of convergence") cycles  or a convergence criterion when robustness is required, in this case when `robust` is true computations stop when convergence is achieved in trend and seasonality.

for `qsmp` the authors do not adivise a default but they use a value close to div(`np`,7).

# Arguments
- `np`: Seasonality.
- `robust`: Robust stl.
- `nl`: Smoothing parameter of the low-pass filter.
- `ns`: Smoothing parameter for the seasonal component.
- `nt`: Smoothing parameter for the trend decomposition.
- `ni`: Number of inner loop cycles.
- `no`: Number of outer loop cycles.
- `spm`: Seasonal post-smoothing.
- `qsmp`: Loess q window for Seasonal post-smoothing.
- `verbose`: If true shows updates for the Seasonal and Trend convergence.
- `cth`: Corvengence threshold for Seasonal and Trend.

# Returns
An `STL` object with the seasonal, trend and remainder components if Yv is an Array and a TimeSeries object with the same components if Yv is a TimeSeries.

* STL: A Seasonal, Trend Decomposition Procedure Based on Loess" Robert B. Cleveland, William S. Cleveland, Jean E. McRae, and Irma Terpenning. Journal of Official Statistics Vol. 6. No. 1, 1990, pp. 3-73 (c) Statistics Sweden.

# Examples
```julia-repl
julia> stl(co2(),365; robust=true, spm=true)
[ Info: Dataset used in Cleveland et al. paper
[ Info: Corvengence achieved (< 0.01); Stopping computation...
STL{TimeSeries.TimeArray{Union{Missing, Float64},2,
Dates.Date,Array{Union{Missing, Float64},2}}}
(4609×3 TimeSeries.TimeArray{Union{Missing, Float64},2,
Dates.Date,Array{Union{Missing, Float64},2}} 1974-05-17 to 1986-12-31,
"stl(Yn, np=365; nl=365, ns=46091, nt=549, ni=1, no=0, spm=true, qsmp=52)")
```
"""
function stl(Yv::TimeArray,np::Integer;
             robust=false,
             nl=nextodd(np),
             ns=10*length(Yv)+1,
             nt=nextodd(1.5*np/(1-1.5/ns)),
             ni=robust ? 1 : 2,
             no=0,
             spm=false,
             qsmp=max(div(np,7),2),
             verbose=false,
             cth = 0.01)

    stl_obj = stl(TimeSeries.values(Yv),np;
              robust=robust, nl=nl, ns=ns, nt=nt,
              ni=ni, no=no, spm=spm, qsmp=qsmp, verbose=verbose)

    STL(TimeArray(stl_obj.ta; timestamp = timestamp(Yv)), stl_obj.call)

end

function stl(Yv,
             np::Integer;
             robust=false,
             nl=nextodd(np),
             ns=10*length(Yv)+1,
             nt=nextodd(1.5*np/(1-1.5/ns)),
             ni=robust ? 1 : 2,
             no=0,
             spm=false,
             qsmp=max(div(np,7),2),
             verbose = false,
             cth = 0.01) #where T<:Union{Missing, Number}

    @assert mod(ns,2)==1 & (ns>=7) "`ns` is chosen of the basis of knowledge of the time series and on the basis of diagnostic methods; must always be odd and at least 7"

    function B(u)
        (u < 1) & !ismissing(u) ? (1.0-u^2)^2 : 0.0
    end

    N = length(Yv)
    # initial robustness weights
    rhov = ones(N)
    # intial trend
    Tv = Tv0 = zeros(N)
    Sv = Sv0 = zeros(N)
    Rv = Array{Float64,1}(undef,N)
    Cv = Array{Float64,1}(undef,N+2*np)
    scnv = false # seasonal convergence flag
    tcnv = false # trend convergence flag
    #for o in 0:no
    o = 0
    while robust | (o <= no)
        for k in 1:ni
            # Updating sesonal and trend components
            ## 1. Detrending (Yv = Tv + Sv)
            Sv = Yv - Tv

            ### Seasonal convergence criterion
            Md = maximum(abs.(skipmissing(Sv-Sv0)))
            M0 = maximum(skipmissing(Sv0))
            m0 = minimum(skipmissing(Sv0))
            scnv = (Md/(M0-m0) < cth)
            if verbose
                println("Outer loop: " * string(o) * " - " * "Inner loop: " * string(k))
                println("Seasonal Convergence: " * string(Md/(M0-m0)))
            end
            Sv0 = Sv

            # Seasonal Smoothing 2,3,4
            ## 2. Cycle-subseries Smoothing
            for csi in 1:np
                Cv[csi:np:N+2*np] = loess(csi:np:N,
                                          Sv[csi:np:N];
                                          q=ns,d=1,rho=rhov[csi:np:N],
                                          predict=(1.0*csi-np):np:(N+np))
            end
            ## 3. Low-Pass Filtering of Smoothed Cycle-Subseries
            ### centered support instead 1:N to balance out machine error
            Lv = loess(1-ceil(N/2):N-ceil(N/2),
                       collect(skipmissing(sma(sma(sma(Cv,np),np),3))),
                       d=1,q=nl,rho=rhov)
            ## 4. Detreending of Smoothed Cycle-Subseries
            ### Lv is substracted to prevent low-frenquency power
            ### from entering the seasonal component.
            Sv = Cv[np+1:end-np] - Lv

            ## 5. Deseasonalizing
            Dv = Yv - Sv

            # Trend Smoothing
            ## 6. Trend Smoothing
            ### centered support instead 1:N to balance out machine error
            ### (floor isntead ceil like in Lv to balance out even lengths)
            Tv = loess(1-floor(N/2):N-floor(N/2),Dv,q=nt,d=1,rho=rhov)

            ### Trend convergence criterion
            Md = maximum(abs.(skipmissing(Tv-Tv0)))
            M0 = maximum(skipmissing(Tv0))
            m0 = minimum(skipmissing(Tv0))
            tcnv = (Md/(M0-m0) < cth)
            if verbose
                println("Trend    Convergence: " * string(Md/(M0-m0)) * "\n")
            end
            Tv0 = Tv

            scnv & tcnv ? break : nothing
        end
        # Computation of robustness weights
        ## These new weights will reduce the influence of transient behaviour
        ## on the trend and seasonal components.

        # Tv and Sv defined everywhere but Rv is not defined
        # where Yv has missing values

        Rv = Yv - Tv - Sv

        if scnv & tcnv
            @info "Corvengence achieved (< " * string(cth) * "); Stopping computation..."
            break
        end

        if 0 < o <= no
            smRv = skipmissing(Rv)
            h = 6*median(abs.(smRv))
            rhov = B.(abs.(Rv)/h)
        end
        o += 1
    end

    if spm
        # Using div(np,7) as default to approximate
        # the q=51 for a np=365 chosen in the original paper
        Sv = loess(1-ceil(N/2):N-ceil(N/2),Sv,q=qsmp)
        Rv = Yv - Tv - Sv
    end

    if !scnv
        @warn "Seasonal convergence not achieved (>= " * string(cth) * ")
         Consider a robust estimation"
    end

    if !tcnv
        @warn "Trend convergence not achieved (>= " * string(cth) * ")
         Consider a robust estimation"
    end

    call = "stl(Yn, np="*string(np)*
        "; nl="*string(nl)*
        ", ns="*string(ns)*
        ", nt="*string(nt)*
        ", ni="*string(ni)*
        ", no="*string(no)*
        ", spm="*string(spm)*
        ", qsmp="*string(qsmp)*")"

    # Dates.datetime2unix(DateTime(1970,1,1,0,0,1)) = 1.0
    stl_ts = (Timestamp = DateTime(1970,1,1,0,0,1):Second(1):DateTime(1970,1,1,0,0,1)+Second(N-1),
              Seasonal = Sv,
              Trend = Tv,
              Remainder = Rv)

    STL(TimeArray(stl_ts; timestamp = :Timestamp), call)
end


text:
for i = (x = 16; 0) : 0 end
println(x + 26)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using ModelingToolkit, Catalyst, Symbolics, Statistics, DataFrames, DiffEqJump, Plots
function create_neighbour(N;dims=1)
    dims_=Tuple(fill(N,dims))
    dict = Dict{CartesianIndex{dims},Vector{CartesianIndex{dims}}}()
    cart = CartesianIndices(dims_)
    for ind in cart
        a = Tuple(ind)
        neighbour = []
        for shift in [1,-1]
            for i in 1:dims
                b = Int.(zeros(dims))
                b[i] = shift
                c = CartesianIndex(a .+ Tuple(b))
                if c in cart
                    push!(neighbour,c)
                end
            end
        end
        push!(dict,ind=>neighbour)
    end
    return dict
end

N = 128
dict_neighbour = @time create_neighbour(N,dims=1)
@parameters  D t
X = [Num(Symbolics.variable(:X,i; T = ModelingToolkit.FnType{Tuple{Any},Real}))(t) for i in 1:N]
begin # construct reaction network
    rxs_tran = []
    for ind in CartesianIndices(X)
        for ind_neig in dict_neighbour[ind]
            push!(rxs_tran,Reaction(D,[X[ind]],[X[ind_neig]]))
        end
    end
    rxs = vcat(rxs_tran)
    @named rs = ReactionSystem(rxs,t,vec(X),[D])
end

h = 1 / N
u0 = 10*ones(Int64, N)
tf = .01
fmt = :png
methods = (Direct(),MNRM(),NRM(),DirectCR(),RSSA())
shortlabels = [string(leg)[12:end-2] for leg in methods]
prob    = DiscreteProblem(rs, u0, (0.0, tf), [1 / (h*h)])
ploth   = plot(reuse=false)
for (i,method) in enumerate(methods)
    println("Benchmarking method: ", method)
    jump_prob = JumpProblem(rs, prob, method, save_positions=(false, false))
    sol = solve(jump_prob, SSAStepper(), saveat=tf/1000.)
    plot!(ploth,sol.t,sol[Int(N//2),:],label=shortlabels[i], format=fmt)
end
plot!(ploth, title="Population at middle lattice site", xlabel="time",format=fmt)

function run_benchmark!(t, jump_prob, stepper)
    sol = solve(jump_prob, stepper)
    @inbounds for i in 1:length(t)
        t[i] = @elapsed (sol = solve(jump_prob, stepper))
    end
end

nsims = 50
benchmarks = Vector{Vector{Float64}}()
for method in methods
    jump_prob = JumpProblem(rs, prob, method, save_positions=(false, false))
    stepper = SSAStepper()
    t = Vector{Float64}(undef,nsims)
    println("Benchmarking method: ", method)
    run_benchmark!(t, jump_prob, stepper)
    push!(benchmarks, t)
end

medtimes = Vector{Float64}(undef,length(methods))
stdtimes = Vector{Float64}(undef,length(methods))
avgtimes = Vector{Float64}(undef,length(methods))
for i in 1:length(methods)
    medtimes[i] = median(benchmarks[i])
    avgtimes[i] = mean(benchmarks[i])
    stdtimes[i] = std(benchmarks[i])
end

df = DataFrame(names=shortlabels,medtimes=medtimes,relmedtimes=(medtimes/medtimes[1]),
                avgtimes=avgtimes, std=stdtimes, cv=stdtimes./avgtimes)

sa = [string(round(mt,digits=4),"s") for mt in df.medtimes]
bar(df.names,df.relmedtimes,legend=:false, fmt=fmt)
# df.relmedtimes
scatter!(df.names, .05 .+ df.relmedtimes, markeralpha=0, series_annotations=sa, fmt=fmt)
ylabel!("median relative to Direct")
title!("256 Site 1D Diffusion CTRW")


text:
module GameOfLife

using CircularArrays

@enum Cell::UInt8 Live=true Dead=false

using Base: show, convert

Base.convert(::Type{Bool}, cell::Cell) = cell == Live

function Base.show(io::IO, ::MIME"text/html", v::Cell)
	displayWidth, displayHeight = get(io, :displaysize, (3, 3))
	style = string("background: ", v == Live ? "black; " : "white; ")
	compact = displayWidth < 3 || displayHeight < 3
	if compact
		style = string(
			style,
			"height: $(displayHeight)px; width: $(displayWidth)px; "
		)
	else
		style = string(
			style,
			"height: $(displayHeight)px; width: $(displayWidth)px; ",
			"border: 1px solid ",
			v == Live ? "white; " : "black; "
		)
	end
	write(io, """""")
end

function Base.show(io::IO, mime::MIME"text/html", vec::AbstractVector{Cell})
	(displayWidth, displayHeight) = get(io, :displaysize, (100, 10))
	cellWidth = max(1, displayWidth ÷ size(vec, 1))
	displayHeight = min(displayHeight, cellWidth)
	cellIO = IOContext(io, :displaysize => (cellWidth, displayHeight))

	write(io, """
""")
	for cell in vec
		Base.show(cellIO, mime, cell)
	end
	write(io, """
""")
end

"""
	Handle displaying the cell grid in a HTML context. Does this by producing a SVG.
	Originally I did this as HTML (see the AbstractVector/Cell show functions above), but
	the layout kept crashing my browser at higher cell sizes.
"""
function Base.show(io::IO, mime::MIME"text/html", cell_grid::AbstractMatrix{Cell})
	displayHeight, displayWidth = get(io, :displaysize, (600,600))
	cellHeight = displayHeight / size(cell_grid, 1)
	cellWidth = displayHeight / size(cell_grid, 2)

	# Calculate a border size
	borderSize = 0.5
	if min(cellHeight / 10, cellWidth / 10) < borderSize
		borderSize = 0
	end

	write(io, """<svg width="$(displayWidth)" height="$(displayHeight)">

		""")

	# CartesianIndices are lazy, so this is efficient
	for i in CartesianIndices(cell_grid)
		@inbounds cell = cell_grid[i]
		# Iteration (and thus destructuring) of CartesianIndex is specifically not supported
		y, x = Tuple(i)
		write(io, """
			<rect x="$((x-1)*cellWidth)" y=$((y-1)*cellHeight) width="$(cellWidth)" height="$(cellHeight)" class="$(
			cell === GameOfLife.Live ? "live" : "dead"
			)" />
			""")
	end

	write(io, "</svg>")
end

"""
Generates a cell that's the updated state for cell_grid(2,2).

"""
function step_cell(cell_grid::AbstractMatrix{Cell})
	# cell_grid is a 1:3,1:3 per Julia 1 indexing. Remove the center cell (2,2)
	centerCellIndex = CartesianIndex(2,2)
	neighbour_indices = filter(
		!=(centerCellIndex),
		CartesianIndices(cell_grid)
	)
	neighbours = view(cell_grid, neighbour_indices)

	# Doesn't call convert(Bool)
	# num_live_neighbours = count(neighbours)
	num_live_neighbours = count(==(Live), neighbours)
	@inbounds if cell_grid[centerCellIndex] == Live
		if num_live_neighbours < 2
			return Dead
		# elseif num_live_neighbours == 2 || num_live_neighbours == 3
		elseif num_live_neighbours < 4
			return Live
		else
			return Dead
		end
	else
		if num_live_neighbours == 3
			return Live
		end
		return Dead
	end
end

"""
	step_grid(grid)

Produce the next iteration of the grid.

Using CircularArray to handle periodic overflow, because we're going to do more complicated things
"""
function step_grid(cell_grid::AbstractMatrix{Cell})
	out = similar(cell_grid)

	circular_cell_grid = CircularArray(cell_grid)

	range = CartesianIndices(cell_grid)
	it_first = first(range)

	range_one = oneunit(it_first)

	Threads.@threads for it in range
	#for it in range
		it_lower = it-range_one
		it_higher = it+range_one
		local_region = it_lower:it_higher

		region = view(circular_cell_grid,local_region)

		@inbounds out[it] = step_cell(region)
	end

	out
end

export step_grid, Cell

end # module


text:
using Bayesian_EMR_MSM
using Test

using Distributions
using LinearAlgebra
import Random

@testset "emr_msm_estimate" begin

    @testset "EMR_MSM_PointEstimate" begin
        Bayesian_EMR_MSM.cmdstan_home!("/cmdstan")
        delta = 1/50

        seed = 2
        rand_gen = Normal(0,5)
        num_obs    = 7
        num_params = 2
        num_layers = 3
        num_est_layers = 2
        num_ts = 2
        σ = 0.001
        Σ = Diagonal(repeat([σ],num_params))
        μ = zeros(num_params)
        num_pred = floor(Int64, delta*500)
        timestep = 0.0015/delta
        F_zero = false
        A_zero = false
        B_zero = false

        num_samples = 10
        num_chains = 2

        Random.seed!(seed)

        # actually only 1 needed
        x_start   = [rand(rand_gen, num_obs, num_params) for ts in 1:num_ts]
        res_start = [rand(rand_gen, num_obs, num_params, num_layers+1)/100
                        for ts in 1:num_ts]

        F = ifelse(F_zero,zeros(num_params),rand(rand_gen, num_params))
        A = ifelse(A_zero,zeros(num_params,num_params),rand(rand_gen, num_params, num_params))
        B = ifelse(B_zero,zeros(num_params,num_params,num_params),rand(rand_gen, num_params, num_params, num_params))

        for i in 1:num_params
            B[:,:,i] = 0.5*(B[:,:,i] + transpose(B[:,:,i]))
        end

        L = Bayesian_EMR_MSM.ResCorrs(rand(rand_gen,
        length(Bayesian_EMR_MSM.ResCorrs,
        num_params, num_layers)), num_params, num_layers)

        model = Bayesian_EMR_MSM.EMR_MSM_Model_PointEstimate(F, A, B, L, μ, Σ)

        timeseries_start = [Bayesian_EMR_MSM.MSM_Timeseries_Point{Float64}(
        x_start[ts],
        res_start[ts],
        repeat([1.0], num_obs)
        ) for ts in 1:num_ts]

        test_point_est = Bayesian_EMR_MSM.EMR_MSM_Estimate(model, timeseries_start,
            timeseries_start
        )

        Bayesian_EMR_MSM.write("point_est", test_point_est)

        pred = Bayesian_EMR_MSM.EMR_MSM_Prediction(test_point_est, num_pred, repeat([timestep],num_ts);
                use_last_layer=false)

        Bayesian_EMR_MSM.write("pred_point", pred)

        print(Bayesian_EMR_MSM.DataFrame(pred.pred_timeseries[1]))

        tau0 = 100.0
        num_pred_samples = 100

        dist_est = Bayesian_EMR_MSM.EMR_MSM_Estimate(pred.pred_timeseries, num_est_layers,
        num_samples, num_chains, tau0, 456)

        Bayesian_EMR_MSM.write("dist_est", dist_est)

        mean_est = Bayesian_EMR_MSM.EMR_MSM_Model_PointEstimate(
                    dist_est.model,
                    num_params,
                    num_est_layers,
                    mean)

        test_pred_ts_med = Bayesian_EMR_MSM.MSM_Timeseries_Point{Float64}(
            dist_est.pred_timeseries,
            x -> median(x))

        dist_pred = Bayesian_EMR_MSM.EMR_MSM_Prediction(
            dist_est, num_pred, repeat([timestep],num_ts), num_pred_samples,
            repeat([1],num_ts);
            use_last_layer=false
        )

        Bayesian_EMR_MSM.write("pred_dist", dist_pred)

        @test true == true
    end

end


text:
#################################################################
# maximum bearing method
#
# Assumes gp only has one sample (looks at meangains)
# Doesn't handle null obs... does it have to?
function bearing_max(gp::GainPattern)
	num_angles = length(gp.angles)
	temp_max = -Inf
	max_bearing = 0
	for i = 1:num_angles
		if validgain(gp.meangains[i])
			if gp.meangains[i] > temp_max
				temp_max = gp.meangains[i]
				max_bearing = gp.angles[i]
			end
		end
	end
	return max_bearing
end
bearing_max(filename::AbstractString) = bearing_max(GainPattern(filename))


text:
@testset "solvecubic.jl" begin
    mix = load(BrusilovskyEoSMixture; names=("methane",))
    solve_cubic = CubicEoS.solve_cubic

    @testset "Unstable 1-phase state" begin
        # definitely inside "bulk" of two-phase region
        @test !first(vt_stability(mix, [0.01], 1e-6, CubicEoS.GAS_CONSTANT_SI * 140.0))
    end

    @testset "Degenerate equation" begin
        @test first(vt_stability(mix, [0.0109], 1e-6, 1588.0623600672689))
        @test first(vt_stability(mix, [0.0109], 1e-6, CubicEoS.GAS_CONSTANT_SI * 190.6))
    end

    @testset "Triple root" begin
        @test solve_cubic(1, -3, 3, -1) == (1.0, 1.0, 1.0)
        @test solve_cubic(1, 3, 3, 1) == (-1.0, -1.0, -1.0)
    end

    @testset "One real root" begin
        @test let
            p, q, r = solve_cubic(1, -5, 1, -5)
            p ≈ 5.0 && isnan(q) && isnan(r)
        end
        @test let
            p, q, r = solve_cubic(1, -5, 9, -5)
            p ≈ 1.0 && isnan(q) && isnan(r)
        end
    end

    @testset "solve_cubic!" begin
        solve_cubic! = CubicEoS.solve_cubic!
        @test let
            coefs = 1, -5, 1, -5
            p, q, r = solve_cubic(coefs...)
            buf = Vector{Float64}(undef, 5)
            solve_cubic!(buf, coefs...)
            any(buf[1:3] .== (p, q, r))
        end
    end
end

max_stars_repo_name: vvpisarev/CubicEoS

max_issues_repo_name: vvpisarev/CubicEoS

max_forks_repo_name: vvpisarev/CubicEoS


text:
using Elf
using Test
using Serialization

@testset "Elf.jl" begin
    # Write your tests here.
    bin = deserialize("testbin")
    ehdr = Elf64_Ehdr(bin)
    @test iself(ehdr)
    @test elfclass(ehdr) == ELFCLASS64
    @test endian(ehdr) == ELFDATA2MSB
    @test sections(bin) == deserialize("sections")
    @test segments(bin) == deserialize("segments")
    @test symbols(bin) == deserialize("symbols")
end


text:
rand_uniform(T::DataType, lb,ub,dims...) = T.(rand(Uniform(T(lb),T(ub)),dims...))
rand_uniform(T::DataType, lb,ub) = rand_uniform(T, lb,ub,1)[1]
rand_uniform(lb,ub,dims...) = rand_uniform(Float32, lb,ub,dims...)
rand_uniform(lb,ub) = rand_uniform(Float32, lb,ub, 1)[1]
Zygote.@nograd rand_uniform
add_dim(x::Array{T, N}) where {T,N} = reshape(x, Val(N+1))

mysigm(x) = one(x) / (one(x) + exp(-x))
#mysigm(x) = NNlib.σ(x)
#ModelingToolkit.@register NNlib.σ

struct Mapper{W,B,F,P}
  W::W
  b::B
  σ::F
  p::P
  function Mapper(W, b, σ, p)
    new{typeof(W),typeof(b),typeof(σ),typeof(p)}(W, b, σ, p)
  end
end
function Mapper(n::Integer,σ::Function=identity; T=Float32, init=dims->ones(T,n), bias=dims->zeros(T,n))
  W = init(n)
  b = bias(n)
  p = vcat(W,b)
  Mapper(W, b, σ, p)
end
function (m::Mapper)(x::AbstractMatrix, p=m.p)
  Wl = length(m.W)
  W = @view p[1 : Wl]
  b = @view p[Wl + 1 : end]
  m.σ.(W .* x .+ b)
end
Base.show(io::IO, m::Mapper) = print(io, "Mapper(", m.W, ", ", m.σ, ")")
initial_params(m::Mapper) = m.p
paramlength(m::Mapper) = length(m.p)
Flux.functor(::Type{<:Mapper}, m) = (m.p,), re -> Mapper(m.W, m.b, m.σ, re...)
Flux.trainable(m::Mapper) = (m.p,)
function get_bounds(m::Mapper, T::DataType=eltype(m.p))
  lb = T.(vcat(fill(-100.1, length(m.W)), fill(-100.1, length(m.b))))
  ub = T.(vcat(fill(100.1, length(m.W)), fill(100.1, length(m.b))))
  lb, ub
end

paramlength(m::Flux.Dense) = length(m.weight) + (m.bias != Flux.Zeros() ? length(m.bias) : 0)
paramlength(m::Union{Flux.Chain,FastChain}) = sum([paramlength(l) for l in m.layers])
#paramlength(m::AbstractFluxLayer) = length(Flux.destructure(m)[1])

function get_bounds(l, T::DataType)
  p, _ = Flux.destructure(l)
  fill(T(-Inf), length(p)), fill(T(Inf), length(p))
end
function get_bounds(m::Union{Flux.Chain, FastChain}, T::DataType)
  lb = vcat([get_bounds(layer, T)[1] for layer in m.layers]...)
  ub = vcat([get_bounds(layer, T)[2] for layer in m.layers]...)
  lb, ub
end
function get_bounds(m::FastDense, T::DataType=nothing)
  lb = T.(vcat(fill(-100.1, m.out*m.in), fill(-100.1, m.out)))
  ub = T.(vcat(fill(100.1, m.out*m.in), fill(100.1, m.out)))
  if !m.bias
    lb = T.(vcat(fill(-100.1, m.out*m.in)))
    ub = T.(vcat(fill(100.1, m.out*m.in)))
  end
  lb, ub
end

function get_bounds(m::Flux.Dense, T=Float32)
  lb = T.(vcat(fill(-100.1, length(m.weight)), fill(-100.1, length(m.bias))))
  ub = T.(vcat(fill(100.1, length(m.weight)), fill(100.1, length(m.bias))))
  if m.bias == Flux.Zeros()
    lb = T.(vcat(fill(-100.1, length(m.weight))))
    ub = T.(vcat(fill(100.1, length(m.weight))))
  end
  lb, ub
end

function reset_state!(m::Union{Flux.Chain, FastChain}, p)
  start_idx = 1
  for l in m.layers
    pl = paramlength(l)
    p_layer = @view p[start_idx:start_idx+pl-1]
    reset_state!(l, p_layer)
    start_idx += pl
  end
end

reset_state!(m,p) = nothing


text:
#= ========== RocksDB ==================
RocksDb specific code.

Each DB should export the symbols in KVS.jl
=#
using RocksDB

struct KVSRocksDB <: AbstractKVS
    name::String     # Name of db
    handle           # Handle to db object
end
KVSRocksDB(name::String) = KVSRocksDB(name, rocksDB_open(name))

KVSRocksDB(name::String, path::String) = KVSRocksDB(name, rocksDB_open(name, path))

function rocksDB_open(name::String, db_dir::String)
    if stat(db_dir).inode == 0
        Base.mkdir(db_dir)
    end
    handle = RocksDB.open_db(db_dir * name, true)
    handle
end

rocksDB_open(name::String) = rocksDB_open(name, kvs_get_db_root())

function kvs_get_db_root()
    getindex(ENV, "HOME") * "/kinant/"
end

function kvs_get_db_root(fs_id)
    getindex(ENV, "HOME") * "/kinant/" * hex(fs_id, sizeof(fs_id))
end

function kvs_open(hdl::KVSRocksDB)
    hdl.handle = rocksDB_open(hdl.name)
    hld
end

function kvs_close(hdl::KVSRocksDB)
    RocksDB.close_db(hdl.handle)
    hdl.name = ""
    hdl.handle = nothing
end

function kvs_put(hdl::KVSRocksDB, key, value)
    RocksDB.db_put(hdl.handle, key, value)
end

function kvs_put_sync(hdl::KVSRocksDB, key, value)
    RocksDB.db_put_sync(hdl.handle, key, value)
end

function kvs_delete(hdl::KVSRocksDB, key)
    RocksDB.db_delete(hdl.handle, key)
end

function kvs_delete_range(hdl::KVSRocksDB, first, last)
    RocksDB.db_delete_range(hdl.handle, first, last)
end

function kvs_get(hdl::KVSRocksDB, key)
    RocksDB.db_get(hdl.handle, key)
end

"""
    kvs_get_many(db::KVSRocksDB, first, last, max::Int; inc_first=false)
Return *max* entries beginning from *first* and upto but not including *last*.
Optionally specify whether the first entry needs to be included.
"""
function kvs_get_many(db::KVSRocksDB, first, last, max::Integer; inc_first=true, raw_read=false)
    #first = "\u00000000000000000000000000000000"
    #last  = "\uffffffffffffffffffffffffffffffff"
    @info "db: $(db.name) $(first) $(last) $(inc_first) $(raw_read)"
    keys = Vector{Any}(undef, max)
    vals = Vector{Any}(undef, max)
    i::Int = 0
    for (k, v) in RocksDB.db_range(db.handle, first, last; raw=raw_read)
        #if k == last break end # Don't include last
        !inc_first && (k == first) && continue # Don't include first
        i += 1
        if i > max # Fill only max entries
            i -= 1
            break
        end
        keys[i] = k
        vals[i] = v
    end
    (keys, vals, i)
end

function kvs_write_batch(db::KVSRocksDB, keys, vals; raw_write=false)
    batch = RocksDB.create_write_batch()
    for i = 1:length(keys)
        RocksDB.batch_put(batch, keys[i], vals[i]; raw=raw_write)
    end
    RocksDB.write_batch(db.handle, batch)
end

function kvs_create_checkpoint(db::KVSRocksDB, path::String)
    RocksDB.db_create_checkpoint(db.handle, path)
end

"""
Update the stream in the db atomically. For Rocks Db this is just a
plain update as each process holds a lock on the db.
"""
function kvs_atomic_update(db::KVSRocksDB, cid, s, op::atomic_op)
    @assert db.name == "recent"
    if op == ATOMIC_INC
        # if s.state != S_OPEN throw(StreamStateException("State != S_OPEN")) end
        s.ref = s.ref + 1
        s.seq = s.seq + 1
    elseif op == ATOMIC_DEC
        s.ref = s.ref - 1
    end
    # TODO: for consul, need to check if someone else is commiting (S_COMMITTING)
    # the stream. Then this process should back off. The best way to do that
    # would be to have a state transition checker throw an exception if the
    # transition is not right.

    RocksDB.db_put(db.handle, cid, s)
    s
end

"""
Try to atomically update db. For Rocks Db this is just a
plain update as each process holds a lock on the db.
"""
function kvs_atomic_try_update(db::KVSRocksDB, cid, s)
    @assert db.name == "recent"
    RocksDB.db_put(db.handle, cid, s)
    true
end

"""
Warning! Extreme caution! Deletes db!
"""
function kvs_clobber!(db::KVSRocksDB)
    db_root::String = getindex(ENV, "HOME")
    db_path = db_root * "/thimble/" * db.name
    kvs_close(db)
    run(`mkdir -p $db_root/bak`)
    run(`mv $db_path $db_root/bak/`)
    run(`unlink $db_path`)
end

"""
    db_chunk_exists(fp)

Returns true if chunk exists in db.
"""
function chunk_exists(db::KVSRocksDB, fp)
    c = RocksDB.db_get(db.handle, fp)
    c == nothing && return false
    return true
end

"""
return true if chunk already exists, else put chunk and return false
"""
function kvs_put_chunk(db::KVSRocksDB, input::Vector{UInt8}, fp, first::Int64, last::Int64)
    if chunk_exists(db, fp)
        return true
    else
        RocksDB.db_put(db.handle, fp, input[first:last])
    end
end

# Dummy function on a single node system
function kvs_put_atomic(db::KVSRocksDB, key, value, index)
    kvs_put(db, key, value) # Ignore index
    true
end

# Dummy function on a single node system
function kvs_get_index(db::KVSRocksDB, key)
    v = kvs_get(db, key)
    return (v, rand(Int64)) # Return a random index
end

"""
Atomically increment the chunk ref counts for all chunks in array.
The incrementing of each chunk's ref is atomic.

For RocksDB increment is implemented by a get and put.
"""
function kvs_atomic_batch_increment(db::KVSRocksDB, chunks)
    @assert db.name == "chunk_ref"
    saved::UInt64 = 0        # Bytes saved for this extent
    index::UInt32 = 0
    exists::BitVector =  falses(length(chunks))
    for i in chunks
        ref::UInt64 = 0
        index += 1
        ret = RocksDB.db_get(db.handle, i.cksum) # Get ref
        (ret != nothing) && (ref = ret)

        # Compute return values
        if (ref > 0)
            exists[index] = true
            saved += i.range.stop - i.range.start
        end

        # Increment ref and save
        ref = ref + 1
        RocksDB.db_put(db.handle, i.cksum, ref) # Put ref
    end
    (exists, saved)
end

# helper functions


text:
"""
`module  JAC.TestFrames`
    ... a submodel of JAC that contains all methods for testing individual methods, modules, ....
"""
module TestFrames

    using Printf, JLD, JAC, ..AngularMomentum, ..Basics, ..Continuum, ..Defaults, ..ManyElectron, ..Nuclear, ..Radial, ..TableStrings, ..SpinAngular

    export testDummy

    function testPrint(sa::String, success::Bool)
        printTest, iostream = Defaults.getDefaults("test flag/stream")
        ok(succ) =  succ ? "[OK]" : "[Fail]"
        sb = sa * TableStrings.hBlank(110);   sb = sb[1:100] * ok(success);    println(iostream, sb)
        return( nothing )
    end

    """
    `TestFrames.testEvaluation_rcsfValues(; short::Bool=true)`
        ... tests on rcsf value.
    """
    function testEvaluation_rcsfValues(; short::Bool=true)
        success = true
        etal = 1.4928400545843579
        printTest, iostream = Defaults.getDefaults("test flag/stream")
        aT = SpinAngular.getTermNumber(3.5, 0, 4, 0)
        bT = SpinAngular.getTermNumber(3.5, 0.5, 1.5, 0)
        coeff = SpinAngular.completlyReducedCfpByIndices(aT, bT)
        if abs(coeff - etal) >  0.000000000000001
            success = false
            if printTest   info(iostream, "$coeff:   $coeff != $etal")   end
        end

        testPrint("testEvaluation_rcsfValues()::", success)
        return(success)
    end

    """
    `TestFrames.testEvaluation_angularCoefientsValues(; short::Bool=true)`
        ... tests on rcsf value.
    """
    function testEvaluation_angularCoefientsValues(; short::Bool=true)
        success = true
        etal = 1.000000000000000

        # First step: Generation of the relativistic CSFs
        configs = [Configuration("1s^2 2s"), Configuration("1s^2 2p")]
        relconfList = ConfigurationR[]
        for conf in configs
        wax = Basics.generate("configuration list: relativistic", conf)
        append!( relconfList, wax)
        end
        subshellList = Basics.generate("subshells: ordered list for relativistic configurations", relconfList)
        Defaults.setDefaults("relativistic subshell list", subshellList; printout=false)
        csfList = CsfR[]
        for relconf in relconfList
        newCsfs = Basics.generate("CSF list: from single ConfigurationR", relconf, subshellList)
        append!( csfList, newCsfs)
        end

        # Second step: Calculation of spin-angular coefficients for a none scalar one-particle operator"
        leftCsf = csfList[1]
        rightCsf1 = csfList[2]
        rightCsf2 = csfList[3]
        op = SpinAngular.OneParticleOperator(1, plus, true)
        coeff = SpinAngular.computeCoefficients(op, leftCsf, rightCsf1, subshellList)
        if abs(coeff[1].T - etal) >  0.000000000000001
            success = false
            if printTest   info(iostream, "$coeff:   $coeff != $etal")   end
        end
        coeff = SpinAngular.computeCoefficients(op, leftCsf, rightCsf2, subshellList)
        if abs(coeff[1].T - etal) >  0.000000000000001
            success = false
            if printTest   info(iostream, "$coeff:   $coeff != $etal")   end
        end

        testPrint("testEvaluation_angularCoefientsValues()::", success)
        return(success)
    end

end # module


text:
using Distributions, Random

export PathDistribution, PathPseudoUniform, PathCooccurrence

abstract type PathDistribution end

struct PathPseudoUniform <: PathDistribution
    vertex_set::UnitRange{Int}
    length_dist::DiscreteUnivariateDistribution
end

function Base.rand(d::PathPseudoUniform)
    m = rand(d.length_dist)
    return rand(d.vertex_set, m)
end

function Random.rand!(p::Vector{Int}, d::PathPseudoUniform)

    m = rand(d.length_dist)
    resize!(p, m)
    sample!(d.vertex_set, p)

end

struct PathConditionalCategorial <: PathDistribution
    α::Float64
    V::UnitRange
end

PathCondCat = PathConditionalCategorial

function get_cond_dist(
    S::InteractionSequence{Int};
    α::Float64=0.0,
    V::Int=length(unique(vcat(S...)))
    )
    tot_len = 0
    c = fill(α, V)
    for path in S
        for v in path
            c[v] += 1.0
            tot_len += 1
        end
    end
    c ./= (tot_len + α * V)
    return c
end

function Base.rand(d::PathCondCat, curr::InteractionSequence{Int})
    m = rand(Poisson(mean(length.(curr))))
    p = get_cond_dist(curr, α=d.α, V=length(d.V))
    return rand(Categorical(p), m)
end

function Random.rand!(x::Vector{Int}, d::PathCondCat, curr::InteractionSequence{Int})
    m = rand(Poisson(mean(length.(curr))))
    resize!(x, m)
    p = get_cond_dist(S, α=d.α, V=length(d.V))
    sample!(Categorical(p), x)
end

struct PathCooccurrence{T<:Union{Int,String}} <: PathDistribution
    vertex_set::Vector{T}
    length_dist::DiscreteUnivariateDistribution
    initial_dist::DiscreteUnivariateDistribution
    P::Matrix{Float64}
    P_dists::Dict{Int,Categorical}
    function PathCooccurrence(
        vertex_set::Vector{S},
        length_dist::DiscreteUnivariateDistribution,
        initial_dist::DiscreteUnivariateDistribution,
        P::Matrix{Float64}
        ) where {S<:Union{Int,String}}
        @assert length(vertex_set) == size(P)[1] "Dimension missmatch. Check vertex set and transition matrix dimensions agree."
        dists = Dict{Int, Categorical}()
        for i in 1:length(vertex_set)
            dists[i] = Categorical(P[:,i])
        end
        new{S}(vertex_set, length_dist, initial_dist, P, dists)
    end
end

# When no vertex set it given we assume integer vertices
function PathCooccurrence(
    length_dist::DiscreteUnivariateDistribution,
    initial_dist::DiscreteUnivariateDistribution,
    P::Matrix{Float64}
    )
    return PathCooccurrence(collect(1:size(P)[1]), length_dist, initial_dist, P)
end

function Base.rand(d::PathCooccurrence{String})
    m = rand(d.length_dist) # Sample length
    log_prob = logpdf(d.length_dist, m) # Intiialise log probability
    unique_vals = Set{Int}() # To store unique vertices observed in path
    z = Vector{String}(undef, m) # Initialisation of storage

    tmp = rand(d.initial_dist) # Sample initial vertex index
    z[1] = d.vertex_set[tmp] # Index the corresponding vertex
    log_prob += logpdf(d.initial_dist, tmp) # Increment log probability
    push!(unique_vals, tmp) # Store value sampled

    for i in 2:m
        component = rand(unique_vals) # Sample component
        tmp = rand(d.P_dists[component]) # Sample vertex index
        z[i] = d.vertex_set[tmp] # Store vertex value
        # Increment log probability (sum over unique values and normalise by number of unique values)
        for val in unique_vals
            log_prob += d.P[val,tmp]/length(unique_vals)
        end
        push!(unique_vals, tmp) # Update unqiue values seen
    end
    return z, log_prob
end

function Base.rand(d::PathCooccurrence{Int})
    m = rand(d.length_dist) # Sample length
    log_prob = logpdf(d.length_dist, m)  # Inistialise log probability

    unique_vals = Set{Int}() # Storage unique vertices observed in path
    z = Path(Vector{Int}(undef, m)) # Storage for sampled observation

    # Draw initial value
    tmp = rand(d.initial_dist) # Sample value
    z[1] = tmp # Store it
    log_prob += logpdf(d.initial_dist, tmp) # Increment log probability
    push!(unique_vals, tmp) # Store value sampled

    for i in 2:m
        component = rand(unique_vals)  # Sample component (which row/col of coocurrence matrix)
        tmp = rand(d.P_dists[component]) # Sample value from component distribution
        z[i] = tmp # Store value
        # Increment log probability (sum over unique values)
        for val in unique_vals
            log_prob += d.P[val,tmp]/length(unique_vals)
        end
        push!(unique_vals, tmp) # Store value sampled
    end
    return z, log_prob
end

function Distributions.logpdf(d::PathPseudoUniform, p::Path{Int})
    m = length(p)
    return logpdf(d.length_dist, m) - m*log(length(d.vertex_set))
end


text:
"""
PSFModels.PSFModel{T} <: AbstractMatrix{T}

Abstract type for PSF models.

In general, all `PSFModel`s have a set of pre-determined axes (the size is set upon creation) but they are lazy. That is, no memory is allocated and the values are calculated on the fly.

# Interface
The interface to define a model is as follows (for an example model `Model`)

| method | description |
|:-------|:------------|
| `Model()` | constructor(s) |
| `Base.size(m::Model)` | size, necessary for `AbstractArray` interface |
| `Base.axes(m::Model)` | axes, necessary for `AbstractArray` interface |
| `(m::Model)(point::AbstractVector)` | evaluate the model at the point in 2d space (x, y) |

browsing through the implementation of [`PSFModels.Gaussian`](@ref) should give a good idea of how to create a model
"""
abstract type PSFModel{T} <: AbstractMatrix{T} end

# in general, parse to static vector
(model::PSFModel)(point...) = model(SVector(point))
(model::PSFModel)(point::Tuple) = model(SVector(point))
# disallow index to avoid confusion
(model::PSFModel)(::CartesianIndex) = error("PSF models should be indexed using `getindex`, (equivalently `[]`)")

# getindex just calls model with reversed indices
Base.getindex(model::PSFModel, idx::Vararg{<:Integer,2}) = model(reverse(idx))
# always inbounds
Base.checkbounds(::Type{Bool}, ::PSFModel, inds...) = true
Base.checkbounds(::Type{Bool}, ::PSFModel, inds::CartesianIndex) = true

# broadcasting hack to slurp other axes (doesn't work for numbers)
Broadcast.combine_axes(psf::PSFModel, other::AbstractMatrix) = axes(other)
Broadcast.combine_axes(other::AbstractMatrix, psf::PSFModel) = axes(other)

## get the position depending on the keyword inputs
function _position(nt::NamedTuple{(:r, :theta)})
    # note: theta is in degrees
    th_rad = deg2rad(nt.theta)
    return CartesianFromPolar()(Polar(nt.r, th_rad))
end
_position(nt::NamedTuple{(:pos,)}) = _position(nt.pos)
_position(pos::AbstractVector) = SVector{length(pos)}(pos)
_position(pos::SVector) = pos
_position(pos::Tuple) = SVector(pos)
_position(nt::NamedTuple{(:x, :y)}) = SA[nt.x, nt.y]
_position(nt::NamedTuple{()}) = SA[0, 0]

@kwcall _position(r, theta)
@kwcall _position(pos)
@kwcall _position(x=0, y=0)
# theta is in degrees
@kwalias _position [
    ρ => r,
    θ => theta
]

function indices_from_extent(pos, extent)
    half_extent = @. 0.5 * extent
    lower = @. floor(Int, pos - half_extent)
    upper = @. ceil(Int, pos + half_extent)
    return last(lower):last(upper), first(lower):first(upper)
end


text:
## Example for documentation

using LearnConvection
using Plots

# problem
problem  = Slack("TKE"; parameters=TKEMassFlux.TKEParameters())

# data
train    = ["general_strat_32_profiles.jld2"]
test     = ["general_strat_32_profiles.jld2"]
D        = 32
N        = 4

# kernel
k        = 1
logγ     = 0.4
logσ     = 0.0
distance = euclidean_distance
kernel   = get_kernel(k, logγ, logσ, distance)

# data
𝒟_train  = LearnConvection.Data.data(train, problem; D=D, N=N);
𝒟_test   = LearnConvection.Data.data(test, problem; D=D, N=N);

# 𝒢 is trained on 𝒟_train
𝒢 = LearnConvection.GaussianProcess.model(𝒟_train; kernel = kernel, stencil_size=8)

# animate the mean GP prediction, where 𝒢 is tested on 𝒟_test
anim = animate_profile_and_model_output(𝒢, 𝒟_test)
gif(anim, "basic_example_8.gif", fps=15)

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
const FloatBits = NamedTuple{(:all, :sign, :exponent, :significand)}

"""
     StandardFloat

A `StandardFloat` is a floating point type that respects IEEE754-2019.

This registers the type so it is available for use with pkg functions.
"""
struct StandardFloat{F} # item Type, is Registered
    type::Union{Nothing, T} where T
    bits::FloatBits
end

function StandardFloat(::Type{T}, bits::FloatBits) where T
   bitsof(T) == sum(bits)-bits.all || throw(DomainError(T, "sum(bits) != bitsof(T)")
   @eval begin
    if isdefined(Base, :exponent_bits)
        Base.exponent_bits(::Type{$T}) = $bits.exponent
        Base.significand_bits(::Type{$T}) = $bits.significand
    else
        exponent_bits(::Type{$T}) = $bits.exponent
        significand_bits(::Type{$T}) = $bits.significand
    end
    sign_bits(::Type{$T}) = $bits.sign
   end
   return StandardFloat{Nothing}(T, bits)
end

function register_ieeefloat(stdfloat::StandardFloat{F}) where {F}
end

"""
    Float128

This placeholder is used to reference related constants.
""" Float128

"""
    Float8

This placeholder is used to reference related constants.
""" Float8

"""
    Binary8

This placeholder is used to reference related constants.
""" Binary8

if !isdefined(Main, :Float128)
    primitive type Float128 <: AbstractFloat 128 end
end
if !isdefined(Main, :Float128)
    primitive type Float8   <: AbstractFloat 8 end
end
if !isdefined(Main, :Binary8)
    primitive type Binary8  <: AbstractFloat 8 end
end

sign_bits(::Type{Float128}) = 1
sign_filter(::Type{Float128}) = ~zero(UInt128) >> 1
sign_mask(::Type{Float128}) = ~sign_filter(Float128)

sign_bits(::Type{Float8}) = 1
sign_filter(::Type{Float8}) = ~zero(UInt8) >> 1
sign_mask(::Type{Float8}) = ~sign_filter(Float8)

if isdefined(Base, :exponent_bits)
    Base.exponent_bits(::Type{Float128}) = 15
    Base.exponent_bits(::Type{Float8}) = 3
else
    exponent_bits(::Type{Float128}) = 15
    exponent_bits(::Type{Float8}) = 3
end

if isdefined(Base, :exponent_mask)
    Base.exponent_mask(::Type{Float128}) = ((UInt128(1) << Base.significand_bits(Float128)) - UInt128(1)) | sign_mask(Float128)
    Base.exponent_mask(::Type{Float8}) = ((UInt8(1) << Base.significand_bits(Float8)) - UInt8(1)) | sign_mask(Float8)
else
    exponent_mask(::Type{Float128}) = ((UInt128(1) << Base.significand_bits(Float128)) - UInt128(1)) | sign_mask(Float128)
    exponent_mask(::Type{Float8}) = ((UInt8(1) << Base.significand_bits(Float8)) - UInt8(1)) | sign_mask(Float8)
end

if isdefined(Base, :significand_bits)
    Base.significand_bits(::Type{Float128}) = 112
    Base.significand_bits(::Type{Float8}) = 4
else
    significand_bits(::Type{Float128}) = 112
    significand_bits(::Type{Float8}) = 4
end

if isdefined(Base, :significand_mask)
    Base.significand_mask(::Type{Float128}) = (UInt128(1) << Base.significand_bits(Float128)) - UInt128(1)
    Base.significand_mask(::Type{Float8}) = (UInt8(1) << Base.significand_bits(Float8)) - UInt8(1)
else
    significand_mask(::Type{Float128}) = (UInt128(1) << Base.significand_bits(Float128)) - UInt128(1)
    significand_mask(::Type{Float8}) = (UInt8(1) << Base.significand_bits(Float8)) - UInt8(1)
end

exponent_filter(::Type{Float128}) = ~exponent_mask(Float128)
exponent_filter(::Type{Float8}) = ~exponent_mask(Float8)

significand_filter(::Type{Float128}) = ~significand_mask(Float128)
significand_filter(::Type{Float8}) = ~significand_mask(Float8)


text:
title = "LyceumAI"
weight = 15


text:
### INFO
# Try to test parameter update.
# Taking a subsample of data
# Increasing the observation noise to Σ= 10^(-2)

using Distributions # to define priors
using Random        # to seed the random number generator
using DataFrames
using CSV
using BridgeSDEInference

HOME_DIR = joinpath(Base.source_dir(), "..", "..", "..", "..")
AUX_DIR = joinpath(pwd(), "src", "auxiliary")
OUT_DIR = joinpath(HOME_DIR, "output")
include(joinpath(AUX_DIR, "read_and_write_data.jl"))
include(joinpath(AUX_DIR, "read_JR_data.jl"))
include(joinpath(AUX_DIR, "transforms.jl"))

# Using BridgeSDEInference

# Fetch the data
fptObsFlag = false
filename = "jr_path_part_obs_3n.csv"
(df, x0, obs, obs_time, fpt,
      fptOrPartObs) = readDataJRmodel(Val(fptObsFlag), joinpath(OUT_DIR, filename))
#take sub set of data: 1 second
obs = obs[1:128]
obs_time = obs_time[1:128]

# Initial parameter guess
θ_true = [3.25, 100.0, 22.0, 50.0 , 135.0, 5.0, 6.0, 0.56, 0.0, 220.0, 0.0, 0.01, 2000.0, 1.0]
θ₀ = [3.25, 100.0, 22.0, 70.0 , 175.0, 5.0, 6.0, 0.56, 0.0, 150.0, 0.0, 0.01, 1500.0, 1.0]
# Boolean for parameters to update
#param_bool = [0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0]
#θ₀ = θ_true.*(1 .+ 0.5.*param_bool.*randn(length(θ_true)))

# P_Target
Pˣ = JRNeuralDiffusion3n(θ₀...)

#

#P_auxiliary
sigmas = θ₀[12:14]
P̃ = [JRNeuralDiffusion3nAux(sigmas..., t₀, u[1], T, v[1]) for (t₀,T,u,v)
     in zip(obs_time[1:end - 1], obs_time[2:end], obs[1:end - 1], obs[2:end])]

model_setup = DiffusionSetup(Pˣ,P̃, PartObs())

# Observation operator and noise
L = @SMatrix [0. 1. -1. 0. 0. 0.]
Σdiagel = 10^(-2)
Σ = @SMatrix [Σdiagel]
set_observations!(model_setup, [L for _ in P̃], [Σ for _ in P̃], obs, obs_time, fpt)

# Imputation grid < observation frequency
dt = 0.0005
set_imputation_grid!(model_setup, dt)

# Set prior on initial observation
x0 = ℝ{6}(0.08, 18, 15, -0.5, 0, 0)
set_x0_prior!(model_setup, GsnStartingPt(x0.*(1 .+ 0.001.*randn(6)), 0.001.*one(SMatrix{6,6,Float64})), x0.*(1 .+ 0.001.*randn(6)))
#set_x0_prior!(model_setup, KnownStartingPt(x0))

initialise!(eltype(x0), model_setup, Vern7(), false, NoChangePt(100))
set_auxiliary!(model_setup; skip_for_save=10^0, adaptive_prop=NoAdaptation())

pCN = 0.99
imputation = Imputation(NoBlocking(), pCN, Vern7())
readj = (100, 0.0005, 0.01, 999.9, 0.234, 50) # readjust_param
readj2 = (100, 0.0001, 0.01, 999.9, 0.234, 50)
p1 = ParamUpdate(MetropolisHastingsUpdt(), 4, θ₀,
            UniformRandomWalk(0.1, true), ImproperPosPrior(),
            UpdtAuxiliary(Vern7(), false) , readj
            )
p2 = ParamUpdate(MetropolisHastingsUpdt(), 10, θ₀,
            UniformRandomWalk(0.2, true), ImproperPosPrior(),
            UpdtAuxiliary(Vern7(), false) , readj
            )
p3 = ParamUpdate(MetropolisHastingsUpdt(), 5, θ₀,
            UniformRandomWalk(0.05, true), ImproperPosPrior(),
            UpdtAuxiliary(Vern7(), false), readj2
            )
p4 = ParamUpdate(MetropolisHastingsUpdt(), 13, θ₀,
            UniformRandomWalk(0.05, true), ImproperPosPrior(),
            UpdtAuxiliary(Vern7(), true), readj2
            )

mcmc_setup = MCMCSetup(imputation, p1, p2, p3, p4)
#mcmc_setup = MCMCSetup(imputation, p1, p2, p3, p4)

#mcmc_setup = MCMCSetup(imputation, p4)
#readjust=(x->x%100==0)
schedule = MCMCSchedule(5*10^4, [[1], [2,3,4], [5]],
                  (save=10^2, verbose=5*10^3, warm_up=100,
                   readjust=(x->x%100==0), fuse=(x->false)))

Random.seed!(4)
out, chains = mcmc(mcmc_setup, schedule, model_setup)
# Print out statement of the acceptance rate
println("")
println("acceptance rates for the for parameters are")
println("b: ", mean(chains.updates[2].accpt_history), ", μy: ",mean(chains.updates[3].accpt_history),
        " C: ", mean(chains.updates[4].accpt_history), " σy: ", mean(chains.updates[5].accpt_history))

# length(0:0.0001:obs_time[100])
error("STOP HERE")

using Makie

p1 = scatter(obs_time, [obs[i][1] for i in 1:length(obs)],  markersize = 0.05)
lines!(p1, out.time, [out.paths[end][i][2] - out.paths[end][i][3] for i in 1:length(out.paths[1])], color = (:red, 0.1))

p2 = Scene()
lines!(p2, out.time, [out.paths[end][i][1]  for i in 1:length(out.paths[1])], color = (:red, 0.5))

p3 = Scene()
lines!(p3, out.time, [out.paths[end][i][4]  for i in 1:length(out.paths[1])], color = (:red, 0.5))
p4 = Scene()
lines!(p4, out.time, [out.paths[end][i][5]  for i in 1:length(out.paths[1])], color = (:red, 0.5))

p5 = Scene()
lines!(p5, out.time, [out.paths[end][i][6]  for i in 1:length(out.paths[1])], color = (:red, 0.5))
out.paths

pscene1 = hbox(
    vbox(p3, p5),
    p4,
    vbox(p1, p2)
)
Makie.save(joinpath(Base.source_dir(), "output","reconstructed_data_with_parameter_estimation1234.png"), pscene1)

using Makie
chain = chains.θ_chain
param_b = Makie.lines([chain[i][4] for i in 1:length(chain)])
Makie.lines!(param_b, [0.0, length(chain)],[50, 50], color = (:red, 0.5))
param_b = Makie.title(param_b, "param b")
param_μy = Makie.lines([chain[i][10] for i in 1:length(chain)])
Makie.lines!(param_μy, [0.0, length(chain)],[220, 220], color = (:red, 0.5))
param_μy = Makie.title(param_μy, "param μy")
param_C = Makie.lines([chain[i][5] for i in 1:length(chain)])
Makie.lines!(param_C, [0.0, length(chain)],[135, 135], color = (:red, 0.5))
param_C =  Makie.title(param_C, "param C")
param_σy = Makie.lines([chain[i][13] for i in 1:length(chain)])
Makie.lines!(param_σy, [0.0, length(chain)],[2000, 2000], color = (:red, 0.5))
param_σy =  Makie.title(param_σy, "param σy")

pscene2 = hbox(
    vbox(param_μy, param_σy),
    vbox(param_b, param_C)
    )

Makie.save(joinpath(Base.source_dir(), "output", "parameter_estimation1234.png"), pscene2)


text:
const A_0 = 4.5#6          # Growth rate adjustment parameter
const α = 4.15e-5 * 24 * 10^6 / (24 * 60 * 60)#3.75e-5 * 24 * 10^6 / (24 * 60 * 60) # photosynthetic efficiency
const C_min = 0.01     # Minimal carbon reserve
const C_struct = 0.2   # Amount of carbon per unit dry weight of structural mass
const γ = 0.5# Exudation parameter
const ϵ = 0.22   # Frond erosion parameter
const I_sat = 90 * 24*60*60/(10^6)#200 * 24 * 60 * 60 / (10^6)  # Irradiance for maximal photosynthesis
const J_max = 1.4e-4 * 24# Maximal nitrate uptake (gN/dm^2/h converted to gN/dm^2/day)
const K_A = .5#0.6        # Structural dry weight per unit area
const K_DW = 0.0785    # Dry weight to wet weight ratio of structural mass
const K_C = 2.1213     # Mass of carbon reserves per gram carbon
const K_N = 2.72       # Mass of nitrogen reserves per gram nitrogen
const N_min = 0.0126#0.01     # Minimal nitrogen reserve
const N_max = 0.0216#0.022    # Maximal nitrogen reserve
const m_2 = 0.039/(2*(1-N_min/N_max))#0.036# 0.03       # Growth rate adjustment parameter
const m_1 = 0.18/(2*(1-N_min/N_max))-m_2# 0.1085     # Growth rate adjustment parameter
const μ_max = 0.18    # Maximal area specific growth ratio
const N_struct = 0.0146#0.01  # Amount of nitrogen per unit dry weight of structural mass
const P_1 = 1.22e-3 * 24    # Maximal photosynthetic rate at T = T?P1K converted to day^-1
const P_2 = 1.3e-3 * 24    # Maximal photosynthetic rate at T = T?P2K converted to day^-1
const a_1 = 0.85       # Photoperiod parameter
const a_2 = 0.3        # Photoperiod parameter
const R_1 =  2.785e-4 * 24   # Respiration rate at T = TR1,
const R_2 = 5.429e-4 * 24   # Respiration rate at T = TR2,
const T_R1 = 285       # Reference temperature for respiration (K)
const T_R2 = 290       # Reference temperature for respiration (K)
const T_P1 = 285       # Reference temperature for photosynthesis (K)
const T_P2 = 288       # Reference temperature for photosynthesis (K)
const T_AP = (1 / T_P1 - 1 / T_P2)^(-1) * log(P_2 / P_1)# 1694.4    # Arrhenius temperature for photosynthesis (K)
const T_PL = 271
const T_PH = 296
const T_APH = 1414.87# 1516.7#25924    # Arrhenius temperature for photosynthesis at high end of range (K)
const T_APL = 4547.89# 4388.5#27774    # Arrhenius temperature for photosynthesis at low end of range (K)
const T_AR = (1 / T_R1 - 1 / T_R2)^(-1) * log(R_2 / R_1)     # Arrhenius temperature for respiration (K)
const U_0p65 = 0.03    # Current speed at which J = 0.65Jmax (m/s)
const K_X = 4          # Nitrate uptake half saturation constant, converted to mmol/m^3

const R_A=1.11e-4*24
const R_B=5.57e-5*24


text:
# This fille contains functions for loading atmosphere model
# curves from the Sonora Bobcat models

using DelimitedFiles
using NamedTupleTools
using BasicInterpolators
using Interpolations

"""
    sonora_photometry_interpolator(:Keck_L′, [metalicity="+0.0"])

Given a supported photometric band and [M/H] metalicity (default=solar),
return a function of temperature (K) and mass (M_jup) that gives the
absolute magnitude of the planet in that bandpass.

Supported bands:
:MKO_Y, :MKO_Z, :MKO_J, :MKO_H, :MKO_K, :MKO_L′, :MKO_M′, :TwoMASS_J, :TwoMASS_H, :TwoMASS_Ks, :Keck_Ks, :Keck_L′, :Keck_Ms, :SDSS_g′, :SDSS_r′, :SDSS_i′, :SDSS_z′, :IRAC_36, :IRAC_45, :IRAC_57, :IRAC_79, :WISE_W1, :WISE_W2, :WISE_W3, :WISE_W4

Supported metalicities:
"+0.0", "-0.5", "+0.5"
"""
function sonora_photometry_interpolator(band, metalicity="+0.0";catalog=datadep"SonoraBobcatEvoPhot")

    #  Load Sonora magnitude table
    mag_table = load_table(joinpath(catalog, "photometry_tables", "mag_table"*metalicity))

    # We first use BasicInterpolators to grid the sparsely sampled models
    # Then we use Interpolations to make a fast, autodiff compatible
    # linear interpolator of the resulting data.
    points = [
        mag_table.Teff./10 mag_table.mass
    ]
    samples = collect(mag_table[band])

    sitp = RBFInterpolator(points, samples, 2)
    minmass, maxmass = extrema(mag_table.mass)
    minT, maxT = extrema(mag_table.Teff)
    temp_mass_to_abs_mag = (teffk, mass_mjup) -> sitp(teffk/10, mass_mjup)
    teff = range(minT, maxT,length=200)
    mass = range(minmass, maxmass,length=200)
    gridded = temp_mass_to_abs_mag.(teff, mass')

    # Now build linear interpolator
    litp = LinearInterpolation((teff./10, mass), gridded, extrapolation_bc=NaN)
    function model_interpolator(teffk, mass)
        if minT <= teffk <= maxT && minmass < mass < maxmass
            return litp(teffk/10, mass)
        else
            return NaN
        end
    end

    return model_interpolator
end
export sonora_photometry_interpolator

function sonora_cooling_interpolator(metalicity="+0.0";catalog=datadep"SonoraBobcatEvoPhot")

    # Load Sonora cooling track
    # valid_lines = filter(readlines("./sonora/evolution_tables/evo_tables+0.0/nc+0.0_co1.0_age")) do line
    valid_lines = filter(readlines(joinpath(catalog, "evolution_tables", "evo_tables"*metalicity, "nc$(metalicity)_co1.0_age"))) do line
        length(line) > 10
    end
    headers = split(valid_lines[1], r"  +")[2:7]
    headers = lowercase.(replace.(headers, r"\W"=>""))
    grid= mapreduce(vcat, valid_lines[2:end]) do line
        permutedims(parse.(Float64, split(line, r"  +")[2:7]))
    end
    # cooling_grid = DataFrame(grid, headers)
    # return grid, headers

    iteffk = findfirst(==("teffk"), headers)
    immsun = findfirst(==("mmsun"), headers)
    iagegyr = findfirst(==("agegyr"), headers)
    agegyr = grid[:,iagegyr]
    agemyr = agegyr .* 1e3
    mmsun = grid[:,immsun]
    mmjup = mmsun ./ mjup2msol
    teffk = grid[:,iteffk]

    points = [
        log.(agemyr) mmjup
    ]
    samples = collect(teffk)

    sitplog = RBFInterpolator(points, samples, 0.5)

    sitp = (agemyr, mmjup) -> sitplog(log(agemyr), mmjup)
    # return (;sitp, agemyr, mmjup, teffk)

    # Now build a fast linear interpolator over this grid.
    minmmjup, maxmmjup = extrema(mmjup)
    minagemyr, maxagemyr = extrema(agemyr)
    agemyrrange = range(minagemyr, maxagemyr,length=2000)
    mmjuprange = range(minmmjup, maxmmjup,length=500)
    gridded = sitp.(agemyrrange, mmjuprange')

    # Now build linear interpolator
    litp = LinearInterpolation((agemyrrange, mmjuprange), gridded, extrapolation_bc=NaN)
    function model_interpolator(agemyr, mmjup)
        if minagemyr <= agemyr <= maxagemyr && minmmjup < mmjup < maxmmjup
            return litp(agemyr, mmjup)
        else
            return NaN
        end
    end

    return model_interpolator
    # return sitp

    # Want (age, mass) -> teffk. This can be combined with the others to get specific fluxes.

end
export sonora_cooling_interpolator

function load_table(fname)

    headers = open(fname, lock = false, read = true) do f
        readline(f)
        readline(f)
        readline(f)
        readline(f)
        readline(f)
        readline(f)
        readline(f)
        readline(f)
        h1 = readline(f)
        h2 = readline(f)

        headers_1 = strip.(split(h1, '|'))
        headers_2 = strip.(split(h2, r"  +"))

        # Starting after 6
        headers_2_1_indices =
            [2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7]

        combined = vcat(
            headers_2[1:6],
            headers_1[headers_2_1_indices] .* '_' .* headers_2[7:end],
        )
        # Normalize headers
        combined = replace.(combined, '\'' => '′')
        combined = replace.(combined, ' ' => "")
        combined = replace.(combined, "/" => '_')
        combined = replace.(combined, "2MASS" => "TwoMASS")
        combined = replace.(combined, r"[^\w′_]" => "")
        return combined
    end

    data = readdlm(fname, String, skipstart = 10, header = false)
    data = [
        try
            typeof(d) <: AbstractString ? parse(Float64, replace(d, '*' => "")) : d
        catch err
            NaN
        end
        for d in data
    ]
    # Return simple table
    return namedtuple(headers, eachcol(data))
end


text:
using PyPlot, Statistics, Printf, Random, Clustering

""" Preconditioner based on maximum entropy method. """
function maxentropy(r::Int,data::Array{Float64,1})

    N = length(data)
    n = N ÷ r

    # throw away random data for equally sized chunks of data
    data = shuffle(data[:])[1:n*r]
    sort!(data)

    nodes = Array{Float64,1}(undef,r)
    nodes[1] = minimum(data)
    nodes[2] = (2*data[1] + data[n] + data[n+1])/4

    for i in 2:r-1
        nodes[i+1] = (data[(i-1)*n] + data[(i-1)*n + 1] + data[i*n] + data[i*n + 1])/4
    end

    nodes[end-1] = (2*data[r*n] + data[(r-1)*n] + data[(r-1)*n-1])/4
    nodes[end] = maximum(data)

    return nodes
end

# """ Sign preserving 1D distance between nodes and a data sample x. """
# function distance!(dist::Array{Float64,1},nodes::Array{Float64,1},x::Float64)
#     n = length(nodes)
#     @boundscheck n == length(dist) || throw(BoundsError())
#     @inbounds for i in 1:n
#         dist[i] = x - nodes[i]
#     end
# end

""" Sign preserving 1D distance between nodes and a data sample x. """
function distance(nodes::Array{Float64,1},x::Float64)
    return x .- nodes
end

""" Find the Best Matching Unit (BMU) given an array dist of distances."""
function bestMatchingUnit(dist::Array{Float64,1})
    return argmin(abs.(dist))
end

function neighbourhood(bmu::Int,s::Int,smax::Int,nodes::Array{Float64,1},w::Float64)

    decay = αexp(s,1e-3)
    shrink = αexp(s,5e-2)
    # decay = αlin(s,smax)
    # shrink = αlin(s,smax)

    # Gaussian neighbourhood function
    #d = decay*exp.(-(((nodes .- nodes[bmu])/(w*shrink)).^2)./2)

    # step function
    d = decay*crenel.(nodes .- nodes[bmu],w*shrink)

    return d
end

function αexp(s::Int,r::Float64=0.05)
    return exp(-r*s)
end

function αlin(s::Int,smax::Int)
    return 1-(s-1)/smax
end

function update!(s::Int,smax::Int,data::Array{Float64,1},nodes::Array{Float64,1},w::Float64=1.0)
    x = rand(data)
    dist = distance(nodes,x)
    bmu = bestMatchingUnit(dist)
    θ = neighbourhood(bmu,s,smax,nodes,w)
    for i in 1:length(nodes)
        nodes[i] += θ[i]*sign(dist[i])
    end
end

function average_rounding_error(data::Array{Float64,1},nodes::Array{Float64,1})
    n = length(data)
    S = 0.0
    for i = 1:n
        S += minimum(abs.(distance(nodes,data[i])))/n
    end
    return S
end

data = vcat(randn(500),randn(500).+4.0)

Nnodes = 32
nodesmaxent = maxentropy(Nnodes,data)
nodes = collect(range(0,stop=1,length=Nnodes))

w = maxspread(data)/length(nodes)

# smax = 600
# er = average_rounding_error(data,nodes)
# ser = @sprintf("%.8f",er)
# println("\nStart: $ser")
# global ermin = er
#
# for s in 1:smax
#     update!(s,smax,data,nodes,w)
#     er = average_rounding_error(data,nodes)
#
#     global ermin
#     if er < ermin
#         ermin = er
#     end
#
#     ser = @sprintf("%.6f",er)
#     print("\r\u1b[K")
#     progress = Int(round(s/smax*100))
#     print("$progress%, $ser")
# end
# er = average_rounding_error(data,nodes)
# ser = @sprintf("%.8f",er)
# println(", End: $ser")

##
R = kmeans(vcat(data',data'),Nnodes)
M = R.centers

nodes = sort(M[1,:])

##
s=0
fig,ax = subplots(1,1,figsize=(8,2))
ax.plot(data,ones(length(data)),"|",lw=0.1)
l1, = ax.plot(nodes,zeros(length(nodes)),".-")
l2, = ax.plot(nodesmaxent,zeros(length(nodes)) .- 0.5,".-")
yticks([])
ylim(-1,2)
title("s = 0", loc="left")
tight_layout()

er = average_rounding_error(data,nodesmaxent)
ser = @sprintf("%.8f",er)
println("\nStart: $ser")
title("s = $s, r = $ser", loc="left")

##
smax = 600
for s in 1:smax
    update!(s,smax,data,nodes,w)
    er = average_rounding_error(data,nodes)
    ser = @sprintf("%.4f",er)
    title("s = $s, r = $ser", loc="left")
    l1.set_data(nodes,zeros(length(nodes)))
    pause(αexp(s,0.5)+0.05)
    #savefig(@sprintf "frames/frame%04d.jpg" s)
    print("\r\u1b[K")
    progress = Int(round(s/smax*100))
    print("$progress%, $ser")
end


text:
# This file is not used; it only exists as a placeholder so that `process_file` can
# be called on it in order to read `CoverageCore.jl.cov`.


text:
abstract type AbstractCFRPolicy <: AbstractPolicy end

function Base.run(p::AbstractCFRPolicy, env::AbstractEnv, stop_condition=StopAfterStep(1), hook=EmptyHook())
    @assert NumAgentStyle(env) isa MultiAgent
    @assert DynamicStyle(env) === SEQUENTIAL
    @assert RewardStyle(env) === TERMINAL_REWARD
    @assert ChanceStyle(env) === EXPLICIT_STOCHASTIC
    @assert DefaultStateStyle(env) isa Information

    RLBase.reset!(env)

    while true
        update!(p, env)
        hook(POST_ACT_STAGE, p, env)
        stop_condition(p, env) && break
    end
    update!(p)
end


text:
"""
`KernelSum(kernels::Array{Kernel};weights::Array{Real}=ones(length(kernels)))`
Create a positive weighted sum of kernels.
One can also use the operator `+`
```
k1 = SqExponentialKernel()
k2 = LinearKernel()
k = KernelSum([k1,k2])
kernelmatrix(k,X) == kernelmatrix(k1,X).+kernelmatrix(k2,X)
kernelmatrix(k,X) == kernelmatrix(k1+k2,X)
kweighted = 0.5*k1 + 2.0*k2
```
"""
struct KernelSum <: Kernel
    kernels::Vector{Kernel}
    weights::Vector{Real}
end

function KernelSum(
    kernels::AbstractVector{<:Kernel};
    weights::AbstractVector{<:Real} = ones(Float64, length(kernels)),
)
    @assert length(kernels) == length(weights) "Weights and kernel vector should be of the same length"
    @assert all(weights .>= 0) "All weights should be positive"
    return KernelSum(kernels, weights)
end

Base.:+(k1::Kernel, k2::Kernel) = KernelSum([k1, k2], weights = [1.0, 1.0])
Base.:+(k1::ScaledKernel, k2::ScaledKernel) = KernelSum([kernel(k1), kernel(k2)], weights = [first(k1.σ²), first(k2.σ²)])
Base.:+(k1::KernelSum, k2::KernelSum) =
    KernelSum(vcat(k1.kernels, k2.kernels), weights = vcat(k1.weights, k2.weights))
Base.:+(k::Kernel, ks::KernelSum) =
    KernelSum(vcat(k, ks.kernels), weights = vcat(1.0, ks.weights))
Base.:+(k::ScaledKernel, ks::KernelSum) =
        KernelSum(vcat(kernel(k), ks.kernels), weights = vcat(first(k.σ²), ks.weights))
Base.:+(k::ScaledKernel, ks::Kernel) =
        KernelSum(vcat(kernel(k), ks), weights = vcat(first(k.σ²), 1.0))
Base.:+(ks::KernelSum, k::Kernel) =
    KernelSum(vcat(ks.kernels, k), weights = vcat(ks.weights, 1.0))
Base.:+(ks::KernelSum, k::ScaledKernel) =
        KernelSum(vcat(ks.kernels, kernel(k)), weights = vcat(ks.weights, first(k.σ²)))
Base.:+(ks::Kernel, k::ScaledKernel) =
        KernelSum(vcat(ks, kernel(k)), weights = vcat(1.0, first(k.σ²)))
Base.:*(w::Real, k::KernelSum) = KernelSum(k.kernels, weights = w * k.weights) #TODO add tests

Base.length(k::KernelSum) = length(k.kernels)

kappa(κ::KernelSum, x, y) = sum(κ.weights[i] * kappa(κ.kernels[i], x, y) for i in 1:length(κ))

function kernelmatrix(κ::KernelSum, X::AbstractMatrix; obsdim::Int = defaultobs)
    sum(κ.weights[i] * kernelmatrix(κ.kernels[i], X, obsdim = obsdim) for i in 1:length(κ))
end

function kernelmatrix(
    κ::KernelSum,
    X::AbstractMatrix,
    Y::AbstractMatrix;
    obsdim::Int = defaultobs,
)
    sum(κ.weights[i] * _kernelmatrix(κ.kernels[i], X, Y, obsdim) for i in 1:length(κ))
end

function kerneldiagmatrix(
    κ::KernelSum,
    X::AbstractMatrix;
    obsdim::Int = defaultobs,
)
    sum(κ.weights[i] * kerneldiagmatrix(κ.kernels[i], X, obsdim = obsdim) for i in 1:length(κ))
end

function Base.show(io::IO,κ::KernelSum)
    printshifted(io,κ,0)
end

function printshifted(io::IO,κ::KernelSum, shift::Int)
    print(io,"Sum of $(length(κ)) kernels:")
    for i in 1:length(κ)
        print(io,"\n"*("\t"^(shift+1))*"- (w=$(κ.weights[i])) ")
        printshifted(io,κ.kernels[i],shift+2)
    end
end


text:
Pkg.add("DataFrames")
Pkg.add("Gadfly")
Pkg.add("Cairo")
Pkg.add("Fontconfig")
using DataFrames
using Gadfly
using Cairo
using Fontconfig

# Plot a Gaussian curve

style = Theme(
    major_label_font_size=35pt,
    minor_label_font_size=30pt,
    background_color="white",
    major_label_font="CMU Serif",
    minor_label_font="CMU Serif",
    line_width=5px
)

mus = [0.807033, 1.1918, 0.667843, 2.32305, 2.0586];
sigmas = [0.416768, 1.15322, 1.22268, 0.458837, 0.896344];
numbers = ["377", "11177", "127413", "140054", "7506"];

function getGaussian(sigma, mu, pos)
    expo = ((-1./2.0) .* (((pos-mu)./sigma).^2));
    return (1./(sigma.*sqrt(2.*3.14159265359))).*exp(expo);
end

for i = 1:5
    println("Creating Gaussian " , numbers[i] , "...");
    mu = mus[i];
    sigma = sigmas[i];

    x1=collect(0:.01:3);
    labels = ["                 Not Similar", "Slightly Similar", "Very Similar", "Nearly Identical                        "];

    output = plot(
        x=x1,
        y=getGaussian(sigma, mu, x1),
        Geom.line,
        Scale.x_continuous(labels = x -> get(labels, round(Int, x)+1, "")),
        Guide.xlabel("Answer"),
        Guide.ylabel("Gaussian"),
        Coord.cartesian(xmin=0),
        style
    );
    draw(PDF(string("graphs/" , numbers[i] , "-gaussian.pdf"), 1600px, 800px), output);
end


text:
@inline function _request_http(verb::AbstractString,
                               full_url::HTTP.URI,
                               headers::AbstractDict,
                               query::Nothing,
                               body::Nothing)
    response = HTTP.request(
        verb,
        full_url,
        headers,
    )
    return response
end

@inline function _request_http(verb::AbstractString,
                               full_url::HTTP.URI,
                               headers::AbstractDict,
                               query::Nothing,
                               body::AbstractString)
    response = HTTP.request(
        verb,
        full_url,
        headers,
        body,
    )
    return response
end

@inline function _request_http(verb::AbstractString,
                               full_url::HTTP.URI,
                               headers::AbstractDict,
                               query::AbstractDict,
                               body::Nothing)
    response = HTTP.request(
        verb,
        full_url,
        headers;
        query = query,
    )
    return response
end

@inline function _request_http(verb::AbstractString,
                               full_url::HTTP.URI,
                               headers::AbstractDict,
                               query::AbstractDict,
                               body::AbstractString)
    response = HTTP.request(
        verb,
        full_url,
        headers,
        body;
        query = query,
    )
    return response
end

@inline function _add_trailing_slash(url::HTTP.URI)::HTTP.URI
    _url_string = _get_http_uri_string(url)
    if endswith(_url_string, "/")
        return url
    end
    return HTTP.URI(string(_url_string, "/"))
end

@inline function _generate_full_url(client::Client, path::AbstractString)::HTTP.URI
    base_url = get_base_url(client)
    result = _generate_full_url(base_url, path)
    return result
end

@inline function _generate_full_url(base_url::BaseURL, path::AbstractString)::HTTP.URI
    base_url_uri = _get_http_uri(base_url)
    result = _generate_full_url(base_url_uri, path)
    return result
end

@inline function _generate_full_url(base_url_uri::HTTP.URI, path::AbstractString)::HTTP.URI
    base_url_uri_with_trailing_slash = _add_trailing_slash(base_url_uri)
    base_url_uri_string = _get_http_uri_string(base_url_uri_with_trailing_slash)
    full_url_uri_string = string(base_url_uri_string, path)
    full_url_uri = HTTP.URI(full_url_uri_string)
    return full_url_uri
end

"""
    request_raw(client, verb, path; query, body, headers)
    request_raw(client, verb, path; query, headers)
    request_raw(client, verb, path; body, headers)
    request_raw(client, verb, path; headers)
"""
@inline function request_raw(client::Client,
                             verb::AbstractString,
                             path::AbstractString;
                             body::Union{AbstractString, Nothing} = nothing,
                             headers::AbstractDict = Dict{String, String}(),
                             query::Union{AbstractDict, Nothing} = nothing)::String
    full_url = _generate_full_url(client,
                                  path)
    _new_headers = Dict{String, String}()
    xml_headers!(_new_headers)
    authentication_headers!(_new_headers, client)
    merge!(_new_headers, headers)
    response = _request_http(verb,
                             full_url,
                             _new_headers,
                             query,
                             body)
    empty!(_new_headers)
    response_body_string::String = String(response.body)::String
    return response_body_string
end

@inline function _write_xml_request_body(body::Nothing)::Nothing
    return nothing
end

@inline function _write_xml_request_body(body::AbstractDict)::String
    return XMLDict.dict_xml(body)
end

"""
    request_xml(client, verb, path; query, body, headers)
    request_xml(client, verb, path; query, headers)
    request_xml(client, verb, path; body, headers)
    request_xml(client, verb, path; headers)
"""
@inline function request_xml(client::Client,
                             verb::AbstractString,
                             path::AbstractString;
                             body::Union{AbstractDict, Nothing} = nothing,
                             headers::AbstractDict = Dict{String, String}(),
                             query::Union{AbstractDict, Nothing} = nothing)
    _new_request_body = _write_xml_request_body(body)
    response_body::String = request_raw(client,
                                        verb,
                                        path;
                                        body = _new_request_body,
                                        headers = headers,
                                        query = query)::String
    response_xml = XMLDict.xml_dict(response_body)
    return response_xml
end


text:
################################################################################
#
# Data filtering utilities
#
################################################################################

export cfreqfilt, cfreq_wavelet, applyfilt

### Filters

function cfreqfilt(nt::Array{Any, 1}, dt::Array{Any, 1}, cfreqs::NTuple{4, R}; type::String = "bandpass", edge = hamming, flag_fft::Bool = true, norm_fact::Union{Nothing, Tuple{judiVector, NTuple{2, R}}} = nothing)

	if flag_fft
		Filter = Array{Array{C, 1}, 1}(undef, length(dt))
	else
		Filter = Array{Array{R, 1}, 1}(undef, length(dt))
	end
	if norm_fact != nothing
		fourNormDat = computeFourierNorm(nt, norm_fact[1], norm_fact[2])
	end
	for i = 1:length(dt)
		Filter[i] = cfreqfilt(nt[i], dt[i], cfreqs; type = type, edge = edge, flag_fft = true)
		if norm_fact != nothing
			Filter[i] ./= fourNormDat[i]
		end
		if !flag_fft
			Filter[i] = real.(ifft(Filter[i]))
		end
	end
	return Filter

end

function cfreqfilt(nt::Int64, dt::R, cfreqs::NTuple{4, R}; type::String = "bandpass", edge = hamming, flag_fft::Bool = true)

	# Initialization
	Filter = zeros(C, nt)

	# Index corresponding to corner frequencies
	i1 = Int64(ceil(cfreqs[1]*nt*dt)+1)
	i2 = Int64(ceil(cfreqs[2]*nt*dt)+1)
	i3 = Int64(floor(cfreqs[3]*nt*dt)+1)
	i4 = Int64(floor(cfreqs[4]*nt*dt)+1)
	if type == "lowpass"
		i1 = 1
		i2 = 1
	elseif type == "highpass"
		i3 = nt
		i4 = nt
	end

	# Populating filter in the frequency domain
	Filter[1:i1-1] .= 0
	k = collect(i1:i2-1); Filter[k] .= (k.-i1)./(i2-i1)
	Filter[i2:i3] .= 1
	k = collect(i3+1:i4); Filter[k] .= R(1).-(k.-i3)./(i4-i3)
	Filter[i4+1:end] .= 0
	SL = Int64((nt-1)/2)
	Filter[SL+2:end] = Filter[SL+1:-1:2]

	# Zeroing edge values in time domain
	filter = ifftshift(R.(edge(length(Filter))).*fftshift(real.(ifft(Filter))))

	# Final Fourier transform
	if flag_fft
		return fft(filter)
	else
		return filter
	end

end

function cfreq_wavelet(SL::Int64, nt::Int64, dt::R, cfreqs::NTuple{4, R}; type::String = "bandpass", edge = hamming, flag_fft::Bool = true)

	return [cfreqfilt(2*SL+1, dt, cfreqs; type = type, edge = edge, flag_fft = false); zeros(R, nt-(2*SL+1))]

end

function whitening_filter(dat::Array{R, 2}, dt::R, B::NTuple{2, R})

	nt = size(dat, 1)
	Dat = fft(dat, 1)
	normDat = sqrt.(sum(abs.(Dat).^2, 2))
	df = 1/(nt*dt)
	freqs = ((1:nt).-1)*df
	idxB = findall(freqs .>= B[1] .& freqs .<= B[2])
	Filter = zeros(R, nt)
	Filter[idxB] = 1f0./normDat[idxB]

end

### Filter application

function applyfilt(dat::judiVector, Filter::Union{Array{Array{C, 1}, 1}, Array{C, 1}})

	dat_filt = Array{Array, 1}(undef, dat.nsrc)
	flag_arr = isa(Filter, Array{C, 1})
	for i = 1:dat.nsrc
		if flag_arr
			dat_filt[i] = applyfilt(dat.data[i], Filter)
		else
			dat_filt[i] = applyfilt(dat.data[i], Filter[i])
		end
	end
	return judiVector(dat.geometry, dat_filt)

end

function applyfilt(dat::Array{R, 2}, Filter::Array{C, 1})

	Dat = fft(zeropadding(dat, length(Filter)), 1)
	return real.(ifft(Dat.*Filter, 1)[1:size(dat, 1), :])

end

function zeropadding(dat::Array{R, 2}, pad::Int64)

	dat_ = zeros(R, (pad, size(dat, 2)))
	dat_[1:size(dat, 1), :] .= dat
	return dat_

end

###

function computeFourierNorm(nt::Array{Any, 1}, dat::judiVector, B::NTuple{2, R}; norm_src::Bool = true)

	fourNorm = Array{Array{R, 1}, 1}(undef, dat.nsrc)
	for i = 1:dat.nsrc
		fourNorm[i] = computeFourierNorm(nt[i], dat.geometry.dt[i], dat.data[i], B)
	end
	return fourNorm

end

function computeFourierNorm(nt::Int64, dt::R, dat::Array{R, 2}, B::NTuple{2, R})

	df = 1f0/(nt*dt)
	freqs = ifftshift(R.(-(nt-1)/2:(nt-1)/2))*df
	Dat = fft(zeropadding(dat, nt), 1)
	fourNorm = vec(sqrt.(sum(abs.(Dat).^2, dims = 2)))
	fourNorm[findall(.!((abs.(freqs) .>= B[1]) .& (abs.(freqs) .<= B[2])))] .= 1f0
	return fourNorm

end