text:
using Morbit
M = Morbit
MOI = M.MOI
#%%

Morbit.print_all_logs()

mop = M.MixedMOP(2)
f1 = x -> sum( (x .- 1).^2 )
f2 = x -> sum( (x .+ 1).^2 )

#M.add_objective!(mop, f1, M.ExactConfig())
#M.add_objective!(mop, f2, M.ExactConfig())
#M.add_objective!(mop, f1, M.RbfConfig(;use_max_points = false))
#M.add_objective!(mop, f2, M.RbfConfig(kernel = :gaussian))
#M.add_objective!(mop, f2, M.RbfConfig())
#M.add_objective!(mop, f1, M.TaylorConfig(; gradients = M.RFD.CFDStamp(1,3), hessians = M.RFD.CFDStamp(1,5)))
M.add_objective!(mop, f1, M.TaylorConfig())
M.add_objective!(mop, f2, M.TaylorConfig())

algo_config = AlgoConfig(;
	max_iter = 10
)
populated_db = nothing

x0 = rand(2)
fx0 = []
#%%
#mop, id, db, sc, ac = M.initialize_data(mop,x0);

X, FX, ret, db = Morbit.optimize( mop, x0; algo_config )


text:
# Error codes used in the device API.

# Convert a error code to a string for printing purposes.
# If the error code is unrecognized, errToStr returns "UNKNOWN".
# param errorCode a negative integer return code
# return a pointer to a string representing the error
function SoapySDR_errToStr(errorCode)
    ccall((:SoapySDR_errToStr, lib), Cstring, (Cint,), errorCode)
end

struct SoapySDRAPIError
    code::Cint
end
Base.show(io::IO, se::SoapySDRAPIError) = print(io, unsafe_string(SoapySDR_errToStr(se.code)))

# Returned when read has a timeout.
const SOAPY_SDR_TIMEOUT = -1

# Returned for non-specific stream errors.
const SOAPY_SDR_STREAM_ERROR = -2

# Returned when read has data corruption.
# For example, the driver saw a malformed packet.
const SOAPY_SDR_CORRUPTION = -3

# Returned when read has an overflow condition.
# For example, and internal buffer has filled.
const SOAPY_SDR_OVERFLOW = -4

# Returned when a requested operation or flag setting
# is not supported by the underlying implementation.
const SOAPY_SDR_NOT_SUPPORTED = -5

# Returned when a the device encountered a stream time
# which was expired (late) or too early to process.
const SOAPY_SDR_TIME_ERROR = -6

# Returned when write caused an underflow condition.
# For example, a continuous stream was interrupted.
const SOAPY_SDR_UNDERFLOW = -7

ext: jl


text:
export FileHeader

struct FileHeader
    th::String
    bfh::BinaryFileHeader
end

function show(io::IO, fh::FileHeader)
    show(fh.bfh)
end

FileHeader() = FileHeader(" "^3200, BinaryFileHeader())
FileHeader(fh::FileHeader) = FileHeader(fh.th, BinaryFileHeader(fh.bfh))
FileHeader(bfh::BinaryFileHeader) =
    FileHeader(" "^3200, BinaryFileHeader(fh.bfh))


text:
# # [Utilities (util.jl)](@id util)
# Unlike every other file in the pipeline, this file is not intended to be run directly.
# Instead, include this in other files. These utilities provide an interface to the Planck
# data products, namely
# 1. binning matrix
# 2. beam ``W_{\ell}^{XY} = B_{\ell}^X B_{\ell}^{Y}``
# 3. foreground model cross-spectra
# 4. ``\texttt{plic}`` reference covariance matrix and reference spectra, for comparison plots

using PowerSpectra
using DataFrames, CSV
using DelimitedFiles
using LinearAlgebra
using FITSIO

# ## Planck Binning
# Planck bins the spectra at the very end, and applies an ``\ell (\ell+1)`` relative
# weighting inside the bin. This utility function generates the binning operator
# ``P_{b\ell}`` such that ``C_b = P_{b \ell} C_{\ell}``. It also returns the mean of the
# left and right bin edges, which is what is used when plotting the Planck spectra.

"""
    util_planck_binning(binfile; lmax=6143)

Obtain the Planck binning scheme.

### Arguments:
- `binfile::String`: filename of the Planck binning, containing left/right bin edges

### Keywords
- `lmax::Int=6143`: maximum multipole for one dimension of the binning matrix

### Returns:
- `Tuple{Matrix{Float64}, Vector{Float64}`: returns (binning matrix, bin centers)
"""
function util_planck_binning(binfile; lmax=6143)
    bin_df = DataFrame(CSV.File(binfile;
        header=false, delim=" ", ignorerepeated=true))
    lb = (bin_df[:,1] .+ bin_df[:,2]) ./ 2
    P = binning_matrix(bin_df[:,1], bin_df[:,2], ℓ -> ℓ*(ℓ+1) / (2π); lmax=lmax)
    return P, lb[1:size(P,1)]
end

# ## Planck Beam
# The Planck effective beams are azimuthally-averaged window functions induced by the
# instrumental optics. This utility function reads the Planck beams from the RIMO, which
# are of the form `TT_2_TT`, `TT_2_EE` etc. Conventionally, the Planck spectra are stored
# with the diagonal of the beam-mixing matrix applied.
#
# ```math
# C_{\ell} = W^{-1}_{\ell} \hat{C}_{\ell}
# ```

"""
    util_planck_beam_Wl([T::Type=Float64], freq1, split1, freq2, split2, spec1_, spec2_;
        lmax=4000, beamdir=nothing)

Returns the Planck beam transfer of [spec1]_to_[spec2], in Wl form.

### Arguments:
- `T::Type=Float64`: optional first parameter specifying numerical type
- `freq1::String`: frequency of first field
- `split1::String`: split of first field (i.e. hm1)
- `freq2::String`: frequency of first field
- `split2::String`: split of second field (i.e. hm2)
- `spec1_`: converted to string, source spectrum like TT
- `spec2_`: converted to string, destination spectrum

### Keywords
- `lmax=4000`: maximum multipole
- `beamdir=nothing`: directory containing beam FITS files. if nothing, will fall back to
        the PowerSpectra.jl beam files.

### Returns:
- `SpectralVector`: the beam Wl, indexed 0:lmax
"""
function util_planck_beam_Wl(T::Type, freq1, split1, freq2, split2, spec1_, spec2_;
                        lmax=4000, beamdir=nothing)
    if isnothing(beamdir)
        @warn "beam directory not specified. switching to PowerSpectra.jl fallback"
        beamdir = PowerSpectra.planck256_beamdir()
    end
    spec1 = String(spec1_)
    spec2 = String(spec2_)

    if parse(Int, freq1) > parse(Int, freq2)
        freq1, freq2 = freq2, freq1
        split1, split2 = split2, split1
    end
    if (freq1 == freq2) && ((split1 == "hm2") && (split2 == "hm1"))
        split1, split2 = split2, split1
    end

    fname = "Wl_R3.01_plikmask_$(freq1)$(split1)x$(freq2)$(split2).fits"
    f = FITS(joinpath(beamdir, "BeamWf_HFI_R3.01", fname))
    Wl = convert(Vector{T}, read(f[spec1], "$(spec1)_2_$(spec2)")[:,1])
    if lmax < 4000
        Wl = Wl[1:lmax+1]
    else
        Wl = vcat(Wl, last(Wl) * ones(T, lmax - 4000))
    end
    return SpectralVector(Wl)
end
util_planck_beam_Wl(T::Type, freq1, split1, freq2, split2, spec1; kwargs...) =
    util_planck_beam_Wl(T, freq1, split1, freq2, split2, spec1, spec1; kwargs...)
util_planck_beam_Wl(freq1::String, split1, freq2, split2, spec1, spec2; kwargs...) =
    util_planck_beam_Wl(Float64, freq1, split1, freq2, split2, spec1, spec2; kwargs...)

# ## Planck Likelihood Specifics
# The Planck likelihood uses a specific choice of spectra and multipole ranges for those
# spectra. We provide some utility functions to retrieve a copy of the spectra order and
# the multipole minimum and maximum for those spectra.

plic_order() = (
    (:TT,"100","100"), (:TT,"143","143"), (:TT,"143","217"), (:TT,"217","217"),
    (:EE,"100","100"), (:EE,"100","143"), (:EE,"100","217"), (:EE,"143","143"),
    (:EE,"143","217"), (:EE,"217","217"),
    (:TE,"100","100"), (:TE,"100","143"), (:TE,"100","217"), (:TE,"143","143"),
    (:TE,"143","217"), (:TE,"217","217")
)

const plic_ellranges = Dict(
    (:TT, "100", "100") => (30, 1197),
    (:TT, "143", "143") => (30, 1996),
    (:TT, "143", "217") => (30, 2508),
    (:TT, "217", "217") => (30, 2508),
    (:EE, "100", "100") => (30, 999),
    (:EE, "100", "143") => (30, 999),
    (:EE, "100", "217") => (505, 999),
    (:EE, "143", "143") => (30, 1996),
    (:EE, "143", "217") => (505, 1996),
    (:EE, "217", "217") => (505, 1996),

    (:TE, "100", "100") => (30, 999),
    (:TE, "100", "143") => (30, 999),
    (:TE, "100", "217") => (505, 999),
    (:TE, "143", "143") => (30, 1996),
    (:TE, "143", "217") => (505, 1996),
    (:TE, "217", "217") => (505, 1996),

    (:TT, "100", "143") => (30, 999),   # not used
    (:TT, "100", "217") => (505, 999),  # not used
)

function get_plic_ellrange(spec::Symbol, freq1, freq2)
    if spec ∈ (:TE, :ET)
        if parse(Float64, freq1) > parse(Float64, freq2)
            freq1, freq2 = freq2, freq1
        end
        return plic_ellranges[:TE, freq1, freq2]
    end
    return plic_ellranges[spec, freq1, freq2]
end

# ## Signal Spectra
# The covariance matrix calculation and and signal simulations require an assumed signal
# spectra. We use the same foreground spectra as used in the ``\text{plic}`` likelihood.
# This returns dictionaries for signal and theory ``C_{\ell}`` in ``\mu\mathrm{K}``
# between two frequencies. The data is stored in the `plicref` directory in the config.

function signal_and_theory(freq1, freq2, config::Dict)
    likelihood_data_dir = joinpath(config["scratch"], "plicref")
    th = read_commented_header(joinpath(likelihood_data_dir,"theory_cl.txt"))
    fg = read_commented_header(joinpath(likelihood_data_dir,
        "base_plikHM_TTTEEE_lowl_lowE_lensing.minimum.plik_foregrounds"))

    for (spec, f1, f2) in plic_order()
        lmin, lmax = get_plic_ellrange(spec, f1, f2)
        const_val = fg[lmax-1,"$(spec)$(f1)X$(f2)"]
        ## constant foreground level after lmax -- there are fitting artifacts otherwise
        fg[(lmax-1):end, "$(spec)$(f1)X$(f2)"] .= const_val
    end

    ## loop over spectra and also fill in the flipped name
    freqs = ("100", "143", "217")
    specs = ("TT", "TE", "ET", "EE")
    for f1 in freqs, f2 in freqs, spec in specs
        if "$(spec)$(f1)X$(f2)" ∉ names(fg)
            if "$(reverse(spec))$(f2)X$(f1)" ∈ names(fg)
                fg[!, "$(spec)$(f1)X$(f2)"] = fg[!, "$(reverse(spec))$(f2)X$(f1)"]
            else
                fg[!, "$(spec)$(f1)X$(f2)"] = zeros(nrow(fg))
            end
        end
    end

    ap(v) = vcat([0., 0.], v)
    ell_fac = fg[!, "l"] .* (fg[!, "l"] .+ 1) ./ (2π);
    signal_dict = Dict{String,Vector{Float64}}()
    theory_dict = Dict{String,Vector{Float64}}()

    for XY₀ in specs  ## XY₀ is the spectrum to store
        f₁, f₂ = parse(Int, freq1), parse(Int, freq2)
        if f₁ <= f₂
            XY = XY₀
        else  ## swap what we're looking for, as fg data only has those cross-spectra
            XY = XY₀[2] * XY₀[1]
            f₁, f₂ = f₂, f₁
        end
        if XY == "ET"
            theory_cl_XY = th[!, "TE"] ./ (th[!, "L"] .* (th[!, "L"] .+ 1) ./ (2π))
        else
            theory_cl_XY = th[!, XY] ./ (th[!, "L"] .* (th[!, "L"] .+ 1) ./ (2π))
        end
        fg_cl_XY = fg[!, "$(XY)$(f₁)X$(f₂)"] ./ (fg[!, "l"] .* (fg[!, "l"] .+ 1) ./ (2π))

        signal_dict[XY₀] =  ap(theory_cl_XY .+ fg_cl_XY[1:2507])
        theory_dict[XY₀] =  ap(theory_cl_XY)
    end
    return signal_dict, theory_dict
end

# ## Planck ``\texttt{plic}`` Reference
# In various parts of the pipeline, we want to compare our results to the official 2018
# data release. These routines load them from disk. They're automatically downloaded to the
# `plicref` directory specified in the configuration TOML.

"""
    PlanckReferenceCov(plicrefpath::String)

Stores the spectra and covariances of the reference plic analysis.
"""
PlanckReferenceCov

struct PlanckReferenceCov{T}
    cov::Array{T,2}
    ells::Vector{Int}
    cls::Vector{T}
    keys::Vector{String}
    sub_indices::Vector{Int}
    key_index_dict::Dict{String,Int}
end

function PlanckReferenceCov(plicrefpath)
    ellspath = joinpath(plicrefpath, "vec_all_spectra.dat")
    clpath = joinpath(plicrefpath, "data_extracted.dat")
    covpath = joinpath(plicrefpath, "covmat.dat")
    keys = ["TT_100x100", "TT_143x143", "TT_143x217", "TT_217x217", "EE_100x100",
        "EE_100x143", "EE_100x217", "EE_143x143", "EE_143x217", "EE_217x217", "TE_100x100",
        "TE_100x143", "TE_100x217", "TE_143x143", "TE_143x217", "TE_217x217"]
    cov = inv(readdlm(covpath))
    ells = readdlm(ellspath)[:,1]
    cls = readdlm(clpath)[:,2]

    subarray_indices = collect(0:(size(cov,1)-2))[findall(diff(ells) .< 0) .+ 1] .+ 1
    sub_indices = [1, subarray_indices..., length(cls)+1]
    key_ind = 1:length(keys)
    key_index_dict = Dict(keys .=> key_ind)

    return PlanckReferenceCov{Float64}(cov, ells, cls, keys, sub_indices, key_index_dict)
end

"""
    get_subcov(pl::PlanckReferenceCov, spec1, spec2)

Extract the sub-covariance matrix corresponding to spec1 × spec2.

### Arguments:
- `pl::PlanckReferenceCov`: data structure storing reference covmat and spectra
- `spec1::String`: spectrum of form i.e. "TT_100x100"
- `spec2::String`: spectrum of form i.e. "TT_100x100"

### Returns:
- `Matrix{Float64}`: subcovariance matrix
"""
function get_subcov(pl::PlanckReferenceCov, spec1, spec2)
    i = pl.key_index_dict[spec1]
    j = pl.key_index_dict[spec2]
    return pl.cov[
        pl.sub_indices[i] : (pl.sub_indices[i + 1] - 1),
        pl.sub_indices[j] : (pl.sub_indices[j + 1] - 1),
    ]
end

"""
    extract_spec_and_cov(pl::PlanckReferenceCov, spec1)

Extract the reference ells, cl, errorbar, and sub-covariance block for a spectrum × itself.

### Arguments:
- `pl::PlanckReferenceCov`: data structure storing reference covmat and spectra
- `spec1::String`: spectrum of form i.e. "TT_100x100"

### Returns:
- `(ells, cl, err, this_subcov)`
"""
function extract_spec_and_cov(pl::PlanckReferenceCov, spec1)
    i = pl.key_index_dict[spec1]
    this_subcov = get_subcov(pl, spec1, spec1)
    ells = pl.ells[pl.sub_indices[i]:(pl.sub_indices[i + 1] - 1)]
    cl = pl.cls[pl.sub_indices[i]:(pl.sub_indices[i + 1] - 1)]
    err = sqrt.(diag(this_subcov))
    return ells, cl, err, this_subcov
end

max_forks_repo_name: sgiardie/PSpipe


text:
export ft_close

function ft_close(ft_handle :: Culong)
  ft_status = ccall((:FT_Close, d2xx),
                      Cuint,
                      (Culong,),
                      ft_handle)
  checkstatus(ft_status)
  return nothing
end


text:
using JLD
using Plots

x = load("ct_results.jld")
max_σ = x["max_σ"]
σ_vec1 = x["σ_vec"]
x = load("sim_results.jld")
max_σ_mat = x["max_σ_mat"]
σ_vec2 = x["σ_vec"]

histogram(σ_vec1, max_σ, bins = σ_vec1, title = "ODE Tsitouras", xlabel = "sigma", ylabel = "Freq")
histogram(σ_vec2, max_σ_mat[:, 1], bins = σ_vec2, title = "2 Subintervals", xlabel = "sigma", ylabel = "Freq")
histogram(σ_vec2, max_σ_mat[:, 2], bins = σ_vec2, title = "3 Subintervals", xlabel = "sigma", ylabel = "Freq")
histogram(σ_vec2, max_σ_mat[:, 3], bins = σ_vec2, title = "12 Subintervals", xlabel = "sigma", ylabel = "Freq")


text:
module UnstructuredGrids

include("Helpers.jl")
include("Kernels.jl")
include("Core.jl")
include("RefCellGallery.jl")
include("Factories.jl")
include("VTK.jl")

export Connections
export RefCell
export UGrid
export VERTEX
export connections
export coordinates
export writevtk
export list, ptrs
export generate_dual_connections
export generate_cell_to_faces
export find_cell_to_faces
export generate_grid_graph
export generate_full_grid_graph

using UnstructuredGrids.Core: Connections
using UnstructuredGrids.Core: UGrid
using UnstructuredGrids.Core: RefCell
using UnstructuredGrids.Core: list, ptrs
using UnstructuredGrids.Core: connections
using UnstructuredGrids.Core: coordinates
using UnstructuredGrids.Core: VERTEX
using UnstructuredGrids.VTK: writevtk
using UnstructuredGrids.Kernels: generate_dual_connections
using UnstructuredGrids.Kernels: generate_cell_to_faces
using UnstructuredGrids.Kernels: find_cell_to_faces
using UnstructuredGrids.Core: generate_grid_graph
using UnstructuredGrids.Core: generate_full_grid_graph

end # module UnstructuredGrids


text:
abstract type AbstractNewtonSolver{T} <: NonlinearSolver{T} end

@define newton_solver_variables begin
    x::Vector{T}
    J::Matrix{T}

    x₀::Vector{T}
    x₁::Vector{T}
    y₀::Vector{T}
    y₁::Vector{T}
    δx::Vector{T}
    δy::Vector{T}

    F!::FT
    Jparams::TJ

    linear::TL

    params::NonlinearSolverParameters{T}
    status::NonlinearSolverStatus{T}
end

function setInitialConditions!(s::AbstractNewtonSolver{T}, x₀::Vector{T}) where {T}
    s.x[:] = x₀
end

function computeJacobian(s::AbstractNewtonSolver)
    computeJacobian(s.x, s.J, s.Jparams)
end

function check_jacobian(s::AbstractNewtonSolver)
    println("Condition Number of Jacobian: ", cond(s.J))
    println("Determinant of Jacobian:      ", det(s.J))
    println("minimum(|Jacobian|):          ", minimum(abs.(s.J)))
    println("maximum(|Jacobian|):          ", maximum(abs.(s.J)))
    println()
end

function print_jacobian(s::AbstractNewtonSolver)
    display(s.J)
    println()
end


text:
macro test_zip(input...)
    n = length(input)
    x = Expr(:tuple, ntuple(i->gensym(), n)...)
    v = Expr(:tuple, map(esc, input)...)
    w = :(zip($(map(esc, input)...)))
    quote
    br = Any[]
    for $x in zip($v...)
        push!(br, $x)
    end
    mr = Any[]
    @itr for $x in $w
        push!(mr, $x)
    end
    @test br == mr
    end
end

macro test_enumerate(input)
    i = gensym()
    x = gensym()
    v = esc(input)
    quote
    br = Any[]
    for ($i,$x) in enumerate($v)
        push!(br, ($i,$x))
    end
    mr = Any[]
    @itr for ($i,$x) in enumerate($v)
        push!(mr, ($i,$x))
    end
    @test br == mr
    end
end

macro test_take(input, n)
    x = gensym()
    v = esc(input)
    quote
    br = Any[]
    for $x in take($v, $n)
        push!(br, $x)
    end
    mr = Any[]
    @itr for $x in take($v, $n)
        push!(mr, $x)
    end
    @test br == mr
    end
end

macro test_takestrict(input, n)
    x = gensym()
    v = esc(input)
    quote
    br = Any[]
    bfailed = false
    try
        for $x in takestrict($v, $n)
            push!(br, $x)
        end
    catch
        bfailed = true
    end

    mr = Any[]
    mfailed = false
    try
        @itr for $x in takestrict($v, $n)
            push!(mr, $x)
        end
    catch
        mfailed = true
    end
    @test br == mr
    @test bfailed == mfailed
    end
end

macro test_drop(input, n)
    x = gensym()
    v = esc(input)
    quote
    br = Any[]
    for $x in drop($v, $n)
        push!(br, $x)
    end
    mr = Any[]
    @itr for $x in drop($v, $n)
        push!(mr, $x)
    end
    @test br == mr
    end
end

macro test_chain(input...)
    x = gensym()
    v = Expr(:tuple, map(esc, input)...)
    w = :(chain($(map(esc, input)...)))
    quote
    br = Any[]
    for $x in chain($v...)
        push!(br, $x)
    end
    mr = Any[]
    @itr for $x in $w
        push!(mr, $x)
    end
    @test br == mr
    end
end


text:
using LinearAlgebra
using Statistics
using PyCall
using Plots
using NCFlex
using Optim
using JuLIP
using ASE
using FiniteDifferences

@pyimport ase.units as ase_units
@pyimport ase.lattice.cubic as ase_lattice_cubic
@pyimport matscipy.fracture_mechanics.crack as crack
@pyimport matscipy.fracture_mechanics.clusters as clusters
@pyimport matscipy.elasticity as elasticity

# Parameters

r_III = 100.0 # radius of region III. Other regions are derived from this
crack_surface = [1, 1, 1] # y-diretion, the surface which is opened
crack_front = [1, -1, 0] # z-direction, crack front line
relax_elasticity = false # if true, C_ij matrix computed with internal relaxation
relax_surface = false # if true, surface energy computed with internal relaxation

imsw = StillingerWeber(brittle=true)

unitcell = bulk(:Si, cubic=true)
variablecell!(unitcell)
set_calculator!(unitcell, imsw)

res = optimize(x -> energy(unitcell, x), x -> gradient(unitcell, x),
               dofs(unitcell), inplace=false)

##

# relaxed lattice constant
alat = unitcell.cell[1, 1]
@show alat
unitcell = bulk(:Si, cubic=true)
set_cell!(unitcell, alat * I(3))
set_calculator!(unitcell, imsw);

##

# 6x6 elastic constant matrix
C = voigt_moduli(unitcell)
C11, C12, C44 = cubic_moduli(C)
@show C11, C12, C44;
##

shift = 2.0

make_bulk() = Atoms(ASEAtoms(ase_lattice_cubic.Diamond(symbol="Si",
                            latticeconstant=alat,
                            directions=[[1,-1,0],[1,1,-2],[1,1,1]]) * (1, 1, 10)))
bulk_at = make_bulk()
set_calculator!(bulk_at, imsw)
surface = make_bulk() #copy(bulk_at) # must be a better way to do this in JuLIP
X = positions(surface) |> mat
X[3, :] .+= shift
set_positions!(surface, X)
wrap_pbc!(surface)
c = Matrix(surface.cell)
c[3, :] += [0.0, 0.0, 10.0]
set_cell!(surface, c)
fixedcell!(surface)
set_calculator!(surface, imsw)
area = norm(cross(bulk_at.cell[:, 1], bulk_at.cell[:, 2]))

if relax_surface
    rattle!(surface, 0.1)
    res = optimize(x -> energy(surface, x), x -> gradient(surface, x),
                dofs(surface), inplace=false)
    print(res)
end
γ = (energy(surface) - energy(bulk_at)) / (2 * area)

@show γ, γ / (ase_units.J / ase_units.m^2)

##

# build the final crystal in correct orientation
r_cut = cutoff(imsw) # cutoff distance for the potential
r_I = r_III - 3 * r_cut

n = [2 * floor(Int, (r_III + r_cut)/ alat),
     2 * floor(Int, (r_III + r_cut)/ alat) - 1, 1]
@show n

cryst = clusters.diamond("Si", alat, n,
                         crack_surface=crack_surface,
                         crack_front=crack_front)

# mark regions and convert to JuLIP
clust = clusters.set_regions(cryst, r_I, r_cut, r_III)
region = clust.arrays["region"]
clust = Atoms(ASEAtoms(clust))

X = positions(clust) |> mat

# object for computing the CLE displacments and gradients
C = voigt_moduli(C11, C12, C44)
Crot = rotate_elastic_moduli(C, rotation_matrix(crack_surface, crack_front))

# old Python code - we fix the rotated elastic constants
# to avoid the disrecpancy in Voigt conversion
crk = crack.CubicCrystalCrack(crack_surface, crack_front,
                              Crot=Crot)

# new Julia code
rac = RectilinearAnisotropicCrack(PlaneStrain(), C11, C12, C44,
                                   crack_surface, crack_front)

k_G1 = crk.k1g(γ)
k_G2 = k1g(rac, γ)
@assert abs(k_G1 - k_G2) < 1e-8

x0, y0, _ = diag(clust.cell) ./ 2

# displacement fields with both approaches
u1, v1 = crk.displacements(X[1, :], X[2, :], x0, y0, 1.0)
u2, v2 = displacements(rac, X[1, :] .- x0, X[2, :] .- y0)

@assert maximum(abs.(u1 - u2)) < 1e-8
@assert maximum(abs.(v1 - v2)) < 1e-8

# check gradient wrt finite differnces
u, ∇u = u_CLE(rac, clust, x0, y0)

@assert maximum((central_fdm(2, 1))(α -> u(1.0, α), 1.2) - ∇u(1.0, 1.2)) < 1e-6

X0 = copy(X)

##

U = [u(1.0, 0.0); zeros(length(clust))']

X = X0 + U
set_positions!(clust, X)

# scatter(X[1, :], X[2, :],
#         color=region, aspect_ratio=:equal, label=nothing)

@pyimport quippy.potential as quippy_potential
@pyimport atomistica

imswq = quippy_potential.Potential("IP SW", param_str="""
<SW_params n_types="1">
<per_type_data type="1" atomic_num="14" />
<per_pair_data atnum_i="14" atnum_j="14" AA="7.049556277" BB="0.6022245584"
      p="4" q="0" a="1.80" sigma="2.0951" eps="2.1675" />
<per_triplet_data atnum_c="14" atnum_j="14" atnum_k="14"
      lambda="42.0" gamma="1.20" eps="2.1675" />
</SW_params>
""")

r = sqrt.((X0[1,:] .- x0).^2 + (X0[2,:] .- y0).^2)
p = plot()
for (label, calculator) in zip(("JuLIP SW", "QUIP SW", "TersoffScr", "KumagaiScr"),
                                (imsw, ASECalculator(imswq), ASECalculator(atomistica.TersoffScr()),
                                ASECalculator(atomistica.KumagaiScr())))
    @show label
    set_calculator!(clust, calculator)
    scatter!(r, norm.(forces(clust)), yscale=:log10, xscale=:log10, label=label, legend=:bottomleft)
end
p


text:
# # File Loaders

function loadtoml(path::AbstractPath, globals)
    defaults = Dict{String,Any}(
        "deps" => Dict{String,String}(),
        "publish" => Dict(
            "ignore" => [],
            "theme" => "Publish.Themes.default",
            "config" => "Project.toml",
            "pages" => ["README.md"],
            "toc" => "toc.md",
            "modules" => [],
            "references" => [],
        ),
    )
    return rmerge(defaults, open(TOML.parse, path), globals)
end

function loadrefs(env)
    # Loading of references/bibliography. We can either load a `.bib` file for
    # use in LaTeX/PDF output, which turns on `citations = biblatex` for
    # control of how citations are printed out to LaTeX.
    #
    # Our other option is to print out basic citations using a JSON or TOML
    # source. TODO: write conversion from JSON/TOML to BIB and back so that we
    # can use either source for biblatex bibliographies.
    if haskey(env["publish"], "latex") && haskey(env["publish"]["latex"], "bibliography")
        env["publish"]["citations"] = "biblatex"
    else
        refs = env["publish"]["references"]
        if isa(refs, String) && isfile(refs)
            # TODO: support YAML and XML as well?
            env["publish"]["references"] =
                endswith(refs, ".toml") ? TOML.parsefile(refs) :
                endswith(refs, ".json") ? JSON.Parser.parsefile(refs) : []
        end
    end
    return env
end

rtree(f, d) = isdir(d) ? (d => [cd(()->rtree(f, x), d) for x in readdir(d) if f(x)]) : d
mktree(dir::Union{AbstractPath,AbstractString}, f=x->true) = FileTrees.maketree(rtree(f, abspath(dir)))

function loadtree(env::AbstractDict, p::AbstractPath)
    ## Custom ignore patterns from configuration are regular expressions.
    ignore = map(Regex, env["publish"]["ignore"])
    fn = function (path::AbstractString)
        ## Always ignore anything starting with '.' or '_'.
        startswith(path, ['.', '_']) && return true
        for each in ignore
            occursin(each, path) && return true
        end
        return false
    end
    ## Construct the initial FileTree manually since this avoids reading in
    ## huge directories that exceed open file limit on default OSX.
    tree = mktree(string(isfile(p) ? dirname(p) : p), !fn)
    return FileTrees.load(tree; lazy=LAZY[]) do file
        loadfile(env, joinpath(basename(tree), path(file)))
    end
end

"""
A dispatch type used to make file loading extensible by extension name.
"""
struct Extension{E} end
Extension(path::AbstractPath) = Extension{Symbol(extension(path))}()

"""
    loadfile(env, path)

Loads a file. Extended by `loadfile` methods that dispatch based on the
[extension](# "`Extension`") of the file.
"""
loadfile(env::AbstractDict, path::AbstractPath) = loadfile(Extension(path), env, path)

loadfile(::Extension{:md}, env, path) = open(path) do io
    parser = init_markdown_parser(env["publish"])
    load_markdown(io, parser)
end

function loadfile(::Extension{:jl}, env, path)
    io = IOBuffer()
    code = String[]
    state = :text
    ## Helper function the reduces code duplication below.
    code_block_helper = function (state)
        if state === :code
            first = findfirst(l -> any(!isspace, l), code)
            last  = findlast(l -> any(!isspace, l), code)
            (first === last === nothing) || join(io, code[first:last])
            empty!(code)
            println(io, CODE_FENCE, "\n")
        end
    end
    open(path) do handle
        for line in eachline(handle)
            occursin(r"#src\b", line) && continue
            m = match(r"^(\s*)([#]*)(.*)", line)
            if m !== nothing
                ws, comments, rest = m[1], m[2], m[3]
                count = length(comments)
                if count == 1
                    ## Remove single whitespace after the comment.
                    line = chop(rest; head=1, tail=0)
                    code_block_helper(state)
                    println(io, line)
                    state = :text
                else
                    ## Start a new code block.
                    state === :text && println(io, CODE_FENCE, "julia")
                    push!(code, string(ws, count === 0 ? "" : '#'^(count-1), rest, '\n'))
                    state = :code
                end
            end
        end
    end
    ## Clean up last code block.
    code_block_helper(state)
    parser = init_markdown_parser(env["publish"])
    return load_markdown(io, parser)
end

function loadfile(::Extension{:ipynb}, env, path)
    dict = open(JSON.parse, path)
    io = IOBuffer()
    if haskey(dict, "cells")
        for cell in dict["cells"]
            if haskey(cell, "cell_type")
                type = cell["cell_type"]
                source = get(cell, "source", "")
                if type == "markdown"
                    join(io, source)
                    println(io)
                elseif type == "code"
                    println(io, CODE_FENCE, "julia")
                    join(io, source)
                    println(io)
                    println(io, CODE_FENCE)
                end
            end
        end
    end
    parser = init_markdown_parser(env["publish"])
    return load_markdown(io, parser)
end

loadfile(::Extension, env, path) = read(path)

"""
Loads "virtual" files for each docstring that is defined within the given
project. Merges these files into the given `tree` as well as appending them to
`pages`.
"""
function loaddocs(tree::FileTree, env::AbstractDict, pages::Vector)
    docs_dir, docs_index = "docstrings", "docstrings.md"
    roots = findmodules(env)
    docs = DataStructures.SortedDict{String,Tuple{Module,Docs.Binding,Docs.MultiDoc}}()
    for mod in visible_modules(env)
        if Base.moduleroot(mod) in roots
            for (k, v) in Docs.meta(mod)
                docs["$k.md"] = (mod, k, v)
            end
        end
    end
    isempty(docs) && return tree, pages # Bail early if no docstrings are found.
    append!(pages, sort!([joinpath(Path(docs_dir), k) for k in keys(docs)]; by=string))
    dtree = maketree(basename(tree) => [docs_index, docs_dir => keys(docs)])
    dtree = FileTrees.load(dtree; lazy=LAZY[]) do file
        name = basename(path(file))
        io = IOBuffer()
        visibility(binding) = Base.isexported(binding.mod, binding.var) ? "public" : "private"
        if name == docs_index
            println(io, "{#docstring-index}")
            println(io, "| Name | Module | Visibility | Category |")
            println(io, "|------|--------|------------|----------|")
            for (k, (mod, bind, doc)) in docs
                vis, cat = visibility(bind), categorise(bind)
                println(io, "| [`$(bind.var)`](docstrings/$k) | `$mod` | `$vis` | `$cat` |")
            end
        else
            mod, bind, doc = docs[name]
            println(io, "```{=html}\n<div class='docs' id='$(bind.var)'>\n```") # TODO: raw latex.
            println(io, "`$(visibility(bind))` `$(bind.var)` --- `$(categorise(bind))`")
            for (n, sig) in enumerate(doc.order)
                println(io)
                println(io, "```{=html}\n<div class='doc' id='$n'>\n```")
                printdoc(io, doc.docs[sig])
                println(io, "```{=html}\n</div>\n```")
            end
            println(io, "```{=html}\n</div>\n```")
        end
        parser = init_markdown_parser(env["publish"])
        return load_markdown(io, parser)
    end
    return merge(tree, dtree), pages
end

"""
    loadpages(tree, env)

Finds all files defined by the project's table of contents and loads them into
the `tree`.
"""
function loadpages(tree::FileTree, env::AbstractDict)
    toc = env["publish"]["toc"]
    tree = try_touch(tree, toc) do
        io = IOBuffer()
        for page in env["publish"]["pages"]
            path = Path(page)
            println(io, "  - [$(filename(path))]($path)")
        end
        return load_markdown(io)
    end
    pages = []
    for (node, enter) in exec(tree[toc][])
        if enter && node.t isa CommonMark.Link
            push!(pages, Path(node.t.destination))
        end
    end
    return tree, pages
end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libopk"], :libopk),
    LibraryProduct(prefix, ["libcobyla"], :libcobyla),
    LibraryProduct(prefix, ["libbobyqa"], :libbobyqa),
    LibraryProduct(prefix, ["libnewuoa"], :libnewuoa),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/emmt/OptimPack/releases/download/v3.1.0"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc)                  => ("$bin_prefix/OptimPack.v3.1.0.aarch64-linux-gnu.tar.gz",          "9b3d559c9021c58b710ee240ef8f2e6a12b935517c72ac1e9f73cfb924999d10"),
    Linux(:aarch64, libc=:musl)                   => ("$bin_prefix/OptimPack.v3.1.0.aarch64-linux-musl.tar.gz",         "8c427f5574fe62d1c393568b82c75da151058ea82580a77de2fe2f4430943f7d"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf) => ("$bin_prefix/OptimPack.v3.1.0.arm-linux-gnueabihf.tar.gz",        "dbcb477fb7cceeb46c08fb21ac6a0367bdba876b965f72833580e2834a672627"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf)  => ("$bin_prefix/OptimPack.v3.1.0.arm-linux-musleabihf.tar.gz",       "bddb5f9ca1f061f9a8c5202cb9807c028cab59bac61f34f092c626ef0a756041"),
    Linux(:i686, libc=:glibc)                     => ("$bin_prefix/OptimPack.v3.1.0.i686-linux-gnu.tar.gz",             "41c293fe0de49d37306ae6219e1de63a1d633ea451cbe62c0d47142759c7a5d9"),
    Linux(:i686, libc=:musl)                      => ("$bin_prefix/OptimPack.v3.1.0.i686-linux-musl.tar.gz",            "3691a1179ea9cea04eaf1e6afe41c5d6f85d02d75da05adc0d021f500ae28c1d"),
    Windows(:i686)                                => ("$bin_prefix/OptimPack.v3.1.0.i686-w64-mingw32.tar.gz",           "ace0885f9c0c614df861688821ee5fbef3de0a6db9d616a0be59970a44ae4151"),
    Linux(:powerpc64le, libc=:glibc)              => ("$bin_prefix/OptimPack.v3.1.0.powerpc64le-linux-gnu.tar.gz",      "7d7e912fb028e129fdcc2a91768b942e4788239b2bcc9421ed11affd3bd7a918"),
    MacOS(:x86_64)                                => ("$bin_prefix/OptimPack.v3.1.0.x86_64-apple-darwin14.tar.gz",      "91914817902f2ce48ff7ecabbdfc303d51e2a65926fa1dd0823579257bc545fa"),
    Linux(:x86_64, libc=:glibc)                   => ("$bin_prefix/OptimPack.v3.1.0.x86_64-linux-gnu.tar.gz",           "4738864c8166b69f128ea1497a53cdc46098e89472b14715d31b25a5a2a5164b"),
    Linux(:x86_64, libc=:musl)                    => ("$bin_prefix/OptimPack.v3.1.0.x86_64-linux-musl.tar.gz",          "2acbcbda8478a7b63ba96d1d8d29848684555a4de5ed2912240f7c9660039142"),
    FreeBSD(:x86_64)                              => ("$bin_prefix/OptimPack.v3.1.0.x86_64-unknown-freebsd11.1.tar.gz", "268f36dddd8e6e71a4ef6f2badb9329315750e7b26fc5f5f27193ccab0b80015"),
    Windows(:x86_64)                              => ("$bin_prefix/OptimPack.v3.1.0.x86_64-w64-mingw32.tar.gz",         "2a914c7bcc0d05875569878b9d0188012d591026239263d27f6fb46b17484c66"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMCoupling.jl/blob/master/LICENSE

# using JuliaFEM
using FEMCoupling

# Validation of kinematic coupling: we have 1x1x1 brick which is fixed at z=0
# Reference node is 1 away from the brick in z direction and twist is applied

# This test is decoupled from JuliaFEM by assembling stiffness matrix and
# saving it as a test resource. Stiffness matrix can be read back using
# readdlm.

X = Dict(
    1 => [0.0, 0.0, 0.0],
    2 => [1.0, 0.0, 0.0],
    3 => [1.0, 1.0, 0.0],
    4 => [0.0, 1.0, 0.0],
    5 => [0.0, 0.0, 1.0],
    6 => [1.0, 0.0, 1.0],
    7 => [1.0, 1.0, 1.0],
    8 => [0.0, 1.0, 1.0],
    9 => [0.5, 0.5, 2.0])

element = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8])
update!(element, "geometry", X)
update!(element, "youngs modulus", 2.080e5)
update!(element, "poissons ratio", 0.3)
update!(element, "density", 7.80e-9)

# problem = Problem(Elasticity, "1 hex8 element", 3)
# add_elements!(problem, [element])
# assemble!(problem, 0.0)
# K = full(problem.assembly.K)
# fn = Pkg.dir("FEMCoupling", "test", "test_dcoupling_3d", "stiffness_matrix")
# writedlm(fn, K)

bc_element = Element(Quad4, [1, 2, 3, 4])
update!(bc_element, "geometry", X)
update!(bc_element, "displacement 1", 0.0)
update!(bc_element, "displacement 2", 0.0)
update!(bc_element, "displacement 3", 0.0)

# bc = Problem(Dirichlet, "fixed", 3, "displacement")
# add_elements!(bc, [bc_element])

# Distributed coupling

coupling_node_ids = [5, 6, 7, 8]
coupling_nodes = [Element(Poi1, [j]) for j in coupling_node_ids]
reference_node_id = 9
reference_node = Element(Poi1, [reference_node_id])
update!(coupling_nodes, "geometry", X)
update!(reference_node, "geometry", X)
update!(reference_node, "point moment 3", 1.0e3)

coupling = Problem(Coupling, "brick", 3, "displacement")
add_coupling_nodes!(coupling, coupling_nodes)
add_reference_node!(coupling, reference_node)

# Analysis

# step = Analysis(Nonlinear)
# add_problems!(step, [problem, bc, coupling])
# xdmf = Xdmf("3d_coupling_results"; overwrite=true)
# add_results_writer!(step, xdmf)
# run!(step)
# close(xdmf.hdf)

K = readdlm(Pkg.dir("FEMCoupling", "test", "test_dcoupling_3d", "stiffness_matrix"))
assemble!(coupling, 0.0)
f = full(sparse(coupling.assembly.f, 24, 1))
free_dofs = (4*3+1):(8*3)
u = zeros(3*8)
u[free_dofs] = K[free_dofs, free_dofs] \ f[free_dofs]
d = Dict(j => u[3*(j-1)+1:3*(j-1)+3] for j=1:8)

# Results

using Base.Test

# Because the twisting is symmetric, we should have equal magnitude of displacement
# in all free nodes
@test isapprox(norm(d[5]), norm(d[6]))
@test isapprox(norm(d[5]), norm(d[7]))
@test isapprox(norm(d[5]), norm(d[8]))

# Comparing calculated results with ABAQUS results. The model used to get
# results is given in directory test_coupling_3d/model.inp

# for i=5:8
#     println("node $i: $(u[i])")
# end

a = 3.75e-2
@test isapprox(d[5], [ a, -a, 0.0])
@test isapprox(d[6], [ a,  a, 0.0])
@test isapprox(d[7], [-a,  a, 0.0])
@test isapprox(d[8], [-a, -a, 0.0])


text:
#### Riemann Matrix
import Base: size
export Riemann

struct Riemann{Int} <: AbstractMatrix{Int}
    n::Int
end # immutable

# Define its size

size(A::Riemann, dim::Integer) = A.n
size(A::Riemann)= size(A,1), size(A,1)

# Index into a Riemann
function getindex(A::Riemann,i::Integer,j::Integer)
#    return (i+1)%(j+1)==0 ? i : -1
    if i<=j && (j+1)%(i+1)==0
        return i
    else
        return -1
    end
end # getindex

# Dense version of Riemann
Matrix(A::Riemann) =[A[i,j] for i=1:size(A,1), j=1:size(A,2)]


text:
mutable struct PartitionTS <: BanditPolicy
    t::Int64
    arms::Vector{PartitionModel}
    retrain::Vector{Int64}
    traj::Trajectory

    # Model parameters
    lims::Matrix{Float64}
    Jmax::Int64
    p::Int64
    pconst::Function
    dconst::Int64
    normalise::Bool
    shape::Float64
    scale::Float64
    prior_init::Function
    Kmax::Int64
    verbose::Bool
end

function PartitionTS(
    A::Int64,
    retrain::Vector{Int64},
    traj::Trajectory,
    lims::Matrix{Float64},
    Jmax::Int64,
    p::Int64,
    pconst::Function,
    dconst::Int64,
    normalise::Bool;
    shape::Float64=0.001,
    scale::Float64=0.001,
    prior_init::Function=identity_prior,
    Kmax::Int64=300,
    verbose::Bool=false,
)
    d = size(lims, 1)
    t = 0
    basis = [tpbasis(d, 0)]
    arms = [
        PartitionModel(
            basis,
            Partition(lims),
            normalise;
            shape=shape,
            scale=scale,
            prior_init=prior_init,
        ) for _ in 1:A
    ]
    return PartitionTS(
        t,
        arms,
        retrain,
        traj,
        lims,
        Jmax,
        p,
        pconst,
        dconst,
        normalise,
        shape,
        scale,
        prior_init,
        Kmax,
        verbose,
    )
end

function update!(pol::PartitionTS, x::Vector{Float64}, a::Int64, r::Float64)
    pol.t += 1
    if pol.t != pol.traj.t
        throw(ErrorException("Policy out of sync with trajectory"))
    end

    if pol.t in pol.retrain
        # Retrain each model from scratch
        for i in 1:length(pol.arms)
            Xfull = pol.traj.X[1:pol.t, :]
            rfull = pol.traj.r[1:pol.t]
            afull = pol.traj.a[1:pol.t]
            Xarm = Xfull[afull .== i, :]
            rarm = rfull[afull .== i]
            pol.arms[i] = partition_search(
                Xarm,
                rarm,
                pol.lims,
                pol.Jmax,
                pol.p,
                pol.pconst,
                pol.dconst,
                pol.normalise;
                shape=pol.shape,
                scale=pol.scale,
                prior_init=pol.prior_init,
                Kmax=pol.Kmax,
                verbose=pol.verbose
            )
        end
    else
        fit!(pol.arms[a], x, r)
    end
end

function choose(x::Vector{Float64}, pol::PartitionTS)
    μ = Vector{Float64}(undef, length(pol.arms))
    for a in 1:length(pol.arms)
        shape = pol.arms[a].shape
        scale = pol.arms[a].scale
        varsim = rand(InverseGamma(shape, scale))

        # Now simulate coefficients
        k = locate(pol.arms[a].P, x)
        coeff = pol.arms[a].models[k].post.coeff
        Σ = varsim * pol.arms[a].models[k].post.cov
        coeffsim = rand(MvNormal(coeff, Σ))
        ψ = expand(x, pol.arms[a].models[k].basis, pol.arms[a].models[k].lims, pol.normalise)
        μ[a] = coeffsim' * ψ
    end
    return findmax(μ)[2]
end


text:
import MitosisStochasticDiffEq as MSDE
using Mitosis
using Test, Random
using LinearAlgebra

@testset "backward information filtering tests" begin
  # define SDE function
  f(u,p,t) = p[1]*u .+ p[2]
  g(u,p,t) = p[3]

  # time span
  tstart = 0.0
  tend = 1.0
  dt = 0.001
  trange = tstart:dt:tend

  # set of linear parameters Eq.~(2.2)
  B, β, σ̃ = -0.1, 0.2, 1.3
  plin = [B, β, σ̃]

  # initial values for ODE
  mynames = (:logscale, :μ, :Σ);
  c = 0.0
  ν = 0.0
  P = 10.0
  myvalues = [c, ν, P]
  NT = NamedTuple{mynames}(myvalues)

  # covariance filter
  kernel = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, plin)
  message, solend = MSDE.backwardfilter(kernel, NT)

  # information filter
  H = inv(P)
  F = H*ν
  myvalues2 = [c, F, H]
  NT2 = NamedTuple{mynames}(myvalues2)
  message2, solend2 = MSDE.backwardfilter(kernel, NT2, filter=MSDE.InformationFilter())

  @test !isapprox(solend, solend2, rtol=1e-1)
  invH = inv(solend2[2])
  invHF = invH*solend2[1]
  @test isapprox(solend[1], invHF, rtol=1e-3)
  @test isapprox(solend[2], invH, rtol=1e-3)
  @test isapprox(solend[3], solend2[3], rtol=1e-10)

  # multivariate tests
  dim = 5
  Random.seed!(123)
  c = randn()
  ν = randn(dim)
  P = randn(dim,dim)
  myvalues = [c, ν, P]
  NT = NamedTuple{mynames}(myvalues)
  H = inv(P)
  F = H*ν
  myvalues2 = [c, F, H]
  NT2 = NamedTuple{mynames}(myvalues2)

  # covariance filter
  kernel = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, plin)
  message, solend = MSDE.backwardfilter(kernel, NT)

  # information filter
  message2, solend2 = MSDE.backwardfilter(kernel, NT2,
    filter=MSDE.InformationFilter())

  @test !isapprox(solend, solend2, rtol=1e-1)
  invH = inv(solend2.x[2])
  invHF = invH*solend2.x[1]
  @test isapprox(solend.x[1], invHF, rtol=1e-3)
  @test isapprox(solend.x[2], invH, rtol=1e-2)
  @test isapprox(solend.x[3], solend2.x[3], rtol=1e-10)

  # test inplace version
  message2, solend2 = MSDE.backwardfilter(kernel, NT2,
    filter=MSDE.InformationFilter(), inplace=true)

  @test !isapprox(solend, solend2, rtol=1e-1)
  invH = inv(solend2.x[2])
  invHF = invH*solend2.x[1]
  @test isapprox(solend.x[1], invHF, rtol=1e-3)
  @test isapprox(solend.x[2], invH, rtol=1e-2)
  @test isapprox(solend.x[3], solend2.x[3], rtol=1e-10)
end
#
# # define SDE function
# f(u,p,t) = p[1]*u .+ p[2]
# g(u,p,t) = p[3]
#
# # time span
# tstart = 0.0
# tend = 1.0
# dt = 0.001
# trange = tstart:dt:tend
#
# m = 3 # some number of Brownian processes
# Random.seed!(123)
# plin = [randn(dim,dim), randn(dim), randn(dim,m)] # B, β, σtil
#
# # multivariate tests
# dim = 5
# Random.seed!(123)
# c = randn()
# ν = randn(dim)
# P = randn(dim,dim)
# myvalues = [c, ν, P]
# NT = NamedTuple{mynames}(myvalues)
# H = inv(P)
# F = H*ν
# myvalues2 = [c, F, H]
# NT2 = NamedTuple{mynames}(myvalues2)
#
# # covariance filter
# kernel = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, plin)
# message, solend = MSDE.backwardfilter(kernel, NT)
#
# # information filter
# message2, solend2 = MSDE.backwardfilter(kernel, NT2,
#   filter=MSDE.InformationFilter())
#
#
# @test !isapprox(solend, solend2, rtol=1e-1)
# invH = inv(solend2.x[2])
# invHF = invH*solend2.x[1]
# @test isapprox(solend.x[1], invHF, rtol=1e-3)
# @test isapprox(solend.x[2], invH, rtol=1e-2)
# @test isapprox(solend.x[3], solend2.x[3], rtol=1e-10)

# test inplace version
# message2, solend2 = MSDE.backwardfilter(kernel, NT2, inplace=true,
#   filter=MSDE.InformationFilter())
#
# @test !isapprox(solend, solend2, rtol=1e-1)
# invH = inv(solend2.x[2])
# invHF = invH*solend2.x[1]
# @test isapprox(solend.x[1], invHF, rtol=1e-3)
# @test isapprox(solend.x[2], invH, rtol=1e-2)
# @test isapprox(solend.x[3], solend2.x[3], rtol=1e-10)


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

"""
    ydn2md(year, day) -> date

### Purpose ###

Convert from year and day number of year to a date.

### Explanation ###

Returns the date corresponding to the `day` of `year`.

### Arguments ###

* `year`: the year, as an integer.
* `day`: the day of `year`, as an integer.

### Output ###

The date, of `Date` type, of \$\\text{day} - 1\$ days after January 1st of
`year`.

### Example ###

Find the date of the 60th and 234th days of the year 2016.

```jldoctest
julia> using AstroLib

julia> ydn2md.(2016, [60, 234])
2-element Array{Dates.Date,1}:
 2016-02-29
 2016-08-21
```

### Note ###

`ymd2dn` converts from a date to day of the year.
"""
function ydn2md(year::Integer, day::Integer)
    return Dates.firstdayofyear(Date(year)) + Dates.Day(day - 1)
end


text:
@testset "notice.jl" begin
    # notice.jl is tested in control.jl and process.jl.
end


text:
using Test
using CellListsBenchmarks

n = 100
d = 2
r = 0.1
seed = 1
iterations = 1
seconds = 1.0
Benchmark("cell_list_serial", n, d, r, seed, iterations, seconds)
Benchmark("cell_list_threads", n, d, r, seed, iterations, seconds)
Benchmark("brute_force", n, d, r, seed, iterations, seconds)
Benchmark("cell_lists", n, d, r, seed, iterations, seconds)
Benchmark("cell_lists", n, d, Float32(r), seed, iterations, seconds)
Benchmark("near_neighbors_serial", n, d, r, seed, iterations, seconds)
Benchmark("near_neighbors_threads", n, d, r, seed, iterations, seconds)
@test true


text:
module Benchmark
  using CxxWrap, MPI
  import ..libjltrilinos

  @wrapmodule(libjltrilinos, :register_benchmark)
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

module InverseDistanceWeighting

using GeoStatsBase

using NearestNeighbors
using StaticArrays
using Distances

import GeoStatsBase: solve

export InvDistWeight

"""
    InvDistWeight(var₁=>param₁, var₂=>param₂, ...)

Inverse distance weighting estimation solver.

## Parameters

* `neighbors` - Number of neighbors (default to all data locations)
* `distance`  - A distance defined in Distances.jl (default to Euclidean()
"""
@estimsolver InvDistWeight begin
  @param neighbors = nothing
  @param distance = Euclidean()
end

function solve(problem::EstimationProblem, solver::InvDistWeight)
  # retrieve problem info
  pdata = data(problem)
  pdomain = domain(problem)

  # result for each variable
  μs = []; σs = []

  for (var,V) in variables(problem)
    # get user parameters
    if var ∈ keys(solver.params)
      varparams = solver.params[var]
    else
      varparams = InvDistWeightParam()
    end

    # get valid data for variable
    X, z = valid(pdata, var)

    # number of data points for variable
    ndata = length(z)

    @assert ndata > 0 "estimation requires data"

    # allocate memory
    varμ = Vector{V}(undef, npoints(pdomain))
    varσ = Vector{V}(undef, npoints(pdomain))

    # fit search tree
    kdtree = KDTree(X, varparams.distance)

    # keep track of estimated locations
    estimated = falses(npoints(pdomain))

    # consider data locations as already estimated
    for (loc, datloc) in datamap(problem, var)
      estimated[loc] = true
      varμ[loc] = pdata[datloc,var]
      varσ[loc] = zero(V)
    end

    # determine number of nearest neighbors to use
    k = varparams.neighbors == nothing ? ndata : varparams.neighbors

    @assert k ≤ ndata "number of neighbors must be smaller or equal to number of data points"

    # pre-allocate memory for coordinates
    coords = MVector{ndims(pdomain),coordtype(pdomain)}(undef)

    # estimation loop
    for location in LinearPath(pdomain)
      if !estimated[location]
        coordinates!(coords, pdomain, location)

        idxs, dists = knn(kdtree, coords, k)

        weights = one(V) ./ dists
        weights /= sum(weights)

        values = view(z, idxs)

        varμ[location] = sum(weights[i]*values[i] for i in eachindex(values))
        varσ[location] = minimum(dists)
      end
    end

    push!(μs, var => varμ)
    push!(σs, var => varσ)
  end

  EstimationSolution(pdomain, Dict(μs), Dict(σs))
end

end


text:
import LinearAlgebra.det
import Statistics.var
import Statistics.std
using Printf
using StatsFuns
#using SortingAlgorithms
#using Distributions
using StatsBase
#using SharedArrays

################################################################################
# Gibbs algorithms
################################################################################

"""
Gibbs sampler using Kalman filter for the OU SDEMEM model.
"""
function gibbs_exact(R::Int,
                     y::Vector,
                     dt::Real,
                     Σ_i_σ_ϵ::Real,
                     Σ_i_ϕ::Array,
                     γ_ϕ_0::Real,
                     γ_σ_ϵ_0::Real,
                     μ_i_ϕ::Array,
                     μ_i_σ_ϵ::Real,
                     α_star_ϕ::Real,
                     α_star_σ_ϵ::Real,
                     log_λ_i_ϕ::Vector,
                     log_λ_i_σ_ϵ::Real,
                     update_interval::Real,
                     start_update::Real,
                     α_power::Real,
                     startval_ϕ::Array,
                     startval_log_σ_ϵ::Real,
                     startval_η::Vector,
                     prior_parameters_η::Array,
                     prior_parameters_σ_ϵ::Array)

    println("Starting Gibbs sampler using Kalman.")

    # dim for problem
    M = length(y)

    # initilization
    chain_ϕ = zeros(M, 3, R)
    chain_log_σ_ϵ = zeros(1, R)
    chain_η = zeros(6, R)
    ϕ_star = zeros(size(startval_ϕ))
    accept_vec = zeros(3,R)
    accept = false

    loglik_vec = zeros(R)

    α_prob_accept_σ_ϵ = 0.
    α_prob_accept_ϕ = zeros(M)

    # first iteration
    chain_ϕ[:,:,1] = startval_ϕ
    chain_log_σ_ϵ[1] = startval_log_σ_ϵ
    print_interval = 1000

    # set star values for η
    chain_η[:,1] = startval_η

    #=
    for j = 1:3

        # get paramaters
        μ_0_j, M_0_j, α_j, β_j = prior_parameters_η[j,:]

        chain_η[j,1] = μ_0_j # set start value to mean
        chain_η[j+3,1] = (α_j - 1)/β_j

    end
    =#

    accept_vec[:,1] = [M;1;3] # we accapt the start values

    # compute likelihood value
    loglik_old_y = kalman(y, exp(startval_log_σ_ϵ), startval_ϕ, dt)

    loglik_vec[1] = sum(loglik_old_y)

    # main loop
    for r = 2:R

        # first gibbs stage, update \phi
        for m = 1:M
            ϕ_star[m,:] = chain_ϕ[m,:,r-1] + cholesky(exp(log_λ_i_ϕ[m])*Σ_i_ϕ[m]).U*randn(3)
        end

        loglik_star_y = kalman(y, exp(chain_log_σ_ϵ[r-1]), ϕ_star, dt)

        loglik_old_ϕ = calc_loglik_ϕ(chain_ϕ[:,:,r-1], chain_η[:,r-1])
        loglik_star_ϕ = calc_loglik_ϕ(ϕ_star, chain_η[:,r-1])

        for m = 1:M

            # TODO remove jacobian!
            #jacobian_old =  sum(chain_ϕ[m,:,r-1])
            #jacobian_star = sum(ϕ_star[m,:])

            #log_α = (loglik_star_y[m] + loglik_star_ϕ[m] + jacobian_star) - (loglik_old_y[m] + loglik_old_ϕ[m] + jacobian_old)
            log_α = (loglik_star_y[m] + loglik_star_ϕ[m]) - (loglik_old_y[m] + loglik_old_ϕ[m])

            α_prob_accept_ϕ[m] = min(1, exp(log_α))

            # correct if we have  NaNs
            if any(isnan,ϕ_star[m,:]) == true || loglik_star_y[m] == NaN
                log_α = log(0)
                α_prob_accept_ϕ[m] = min(1, exp(log_α))
            end

            accept = log(rand()) < log_α

            if accept
                chain_ϕ[m,:,r] = ϕ_star[m,:]
                accept_vec[1,r] = accept_vec[1,r] +1
                loglik_old_y[m] = loglik_star_y[m]
            else
                chain_ϕ[m,:,r] = chain_ϕ[m,:,r-1]
            end

        end

        # secound gibbs stage, update σ_ϵ
        log_σ_ϵ_star = chain_log_σ_ϵ[r-1] + sqrt(exp(log_λ_i_σ_ϵ)*Σ_i_σ_ϵ)*randn()

        # TODO remove jacobian!
        #jacobian_old =  chain_log_σ_ϵ[r-1]
        #jacobian_star = log_σ_ϵ_star

        loglik_star_y = kalman(y, exp(log_σ_ϵ_star), chain_ϕ[:,:,r], dt)

        prior_old_σ_ϵ = normlogpdf(prior_parameters_σ_ϵ[1],1/prior_parameters_σ_ϵ[2],chain_log_σ_ϵ[r-1])
        prior_star_σ_ϵ = normlogpdf(prior_parameters_σ_ϵ[1],1/prior_parameters_σ_ϵ[2],log_σ_ϵ_star)

        #log_α = (sum(loglik_star_y) + prior_star_σ_ϵ + jacobian_star) - (sum(loglik_old_y) + prior_old_σ_ϵ + jacobian_old)
        log_α = (sum(loglik_star_y) + prior_star_σ_ϵ) - (sum(loglik_old_y) + prior_old_σ_ϵ)

        α_prob_accept_σ_ϵ = min(1, exp(log_α))

        # correct if we have  NaNs
        if any(isnan,loglik_star_y) == true || log_σ_ϵ_star == NaN
            log_α = log(0)
            α_prob_accept_σ_ϵ = min(1, exp(log_α))
        end

        accept = log(rand()) < log_α

        if accept
            chain_log_σ_ϵ[r] = log_σ_ϵ_star
            accept_vec[2,r] = accept_vec[2,r] + 1
            loglik_old_y = loglik_star_y
        else
            chain_log_σ_ϵ[r] = chain_log_σ_ϵ[r-1]
        end

        # third stage, update η
        for j = 1:3

            # get paramaters
            μ_0_j, M_0_j, α_j, β_j = prior_parameters_η[j,:]

            ϕ_j = chain_ϕ[:,j,r-1]
            ϕ_j_mean = 1/M*sum(ϕ_j)

            α = α_j + div(M,2)
            β = β_j + 0.5*sum((ϕ_j .- ϕ_j_mean).^2) + (M*M_0_j)/(2*(M+M_0_j))*(ϕ_j_mean-μ_0_j)^2

            τ_j = rand_gamma(α)
            τ_j = τ_j/β

            μ = (M*τ_j*ϕ_j_mean)/(M*τ_j + M_0_j*τ_j) + (M_0_j*τ_j*μ_0_j)/(M*τ_j + M_0_j*τ_j)
            τ = M*τ_j + M_0_j*τ_j

            μ_j = μ + sqrt(1/τ)*randn()

            chain_η[j,r] = μ_j
            chain_η[j+3,r] = τ_j

            accept_vec[3,r] = accept_vec[3,r] +1

        end

        loglik_vec[r] = sum(loglik_old_y)

        print_progress(r,print_interval,accept_vec,M)

        log_λ_i_σ_ϵ,μ_i_σ_ϵ,Σ_i_σ_ϵ = adaptive_update_σ_ϵ(r,update_interval,start_update,log_λ_i_σ_ϵ,μ_i_σ_ϵ,Σ_i_σ_ϵ,γ_σ_ϵ_0,α_power,chain_log_σ_ϵ,α_prob_accept_σ_ϵ)
        adaptive_update_ϕ!(r,update_interval,start_update,log_λ_i_ϕ,μ_i_ϕ,Σ_i_ϕ,γ_ϕ_0,α_power,M,chain_ϕ,α_prob_accept_ϕ)

    end

    println("Ending Gibbs sampler using Kalman.")

    return chain_ϕ, chain_log_σ_ϵ, chain_η, accept_vec, loglik_vec

end

"""
Gibbs sampler using CPMMH for the OU SDEMEM model.

Gibbs scheme #2 where we integrate out the random numbers.
"""
function gibbs_cpmmh(R::Int,
                     y::Array,
                     dt::Real,
                     Σ_i_σ_ϵ::Real,
                     Σ_i_ϕ::Array,
                     γ_ϕ_0::Real,
                     γ_σ_ϵ_0::Real,
                     μ_i_ϕ::Array,
                     μ_i_σ_ϵ::Real,
                     α_star_ϕ::Real,
                     α_star_σ_ϵ::Real,
                     log_λ_i_ϕ::Vector,
                     log_λ_i_σ_ϵ::Real,
                     update_interval::Real,
                     start_update::Real,
                     α_power::Real,
                     startval_ϕ::Array,
                     startval_log_σ_ϵ::Real,
                     startval_η::Vector,
                     prior_parameters_η::Array,
                     prior_parameters_σ_ϵ::Array,
                     nbr_particles::Vector,
                     T_vec::Vector,
                     ρ::Real,
                     particel_filter::Function)

    println("Starting Gibbs sampler using CPMMH.")

    # dim for problem
    M = length(y)

    # initilization
    chain_ϕ = zeros(M, 3, R)
    chain_log_σ_ϵ = zeros(1, R)
    chain_η = zeros(6, R)
    ϕ_star = zeros(size(startval_ϕ))
    accept_vec = zeros(3,R)
    accept = false
    loglik_vec = zeros(R)

    α_prob_accept_σ_ϵ = 0.
    α_prob_accept_ϕ = zeros(M)

    # first iteration
    chain_ϕ[:,:,1] = startval_ϕ
    chain_log_σ_ϵ[1] = startval_log_σ_ϵ
    print_interval = 1000

    if ρ > 0
        run_sort = true
    else
        run_sort = false
    end

    # set star values for η
    chain_η[:,1] = startval_η

    #=
    # sample star values for η
    for j = 1:3

        # get paramaters
        μ_0_j, M_0_j, α_j, β_j = prior_parameters_η[j,:]

        chain_η[j,1] = μ_0_j # set start value to mean
        chain_η[j+3,1] = (α_j - 1)/β_j # set start value to mode

    end
    =#

    accept_vec[:,1] = [M;1;3] # we accapt the start values

    # set random numbers
    #u_prop_old = randn(nbr_particles,N+1,M)
    #u_resample_old = randn(N,2,M)

    u_prop_old = [randn(nbr_particles[1], T_vec[1]+1)]
    u_resample_old = [randn(T_vec[1], 2)]

    for i in 2:M_subjects;
        append!(u_prop_old, [randn(nbr_particles[i],T_vec[i]+1)])
        append!(u_resample_old, [randn(T_vec[i],2)])
    end

    # estimate likelihood valie
    loglik_old_y = particel_filter(y, exp(startval_log_σ_ϵ), startval_ϕ, dt, u_prop_old, u_resample_old,  nbr_particles, run_sort)

    loglik_vec[1] = sum(loglik_old_y)

    #nbr_nan  = 0
    # main loop
    for r = 2:R

        # first gibbs stage, update \phi
        for m = 1:M
            ϕ_star[m,:] = chain_ϕ[m,:,r-1] + cholesky(exp(log_λ_i_ϕ[m])*Σ_i_ϕ[m]).U*randn(3)
        end

        u_prop_new = [randn(nbr_particles[1], T_vec[1]+1)]
        u_resample_new = [randn(T_vec[1], 2)]

        for i in 2:M_subjects;
            append!(u_prop_new, [randn(nbr_particles[i],T_vec[i]+1)])
            append!(u_resample_new, [randn(T_vec[i],2)])
        end

        u_prop_star = ρ*u_prop_old + sqrt(1-ρ^2)*u_prop_new # TODO check this!
        u_resample_star = ρ*u_resample_old + sqrt(1-ρ^2)*u_resample_new

        loglik_star_y = particel_filter(y, exp(chain_log_σ_ϵ[r-1]), ϕ_star, dt, u_prop_star, u_resample_star,  nbr_particles, run_sort)

        loglik_old_ϕ = calc_loglik_ϕ(chain_ϕ[:,:,r-1], chain_η[:,r-1])
        loglik_star_ϕ = calc_loglik_ϕ(ϕ_star, chain_η[:,r-1])

        for m = 1:M

            #jacobian_old =  sum(chain_ϕ[m,:,r-1])
            #jacobian_star = sum(ϕ_star[m,:])

            #log_α = (loglik_star_y[m] + loglik_star_ϕ[m] + jacobian_star) - (loglik_old_y[m] + loglik_old_ϕ[m] + jacobian_old)
            log_α = (loglik_star_y[m] + loglik_star_ϕ[m]) - (loglik_old_y[m] + loglik_old_ϕ[m])

            α_prob_accept_ϕ[m] = min(1, exp(log_α))

            # correct if we have  NaNs
            if loglik_star_y[m] == NaN
                log_α = log(0)
                α_prob_accept_ϕ[m] = min(1, exp(log_α))
                nbr_nan = nbr_nan +1
            end

            accept = log(rand()) < log_α

            if accept
                chain_ϕ[m,:,r] = ϕ_star[m,:]
                accept_vec[1,r] = accept_vec[1,r] +1
                loglik_old_y[m] = loglik_star_y[m]
                u_prop_old[m] = u_prop_star[m]
                u_resample_old[m] = u_resample_star[m]
            else
                chain_ϕ[m,:,r] = chain_ϕ[m,:,r-1]
            end

        end

        # secound gibbs stage, update log_σ_ϵ
        log_σ_ϵ_star = chain_log_σ_ϵ[r-1] + sqrt(exp(log_λ_i_σ_ϵ)*Σ_i_σ_ϵ)*randn()

        #jacobian_old =  chain_log_σ_ϵ[r-1]
        #jacobian_star = log_σ_ϵ_star

        loglik_star_y = particel_filter(y, exp(log_σ_ϵ_star), chain_ϕ[:,:,r], dt, u_prop_old, u_resample_old,  nbr_particles, run_sort)

        prior_old_log_σ_ϵ = normlogpdf(prior_parameters_σ_ϵ[1],1/prior_parameters_σ_ϵ[2],chain_log_σ_ϵ[r-1])
        prior_star_log_σ_ϵ = normlogpdf(prior_parameters_σ_ϵ[1],1/prior_parameters_σ_ϵ[2],log_σ_ϵ_star)

        #log_α = (sum(loglik_star_y) + prior_star_log_σ_ϵ + jacobian_star) - (sum(loglik_old_y) + prior_old_log_σ_ϵ + jacobian_old)
        log_α = (sum(loglik_star_y) + prior_star_log_σ_ϵ) - (sum(loglik_old_y) + prior_old_log_σ_ϵ)

        α_prob_accept_σ_ϵ = min(1, exp(log_α))

        # correct if we have  NaNs
        if any(isnan,loglik_star_y) == true
            #nbr_nan = nbr_nan +1
            #println(r)
            #println(log_σ_ϵ_star)
            #println(log_λ_i_σ_ϵ)
            #println(Σ_i_σ_ϵ)
            #println(chain_ϕ[:,:,r])
            #println(loglik_star_y)
            #println(particel_filter(y, exp(log_σ_ϵ_star), chain_ϕ[:,:,r], dt, u_prop_old, u_resample_old,  nbr_particles, run_sort))
            log_α = log(0)
            α_prob_accept_σ_ϵ = min(1, exp(log_α))
        end

        accept = log(rand()) < log_α

        if accept
            chain_log_σ_ϵ[r] = log_σ_ϵ_star
            accept_vec[2,r] = accept_vec[2,r] + 1
            loglik_old_y = loglik_star_y
        else
            chain_log_σ_ϵ[r] = chain_log_σ_ϵ[r-1]
        end

        # third stage, update η
        for j = 1:3

            # get paramaters
            μ_0_j, M_0_j, α_j, β_j = prior_parameters_η[j,:]

            ϕ_j = chain_ϕ[:,j,r-1]
            ϕ_j_mean = 1/M*sum(ϕ_j)

            α = α_j + div(M,2)
            β = β_j + 0.5*sum((ϕ_j .- ϕ_j_mean).^2) + (M*M_0_j)/(2*(M+M_0_j))*(ϕ_j_mean-μ_0_j)^2

            τ_j = rand_gamma(α) # TODO check scaling for Gamma
            τ_j = τ_j/β

            μ = (M*τ_j*ϕ_j_mean)/(M*τ_j + M_0_j*τ_j) + (M_0_j*τ_j*μ_0_j)/(M*τ_j + M_0_j*τ_j)
            τ = M*τ_j + M_0_j*τ_j

            μ_j = μ + sqrt(1/τ)*randn()

            chain_η[j,r] = μ_j
            chain_η[j+3,r] = τ_j

            accept_vec[3,r] = accept_vec[3,r] +1

        end

        loglik_vec[r] = sum(loglik_old_y)

        print_progress(r,print_interval,accept_vec,M)

        log_λ_i_σ_ϵ,μ_i_σ_ϵ,Σ_i_σ_ϵ = adaptive_update_σ_ϵ(r,update_interval,start_update,log_λ_i_σ_ϵ,μ_i_σ_ϵ,Σ_i_σ_ϵ,γ_σ_ϵ_0,α_power,chain_log_σ_ϵ,α_prob_accept_σ_ϵ)
        adaptive_update_ϕ!(r,update_interval,start_update,log_λ_i_ϕ,μ_i_ϕ,Σ_i_ϕ,γ_ϕ_0,α_power,M,chain_ϕ,α_prob_accept_ϕ)

    end

    #println(nbr_nan)
    println("Ending Gibbs sampler using CPMMH.")

    return chain_ϕ, chain_log_σ_ϵ, chain_η, accept_vec, loglik_vec

end

################################################################################
# help functions for Gibbs algorithms
################################################################################

# print progress to consol
function print_progress(r::Int,print_interval::Int,accept_vec::Array,M::Int)

    if mod(r-1,print_interval) == 0
        println("-------------")
        @printf "Percentage done:  %.2f %%\n" r/R*100
        @printf "Accaptance rate on %.0f to %.0f for ϕ: %.2f %%:\n" r-print_interval r sum(accept_vec[1,r-print_interval:r])/(M*print_interval)*100
        @printf "Accaptance rate on %.0f to %.0f for σ_ϵ: %.2f %%:\n" r-print_interval r sum(accept_vec[2,r-print_interval:r])/print_interval*100
    end

end

# Generlized AM adaptaton for the sigma_epsilon parameter,
# algorithm 4 in https://people.eecs.berkeley.edu/~jordan/sail/readings/andrieu-thoms.pdf)
function adaptive_update_σ_ϵ(r::Int,update_interval::Int,start_update::Int,
    log_λ_i_σ_ϵ::Real,μ_i_σ_ϵ::Real,Σ_i_σ_ϵ::Real,γ_σ_ϵ_0::Real,α_power::Real,
    chain_σ_ϵ::Array,α_prob_accept_σ_ϵ::Real)

    if mod(r-1,update_interval) == 0 && r-1 >= start_update

        γ = γ_σ_ϵ_0/r^α_power

        log_λ_i_σ_ϵ = log_λ_i_σ_ϵ + γ*(α_prob_accept_σ_ϵ - α_star_σ_ϵ)
        μ_i_σ_ϵ = μ_i_σ_ϵ + γ*(chain_σ_ϵ[r] - μ_i_σ_ϵ)
        Σ_i_σ_ϵ = Σ_i_σ_ϵ + γ*((chain_σ_ϵ[r]-μ_i_σ_ϵ)^2 - Σ_i_σ_ϵ)
        # perhaps we need the scalar version of Σ_i_ϕ[m] = (Σ_i_ϕ_tmp + Σ_i_ϕ_tmp')/2 here # add the transpose and diviade by 2 for numerical stability (to make sure that Σ_i_ϕ[m] is symmertic)

    end

    return log_λ_i_σ_ϵ,μ_i_σ_ϵ,Σ_i_σ_ϵ
end

# Generlized AM adaptaton for the random effects ϕ,
# algorithm 4 in https://people.eecs.berkeley.edu/~jordan/sail/readings/andrieu-thoms.pdf)
function adaptive_update_ϕ!(r::Int,update_interval::Int,start_update::Int,
    log_λ_i_ϕ::Vector,μ_i_ϕ::Array,Σ_i_ϕ::Array,γ_ϕ_0::Real,α_power::Real,
    M::Int,chain_ϕ::Array,α_prob_accept_ϕ::Vector)

    if mod(r-1,update_interval) == 0 && r-1 >= start_update

        γ = γ_ϕ_0/r^α_power

        for m in 1:M
            Σ_i_ϕ_tmp = Σ_i_ϕ[m]
            log_λ_i_ϕ[m] =  log_λ_i_ϕ[m]  + γ*(α_prob_accept_ϕ[m] - α_star_ϕ)
            μ_i_ϕ[m,:] = μ_i_ϕ[m,:] + γ*(chain_ϕ[m,:,r] - μ_i_ϕ[m,:])
            Σ_i_ϕ_tmp = Σ_i_ϕ_tmp + γ*((chain_ϕ[m,:,r]-μ_i_ϕ[m,:])*(chain_ϕ[m,:,r]-μ_i_ϕ[m,:])' - Σ_i_ϕ_tmp)
            Σ_i_ϕ[m] = (Σ_i_ϕ_tmp + Σ_i_ϕ_tmp')/2 # add the transpose and diviade by 2 for numerical stability (to make sure that Σ_i_ϕ[m] is symmertic)
        end

    end

end

function calc_loglik_ϕ(ϕ, η)

    loglik = zeros(size(ϕ,1))

    for m = 1:size(ϕ,1)
        for i in 1:size(ϕ,2)
            loglik[m] = loglik[m] + normlogpdf(η[i], sqrt(1/η[i+3]), ϕ[m,i])
        end
    end

    return loglik

end

# sample from Gamma(α,1)
# code adapted from https://www.cs.toronto.edu/~radford/csc2541.F04/gamma.html
function rand_gamma(α::Real)

    x = 0

    # TODO: We should have some check here such that alpha is a natrual number
    for i = 1:Int64(α)
        x = x + log(rand())
    end

    return -x

end

# parallel kalman/pf filter code (without diagnostics)

################################################################################
# Kalman/pf
################################################################################

"""
Parallel Kalman filter
"""
function kalman(y_full::Array, σ_ϵ::Real, ϕ::Array, τ_s::Real, diagnostics::Bool=false)

    M = length(y_full)

    #loglik_est = zeros(M)

    if diagnostics; path = Vector{Float64}[]; end #[zeros(length(100))]; end

    #loglik_est = SharedArray{Float64}(M)
    loglik_est = zeros(M)

    #@sync @distributed for m in 1:M

    Threads.@threads for m in 1:M

        if diagnostics
            loglik_est[m],x = kalman_filter(y_full[m].mV, σ_ϵ, ϕ[m,:], τ_s, diagnostics)
            append!(path,[x])
        else
            loglik_est[m] = kalman_filter(y_full[m].mV, σ_ϵ, ϕ[m,:], τ_s, diagnostics)
        end

    end

    #diagnostics == true ? (return Array(loglik_est), path) : (return Array(loglik_est))
    diagnostics == true ? (return loglik_est, path) : (return loglik_est)

end

# Kalman filter. Estimate the likelihood for the OU process (LTI) system. The code follows the notation in "Grey-box pharmacokinetic/pharmacodynamic modelling of euglycaemic clamp study".
function kalman_filter(y::Vector, σ_ϵ::Real, ϕ::Vector, τ_s::Real, diagnostics::Bool=false)

    T = length(y)

    if diagnostics; x_hat_kalman = zeros(T); end

    #println(loglik_est[m])
    λ = exp(ϕ[1])
    ν = exp(ϕ[2])
    σ = exp(ϕ[3]) # this is actually the correct parameterization for the Kalman filter

    B = ν
    A = -λ
    C = 1
    S = σ_ϵ^2

    # start values for the kalman filter
    P_start = var(y)
    x_hat_start = 0 # TODO: This is not correct we can fix this by sampling the start values form N(0, σ_ϵ)

    P_k = P_start
    x_k = x_hat_start

    loglik_est = 0.

    # main loop
    for k = 1:T

        x_k = exp(A*τ_s)*x_k + (1/A)*(exp(τ_s*A)-1)*B
        P_k = exp(A*τ_s)*P_k*exp(A*τ_s) + σ^2*(1/(2*A))*(exp(2*A*τ_s)-1)

        ϵ_k = y[k]-C*x_k
        R_k = C*P_k*C + S
        K = P_k*C*inv(R_k)

        #loglik_est = loglik_est - 0.5*(log(det(R_k)) + ϵ_k*inv(R_k)*ϵ_k)
        loglik_est = loglik_est + log_normalpdf(C*x_k, sqrt(R_k), y[k])

        x_k = x_k + K*ϵ_k
        P_k = P_k - K*R_k*K

        if diagnostics; x_hat_kalman[k] = x_k; end

    end

    diagnostics == true ? (return loglik_est, x_hat_kalman) : (return loglik_est)

end

"""
Parallel cpf
"""
function cpf(y_full::Array, σ_ϵ::Real, ϕ::Array, dt::Real, u_prop::Vector, u_resample::Vector, N_vec::Vector, run_sort::Bool)

    M = length(y_full)

    loglik_est = zeros(M)

    Threads.@threads for m in 1:M

        loglik_est[m] = particlefilter(y_full[m].mV, σ_ϵ, ϕ[m,:], dt, u_prop[m], u_resample[m], N_vec[m], run_sort, false)

    end

    return loglik_est

end

function cpf(y_full::Array, σ_ϵ::Real, ϕ::Array, dt::Real, u_prop::Vector, u_resample::Vector, N_vec::Vector, run_sort::Bool,diagnostics::Bool)

    M = length(y_full)

    if diagnostics;
        path = Matrix{Float64}[];
        nbr_resamlping = zeros(M)
    end

    loglik_est = zeros(M)

    for m in 1:M

        loglik_est[m],x,nbr_resamlping[m] = particlefilter(y_full[m].mV, σ_ϵ, ϕ[m,:], dt, u_prop[m], u_resample[m], N_vec[m], run_sort, diagnostics)
        push!(path, x)

    end

    diagnostics == true ? (return loglik_est, path, nbr_resamlping) : (return loglik_est)

end

# Estimates the loglikelihood using the correlated pf.
function particlefilter(y::Vector, σ_ϵ::Real, ϕ::Array, dt::Real, u_prop::Array, u_resample::Array, N::Real, run_sort::Bool, diagnostics::Bool=false)

    # pre-allocation

    loglik = 0.

    T = length(y)
    u_resample_calc = zeros(T)

    λ = exp(ϕ[1])
    ν = exp(ϕ[2])
    σ = exp(ϕ[3]) # this is actually the correct parameterization for the Kalman filter

    θ_1 = λ
    θ_2 = 1/λ*ν
    θ_3 = σ

    # convert standard normal to standard uniforms
    u_resample_temp = u_resample
    u_prop_temp = u_prop[:,2:end]
    u_prop_temp_init = u_prop[:,1]

    for i in 1:T
      u_resample_calc[i] = exp(-(u_resample_temp[i,1]^2+u_resample_temp[i,2]^2)/2)
    end

    # pre-allocation
    x = zeros(N,T) # particels
    w = zeros(N,T) # weigts

    nbr_resamlping = T

    # set start values
    xinit = zeros(N) + std(y)/10*u_prop_temp_init #(std(y[m,:]) + σ_ϵ)*randn(N)

    for t in 1:T

        if t == 1 # first iteration

          # propagate particels
          x[:,t] = stateprop(xinit, θ_1, θ_2, θ_3, dt, u_prop_temp[:,t]) #r*x_resample.*exp.(-x_resample .+ e[:,t])

        else

            # sort particles and weigths
            # alg=RadixSort
            if run_sort == true
                sorted_idx = sortperm(x[:,t-1]; alg = QuickSort) # ensure that we use QuickSort
                w[:,t-1] = w[sorted_idx,t-1]
                x[:,t-1] = x[sorted_idx,t-1]
            end

            # resample particels using systematic resampling
            ind = sysresample2(w[:,t-1],N,u_resample_calc[t])
            x_resample = x[ind,t-1]
            w[:,t-1] = 1/N*ones(N)

            # propagate particels
            x[:,t] = stateprop(x_resample, θ_1, θ_2, θ_3, dt, u_prop_temp[:,t]) #r*x_resample.*exp.(-x_resample .+ e[:,t])

        end

        # calc weigths and update loglik
        loglik = loglik + calc_weigths(w,x[:,t],y[t],σ_ϵ,t)

    end

    diagnostics == true ? (return loglik, x,nbr_resamlping) : (return loglik)

end

function bridge_pf(y_full::Array, σ_ϵ::Real, ϕ::Array, dt::Real, u_prop::Vector, u_resample::Vector, N_vec::Vector, run_sort::Bool)

    M = length(y_full)

    loglik_est = zeros(M)

    Threads.@threads for m in 1:M
        loglik_est[m] = diffusionbridgeparticlefilter(y_full[m].mV, σ_ϵ, ϕ[m,:], dt, u_prop[m], u_resample[m], N_vec[m], run_sort,false)
    end

    return loglik_est

end

function bridge_pf(y_full::Array, σ_ϵ::Real, ϕ::Array, dt::Real, u_prop::Vector, u_resample::Vector, N_vec::Vector, run_sort::Bool, diagnostics::Bool)

    M = length(y_full)

    if diagnostics;
        path = Matrix{Float64}[];
        nbr_resamlping = zeros(M)
    end

    loglik_est = zeros(M)

    for m in 1:M

        if diagnostics
            loglik_est[m],x,nbr_resamlping[m] = diffusionbridgeparticlefilter(y_full[m].mV, σ_ϵ, ϕ[m,:], dt, u_prop[m], u_resample[m], N_vec[m], run_sort,diagnostics)
            push!(path, x)
        else
            loglik_est[m] = diffusionbridgeparticlefilter(y_full[m].mV, σ_ϵ, ϕ[m,:], dt, u_prop[m], u_resample[m], N_vec[m], run_sort,diagnostics)
        end

    end

    diagnostics == true ? (return loglik_est, path, nbr_resamlping) : (return loglik_est)

end

# Estimates the loglikelihood using the correlated pf.
function diffusionbridgeparticlefilter(y::Vector, σ_ϵ::Real, ϕ::Array, dt::Real, u_prop::Array, u_resample::Array, N::Real, run_sort::Bool,diagnostics::Bool=false)

    # pre-allocation

    loglik = 0.

    T = length(y)
    u_resample_calc = zeros(T)

    λ = exp(ϕ[1])
    ν = exp(ϕ[2])
    σ = exp(ϕ[3]) # this is actually the correct parameterization for the Kalman filter

    # convert standard normal to standard uniforms
    u_resample_temp = u_resample
    u_prop_temp = u_prop[:,2:end]
    u_prop_temp_init = u_prop[:,1]

    for i in 1:T
      u_resample_calc[i] = exp(-(u_resample_temp[i,1]^2+u_resample_temp[i,2]^2)/2)
    end

    # pre-allocation
    x = zeros(N,T) # particels
    w = zeros(N,T) # weigts

    nbr_resamlping = T

    # set start values
    xinit = zeros(N) + std(y)/10*u_prop_temp_init #(std(y[m,:]) + σ_ϵ)*randn(N)

    for t in 1:T

        if t == 1 # first iteration

          #x[:,1] = xinit

          #loglik = loglik + calc_weigths(w,x[:,t],y[t],σ_ϵ,t)

          # propagate particels
          x_resample = xinit
          x[:,t], p_hat_mean, p_hat_std = sample_partilces(N, y[t], x_resample, u_prop_temp[:,t], dt, λ, ν, σ, σ_ϵ)

        else

            if run_sort == true
                sorted_idx = sortperm(x[:,t-1]; alg = QuickSort) # ensure that we use QuickSort
                w[:,t-1] = w[sorted_idx,t-1]
                x[:,t-1] = x[sorted_idx,t-1]
            end

            # resample particels using systematic resampling
            ind = sysresample2(w[:,t-1],N,u_resample_calc[t])
            x_resample = x[ind,t-1]
            w[:,t-1] = 1/N*ones(N)

            # propagate particels
            x[:,t], p_hat_mean, p_hat_std = sample_partilces(N, y[t], x_resample, u_prop_temp[:,t], dt, λ, ν, σ, σ_ϵ)

        end

        loglik = loglik + calc_apf_weigths(N, w, x[:,t], x_resample, p_hat_mean, p_hat_std, y[t],σ_ϵ,t, dt, λ, ν, σ)

    end

    diagnostics == true ? (return loglik, x,nbr_resamlping) : (return loglik)

end

################################################################################
## Testing of particel filters
################################################################################

# Kalman filter. Estimate the likelihood for the OU process (LTI) system. The code follows the notation in "Grey-box pharmacokinetic/pharmacodynamic modelling of euglycaemic clamp study".
function kalman_filter_test(y::Vector, σ_ϵ::Real, ϕ::Vector, τ_s::Real, diagnostics::Bool=false)

    T = length(y)

    if diagnostics; x_hat_kalman = zeros(T); end

    #println(loglik_est[m])
    λ = exp(ϕ[1])
    ν = exp(ϕ[2])
    σ = exp(ϕ[3]) # this is actually the correct parameterization for the Kalman filter

    B = ν
    A = -λ
    C = 1
    S = σ_ϵ^2

    # start values for the kalman filter
    P_start = 0 #var(y)
    x_hat_start = 5 # TODO: This is not correct we can fix this by sampling the start values form N(0, σ_ϵ)

    P_k = P_start
    x_k = x_hat_start

    loglik_est = 0.

    # main loop
    for k = 1:T

        x_k = exp(A*τ_s)*x_k + (1/A)*(exp(τ_s*A)-1)*B
        P_k = exp(A*τ_s)*P_k*exp(A*τ_s) + σ^2*(1/(2*A))*(exp(2*A*τ_s)-1)

        ϵ_k = y[k]-C*x_k
        R_k = C*P_k*C + S
        K = P_k*C*inv(R_k)

        #loglik_est = loglik_est - 0.5*(log(det(R_k)) + ϵ_k*inv(R_k)*ϵ_k)
        loglik_est = loglik_est + log_normalpdf(C*x_k, sqrt(R_k), y[k])

        x_k = x_k + K*ϵ_k
        P_k = P_k - K*R_k*K

        if diagnostics; x_hat_kalman[k] = x_k; end

    end

    diagnostics == true ? (return loglik_est, x_hat_kalman) : (return loglik_est)

end

# Estimates the loglikelihood using the correlated pf.
function particlefilter_test(y::Vector, σ_ϵ::Real, ϕ::Array, dt::Real, u_prop::Array, u_resample::Array, N::Real, run_sort::Bool, diagnostics::Bool=false)

    # pre-allocation

    loglik = 0.

    T = length(y)
    u_resample_calc = zeros(T)

    λ = exp(ϕ[1])
    ν = exp(ϕ[2])
    σ = exp(ϕ[3])

    θ_1 = λ
    θ_2 = 1/λ*ν
    θ_3 = σ

    # convert standard normal to standard uniforms
    u_resample_temp = u_resample
    u_prop_temp = u_prop[:,2:end]
    u_prop_temp_init = u_prop[:,1]

    for i in 1:T
      u_resample_calc[i] = exp(-(u_resample_temp[i,1]^2+u_resample_temp[i,2]^2)/2)
    end

    # pre-allocation
    x = zeros(N,T) # particels
    w = zeros(N,T) # weigts

    nbr_resamlping = T

    # set start values
    xinit = 5*ones(N)

    for t in 1:T

        if t == 1 # first iteration

            # propagate particels
            x[:,t] = stateprop(xinit, θ_1, θ_2, θ_3, dt, u_prop_temp[:,t]) #r*x_resample.*exp.(-x_resample .+ e[:,t])

        else

            # sort particles and weigths
            # alg=RadixSort
            if run_sort == true
                sorted_idx = sortperm(x[:,t-1]; alg = QuickSort) # ensure that we use QuickSort
                w[:,t-1] = w[sorted_idx,t-1]
                x[:,t-1] = x[sorted_idx,t-1]
            end

            # resample particels using systematic resampling
            ind = sysresample2(w[:,t-1],N,u_resample_calc[t])
            x_resample = x[ind,t-1]
            w[:,t-1] = 1/N*ones(N)

            # propagate particels
            x[:,t] = stateprop(x_resample, θ_1, θ_2, θ_3, dt, u_prop_temp[:,t]) #r*x_resample.*exp.(-x_resample .+ e[:,t])

        end

        # calc weigths and update loglik
        loglik = loglik + calc_weigths(w,x[:,t],y[t],σ_ϵ,t)

    end

    diagnostics == true ? (return loglik, x,nbr_resamlping) : (return loglik)

end

# The bridge filter
function diffusionbridgeparticlefilter_test(y::Vector, σ_ϵ::Real, ϕ::Array, dt::Real, u_prop::Array, u_resample::Array, N::Real, run_sort::Bool,diagnostics::Bool=false)

    # pre-allocation

    loglik = 0.

    T = length(y)
    u_resample_calc = zeros(T)

    λ = exp(ϕ[1])
    ν = exp(ϕ[2])
    σ = exp(ϕ[3]) # this is actually the correct parameterization for the Kalman filter

    # convert standard normal to standard uniforms
    u_resample_temp = u_resample
    u_prop_temp = u_prop[:,2:end]
    u_prop_temp_init = u_prop[:,1]

    for i in 1:T
      u_resample_calc[i] = exp(-(u_resample_temp[i,1]^2+u_resample_temp[i,2]^2)/2)
    end

    # pre-allocation
    x = zeros(N,T) # particels
    w = zeros(N,T) # weigts

    nbr_resamlping = T

    # set start values
    xinit = 5*ones(N)

    for t in 1:T

        if t == 1 # first iteration

          #x[:,1] = xinit

          #loglik = loglik + calc_weigths(w,x[:,t],y[t],σ_ϵ,t)

          # propagate particels
          x_resample = xinit
          x[:,t], p_hat_mean, p_hat_std = sample_partilces(N, y[t], x_resample, u_prop_temp[:,t], dt, λ, ν, σ, σ_ϵ)

        else

            if run_sort == true
                sorted_idx = sortperm(x[:,t-1]; alg = QuickSort) # ensure that we use QuickSort
                w[:,t-1] = w[sorted_idx,t-1]
                x[:,t-1] = x[sorted_idx,t-1]
            end

            # resample particels using systematic resampling
            ind = sysresample2(w[:,t-1],N,u_resample_calc[t])
            x_resample = x[ind,t-1]
            w[:,t-1] = 1/N*ones(N)

            # propagate particels
            x[:,t], p_hat_mean, p_hat_std = sample_partilces(N, y[t], x_resample, u_prop_temp[:,t], dt, λ, ν, σ, σ_ϵ)

        end

        loglik = loglik + calc_apf_weigths(N, w, x[:,t], x_resample, p_hat_mean, p_hat_std, y[t],σ_ϵ,t, dt, λ, ν, σ)

    end

    diagnostics == true ? (return loglik, x,nbr_resamlping) : (return loglik)

end

function calc_apf_weigths(N::Int, w::Array, x_new::Vector,x_old::Vector, p_hat_mean::Vector, p_hat_std::Real, y::Real, σ_ϵ::Real, t::Int, dt::Real, λ::Real, ν::Real, σ::Real)

    # compute weigths
    x_transtion_std = sqrt((σ^2/(2*λ))*(1-exp(-2*λ*dt)))

    logw = zeros(N)
    w_temp = zeros(N)

    # here we do not need the t > 1 condition since we t > 1 now
    @inbounds for i in 1:N
        x_transtion_mean = 1/λ*ν + (x_old[i] - 1/λ*ν)*exp(-λ*dt)
        logw[i] = log_normalpdf(x_new[i], σ_ϵ, y) + log_normalpdf(x_transtion_mean, x_transtion_std, x_new[i]) - log_normalpdf(p_hat_mean[i], p_hat_std, x_new[i])
    end

    # find largets wegith
    constant = maximum(logw)

    # subtract largets weigth
    @inbounds for i in 1:N; w_temp[i] = exp(logw[i] - constant); end

    # calc sum of weigths
    w_sum = sum(w_temp)

    # normalize weigths
    for i in 1:N; w_temp[i] = w_temp[i]/w_sum; end
    w[:,t] = w_temp

    return constant + log(w_sum) - log(N)

end

function sample_partilces(N::Int, y::Real, x::Vector, u_prop_temp::Vector, dt::Real, λ::Real, ν::Real, σ::Real, σ_ϵ::Real)

    x_new = zeros(N)
    p_hat_mean = zeros(N)
    α_0 = zeros(N)

    @inbounds for i in 1:N; α_0[i] = x[i]*exp(-λ*dt)+ν/λ*(1-exp(-λ*dt)) ; end
    β_0 = σ^2/(2*λ)*(1-exp(-2*λ*dt))

    @inbounds for i in 1:N; p_hat_mean[i] = α_0[i] + β_0*(β_0 + σ_ϵ^2)^(-1)*(y-α_0[i]); end
    p_hat_std = sqrt(β_0*(1-(β_0 + σ_ϵ^2)^(-1)*β_0))

    @inbounds for i in 1:N; x_new[i] =  p_hat_mean[i] + p_hat_std*u_prop_temp[i]; end

    return x_new, p_hat_mean, p_hat_std

end

################################################################################
# help functions for kalman/pf
################################################################################

# sample exactly from the OU process
function stateprop(x_old::Vector, θ_1::Real, θ_2::Real, θ_3::Real, dt::Real, u::Array)

    x_prop = zeros(length(x_old))

    σ_cond = sqrt((θ_3^2/(2*θ_1))*(1-exp(-2*θ_1*dt)))

    for i in 1:length(x_old)
        μ_cond = θ_2 + (x_old[i] - θ_2)*exp(-θ_1*dt)
        x_prop[i] = μ_cond + σ_cond*u[i]
    end

    return x_prop

end

# calc weigths and update loglik
function calc_weigths(w::Array, x::Array, y::Real, σ_ϵ::Real, t::Int)

    N = length(x) # nbr particels
    logw = zeros(N)
    w_temp = zeros(N)

    # calc w
    for i in 1:N; logw[i] = log_normalpdf(x[i], σ_ϵ, y); end

    # find largets wegith
    constant = maximum(logw)

    # subtract largets weigth
    for i in 1:N; w_temp[i] = exp(logw[i] - constant); end

    # calc sum of weigths
    w_sum = sum(w_temp)

    # update loglik
    #loglik =  loglik + constant + log(w_sum) - log(N)

    # normalize weigths
    for i in 1:N; w_temp[i] = w_temp[i]/w_sum; end

    w[:,t] = w_temp

    # return loglik
    return constant + log(w_sum) - log(N)

end

function log_normalpdf(μ::Real, σ::Real, x::Real)

  R = σ^2
  ϵ = x-μ

  return -0.5*(log(2*pi) + log(det(R)) + ϵ*inv(R)*ϵ)

  #return -0.5*(log(det(R)) + ϵ*inv(R)*ϵ)

end

#normlogpdf(1, 1, 1)
#log_normalpdf(1,1,1)

#=
function log_normalpdf(μ::Real, σ::Real, x::Real)

    # normlogpdf(μ::Real, σ::Real, x::Number)
    return normlogpdf(μ, σ, x)

end
=#

# Systematic resampling. Code adapted from Andrew Golightly.
function sysresample2(wts::Array,N::Int64,uni::Real)

  vec = zeros(Int64,N)
  wsum = sum(wts)
  k = 1
  u = uni/N
  wsumtarg = u
  wsumcurr=wts[k]/wsum
  delta = 1/N

  for i = 1:N
    while wsumcurr < wsumtarg
      k = k+1
      wsumcurr=wsumcurr+wts[k]/wsum
    end
    vec[i]=k
    wsumtarg=wsumtarg+delta
  end

  return vec

end

max_stars_repo_path: src/SDEMEM OU neuron data/mcmc.jl

max_issues_repo_path: src/SDEMEM OU neuron data/mcmc.jl

max_forks_repo_path: src/SDEMEM OU neuron data/mcmc.jl


text:
function brent(f::Function, x0::Number, x1::Number, args::Tuple=();
               xtol::AbstractFloat=1e-7, ytol=2eps(Float64),
               maxiter::Integer=50)
    EPS = eps(Float64)
    y0 = f(x0,args...)
    y1 = f(x1,args...)
    if abs(y0) < abs(y1)
        # Swap lower and upper bounds.
        x0, x1 = x1, x0
        y0, y1 = y1, y0
    end
    x2 = x0
    y2 = y0
    x3 = x2
    bisection = true
    for q in 1:maxiter
        # x-tolerance.
        if abs(x1-x0) < xtol
            return x1
        end

        # Use inverse quadratic interpolation if f(x0)!=f(x1)!=f(x2)
        # and linear interpolation (secant method) otherwise.
        if abs(y0-y2) > ytol && abs(y1-y2) > ytol
            x = x0*y1*y2/((y0-y1)*(y0-y2)) +
                x1*y0*y2/((y1-y0)*(y1-y2)) +
                x2*y0*y1/((y2-y0)*(y2-y1))
        else
            x = x1 - y1 * (x1-x0)/(y1-y0)
        end

        # Use bisection method if satisfies the conditions.
        delta = abs(2EPS*abs(x1))
        min1 = abs(x-x1)
        min2 = abs(x1-x2)
        min3 = abs(x2-x3)
        if (x < (3x0+x1)/4 && x > x1) ||
           (bisection && min1 >= min2/2) ||
           (!bisection && min1 >= min3/2) ||
           (bisection && min2 < delta) ||
           (!bisection && min3 < delta)
            x = (x0+x1)/2
            bisection = true
        else
            bisection = false
        end

        y = f(x,args...)
        # y-tolerance.
        if abs(y) < ytol
            return x
        end
        x3 = x2
        x2 = x1
        if sign(y0) != sign(y)
            x1 = x
            y1 = y
        else
            x0 = x
            y0 = y
        end
        if abs(y0) < abs(y1)
            # Swap lower and upper bounds.
            x0, x1 = x1, x0
            y0, y1 = y1, y0
        end
    end
    error("Max iteration exceeded")
end


text:
push!(LOAD_PATH, "/Users/antoinerestivo/Desktop/Permanents.jl-main/src")

using Documenter, Permanents

About = "Introduction" => "intro.md"

GettingStarted = "gettingstarted.md"

Functions = "Functions" => "index.md"

#
# Examples = "Examples" => [
#         "examples/flux.md"
#     ]
#
# License = "License" => "license.md"

PAGES = [About, GettingStarted, Functions]

# makedocs(
#     sitename = "Permanents.jl",
#     modules = [Permanents],
#     authors = "Benoit Seron, Antoine Restivo",
#     format = Documenter.HTML(),
# )

makedocs(
    sitename = "Permanents.jl",
    modules = [Permanents],
    authors = "Benoit Seron, Antoine Restivo",
    format = Documenter.HTML(),
    pages = PAGES
)

# deploydocs(repo = "github.com/Evizero/Augmentor.jl.git")
# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/AntoineRestivo/Permanents.jl.git"
)


text:
# test all solvers with the well known Woods test function
nlp = woods(n = 10)

solver = ALL_solvers[21]
(x, f, gNorm, iter, optimal, tired, status) = solver(nlp, verbose = true)

reset!(nlp);


text:
using OutlierDetectionData
using OutlierDetection: CLASS_NORMAL, CLASS_OUTLIER, normal_count, outlier_count
using Test

function test_tabular_dataset(dataset)
    names = dataset.list()
    always_accept = "DATADEPS_ALWAYS_ACCEPT" => true

    for name in names
        @info "Testing $name"
        X, y = withenv(() -> dataset.load(name), always_accept)
        @test ndims(y) == 1 # make sure that the labels are a vector
        @test ndims(X) == 2 # all tabular datasets are two-dimensional
        inliers, outliers = normal_count(y), outlier_count(y)

        # at least one outlier and inlier
        @test inliers > 0
        @test outliers > 0

        # no other elements than inliers and outliers
        @test inliers + outliers == length(y)
    end
end

@testset "Tabular datasets" begin
    test_tabular_dataset(ODDS)
    test_tabular_dataset(ELKI)
    test_tabular_dataset(TSAD)
end

@testset "Use prefix to list datasets" begin
    # empty prefix should return all datasets
    @test ODDS.list() == ODDS.list("")
    @test ELKI.list() == ELKI.list("")
    @test TSAD.list() == TSAD.list("")

    # prefix should return datasets that start with the prefix
    odds_starting_with_s = ["satellite", "satimage-2", "shuttle", "smtp", "speech"]
    odds_starting_with_sat = ["satellite", "satimage-2"]
    @test ODDS.list("s") == odds_starting_with_s
    @test ODDS.list("sat") == odds_starting_with_sat
    @test length(ELKI.list("Annthyroid_")) == 84 # there are 84 datasets starting with "Annthyroid_"
end


text:
abstract type Simulation{T,Dim} end

"""
    FrequencySimulation([particles::AbstractParticles=[],]
                        source::AbstractSource)

Build a FrequencySimulation. If particles are not provided, an empty array is used.

After building, you can [`run`](@ref) the simulation to get a [`FrequencySimulationResult`](@ref).
"""
mutable struct FrequencySimulation{T<:AbstractFloat,Dim,P<:PhysicalMedium} <: Simulation{T,Dim}
    "Vector of particles, can be of different types."
    particles::AbstractParticles
    "Source wave, where source.medium is the background medium of the simulation."
    source::AbstractSource{T,P}
end

# Constructor which infers parametric types from input arguments, note that we  don't need to do much type checking as the struct will error is inconsistent
function FrequencySimulation(particles::AbstractParticles{T,Dim}, source::AbstractSource{T,P}) where {Dim,T,P<:PhysicalMedium{T,Dim}}
    FrequencySimulation{T,Dim,P}(particles, source)
end

# A simulation with just sources is perfectly reasonable
function FrequencySimulation(source::AbstractSource{T,P}) where {Dim,T,P<:PhysicalMedium{T,Dim}}
    FrequencySimulation{T,Dim,P}(Vector{AbstractParticle{T,Dim}}(undef,0), source)
end

import Base.show
function show(io::IO, mime::MIME"text/plain", s::FrequencySimulation{T}) where {T}
    # text/plain is repl type
    # FrequencySimulation paramaters can be determined entirely from the medium and shape so we do not need to print them
    println(io, "FrequencySimulation{$T}")
    print(io,   "medium    = ")
    show(io, s.source.medium)
    println(io)
    println(io, "particles = ")
    for particle in s.particles
        show(io, particle)
    end
    return
end

import Base.run

run(particles::AbstractParticles, source::AbstractSource, x::Union{Shape,AbstractVector}, ω; kws...) = run(FrequencySimulation(particles,source), x, ω; kws...)

run(source::AbstractSource, x::Union{Shape,Vector}, ω; kws...) = run(FrequencySimulation(source), x, ω; kws...)

function run(sim::FrequencySimulation{T,Dim,P}, x::AbstractVector{T}, ωs::AbstractVector{T}=T[];
        kws...)::(SimulationResult{T,Dim,FieldDim} where FieldDim) where {Dim,P,T}
    run(sim,[x],ωs; kws...)
end

function run(sim::FrequencySimulation{T,Dim,P}, x::AbstractVector{T}, ω::T;
        kws...)::(SimulationResult{T,Dim,FieldDim} where FieldDim) where {Dim,P,T}
    run(sim,[x],[ω]; kws...)
end

# Main run function, all other run functions use this
function run(sim::FrequencySimulation{T,Dim,P}, x_vec::Union{Vector{Vector{T}},Vector{SVector{Dim,T}}}, ω::T;
        basis_order::Int = 5,
        only_scattered_waves::Bool = false, kws...) where {Dim,FieldDim,T,P<:PhysicalMedium{T,Dim,FieldDim}}

    x_vec = [SVector{Dim,T}(x...) for x in x_vec]

    # return just the source if there are no particles
    if isempty(sim.particles)
        f = field(sim.source)
        field_vec = f.(x_vec,ω)
    else
        # Calculate the outgoing basis coefficients around each particle, this is where most of the maths happens
        a_vec = basis_coefficients(sim, ω; basis_order=basis_order)

        if only_scattered_waves # remove source wave
            sim = FrequencySimulation(sim.particles,
                constant_source(sim.source.medium, zero(T)*im)
            )
        end

        # Evaluate the total field at the requested x positions
        field_vec = field(sim, ω, x_vec, a_vec)
    end

    # Construct results object
    field_vec = reshape(map(f->SVector{FieldDim,Complex{T}}(f), field_vec), :, 1)
    return FrequencySimulationResult{T,Dim,FieldDim}(field_vec, x_vec, Vector([ω]))

end

function run(sim::FrequencySimulation{T,Dim,P}, x_vec::Union{Vector{Vector{T}},Vector{SVector{Dim,T}}}, ωs::AbstractArray{T}=T[];
        ts::AbstractArray{T} = T[], result_in_time = !isempty(ts),
        basis_order::Int = 5,
        min_basis_order::Int = basis_order,
        basis_order_vec::AbstractVector{Int} = [-1],
        kws...)::(SimulationResult{T,Dim,FieldDim} where FieldDim)  where {Dim,P,T}

    if isempty(ωs) ωs = t_to_ω(ts) end
    x_vec = [SVector{Dim,T}(x...) for x in x_vec]

    # Considering basis_order to be the maximum basis order, then to avoid too much truncation error we use smaller basis orders on the smaller frequencies.
    if basis_order_vec == [-1]
        max_basis_order = max(basis_order,min_basis_order)
        basis_order_vec = Int.(round.(
            LinRange(min_basis_order, max_basis_order, length(ωs))
        ))
        basis_order_vec = basis_order_vec[sortperm(ωs)]
    end

    # if user asks for ω = 0, then we provide
    if first(ωs) == zero(T)
        # Compute for each angular frequency, then join up all the results
        fields = mapreduce(
            i -> run(sim,x_vec,ωs[i]; basis_order = basis_order_vec[i], kws...).field,
        hcat, eachindex(ωs)[2:end])

        # extrapolate field at ω = 0, which should be real when the time signal is real
        zeroresponse = real(ωs[3].*fields[:,1] - ωs[2].*fields[:,2]) ./ (ωs[3]-ωs[2])
        fields = reshape([zeroresponse; fields[:]], length(zeroresponse), size(fields,2)+1)
    else
        fields = mapreduce(
            i -> run(sim,x_vec,ωs[i]; basis_order = basis_order_vec[i], kws...).field,
        hcat, eachindex(ωs))
    end

    if !result_in_time
        FrequencySimulationResult(fields,x_vec,Vector(ωs))
    else
        @error "Using `result_in_time=true` in `run` is depricated and will be removed in later versions. Please instead use `frequency_to_time` on the results of `run`."
        # if isempty(ts) ts = ω_to_t(ωs) end
        #
        # # better to use the defaults of TimeSimulationResult's Constructor.
        # frequency_to_time(FrequencySimulationResult(fields,x_vec,Vector(ωs)); t_vec = reshape(ts,length(ts)), time_kws...)
    end
end

# Add docstring to run functions
"""
    run(sim::FrequencySimulation, x, ω; basis_order=5)

Run the simulation `sim` for the position `x` and angular frequency `ω`.

Position can be an SVector or Vector{SVector} and frequency can be a float or
vector of floats.
"""
function run(s::FrequencySimulation) throw(MethodError(run, (s,))) end

"""
    run(sim::FrequencySimulation, region::Shape;
        res=20, xres=res, yres=res, basis_order=5)

Run the simulation `sim` for a grid of positions in region and for angular frequency `ω`.

Frequency can be a float or vector of floats. The resolution of the grid points is defined
by xres and yres.
"""
function run(sim::FrequencySimulation{T,Dim}, region::Shape, ω_vec::AbstractVector;
            kws...) where {T,Dim}

    x_vec, inds = points_in_shape(region; kws...)
    # x_vec is a square grid of points and x_vec[inds] are the points in the region.

    result = run(sim, x_vec[inds], ω_vec; kws...)

    field_mat = zeros(typeof(result.field[1]), length(x_vec), length(ω_vec))
    field_mat[inds,:] = result.field

    return FrequencySimulationResult(field_mat, x_vec, ω_vec)
end

"""
    basis_coefficients(sim::FrequencySimulation, ω::AbstractFloat; basis_order::Int=5)::Matrix{Complex}

Return coefficients for bases around each particle for a given simulation and angular frequency (ω).
"""
function basis_coefficients(sim::FrequencySimulation{T,Dim,P}, ω::T; basis_order::Int = 5) where {Dim,P,T}

    # Precompute T-matrices for these particles
    t_matrices = get_t_matrices(sim.source.medium, sim.particles, ω, basis_order)

    # Compute scattering matrix for all particles
    S = scattering_matrix(sim.source.medium, sim.particles, t_matrices, ω, basis_order)

    # Get forcing vector from source, forms the right hand side of matrix equation to find basis_coefficients
    source_coefficient = regular_spherical_coefficients(sim.source)
    forcing = reduce(vcat, [source_coefficient(basis_order,origin(p),ω) for p in sim.particles])

    # Find scattering coefficients by solving this forcing
    a = (S + I) \ forcing

    # reshape and multiply by t-matrix to get the scattering coefficients
    a = reshape(a,basisorder_to_basislength(P,basis_order),length(sim.particles))
    for i in axes(a,2)
        a[:,i] = t_matrices[i] * a[:,i]
    end
    return a
end

function field(sim::FrequencySimulation{T,Dim,P}, ω::T, x_vec::Vector{v}, a_vec) where {Dim,P,T,v <: AbstractArray{T}}

    N = basislength_to_basisorder(P,size(a_vec,1))
    num_particles = length(sim.particles)
    basis = outgoing_basis_function(sim.source.medium, ω)

    function sum_basis(x)
        sum(eachindex(sim.particles)) do i
            p = sim.particles[i]
            sum(a_vec[:,i] .* basis(N, x-origin(p)))
        end
    end
    map(x_vec) do x
        j = findfirst(p -> x ∈ p, sim.particles)
        if typeof(j) === Nothing
            field(sim.source)(x,ω) + sum_basis(x)
        else
            p = sim.particles[j]
            internal_field(x, p, sim.source, ω, collect(a_vec[:,j]))
        end
    end
end