text:
module PlottingTools

using Colors
using RecipesBase
using LaTeXStrings
using Printf: @sprintf

export TOMPlot,
       Spectrum1D,
       SpectrumAngle1D,
       Spectrum2D,
       EMF2D,
       EMFAngle2D,
       FitSpectrum,
       FitSpectrumEllip,
       SpaceSolution,
       RIprofile,
       PBGDispersion1D,
       PBGDispersion1Dalt,
       PBGDispersion1Dimre,
       PBGDispersion2D,
       PBGDispersion2Dalt

"""

    All the recipes included in this file works with Plots.jl

"""

"""

    Recipe for producing the ThreeOmegaMethod solution plot.

        plot(TOMPlot(), solution; s=(640,480))
        gui()

            solution: structure from the ThreeOmegaMethod
                s: size of the figure

"""
struct TOMPlot end
@recipe function f(::TOMPlot, solution; s=(640,480))
    seriestype  :=  :path
    linestyle --> [:solid :dashdot]
    xlabel --> "Frequency [Hz]"
    ylabel --> "Temperature rise [°C]"
    label --> ["Real" "Imag"]
    xscale --> :log10
    # tickfont --> font(12)
    # legendfont --> font(10)
    size --> s
    solution.f, [solution.ΔThreal, solution.ΔThimag]
end

"""

    Plot the spectrum of Reflectance, Transmittance and Absorbance as input.

        plot(Spectrum1D(), λ, S; s=(640,480))
        gui()

            λ: wavelength range
            S: spectrum
                s: size of the figure

"""
struct Spectrum1D end
@recipe function f(::Spectrum1D, λ, S; s=(640,480))
    seriestype  :=  :path
    linestyle --> :solid
    xlabel --> "Wavelength [nm]"
    # ylabel --> "Reflectance"
    # xlim --> (λ[1], λ[end])
    # ylim --> (0.0, 1.0)
    # tickfont --> font(12)
    # legendfont --> font(10)
    size --> s
    vec(λ), S
end

"""

    Plot the spectrum of Reflectance, Transmittance and Absorbance as input for
    angle of incidence.

        plot(SpectrumAngle1D(), θ, S, s=(640,480))
        gui()

            θ: angle of incidence
            S: spectrum
                s: size of the figure

"""
struct SpectrumAngle1D end
@recipe function f(::SpectrumAngle1D, θ, S; s=(640,480))
    seriestype  :=  :path
    linestyle --> :solid # :dash :dashdot]
    xlabel --> "Angle of incidence [°]"
    # ylabel --> "Reflectance"
    # xlim --> (θ[1], θ[end])
    # ylim --> (0.0, 1.0)
    # tickfont --> font(12)
    # legendfont --> font(10)
    size --> s
    vec(θ), S
end

"""

    Plot the spectra input (Reflectance, Transmittance and Absorbance) for λ and θ.

        plot(Spectrum2D(), λ, θ, S; num_levels=80, s=(640,480))
        gui()

            λ: wavelength range
            θ: angle of incidence
            S: spectrum
                num_levels: number of levels for the contour plot
                s: size of the figure

"""
struct Spectrum2D end
@recipe function f(::Spectrum2D, λ, θ, S; num_levels=80, s=(640,480))
    seriestype  :=  :contour
    fill --> true
    levels --> num_levels
    # color --> colormap("RdBu")
    xlabel --> "Wavelength [nm]"
    ylabel --> "Angle of incidence [°]"
    # xlim --> (λ[1], λ[end])
    # ylim --> (θ[1], θ[end])
    # tickfont --> font(12)
    # legendfont --> font(10)
    size --> s
    vec(λ), vec(θ), Matrix(S')
end

"""

    Plot the EMF for depth and λ, at a given angle of incidence θ.

        plot(EMF2D(), λ, ℓ, emf; num_levels=80, s=(640,480))
        gui()

            λ: wavelength range
            ℓ: multilayer depth
            emf: electromagnetic field
                num_levels: number of levels for the contour plot
                s: size of the figure

"""
struct EMF2D end
@recipe function f(::EMF2D, λ, ℓ, emf; num_levels=80, s=(640,480))
    seriestype  :=  :contour
    fill --> true
    levels --> num_levels
    # color --> colormap("RdBu")
    xlabel --> "Wavelength [nm]"
    ylabel --> "Depth profile [nm]"
    # title --> "EMF Intensity for s-wave"
    # tickfont --> font(12)
    # size --> s
    vec(λ), vec(ℓ), Matrix(emf')
end

"""

    Plot the EMF for depth and θ, at a given λ.

        plot(EMFAngle2D(), θ, ℓ, emf; num_levels=80, s=(640,480))
        gui()

            θ: angle of incidence
            ℓ: multilayer depth
            emf: electromagnetic field
                num_levels: number of levels for the contour plot
                s: size of the figure

"""
struct EMFAngle2D end
@recipe function f(::EMFAngle2D, θ, ℓ, emf; num_levels=80, s=(640,480))
    seriestype  :=  :contour
    fill --> true
    levels --> num_levels
    # color --> colormap("RdBu")
    xlabel --> "Angle of incidence [°]"
    ylabel --> "Depth profile [nm]"
    # title --> "EMF Intensity"
    # tickfont --> font(12)
    # size --> s
    vec(θ), vec(ℓ), Matrix(emf')
end

"""

    Recipe for plotting a comparison of the model and experimental spectra.

        plot(FitSpectrum(), x, Xexp, Xmodel; s=(640,480))
        gui()

            x: range of variable (λ or θ)
            Xexp: experimental spectrum
            Xmodel: model spectrum
                s: size of the figure

"""
struct FitSpectrum end
@recipe function f(::FitSpectrum, x, Xexp, Xmodel; s=(640,480))
    # linecolor --> customcolors[1]
    seriestype := :path
    linestyle --> :solid
    @series begin
        seriestype := :scatter
        markershape --> :circle
        markersize --> 5
        markeralpha --> 0.5
        # markercolor --> customcolors[2]
        markerstrokewidth --> 0.5
        label --> "Data"
        x, Xexp
    end
    label --> "Model"
    # tickfont --> font(12)
    # legendfont --> font(10)
    size --> s
    x, Xmodel
end

"""

    Plot the objective function error obtained from SpaceSolutionRef2D and
    SpaceSolutionTra2D for the space solution given.

        plot(SpaceSolution(),
                x1, x2, S;
                lims=[x1[1], x1[end], x2[1], x2[end]], num_levels=80, s=(640,480)
        )
        gui()

            x1: range of first dimension
            x2: range of second dimension
            S: matrix with the objective function values
            Optional:
                lims: limits of the two axes
                num_levels: number of levels of the contour plot
                s: size of the figure

"""
struct SpaceSolution end
@recipe function f(
    ::SpaceSolution, x1, x2, S;
    lims=[vec(x1)[1], vec(x1)[end], vec(x2)[1], vec(x2)[end]], num_levels=80, s=(640,480),
)
    seriestype := :contour
    fill --> true
    levels --> num_levels
    # color --> colormap("RdBu")
    # tickfont --> font(12)
    size --> s
    xlims --> (lims[1], lims[2])
    ylims --> (lims[3], lims[4])
    vec(x1), vec(x2), Matrix(S')
end

"""

    Recipe for plotting a comparison of the model and experimental spectra of ellipsometry data.

        plot(FitSpectrumEllip(), x, Xexp, Xmodel; s=(640,480))
        gui()

            x: range of variable (λ or θ)
            Xexp: experimental spectrum
            Xmodel: model spectrum
                s: size of the figure

"""
struct FitSpectrumEllip end
@recipe function f(::FitSpectrumEllip, x, Xexp, Xmodel; s=(640,480))
    linecolor --> [RGBA(0.0,0.0,0.0,0.6) RGBA(0.0,0.0,0.0,0.6)]
    seriestype := :path
    linestyle --> :solid
    @series begin
        seriestype := :scatter
        markershape --> [:circle :square]
        markersize --> 5
        markeralpha --> 0.5
        markercolor --> [RGBA(0.9019607843137255, 0.6235294117647059, 0.0, 0.7) RGBA(0.33725490196078434, 0.7058823529411765, 0.9137254901960784, 0.7)]
        markerstrokewidth --> 0.5
        markerstrokecolor --> [RGB(0.9019607843137255, 0.6235294117647059, 0.0) RGB(0.33725490196078434, 0.7058823529411765, 0.9137254901960784)]
        label --> ["ℜ{ρ}" "ℑ{ρ}"]
        x, Xexp
    end
    label --> ["Model" ""]
    size --> s
    x, Xmodel
end

"""

    Plots the index of refraction ath certain wavelength (usually λ0) of the multilayer
    structure.

        plot(RIprofile(),
             solution;
             wave=:b, λ=solution.Misc.λ0, θ=solution.Beam.θ[1], s=(640,480),
        )
        gui()

            solution: structure solution from TMMOptics
                wave = :b (both, default), :p (p-wave), :s (s-wave) of the EMF to overlap
                θ: angle for which to overlap the EMF, by default is taken the first one
                λ: wavelength for which to overlap the EMF, by default is taken the reference one
                s: size of the figure

"""
struct RIprofile end
@recipe function f(
    ::RIprofile, solution;
    plotemf=false, wave=:b, λ=[solution.Misc.λ0], θ=[solution.Beam.θ[1]], s=(640,480),
)
    d = solution.Misc.d[2:end-1] # remove incident and emergent media
    doffset = 0.05*sum(d)
    new_d = [-doffset; cumsum([0; d; doffset], dims=1)]
    w = diff(new_d)
    h = solution.Misc.nλ0
    x = new_d[1:end-1]
    y = zeros(length(x))
    block1 = [x[i] .+ vcat(0.0, w[i], w[i], 0.0) for i=1:length(x)]
    block2 = [y[i] .+ vcat(0.0, 0.0, h[i], h[i]) for i=1:length(x)]
    xlabel --> "Thickness profile [nm]"
    ylabel --> "Refractive index at λ₀"
    @series begin
        seriestype := :shape
        linewidth --> 0.0
        legend --> false
        size --> s
        label --> ""
        # Generate colors for different layers
        color --> Array(colorsUniqueLayers(solution.Misc.d)')
        alpha --> 0.6
        block1, block2
    end
    if plotemf # overlap EMF
        if ~isempty(solution.Field.emfp) # Check is not empty
            aux1 = findmin(abs.(solution.Beam.λ .- λ[1]))[2][1]
            aux2 = findmin(abs.(solution.Beam.θ .- θ[1]))[2][1]
            λaux1 = solution.Beam.λ[aux1]
            θaux2 = solution.Beam.θ[aux2]
            if isequal(wave, :b)
                @series begin
                    seriestype := :path
                    linestyle --> :solid
                    label --> "EMF-p"
                    solution.Misc.ℓ, solution.Field.emfp[aux1, aux2, :]
                end
                @series begin
                    seriestype := :path
                    linestyle --> :dash
                    label --> "EMF-s"
                    solution.Misc.ℓ, solution.Field.emfs[aux1, aux2, :]
                end
            elseif isequal(wave, :p)
                @series begin
                    seriestype := :path
                    linestyle --> :solid
                    label --> "EMF-p"
                    solution.Misc.ℓ, solution.Field.emfp[aux1, aux2, :]
                end
            elseif isequal(wave, :s)
                @series begin
                    seriestype := :path
                    linestyle --> :solid
                    label --> "EMF-s"
                    solution.Misc.ℓ, solution.Field.emfs[aux1, aux2, :]
                end
            end # if isequal(wave, :b)
        end # if ~isempty(solution.Field.emfp) # Check is not empty
    end # if plotemf
end

"""

    Check the unique layers in the system.

"""
function colorsUniqueLayers(d)
    # Generate colors for different layers, using the thickness vector
    unique_layers = unique(d)
    # cols = distinguishable_colors((num_layers+1), [RGB(1,1,1)])[2:end]
    cols0 = distinguishable_colors(length(unique_layers)*2, [RGB(1,1,1)])[2:end]
    cols = cols0[1:2:end]
    # pcols = map(col -> (red(col), green(col), blue(col)), cols)
    # Match one color for one layer in the structure
    assigned_cols = Array{RGB}(undef, length(d))
    for i ∈ eachindex(cols)
        assigned_cols[unique_layers[i] .== d] .= cols[i]
    end
    # # reassign incident and substrate to differ from possible same layers
    assigned_cols[1] = cols0[1]
    assigned_cols[end] = cols0[end-1]
    return assigned_cols
end

"""

    Plots the photonic dispersion (Bloch wavevector) of a photonic crystal structure,
    computed for a range of frequencies (wavelengths) and one angle of incidence.

    Takes both polarisation types and renders them into one piece. You need to pass either
    the real or imaginary parts of them.

        plot(PBGDispersion1Dalt(), Bloch; kpart=:real, s=(640,480))
        gui()

            solution.Bloch: Bloch structure of the solution
            kpart: either real (:real, default) or imginary (:imag)
                s: size of the figure

"""
struct PBGDispersion1Dalt end
@recipe function f(::PBGDispersion1Dalt, Bloch; kpart=:real, s=(640,480))
    ω = Bloch.ω .* Bloch.Λ ./ 2.0 ./ π # frequency range normalized
    xlabel --> "KΛ/π"
    ylabel --> "ωΛ/(2π)"
    # title --> "p/TM     s/TE"
    xticks --> (-1:0.5:1, string.([1.0, 0.5, 0.0, 0.5, 1.0]))
    yticks --> ([0.], [" "])
    size --> s
    framestyle --> :zerolines
    @series begin
        seriestype := :path
        linestyle --> :solid
        label --> "p/TM"
        κ = Bloch.κp .* Bloch.Λ ./ π
        isequal(kpart, :real) ? -real.(κ) : -imag.(κ), vec(ω)
    end
    @series begin
        seriestype := :path
        linestyle --> :dashdot
        label --> "s/TE"
        κ = Bloch.κs .* Bloch.Λ ./ π
        isequal(kpart, :real) ? real.(κ) : imag.(κ), vec(ω)
    end
end

"""

    Plots the photonic dispersion (Bloch wavevector) of a photonic crystal structure,
    computed for a range of frequencies (wavelengths) and one angle of incidence.

    Takes one polarisation type in complex format, and plots on the left the imaginary
    part and on the right the real part.

        plot(PBGDispersion1Dimre(), Bloch; s=(640,480))
        gui()

            solution.Bloch: Bloch structure of the solution
            wave: either p-wave (:p, default) or s-wave (:s)
                s: size of the figure

"""
struct PBGDispersion1Dimre end
@recipe function f(::PBGDispersion1Dimre, Bloch; wave=:p, s=(640,480))
    ω = Bloch.ω .* Bloch.Λ ./ 2.0 ./ π # frequency range normalized
    xlabel --> "KΛ/π"
    ylabel --> "ωΛ/(2π)"
    xticks --> (-1:0.5:1, string.([1.0, 0.5, 0.0, 0.5, 1.0]))
    yticks --> ([0.], [" "])
    size --> s
    framestyle --> :zerolines
    if isequal(wave, :p)
        κ = Bloch.κp .* Bloch.Λ ./ π
    elseif isequal(wave, :s)
    else
        throw("The wave parameter must be equal to :p or :s.")
    end
    @series begin
        seriestype := :path
        linestyle --> :solid
        label --> "Imag"
        -abs.(imag.(κ)), vec(ω)
    end
    @series begin
        seriestype := :path
        linestyle --> :dashdot
        label --> "Real"
        real.(κ), vec(ω)
    end
end

"""

    Plots the photonic dispersion (Bloch wavevector) of a photonic crystal structure,
    computed for a range of frequencies (wavelengths) and one angle of incidence.

    Takes only one polarisation. If you want you can pass optionally the part of the
    wavevector to plot. By default, plots the real one (as the imaginary does not represent
    much).

        plot(PBGDispersion1D(), solution.Bloch; wave=:p, kpart=:real, s=(640,480))
        gui()

            solution.Bloch: Bloch structure of the solution
                wave: either p-wave (:p, default) or s-wave (:s)
                kpart: either real (:real, default) or imginary (:imag)
                s: size of the figure

"""
struct PBGDispersion1D end
@recipe function f(::PBGDispersion1D, Bloch; kpart=:real, wave=:p, s=(640,480))
    ω = Bloch.ω .* Bloch.Λ ./ 2.0 ./ π # frequency range normalized
    if isequal(wave, :p)
        κ = Bloch.κp .* Bloch.Λ ./ π
    elseif isequal(wave, :s)
        κ = Bloch.κs .* Bloch.Λ ./ π
    else
        throw("The wave parameter should be either :p or :s.")
    end
    seriestype := :path
    linestyle --> :solid
    xlabel --> "KΛ/π"
    ylabel --> "ωΛ/(2π)"
    label --> ""
    size --> s
    isequal(kpart, :real) ? real.(κ) : imag.(κ), ω
end

"""

    Plots the photonic dispersion (Bloch wavevector) of a photonic crystal structure,
    computed for a range of frequencies (wavelengths) and a range of angle of incidences.

        plot(PBGDispersion2Dalt(), solution.Bloch; s=(640,480), num_levels=90)
        gui()

            solution: solution structure from TMMOptics
                s: size of the figure

"""
struct PBGDispersion2Dalt end
@recipe function f(::PBGDispersion2Dalt, Bloch; s=(640,480), num_levels=90)
    qz = Bloch.qz; ωh = Bloch.ωh; ωl = Bloch.ωl
    ω = Bloch.ω .* Bloch.Λ ./ 2.0 ./ π # frequency range normalized
    # k normalized for angle-frequency dependence
    # logical matrices, used to select points which belong to the forbidden bands
    κp = cos.(Bloch.κp .* Bloch.Λ) # rhs equation 45 in the paper, π comes from previous normalization
    κs = cos.(Bloch.κs .* Bloch.Λ) # rhs equation 45 in the paper
    κp[abs.(κp) .< 1.0] .= 1.0 # propagating p wave
    κp[abs.(κp) .> 1.0] .= 0.0 # evanescent p wave
    κs[abs.(κs) .< 1.0] .= 1.0 # propagating s wave
    κs[abs.(κs) .> 1.0] .= 0.0 # evanescent s wave
    size --> s
    xlabel --> "Parallel wavevector, qz (2π/Λ)"
    ylabel --> "ωΛ/(2π)"
    title --> "p/TM-wave                        s/TE-wave"
    grid --> true
    # legend --> false
    gridstyle --> :dot
    gridlinewidth --> 1.5
    xlims --> (-qz[end], qz[end])
    ylims --> (0.0, ω[1])
    @series begin
        seriestype := :contour
        color --> colormap("Blues")
        # linewidth --> 2.0
        colorbar_entry --> false
        fill --> true
        levels --> num_levels
        # alpha --> 0.3
        label --> ""
        -vec(qz), vec(ω), real.(κp)
    end
    @series begin
        seriestype := :contour
        color --> colormap("Reds")
        # linewidth --> 2.0
        colorbar_entry --> false
        fill --> true
        levels --> num_levels
        # alpha --> 0.3
        label --> ""
        vec(qz), vec(ω), real.(κs)
    end
    @series begin
        seriestype := :path
        linestyle --> :solid
        linewidth --> 1.5
        linecolor --> RGBA(0.0, 0.0, 0.0, 0.7)
        label --> "Light line"
        [0.0, ω[1]], [0.0, ω[1]]
    end
    @series begin
        seriestype := :path
        linestyle --> :solid
        linewidth --> 1.5
        linecolor --> RGBA(0.0, 0.0, 0.0, 0.7)
        label --> ""
        [0.0, -ω[1]], [0.0, ω[1]]
    end
    @series begin
        seriestype := :path
        linestyle --> :dash
        linewidth --> 1.5
        linecolor --> RGBA(0.0, 0.0, 0.0, 0.7)
        label --> "ωₕ"
        [-qz[end], qz[end]], [ωh, ωh]
    end
    @series begin
        seriestype := :path
        linestyle --> :dashdot
        linewidth --> 1.5
        linecolor --> RGBA(0.0, 0.0, 0.0, 0.7)
        label --> "ωₗ"
        [-qz[end], qz[end]], [ωl, ωl]
    end
end

"""

    Plots the photonic dispersion (Bloch wavevector) of a photonic crystal structure,
    computed for a range of frequencies (wavelengths) and a range of angle of incidences.

    This function plots the Bloch wavevector for only one of the polarisation types.

        plot(PBGDispersion2D(), solution.Bloch; wave=:p, s=(640,480), num_levels=90)
        gui()

            solution.Bloch: solution.Bloch structure from TMMOptics
                wave: either :p (p-wave, default) or :s (s-wave)
                s: size of the figure

"""
struct PBGDispersion2D end
@recipe function f(::PBGDispersion2D, Bloch; wave=:p, s=(640,480), num_levels=90)
    qz = Bloch.qz; ωh = Bloch.ωh; ωl = Bloch.ωl
    Λ = Bloch.Λ # periodicity
    ω = Bloch.ω .* Λ ./ 2.0 ./ π # frequency range normalized
    if isequal(wave, :p)
        kbloch = Bloch.κp .* Λ# ./ π
        cmap = colormap("Blues")
    elseif isequal(wave, :s)
        kbloch = Bloch.κs .* Λ# ./ π
        cmap = colormap("Reds")
    else
        throw("The wave type must be either :p or :s.")
    end
    # k normalized for angle-frequency dependence
    # logical matrices, used to select points which belong to the forbidden bands
    κ = cos.(kbloch)#.*π) # rhs equation 45 in the paper, π comes from previous normalization
    κ[abs.(κ) .< 1.0] .= 1.0 # propagating wave
    κ[abs.(κ) .> 1.0] .= 0.0 # evanescent wave
    size --> s
    xlabel --> "Parallel wavevector, qz (2π/Λ)"
    ylabel --> "ωΛ/(2π)"
    grid --> true
    # legend --> false
    gridstyle --> :dot
    gridlinewidth --> 1.5
    xlims --> (0.0, qz[end])
    ylims --> (0.0, ω[1])
    @series begin
        seriestype := :contour
        fill --> true
        levels --> num_levels
        colorbar_entry --> false
        color --> cmap
        # alpha --> 0.9
        label --> ""
        vec(qz), vec(ω), real.(κ)
    end
    @series begin
        seriestype := :path
        linestyle --> :solid
        linewidth --> 1.5
        linecolor --> RGBA(0.0, 0.0, 0.0, 0.7)
        label --> "Light line"
        [0.0, ω[1]], [0.0, ω[1]]
    end
    @series begin
        seriestype := :path
        linestyle --> :dash
        linewidth --> 1.5
        linecolor --> RGBA(0.0, 0.0, 0.0, 0.7)
        label --> "ωₕ"
        [0.0, qz[end]], [ωh, ωh]
    end
    @series begin
        seriestype := :path
        linestyle --> :dashdot
        linewidth --> 1.0
        linecolor --> RGBA(0.0, 0.0, 0.0, 0.7)
        label --> "ωₗ"
        [0.0, qz[end]], [ωl, ωl]
    end
end

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule AMGX_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("AMGX")
JLLWrappers.@generate_main_file("AMGX", UUID("656d14af-56e4-5275-8e68-4e861d7b5043"))
end  # module AMGX_jll

ext: jl


text:
import FileIO
import CSV
import StatsBase
import AbstractPlotting
import Makie
import DataFrames

using Porta

α = 40 / 180 * pi
solidtop = U1(pi)
solidbottom = U1(-pi + 2α)
ghosttop = U1(-pi + 2α)
ghostbottom = U1(-pi)

"""
sample(dataframe, part, max)

Samples points from a dataframe with the given `dataframe`, `part` id, and the `max` number
of samples limit.
"""
function sample(dataframe, part, max)
    groupdataframe = DataFrames.groupby(dataframe, :partid)
    subdataframe = groupdataframe[(partid=part,)]
    ϕ = subdataframe.x ./ 180 .* pi
    θ = subdataframe.y ./ 180 .* pi
    coordinates = map(x -> Geographic(x[1], x[2]), eachrow([ϕ θ]))
    sampledpoints = Array{Geographic,1}(undef, max)
    count = length(coordinates)
    if count > max
        StatsBase.sample!(coordinates, sampledpoints, replace = false, ordered = true)
    else
        sampledpoints = coordinates
    end
    sampledpoints
end

# Made with Natural Earth.
# Free vector and raster map data @ naturalearthdata.com.
countries = Dict("iran" => [0.800, 0.867, 0.090],
                 "us" => [0.788, 0.200, 0.067],
                 "china" => [0.882, 0.882, 0.431],
                 "ukraine" => [0.776, 0.447, 0.788],
                 "australia" => [0.827, 0.145, 0.510],
                 "germany" => [0.082, 0.784, 0.573],
                 "israel" => [0.522, 0.333, 0.875],
                 "canada" => [0.212, 0.933, 0.369],
                 "india" => [0.525, 0.902, 0.584],
                 "southkorea" => [0.180, 0.827, 0.784],
                 "france" => [0.518, 0.812, 0.184],
                 "antarctica" => [0.875, 0.651, 0.357])
# The path to the dataset
path = "test/data/natural_earth_vector"
# The scene object that contains other visual objects
scene = Makie.Scene(backgroundcolor = :white, show_axis = false, resolution = (360, 360))
maxsamples = 360
segments = 36
s3rotation = Quaternion(0, ℝ³(0, 0, 1))
config = Biquaternion(ℝ³(0, 0, 0))
solidwhirls = []
ghostwhirls = []
parts = 3
for country in countries
    countryname = country[1]
    dataframe = DataFrames.DataFrame(CSV.File(joinpath(path, "$(countryname)-nodes.csv")))
    solidcolor = Makie.RGBAf0(country[2]..., 0.9)
    ghostcolor = Makie.RGBAf0(country[2]..., 0.1)
    part = 0
    if countryname == "antarctica"
        points = sample(dataframe, part, maxsamples)
        ghostwhirl = Whirl(scene,
                           points,
                           top = ghosttop,
                           bottom = ghostbottom,
                           s3rotation = s3rotation,
                           config = config,
                           segments = segments,
                           color = ghostcolor,
                           transparency = true)
        push!(ghostwhirls, ghostwhirl)
    elseif countryname in ["us", "france", "iran", "canada", "china"]
        for part in 0:parts-1
            # Sample a random subset of the points
            points = sample(dataframe, part, maxsamples)
            solidwhirl = Whirl(scene,
                               points,
                               top = solidtop,
                               bottom = solidbottom,
                               s3rotation = s3rotation,
                               config = config,
                               segments = segments,
                               color = solidcolor,
                               transparency = false)
            push!(solidwhirls, solidwhirl)
            ghostwhirl = Whirl(scene,
                               points,
                               top = ghosttop,
                               bottom = ghostbottom,
                               s3rotation = s3rotation,
                               config = config,
                               segments = segments,
                               color = ghostcolor,
                               transparency = true)
            push!(ghostwhirls, ghostwhirl)
        end
    else
        # Sample a random subset of the points
        points = sample(dataframe, part, maxsamples)
        solidwhirl = Whirl(scene,
                           points,
                           top = solidtop,
                           bottom = solidbottom,
                           s3rotation = s3rotation,
                           config = config,
                           segments = segments,
                           color = solidcolor,
                           transparency = false)
        push!(solidwhirls, solidwhirl)
    end
end
# Use QGIS to design a geo map
framecolor = FileIO.load("test/data/basemap90grid.png")
framesprites = []
framesprite1 = Frame(scene,
                     solidtop,
                     framecolor,
                     s3rotation = s3rotation,
                     config = config,
                     segments = segments,
                     transparency = false)
framesprite2 = Frame(scene,
                     solidbottom,
                     framecolor,
                     s3rotation = s3rotation,
                     config = config,
                     segments = segments,
                     transparency = false)
push!(framesprites, framesprite1)
push!(framesprites, framesprite2)
frames = 360

"""
    animate(i)

Update the state of observables with the given frame number `i`.
"""
function animate(i)
    step = (i - 1) / frames
    println("Step: ", 100step)
    τ = step * 2pi - pi
    q = Quaternion(τ, ℝ³(0, 0, 1))
    for whirl in solidwhirls
        update(whirl, q)
    end
    for whirl in ghostwhirls
        update(whirl, q)
    end
    update(framesprite1, q)
    update(framesprite2, q)
end

# update eye position
# scene.camera.eyeposition.val
upvector = Makie.Vec3f0(0, 0, 1)
eyeposition = Makie.Vec3f0(sqrt(3) / 3, -sqrt(3) / 3, sqrt(3) / 3) .* 4
lookat = Makie.Vec3f0(sqrt(3) / 3, sqrt(3) / 3, sqrt(3) / 3)
Makie.update_cam!(scene, eyeposition, lookat, upvector)
scene.center = false # prevent scene from recentering on display
# Makie.save("gallery/drorbarnatan2010.jpg", scene)
Makie.record(scene, "gallery/drorbarnatan2010.gif") do io
    for i in 1:frames
        animate(i)
        Makie.recordframe!(io) # record a new frame
    end
end


text:
import RDates: WeekendCalendar, @rd_str, setcachedcalendar!
import Dates: Date

setcachedcalendar!(CALENDARS, "TARGET", WeekendCalendar() + RuleBasedCalendar([
    # New Years Day
    PeriodicCalendarRule(rd"1y", rd"1JAN"),
    # Labour Day
    WindowCalendarRule(PeriodicCalendarRule(rd"1y", rd"1MAY"); from=Date(2000,1,1)),
    # Good Friday
    WindowCalendarRule(PeriodicCalendarRule(rd"1y", rd"0E-1FRI"); from=Date(2000,1,1)),
    # Easter Monday
    WindowCalendarRule(PeriodicCalendarRule(rd"1y", rd"0E+1MON"); from=Date(2000,1,1)),
    # Christmas Day
    PeriodicCalendarRule(rd"1y", rd"25DEC"),
    # Day of Goodwill
    WindowCalendarRule(PeriodicCalendarRule(rd"1y", rd"26DEC"); from=Date(2000,1,1)),
    # End of Year
    ExplicitDateRule(inclusions=[Date(1998,12,31), Date(1999,12,31), Date(2001,12,31)]),
]))


text:
import Backtest
import Test: @test, @testset, @test_throws
import NamedArrays: NamedArray
import Statistics: mean, std
import DataFrames: DataFrame
import Dates: Minute, Day, Date, DateTime
import CSV
import Random: MersenneTwister, randperm

# Type definitions (defined here to avoid local-type-definition error)
struct  UndefWindowFieldOp <: Backtest.AbstractFields.AbstractWindowFieldOperation end
struct UndefCSFieldOp <: Backtest.AbstractFields.AbstractCrossSectionalFieldOperation end

@testset "Backtest.jl" begin

    # Tests the AbstractFields module
    @testset "AbstractFields" begin
        # Test that exceptions are thrown when calling dofieldop on abstract types
        windowfieldop = UndefWindowFieldOp()
        crossectionalfieldop = UndefCSFieldOp()

        @test_throws Exception Backtest.AbstractFields.dofieldop(windowfieldop, [])
        @test_throws Exception Backtest.AbstractFields.dofieldop(crossectionalfieldop, Dict{String, Any}())
    end

    # Tests the ConcreteFields module
    @testset "ConcreteFields" begin
        # Test default constructor for OHLCV
        @test Backtest.ConcreteFields.Open() == Backtest.ConcreteFields.Open("Open")
        @test Backtest.ConcreteFields.High() == Backtest.ConcreteFields.High("High")
        @test Backtest.ConcreteFields.Low() == Backtest.ConcreteFields.Low("Low")
        @test Backtest.ConcreteFields.Close() == Backtest.ConcreteFields.Close("Close")
        @test Backtest.ConcreteFields.Volume() == Backtest.ConcreteFields.Volume("Volume")

        # Test default window field operation constructors
        @test Backtest.ConcreteFields.Returns("", "", 2) == Backtest.ConcreteFields.Returns("", "")
        @test Backtest.ConcreteFields.LogReturns("", "", 2) == Backtest.ConcreteFields.LogReturns("", "")

        function windows(vec, n)
            return [windowdata[i:i+n] for i = 1:length(vec)-n]
        end

        # Data to be used for window field operations
        windowdata = [i*1.0 for i = 1:15]
        n_to_windows = Dict{Integer, Vector}()
        n_to_windows[2]=windows(windowdata, 2)
        n_to_windows[3]=windows(windowdata, 3)
        n_to_windows[10]=windows(windowdata, 10)

        makeexpectedreturns(n) = (x -> (x[n]-x[1]) / x[1])
        expectedreturns = makeexpectedreturns(2)

        makeexpectedlogreturns(n) = (x -> log(x[n]/x[1]))
        expectedlogreturns = makeexpectedlogreturns(2)

        makeexpectedsma(n) = (x->sum(x)/length(x))

        for n in keys(n_to_windows)
            # make field operations
            returnsop = Backtest.ConcreteFields.Returns("", "", n)
            logreturnsop =Backtest.ConcreteFields.LogReturns("", "", n)
            smaop = Backtest.ConcreteFields.SMA("", "", n)

            # test dofieldop
            windows = n_to_windows[n]
            for window in windows
                @test Backtest.ConcreteFields.dofieldop(returnsop, window) == makeexpectedreturns(n)(window)
                @test Backtest.ConcreteFields.dofieldop(logreturnsop, window) == makeexpectedlogreturns(n)(window)
                @test Backtest.ConcreteFields.dofieldop(smaop, window) == makeexpectedsma(n)(window)
            end
        end

        # Create data to be used for cross sectional field operations
        names = ["A", "B", "C", "D"]
        csdata = [
            NamedArray([1.0, 2, 3, 4], (names,) ),
            NamedArray([3, 1, 40, -2], (names,))
        ]

        # Test Z-Score
        means = mean(csdata)
        means = (x -> sum(x)/length(x)).(csdata)
        stddevs = std.(csdata)
        function makezscorefn(mu, sigma)
            function zscorefn(x)
                return (x-mu)/sigma
            end
            return zscorefn
        end
        zscorefns = [makezscorefn(mean, stddev) for (mean, stddev) in zip(means, stddevs)]

        zsc = Backtest.ConcreteFields.ZScore("", "")
        for (i, csdatum) in enumerate(csdata)
            expectedresults = zscorefns[i].(csdatum)
            actualresults = Backtest.ConcreteFields.dofieldop(zsc, csdatum)
            @test expectedresults == actualresults
        end

        # Test rank
        expectedranks = [
            NamedArray([4, 3, 2, 1], (names,)),
            NamedArray([2, 3, 1, 4], (names,))
        ]
        makerankdofieldop() = (datum -> Backtest.ConcreteFields.dofieldop(Backtest.ConcreteFields.Rank("", ""), datum))
        actualranks = makerankdofieldop().(csdata)
        @test expectedranks == actualranks
    end

    # Tests the Engine module
    @testset "Engine" begin
        nbarsstoredoptions =
        assetids = ["A", "B", "C"]
        (O, H, L, C, V) = ("Open", "High", "Low", "Close", "Volume")
        function testengine(nbarstostore)
            lattice = Backtest.Engine.CalcLattice(nbarstostore, assetids)
            @test Backtest.Engine.numbarsavailable(lattice) == 0
            @test_throws Exception Backtest.Engine.data(lattice)

            # Make field operations and add them to the backtest
            fieldoperations = [
                Backtest.ConcreteFields.Open(O),
                Backtest.ConcreteFields.High(H),
                Backtest.ConcreteFields.Low(L),
                Backtest.ConcreteFields.Close(C),
                Backtest.ConcreteFields.Volume(V),
                Backtest.ConcreteFields.SMA("SMA1High", H, 1),
                Backtest.ConcreteFields.SMA("SMA2Open", O, 2),
                Backtest.ConcreteFields.Rank("RankLow", L),
                Backtest.ConcreteFields.Rank("RankSMA1High", "SMA1High")
            ]
            Backtest.Engine.addfields!(lattice, fieldoperations)

            # Make bars of data and run them
            allbars = [
                Dict(
                    "A" => Dict(O=>10, H=>15, L=>8, C=>11, V=>10000),
                    "B" => Dict(O=>100, H=>101, L=>90, C=>93, V=>101),
                    "C" => Dict(O=>60, H=>80, L=>60, C=>80, V=>10000),
                ),
                Dict(
                    "A" => Dict(O=>11, H=>11, L=>3, C=>6, V=>8000),
                    "B" => Dict(O=>93, H=>100, L=>90, C=>99, V=>101),
                    "C" => Dict(O=>80, H=>80, L=>60, C=>80, V=>10000),
                )
            ]

            # Test that the first bar of data gets loaded properly and
            # computes dependent fields
            Backtest.Engine.newbar!(lattice, allbars[1])
            @test Backtest.Engine.data(lattice, "A", O) == 10
            @test Backtest.Engine.data(lattice, "B", "SMA1High") == 101
            @test Backtest.Engine.data(lattice, "B", "RankLow") == 1
            @test Backtest.Engine.data(lattice, "C", "RankLow") == 2
            @test Backtest.Engine.data(lattice, "B", "RankSMA1High") == 1

            @test Backtest.Engine.numbarsavailable(lattice) == 1

            # Test that the second bar of data gets load
            Backtest.Engine.newbar!(lattice, allbars[2])
            Backtest.Engine.data(lattice, "A", O)
            @test Backtest.Engine.data(lattice, "A", O) == 11
            @test Backtest.Engine.data(lattice, "B", "SMA1High") == 100
            @test Backtest.Engine.data(lattice, "B", "RankLow") == 1
            @test Backtest.Engine.data(lattice, "C", "RankLow") == 2
            @test Backtest.Engine.data(lattice, "B", "RankSMA1High") == 1
            @test Backtest.Engine.data(lattice, "A", "SMA2Open") == (10+11)/2
            @test Backtest.Engine.data(lattice, "B", "SMA2Open") == (93+100)/2
            @test Backtest.Engine.data(lattice, "C", "SMA2Open") == (60+80)/2

            @test Backtest.Engine.numbarsavailable(lattice) == 2
        end

        testengine(-1)
        testengine(5)
        testengine(1_000_000_000)
    end

    # Tests the DataReaders module
    @testset "DataReaders" begin
        # Generate some data that spans 10 files with 1000 entries each
        nfiles = 10
        barsperfile = 5

        basedate = DateTime(1)
        dates = [basedate + Day(i) for i = 1:nfiles]
        datadir = "____tempdata___"
        mkdir(datadir)
        try
            filenames = ["$datadir/$(Date(date)).csv" for date in dates]

            (openbase, highbase, lowbase, closebase, volumebase) = (100, 120, 79.3, 98.1, 10000)
            oval = ((i, j) -> i+j)
            hval = ((i, j) -> i+j+10)
            lval = ((i, j) -> i+j-20)
            cval = ((i, j) -> i+j)
            vval = ((i, j) -> 100)
            for (i, filename) in enumerate(filenames)
                df = DataFrame(
                        DateTime = [dates[i] + Minute(j-1) for j in 1:barsperfile],
                        Open = [oval(i, j) for j in 1:barsperfile],
                        High = [hval(i, j) for j in 1:barsperfile],
                        Low = [lval(i, j) for j in 1:barsperfile],
                        Close = [cval(i, j) for j in 1:barsperfile],
                        Volume = [vval(i, j) for j in 1:barsperfile]
                    )
                CSV.write(filename, df)
            end

            function makedatareaders()
                # Make all of the datareaders
                return [
                    Backtest.DataReaders.InMemoryDataReader(
                        "",
                        filenames,
                        datetimecol="DateTime",
                        dtfmt="yyyy-mm-ddTHH:MM:SS"),
                    Backtest.DataReaders.PerFileDataReader(
                        "",
                        filenames,
                        datetimecol="DateTime",
                        dtfmt="yyyy-mm-ddTHH:MM:SS"
                    )
                ]
            end

            # end

            function runthroughdatareader(datareader)
                for i in 1:nfiles
                    for j in 1:barsperfile
                        peekval =  Backtest.DataReaders.peek(datareader)
                        popval =  Backtest.DataReaders.popfirst!(datareader)
                        @assert peekval == popval
                        @assert peekval["Open"] == oval(i, j)
                        @assert peekval["High"] == hval(i, j)
                        @assert peekval["Low"] == lval(i, j)
                        @assert peekval["Close"] == cval(i, j)
                        @assert peekval["Volume"] == vval(i, j)
                    end
                end
                return true
            end

            for datareader in makedatareaders()
                # Test that that fastforward throws an exception when the given date is too early
                @test_throws Backtest.Exceptions.DateTooEarlyError Backtest.DataReaders.fastforward!(datareader, DateTime(0))
                @test_throws Backtest.Exceptions.DateTooEarlyError Backtest.DataReaders.fastforward!(datareader, Date(0, 1, 1))

                # Test that fastforward makes the underlying DataReader's current time
                # be greater than the input time.
                @test Backtest.DataReaders.peek(datareader)["DateTime"] > basedate

                # Test that peek gives the expected result
                @test runthroughdatareader(datareader)

            end

            for datareader in makedatareaders()
                # Test that that fastforward throws an exception when the given date is too late
                @test_throws Backtest.Exceptions.DateTooFarOutError Backtest.DataReaders.fastforward!(datareader, DateTime(2030))
            end
        finally
            rm(datadir, recursive=true) # clean up
        end
    end

    # Tests the Events module
    @testset "Events" begin
        exampleorder = Backtest.Orders.MarketOrder("A", 10)
        sortedevents = [
            Backtest.Events.NewBarEvent(DateTime(0), Dict(""=>Dict(""=>nothing))),
            Backtest.Events.FieldCompletedProcessingEvent(DateTime(1)),
            Backtest.Events.OrderAckEvent(DateTime(1, 2), "__an_order_id"),
            Backtest.Events.OrderFillEvent(DateTime(1, 3), exampleorder, 100, -1),

            Backtest.Events.NewBarEvent(DateTime(2), Dict(""=>Dict(""=>nothing))),
            Backtest.Events.FieldCompletedProcessingEvent(DateTime(3)),
            Backtest.Events.OrderAckEvent(DateTime(3, 2), "__an_order_id"),
            Backtest.Events.OrderFillEvent(DateTime(3, 3), exampleorder, 100, -1),

            Backtest.Events.NewBarEvent(DateTime(4), Dict(""=>Dict(""=>nothing))),
            Backtest.Events.FieldCompletedProcessingEvent(DateTime(5)),
            Backtest.Events.OrderAckEvent(DateTime(5, 2), "__an_order_id"),
            Backtest.Events.OrderFillEvent(DateTime(5, 3), exampleorder, 100, -1)
        ]

        seed = MersenneTwister(1234)
        indexpermutations = randperm(seed, length(sortedevents))

        eventq = Backtest.Events.EventQueue()
        @test Backtest.Events.empty(eventq)

        for index in indexpermutations
            Backtest.Events.push!(eventq, sortedevents[index])
            @test !Backtest.Events.empty(eventq)
        end

        for event in sortedevents
            @test event == Backtest.Events.peek(eventq)
            @test event == Backtest.Events.pop!(eventq)
        end

        @test Backtest.Events.empty(eventq)
    end
    #
    # Tests the Backtest.jl API
    @testset "API" begin
        # TODO: Write tests for the backtest API
    end

    @testset "Utils" begin

    end
end

max_stars_repo_name: maxholloway/backtest

max_issues_repo_name: maxholloway/backtest

max_forks_repo_name: maxholloway/backtest


text:
############
## TlsIID ##
############

type TlsIID <: AbstrUnivarModel
    ν::Float64
    μ::Float64
    a::Float64
    function TlsIID(ν::Float64, μ::Float64, a::Float64)
        if a <= 0
            error("scaling parameter a must be larger than 0")
        end
        if ν <= 0
            error("ν parameter must be larger than 0")
        end
        if ν <= 2
            warn("ν parameter is smaller than 2")
        end
        return new(ν, μ, a)
    end
end

function TlsIID(params::Array{Float64, 1})
    return TlsIID(params...)
end

import Base.Multimedia.display
function display(mod::TlsIID)
    nDigits = 3
    println("\nIID, constant univariate t-location scale distribution:")
    println("    ν: $(round(mod.ν, nDigits))")
    println("    μ: $(round(mod.μ, nDigits))")
    println("    a: $(round(mod.a, nDigits))")
end

function description(mod::TlsIID)
    return "IID, Student's t location scale distributed"
end

## getParams
##----------

function getParams(mod::TlsIID)
    return (mod.ν, mod.μ, mod.a)
end

## simulate
##---------

function simulate(mod::TlsIID, nObs::Int)
    ## simulate from TlsIID model
    t = TDist(mod.ν)
    return rand(t, nObs)*mod.a + mod.μ
end

## estimate
##---------

function estimate(dt::Type{TlsIID}, data::Array{Float64, 1})
    tlsFit = fit(TLSDist, data)
    return TlsIID(getParams(tlsFit)...)
end

## resimulate
##-----------

@doc doc"""
Fit model and resimulate data from it.
"""->
function resimulate(data::AbstractTimenum,
                    dt::Type{Econometrics.TlsIID},
                    nPaths::Int = 1)
    ## get returns without NAs
    pureRets = dropna(data.vals[1])
    nRets = size(pureRets, 1)

    ## fit model to pure returns
    modFit = estimate(dt, pureRets)

    return resimulate(data, modFit, nPaths)
end


text:
"""
Return a new processing key with the number incremented.
It checks for existing keys and returns a string with the next key to be used.

#### Arguments

* `d`: Dictionary containing existing keys
* `key_name`: Base of the

#### Returns

* AbstractString with new key name

#### Returns

```julia
results_storage = Dict()
results_storage[new_processing_key(results_storage, "FTest")] = 4
results_storage[new_processing_key(results_storage, "FTest")] = 49

# Dict(Any, Any) with 2 entries
#   "FTest1" => 4
#   "FTest2" => 49
```
"""
function new_processing_key(d::Dict, key_name::AbstractString)
    key_numb = 1
    key = string(key_name, key_numb)
    while haskey(d, key)
        key_numb += 1
        key = string(key_name, key_numb)
    end
    return key
end

"""
Find dictionary keys containing a string.

#### Arguments

* `d`: Dictionary containing existing keys
* `partial_key`: AbstractString you want to find in key names

#### Returns

* Array containg the indices of dictionary containing the partial_key

#### Returns

```julia
results_storage = Dict()
results_storage[new_processing_key(results_storage, "FTest")] = 4
results_storage[new_processing_key(results_storage, "Turtle")] = 5
results_storage[new_processing_key(results_storage, "FTest")] = 49

find_keys_containing(results_storage, "FTest")

# 2-element Array{Int64,1}:
#  1
#  3
```
"""
function find_keys_containing(d::Dict, partial_key::AbstractString)
    valid_keys = [startswith(i, partial_key) for i in collect(keys(d))]
    findall((in)(true), valid_keys)
end

"""
Extract the path, filename and extension of a file

#### Arguments

* `fname`: AbstractString with the full path to a file

#### Output

* Three strings containing the path, file name and file extension

#### Returns

```julia
fileparts("/Users/test/subdir/test-file.bdf")

# ("/Users/test/subdir/","test-file","bdf")
```
"""
function fileparts(fname::AbstractString)
    if fname == ""
        pathname = ""
        filename = ""
        extension = ""
    else

        pathname = dirname(fname)
        if pathname == ""
            #nothing
        else
            pathname = string(pathname, "/")
        end
        filename = splitext(basename(fname))[1]
        extension = splitext(basename(fname))[2][2:end]
    end

    return pathname, filename, extension
end

"""
Find the closest number to a target in an array and return the index

#### Arguments

* `list`: Array containing numbers
* `target`: Number to find closest to in the list

#### Output

* Index of the closest number to the target

#### Returns

```julia
_find_closest_number_idx([1, 2, 2.7, 3.2, 4, 3.1, 7], 3)

# 6
```
"""
function _find_closest_number_idx(
    list::AbstractArray{T,1},
    target::Number,
) where {T<:Number}
    diff_array = abs.(list .- target)
    targetIdx = something(findfirst(isequal(minimum(diff_array)), diff_array), 0)
end

#######################################
#
# DataFrame manipulation
#
#######################################

function add_dataframe_static_rows(a::DataFrame, args...)
    @debug("Adding column(s)")
    for kwargs in args
        @debug(kwargs)
        for k in kwargs
            name = Symbol(k[1])
            code = k[2]
            expanded_code = vec(repeat([k[2]], size(a, 1), 1))
            @debug("Name: $name  Code: $code")
            a[!, name] = expanded_code
            # insert!(a, size(a, 2) + 1, expanded_code, name)
            # rename!(a, Symbol(string("x", size(a, 2))) => name)
        end
    end
    return a
end


text:
using Gadfly
using DataFrames

function rmean(v)
cumsum(v)./[1:length(v)]
end

theme = Theme(panel_fill=color("light grey"), grid_color=color("white"))
geo = Geom.line

xobsdf = DataFrame(hcat(0:49,readcsv("autoreg50fo.csv"; has_header=false)))
names!(xobsdf, [:t,:RNA, :P, :P2, :DNA])
p = plot(melt(xobsdf, "t"), x="t", y="value", color="variable",
    Guide.ylabel("#"), Guide.title(""),
    theme )
draw(PNG("autoreg50fo.png",18cm, 9cm), p)
draw(PDF("autoreg50fo.pdf",18cm, 9cm), p)


text:
using MicrobiomeAnalysis
using Test

@testset "MicrobiomeAnalysis.jl" begin

    @test @isdefined HintikkaXOData
    @test @isdefined OKeefeDSData

end


text:
#################################Test1 : Array Element types
println("==Test 1: Array Element Types==")
println("One is an Integer : ",isa(1, Int64))
println("One is an Number  : ",isa(1, Number))          #At this case it is true, that 1 is a number
println("Array of One is an Array of Integer : ",isa([1],Array{Int64}))
println("Array of One is an Array of Number  : ",isa([1],Array{Number}))  #But isn't [1] an array of number?
println("Make sense \n Parametric type behaves differently than just a subtype matching.  This is to improve the performance such that an Array of Int cannot be confused with a Number type more complex and thus impairing the performance of the algorithms at the lower level.")
println("\n\n\n\n")
#################################Test2 : Union Argument
println("==Test 2: Union Arguments")
x = [1,2,3]
thatfunction(i::Union{Colon,Int64,Array{Int64}}) = x[i]
println("Int?",thatfunction(1))
println("ArrayInt?",thatfunction([1,2]))
println("Colon?",thatfunction(:))
println("Oh wait... It works...")


text:
@with_kw type RadiusRandom <: Policy
    radius::Float64     = 1.0
    rng::AbstractRNG    = Base.GLOBAL_RNG
end

POMCPOW.next_action(gen::RadiusRandom, pomdp::POMDP, b, h) = rand_in_radius(gen.rng, gen.radius)
function rand_in_radius(rng::AbstractRNG, radius::Float64)
    r = radius*rand(rng)
    angle = 2*pi*rand(rng)
    return r*@SVector([cos(angle), sin(angle)])
end

POMDPs.action(p::RadiusRandom, b) = rand_in_radius(p.rng, p.radius)
POMDPs.updater(p::RadiusRandom) = VoidUpdater()

Base.rand(rng::AbstractRNG, rr::RadiusRandom) = rand_in_radius(rng, rr.radius)

@with_kw type AdaptiveRadiusRandom <: Policy
    max_radius::Float64     = 1.0
    to_zero_first::Bool     = true
    to_light_second::Bool   = true
    rng::AbstractRNG        = Base.GLOBAL_RNG
end

function POMCPOW.next_action(gen::AdaptiveRadiusRandom, pomdp::POMDP, b, h::BeliefNode)
    if gen.to_zero_first && n_children(h) < 1
        m = mean(b)
        n = norm(m)
        return clip_to_max(-m, gen.max_radius)
    elseif gen.to_light_second && n_children(h) < 2
        m = mean(b)
        n = norm(m)
        return clip_to_max(Vec2(pomdp.min_noise_loc-m[1], 0.0), gen.max_radius)
    else
        rad = min(gen.max_radius)
        return rand_in_radius(gen.rng, rad)
    end
    #=
    @show h.node
    @show h.tree.tried[h.node]
    @show h.tree.total_n[h.node]
    @show h.tree.n[first(h.tree.tried[h.node])]
    =#
end

@with_kw type SimpleFeedback <: Policy
    gain::Float64 = 1.0
    max_radius::Float64 = 1.0
end

function clip_to_max(s::AbstractVector{Float64}, max::Float64)
    nrm = norm(s)
    if nrm > max
        return max/nrm*s
    else
        return s
    end
end

POMDPs.action(p::SimpleFeedback, b::SymmetricNormal2) = clip_to_max(-p.gain*b.mean, p.max_radius)
POMDPs.action(p::SimpleFeedback, b) = clip_to_max(-p.gain*mean(b), p.max_radius)
POMDPs.action(p::SimpleFeedback, o::Vec2) = clip_to_max(-p.gain*o, p.max_radius)

immutable LightDarkLQRSolver <: Solver end

function POMDPs.solve(solver::LightDarkLQRSolver, pomdp::LightDark2D)
    A = eye(2)
    B = eye(2)
    K = dlqr(A, B, pomdp.Q, pomdp.R)
    return LightDarkLQRPolicy(K)
end

immutable LightDarkLQRPolicy <: Policy
    K::AbstractMatrix{Float64}
end

POMDPs.action(p::LightDarkLQRPolicy, s::Vec2) = -p.K*s
POMDPs.action(p::LightDarkLQRPolicy, b) = -p.K*mean(b)

immutable InfoGatherHeur <: Policy
    info_x::Float64
    std_thresh::Float64
    homing::Policy
    exploring_up::Updater
    homing_up::Updater
end

# for now, mode can be :exploring or :homing
immutable ModeAugmentedBelief{B}
    mode::Symbol
    b::B
end

@recipe function f(b::ModeAugmentedBelief) b.b end

function POMDPs.action(h::InfoGatherHeur, b::ModeAugmentedBelief)
    if b.mode == :homing
        return action(h.homing, b.b)
    else
        target = Vec2(h.info_x, 0.0)
        diff = mean(b.b) - target
        return action(h.homing, diff)
    end
end

POMDPs.updater(h::InfoGatherHeur) = InfoGatherUpdater(h.std_thresh, h.exploring_up, h.homing_up)

immutable InfoGatherUpdater <: Updater
    std_thresh::Float64
    exploring_up::Updater
    homing_up::Updater
end

POMDPs.initialize_belief(up::InfoGatherUpdater, b) = ModeAugmentedBelief(:exploring, b)

function POMDPs.update(up::InfoGatherUpdater, b::ModeAugmentedBelief, a, o)
    if b.mode != :homing
        bnew = update(up.exploring_up, b.b, a, o)
        std = det(cov(cat(2, particles(bnew)...)'))^(1/4)
        if std < up.std_thresh
            return ModeAugmentedBelief(:homing, bnew)
        else
            return ModeAugmentedBelief(:exploring, bnew)
        end
    else
        bnew = update(up.homing_up, b.b, a, o)
        return ModeAugmentedBelief(:homing, bnew)
    end
end


text:
export print_tree

struct BlockTreeCharSet
    mid
    terminator
    skip
    dash
end

# Default Charset
BlockTreeCharSet() = BlockTreeCharSet('├','└','│','─')

_charwidth(c::Char) = textwidth(c)
_charwidth(s) = sum(map(textwidth, collect(s)))

"""
    print_prefix(io, depth, charset, active_levels)

print prefix of a tree node in a single line.
"""
function print_prefix(io::IO, depth, charset, active_levels)
    for current_depth in 2:depth
        if current_depth in active_levels
            print(io, charset.skip, " "^(_charwidth(charset.dash) + 1))
        else
            print(io, " "^(_charwidth(charset.skip) + _charwidth(charset.dash) + 1))
        end
    end
end

"""
    print_title(io, block)

Print the title of given `block` of an [`AbstractBlock`](@ref).
"""
function print_title(io::IO, x::AbstractBlock)
    printstyled(io, "nqubits: ", nqubits(x); bold=false, color=:cyan)
end

"""
    print_annotation(io, root, node, child, k)

Print the annotation of `k`-th `child` of node, aka the `k`-th element of
`subblocks(node)`.
"""
print_annotation(io::IO, root, node, child, k) = nothing

"""
    print_tree([io=stdout], root)

Print the block tree.
"""
print_tree(root; kwargs...) = print_tree(stdout, root; kwargs...)

Base.show(io::IO, blk::AbstractBlock) = show(io, "plain/text", blk)
Base.show(io::IO, ::MIME"plain/text", blk::AbstractBlock) = print_tree(io, blk)

function Base.show(io::IO, ::MIME"plain/text",
        blk::TagBlock{<:PrimitiveBlock})
    return print_tree(io, blk; title=false, compact=false)
end

Base.show(io::IO, ::MIME"plain/text", blk::PrimitiveBlock) = print_tree(io, blk; title=false, compact=true)

print_tree(io::IO, root::AbstractBlock; kwargs...) = print_tree(io, root, root; kwargs...)

"""
    print_tree(io, root, node[, depth=1, active_levels=()]; kwargs...)

Print the block tree.

# Keywords

- `maxdepth`: max tree depth to print
- `charset`: default is ('├','└','│','─'). See also [`BlockTreeCharSet`](@ref).
- `title`: control whether to print the title, `true` or `false`, default is `true`
"""
function print_tree(
    io::IO,
    root::AbstractBlock,
    node::AbstractBlock,
    depth::Int=1,
    islast::Bool=false,
    active_levels=();
    maxdepth=5, charset=BlockTreeCharSet(), title=true, compact=false)

    if root === node && title
        print_title(io, root)
        println(io)
    end

    if root === node
        islast = true
    end

    print_block(io, node)
    _println(io, node, islast, compact)

    kwargs = (maxdepth=maxdepth, charset=charset, title=false, compact=compact)
    for (k, each_node) in enumerate(subblocks(node))
        if k == lastindex(subblocks(node))
            print_prefix(io, depth, charset, active_levels)
            print(io, charset.terminator, charset.dash)
            print(io, " ")
            print_annotation(io, root, node, each_node, k)
            print_tree(io, root, each_node, depth+1, islast && isempty(subblocks(each_node)), active_levels; kwargs...)
        else
            print_prefix(io, depth, charset, active_levels)
            print(io, charset.mid, charset.dash)
            print(io, " ")
            print_annotation(io, root, node, each_node, k)
            print_tree(io, root, each_node, depth+1, false, (active_levels..., depth+1); kwargs...)
        end
    end

    return nothing
end

function print_tree(
    io::IO,
    root::AbstractBlock,
    node::TagBlock,
    depth::Int=1,
    islast::Bool=false,
    active_levels=();
    maxdepth=5, charset=BlockTreeCharSet(), title=true, compact=false)

    if root === node && title
        print_title(io, root)
        println(io)
    end

    if root === node
        print_annotation(io, root, root, root, 1)
        islast = true
    end

    kwargs = (maxdepth=maxdepth, charset=charset, title=false, compact=compact)
    child = content(node)
    print_annotation(io, root, node, child, 1)
    print_tree(io, root, child, depth+1, islast && isempty(subblocks(child)), active_levels; kwargs...)

    return nothing
end

_println(io::IO, node::AbstractBlock, islast, compact) = !compact && println(io)
function _println(io::IO, node::Union{TagBlock, PrimitiveBlock}, islast, compact)
    if !compact && !islast
        println(io)
    end
end

# Custom layouts
color(m::AbstractBlock) = color(typeof(m))
color(::Type{<:Swap}) = :magenta
color(::Type{<:ControlBlock}) = :red
color(::Type{<:ChainBlock}) = :blue
color(::Type{<:MathGate}) = :red
color(::Type{<:PutBlock}) = :cyan
color(::Type{T}) where {T <: PauliString} = :cyan
color(::Type{<:RepeatedBlock}) = :cyan
color(::Type{<:GeneralMatrixBlock}) = :red

# color(::Type{T}) where {T <: PauliString} = :cyan
# color(::Type{T}) where {T <: Sequential} = :blue

print_block(io::IO, g::PhaseGate) = print(io, "phase(", g.theta, ")")
print_block(io::IO, S::ShiftGate) = print(io, "shift(", S.theta, ")")
print_block(io::IO, R::RotationGate) = print(io, "rot(", content(R), ", ", R.theta, ")")
print_block(io::IO, swap::Swap) = printstyled(io, "swap", swap.locs; bold=true, color=color(Swap))
print_block(io::IO, x::KronBlock) = printstyled(io, "kron"; bold=true, color=color(KronBlock))
print_block(io::IO, x::ChainBlock) = printstyled(io, "chain"; bold=true, color=color(ChainBlock))
print_block(io::IO, x::ReflectGate{N}) where N = print(io, "reflect($(summary(x.psi)))")
print_block(io::IO, c::Concentrator) = print(io, "Concentrator: ", occupied_locs(c))
print_block(io::IO, c::CachedBlock) = print_block(io, content(c))
print_block(io::IO, c::Sum) = printstyled(io, "sum"; bold=true, color=color(ChainBlock))
print_block(io::IO, c::TagBlock) = nothing
print_block(io::IO, c::GeneralMatrixBlock) = printstyled(io, "matblock(...)"; color=color(GeneralMatrixBlock))

function print_block(io::IO, c::Measure{N, K, OT}) where {N, K, OT}
    strs = String[]
    if c.operator != ComputationalBasis()
        push!(strs, "operator=$(repr(c.operator))")
    end

    if c.locations != AllLocs()
        push!(strs, "locs=$(repr(c.locations))")
    end

    if c.collapseto !== nothing
        push!(strs, "collapseto=$(c.collapseto)")
    end

    if c.remove
        push!(strs, "remove=true")
    end

    out = join(strs, ", ")
    if !isempty(strs)
        out = "Measure($N;" * out
    else
        out = "Measure($N" * out
    end

    return print(io, out, ")")
end

# TODO: use OhMyREPL's default syntax highlighting for functions
function print_block(io::IO, m::MathGate{N, <:LegibleLambda}) where N
    printstyled(io, "mathgate($(m.f); nbits=$N, bview=$(nameof(m.v)))"; bold=true, color=color(m))
end

function print_block(io::IO, m::MathGate{N, <:Function}) where N
    printstyled(io, "mathgate($(nameof(m.f)); nbits=$N, bview=$(nameof(m.v)))"; bold=true, color=color(m))
end

function print_block(io::IO, te::TimeEvolution)
    println(io, "Time Evolution Δt = $(te.dt), tol = $(te.tol)")
    print_tree(io, content(te.H); title=false)
end

function print_block(io::IO, x::ControlBlock)
    printstyled(io, "control("; bold=true, color=color(ControlBlock))

    for i in eachindex(x.ctrl_locs)
        printstyled(io, x.ctrl_locs[i]; bold=true, color=color(ControlBlock))

        if i != lastindex(x.ctrl_locs)
            printstyled(io, ", "; bold=true, color=color(ControlBlock))
        end
    end
    printstyled(io, ")"; bold=true, color=color(ControlBlock))
end

function print_block(io::IO, pb::PutBlock{N}) where N
    printstyled(io, "put on ("; bold=true, color=color(PutBlock))
    for i in eachindex(pb.locs)
        printstyled(io, pb.locs[i]; bold=true, color=color(PutBlock))
        if i != lastindex(pb.locs)
            printstyled(io, ", "; bold=true, color=color(PutBlock))
        end
    end
    printstyled(io, ")"; bold=true, color=color(PutBlock))
end

function print_block(io::IO, rb::RepeatedBlock{N}) where N
    printstyled(io, "repeat on ("; bold=true, color=color(RepeatedBlock))
    for i in eachindex(rb.locs)
        printstyled(io, rb.locs[i]; bold=true, color=color(RepeatedBlock))
        if i != lastindex(rb.locs)
            printstyled(io, ", "; bold=true, color=color(RepeatedBlock))
        end
    end
    printstyled(io, ")"; bold=true, color=color(RepeatedBlock))
end

function print_block(io::IO, x::PauliString)
    printstyled(io, "PauliString"; bold=true, color=color(PauliString))
end

# forward to simplify interfaces
function print_annotation(io::IO,
    root::AbstractBlock,
    node::AbstractBlock,
    child::AbstractBlock, k=1)
    print_annotation(io, child)
end

print_annotation(io::IO, node::AbstractBlock) = nothing # skip
print_annotation(io::IO, c::Daggered) = printstyled(io, " [†]"; bold=true, color=:yellow)
print_annotation(io::IO, c::CachedBlock) =
    printstyled(io, "[cached] "; bold=true, color=:yellow)

function print_annotation(io::IO, x::Scale)
    if x.alpha == im
        printstyled(io, "[+im] "; bold=true, color=:yellow)
    elseif x.alpha == -im
        printstyled(io, "[-im] "; bold=true, color=:yellow)
    elseif x.alpha == 1
        printstyled(io, "[+] "; bold=true, color=:yellow)
    elseif x.alpha == -1
        printstyled(io, "[-] "; bold=true, color=:yellow)
    elseif real(x.alpha) == 0
        printstyled(io, "[scale: ", imag(x.alpha), "im] "; bold=true, color=:yellow)
    else
        printstyled(io, "[scale: ", x.alpha, "] "; bold=true, color=:yellow)
    end
end

function print_annotation(io::IO, x::Scale{Val{S}}) where S
    print_annotation(io, Scale(S, x))
end

function print_annotation(
    io::IO,
    root::AbstractBlock,
    node::KronBlock,
    child::AbstractBlock,
    k::Int)

    printstyled(io, node.locs[k]; bold=true, color=:white)
    print(io, "=>")
    print_annotation(io, child)
end

function print_annotation(
    io::IO,
    root::AbstractBlock,
    node::ControlBlock,
    child::AbstractBlock,
    k::Int)

    printstyled(io, node.locs; bold=true, color=:white)
    print(io, " ")
    print_annotation(io, child)
end


text:
@testset "Test Function Tests" begin

  @test isdefined(Julz, :test) == true

  dummy_package = "$(dirname(@__FILE__))/../../../../JulzDummy"

  cp("dummy", dummy_package, remove_destination=true)

  err_data = "nil"
  out_data = "nil"
  did_break = false

  cd(dummy_package) do

    originalSTDERR = STDERR
    originalSTDOUT = STDOUT

    (errRead, errWrite) = redirect_stderr()
    (outRead, outWrite) = redirect_stdout()

    println("\ndummy start\n")
    try
      Julz.test()
    catch
      did_break = true
    end
    println("dummy end")

    close(errWrite)
    close(outWrite)

    err_data = readavailable(errRead)
    out_data = readavailable(outRead)

    close(errRead)
    close(outRead)

    redirect_stderr(originalSTDERR)
    redirect_stdout(originalSTDOUT)

  end

  rm(dummy_package, force=true, recursive=true)

  @test contains(String(err_data), "JulzDummy tests passed")

  if did_break

    println(" ================================ ")
    println("  dummy broken test output (out)  ")
    println(" ================================ ")

    println(String(out_data))

    println(" ================================ ")
    println("  dummy broken test output (err)  ")
    println(" ================================ ")

    println(String(err_data))

    println(" ================================ ")
    println("  dummy broken test output (end)  ")
    println(" ================================ ")

  end

end


text:
using Test

@testset "all" begin
    @testset "FED" begin include("FED.jl") end
end


text:
"""
```
AbstractDSGEModel{T} <: AbstractModel{T}
```

The AbstractDSGEModel is defined as a subtype of AbstractModel
"""
abstract type AbstractDSGEModel{T} <: ModelConstructors.AbstractModel{T} end

function Base.show(io::IO, m::AbstractDSGEModel)
    @printf io "Dynamic Stochastic General Equilibrium Model\n"
    @printf io "no. states:             %i\n" n_states(m)
    @printf io "no. anticipated shocks: %i\n" n_anticipated_shocks(m)
    @printf io "data vintage:           %s\n" data_vintage(m)
    @printf io "description:\n %s\n"          description(m)
end

# Number of anticipated policy shocks
n_anticipated_shocks(m::AbstractDSGEModel) = get_setting(m, :n_anticipated_shocks)
n_anticipated_shocks_padding(m::AbstractDSGEModel) = get_setting(m, :n_anticipated_shocks_padding)

# Dates, indices, number of periods for each regime
date_presample_start(m::AbstractDSGEModel) = get_setting(m, :date_presample_start)
date_mainsample_start(m::AbstractDSGEModel) = get_setting(m, :date_mainsample_start)
date_zlb_start(m::AbstractDSGEModel) = get_setting(m, :date_zlb_start)
date_zlb_end(m::AbstractDSGEModel) = get_setting(m, :date_zlb_end)
date_presample_end(m::AbstractDSGEModel) = Dates.lastdayofquarter(get_setting(m, :date_mainsample_start) - Dates.Month(3))
date_prezlb_end(m::AbstractDSGEModel) = Dates.lastdayofquarter(get_setting(m, :date_zlb_start) - Dates.Month(3))
date_mainsample_end(m::AbstractDSGEModel) = Dates.lastdayofquarter(get_setting(m, :date_forecast_start) - Dates.Month(3))
date_conditional_end(m::AbstractDSGEModel) = get_setting(m, :date_conditional_end)

index_presample_start(m::AbstractDSGEModel) = 1
index_mainsample_start(m::AbstractDSGEModel) = subtract_quarters(date_mainsample_start(m), date_presample_start(m)) + 1
index_zlb_start(m::AbstractDSGEModel) = subtract_quarters(date_zlb_start(m), date_presample_start(m)) + 1
index_forecast_start(m::AbstractDSGEModel) = subtract_quarters(date_forecast_start(m), date_presample_start(m)) + 1

"""
```
index_shockdec_start(m::AbstractDSGEModel)
```

Returns the index starting from which the shock decomposition is saved, where 1 is the index corresponding to date_mainsample_start(m).
"""
index_shockdec_start(m::AbstractDSGEModel) = subtract_quarters(date_shockdec_start(m), date_mainsample_start(m)) + 1

"""
```
index_shockdec_end(m::AbstractDSGEModel)
```

Returns the last index for which the shock decomposition is saved, where 1 is the index corresponding to date_mainsample_start(m).
"""
index_shockdec_end(m::AbstractDSGEModel) = subtract_quarters(date_shockdec_end(m), date_mainsample_start(m)) + 1

n_presample_periods(m::AbstractDSGEModel)   = subtract_quarters(date_mainsample_start(m), date_presample_start(m))
n_prezlb_periods(m::AbstractDSGEModel)      = subtract_quarters(date_zlb_start(m), date_mainsample_start(m))
n_zlb_periods(m::AbstractDSGEModel)         = subtract_quarters(date_forecast_start(m), date_zlb_start(m))
n_mainsample_periods(m::AbstractDSGEModel)  = subtract_quarters(date_forecast_start(m), date_mainsample_start(m))
n_conditional_periods(m::AbstractDSGEModel) = subtract_quarters(date_conditional_end(m), date_mainsample_end(m))

inds_presample_periods(m::AbstractDSGEModel)  = collect(index_presample_start(m):(index_mainsample_start(m)-1))
inds_prezlb_periods(m::AbstractDSGEModel)     = collect(index_mainsample_start(m):(index_zlb_start(m)-1))
inds_zlb_periods(m::AbstractDSGEModel)        = collect(index_zlb_start(m):(index_forecast_start(m)-1))
inds_mainsample_periods(m::AbstractDSGEModel) = collect(index_mainsample_start(m):(index_forecast_start(m)-1))

# Convenience functions for working with heterogeneous agent models
# that differentiate between backward looking "state" variables and "jump" variables
n_backward_looking_states(m::AbstractDSGEModel) = get_setting(m, :n_backward_looking_states)
n_jumps(m::AbstractDSGEModel) = get_setting(m, :n_jumps)
n_model_states(m::AbstractDSGEModel) = get_setting(m, :n_model_states)

# The numbers for n_states, and n_jumps assumes normalization
# There is a procedure in klein_solve that normalizes the state variable grids
# by removing an entry from them. Hence the number of states and jumps being
# tracked should be 1 less if normalized
# However, the number of jumps and states unnormalized are required for
# the construction of the Jacobian, hence the reason for these helpers
n_jumps_unnormalized(m::AbstractDSGEModel) = n_jumps(m) + get_setting(m, :jumps_normalization_factor)
n_backward_looking_states_unnormalized(m::AbstractDSGEModel) = n_backward_looking_states(m) + get_setting(m, :backward_looking_states_normalization_factor)
n_model_states_unnormalized(m::AbstractDSGEModel) = n_jumps_unnormalized(m) + n_backward_looking_states_unnormalized(m)
n_model_states_original(m::AbstractDSGEModel) = get_setting(m, :n_model_states_original)

"""
```
AbstractCTModel{T} <: AbstractDSGEModel{T}
```

The AbstractCTModel is defined as a subtype of AbstractDSGEModel to accomodate the
numerical methods and procedures specific to continuous time models.
"""
abstract type AbstractCTModel{T} <: AbstractDSGEModel{T} end

n_states(m::AbstractCTModel) = sum(map(i -> length(collect(m.endogenous_states)[i][2]), 1:length(keys(m.endogenous_states))))
n_shocks_expectational(m::AbstractCTModel) = sum(map(i -> length(collect(m.expected_shocks)[i][2]), 1:length(keys(m.expected_shocks))))

"""
```
AbstractHetModel{T} <: AbstractDSGEModel{T}
```

The AbstractHetModel is defined as a subtype of AbstractDSGEModel to accomodate a bunch of stuff, but for now just different impulse response functions.
"""
abstract type AbstractHetModel{T} <: AbstractDSGEModel{T} end

"""
```
AbstractRepModel{T} <: AbstractDSGEModel{T}
```

The AbstractRepresentativeModel is defined as a subtype of AbstractDSGEModel to accomodate a bunch of stuff, but for now just different impulse response functions.
"""
abstract type AbstractRepModel{T} <: AbstractDSGEModel{T} end

"""
```
get_dict(m, class, index)
```
"""
function get_dict(m::AbstractDSGEModel, class::Symbol)
    if class == :states
        m.endogenous_states
    elseif class == :obs
        m.observables
    elseif class == :pseudo
        m.pseudo_observables
    elseif class in [:shocks, :stdshocks]
        m.exogenous_shocks
    else
        throw(ArgumentError("Invalid class: $class. Must be :states, :obs, :pseudo, :shocks, or :stdshocks"))
    end
end

"""
```
get_key(m, class, index)
```
Returns the name of the state (`class = :states`), observable (`:obs`),
pseudo-observable (`:pseudo`), or shock (`:shocks` or `:stdshocks`)
corresponding to the given `index`.
"""
function get_key(m::AbstractDSGEModel, class::Symbol, index::Int)
    dict = get_dict(m, class)
    out = Base.filter(key -> dict[key] == index, collect(keys(dict)))
    if length(out) == 0
        error("Key corresponding to index $index not found for class: $class")
    elseif length(out) > 1
        error("Multiple keys corresponding to index $index found for class: $class")
    else
        return out[1]
    end
end

# Parse population mnemonic into 2 Nullable{Symbol}s from one
function parse_population_mnemonic(m::AbstractDSGEModel)
    mnemonic = get_setting(m, :population_mnemonic)
    if isnull(mnemonic)
        return [Nullable{Symbol}(), Nullable{Symbol}()]
    else
        return map(s -> Nullable(Symbol(s)), split(string(get(mnemonic)), DSGE_DATASERIES_DELIM))
    end
end

# From an augmented state space with anticipated policy shocks, get indices
# corresponding to pre-ZLB states, shocks, and observables
function inds_states_no_ant(m::AbstractDSGEModel)
    if n_anticipated_shocks(m) > 0
        ind_ant1 = m.endogenous_states[:rm_tl1]
        ind_antn = m.endogenous_states[Symbol("rm_tl$(n_anticipated_shocks(m))")]
        return [1:(ind_ant1-1); (ind_antn+1):n_states_augmented(m)]
    else
        return collect(1:n_states_augmented(m))
    end
end

function inds_shocks_no_ant(m::AbstractDSGEModel)
    if n_anticipated_shocks(m) > 0
        ind_ant1 = m.exogenous_shocks[:rm_shl1]
        ind_antn = m.exogenous_shocks[Symbol("rm_shl$(n_anticipated_shocks(m))")]
        return [1:(ind_ant1-1); (ind_antn+1):n_shocks_exogenous(m)]
    else
        return collect(1:n_shocks_exogenous(m))
    end
end

function inds_obs_no_ant(m::AbstractDSGEModel)
    if n_anticipated_shocks(m) > 0
        ind_ant1 = m.observables[:obs_nominalrate1]
        ind_antn = m.observables[Symbol("obs_nominalrate$(n_anticipated_shocks(m))")]
        return [1:(ind_ant1-1); (ind_antn+1):n_observables(m)]
    else
        return collect(1:n_observables(m))
    end
end

# Interface for accessing parameters
get_parameters(m::AbstractDSGEModel) = m.parameters

# Interface for accessing rng
get_rng(m::AbstractDSGEModel) = m.rng

# Interface for accessing settings dictionary
get_settings(m::AbstractDSGEModel) = m.settings

# Interface for accessing observables dictionary
get_observables(m::AbstractDSGEModel) = m.observables
get_pseudo_observables(m::AbstractDSGEModel) = m.pseudo_observables
get_exogenous_shocks(m::AbstractDSGEModel) = m.exogenous_shocks

# Interface for data
cond_vintage(m::AbstractDSGEModel)    = get_setting(m, :cond_vintage)
cond_id(m::AbstractDSGEModel)         = get_setting(m, :cond_id)
cond_full_names(m::AbstractDSGEModel) = get_setting(m, :cond_full_names)
cond_semi_names(m::AbstractDSGEModel) = get_setting(m, :cond_semi_names)
use_population_forecast(m::AbstractDSGEModel) = get_setting(m, :use_population_forecast)
hpfilter_population(m::AbstractDSGEModel)     = get_setting(m, :hpfilter_population)

# Interface for general computation settings
use_parallel_workers(m::AbstractDSGEModel)    = get_setting(m, :use_parallel_workers)

# Interface for estimation settings
reoptimize(m::AbstractDSGEModel)          = get_setting(m, :reoptimize)
calculate_hessian(m::AbstractDSGEModel) = get_setting(m, :calculate_hessian)
hessian_path(m::AbstractDSGEModel)      = get_setting(m, :hessian_path)
n_hessian_test_params(m::AbstractDSGEModel) = get_setting(m, :n_hessian_test_params)

# Interface for Metropolis-Hastings settings
n_mh_blocks(m::AbstractDSGEModel)       =  get_setting(m, :n_mh_blocks)
n_mh_param_blocks(m::AbstractDSGEModel) =  get_setting(m, :n_mh_param_blocks)
n_mh_simulations(m::AbstractDSGEModel)  =  get_setting(m, :n_mh_simulations)
n_mh_burn(m::AbstractDSGEModel)         =  get_setting(m, :n_mh_burn)
mh_thin(m::AbstractDSGEModel)           =  get_setting(m, :mh_thin)

# Interface for forecast settings
date_forecast_start(m::AbstractDSGEModel)   = get_setting(m, :date_forecast_start)
forecast_block_size(m::AbstractDSGEModel)   = get_setting(m, :forecast_block_size)
forecast_start_block(m::AbstractDSGEModel)  = get_setting(m, :forecast_start_block)
forecast_input_file_overrides(m::AbstractDSGEModel) = get_setting(m, :forecast_input_file_overrides)
forecast_uncertainty_override(m::AbstractDSGEModel) = get_setting(m, :forecast_uncertainty_override)
forecast_smoother(m::AbstractDSGEModel)     = get_setting(m, :forecast_smoother)
forecast_tdist_df_val(m::AbstractDSGEModel) = get_setting(m, :forecast_tdist_df_val)
forecast_tdist_shocks(m::AbstractDSGEModel) = get_setting(m, :forecast_tdist_shocks)
forecast_zlb_value(m::AbstractDSGEModel)    = get_setting(m, :forecast_zlb_value)
impulse_response_horizons(m::AbstractDSGEModel) = get_setting(m, :impulse_response_horizons)
n_shockdec_periods(m::AbstractDSGEModel)    = index_shockdec_end(m) - index_shockdec_start(m) + 1

# Interface for alternative policy settings
alternative_policy(m::AbstractDSGEModel) = get_setting(m, :alternative_policy)

function date_forecast_end(m::AbstractDSGEModel)
    date = date_forecast_start(m) + Dates.Month(3 * (forecast_horizons(m)-1))
    return Dates.lastdayofquarter(date)
end

function forecast_horizons(m::AbstractDSGEModel; cond_type::Symbol = :none)
    horizons = get_setting(m, :forecast_horizons)
    if cond_type == :none
        return horizons
    else
        return horizons - n_conditional_periods(m)
    end
end

function date_shockdec_start(m::AbstractDSGEModel)
    startdate = get_setting(m, :shockdec_startdate)
    if (Nullables.isnull(startdate)) | (startdate===nothing)
        return DSGE.date_mainsample_start(m)
    else
        return get(startdate)
    end
end

function date_shockdec_end(m::AbstractDSGEModel)
    enddate =  get_setting(m, :shockdec_enddate)
    if !isnull(enddate)
        return get(enddate)
    else
        return date_forecast_end(m)
    end
end

"""
```
load_parameters_from_file(m::AbstractDSGEModel,path::String)
```
Returns a vector of parameters, read from a file, suitable for updating `m`.
"""
function load_parameters_from_file(m::AbstractDSGEModel, path::String)

    if isfile(path) && splitext(path)[2] == ".h5"
        x  = h5open(path, "r") do file
            try
                read(file, "params")
            catch
                error("$path does not contain variable params")
            end
        end
    else
        error("$path is not a valid HDF5 file.")
    end

    if m.spec=="smets_wouters" && length(x)==42
        x = vcat(x, zeros(7))
    end

    @assert length(x) == length(m.parameters)
    @assert eltype(x) == typeof(m.parameters[1].value)
    return x
end

"""
```
specify_mode!(m::AbstractDSGEModel, mode_file::String=""; verbose=:low)
```

Updates the values of `m.parameters` with the values from
`mode_file`. Sets `reoptimize` setting to `false`.

Usage: should be run before calling `estimate(m)`, e.g.:

    m = Model990()
    specify_mode!(m, modefile)
    estimate(m)
"""
function specify_mode!(m::AbstractDSGEModel, mode_file::String = ""; verbose=:low)

    m <= Setting(:reoptimize, false)

    if mode_file == ""
        mode_file = inpath(m, "user", "paramsmode.h5")
    end
    mode_file = normpath(mode_file)

    DSGE.update!(m,load_parameters_from_file(m,mode_file))

    println(verbose, :low, "Loaded previous mode from $mode_file.")
end

"""
```
specify_hessian!(m::AbstractDSGEModel, path::String=""; verbose=:low)
```

Specify a Hessian matrix calculated at the posterior mode to use in the model estimation. If
no path is provided, will attempt to detect location.
"""
function specify_hessian!(m::AbstractDSGEModel, path::String=""; verbose=:low)
    if isempty(path)
        path = inpath(m, "user", "hessian.h5")
    end

    if isfile(path) && splitext(path)[2] == ".h5"
        m <= Setting(:hessian_path, normpath(abspath(path)))
    else
        error("Invalid input Hessian file: $path",)
    end

    m <= Setting(:calculate_hessian, false)

    println(verbose, :low, "Specified hessian from $path.")
end

"""
```
transform_to_model_space!(m::AbstractDSGEModel, values::Vector{T}) where T<:AbstractFloat
```

Transforms `values` from the real line to the model space, and assigns `values[i]` to
`m.parameters[i].value` for non-steady-state parameters. Recomputes the steady-state
paramter values.

### Arguments
- `m`: the model object
- `values`: the new values to assign to non-steady-state parameters.
"""
function transform_to_model_space!(m::AbstractDSGEModel, values::Vector{T}) where {T<:AbstractFloat}
    new_values = transform_to_model_space(m.parameters, values)
    DSGE.update!(m, new_values)
    steadystate!(m)
end

"""
```
update!(m::AbstractDSGEModel, values::Vector{T}) where T<:AbstractFloat
```

Update `m.parameters` with `values`, recomputing the steady-state parameter values.

### Arguments:
- `m`: the model object
- `values`: the new values to assign to non-steady-state parameters.
"""
function update!(m::AbstractDSGEModel, values::Vector{T}) where T<:AbstractFloat
    ModelConstructors.update!(m.parameters, values)
    steadystate!(m)
end

"""
```
update!(m::AbstractDSGEModel, values::ParameterVector{T}) where T
```
Update `m.parameters` with `values`, recomputing the steady-state parameter values.
### Arguments:
- `m`: the model object
- `values`: the new values to assign to non-steady-state parameters.
"""
function update!(m::AbstractDSGEModel, values::ParameterVector{T}) where T
    ModelConstructors.update!(m.parameters, [θ.value for θ in values])
    steadystate!(m)
end

"""
```
mutable struct ShockGroup
```

The `ShockGroup` mutable struct is used in `prepare_means_table_shockdec` and
`plot_shock_decompositions`. When plotting shock decompositions, we usually want
to group the shocks into categories (financial, monetary policy, etc.) so that
the resulting grouped bar plot is legible.

### Fields

- `name::String`
- `shocks::Vector{Symbol}`
- `color::Colorant`
"""
mutable struct ShockGroup
    name::String
    shocks::Vector{Symbol}
    color::Colorant
end

function ShockGroup(name::String, shocks::Vector{Symbol}, color_name::Symbol)
    color = parse(Colorant, color_name)
    return ShockGroup(name, shocks, color)
end

mutable struct SteadyStateConvergenceError <: Exception
    msg::String
end
SteadyStateConvergenceError() = SteadyStateConvergenceError("SteadyState didn't converge")
Base.showerror(io::IO, ex::SteadyStateConvergenceError) = print(io, ex.msg)


text:
# ---
# title: 459. Repeated Substring Pattern
# id: problem459
# author: Indigo
# date: 2021-02-04
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/repeated-substring-pattern/description/>
# hidden: true
# ---
#
# Given a non-empty string check if it can be constructed by taking a substring
# of it and appending multiple copies of the substring together. You may assume
# the given string consists of lowercase English letters only and its length
# will not exceed 10000.
#
#
#
# **Example 1:**
#
#
#
#     Input: "abab"
#     Output: True
#     Explanation: It's the substring "ab" twice.
#
#
# **Example 2:**
#
#
#
#     Input: "aba"
#     Output: False
#
#
# **Example 3:**
#
#
#
#     Input: "abcabcabcabc"
#     Output: True
#     Explanation: It's the substring "abc" four times. (And the substring "abcabc" twice.)
#
#
#
## @lc code=start
using LeetCode

function repeated_substring_pattern(s::String)
    findnext(s, s^2, 2) != length(s) + 1 : 2 * length(s)
end
## @lc code=end


text:
abstract Configuration

immutable Particle
    index::Int
    species::Int
end

abstract ParticlesConfiguration <: Configuration

function get_Nspecies(config::ParticlesConfiguration)
    @assert length(config.particle_positions) == length(config.lattice_status)
    return length(config.particle_positions)
end

function get_Nparticles(config::ParticlesConfiguration)
    return map(length, config.particle_positions)
end

function get_Ntotal(config::ParticlesConfiguration)
    return sum(get_Nparticles(config))
end

type LatticeParticlesConfiguration <: ParticlesConfiguration
    particle_positions::Vector{Vector{Int}}
    lattice_status::Vector{Vector{Int}}

    function LatticeParticlesConfiguration(particle_positions::Vector{Vector{Int}}, Nsites::Int)
        Nparticles = map(length, particle_positions)
        lattice_status = Vector{Int}[]
        for j in 1:length(particle_positions) # = Nspecies
            @assert findmin(particle_positions[j])[1] >= 1
            @assert findmax(particle_positions[j])[1] <= Nsites
            push!(lattice_status, zeros(Int, Nsites))
            for i in 1:Nparticles[j]
                lattice_status[j][particle_positions[j][i]] = i
            end
        end
        return new(particle_positions, lattice_status)
    end
end

function get_Nsites(config::LatticeParticlesConfiguration)
    @assert length(config.lattice_status) != 0
    # @assert all(length(config.lattice_status[1]) .== map(length, config.lattice_status[2:end]))
    return length(config.lattice_status[1])
end

# this is probably unnecessary?
function is_valid_configuration(config::LatticeParticlesConfiguration)
    guy = true
    #FIXME: loop over species and do some checks..
    return guy
end

function is_valid_particle(config::LatticeParticlesConfiguration, particle::Particle)
    return (particle.species <= get_Nspecies(config) &&
            particle.index <= get_Nparticles(config)[particle.species])
end

function is_occupied(config::LatticeParticlesConfiguration, position::Int, species::Int)
    return config.lattice_status[species][position] != 0
end

function GutzSpinhalfConfiguration(positions_up::Vector{Int}, Nsites)
    @assert all(positions_up .<= Nsites)
    positions_down = setdiff(collect(1:Nsites), positions_up)
    return LatticeParticlesConfiguration(Vector{Int}[positions_up, positions_down], Nsites)
end

function is_valid_Gutz_spinhalf_configuration(config::LatticeParticlesConfiguration)
    guy = true
    # guy &= get_Nspecies(config) == 2 # spin-1/2 spinons
    # guy &= dot(config.lattice_status[1], config.lattice_status[2]) == 0 # no double occupancy [this is an O(N) operations check..]
    # guy &= isempty(setdiff([1:config.Nsites], union(config.particle_positions[1], config.particle_positions[1]))) # particles cover the lattice
    # the above is fairly time-consuming, so use this somewhat weaker check instead
    # guy &= get_Ntotal(config) == get_Nsites(config)
    return guy
end

# takes config_alpha1 and config_alpha2 and returns a tuple of the 2 swapped configurations (beta1, beta2)
# specific to swaps with equivalent particle numbers in the subsystem (for each species), which should be all we ever need..
function swap(config_alpha1::LatticeParticlesConfiguration, config_alpha2::LatticeParticlesConfiguration,
              subsystemA::Vector{Int})
    @assert get_Nspecies(config_alpha1) == get_Nspecies(config_alpha2)
    @assert get_Nparticles(config_alpha1) == get_Nparticles(config_alpha1)
    @assert get_Nsites(config_alpha1) == get_Nsites(config_alpha2)
    Nsites = get_Nsites(config_alpha1)

    @assert 0 < length(subsystemA) < Nsites # subsystemA is not the whole system nor empty
    @assert minimum(subsystemA) >= 1 && maximum(subsystemA) <= Nsites

    # check that there are nA particles in subsystem A for both configs
    nA = get_Nparticles_in_subsystem(config_alpha1, subsystemA)
    @assert nA == get_Nparticles_in_subsystem(config_alpha2, subsystemA)

    swapped_particle_positions_beta1 = Vector{Int}[]
    swapped_particle_positions_beta2 = Vector{Int}[]
    for s in 1:get_Nspecies(config_alpha1)
        push!(swapped_particle_positions_beta1, [config_alpha2.particle_positions[s][1:nA[s]]; config_alpha1.particle_positions[s][nA[s]+1:end]])
        push!(swapped_particle_positions_beta2, [config_alpha1.particle_positions[s][1:nA[s]]; config_alpha2.particle_positions[s][nA[s]+1:end]])
    end
    @assert all(map(length, swapped_particle_positions_beta1) .== map(length, swapped_particle_positions_beta2) .== get_Nparticles(config_alpha1))

    return LatticeParticlesConfiguration(swapped_particle_positions_beta1, Nsites),
           LatticeParticlesConfiguration(swapped_particle_positions_beta2, Nsites)
end

function get_Nparticles_in_subsystem(config::LatticeParticlesConfiguration, subsystem::Vector{Int})
    @assert 0 <= length(subsystem) <= get_Nsites(config)
    if length(subsystem) == 0
        return 0
    else
        return map(v -> length(find(v[subsystem] .!= 0)), config.lattice_status)
    end
end

function initialize_random_hardcore_two_copy_configuration(Nsites::Int, Nparticles::Int, subsystemA::Vector{Int}, nA::Int)
    Nsites_subsystemA = length(subsystemA)

    @assert Nparticles <= Nsites
    @assert 0 < Nsites_subsystemA < Nsites # subsystemA is not the whole system nor empty
    @assert minimum(subsystemA) >= 1 && maximum(subsystemA) <= Nsites
    @assert 0 <= nA <= Nparticles
    @assert nA <= Nsites_subsystemA && Nparticles - nA <= Nsites - Nsites_subsystemA # since we have hard-core particles here

    subsystemB = setdiff(collect(1:Nsites), subsystemA)

    positions_alpha1_subsystemA = subsystemA[random_combination(Nsites_subsystemA, nA)]
    positions_alpha1_subsystemB = subsystemB[random_combination(Nsites - Nsites_subsystemA, Nparticles - nA)]
    positions_alpha1 = [positions_alpha1_subsystemA; positions_alpha1_subsystemB]

    positions_alpha2_subsystemA = subsystemA[random_combination(Nsites_subsystemA, nA)]
    positions_alpha2_subsystemB = subsystemB[random_combination(Nsites - Nsites_subsystemA, Nparticles - nA)]
    positions_alpha2 = [positions_alpha2_subsystemA; positions_alpha2_subsystemB]

    config_alpha1 = LatticeParticlesConfiguration(Vector{Int}[positions_alpha1], Nsites)
    config_alpha2 = LatticeParticlesConfiguration(Vector{Int}[positions_alpha2], Nsites)
    @assert get_Nparticles_in_subsystem(config_alpha1, subsystemA)[1] == get_Nparticles_in_subsystem(config_alpha2, subsystemA)[1] == nA
    @assert get_Nparticles_in_subsystem(config_alpha1, subsystemB)[1] == get_Nparticles_in_subsystem(config_alpha2, subsystemB)[1] == Nparticles - nA

    config_beta1, config_beta2 = swap(config_alpha1, config_alpha2, subsystemA)
    @assert get_Nparticles_in_subsystem(config_beta1, subsystemA)[1] == get_Nparticles_in_subsystem(config_beta2, subsystemA)[1] == nA
    @assert get_Nparticles_in_subsystem(config_beta1, subsystemB)[1] == get_Nparticles_in_subsystem(config_beta2, subsystemB)[1] == Nparticles - nA

    return config_alpha1, config_alpha2, config_beta1, config_beta2, subsystemB
end

function exchange_particles!(config::LatticeParticlesConfiguration, particle1_ind::Int, particle2_ind::Int, species::Int)
    @assert species <= get_Nspecies(config)
    @assert particle1_ind <= get_Nparticles(config)[species]
    @assert particle2_ind <= get_Nparticles(config)[species]
    @assert particle1_ind != particle2_ind

    particle1_orig_pos = config.particle_positions[species][particle1_ind]
    particle2_orig_pos = config.particle_positions[species][particle2_ind]

    config.particle_positions[species][particle1_ind], config.particle_positions[species][particle2_ind] = particle2_orig_pos, particle1_orig_pos
    config.lattice_status[species][particle1_orig_pos], config.lattice_status[species][particle2_orig_pos] = particle2_ind, particle1_ind

    return config
end


text:
# LibTOAST.jl: interface to the TOAST++ library
# Copyright (C) 2019 Samuel Powell

@testset "Assembly" begin

  # Check we can create system matrices and convert to Julia
  S1 = SystemMatrix(mesh3D)
  assemble!(S1, FF)
  @test norm(sparse(S1),1)-norm(sparse(assemble(mesh3D, FF)),1) < 1e-14

  # Check BNDFF integral
  S1 = SystemMatrix(mesh3D)
  assemble!(S1, BNDFF)
  @test norm(sparse(S1),1)-norm(sparse(assemble(mesh3D, BNDFF)),1) < 1e-14

  # Check FF integral
  @test sum(sparse(assemble(mesh3D, FF))) ≈ fullsize(mesh3D)

  # Check DD integral
  @test sum(sparse(assemble(mesh3D, DD))) < 1e-12

  # Check we can add to system matrices in place
  assemble!(S1, DD)
  S2 = sparse(assemble(mesh3D, FF)) + sparse(assemble(mesh3D, DD))
  @test norm(sparse(S1),1)-norm(sparse(S2),1) < 1e-14

  # Check we can assemble with parameters
  p1 = NodalCoeff(mesh3D,rand(nodecount(mesh3D)))
  p2 = NodalCoeff(mesh3D,rand(nodecount(mesh3D)))
  p3 = NodalCoeff(mesh3D,rand(nodecount(mesh3D)))

  S1 = SystemMatrix(mesh3D)
  assemble!(S1, PFF, p1)
  assemble!(S1, PDD, p2)
  assemble!(S1, BNDPFF, p3)

  S2 = sparse(assemble(mesh3D, PFF, p1)) +
       sparse(assemble(mesh3D, PDD, p2)) +
       sparse(assemble(mesh3D, BNDPFF, p3))

  @test norm(sparse(S1),1)-norm(S2,1) < 1e-14

end


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ d627801e-3479-11eb-1fab-1dfaf4af1144
begin
	using Random
	Random.seed!()

	passLength, numMatchesForLog = 8, 4
	possibleChars = ['a':'z' ; 'A':'Z' ; '0':'9']

	correctPassword = "3xsfe6p3"

	numMatch(loginPassword) = sum([loginPassword[i] == correctPassword[i] for i in 1:passLength])

	Nₚ = 10^7

	passwords = [String(rand(possibleChars, passLength)) for _ in 1:Nₚ]
	numLogs = sum([numMatch(p) >= numMatchesForLog for p in passwords])
end

# ╔═╡ 81248f92-347b-11eb-046d-1327fefd71be
begin
	using StatsBase, Combinatorics, Plots; gr()
	matchExists1(n) = 1 - prod([k/365 for k in 365:-1:365-n+1])
	matchExists2(n) = 1 - factorial(365, 365 - big(n))/365^big(n)
	function bdEvent(n)
		birthdays = rand(1:365, n)
		dayCounts = counts(birthdays, 1:365)
		return maximum(dayCounts) > 1
	end

	probEst(n) = sum([bdEvent(n) for _ in 1:Nᵦ]) / Nᵦ

	xGrid = 1:50
	analyticSolution1 = [matchExists1(n) for n in xGrid]
	analyticSolution2 = [matchExists2(n) for n in xGrid]

	Nᵦ = 10^3
	mcEstimates = [probEst(n) for n in xGrid]
	plot(xGrid, analyticSolution1, c=:blue, label="Analytic solution")
	scatter!(xGrid, mcEstimates, c=:red, ms=4, msw=1, shape=:xcross, label="MC estimate", xlims=(0, 50), ylim=(0, 1), xlabel="Number of people in room", ylabel="Probaility of birthday match", leg=:topleft)

end

# ╔═╡ fcc139c4-3481-11eb-3765-6fe44d8b9951
begin
	using LaTeXStrings

	nₗ, Nₗ = 5, 10^5

	function isUpperLattice(v)
		for i in 1:Int(length(v)/2)
			sum(v[1:2*i-1]) >= i ? continue : return false
		end
		return true
	end

	omega = unique(permutations([zeros(Int, nₗ);ones(Int, nₗ)]))
	A = omega[isUpperLattice.(omega)]
	pA_modelI = length(A)/length(omega)

	function randomWalkPath(n)
		x, y = 0, 0
		path = []
		while x<n && y<n
			if rand() < 0.5
				x += 1
				push!(path, 0)
			else
				y += 1
				push!(path, 1)
			end
		end
		append!(path, x<n ? zeros(Int64, n-x) : ones(Int64, n-y))
		return path
	end

	PA_modelIIest = sum([isUpperLattice(randomWalkPath(nₗ)) for _ in 1:Nₗ]) / Nₗ

	function plotPath(v, l, c)
		x, y = 0, 0
		graphX, graphY = [x], [y]
		for i in v
			if i == 0
				x += 1
			else
				y += 1
			end
			push!(graphX, x), push!(graphY, y)
		end
		plot!(graphX, graphY,
			la=0.8, lw = 2, label=l, c=c, ratio = :equal, legend=:topleft,xlims=(0,nₗ), ylims=(0,nₗ), xlabel="East rightarrow", ylabel="North rightarrow")
	end
	plot()
	plotPath(rand(A), "Upper lattice path", :blue)
	plotPath(rand(setdiff(omega, A)), "Non-Upper lattice path", :red)
	plot!([0, nₗ], [0,nₗ],ls=:dash,c=:black, label="")

end

# ╔═╡ 2462b3a2-3478-11eb-2047-d17d9419ef0a
md"# _Probaility_"

# ╔═╡ 5235dc84-3478-11eb-2f71-11b5b2f75b44
md"
1. A sample space
2. A collection of events
3. A probaility measure
"

# ╔═╡ 8e828a70-3478-11eb-2539-f5be394fbbdf
md"## Event sum of dice"

# ╔═╡ a8ab5e90-3478-11eb-3f35-339109ff8c54
begin
	Nₜ, faces = 10^6, 1:6
	numSol = sum([iseven(i+j) for i in faces, j in faces]) / length(faces)^2
	mcest = sum([iseven(rand(faces) + rand(faces)) for i in 1:Nₜ]) / Nₜ
	println("Numerical solution = $numSol \n Monte Carlo estimate = $mcest")
end

# ╔═╡ 7c6ca8f6-3479-11eb-3f7c-d35794372677
numSol

# ╔═╡ 826c806e-3479-11eb-022d-53c8acda5a74
mcest

# ╔═╡ c95a38d6-3479-11eb-2aed-d9ce01229a1f
md"## Password matching"

# ╔═╡ 5b1c6338-347b-11eb-043a-49549cd35c32
md"## The Birthday Problem"

# ╔═╡ 592036ba-347e-11eb-07da-1134d98b2a0f
md"## Fishing with and without replacement"

# ╔═╡ 5f033766-347f-11eb-3dce-9d840202929b
begin
	function proportionFished(gF, sF, n, N, withReplacement = false)
		function fishing()
			fishInPond = [ones(Int64, gF); zeros(Int64, sF)]
			fishCaught = Int64[]

			for fish in 1:n
				fished = rand(fishInPond)
				push!(fishCaught, fished)
				if withReplacement == false
					deleteat!(fishInPond, findfirst(x->x==fished, fishInPond))
				end
			end
			sum(fishCaught)
		end
		simulations = [fishing() for _ in 1:N]
		proportions = counts(simulations, 0:n)/N

		if withReplacement
			plot!(0:n, proportions,line=:stem, marker=:circle, c=:blue, ms=4, msw=1, label="With Replacement", xlabel="n", ylabel="Probaility", ylims=(0, 0.6))
		else
			plot!(0:n, proportions, line=:stem, marker=:xcross, c=:red, msw=4, label="Without replacement")
		end
	end

	N₄ = 10^6
	goldFish, silverFish, n = 3, 4, 3
	plot()
	proportionFished(goldFish, silverFish, n, N₄)
	proportionFished(goldFish, silverFish, n, N₄, true)
end

# ╔═╡ ac4a60ee-347f-11eb-2580-659176e69bb2
md"## Lattice paths"

# ╔═╡ abacb3da-347f-11eb-06fd-df1a2f6812b6

# ╔═╡ ab179624-347f-11eb-141e-815098dddcd0

# ╔═╡ Cell order:
# ╟─2462b3a2-3478-11eb-2047-d17d9419ef0a
# ╟─5235dc84-3478-11eb-2f71-11b5b2f75b44
# ╟─8e828a70-3478-11eb-2539-f5be394fbbdf
# ╠═a8ab5e90-3478-11eb-3f35-339109ff8c54
# ╠═7c6ca8f6-3479-11eb-3f7c-d35794372677
# ╠═826c806e-3479-11eb-022d-53c8acda5a74
# ╟─c95a38d6-3479-11eb-2aed-d9ce01229a1f
# ╠═d627801e-3479-11eb-1fab-1dfaf4af1144
# ╟─5b1c6338-347b-11eb-043a-49549cd35c32
# ╠═81248f92-347b-11eb-046d-1327fefd71be
# ╟─592036ba-347e-11eb-07da-1134d98b2a0f
# ╠═5f033766-347f-11eb-3dce-9d840202929b
# ╟─ac4a60ee-347f-11eb-2580-659176e69bb2
# ╠═fcc139c4-3481-11eb-3765-6fe44d8b9951
# ╠═abacb3da-347f-11eb-06fd-df1a2f6812b6
# ╠═ab179624-347f-11eb-141e-815098dddcd0

max_stars_repo_name: ZhouZhuofei/juliaDL

max_issues_repo_name: ZhouZhuofei/juliaDL

max_forks_repo_name: ZhouZhuofei/juliaDL