text:
#versioninfo()

#import Pkg; Pkg.add("Arpack")
#import Pkg; Pkg.add("ArgParse")

using LinearAlgebra
using Arpack
using SparseArrays
using ArgParse

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--Lx"
            help = "set Nsize (should be >=2)"
            arg_type = Int
            default = 4
        "--Ly"
            help = "set Nsize (should be >=2)"
            arg_type = Int
            default = 4
        "--J1"
            help = "set J1"
            arg_type = Float64
            default = 1.0
        "--J2"
            help = "set J2"
            arg_type = Float64
            default = 0.0
        "--Sz"
            help = "set Sz"
            arg_type = Int
            default = 0
    end
    return parse_args(s)
end

function snoob(x)
    next = 0
    if(x>0)
        smallest = x & -(x)
        ripple = x + smallest
        ones = xor(x,ripple)
        ones = div((ones >> 2),smallest)
        next = ripple | ones
    end
    return next
end

function count_bit(n)
    count = 0
    while (n>0)
        count += n & 1
        n >>= 1
    end
    return count
end

function init_parameters(N,Sz)
    Nup = div(N,2) + Sz
    Nhilbert = binomial(N,Nup)
    ihfbit = 1 << (div(N,2))
    irght = ihfbit-1
    ilft = xor(((1<<N)-1),irght)
    iup = (1<<(N-Nup))-1
    return Nup, Nhilbert, ihfbit, irght, ilft, iup
end

function make_list(N,Nup,Nhilbert,ihfbit,irght,ilft,iup)
    list_1 = zeros(Int,Nhilbert)
    list_ja = zeros(Int,ihfbit)
    list_jb = zeros(Int,ihfbit)
    ii = iup
    ja = 0
    jb = 0
    ia_old = ii & irght
    ib_old = div((ii & ilft),ihfbit)
    list_1[1] = ii
    list_ja[ia_old+1] = ja
    list_jb[ib_old+1] = jb
    ii = snoob(ii)
    for i in 1:Nhilbert-1
        ia = ii & irght
        ib = div((ii & ilft),ihfbit)
        if ib == ib_old
            ja += 1
        else
            jb += ja+1
            ja = 0
        end
        list_1[i+1] = ii
        list_ja[ia+1] = ja
        list_jb[ib+1] = jb
        ia_old = ia
        ib_old = ib
        ii = snoob(ii)
    end
    return list_1, list_ja, list_jb
end

function get_ja_plus_jb(ii,irght,ilft,ihfbit,list_ja,list_jb)
    ia = ii & irght
    ib = div((ii & ilft),ihfbit)
    ja = list_ja[ia+1]
    jb = list_jb[ib+1]
    return ja+jb
end

function make_hamiltonian(Jxx,Jzz,list_isite1,list_isite2,N,Nint,Nhilbert,irght,ilft,ihfbit,list_1,list_ja,list_jb)
    listki = ones(Int,(Nint+1)*Nhilbert)
    loc = ones(Int,(Nint+1)*Nhilbert)
    elemnt = zeros((Nint+1)*Nhilbert)
    for k in 0:Nint
        for i in 0:Nhilbert-1
            listki[k*Nhilbert+i+1] = i+1
        end
    end
    for k in 0:Nint-1
        isite1 = list_isite1[k+1]
        isite2 = list_isite2[k+1]
        is1 = 1<<isite1
        is2 = 1<<isite2
        is12 = is1 + is2
        wght = 2.0*Jxx[k+1]
        diag = Jzz[k+1]
        for i in 0:Nhilbert-1
            ii = list_1[i+1]
            ibit = ii & is12
            loc[Nint*Nhilbert+i+1] = i+1
            if (ibit==0 || ibit==is12)
                elemnt[Nint*Nhilbert+i+1] += diag
            else
                elemnt[Nint*Nhilbert+i+1] -= diag
                iexchg = xor(ii,is12)
                newcfg = get_ja_plus_jb(iexchg,irght,ilft,ihfbit,list_ja,list_jb)
                elemnt[k*Nhilbert+i+1] = wght
                loc[k*Nhilbert+i+1] = newcfg+1
            end
        end
    end
    HamCSC = dropzeros(sparse(listki,loc,elemnt))
    return HamCSC
end

function calc_zcorr(Nhilbert,Ncorr,list_corr_isite1,list_corr_isite2,psi,list_1)
    szz = zeros(Float64,Ncorr)
    for k in 0:Ncorr-1 # loop for all bonds for correlations
        isite1 = list_corr_isite1[k+1]
        isite2 = list_corr_isite2[k+1]
        is1 = 1<<isite1
        is2 = 1<<isite2
        is12 = is1 + is2
        corr = 0.0
        for i in 0:Nhilbert-1 # loop for all spin configurations with fixed Sz
            ii = list_1[i+1]
            ibit = ii & is12 # find sgmz.sgmz|uu> = |uu> or sgmz.sgmz|dd> = |dd>
            if (ibit==0 || ibit==is12) # if (spin1,spin2) = (00) or (11): factor = +1
                factor = +1.0
            else # if (spin1,spin2) = (01) or (10): factor = -1
                factor = -1.0
            end
            corr += factor*psi[i+1]^2 # psi[i]: real
        end
        szz[k+1] = 0.25 * corr
        if (isite1==isite2)
            szz[k+1] = 0.25
        end
    end
    return szz
end

function calc_xcorr(Nhilbert,Ncorr,list_corr_isite1,list_corr_isite2,psi,irght,ilft,ihfbit,list_1,list_ja,list_jb)
    sxx = zeros(Float64,Ncorr)
    for k in 0:Ncorr-1 # loop for all bonds for correlations
        isite1 = list_corr_isite1[k+1]
        isite2 = list_corr_isite2[k+1]
        is1 = 1<<isite1
        is2 = 1<<isite2
        is12 = is1 + is2
        corr = 0.0
        for i in 0:Nhilbert-1 # loop for all spin configurations with fixed Sz
            ii = list_1[i+1]
            ibit = ii & is12 # find sgmz.sgmz|ud> = -|ud> or sgmz.sgmz|du> = -|du>
            if (ibit==is1 || ibit==is2) # if (spin1,spin2) = (10) or (01)
                iexchg = xor(ii,is12) # find S+.S-|du> = |ud> or S-.S+|ud> = |du>
                newcfg = get_ja_plus_jb(iexchg,irght,ilft,ihfbit,list_ja,list_jb)
                corr += psi[i+1]*psi[newcfg+1] # psi[i]: real
            end
        end
        sxx[k+1] = 0.25 * corr
        if (isite1==isite2)
            sxx[k+1] = 0.25
        end
    end
    return sxx
end

function make_lattice(Lx,Ly,J1,J2)
#    Jxx = []
#    Jzz = []
#    list_isite1 = []
#    list_isite2 = []
    Jxx = Float64[]
    Jzz = Float64[]
    list_isite1 = Int[]
    list_isite2 = Int[]
    Nint = 0
    for iy in 0:Ly-1
        for ix in 0:Lx-1
            site1 = ix + Lx*iy
            site1x = mod((ix+1),Lx) + Lx*iy
            site1y = ix + Lx*(mod((iy+1),Ly))
            site1xpy = mod((ix+1),Lx) + Lx*(mod((iy+1),Ly))
            site1xmy = mod((ix+1),Lx) + Lx*(mod((iy-1+Ly),Ly))
#
            push!(list_isite1,site1)
            push!(list_isite2,site1x)
            push!(Jxx,J1)
            push!(Jzz,J1)
            Nint += 1
#
            push!(list_isite1,site1)
            push!(list_isite2,site1y)
            push!(Jxx,J1)
            push!(Jzz,J1)
            Nint += 1
#
            push!(list_isite1,site1)
            push!(list_isite2,site1xpy)
            push!(Jxx,J2)
            push!(Jzz,J2)
            Nint += 1
#
            push!(list_isite1,site1)
            push!(list_isite2,site1xmy)
            push!(Jxx,J2)
            push!(Jzz,J2)
            Nint += 1
        end
    end
    return Jxx, Jzz, list_isite1, list_isite2, Nint
end

function main()
    parsed_args = parse_commandline()
#    println("Parsed args:")
#    for (arg,val) in parsed_args
#        println("  $arg  =>  $val")
#    end

#    Lx = 4 # should be Lx>=2
#    Ly = 4 # should be Ly>=2
#    Sz = 0
#    J1 = 1.0
#    J2 = 0.0

    Lx = parsed_args["Lx"]
    Ly = parsed_args["Ly"]
    Sz = parsed_args["Sz"]
    J1 = parsed_args["J1"]
    J2 = parsed_args["J2"]
    N = Lx*Ly

    Nup, Nhilbert, ihfbit, irght, ilft, iup = init_parameters(N,Sz)

    binirght = string(irght,base=2,pad=N)
    binilft = string(ilft,base=2,pad=N)
    biniup = string(iup,base=2,pad=N)

    println("J1=",J1)
    println("J2=",J2)
    println("Lx=",Lx)
    println("Ly=",Ly)
    println("N=",N)
    println("Sz=",Sz)
    println("Nup=",Nup)
    println("Nhilbert=",Nhilbert)
    println("ihfbit=",ihfbit)
    println("irght,binirght=",irght," ",binirght)
    println("ilft,binilft=",ilft," ",binilft)
    println("iup,biniup=",iup," ",biniup)

    println("")
    @time list_1, list_ja, list_jb = make_list(N,Nup,Nhilbert,ihfbit,irght,ilft,iup)

    #println("list_1=",list_1)
    #println("list_ja=",list_ja)
    #println("list_jb=",list_jb)

    #println("")
    #println("i ii binii ja+jb")
    #for i in 0:Nhilbert-1
    #    ii = list_1[i+1]
    #    binii = string(ii,base=2,pad=N)
    #    ind = get_ja_plus_jb(ii,irght,ilft,ihfbit,list_ja,list_jb)
    #    println(i," ",ii," ",binii," ",ind)
    #end

    @time Jxx, Jzz, list_isite1, list_isite2, Nint = make_lattice(Lx,Ly,J1,J2)
    println(Jxx)
    println(Jzz)
    println(list_isite1)
    println(list_isite2)
    println("Nint=",Nint)

    println("")
    @time Ham = make_hamiltonian(Jxx,Jzz,list_isite1,list_isite2,N,Nint,Nhilbert,irght,ilft,ihfbit,list_1,list_ja,list_jb)
    #println(Ham)

    #@time ene,vec = eigs(Ham, nev=5, which=:LM)
    @time ene,vec = eigs(Ham, nev=5, which=:SR)
    ene = ene/N/4.0
    #println("# GS energy:",ene[1])
    println("# energy:",ene[1]," ",ene[2]," ",ene[3]," ",ene[4]," ",ene[5])
    #vec_sgn = sign(argmax(vec[:,1]))
    #println("# GS wave function:")
    #for i in 0:Nhilbert-1
    #    ii = list_1[i+1]
    #    binii = string(ii,base=2,pad=N)
    #    println(i," ",vec[i+1,1]*vec_sgn," ",binii)
    #end

    println("")
    Ncorr = N # number of total correlations
    list_corr_isite1 = [0 for k in 0:Ncorr-1] # site 1
    list_corr_isite2 = [k for k in 0:Ncorr-1] # site 2
    println(list_corr_isite1)
    println(list_corr_isite2)
    psi = vec[:,1] # choose the ground state
    @time szz = calc_zcorr(Nhilbert,Ncorr,list_corr_isite1,list_corr_isite2,psi,list_1)
    @time sxx = calc_xcorr(Nhilbert,Ncorr,list_corr_isite1,list_corr_isite2,psi,irght,ilft,ihfbit,list_1,list_ja,list_jb)
    ss = szz+sxx+sxx
    stot2 = N*sum(ss)
    println("# szz:",szz)
    println("# sxx:",sxx)
    println("# ss:",ss)
    println("# stot(stot+1):",stot2)

end

main()


text:
include("operator.jl")

function k_nearest_neighbor(coords1, coords2, k)
    """Compute k nearest neighbors between coords1 and coords2.

    Args:
        coords1: coordinates of centroid points (N, C)
        coords2: coordinates of all points (M, C)
        k: number of nearest neighbors

    Returns:
        idxs: indices of k nearest neighbors (N, k)
        square distances: square distance for kNN (N, k)
    """

    # batch proc.
    point_pairwise_distances = square_distance(coords1, coords2)
    # idxs = sortperm(point_pairwise_distances, dim=2)[:, 1:k]

    # square_dists = np.take_along_axis(point_pairwise_distances, idxs, axis=-1)

    return idxs
end

ext: jl


text:
using DataFrames, PyPlot, Statistics, FeatherFiles

AREA = 0.25445

dfw = DataFrame(load("./data/windyear.feather"))

plot(dfw.Power_Generation*AREA)
title("Wind power, 0.75 m² Wind Drone, Gouda, NL")


text:
# ---
# title: 82. Remove Duplicates from Sorted List II
# id: problem82
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Linked List
# link: <https://leetcode.com/problems/remove-duplicates-from-sorted-list-ii/description/>
# hidden: true
# ---
#
# Given a sorted linked list, delete all nodes that have duplicate numbers,
# leaving only _distinct_ numbers from the original list.
#
# Return the linked list sorted as well.
#
# **Example 1:**
#
#
#
#     Input: 1->2->3->3->4->4->5
#     Output: 1->2->5
#
#
# **Example 2:**
#
#
#
#     Input: 1->1->1->2->3
#     Output: 2->3
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module GridPortions

using Gridap
using Gridap.CellValuesGallery

export GridPortion

import Gridap: points
import Gridap: cells
import Gridap: celltypes
import Gridap: cellorders

struct GridPortion{D,Z} <: Grid{D,Z}
  oldgrid::Grid{D,Z}
  newcell_to_oldcell::IndexCellValue{Int}
  _newcells
  _newcelltypes
  _newcellorders
end

function GridPortion(oldgrid::Grid, newcell_to_oldcell)
  _newcell_to_oldcell = _setup_newcell_to_oldcell(newcell_to_oldcell)
  oldcells = cells(oldgrid)
  oldcelltypes = celltypes(oldgrid)
  oldcellorders = cellorders(oldgrid)
  _newcells = reindex(oldcells, _newcell_to_oldcell)
  _newcelltypes = reindex(oldcelltypes, _newcell_to_oldcell)
  _newcellorders = reindex(oldcellorders, _newcell_to_oldcell)
  GridPortion(
    oldgrid,_newcell_to_oldcell,_newcells,_newcelltypes,_newcellorders)
end

points(g::GridPortion) = points(g.oldgrid)

cells(g::GridPortion) = g._newcells

celltypes(g::GridPortion) = g._newcelltypes

cellorders(g::GridPortion) = g._newcellorders

_setup_newcell_to_oldcell(a::IndexCellValue) = a

_setup_newcell_to_oldcell(a::AbstractVector) = CellValueFromArray(a)

end # module


text:
@doc raw"""
    Coeffgamm2Simple(k, T, m)

The coefficient before the two point function flow.
"""
function Coeffgamm2Simple(k, T, m)
    (
        k * (
            -(coth(Epi(k, m) / (2 * T)) / Epi(k, m)^3) -
            csch(Epi(k, m) / (2 * T))^2 / (2 * T * Epi(k, m)^2)
        )
    ) / (16 * pi^2)
end

@doc raw"""
    VImSimple(p0, ps, q0, k, m, T, Npi, lam4pik)

compute $\mathrm{Im}V_k(q_0)$, the `k` dependence of $\lambda_{4\pi}$ and $m$ are neglected.

`VImSimple` only contains $V(q_0)$, no $V(-q_0)$

`VImSimple` contains type-1 delta function

`VImSimple` doesn't contains type-2 delta function
"""
function VImSimple(p0, ps, q0, k, m, T, Npi, lam4pik, UVScale)
    lam4pik^2 *
    (2 + Npi) *
    π *
    (
        3 * (
            F1All(p0 - q0, ps, k, m, T) - F1All(p0 - q0, ps, UVScale, m, T) +
            F1All(p0 + q0, ps, k, m, T) - F1All(p0 + q0, ps, UVScale, m, T)
            +
            F2All(p0 - q0, ps, k, m, T) - F2All(p0 - q0, ps, UVScale, m, T) +
            F2All(p0 + q0, ps, k, m, T) - F2All(p0 + q0, ps, UVScale, m, T)
        )
        # +
        # (2 + Npi) * (
        #     F1All(1e-8 - 1e-14, 1e-8, k, m, T) -
        #     F1All(1e-8 - 1e-14, 1e-8, UVScale, m, T) +
        #     F2All(1e-8 - 1e-14, 1e-8, k, m, T) -
        #     F2All(1e-8 - 1e-14, 1e-8, UVScale, m, T)
        # )
    )
end

@doc raw"""
    VImintqsSimple(p0, ps, k, T, Npi, m, lamda,UVScale)

compute $\int_0^kq_s^2dqs\int_{-1}^{1}d\cos\theta\mathrm{Im}V_k(q_0)$, the `k` dependence of $\lambda_{4\pi}$ and $m$ are neglected.

`VImintqsSimple` only contains $V(q_0)$, no $V(-q_0)$

`VImintqsSimple` contains type-1 delta function

`VImintqsSimple` doesn't contains type-2 delta function
"""
function VImintqsSimple(p0, ps, k, T, Npi, m, lamda,UVScale)
    hcubature(
        x ->
            x[1]^2 * VImSimple(
                p0,
                sqrt(x[1]^2 + ps^2 + 2 * x[1] * x[2] * ps),
                Epi(k, m),
                k,
                m,
                T,
                Npi,
                lamda,UVScale,
            ),
        [0.0, -1.0],
        [k, 1.0],atol=1e-6,rtol=1e-6,initdiv=10,maxevals=8000,
    )[1]
end

function VImintqsSimple(p0, ps,q0, k, T, Npi, m, lamda,UVScale)
    hcubature(
        x ->
            x[1]^2 * VImSimple(
                p0,
                sqrt(x[1]^2 + ps^2 + 2 * x[1] * x[2] * ps),
                q0,
                k,
                m,
                T,
                Npi,
                lamda,UVScale,
            ),
        [0.0, -1.0],
        [k, 1.0],atol=1e-2,rtol=1e-2,initdiv=40,maxevals=16000,
    )[1]
end

@doc raw"""
    propImSimple(p0, ps, T,IRScale,UVScale, Npi, m, lamda)
`x[1]` is `qs`, `x[2]` is `costh`, `x[3]` is `k`

# Arguments
- `m`: mass square, it's a constant number.
- `lamda`: $\lambda_{4\pi}$, it's a constant number.
"""
function propImSimple(p0, ps, T, IRScale, UVScale, Npi, m, lamda)
    -hquadrature(
        x ->
            2*Coeffgamm2Simple(x, T, m) *
            VImintqsSimple(p0, ps, x, T, Npi, m, lamda,UVScale),
        IRScale,
        UVScale,
        atol = 1e-4,
        rtol = 1e-4,maxevals=200,
    )[1]
end


text:
"""
Author: Dylan Morris
Date: 10 Feb 2022 at 6:15:09 pm
Description:
Use this script to quickly plot diagnostics of a chain. This could involve scatter-plots
to assess correlations between parameters, traceplots to check for divergent behaviours and
also just assessing the overall reasonableness of the posterior samples.
"""

let
    # for dataframe handling
    using DataFrames
    using CSV
    using Chain

    # nice plotting
    using StatsPlots
    using PlotThemes
    using Plots
    using LaTeXStrings

    # analysis/commonly used
    using Statistics
    using StatsBase
    using Dates
    using KernelDensity
    using Random
    using Distributions
    using BenchmarkTools
    using ProgressBars
end

gr()

default(
	linewidth=1.5,
	label=nothing,
	framestyle=:box,
)
scalefontsizes(0.95)

##

file_date = ARGS[1]

file_name = (
    "results/"
    * file_date
    * "/posterior_sample_"
    * file_date
    * ".csv"
)

samples = CSV.read(file_name, DataFrame)

dir_name = joinpath("figs", "posterior_summary", file_date)
if !ispath(dir_name)
    mkpath(dir_name)
end

let
    num_samples = nrow(samples)
    # we run 4 independent chains and so this just ensures that we get the right number of samples
    # per chain
    num_each_chain = num_samples ÷ 4

    sp = 1
    sp_total = 1
    page = 1

    names_to_plot = names(samples)
    names_to_plot_keep = ones(Bool, length(names_to_plot))

    for (i, name) in enumerate(names_to_plot)
        # this gross chunk of code enables us to avoid plotting of some of the nuisance parameters
        if length(name) >= length("mu_hat") && name[1:length("mu_hat")] == "mu_hat"
            names_to_plot_keep[i] = false
        elseif length(name) >= length("prop_md") && name[1:length("prop_md")] == "prop_md"
            names_to_plot_keep[i] = false
        elseif length(name) >= length("prop_masks") && name[1:length("prop_masks")] == "prop_masks"
            names_to_plot_keep[i] = false
        elseif length(name) >= length("brho") && name[1:length("brho")] == "brho"
            names_to_plot_keep[i] = false
        elseif length(name) >= length("ve") && name[1:length("ve")] == "ve"
            names_to_plot_keep[i] = false
        elseif length(name) >= length("md") && name[1:length("md")] == "md"
            names_to_plot_keep[i] = false
        elseif length(name) >= length("masks") && name[1:length("masks")] == "masks"
            names_to_plot_keep[i] = false
        elseif length(name) >= length("sus") && name[1:length("sus")] == "sus"
            names_to_plot_keep[i] = false
        elseif length(name) >= length("micro") && name[1:length("micro")] == "micro"
            names_to_plot_keep[i] = false
        elseif length(name) >= length("macro") && name[1:length("macro")] == "macro"
            names_to_plot_keep[i] = false
        elseif length(name) >= length("mob") && name[1:length("mob")] == "mob"
            names_to_plot_keep[i] = false
        end
    end

    names_to_plot = names_to_plot[names_to_plot_keep]

    fig = plot(
        layout = (5,5),
        dpi=200,
        size=(1000,1200),
        # link=:x,
        framestyle=:minimal,
        legend = false,
    )

    for name in ProgressBar(names_to_plot)

        if sp == 1
            fig = plot(
                layout = (5,5),
                dpi=200,
                size=(1000,1200),
                # link=:x,
                framestyle=:minimal,
                legend = false,
            )
        end
        # plot!(fig, subplot = sp, samples[!,name])
        for i in 1:4
            plot!(fig, subplot = sp, kde(samples[(1 + (i - 1) * num_each_chain):(i * num_each_chain),name]))
        end
        xlabel!(fig, subplot = sp, name)
        sp += 1
        sp_total += 1
        if sp == 5*5 + 1 || sp_total == length(names_to_plot) + 1
            sp = 1
            savefig(fig, dir_name*"/page"*string(page)*".pdf")
            page += 1
        end
    end

    using PDFmerger

    pdf_filenames = [
        dir_name*"/page"*string(p)*".pdf" for p in 1:page-1
    ]

    merge_pdfs(
        pdf_filenames,
        dir_name*"/kde.pdf",
        cleanup=true,
    )

    sp = 1
    sp_total = 1
    page = 1

    fig = plot(
        layout = (5,5),
        dpi=200,
        size=(1000,1200),
        # link=:x,
        framestyle=:minimal,
        legend = false,
    )

    for name in ProgressBar(names_to_plot)
        if sp == 1
            fig = plot(
                layout = (5,5),
                dpi=200,
                size=(1000,1200),
                # link=:x,
                framestyle=:minimal,
                legend = false,
            )
        end
        # plot!(fig, subplot = sp, samples[!,name])
        for i in 1:4
            plot!(fig, subplot = sp, samples[(1 + (i - 1) * num_each_chain):(i * num_each_chain),name])
        end
        xlabel!(fig, subplot = sp, name)
        sp += 1
        sp_total += 1
        if sp == 5*5 + 1 || sp_total == length(names_to_plot) + 1
            sp = 1
            savefig(fig, dir_name*"/page"*string(page)*".pdf")
            page += 1
        end
    end

    using PDFmerger

    pdf_filenames = [
        dir_name*"/page"*string(p)*".pdf" for p in 1:page-1
    ]

    merge_pdfs(
        pdf_filenames,
        dir_name*"/traceplots.pdf",
        cleanup=true,
    )

end


text:
using Revise
using Test

@test isempty(detect_ambiguities(Revise, Base, Core))

using Pkg, Unicode, Distributed, InteractiveUtils, REPL
using OrderedCollections: OrderedSet

include("common.jl")

to_remove = String[]

throwing_function(bt) = bt[2]

function rm_precompile(pkgname::AbstractString)
    filepath = Base.cache_file_entry(Base.PkgId(pkgname))
    for depot in DEPOT_PATH
        fullpath = joinpath(depot, filepath)
        isfile(fullpath) && rm(fullpath)
    end
end

# A junk module that we can evaluate into
module ReviseTestPrivate
struct Inner
    x::Float64
end
end

function private_module()
    modname = gensym()
    Core.eval(ReviseTestPrivate, :(module $modname end))
end

function compare_sigs(ex, typex)
    mod = private_module()
    f = Core.eval(mod, ex)
    mths = methods(f)
    for (m, tex) in zip(mths, typex)
        t = Core.eval(mod, tex)
        @test m.sig == t
    end
end

function compare_sigs(ex)
    sig = Revise.get_signature(ex)
    typex = Revise.sig_type_exprs(sig)
    compare_sigs(ex, typex)
end

@testset "Revise" begin

    function collectexprs(ex::Revise.RelocatableExpr)
        exs = Revise.RelocatableExpr[]
        for item in Revise.LineSkippingIterator(ex.args)
            push!(exs, item)
        end
        exs
    end

    function get_docstring(ds)
        docstr = ds.content[1]
        while !isa(docstr, AbstractString)
            docstr = docstr.content[1]
        end
        return docstr
    end

    @testset "LineSkipping" begin
        ex = Revise.relocatable!(quote
                                 f(x) = x^2
                                 g(x) = sin(x)
                                 end)
        @test length(ex.args) == 4  # including the line number expressions
        exs = collectexprs(ex)
        @test length(exs) == 2
        @test isequal(exs[1], Revise.relocatable!(:(f(x) = x^2)))
        @test !isequal(exs[2], Revise.relocatable!(:(f(x) = x^2)))
        @test isequal(exs[2], Revise.relocatable!(:(g(x) = sin(x))))
        @test !isequal(exs[1], Revise.relocatable!(:(g(x) = sin(x))))
        @test string(ex) == """
quote
    f(x) = begin
            x ^ 2
        end
    g(x) = begin
            sin(x)
        end
end"""
    end

    @testset "Parse errors" begin
        warnfile = randtmp()
        open(warnfile, "w") do io
            redirect_stderr(io) do
                md = Revise.FileModules(Main)
                @test Revise.parse_source!(md, """
begin # this block should parse correctly, cf. issue #109

end
f(x) = 1
g(x) = 2
h{x) = 3  # error
k(x) = 4
""",
                                            :test, 1, Main) == nothing
                @test haskey(md[Main].defmap, convert(Revise.RelocatableExpr, :(g(x) = 2)))
            end
        end
        @test occursin("parsing error near line 6", read(warnfile, String))
        rm(warnfile)
        @test Revise.is_linenumber(LineNumberNode(5, "foo.jl"))  # issue #100
    end

    @testset "Methods and signatures" begin
        compare_sigs(:(foo(x) = 1))
        compare_sigs(:(foo(x::Int) = 1))
        # where signatures
        compare_sigs(:(foo(x::T) where T<:Integer = T))
        compare_sigs(:(foo(x::V) where V<:Array{T} where T = T))
        # Varargs
        compare_sigs(:(foo(x::Int, y...) = 1))
        compare_sigs(:(foo(x::Int, y::Symbol...) = 1))
        compare_sigs(:(foo(x::T, y::U...) where {T<:Integer,U} = U))
        compare_sigs(:(foo(x::Array{Float64,K}, y::Vararg{Symbol,K}) where K = K))
        # Default args
        compare_sigs(:(foo(x, y=0) = 1))
        compare_sigs(:(foo(x, y::Int=0) = 1))
        compare_sigs(:(foo(x, y="hello", z::Int=0) = 1))
        compare_sigs(:(foo(x::Array{Float64,K}, y::Int=0) where K = K))
        # Keyword args
        compare_sigs(:(foo(x; y=0) = 1))
        compare_sigs(:(foo(x; y::Int=0) = 1))
        compare_sigs(:(foo(x; y="hello", z::Int=0) = 1))
        compare_sigs(:(foo(x::Array{Float64,K}; y::Int=0) where K = K))
        # Default and keyword args
        compare_sigs(:(foo(x, y="hello"; z::Int=0) = 1))
        # Destructured args
        compare_sigs(:(foo(x, (count, name)) = 1))

        # Do it all again for long-form declarations
        compare_sigs(:(function foo(x) return 1 end))
        compare_sigs(:(function foo(x::Int) return 1 end))
        # where signatures
        compare_sigs(:(function foo(x::T) where T<:Integer return T end))
        compare_sigs(:(function foo(x::V) where V<:Array{T} where T return T end))
        # Varargs
        compare_sigs(:(function foo(x::Int, y...) return 1 end))
        compare_sigs(:(function foo(x::Int, y::Symbol...) return 1 end))
        compare_sigs(:(function foo(x::T, y::U...) where {T<:Integer,U} return U end))
        compare_sigs(:(function foo(x::Array{Float64,K}, y::Vararg{Symbol,K}) where K return K end))
        # Default args
        compare_sigs(:(function foo(x, y=0) return 1 end))
        compare_sigs(:(function foo(x, y::Int=0) return 1 end))
        compare_sigs(:(function foo(x, y="hello", z::Int=0) return 1 end))
        compare_sigs(:(function foo(x::Array{Float64,K}, y::Int=0) where K return K end))
        # Keyword args
        compare_sigs(:(function foo(x; y=0) return 1 end))
        compare_sigs(:(function foo(x; y::Int=0) return 1 end))
        compare_sigs(:(function foo(x; y="hello", z::Int=0) return 1 end))
        compare_sigs(:(function foo(x::Array{Float64,K}; y::Int=0) where K return K end))
        # Default and keyword args
        compare_sigs(:(function foo(x, y="hello"; z::Int=0) return 1 end))
        # Destructured args
        compare_sigs(:(function foo(x, (count, name)) return 1 end))

        # Return type annotations
        @test Revise.sig_type_exprs(:(typeinfo_eltype(typeinfo::Type)::Union{Type,Nothing})) ==
              Revise.sig_type_exprs(:(typeinfo_eltype(typeinfo::Type)))

        # Overloading call
        def = :((i::Inner)(::String) = i.x)
        sig = Revise.get_signature(def)
        sigexs = Revise.sig_type_exprs(sig)
        Core.eval(ReviseTestPrivate, def)
        i = ReviseTestPrivate.Inner(3)
        m = @which i("hello")
        @test Core.eval(ReviseTestPrivate, sigexs[1]) == m.sig
        def = :((::Type{Inner})(::Dict) = 17)
        sig = Revise.get_signature(def)
        sigexs = Revise.sig_type_exprs(sig)
        Core.eval(ReviseTestPrivate, def)
        m = @which ReviseTestPrivate.Inner(Dict("a"=>1))
        @test Core.eval(ReviseTestPrivate, sigexs[1]) == m.sig

        # Annotations
        refex =  Revise.relocatable!(:(function foo(x) x^2 end))
        for ex in (:(@inline function foo(x) x^2 end),
                   :(@noinline function foo(x) x^2 end),
                   :(@propagate_inbounds function foo(x) x^2 end))
            @test Revise.get_signature(ex) == :(foo(x))
            @test Revise.relocatable!(Revise.funcdef_expr(ex)) == refex
        end

        # @eval-defined methods
        ex = :(@eval getindex(A::Array, i1::Int, i2::Int, I::Int...) = (@_inline_meta; arrayref($(Expr(:boundscheck)), A, i1, i2, I...)))
        @test Revise.get_signature(ex) == :(getindex(A::Array, i1::Int, i2::Int, I::Int...))
        @test Revise.relocatable!(Revise.funcdef_expr(ex)) == Revise.relocatable!(
            :(getindex(A::Array, i1::Int, i2::Int, I::Int...) = (@_inline_meta; arrayref($(Expr(:boundscheck)), A, i1, i2, I...)))
            )
    end

    @testset "Comparison and line numbering" begin
        fl1 = joinpath(@__DIR__, "revisetest.jl")
        fl2 = joinpath(@__DIR__, "revisetest_revised.jl")
        fl3 = joinpath(@__DIR__, "revisetest_errors.jl")
        include(fl1)  # So the modules are defined
        # test the "mistakes"
        @test ReviseTest.cube(2) == 16
        @test ReviseTest.Internal.mult3(2) == 8
        oldmd = Revise.parse_source(fl1, Main)
        Revise.instantiate_sigs!(oldmd)
        newmd = Revise.parse_source(fl2, Main)
        revmd = Revise.eval_revised(newmd, oldmd)
        @test ReviseTest.cube(2) == 8
        @test ReviseTest.Internal.mult3(2) == 6

        @test length(revmd) == 3
        @test haskey(revmd, ReviseTest) && haskey(revmd, ReviseTest.Internal)

        dvs = collect(revmd[ReviseTest].defmap)
        @test length(dvs) == 3
        (def, val) = dvs[1]
        @test isequal(def,  Revise.relocatable!(:(square(x) = x^2)))
        @test val == (DataType[Tuple{typeof(ReviseTest.square),Any}], 0)
        @test Revise.firstlineno(def) == 5
        m = @which ReviseTest.square(1)
        @test m.line == 5
        @test revmd[ReviseTest].sigtmap[Tuple{typeof(ReviseTest.square),Any}] == def
        (def, val) = dvs[2]
        @test isequal(def, Revise.relocatable!(:(cube(x) = x^3)))
        @test val == (DataType[Tuple{typeof(ReviseTest.cube),Any}], 0)
        m = @which ReviseTest.cube(1)
        @test m.line == 7
        @test revmd[ReviseTest].sigtmap[Tuple{typeof(ReviseTest.cube),Any}] == def
        (def, val) = dvs[3]
        @test isequal(def, Revise.relocatable!(:(fourth(x) = x^4)))
        @test val == (DataType[Tuple{typeof(ReviseTest.fourth),Any}], 0)
        m = @which ReviseTest.fourth(1)
        @test m.line == 9
        @test revmd[ReviseTest].sigtmap[Tuple{typeof(ReviseTest.fourth),Any}] == def

        dvs = collect(revmd[ReviseTest.Internal].defmap)
        @test length(dvs) == 2
        (def, val) = dvs[1]
        @test isequal(def,  Revise.relocatable!(:(mult2(x) = 2*x)))
        @test val == (DataType[Tuple{typeof(ReviseTest.Internal.mult2),Any}], 2)  # 2 because it shifted down 2 lines and was not evaled
        @test Revise.firstlineno(def) == 11
        m = @which ReviseTest.Internal.mult2(1)
        @test m.line == 11
        @test revmd[ReviseTest.Internal].sigtmap[Tuple{typeof(ReviseTest.Internal.mult2),Any}] == def
        (def, val) = dvs[2]
        @test isequal(def, Revise.relocatable!(:(mult3(x) = 3*x)))
        @test val == (DataType[Tuple{typeof(ReviseTest.Internal.mult3),Any}], 0)  # 0 because it was freshly-evaled
        m = @which ReviseTest.Internal.mult3(1)
        @test m.line == 14
        @test revmd[ReviseTest.Internal].sigtmap[Tuple{typeof(ReviseTest.Internal.mult3),Any}] == def

        # Backtraces
        newmd = Revise.parse_source(fl3, Main)
        revmd = Revise.eval_revised(newmd, revmd)
        try
            ReviseTest.cube(2)
            @test false
        catch err
            @test isa(err, ErrorException) && err.msg == "cube"
            bt = throwing_function(stacktrace(catch_backtrace()))
            @test bt.func == :cube && bt.file == Symbol(fl3) && bt.line == 7
        end
        try
            ReviseTest.Internal.mult2(2)
            @test false
        catch err
            @test isa(err, ErrorException) && err.msg == "mult2"
            bt = throwing_function(stacktrace(catch_backtrace()))
            @test bt.func == :mult2 && bt.file == Symbol(fl3) && bt.line == 13
        end
    end

    @testset "Display" begin
        io = IOBuffer()
        show(io, Revise.relocatable!(:(@inbounds x[2])))
        str = String(take!(io))
        @test str == ":(@inbounds x[2])"
        fm = Revise.parse_source(joinpath(@__DIR__, "revisetest.jl"), Main)
        Revise.instantiate_sigs!(fm)
        @test string(fm) == "OrderedCollections.OrderedDict(Main=>FMMaps(<1 expressions>, <0 signatures>),Main.ReviseTest=>FMMaps(<2 expressions>, <2 signatures>),Main.ReviseTest.Internal=>FMMaps(<2 expressions>, <2 signatures>))"
        fmmr = fm[ReviseTest]
        @test string(fmmr) == "FMMaps(<2 expressions>, <2 signatures>)"
        io = IOBuffer()
        print(IOContext(io, :limit=>false), fmmr)
        str = String(take!(io))
        @test str == "FMMaps with the following expressions:\n  :(square(x) = begin\n          x ^ 2\n      end)\n  :(cube(x) = begin\n          x ^ 4\n      end)\n"
    end

    @testset "File paths" begin
        testdir = randtmp()
        mkdir(testdir)
        push!(to_remove, testdir)
        push!(LOAD_PATH, testdir)
        for (pcflag, fbase) in ((true, "pc"), (false, "npc"))  # precompiled & not
            modname = uppercase(fbase)
            pcexpr = pcflag ? :() : :(__precompile__(false))
            # Create a package with the following structure:
            #   src/PkgName.jl   # PC.jl = precompiled, NPC.jl = nonprecompiled
            #   src/file2.jl
            #   src/subdir/file3.jl
            #   src/subdir/file4.jl
            # exploring different ways of expressing the `include` statement
            dn = joinpath(testdir, modname, "src")
            mkpath(dn)
            open(joinpath(dn, modname*".jl"), "w") do io
                println(io, """
$pcexpr
module $modname

export $(fbase)1, $(fbase)2, $(fbase)3, $(fbase)4, $(fbase)5, using_macro_$(fbase)

$(fbase)1() = 1

include("file2.jl")
include("subdir/file3.jl")
include(joinpath(@__DIR__, "subdir", "file4.jl"))
otherfile = "file5.jl"
include(otherfile)

# Update order check: modifying `some_macro_` to return -6 doesn't change the
# return value of `using_macro_` (issue #20) unless `using_macro_` is also updated,
# *in this order*:
#   1. update the `@some_macro_` definition
#   2. update the `using_macro_` definition
macro some_macro_$(fbase)()
    return 6
end
using_macro_$(fbase)() = @some_macro_$(fbase)()

end
""")
            end
            open(joinpath(dn, "file2.jl"), "w") do io
                println(io, "$(fbase)2() = 2")
            end
            mkdir(joinpath(dn, "subdir"))
            open(joinpath(dn, "subdir", "file3.jl"), "w") do io
                println(io, "$(fbase)3() = 3")
            end
            open(joinpath(dn, "subdir", "file4.jl"), "w") do io
                println(io, "$(fbase)4() = 4")
            end
            open(joinpath(dn, "file5.jl"), "w") do io
                println(io, "$(fbase)5() = 5")
            end
            sleep(2.1)   # so the defining files are old enough not to trigger mtime criterion
            @eval using $(Symbol(modname))
            fn1, fn2 = Symbol("$(fbase)1"), Symbol("$(fbase)2")
            fn3, fn4 = Symbol("$(fbase)3"), Symbol("$(fbase)4")
            fn5 = Symbol("$(fbase)5")
            fn6 = Symbol("using_macro_$(fbase)")
            @eval @test $(fn1)() == 1
            @eval @test $(fn2)() == 2
            @eval @test $(fn3)() == 3
            @eval @test $(fn4)() == 4
            @eval @test $(fn5)() == 5
            @eval @test $(fn6)() == 6
            m = @eval first(methods($fn1))
            @test Revise.get_def(m) == convert(Revise.RelocatableExpr, :( $fn1() = 1 ))
            sleep(0.1)  # to ensure that the file watching has kicked in
            # Change the definition of function 1 (easiest to just rewrite the whole file)
            open(joinpath(dn, modname*".jl"), "w") do io
                println(io, """
$pcexpr
module $modname
export $(fbase)1, $(fbase)2, $(fbase)3, $(fbase)4, $(fbase)5, using_mac$(fbase)
$(fbase)1() = -1
include("file2.jl")
include("subdir/file3.jl")
include(joinpath(@__DIR__, "subdir", "file4.jl"))
otherfile = "file5.jl"
include(otherfile)

macro some_macro_$(fbase)()
    return -6
end
using_macro_$(fbase)() = @some_macro_$(fbase)()

end
""")  # just for fun we skipped the whitespace
            end
            yry()
            @eval @test $(fn1)() == -1
            @eval @test $(fn2)() == 2
            @eval @test $(fn3)() == 3
            @eval @test $(fn4)() == 4
            @eval @test $(fn5)() == 5
            @eval @test $(fn6)() == 6      # because it hasn't been re-macroexpanded
            @test revise(eval(Symbol(modname)))
            @eval @test $(fn6)() == -6
            # Redefine function 2
            open(joinpath(dn, "file2.jl"), "w") do io
                println(io, "$(fbase)2() = -2")
            end
            yry()
            @eval @test $(fn1)() == -1
            @eval @test $(fn2)() == -2
            @eval @test $(fn3)() == 3
            @eval @test $(fn4)() == 4
            @eval @test $(fn5)() == 5
            @eval @test $(fn6)() == -6
            open(joinpath(dn, "subdir", "file3.jl"), "w") do io
                println(io, "$(fbase)3() = -3")
            end
            yry()
            @eval @test $(fn1)() == -1
            @eval @test $(fn2)() == -2
            @eval @test $(fn3)() == -3
            @eval @test $(fn4)() == 4
            @eval @test $(fn5)() == 5
            @eval @test $(fn6)() == -6
            open(joinpath(dn, "subdir", "file4.jl"), "w") do io
                println(io, "$(fbase)4() = -4")
            end
            yry()
            @eval @test $(fn1)() == -1
            @eval @test $(fn2)() == -2
            @eval @test $(fn3)() == -3
            @eval @test $(fn4)() == -4
            @eval @test $(fn5)() == 5
            @eval @test $(fn6)() == -6
            open(joinpath(dn, "file5.jl"), "w") do io
                println(io, "$(fbase)5() = -5")
            end
            yry()
            @eval @test $(fn1)() == -1
            @eval @test $(fn2)() == -2
            @eval @test $(fn3)() == -3
            @eval @test $(fn4)() == -4
            @eval @test $(fn5)() == -5
            @eval @test $(fn6)() == -6
            # Check module2files
            files = [joinpath(dn, modname*".jl"), joinpath(dn, "file2.jl"),
                     joinpath(dn, "subdir", "file3.jl"),
                     joinpath(dn, "subdir", "file4.jl"),
                     joinpath(dn, "file5.jl")]
            @test Revise.module2files[Symbol(modname)] == files
        end
        # Remove the precompiled file
        rm_precompile("PC")

        # Test files paths that can't be statically parsed
        dn = joinpath(testdir, "LoopInclude", "src")
        mkpath(dn)
        open(joinpath(dn, "LoopInclude.jl"), "w") do io
            println(io, """
module LoopInclude

export li_f, li_g

for fn in ("file1.jl", "file2.jl")
    include(fn)
end

end
""")
        end
        open(joinpath(dn, "file1.jl"), "w") do io
            println(io, "li_f() = 1")
        end
        open(joinpath(dn, "file2.jl"), "w") do io
            println(io, "li_g() = 2")
        end
        sleep(2.1) # so the defining files are old enough not to trigger mtime criterion
        @eval using LoopInclude
        sleep(0.1) # to ensure file-watching is set up
        @test li_f() == 1
        @test li_g() == 2
        sleep(1.1)  # ensure watching is set up
        yry()
        open(joinpath(dn, "file1.jl"), "w") do io
            println(io, "li_f() = -1")
        end
        yry()
        @test li_f() == -1

        pop!(LOAD_PATH)
    end

    # issue #36
    @testset "@__FILE__" begin
        testdir = randtmp()
        mkdir(testdir)
        push!(to_remove, testdir)
        push!(LOAD_PATH, testdir)
        dn = joinpath(testdir, "ModFILE", "src")
        mkpath(dn)
        open(joinpath(dn, "ModFILE.jl"), "w") do io
            println(io, """
module ModFILE

mf() = @__FILE__, 1

end
""")
        end
        sleep(2.1) # so the defining files are old enough not to trigger mtime criterion
        @eval using ModFILE
        @test ModFILE.mf() == (joinpath(dn, "ModFILE.jl"), 1)
        sleep(0.1)
        open(joinpath(dn, "ModFILE.jl"), "w") do io
            println(io, """
module ModFILE

mf() = @__FILE__, 2

end
""")
        end
        yry()
        @test ModFILE.mf() == (joinpath(dn, "ModFILE.jl"), 2)
        rm_precompile("ModFILE")
    end

    # issue #8
    @testset "Module docstring" begin
        testdir = randtmp()
        mkdir(testdir)
        push!(to_remove, testdir)
        push!(LOAD_PATH, testdir)
        dn = joinpath(testdir, "ModDocstring", "src")
        mkpath(dn)
        open(joinpath(dn, "ModDocstring.jl"), "w") do io
            println(io, """
" Ahoy! "
module ModDocstring

include("dependency.jl")

f() = 1

end
""")
        end
        open(joinpath(dn, "dependency.jl"), "w") do io
            println(io, "")
        end
        sleep(2.1) # so the defining files are old enough not to trigger mtime criterion
        @eval using ModDocstring
        sleep(2)
        @test ModDocstring.f() == 1
        ds = @doc(ModDocstring)
        @test get_docstring(ds) == "Ahoy! "

        sleep(0.1)  # ensure watching is set up
        open(joinpath(dn, "ModDocstring.jl"), "w") do io
            println(io, """
" Ahoy! "
module ModDocstring

include("dependency.jl")

f() = 2

end
""")
        end
        yry()
        @test ModDocstring.f() == 2
        ds = @doc(ModDocstring)
        @test get_docstring(ds) == "Ahoy! "

        open(joinpath(dn, "ModDocstring.jl"), "w") do io
            println(io, """
" Hello! "
module ModDocstring

include("dependency.jl")

f() = 3

end
""")
        end
        yry()
        @test ModDocstring.f() == 3
        ds = @doc(ModDocstring)
        @test get_docstring(ds) == "Hello! "

        pop!(LOAD_PATH)
    end

    @testset "Line numbers" begin
        # issue #27
        testdir = randtmp()
        mkdir(testdir)
        push!(to_remove, testdir)
        push!(LOAD_PATH, testdir)
        modname = "LineNumberMod"
        dn = joinpath(testdir, modname, "src")
        mkpath(dn)
        open(joinpath(dn, modname*".jl"), "w") do io
            println(io, """
module $modname
include("incl.jl")
end
""")
        end
        open(joinpath(dn, "incl.jl"), "w") do io
            println(io, """
0
0
1
2
3
4
5
6
7
8

function foo(x)
    return x+5
end

foo(y::Int) = y-51
""")
        end
        sleep(2.1) # so the defining files are old enough not to trigger mtime criterion
        @eval using LineNumberMod
        lines = Int[]
        files = String[]
        for m in methods(LineNumberMod.foo)
            push!(files, String(m.file))
            push!(lines, m.line)
        end
        @test all(f->endswith(string(f), "incl.jl"), files)
        sleep(0.1)  # ensure watching is set up
        open(joinpath(dn, "incl.jl"), "w") do io
            println(io, """
0
0
1
2
3
4
5
6
7
8

function foo(x)
    return x+6
end

foo(y::Int) = y-51
""")
        end
        yry()
        for m in methods(LineNumberMod.foo)
            @test endswith(string(m.file), "incl.jl")
            @test m.line ∈ lines
        end
    end

    # Issue #43
    @testset "New submodules" begin
        testdir = randtmp()
        mkdir(testdir)
        push!(to_remove, testdir)
        push!(LOAD_PATH, testdir)
        dn = joinpath(testdir, "Submodules", "src")
        mkpath(dn)
        open(joinpath(dn, "Submodules.jl"), "w") do io
            println(io, """
module Submodules
f() = 1
end
""")
        end
        sleep(2.1) # so the defining files are old enough not to trigger mtime criterion
        @eval using Submodules
        @test Submodules.f() == 1
        sleep(0.1)  # ensure watching is set up
        open(joinpath(dn, "Submodules.jl"), "w") do io
            println(io, """
module Submodules
f() = 1
module Sub
g() = 2
end
end
""")
        end
        yry()
        @test Submodules.f() == 1
        @test Submodules.Sub.g() == 2
    end

    @testset "Method deletion" begin
        Core.eval(Base, :(revisefoo(x::Float64) = 1)) # to test cross-module method scoping
        testdir = randtmp()
        mkdir(testdir)
        push!(to_remove, testdir)
        push!(LOAD_PATH, testdir)
        dn = joinpath(testdir, "MethDel", "src")
        mkpath(dn)
        open(joinpath(dn, "MethDel.jl"), "w") do io
            println(io, """
__precompile__(false)   # "clean" Base doesn't have :revisefoo
module MethDel
f(x) = 1
f(x::Int) = 2
g(x::Vector{T}, y::T) where T = 1
g(x::Array{T,N}, y::T) where N where T = 2
g(::Array, ::Any) = 3
h(x::Array{T}, y::T) where T = g(x, y)
k(::Int; badchoice=1) = badchoice
Base.revisefoo(x::Int) = 2
struct Private end
Base.revisefoo(::Private) = 3

dfltargs(x::Int8, y::Int=0, z::Float32=1.0f0) = x+y+z

hasmacro1(@nospecialize(x)) = x
hasmacro2(@nospecialize(x::Int)) = x
hasmacro3(@nospecialize(x::Int), y::Float64) = x

hasdestructure1(x, (count, name)) = name^count
hasdestructure2(x, (count, name)::Tuple{Int,Any}) = name^count

end
""")
        end
        @eval using MethDel
        @test MethDel.f(1.0) == 1
        @test MethDel.f(1) == 2
        @test MethDel.g(rand(3), 1.0) == 1
        @test MethDel.g(rand(3, 3), 1.0) == 2
        @test MethDel.g(Int[], 1.0) == 3
        @test MethDel.h(rand(3), 1.0) == 1
        @test MethDel.k(1) == 1
        @test MethDel.k(1; badchoice=2) == 2
        @test MethDel.hasmacro1(1) == 1
        @test MethDel.hasmacro2(1) == 1
        @test MethDel.hasmacro3(1, 0.0) == 1
        @test MethDel.hasdestructure1(0, (3, "hi")) == "hihihi"
        @test MethDel.hasdestructure2(0, (3, "hi")) == "hihihi"
        @test Base.revisefoo(1.0) == 1
        @test Base.revisefoo(1) == 2
        @test Base.revisefoo(MethDel.Private()) == 3
        @test MethDel.dfltargs(Int8(2)) == 3.0f0
        @test MethDel.dfltargs(Int8(2), 5) == 8.0f0
        @test MethDel.dfltargs(Int8(2), 5, -17.0f0) == -10.0f0
        sleep(0.1)  # ensure watching is set up
        open(joinpath(dn, "MethDel.jl"), "w") do io
            println(io, """
module MethDel
f(x) = 1
g(x::Array{T,N}, y::T) where N where T = 2
h(x::Array{T}, y::T) where T = g(x, y)
k(::Int; goodchoice=-1) = goodchoice
dfltargs(x::Int8, yz::Tuple{Int,Float32}=(0,1.0f0)) = x+yz[1]+yz[2]
end
""")
        end
        yry()
        @test MethDel.f(1.0) == 1
        @test MethDel.f(1) == 1
        @test MethDel.g(rand(3), 1.0) == 2
        @test MethDel.g(rand(3, 3), 1.0) == 2
        @test_throws MethodError MethDel.g(Int[], 1.0)
        @test MethDel.h(rand(3), 1.0) == 2
        @test_throws MethodError MethDel.k(1; badchoice=2)
        @test MethDel.k(1) == -1
        @test MethDel.k(1; goodchoice=10) == 10
        @test_throws MethodError MethDel.hasmacro1(1)
        @test_throws MethodError MethDel.hasmacro2(1)
        @test_throws MethodError MethDel.hasmacro3(1, 0.0)
        @test_throws MethodError MethDel.hasdestructure1(0, (3, "hi"))
        @test_throws MethodError MethDel.hasdestructure2(0, (3, "hi"))
        @test Base.revisefoo(1.0) == 1
        @test_throws MethodError Base.revisefoo(1)
        @test_throws MethodError Base.revisefoo(MethDel.Private())
        @test MethDel.dfltargs(Int8(2)) == 3.0f0
        @test MethDel.dfltargs(Int8(2), (5,-17.0f0)) == -10.0f0
        @test_throws MethodError MethDel.dfltargs(Int8(2), 5) == 8.0f0
        @test_throws MethodError MethDel.dfltargs(Int8(2), 5, -17.0f0) == -10.0f0

        Base.delete_method(first(methods(Base.revisefoo)))
    end

    @testset "Pkg exclusion" begin
        push!(Revise.dont_watch_pkgs, :Example)
        push!(Revise.silence_pkgs, :Example)
        @eval import Example
        for k in keys(Revise.fileinfos)
            if occursin("Example", k)
                error("Should not track files in Example")
            end
        end
        # Ensure that silencing works
        sfile = Revise.silencefile[]  # remember the original
        try
            sfiletemp = tempname()
            Revise.silencefile[] = sfiletemp
            Revise.silence("GSL")
            @test isfile(sfiletemp)
            pkgs = readlines(sfiletemp)
            @test any(p->p=="GSL", pkgs)
            rm(sfiletemp)
        finally
            Revise.silencefile[] = sfile
        end
    end

    @testset "Manual track" begin
        srcfile = joinpath(tempdir(), randtmp()*".jl")
        open(srcfile, "w") do io
            print(io, """
revise_f(x) = 1
""")
        end
        includet(srcfile)
        @test revise_f(10) == 1
        sleep(0.1)
        open(srcfile, "w") do io
            print(io, """
revise_f(x) = 2
""")
        end
        yry()
        @test revise_f(10) == 2
        push!(to_remove, srcfile)

        # Do it again with a relative path
        curdir = pwd()
        cd(tempdir())
        srcfile = randtmp()*".jl"
        open(srcfile, "w") do io
            print(io, """
        revise_floc(x) = 1
        """)
        end
        include(joinpath(pwd(), srcfile))
        @test revise_floc(10) == 1
        Revise.track(srcfile)
        sleep(0.1)
        open(srcfile, "w") do io
            print(io, """
        revise_floc(x) = 2
        """)
        end
        yry()
        @test revise_floc(10) == 2
        push!(to_remove, joinpath(tempdir(), srcfile))
        cd(curdir)
    end

    @testset "Auto-track user scripts" begin
        srcfile = joinpath(tempdir(), randtmp()*".jl")
        push!(to_remove, srcfile)
        open(srcfile, "w") do io
            println(io, "revise_g() = 1")
        end
        # By default user scripts are not tracked
        include(srcfile)
        yry()
        @test revise_g() == 1
        sleep(0.1)
        open(srcfile, "w") do io
            println(io, "revise_g() = 2")
        end
        yry()
        @test revise_g() == 1
        # Turn on tracking of user scripts
        empty!(Revise.included_files)  # don't track files already loaded (like this one)
        Revise.tracking_Main_includes[] = true
        try
            srcfile = joinpath(tempdir(), randtmp()*".jl")
            push!(to_remove, srcfile)
            open(srcfile, "w") do io
                println(io, "revise_g() = 1")
            end
            include(srcfile)
            yry()
            @test revise_g() == 1
            sleep(0.1)
            open(srcfile, "w") do io
                println(io, "revise_g() = 2")
            end
            yry()
            @test revise_g() == 2
        finally
            Revise.tracking_Main_includes[] = false  # restore old behavior
        end
    end

    @testset "Distributed" begin
        allworkers = [myid(); addprocs(2)]
        @everywhere using Revise
        dirname = randtmp()
        mkdir(dirname)
        @everywhere push_LOAD_PATH!(dirname) = push!(LOAD_PATH, dirname)
        for p in allworkers
            remotecall_wait(push_LOAD_PATH!, p, dirname)
        end
        push!(to_remove, dirname)
        modname = "ReviseDistributed"
        dn = joinpath(dirname, modname, "src")
        mkpath(dn)
        open(joinpath(dn, modname*".jl"), "w") do io
            println(io, """
module ReviseDistributed

f() = π
g(::Int) = 0

end
""")
        end
        sleep(2.1)   # so the defining files are old enough not to trigger mtime criterion
        using ReviseDistributed
        @everywhere using ReviseDistributed
        for p in allworkers
            @test remotecall_fetch(ReviseDistributed.f, p)    == π
            @test remotecall_fetch(ReviseDistributed.g, p, 1) == 0
        end
        sleep(0.1)
        open(joinpath(dn, modname*".jl"), "w") do io
            println(io, """
module ReviseDistributed

f() = 3.0

end
""")
        end
        yry()
        sleep(1.0)
        @test_throws MethodError ReviseDistributed.g(1)
        for p in allworkers
            @test remotecall_fetch(ReviseDistributed.f, p)    == 3.0
            p == myid() && continue
            @test_throws RemoteException remotecall_fetch(ReviseDistributed.g, p, 1)
        end
    end

    @testset "Git" begin
        if haskey(ENV, "CI")   # if we're doing CI testing (Travis, Appveyor, etc.)
            # First do a full git checkout of a package (we'll use Revise itself)
            @warn "checking out a development copy of Revise for testing purposes"
            pkg = Pkg.Types.PackageSpec("Revise")
            pkg.repo = Pkg.Types.GitRepo("", "")
            Pkg.API.add_or_develop(Pkg.Types.Context(), [pkg], mode=:develop)
        end
        loc = Base.find_package("Revise")
        if occursin("dev", loc)
            repo, path = Revise.git_repo(loc)
            @test repo != nothing
            files = Revise.git_files(repo)
            @test "README.md" ∈ files
            src = Revise.git_source(loc, "946d588328c2eb5fe5a56a21b4395379e41092e0")
            @test startswith(src, "__precompile__")
            src = Revise.git_source(loc, "eae5e000097000472280e6183973a665c4243b94") # 2nd commit in Revise's history
            @test src == "module Revise\n\n# package code goes here\n\nend # module\n"
        else
            @warn "skipping git tests because Revise is not under development"
        end
    end

    @testset "Recipes" begin
        # Tracking Base
        Revise.track(Base)
        @test any(k->endswith(k, "number.jl"), keys(Revise.fileinfos))
        @test length(filter(k->endswith(k, "file.jl"), keys(Revise.fileinfos))) == 2

        # Determine whether a git repo is available. Travis & Appveyor do not have this.
        repo, path = Revise.git_repo(Revise.juliadir)
        if repo != nothing
            # Tracking Core.Compiler
            Revise.track(Core.Compiler)
            @test any(k->endswith(k, "compiler.jl"), keys(Revise.fileinfos))

            # Tracking stdlibs
            Revise.track(Unicode)
            @test any(k->endswith(k, "Unicode.jl"), keys(Revise.fileinfos))
            @test Revise.get_def(first(methods(Unicode.isassigned))) isa Revise.RelocatableExpr

            # Test that we skip over files that don't end in ".jl"
            logs, _ = Test.collect_test_logs() do
                Revise.track(REPL)
            end
            @test isempty(logs)
        else
            @warn "skipping Core.Compiler and stdlibs tests due to lack of git repo"
        end
    end

    @testset "Cleanup" begin
        warnfile = randtmp()
        open(warnfile, "w") do io
            redirect_stderr(io) do
                for name in to_remove
                    try
                        rm(name; force=true, recursive=true)
                        deleteat!(LOAD_PATH, findall(LOAD_PATH .== name))
                    catch
                    end
                end
                yry()
            end
        end
        if !Sys.isapple()
            @test occursin("is not an existing directory", read(warnfile, String))
        end
        rm(warnfile)
    end
end

GC.gc(); GC.gc(); GC.gc()   # work-around for https://github.com/JuliaLang/julia/issues/28306


text:
module Haversine

using Parameters

R = 6.371e6 # earth's volumetric mean radius in meters

@with_kw struct GeoLocation
    λ::Real
    ϕ::Real
end

function HaversineDistance(p1::GeoLocation, p2::GeoLocation)::Float64
    λ1, ϕ1 = p1.λ, p1.ϕ
    λ2, ϕ2 = p2.λ, p2.ϕ
    Δϕ = ϕ2 - ϕ1
    Δλ = λ2 - λ1
    a = sind(Δϕ / 2)^2 + cosd(ϕ1) * cosd(ϕ2) * sind(Δλ / 2)^2
    c = 2 * atan(sqrt(a), sqrt(1 - a))
    return c * R
end

function HaversineBearing(p1::GeoLocation, p2::GeoLocation)::Float64
    λ1, ϕ1 = p1.λ, p1.ϕ
    λ2, ϕ2 = p2.λ, p2.ϕ
    Δλ = λ2 - λ1
    θ = atand(sind(Δλ) * cosd(ϕ2), cosd(ϕ1) * sind(ϕ2) - sind(ϕ1) * cosd(ϕ2) * cosd(Δλ))
    return θ
end

function HaversineDestination(geopoint::GeoLocation, θ::Real, d::Real)::GeoLocation
    λ1, ϕ1 = geopoint.λ, geopoint.ϕ
    δ = d / R
    ϕ2 = asind(sind(ϕ1) * cos(δ) + cosd(ϕ1) * sin(δ) * cosd(θ))
    λ2 = λ1 + atand(sind(θ) * sin(δ) * cosd(ϕ1), cos(δ) - sind(ϕ1) * sind(ϕ2))
    return GeoLocation(λ=λ2, ϕ=ϕ2)
end

function HaversineDistance(p1::AbstractArray{GeoLocation}, p2::AbstractArray{GeoLocation})::Vector{Float64}
    return map(HaversineDistance, p1, p2)
end

function HaversineBearing(p1::AbstractArray{GeoLocation}, p2::AbstractArray{GeoLocation})::Vector{Float64}
    return map(HaversineBearing, p1, p2)
end

function HaversineDestination(
        p::Union{AbstractArray, GeoLocation},
        θ::Union{AbstractArray, Real},
        d::Union{AbstractArray, Real}
    )::Vector{GeoLocation}
    if isa(p, GeoLocation)
        return broadcast((y, z) -> HaversineDestination(p, y, z), θ, d)
    end
    return broadcast((x, y, z) -> HaversineDestination(x, y, z), p, θ, d)
end

export GeoLocation, HaversineDistance, HaversineBearing, HaversineDestination

end


text:
##### MEAN-VARIANCE OPTIMIZATION
#This is not supposed to produce great results, since mean and covariance
#is supposed to be very exact before good weights are found
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")
#Esben's path
path = "$(homedir())/Documents/GitHub/Thesis/Data/MonthlyReturns"

mainData = loadIndexDataNoDurLOGReturn(path)

mainDataArr = Array(mainData)

startRow = 1
endRow = 100
trainingObservations = 100
mainXarr = mainDataArr[startRow:endRow,1:10]
indexes = size(mainXarr)[2]

#### WITH RISK AVERSION

nRows = size(mainDataArr)[1]
nCols = size(mainDataArr)[2]

##LOOP TO TEST IT
portfolioMVTotal = []
portfolioMVTotal = push!(portfolioMVTotal, 1)
portfolioMVIndividual = []
portfolioMVIndividual = push!(portfolioMVIndividual, 0)
wMVtracker = zeros((nRows-trainingObservations-1),10)

w1N = repeat([0.1], outer = 10)
portfolio1NTotal = []
portfolio1NTotal = push!(portfolio1NTotal, 1)
portfolio1NIndividual = []
portfolio1NIndividual = push!(portfolio1NIndividual, 0)

for rStart = 1:(nRows-trainingObservations-1)
    startRow = rStart
    endRow = (rStart-1)+trainingObservations
    mainXarr = mainDataArr[startRow:endRow,1:10]
    indexes = size(mainXarr)[2]

    #creating the covariance matrix of the training data
    Sigma = cov(mainXarr)

    #creating the mean of the returns for SAA (forecasts should be used instead)
    Mu = mean(mainXarr, 1)'

    #setting risk aversion.
    gamma = 10

    M = JuMP.Model(solver = GurobiSolver(OutputFlag = 0))
    @variables M begin
            w[1:indexes]
    end
    @objective(M,Min, gamma*w'*Sigma*w-(Mu'*w)[1]) #efficient frontier
    @constraint(M, 0 .<= w)
    @constraint(M, sum(w[i] for i=1:indexes) == 1)

    solve(M)
    wStar = getvalue(w)
    for i = 1:length(wStar)
        if wStar[i] >=0.5
            println("Index $i is chosen")
        end
    end
    wMVtracker[rStart,:] = wStar

    #creating the forecastRow in arithmetic returns (not log returns)
    #then we can multiply and sum across the row
    forecastRow = exp(mainDataArr[(endRow+1):(endRow+1),1:10])

    periodReturn = Array(forecastRow*wStar)[1]
    currentLength = length(portfolioMVTotal)
    portfolioMVTotal = push!(portfolioMVTotal, portfolioMVTotal[currentLength]*periodReturn)
    portfolioMVIndividual = push!(portfolioMVIndividual, periodReturn)

    period1NReturn = Array(forecastRow*w1N)[1]
    current1NLength = length(portfolio1NTotal)
    portfolio1NTotal = push!(portfolio1NTotal, portfolio1NTotal[current1NLength]*period1NReturn)
    portfolio1NIndividual = push!(portfolio1NIndividual, period1NReturn)
end

#Checking the total return of the portfolio by checking last index
portfolioMVTotal[length(portfolioMVTotal)]
portfolio1NTotal[length(portfolio1NTotal)]

portfolioMVTotal = convert(Array{Float64,1}, portfolioMVTotal)
portfolioMVIndividual = convert(Array{Float64,1}, portfolioMVIndividual)
(mean(portfolioMVIndividual)-1)*100
std(portfolioMVIndividual)
SharpeRatioMV = (mean(portfolioMVIndividual)-1)/std(portfolioMVIndividual)

portfolio1NTotal = convert(Array{Float64,1}, portfolio1NTotal)
portfolio1NIndividual = convert(Array{Float64,1}, portfolio1NIndividual)
(mean(portfolio1NIndividual)-1)*100
std(portfolio1NIndividual)
SharpeRatio1N = (mean(portfolio1NIndividual)-1)/std(portfolio1NIndividual)

combinedPortfolios = hcat(portfolioMVTotal, portfolio1NTotal, portfolioMVIndividual, portfolio1NIndividual)
writedlm("MV100obsG10versus1N.csv", combinedPortfolios,",")

writedlm("MV100obsPortfolioWeightsG10.csv", wMVtracker,",")

function meanVariancePPDOptimization(fullX, startRow, endRow, gammaValue, meanForecasts)
    mainXarr = fullX[startRow:endRow, 1:10]
    w1N = repeat([0.1], outer = 10)

    #creating the covariance matrix of the training data
    Sigma = cov(mainXarr)

    #using the PPD mean forecasts
    Mu = meanForecasts

    #setting risk aversion.
    gamma = gammaValue

    M = JuMP.Model(solver = GurobiSolver(OutputFlag = 0))
    @variables M begin
            w[1:indexes]
    end
    @objective(M,Min, gamma*w'*Sigma*w-(Mu'*w)[1]) #efficient frontier
    @constraint(M, 0 .<= w)
    @constraint(M, sum(w[i] for i=1:indexes) == 1)

    solve(M)
    wStar = getvalue(w)

    #creating the forecastRow in arithmetic returns (not log returns)
    #then we can multiply and sum across the row
    forecastRow = exp(fullX[(endRow+1):(endRow+1),1:10])

    periodReturn = Array(forecastRow*wStar)[1]
    period1NReturn = Array(forecastRow*w1N)[1]

    return periodReturn, period1NReturn, wStar
end

"""
Keeps track of returns for both the proposed method and 1/N.
Needs the overall array as well as the two next returns to push
"""

function trackPortfolioReturns(ReturnOverall, individualReturn, newPeriodReturn, Return1NOverall, Individual1NReturn, period1NReturn)
    currentLength = length(ReturnOverall)
    ReturnOverall = push!(ReturnOverall, ReturnOverall[currentLength]*newPeriodReturn)
    individualReturn = push!(individualReturn, newPeriodReturn)

    current1NLength = length(Return1NOverall)
    Return1NOverall = push!(Return1NOverall, Return1NOverall[current1NLength]*period1NReturn)
    Individual1NReturn = push!(Individual1NReturn, period1NReturn)
end

function trackPortfolioWeights(weights, periodWeights, iteration)
    weights[iteration, 1:10] = periodWeights
end

"""
Utilizing LASSO to create a classification; 1, it goes up, 0, it goes down
The portfolio is then a 1/K, where K is the amount of indexes predicted to go up
No weighting between them or variance optimization, hence quite basic
"""
function ClassificationPortfolio(fullX, firstRow, indicators)

end

####SOCP FORM OF MEAN-VARIANCE
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")
#Esben's path
path = "$(homedir())/Documents/GitHub/Thesis/Data/MonthlyReturns"

mainData = loadIndexDataNoDurLOGReturn(path)

mainDataArr = Array(mainData)
indexes = 10

#### WITH RISK AVERSION

nRows = size(mainDataArr)[1]
nCols = size(mainDataArr)[2]

##LOOP TO TEST IT
portfolioMVTotal = []
portfolioMVTotal = push!(portfolioMVTotal, 1)
portfolioMVIndividual = []
portfolioMVIndividual = push!(portfolioMVIndividual, 0)
wMVtracker = zeros((nRows-trainingObservations-1),10)

w1N = repeat([0.1], outer = 10)
portfolio1NTotal = []
portfolio1NTotal = push!(portfolio1NTotal, 1)
portfolio1NIndividual = []
portfolio1NIndividual = push!(portfolio1NIndividual, 0)

for rStart = 1:(nRows-trainingObservations-1)
    startRow = rStart
    endRow = (rStart-1)+trainingObservations
    mainXarr = mainDataArr[startRow:endRow,1:10]
    indexes = size(mainXarr)[2]

    #creating the covariance matrix of the training data
    Sigma = cov(mainXarr)
    #A=U^(T)U where U is upper triangular with real positive diagonal entries
    F = lufact(Sigma)
    U = F[:U]  #Cholesky factorization of Sigma
    #creating the mean of the returns for SAA (forecasts should be used instead)
    Mu = mean(mainXarr, 1)'

    #setting risk aversion.
    gamma = 1000

    M = JuMP.Model(solver = GurobiSolver(OutputFlag = 0))
    @variables M begin
            w[1:indexes]
            u[1:indexes]
            z
            y
    end

    @objective(M,Min, gamma*y - (Mu'*w)[1])
    @constraint(M, 0 .<= w)
    @constraint(M, sum(w[i] for i=1:indexes) == 1)
    @constraint(M, norm([2*U'*w;y-1]) <= y+1)
    solve(M)
    wStar = getvalue(w)

    wMVtracker[rStart,:] = wStar

    #creating the forecastRow in arithmetic returns (not log returns)
    #then we can multiply and sum across the row
    forecastRow = (exp10(mainDataArr[(endRow+1):(endRow+1),1:10])-1)*100

    periodReturn = Array(forecastRow*wStar)[1]
    currentLength = length(portfolioMVTotal)
    portfolioMVTotal = push!(portfolioMVTotal, portfolioMVTotal[currentLength]*periodReturn)
    portfolioMVIndividual = push!(portfolioMVIndividual, periodReturn)

    period1NReturn = Array(forecastRow*w1N)[1]
    current1NLength = length(portfolio1NTotal)
    portfolio1NTotal = push!(portfolio1NTotal, portfolio1NTotal[current1NLength]*period1NReturn)
    portfolio1NIndividual = push!(portfolio1NIndividual, period1NReturn)
end

#Checking the total return of the portfolio by checking last index
portfolioMVTotal[length(portfolioMVTotal)]
portfolio1NTotal[length(portfolio1NTotal)]

portfolioMVTotal = convert(Array{Float64,1}, portfolioMVTotal)
portfolioMVIndividual = convert(Array{Float64,1}, portfolioMVIndividual)
(mean(portfolioMVIndividual)-1)*100
std(portfolioMVIndividual)
SharpeRatioMV = (mean(portfolioMVIndividual)-1)/std(portfolioMVIndividual)

portfolio1NTotal = convert(Array{Float64,1}, portfolio1NTotal)
portfolio1NIndividual = convert(Array{Float64,1}, portfolio1NIndividual)
(mean(portfolio1NIndividual)-1)*100
std(portfolio1NIndividual)
SharpeRatio1N = (mean(portfolio1NIndividual)-1)/std(portfolio1NIndividual)

writedlm("MVWeights2.csv", wMVtracker,",")
2.98

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
#####
##### JSON Serialization
#####

import JSON3

for T in [
    # Reports
    Report.Loss, Report.LearningStatus, Report.Evaluation, Report.Checkpoint,
    Report.Learning, Report.Samples, Report.StageSamples, Report.Memory,
    Report.SelfPlay, Report.Perfs, Report.Iteration, Report.Initial,
    # Network Hyperparameters
    NetLib.SimpleNetHP, NetLib.ResNetHP, GraphNeuralNetworkLib.SimpleGraphNetHP,
    # Parameters
    Params, SelfPlayParams, LearningParams, ArenaParams,
    SimParams, MctsParams, MemAnalysisParams,
    # Optimisers
    CyclicNesterov, Adam,
    # Schedules
    ConstSchedule, PLSchedule, StepSchedule
  ]
  @eval JSON3.StructType(::Type{<:$T}) = JSON3.Struct()
end

# Abstract types

JSON3.StructType(::Type{OptimiserSpec}) = JSON3.AbstractType()
JSON3.subtypekey(::Type{OptimiserSpec}) = :type
JSON3.subtypes(::Type{OptimiserSpec}) =
  (adam=Adam, cyclic_nesterov=CyclicNesterov)

JSON3.StructType(::Type{<:AbstractSchedule}) = JSON3.AbstractType()
JSON3.subtypekey(::Type{<:AbstractSchedule}) = :type
JSON3.subtypes(::Type{<:AbstractSchedule}) =
  (piecewise_linear=PLSchedule, step=StepSchedule, constant=ConstSchedule)

# TODO: for the subtypes of abstract types above, the `subtypekey` field
# is missing so deserialization won't work.


text:
"""
    abstract type DatasetRecipe

A recipe that contains configuration for loading a data container.
Calling it with a path returns a data container and the blocks
that each sample is made of.

## Examples

For example implementations, see [`Vision.ImageFolders`](#).

## Extending

### Interface

- [`loadrecipe`](#)`(::DatasetRecipe, args...; kwargs...) -> (data, blocks)`
    This loads a data container the [`Block`]s that each observation corresponds to.
    For most recipes the only argument beside the recipe is a path to a folder on disk.
- [`recipeblocks`](#)`(::Type{DatasetRecipe}) -> TBlocks`
    The type of `blocks` returned by `loadrecipe`. Should be as specific as possible.
    Used for discovery.

### Invariants

Given

```julia
data, blocks = loadrecipe(recipe, args...; kwargs...)
```

the following must hold:

- `∀i ∈ 1:numobs(data): checkblock(blocks, data[i])`, i.e.
    `data` must be a data container of observations that are valid `blocks`.
- `numobs(data) ≥ 1`, i.e. there is at least one observation if the data was loaded
    without error.
"""
abstract type DatasetRecipe end

"""
    loadrecipe(recipe, path)

Load a recipe from a path. Return a data container `data` and concrete
`blocks`.
"""
function loadrecipe end

"""
    recipeblocks(TRecipe) -> TBlocks
    recipeblocks(recipe) -> TBlocks

Return the `Block` _types_ for the data container that recipe
type `TRecipe` creates. Does not return `Block` instances as the exact
configuration may not be known until the dataset is being
loaded.

#### Examples

```julia
recipeblocks(ImageFolders) == Tuple{Image{2}, Label}
```
"""
recipeblocks(::R) where {R <: DatasetRecipe} = recipeblocks(R)

function testrecipe(recipe::Datasets.DatasetRecipe, path::AbstractPath)
    data, blocks = loadrecipe(recipe, path)
    return testrecipe(recipe, data, blocks)
end

function testrecipe(recipe::Datasets.DatasetRecipe, data, blocks)
    # `blocks` must be an instance of `recipeblocks(recipe)`
    @test blocks isa Datasets.recipeblocks(recipe)

    # Observations must be compatible with `blocks`
    @test checkblock(blocks, getobs(data, 1))
end


text:
@testset "Poly" begin

  function _test_sturm()
    s = rand(1:20)
    Zx, x = FlintZZ["x"]
    M = random_symmetric_matrix(s)
    f = charpoly(Zx, M)
    while !issquarefree(f) || iszero(coeff(f, 0))
      M = random_symmetric_matrix(s)
      f = charpoly(Zx, M)
    end
    npos = Hecke.number_positive_roots(f)
    ff=factor(f)
    sgtpos=0
    sgtneg=0
    for (h,v) in ff.fac
      if degree(h)==1
        if coeff(h,0)>0
          sgtneg+=v
        else
          sgtpos+=v
        end
        continue
      end
      p=64
      while p<4096
        sgtposf=0
        sgtnegf=0
        R=AcbField(p, cached = false)
        Rx=AcbPolyRing(R, Symbol("x"), false)
        g=Rx(h)
        l=roots(g)
        for i=1:length(l)
          y=real(l[i])
          if ispositive(y)
            sgtposf+=1
          end
          if isnegative(y)
            sgtnegf+=1
          end
        end
        if sgtposf+sgtnegf==degree(h)
          sgtpos+=sgtposf*v
          sgtneg+=sgtnegf*v
          break
        else
          p*=2
        end
      end
      if p > 4096
        error("Precision issue")
      end
    end
    return sgtpos == npos
  end

  function random_symmetric_matrix(x::Int)
    M = zero_matrix(FlintZZ, x, x)
    for i = 1:x
      for j= i:x
        a = rand(1:5)
        M[i, j] = a
        M[j, i] = a
      end
    end
    return M
  end

  for i = 1:20
    @test _test_sturm()
  end

end


text:
## Quite a few of these are borrowed from Flux.jl

# Misc
nfan() = 1, 1 # fan_in, fan_out
nfan(n) = 1, n # A vector is treated as a n×1 matrix
nfan(n_out, n_in) = n_in, n_out # In case of Dense kernels: arranged as matrices
nfan(dims::Tuple) = nfan(dims...)
nfan(dims...) = prod(dims[1:(end - 2)]) .* (dims[end - 1], dims[end]) # In case of convolution kernels

# Neural Network Initialization
## NOTE: Would be great if these could be moved into its own package and NN frameworks
##       could just import it.
"""
    zeros32(rng::AbstractRNG, size...) = zeros(Float32, size...)

Return an `Array{Float32}` of zeros of the given `size`. (`rng` is ignored)
"""
zeros32(rng::AbstractRNG, args...; kwargs...) = zeros(rng, Float32, args...; kwargs...)

"""
    ones32(rng::AbstractRNG, size...) = ones(Float32, size...)

Return an `Array{Float32}` of ones of the given `size`. (`rng` is ignored)
"""
ones32(rng::AbstractRNG, args...; kwargs...) = ones(rng, Float32, args...; kwargs...)

Base.zeros(rng::AbstractRNG, args...; kwargs...) = zeros(args...; kwargs...)
Base.ones(rng::AbstractRNG, args...; kwargs...) = ones(args...; kwargs...)

"""
    glorot_uniform(rng::AbstractRNG, size...; gain = 1)

Return an `Array{Float32}` of the given `size` containing random numbers drawn from a uniform distribution on the interval ``[-x, x]``, where `x = gain * sqrt(6 / (fan_in + fan_out))`. This method is described in [1] and also known as Xavier initialization.

# References

[1] Glorot, Xavier, and Yoshua Bengio. "Understanding the difficulty of training deep feedforward neural networks." _Proceedings of the thirteenth international conference on artificial intelligence and statistics_. 2010.
"""
function glorot_uniform(rng::AbstractRNG, dims::Integer...; gain::Real=1)
    scale = Float32(gain) * sqrt(24.0f0 / sum(nfan(dims...)))
    return (rand(rng, Float32, dims...) .- 0.5f0) .* scale
end

"""
    glorot_normal(rng::AbstractRNG, size...; gain = 1)

Return an `Array{Float32}` of the given `size` containing random numbers drawn from a normal distribution with standard deviation `gain * sqrt(2 / (fan_in + fan_out))`. This method is described in [1] and also known as Xavier initialization.

# References

[1] Glorot, Xavier, and Yoshua Bengio. "Understanding the difficulty of training deep feedforward neural networks." _Proceedings of the thirteenth international conference on artificial intelligence and statistics_. 2010.
"""
function glorot_normal(rng::AbstractRNG, dims::Integer...; gain::Real=1)
    std = Float32(gain) * sqrt(2.0f0 / sum(nfan(dims...)))
    return randn(rng, Float32, dims...) .* std
end

# PRNG Handling
replicate(rng::AbstractRNG) = copy(rng)
replicate(rng::CUDA.RNG) = deepcopy(rng)

# Training Check
@inline istraining() = false
@inline istraining(::Val{training}) where {training} = training
@inline istraining(st::NamedTuple) = istraining(st.training)

# Linear Algebra
@inline _norm(x; dims=Colon()) = sqrt.(sum(abs2, x; dims=dims))
@inline function _norm_except(x::AbstractArray{T, N}, except_dim=N) where {T, N}
    _norm(x; dims=filter(i -> i != except_dim, 1:N))
end

# Convolution
function convfilter(rng::AbstractRNG, filter::NTuple{N, Integer},
                    ch::Pair{<:Integer, <:Integer};
                    init=glorot_uniform, groups=1) where {N}
    cin, cout = ch
    @assert cin % groups==0 "Input channel dimension must be divisible by groups."
    @assert cout % groups==0 "Output channel dimension must be divisible by groups."
    return init(rng, filter..., cin ÷ groups, cout)
end

expand(N, i::Tuple) = i
expand(N, i::Integer) = ntuple(_ -> i, N)

_maybetuple_string(pad) = string(pad)
_maybetuple_string(pad::Tuple) = all(==(pad[1]), pad) ? string(pad[1]) : string(pad)

# Padding
struct SamePad end

function calc_padding(lt, pad, k::NTuple{N, T}, dilation, stride) where {T, N}
    expand(Val(2 * N), pad)
end

function calc_padding(lt, ::SamePad, k::NTuple{N, T}, dilation, stride) where {N, T}
    # Ref: "A guide to convolution arithmetic for deep learning" https://arxiv.org/abs/1603.07285
    # Effective kernel size, including dilation
    k_eff = @. k + (k - 1) * (dilation - 1)
    # How much total padding needs to be applied?
    pad_amt = @. k_eff - 1
    # In case amount of padding is odd we need to apply different amounts to each side.
    return Tuple(mapfoldl(i -> [cld(i, 2), fld(i, 2)], vcat, pad_amt))
end

# Handling ComponentArrays
## NOTE: We should probably upsteam some of these
function Base.zero(c::ComponentArray{T, N, <:CuArray{T}}) where {T, N}
    ComponentArray(zero(getdata(c)), getaxes(c))
end

Base.vec(c::ComponentArray{T, N, <:CuArray{T}}) where {T, N} = getdata(c)

function Base.:-(x::ComponentArray{T, N, <:CuArray{T}}) where {T, N}
    ComponentArray(-getdata(x), getaxes(x))
end

function Base.similar(c::ComponentArray{T, N, <:CuArray{T}},
                      l::Vararg{Union{Integer, AbstractUnitRange}}) where {T, N}
    return similar(getdata(c), l)
end

function Functors.functor(::Type{<:ComponentArray}, c)
    return NamedTuple{propertynames(c)}(getproperty.((c,), propertynames(c))),
           ComponentArray
end

# Updating a monolithic vector is way faster than chunks of smaller ones
# Also helps in distributed settings
Optimisers.setup(opt, ps::ComponentArray) = Optimisers.setup(opt, getdata(ps))

function Optimisers.update(opt_state, ps::ComponentArray, gs::ComponentArray)
    opt_state, ps_new = Optimisers.update(opt_state, getdata(ps), getdata(gs))
    return opt_state, ComponentArray(ps_new, getaxes(ps))
end

function Optimisers.update!(st, ps::ComponentArray, gs::ComponentArray)
    st, ps_ = Optimisers.update!(st, NamedTuple(ps), NamedTuple(gs))
    return st, ComponentArray(ps_)
end

function ComponentArrays.make_carray_args(nt::NamedTuple)
    data, ax = ComponentArrays.make_carray_args(Vector, nt)
    data = length(data) == 0 ? Float32[] :
           (length(data) == 1 ? [data[1]] : reduce(vcat, data))
    return (data, ax)
end

## For being able to print empty ComponentArrays
function ComponentArrays.last_index(f::FlatAxis)
    nt = ComponentArrays.indexmap(f)
    length(nt) == 0 && return 0
    return ComponentArrays.last_index(last(nt))
end

ComponentArrays.recursive_length(nt::NamedTuple{(), Tuple{}}) = 0

# Return Nothing if field not present
function safe_getproperty(x::Union{ComponentArray, NamedTuple}, k::Symbol)
    k ∈ propertynames(x) && return getproperty(x, k)
    return nothing
end

# Getting typename
get_typename(::T) where {T} = Base.typename(T).wrapper

# For Normalization
@inline @generated safe_copy(x::T) where {T} = hasmethod(copy, (T,)) ? :(copy(x)) : :x

@inline @generated safe_vec(x::T) where {T} = hasmethod(vec, (T,)) ? :(vec(x)) : :x

@inline @inbounds function get_reshape_dims(sx::NTuple{N, <:Int},
                                            ly::Int)::typeof(sx) where {N}
    if ly == sx[N - 1]
        return ntuple(i -> i == N - 1 ? ly : 1, N)
    elseif N > 2 && ly == sx[N - 1] * sx[N - 2]
        return ntuple(i -> i == (N - 1) || i == (N - 2) ? sx[i] : 1, N)
    else
        error("Invalid Dimensions")
    end
end

@inline reshape_into_proper_shape(x::Nothing, y)::Nothing = x
@inline reshape_into_proper_shape(x, y)::typeof(y) = reshape(x,
                                                             get_reshape_dims(size(y),
                                                                              length(x)))

# RNN Utilities
@inline gate(h::Int, n::Int) = (1:h) .+ h * (n - 1)
@inline gate(x::AbstractVector, h::Int, n::Int) = view(x, gate(h, n))
@inline gate(x::AbstractMatrix, h::Int, n::Int) = view(x, gate(h, n), :)

"""
    multigate(x::AbstractArray, ::Val{N})

Split up `x` into `N` equally sized chunks (along dimension `1`).
"""
@inline multigate(x::AbstractArray, ::Val{N}) where {N} = gate.((x,), size(x, 1) ÷ N, 1:N)

# Val utilities
get_known(::Val{T}) where {T} = T


text:
using Distributions, StatsBase, Plots, LaTeXStrings; pyplot()

mu0, mu1, sd, tau  = 15, 18, 2, 17.5
dist0, dist1 = Normal(mu0,sd), Normal(mu1,sd)
grid = 5:0.1:25
h0grid, h1grid = tau:0.1:25, 5:0.1:tau

println("Probability of Type I error: ", ccdf(dist0,tau))
println("Probability of Type II error: ", cdf(dist1,tau))

plot(grid, pdf.(dist0,grid),
	c=:blue, label="Bolt type 15g")
plot!(h0grid, pdf.(dist0, h0grid),
	c=:blue, fa=0.2, fillrange=[0 1], label="")
plot!(grid, pdf.(dist1,grid),
	c=:green, label="Bolt type 18g")
plot!(h1grid, pdf.(dist1, h1grid),
	c=:green, fa=0.2, fillrange=[0 1], label="")
plot!([tau, 25],[0,0],
	c=:red, lw=3, label="Rejection region",
	xlims=(5, 25), ylims=(0,0.25) , legend=:topleft,
    xlabel="x", ylabel="Density")
annotate!([(16, 0.02, text(L"\beta")),(18.5, 0.02, text(L"\alpha")),
            (15, 0.21, text(L"H_0")),(18, 0.21, text(L"H_1"))])

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# inference on Distribution type D not needed any more
# function paramtypes(::Type{D}) where D<:Distribution
#     isconcretetype(D) || error("Expected a concrete distibution type," *
#         " Did you specify all type parameters, e.g. $D{Float64}?")
#     D.types
# end
# function tupleofvectype(::Type{D}) where D<:Distribution
#     # https://discourse.julialang.org/t/allocations-in-comprehensions/21309/3
#     V = Tuple{ntuple(
#         i -> Vector{Union{Missing,paramtypes(D)[i]}}, length(paramtypes(D)))...}
# end

"""
    vectuptotupvec(vectup)

Typesafe convert from Vector of Tuples to Tuple of Vectors.

# Arguments
* `vectup`: A Vector of identical Tuples

# Examples
```jldoctest; output=false, setup = :(using Distributions)
vectup = [(1,1.01, "string 1"), (2,2.02, "string 2")]
vectuptotupvec(vectup) == ([1, 2], [1.01, 2.02], ["string 1", "string 2"])
# output
true
```
"""
function vectuptotupvec(vectup::
    AbstractVector{<:Union{Missing,NTuple{N, Any}}}) where N
    Ti = nonmissingtype(eltype(vectup)).parameters
    npar = length(Ti)
    Tim = ntuple(i -> Union{Missing,Ti[i]}, npar)
    #TiVec = ntuple(i -> Vector{Tim[i]}, npar)
    #vectupT = Tuple{TiVec...}
    # need to help collect both with eltype and return type
    ntuple(i ->
        collect(Tim[i], passmissing(getindex).(vectup, i))::Vector{Tim[i]}, npar)
end

# function vectuptotupvec_old(vectup::AbstractVector{TT}) where
#     {TT<:Union{Missing,NTuple{N, Any}}} where N
#     v1 = first(skipmissing(vectup))
#     types = typeof.(v1)
#     imiss = findall(ismissing, vectup) # unfortunately allocating
#     if length(imiss) != 0
#         # replace missing by a tuple of correct type (replaced later be missings)
#         # because there is no missing of correct type during getindex
#         vectupc = mappedarray((x -> ismissing(x) ? v1 : x), vectup)
#         function f(i)
#             v = allowmissing(
#                 getindex.(vectupc,i))::Vector{Union{Missing,types[i]}}
#             v[imiss] .= missing
#             v
#         end
#         ntuple(f, length(v1))
#     else
#         ntuple((i ->
#             allowmissing(getindex.(vectup,i))::Vector{Union{Missing,types[i]}}
#         ), length(v1))
#     end
# end

"""
    AbstractDistributionVector{D <: Distribution}

Is any type able represent a vector of distribution of the same type.
This corresponds to a sequence of random variables, each characterized
by the same type of distribution but with different parameters.
This allows aggregating functions to work, for
example, computing the distribution of the sum of random variables by
[`sum(dv::AbstractDistributionVector)`](@ref).

It is parametrized by `D <: Distribution` defining the type of the distribution
used for all the random variables.

Items may be missing. Hence the element type of the iterator is
`Union{Missing,D}`.

AbstractDistributionVector
- is iterable
- has length and index access, i.e. `dv[i]::D`
- access to entire parameter vectors: `params(dv,Val(i))`
- conversion to Tuple of Vectors: `params(dv)`
- array of random numbers: `rand(n, dv)`: adding one
  dimension that represents across random variables
- query if entry is missing without needing to construct the
  distribution entry: `ismissing(dv,i)`:

Specific implementations,  need
to implement at minimum methods `length` and `getindex`, and `params`.

There are two standard implementations:
- [`SimpleDistributionVector`](@ref): fast indexing but slower `params` method
- [`ParamDistributionVector`](@ref): possible allocations in indexing but
  faster `params`

# Examples
```jldoctest; output = false, setup = :(using Distributions)
import LinearAlgebra: I
dmn1 = MvNormal([0,0,0], 1 * I)
dmn2 = MvNormal([1,1,1], 2 * I)
dv = SimpleDistributionVector(dmn1, dmn2, missing, missing);
sample = rand(dv,2);
# 4 distr, each 2 samples of length 3
size(sample) == (3,2,4)
# output
true
```
"""
abstract type AbstractDistributionVector{D <: Distribution} <: AbstractVector{Union{Missing,D}} end

size(dv::AbstractDistributionVector) = (length(dv),)

# the following already definded by AbstractVector{D}
#Base.eltype(::Type{<:AbstractDistributionVector{D}}) where D = Union{Missing,D}
#Base.ndims(::Type{<:AbstractDistributionVector}) = 1

# function Base.iterate(dv::AbstractDistributionVector, state=1)
#     state > length(dv) ? nothing : (dv[state], state+1)
# end

# function Base.iterate(rds::Iterators.Reverse{AbstractDistributionVector},
#     state=length(rds.itr))
#     state < 1 ? nothing : (rds.itr[state], state-1)
# end
# Base.firstindex(dv::AbstractDistributionVector) = 1
# Base.lastindex(dv::AbstractDistributionVector) = length(dv)

# function Base.getindex(dv::AbstractDistributionVector, i::Number)
#     @info "getindex(dv::AbstractDistributionVector, i::Number)"
#     dv[convert(Int, i)]
# end

function params(dv::AbstractDistributionVector, ::Val{i}) where i
    # need to help compile to determine the type of tupvec
    T = typeof(params(first(skipmissing(dv)))[i])
    Tm = Union{Missing,T}
    collect(Tm, passmissing(getindex).(passmissing(params).(dv),i))::Vector{Tm}
end

# extends Random
function rand(dv::AbstractDistributionVector, n::Integer)
    x1 = rand(first(skipmissing(dv)), n)
    xm = Fill(missing,size(x1))
    #xm = fill(missing,size(x1))
    fmiss(x)::Union{typeof(xm),typeof(x1)} = (ismissing(x) ? xm : rand(x,n))
    vecarr = convert(Vector{Union{typeof(xm),typeof(x1)}}, fmiss.(dv))::Vector{Union{typeof(xm),typeof(x1)}}
    VectorOfArray(vecarr)
end
function rand(dv::AbstractDistributionVector)
    x1 = rand(first(skipmissing(dv)))
    xm = Fill(missing,size(x1))
    #xm = fill(missing,size(x1))
    fmiss(x)::Union{typeof(xm),typeof(x1)} = (ismissing(x) ? xm : rand(x))
    vecarr = convert(Vector{Union{typeof(xm),typeof(x1)}}, fmiss.(dv))::Vector{Union{typeof(xm),typeof(x1)}}
    nonmissingtype(eltype(dv)) <: UnivariateDistribution ? vecarr : VectorOfArray(vecarr)
end
function rand!(vecarr::AbstractArray{T}, dv::AbstractDistributionVector) where T
    x1 = rand(first(skipmissing(dv)))
    xm = Fill(missing,size(x1))
    #xm = fill(missing,size(x1))
    fmiss(x)::Union{typeof(xm),typeof(x1)} = (ismissing(x) ? xm : rand(x))
    vecarr .= fmiss.(dv)
end
function ismissing(dv::AbstractDistributionVector, i::Int)
    ismissing(dv[i])
end

## SimpleDistributionVector
"""
    SimpleDistributionVector{D <: Distribution, V}

Is an Vector-of-Distribution based implementation of
[`AbstractDistributionVector`](@ref).

Vector of random var can be created by
- specifying the distributions as arguments.
```jldoctest; output = false, setup = :(using Distributions)
d1 = LogNormal(log(110), 0.25)
d2 = LogNormal(log(100), 0.15)
dv = SimpleDistributionVector(d1, d2, missing);
isequal(params(dv, Val(1)), [log(110), log(100), missing])
# output
true
```

- providing the Type of distribution and vectors of each parameter
```jldoctest; output = false, setup = :(using Distributions)
mu = [1.1,1.2,1.3]
sigma = [1.01, 1.02, missing]
dv = SimpleDistributionVector(LogNormal{eltype(mu)}, mu, sigma);
isequal(params(dv, Val(1)), [1.1,1.2,missing])
# output
true
```
Note that if one of the parameters is missing, then the entire entry of
the distribution is marked missing.

Since Distributions are stored directly, indexing passes a reference.
However, getting parameter vectors, required iterating all distributions,
and allocating a new vector.
"""
struct SimpleDistributionVector{D <: Distribution, V <: AbstractVector} <:
    AbstractDistributionVector{D}
    dvec::V
    # inner constructor checking ?
end

function SimpleDistributionVector(::Type{D}, dvec::V) where
{D<:Distribution,  V<:AbstractVector}
    isconcretetype(D) || error("Expected a concrete distibution type," *
        " Did you specify all type parameters, e.g. $D{Float64}?")
    Missing <: eltype(V) || error(
        "Expected type of parameters to allow for missing." *
        " Can you use 'allowmissing' in constructing the " *
        "SimpleDistributionVector?")
    eltype(V) <: Union{Missing, <:D} || error(
        "Expected type of parameters of 'Union{Missing, $(D)}' "*
        " but got $(eltype(V)).")
    SimpleDistributionVector{D, V}(dvec)
end

function SimpleDistributionVector(dv::Vararg{Union{Missing,D},N}) where
    {D<:Distribution, N}
    N == 0 && error(
        "Provide at least one argument, i.e. distribtution," *
        "i n SimpleDistributionVector(x...).")
    d1 = first(skipmissing(dv))
    dvec = collect(Union{Missing, typeof(d1)}, dv)::Vector{Union{Missing, typeof(d1)}}
    SimpleDistributionVector(D, allowmissing(dvec))
end

function SimpleDistributionVector(::Type{D}, pvec::Vararg{Any,N}) where
    {D<:Distribution, N}
    #Union{Missing,typeof(pvec[1])}
    pvecm = ntuple(i ->
        allowmissing(pvec[i])::Vector{Union{Missing,eltype(pvec[i])}}, N)
    # if one parameter has missing, the entire tuple must be set to missing
    dvec = allowmissing(collect(zip(pvecm...)))
    for i in axes(dvec,1)
        if any(ismissing.(dvec[i]))
            dvec[i] = missing
        end
    end
    Tm = Union{Missing, D}
    dv = collect(Tm,
        ismissing(x) ? missing : D(x...) for x in dvec)::Vector{Tm}
    SimpleDistributionVector(D, dv)
end

length(dv::SimpleDistributionVector) = length(dv.dvec)
IndexStyle(::Type{<:SimpleDistributionVector{D,V}}) where {D,V} = IndexStyle(V)
similar(dv::SimpleDistributionVector{D,V}) where {D,V} =
    SimpleDistributionVector{D,V}(similar(dv.dvec))

getindex(dv::SimpleDistributionVector,i::Int) = dv.dvec[i] #::Union{Missing, D}
#function getindex(dv::SimpleDistributionVector, I::Union{AbstractVector{<:Integer},AbstractVector{<:Bool}})
function getindex(dv::SimpleDistributionVector, I) # ambiguous
    # # Catesian indexing?: provided a tuple -> get only a single value in printing
    #length(I) != 1 && return dv[I[1]]
    #dvecsub = dv.dvec[I...]
    dvecsub = dv.dvec[I]
    # #@show I, length(I), typeof(I), typeof(dvecsub)
    # length(I) != 1 && return(dvecsub)
    typeof(dv)(dvecsub)
end
setindex!(dv::SimpleDistributionVector, d, i::Int) = dv.dvec[i] = d
#Base.setindex!(dv::SimpleDistributionVector, d, I::Vararg{Int, N}) where N = dv.dvec[I...] = d

# params(i) already defined as default in AbstractDistributionVector
# function StatsBase.params(dv::SimpleDistributionVector, i::Integer)
#    # mappedarray(e -> passmissing(getindex)(e,i), dv.dvec)
#    # currentl does not work, see
# https://github.com/JuliaArrays/MappedArrays.jl/issues/40
#    passmissing(getindex).(passmissing(params).(dv.dvec),i)
# end

function params(dv::SimpleDistributionVector)
    Ti = typeof(params(first(skipmissing(dv.dvec)))).parameters
    npar = length(Ti)
    Tim = ntuple(i -> Union{Missing, Ti[i]}, npar)
    f(x,i) = ismissing(x) ? missing : getindex(params(x),i)
    ntuple(i -> collect(Tim[i], f.(dv.dvec, i))::Vector{Tim[i]}, npar)
end

## ParamDistributionVector
"""
    ParamDistributionVector{D <: Distribution, V}

Is an Tuple of Vectors based implementation of
[`AbstractDistributionVector`](@ref).

Vector of random var can be created by
- specifying the distributions as arguments with some overhead of converting
  the Distributions to vectors of each parameter
```jldoctest; output = false, setup = :(using Distributions)
d1 = LogNormal(log(110), 0.25)
d2 = LogNormal(log(100), 0.15)
dv = ParamDistributionVector(d1, d2, missing);
isequal(params(dv, Val(1)), [log(110), log(100), missing])
# output
true
```

- providing the Type of distribution and vectors of each parameter
```jldoctest; output = false, setup = :(using Distributions)
mu = [1.1,1.2,1.3]
sigma = [1.01, 1.02, missing]
dv = ParamDistributionVector(LogNormal{eltype(mu)}, mu, sigma);
ismissing(dv[3])
isequal(params(dv, Val(1)), [1.1,1.2,1.3]) # third still not missing here
# output
true
```
Note that if one of the parameters for entry `i` is missing, then `dv[i]`
is missing.

Since distributions are stored by parameter vectors, the acces to these
vectors is just passing a reference.
Indexing, will create Distribution types.
"""
struct ParamDistributionVector{D <: Distribution, V <: Tuple} <:
    AbstractDistributionVector{D}
    params::V
    # inner constructor checking ?
end

function ParamDistributionVector(::Type{D}, params::V) where
{D<:Distribution,  V<:Tuple}
    isconcretetype(D) || error("Expected a concrete distibution type," *
        " Did you specify all type parameters, e.g. $D{Float64}?")
    all(map((x -> x <: AbstractVector), V.parameters)) || error(
        "Expected all entries in Tuple param to be AbstractVectors.")
    all(map((x -> Missing <: eltype(x)), V.parameters)) || error(
        "Expected type of each vector in params to allow for missing.")
    lenparams = map(length, params)
    all(map(x -> x == first(lenparams), lenparams)) || error(
        "Expected all vectors in params tuple to be of the same length.")
    ParamDistributionVector{D, V}(params)
end

function ParamDistributionVector(dtup::Vararg{Union{Missing,D},N}) where
    {D<:Distribution, N}
    N == 0 && error(
        "Provide at least one distribution in ParamDistributionVector(x...).")
    # need to help compile to determine the type of tupvec
    Ti = typeof(params(first(skipmissing(dtup)))).parameters
    npar = length(Ti)
    tupvec = ntuple(i -> begin
        Tim = Union{Missing, Ti[i]}
        collect(Tim, passmissing(getindex).(passmissing(params).(
            dtup),i))::Vector{Tim}
    end, npar)
    ParamDistributionVector(D, tupvec)
end

function allowmissingtuple(pvec::Vararg{Any,N}) where N
    ntuple(i -> allowmissing(pvec[i]), N) # works only typestable from Julia1.6
end
function ParamDistributionVector(::Type{D}, pvec::Vararg{Any,N}) where
    {D<:Distribution, N}
    isMissingAllowed = all(.<:(Missing, eltype.(pvec)))
    pvecm = isMissingAllowed ? pvec : begin
        inferred_type = first(Base.return_types(allowmissingtuple, typeof.(pvec)))
        inferred_type === Tuple && error(
            "All parameter vectors must allow for missings. Use allowmissing(parm).")
        allowmissingtuple(pvec...)
    end
    ParamDistributionVector(D, pvecm)
end

length(dv::ParamDistributionVector) = length(first(dv.params))
IndexStyle(::Type{<:ParamDistributionVector{D,V}}) where {D,V} =
    IndexStyle(first(V.parameters))

similar(dv::ParamDistributionVector{D,V}) where {D,V} =
    typeof(dv)(ntuple(i->similar(dv.params[i]), length(dv.params)))

function ismissing(dv::ParamDistributionVector, i::Int;
    params_i = getindex.(dv.params, Ref(i)))
    ismissingt = ntuple(i->ismissing(params_i[i]), length(params_i))
    any(ismissingt)::Bool
end
function getindex(dv::ParamDistributionVector, i::Int)
    params_i = getindex.(dv.params, Ref(i))
    ismissing(dv,i; params_i = params_i) && return missing
    nonmissingtype(eltype(dv))(params_i...)
end
function getindex(dv::ParamDistributionVector, I)
    tupvec = map(x -> x[I], dv.params)
    typeof(dv)(tupvec)
    #ParamDistributionVector{D,V}(tupvec)
    #Base.typename(DV).wrapper((dv[i] for i in I)...)
end
function setindex!(dv::ParamDistributionVector, d, i::Int)
    p = ismissing(d) ? ntuple(_->missing,length(dv.params)) : params(d)
    for ti in eachindex(dv.params)
        dv.params[ti][i] = p[ti]
    end
end

function params(dv::ParamDistributionVector, ::Val{i}) where i
    # if types of parameters differ, then a union type is returned -> need Val
    dv.params[i]
end

params(dv::ParamDistributionVector) = dv.params