text:
OUTNAME = "instrdef.jl"

instdef = open(readall, ARGS[1])

TERMS  = {}
BINARY = {}
MEMORY = {}
CAST   = {}
OTHER  = {}

parse_line(line, name) = begin
    t = split(lstrip(lstrip(line, collect("$name(")), [' ', '(']), ',')
    (int(t[1]), rstrip(lstrip(t[2])), rstrip(split(lstrip(t[3]), ')')[1]))
end

for line in split(chomp(instdef), '\n')
    if beginswith(line, "HANDLE_TERM_INST")
        push!(TERMS, parse_line(line, "HANDLE_TERM_INST"))
    end
    if beginswith(line, "HANDLE_BINARY_INST")
        push!(BINARY, parse_line(line, "HANDLE_BINARY_INST"))
    end
    if beginswith(line, "HANDLE_MEMORY_INST")
        push!(MEMORY, parse_line(line, "HANDLE_MEMORY_INST"))
    end
    if beginswith(line, "HANDLE_CAST_INST")
        push!(CAST, parse_line(line, "HANDLE_CAST_INST"))
    end
    if beginswith(line, "HANDLE_OTHER_INST")
        push!(OTHER, parse_line(line, "HANDLE_OTHER_INST"))
    end
end

io = open(OUTNAME, "w")
write(io, """
immutable InstructionDef{T}
    opcode::Cint
    capi_name::String
    capi_classname::String
end

""")

write(io, """
const InstructionDefs = [
""")

writeinst(typ, t) = write(io, "\t$(t[1]) => InstructionDef{:$typ}($t...),\n")

map(t -> writeinst(:Term,   t), TERMS)
map(t -> writeinst(:Binary, t), BINARY)
map(t -> writeinst(:Memory, t), MEMORY)
map(t -> writeinst(:Cast,   t), CAST)
map(t -> writeinst(:Other,  t), OTHER)

write(io, "]\n")
close(io)


text:
"""
`Zoom <: ImageTransformation`

Description
============

Multiplies the image height and image width equally by some
constant factor. This means that the size of the output image
depends on the size of the input image. If one wants to resize
each dimension by a different factor, use :class:`Scale` instead.

Usage
======

    Zoom()

    Zoom(factor = 2)

Arguments
==========

- **`factor`** : The constant factor that should be used to scale
both width and height of the output image

Methods
========

- **`transform`** : Applies the transformation to the given Image
or set of images

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
========

    using Augmentor
    using TestImages

    # load an example image
    img = testimage("lena")

    # Zoom the image to twice its size
    img_new = transform(Zoom(2), img)

see also
=========

`ImageTransformation`, `ProbableOperation`, `transform`
"""
immutable Zoom <: ImageTransformation
    factor::Float64

    function Zoom(factor::Real)
        factor > 0 || throw(ArgumentError("factor has to be greater than 0"))
        new(Float64(factor))
    end
end

Zoom(; factor = 2) = Zoom(factor)

Base.showcompact(io::IO, op::Zoom) = print(io, "Zoom by $(op.factor)x.")
multiplier(::Zoom) = 1

function transform{T<:AbstractImage}(op::Zoom, img::T)
    x, y = size(img)
    xn = ceil(Int, x * op.factor)
    yn = ceil(Int, y * op.factor)
    result = copyproperties(img, Images.imresize(img, (xn, yn)))
    _log_operation!(result, op)
end

ext: jl


text:
function ifeast!(A::AbstractMatrix, X₀::AbstractMatrix, nodes::Integer, iter::Integer;
                c=complex(0.0,0.0), r=1.0, debug=false, ϵ=0.05)
    N, m₀ = size(X₀)
    if size(A, 1) != size(A, 2)
        error("Incorrect dimensions of A, must be square")
    elseif size(A,1) != N
        error("Incorrect dimensions of X₀, must match A")
    end
    X = deepcopy(X₀)
    Λ = zeros(m₀)
    res = zeros(m₀)
    resvecs = zeros(ComplexF64, N, m₀)
    temp = complex.(X)
    Aq, Bq, Xq = zeros(ComplexF64, m₀, m₀), zeros(ComplexF64, m₀, m₀), zeros(ComplexF64, m₀, m₀)
    θ = LinRange(π/nodes, 2*π-π/nodes, nodes)
    Q = similar(X, ComplexF64)
    ZmA = similar(A, ComplexF64)
    for nit=1:iter
        Q .= 0.0
        for i=1:nodes
            ZmA .= 0.0
                for j=1:N
                    ZmA[j,j] = r*exp(θ[i]*im)+c
                end
            ZmA .-= A
            for m=1:m₀
                temp[:, m] .= bicgstabl(ZmA, X[:, m])
            end
            Q .+= temp .* exp(θ[i]*im)./nodes
        end
        # Q .= Matrix(qr(Q).Q)
        mul!(resvecs, A, Q) ## why does this one allocate?
        mul!(Aq, Q', resvecs) ### Aq = Q' * A * Q
        mul!(Bq, Q', Q) ### Bq = Q' * Q
        F = eigen!(Aq, Bq)
        Λ = F.values
        Xq = F.vectors
        # Λ, Xq = eigen!(Aq, Bq) ### eigen! allocates
        # for i=1:m₀
        #     Xq[:,i] = Xq[:,i] / vecnorm(Xq)
        # end
        mul!(X, Q, Xq) ### X = Q*Xq
        for i=1:m₀
            vn = norm(@view X[:,i]) ## allocates
            for j=1:N
                X[j,i] /= vn
            end
            # res[i] = vecnorm(A*X[:,i]-X[:,i]*Λ[i])/vecnorm(X[:,i]*Λ[i])
        end
        mul!(resvecs, A, X) ## allocates
        mul!(temp, X, Diagonal(Λ)) ## Diagonal allocates
        resvecs .-= temp ### resvecs = A*X-X*Diagonal(Λ)
        for i in 1:m₀
            res[i]=norm(@view resvecs[:,i])
        end
        if debug println(nit) end
    end
    Λ, X, res
end


text:
using ASEPotential
using DFTK
using Random
using Test
using Unitful
using UnitfulAtomic

@testset "ASEPotential.jl" begin

end


text:
mutable struct Legend
    annotation::Vector{Annotation}
    handle::Vector{Handle}
end

Legend( ; annotation, handle) = Legend(annotation, handle)


text:
"""
Restricted Two-body Model tests.
"""
module R2BPTests

using LinearAlgebra, Unitful, GeneralAstrodynamics, Test

@testset verbose=false "R2BP Determination" begin

    rᵢ = [0.0, 11681.0, 0.0] * u"km"
    vᵢ = [5.134, 4.226, 2.787] * u"km/s"
    corbit = Orbit(CartesianState(rᵢ, vᵢ), Earth)

    @test all(
        keplerian(corbit) .≈ (
            0.723452708202361,
            24509.265399338536u"km",
            151.50460766373865u"°",
            90.0u"°",
            270.0034742609256u"°",
            89.99652573907436u"°"
        )
    )

    @test CartesianState(cartesian(KeplerianState(keplerian(corbit)...), massparameter(Earth))...) ≈ state(corbit)

    korbit = Orbit(
        KeplerianState(
            0.723452708202361,
            24509.265399338536u"km",
            151.50460766373865u"°",
            90.0u"°",
            270.0034742609256u"°",
            89.99652573907436u"°"
        ),
        Earth
    )

    @test all(
        upreferred.(statevector(KeplerianState(keplerian(cartesian(korbit)..., massparameter(Earth))...))) .≈
        upreferred.(statevector(KeplerianState(keplerian(corbit)...)))
    )

end

@testset verbose=false "Kepler's Algorithm" begin

    @testset "Unitful" begin

        rᵢ = [0.0, 11681.0, 0.0] * u"km"
        vᵢ = [5.134, 4.226, 2.787] * u"km/s"
        corbit = Orbit(CartesianState(rᵢ, vᵢ), Earth)

        @test all(upreferred.(statevector(state(kepler(corbit)))) .≈ upreferred.(statevector(state(corbit))))

    end

end

@testset verbose=true "Lambert Solvers" begin

    @testset "Universal" begin

        rᵢ = [0.0, 11681.0, 0.0]u"km"
        vᵢ = [5.134, 4.226, 2.787]u"km/s"
        initial = Orbit(CartesianState(rᵢ, vᵢ), Earth)

        Δt = 1000u"s"
        final = kepler(initial, Δt; tol=1e-12)

        v₁, v₂ = lambert_universal(position(initial), position(final), massparameter(Earth), Δt; trajectory=:short, tolerance=1e-6, max_iter=1000)

        @test all(v₁ .≈ velocity(initial))
        @test all(v₂ .≈ velocity(final))

    end

    @testset "Oldenhuis" begin

        rᵢ = [0.0, 11681.0, 0.0]u"km"
        vᵢ = [5.134, 4.226, 2.787]u"km/s"
        initial = Orbit(CartesianState(rᵢ, vᵢ), Earth)

        Δt = 1000u"s"
        final = kepler(initial, Δt; tol=1e-12)

        m = 0
        v₁, v₂ = lambert(position(initial), position(final), Δt, m, massparameter(Earth))

        @test_skip all(v₁ .≈ velocity(initial))
        @test_skip all(v₂ .≈ velocity(final))

    end

end

end # module


text:
export randErdosRenyi

function randErdosRenyi(n::Int, p = 0.5::Float64 )
    randHermitian(Bernoulli(p),n)
end


text:
"""
Manages the build of a project. Includes multiple documents
"""
type MarbleEnv{T}
    settings::SettingsBundle
    cache::States # Storage of ached files
    documents::Array{T}
    # log

    MarbleEnv(settings, cache, documents) = new(settings, cache, documents)
end

# Construct a MarbleEnv with an array of MarbleDocs
# Oh the price of explicit typeing
function MarbleEnv(settings)
    cache = States("$(settings["cachedir"])/filehash_list.json")
    documents = map(settings["documents"]) do docname
        show(docname)
        return MarbleDoc(
            docname[1:findlast(docname, '.') - 1],
            readstring(docname),
            cache,
            SettingsBundle(settings))
    end

    this = MarbleEnv{MarbleDoc}(
        settings,
        cache,
        documents)
    return this
end

function Base.parse(env::MarbleEnv)
    for doc in env.documents
        parse(doc)
    end
end

function process(env::MarbleEnv)
    for doc in env.documents
        process(doc)
    end
end

function render(env::MarbleEnv)
    for doc in env.documents
        render(doc)
    end
end

function template(env::MarbleEnv)
    for doc in env.documents
        template(doc)
    end
end

function build(env::MarbleEnv)
    for doc in env.documents
        build(doc)
    end
end

function buildall(env::MarbleEnv)
    parse(env)
    process(env)
    render(env)
    template(env)
    build(env)
end


text:
# Create object for making move (Rock, Paper, Scissors)
# One perk of this structure is that we can add new moves with slight changes to this struct anytime!
# e.g. Rock, Paper, Scissors, Spock, Lizard and other derivatives of Rock Paper Scissors
struct Thing
    name::String
    wins_with::String
end

# Create possible move objects
scissors = Thing("scissors", "paper")
rock = Thing("rock", "scissors")
paper = Thing("paper", "rock")

# Define winner check function
check_winner(player1::Thing, player2::Thing) = (player1.name == player2.name) ? 0 : (player1.wins_with == player2.name ?  1 : 2)

# Function to change string to a move object
function to_thing(str)
    arr = [scissors, rock, paper]
    return arr[findfirst(x -> x.name == str, arr)]
end

# Function to play_again
function play_again()
    (ask_play_again("Play again? (y/n): ") == "y") ? run(`julia $PROGRAM_FILE`) : exit()
end

# Input part of play_again
function ask_play_again(prompt)
    print(prompt)
    txt = split(lowercase(readline()))
    ((length(txt) > 1) || !(txt[1] in ["y", "n"])) ? ask_play_again(prompt) : return txt[1]
end

# Get input from user, and scrub input
function get_input(prompt)
    print(prompt)
    txt = split(lowercase(readline()))
    ((length(txt) > 1) || !(txt[1] in ["rock", "scissors", "paper"])) ? get_input(prompt) : return txt[1]
end

# Create points for players
AI_points = 0
human_points = 0

while (AI_points + human_points) < 3
    global AI_points, human_points
    # Ask for input and get move by user
    player1_move = to_thing(get_input("Enter your move (Rock, Paper or Scissors): "))

    # Make AI make random move
    player2_move = rand([rock, scissors, paper])

    # Show user the move the AI made
    println("AI: ", player2_move.name)

    # Check winner
    result = check_winner(player1_move, player2_move)

    # Print result. If tie, then game continues, since game counter is sum of AI_points and human_points, the game goes on one more time
    (result == 1) ? (println("You win"); human_points += 1) : (result == 2) ? (println("AI wins"); AI_points += 1) : (println("Tie. No one gets a point. Play again.\n"); continue)

    # Print current standings
    println("Current points - You: $human_points, AI: $AI_points\n")

    # If any player gets 2 points out of 3, playing game again is pointless, so let the player with 2 points win.
    max(AI_points, human_points) == 2 ? break : ""
end

# Print game result
(AI_points > human_points) ? println("\nGAME! AI Wins.") : println("\nGAME! You Win.")

# Ask if user wants to play again
play_again()

max_stars_repo_name: soumitradev/RockPaperScissors

max_issues_repo_name: soumitradev/RockPaperScissors

max_forks_repo_name: soumitradev/RockPaperScissors


text:
@everywhere function SimulateBattle(ArmiesAttacker::Int64, ArmiesDefender::Int64)
  while ArmiesDefender > 0 && ArmiesAttacker > 0
    DiceAttacker = rand(1:6, min(3, ArmiesAttacker))
    DiceDefender = rand(1:6, min(2, ArmiesDefender))
    sort!(DiceAttacker, rev=true)
    sort!(DiceDefender, rev=true)
    if DiceAttacker[1] > DiceDefender[1]
      ArmiesDefender -= 1
    else
    ArmiesAttacker -= 1
    end
    if ArmiesDefender >1 && ArmiesAttacker>1
      if DiceAttacker[2] > DiceDefender[2]
        ArmiesDefender -=1
      else ArmiesAttacker -=1
    end
  end
  end
  (ArmiesAttacker, ArmiesDefender)
end

@everywhere function SimulateBattleFast(ArmiesAttacker::Int64, ArmiesDefender::Int64)
  while ArmiesDefender > 2 && ArmiesAttacker > 3
    a = rand()
    if a <= 0.372
      ArmiesDefender -= 2
    elseif a > 0.372 && a <= 0.707
      ArmiesDefender -= 1
      ArmiesAttacker -= 1
    else
      ArmiesAttacker -= 2
    end
  end

  while ArmiesDefender > 0 && ArmiesAttacker > 0
    DiceAttacker = rand(1:6, min(3, ArmiesAttacker))
    DiceDefender = rand(1:6, min(2, ArmiesDefender))
    sort!(DiceAttacker, rev=true)
    sort!(DiceDefender, rev=true)
    if DiceAttacker[1] > DiceDefender[1]
      ArmiesDefender -= 1
    else
      ArmiesAttacker -= 1
    end
    if ArmiesDefender > 1 && ArmiesAttacker > 1
      if DiceAttacker[2] > DiceDefender[2]
        ArmiesDefender -= 1
      else ArmiesAttacker -= 1
    end
  end
  end
  (ArmiesAttacker, ArmiesDefender)
end

function SimulateRiskio(ArmiesAttacker, ArmiesDefender, NSim::Int64)
  result = Array{Float32}(length(ArmiesAttacker)
                          , length(ArmiesDefender))
    for (k, AA) in enumerate(ArmiesAttacker)
      for (j, AD) in enumerate(ArmiesDefender)
       battles = SharedArray(Int64,NSim,3)
        @sync @parallel for i in 1:NSim
           b = SimulateBattleFast(AA,AD)
           battles[i,1] = b[1]
           battles[i,2] = b[2]
        end
      AttackerWins = battles[:,1] .> battles[:,2]
      result[k,j] = mean(AttackerWins)
    end
  end
  result
end

d = @time SimulateRiskio(1:30,1:30,10000)

Pkg.add("Gadfly")

using Gadfly
using Colors
using Compose

draw(SVG("d:/probs.svg",15cm, 15cm),spy(d))

spy(d)

max_stars_repo_path: Risiko Simu.jl

max_stars_repo_name: wolfganglederer/risiko

max_issues_repo_path: Risiko Simu.jl

max_issues_repo_name: wolfganglederer/risiko

max_forks_repo_path: Risiko Simu.jl

max_forks_repo_name: wolfganglederer/risiko


text:
using BandedMatrices, Base.Test

for T in (Float64,Complex128,Float32,Complex64)
    A=brand(T,10,10,3,2)
    Q,R=qr(A)
    @test_approx_eq full(Q)*full(R) A
    b=rand(T,10)
    @test_approx_eq A_mul_B!(similar(b),Q,Ac_mul_B!(similar(b),Q,b)) b
    for j=1:size(A,2)
        @test_approx_eq Ac_mul_B(Q,A[:,j]) R[:,j]
    end
    A=brand(T,14,10,3,2)

    Q,R=qr(A)

    for k=1:size(Q,1),j=1:size(Q,2)
        @test_approx_eq Q[k,j] full(Q)[k,j]
    end

    @test_approx_eq full(Q)*full(R) A

    A=brand(T,10,14,3,2)

    Q,R=qr(A)

    for k=1:size(Q,1),j=1:size(Q,2)
        @test_approx_eq Q[k,j] full(Q)[k,j]
    end

    @test_approx_eq full(Q)*full(R) A
    A=brand(T,100,100,3,4)
    Q,R=qr(A)
    b=rand(T,100)
    @test_approx_eq R\(Q'*b) full(A)\b
end

A=brand(10,10,3,2)
b=rand(Complex128,10)
Q,R=qr(A)
@test_approx_eq R\(Q'*b) full(A)\b

A=brand(Complex128,10,10,3,2)
b=rand(10)
Q,R=qr(A)
@test_approx_eq R\(Q'*b) full(A)\b

println("Large dimension QR solve")
A=brand(100000,100000,3,4)
b=rand(size(A,1))
@time Q,R=qr(A)
@time R\(Q'*b)

println("Time should be   0.030149 seconds (23 allocations: 17.548 MB)")
println("Time should be   0.008210 seconds (29 allocations: 8.012 MB)")

A=brand(2000,2000,400,400)
b=rand(size(A,1))
@time Q,R=qr(A)
@time R\(Q'*b)

println("Time should be   0.401671 seconds (24 allocations: 36.665 MB)")
println("Time should be   0.009276 seconds (30 allocations: 12.262 MB)")


text:
"""
    plot(c::AbstractChains, ptype::Vector{Symbol}=[:trace, :density];
       <keyword arguments>)::Array{Plots.Plot}

Function that takes a MCMC chain and creates various different plots (trace &
density by default).

# Arguments
- 'vars::Vector{String}=String[]' : specifies the variables of the chain that
                                    are plotted
- 'filename::String=""' : when given, the plots will be saved to a file
- 'fmt::Symbol=:pdf' : Format of the output file
- 'args...': This includes specific arguments for the different plot types
             , like the number of bins for the contourplot or if the barplots
             bars should be stacked or not. Check the specific plot functions
             to use these arguments.

            Plot attributes from the Plots.jl package can also be passed here:
            I.e. try passing background_color=:black for a black background on
            your plots. List of attributes here,
            https://docs.juliaplots.org/latest/generated/attributes_series/
            https://docs.juliaplots.org/latest/generated/attributes_plot/
            https://docs.juliaplots.org/latest/generated/attributes_subplot/
            https://docs.juliaplots.org/latest/generated/attributes_axis/
            though not all are supported:
            https://docs.juliaplots.org/latest/generated/supported/
"""
function plot(c::AbstractChains, ptype::Vector{Symbol}=[:trace, :density];
              vars::Vector{String}=String[], filename::String="",
              fmt::Symbol=:pdf, fuse::Bool=false, f_layout=nothing,
              fsize::Tuple{Number, Number}=(0, 0), args...
              )::Union{Vector{Plots.Plot}, Tuple{Vector{Plots.Plot}, Plots.Plot}}
  if !isempty(vars)
    indeces = check_vars(c.names, vars)
  else
    indeces = collect(1:length(c.names))
  end # if / else
  ilength = length(indeces)
  if :contour in ptype && ilength == 1
    filter!(e -> e ≠ :contour, ptype)
    if isempty(ptype)
      throw(ArgumentError("Drawing a contourplot requires at least 2 variables."))
    else
      @warn "Drawing a contourplot requires at least 2 variables. Only drawing other plots"
    end # if/else
  end # if
  n = length(ptype)
  p = Array{Plots.Plot}(undef, n)
  for i in 1:n
    ptype[i] == :bar ? p[i] = bar_int(c, indeces; args...) :
    ptype[i] == :mixeddensity ? p[i] = mixeddensityplot(c, indeces; args...) :
    p[i] = Plots.plot(c, indeces; ptype=ptype[i], size=(ilength * 500, 300), args...)
    if !fuse
      if n != 1
        println("Press ENTER to draw next plot")
        readline(stdin)
      end # if
      display(p[i])
    end # if
  end # for
  if fuse
    isnothing(f_layout) && (f_layout = (n, 1))
    fsize == (0, 0) && (fsize = (ilength * 500, n * 300))
    allplots = Plots.plot(p..., layout=f_layout, size=fsize)
    display(allplots)
    filename != "" && check_filename(filename, fmt, allplots)
    return (p, allplots)
  end # if
  filename != "" && check_filename(filename, fmt, [p])
  return p
end # plot

"""
  check_vars(sim_names::Vector{AbstractString},
             vars::Vector{String})::Vector{Int64}

--- INTERNAL ---
Helper function that returns a list of indeces that correspond to specific
variables. Only those variables are then plotted in the following steps.
"""
function check_vars(sim_names::Vector{AbstractString}, vars::Vector{String})::Vector{Int64}
    names = []
    for var in vars
        if endswith(var, r"\[[0-9]*\]")
            for sim_name in sim_names
                if sim_name == var
                    push!(names, sim_name)
                end # if
            end # for
        else
            for sim_name in sim_names
                if sim_name == var || occursin(Regex(var * "\\[[0-9]+\\]"), sim_name)
                  #(occursin(sim_name,r"\[[0-9]*\]") && sim_name[1 : findfirst("[", sim_name)[1]-1] == var))
                    push!(names, sim_name)
                end # if
            end # for
        end # if / else
    end # for
    unique!(names)
    indeces = []
    for name in names
      index = findfirst(isequal(name), sim_names)
      push!(indeces, index)
    end # for
    sort!(indeces)
    return indeces
end # check_vars

"""
    check_filename(filename, fmt, plots)

--- INTERNAL ---
Helper function that checks if a user-given filename is valid, and saves the
created plots to that file.
"""
function check_filename(filename, fmt, plots)
  if !(fmt in [:pdf, :png, :ps, :svg])
      @warn "The given format is not supported. Use :pdf, :png, :ps, or :svg. Defaulting to pdf."
      fmt = :pdf
  end # if
  for (i, plot) in enumerate(plots)
      i == 1 ? savefig(plot[i], "$filename.$fmt") :
               savefig(plot[i], "$filename_$i.$fmt")
  end # for
end # check_filename

#################### Internal Plot Function ####################
@recipe function f(c::AbstractChains, indeces::Vector{Int64}; ptype=:autocor,
                   maxlag=round(Int, 10 * log10(length(c.range))), bins=100,
                   trim=(0.025, 0.975))
  grid --> :dash
  gridalpha --> 0.5
  legend --> false
  legendtitle --> "Chain"
  legendtitlefonthalign := :left
  margin --> 7mm

  arr = []
  ptype == :autocor ? push!(arr, Autocor(c, indeces, maxlag)) :
  ptype == :contour ? push!(arr, Contour(c, indeces, bins)) :
  ptype == :density ? push!(arr, Density(c, indeces, trim)) :
  ptype == :mean    ? push!(arr, Mean(c, indeces)) :
  ptype == :trace   ? push!(arr, Trace(c, indeces)) :
    throw(ArgumentError("unsupported plot type $ptype"))
  return Tuple(arr)
end # recipe

# structs that are used for the recipes
struct Autocor; c; indeces; maxlag; end
struct Contour; c; indeces; bins; end
struct Density; c; indeces; trim; end
struct Mean; c; indeces; end
struct Trace; c; indeces; end

"""
    f(acor::Autocor)

--- INTERNAL ---
Recipe for Autocor plots
"""
@recipe function f(acor::Autocor)
  xguide --> "Lag"
  yguide --> "Autocorrelation"
  xlims --> (0, +Inf)
  layout --> (1, length(acor.indeces))

  nrows, nvars, nchains = size(acor.c.value)
  lags = 0:acor.maxlag
  ac = autocor(acor.c, lags=collect(lags))
  x = repeat(collect(lags * step(acor.c)), outer=[nchains])
  for (index, i) in enumerate(acor.indeces)
    subplot := index
    y = vec(ac.value[i,:,:])
    ac_group = repeat(acor.c.chains, inner=[length(lags)])
    for (index, chain) in enumerate(acor.c.chains)
      idxs = findall(==(chain), ac_group)
      @series begin
        label --> string(index)
        title --> acor.c.names[i]
        seriestype := :line
        x[idxs], y[idxs]
      end # series
    end # for
  end # for
  primary := false
end # recipe

"""
    bar_int(c::AbstractChains, indeces::Vector{Int64}; args...)::Plots.Plot

--- INTERNAL ---
Helper function for creating barplots.
"""
function bar_int(c::AbstractChains, indeces::Vector{Int64}; args...)::Plots.Plot
  nrows, nvars, nchains = size(c.value)
  ilength = length(indeces)
  bar_plots = Array{Plots.Plot}(undef, ilength)
  for (index, i) in enumerate(indeces)
    S = unique(c.value[:, i, :])
    n = length(S)
    x = repeat(S, 1, nchains)
    y = zeros(n, nchains)
    for j in 1:nchains
      m = StatsBase.countmap(c.value[:, i, j])
      for k in 1:n
        if S[k] in keys(m)
          y[k, j] = m[S[k]] / nrows
        end # if
      end # for
    end # for
    bar_plots[index] = groupbar(vec(x), vec(y);
                                group=repeat(c.chains, inner=[n]),
                                title=c.names[i], args...)
  end # for
  p = Plots.plot(bar_plots..., layout=(1, ilength), size=(ilength * 500, 300))
  return p
end # function

"""
    f(cont::Contour)

--- INTERNAL ---
Recipe for contour plots
"""
@recipe function f(cont::Contour)
  layout --> (1, sum(collect(1:length(cont.indeces) - 1)))
  legend := false
  colorbar --> :best
  subplot_index = 0
  nrows, nvars, nchains = size(cont.c.value)
  offset = 1e4 * eps()
  n = nrows * nchains
  for (index, i) in enumerate(cont.indeces[1:end-1])
    X = cont.c.value[:, i, :]
    qx = range(minimum(X) - offset, stop=maximum(X) + offset, length=cont.bins + 1)
    mx = map(k -> mean([qx[k], qx[k + 1]]), 1:cont.bins)
    idx = Int[findfirst(k -> qx[k] <= x < qx[k + 1], 1:cont.bins) for x in X]
    for j in cont.indeces[index+1:end]
      subplot_index += 1
      Y = cont.c.value[:, j, :]
      qy = range(minimum(Y) - offset, stop=maximum(Y) + offset, length=cont.bins + 1)
      my = map(k -> mean([qy[k], qy[k + 1]]), 1:cont.bins)
      idy = Int[findfirst(k -> qy[k] <= y < qy[k + 1], 1:cont.bins) for y in Y]
      density = zeros(cont.bins, cont.bins)
      for k in 1:n
        density[idx[k], idy[k]] += 1.0 / n
      end
      @series begin
        subplot := subplot_index
        xguide --> cont.c.names[i]
        yguide --> cont.c.names[j]
        colorbar_title --> "Density"
        title --> cont.c.names[i]
        seriestype := :contour
        mx, my, density
      end # series
    end # for
  end # for
  primary := false
end # recipe

"""
    f(dens::Density)

--- INTERNAL ---
Recipe for density plots
"""
@recipe function f(dens::Density)
  xguide --> "Value"
  yguide --> "Density"
  ylims --> (0.0, +Inf)
  layout --> (1, length(dens.indeces))

  nrows, nvars, nchains = size(dens.c.value)
  for (index, i) in enumerate(dens.indeces)
    subplot := index
    val = Array{Vector{Float64}}(undef, nchains)
    dens_group = []
    for j in 1:nchains
      qs = quantile(dens.c.value[:, i, j], [dens.trim[1], dens.trim[2]])
      mask = [qs[1] .<= dens.c.value[:, i, j] .<= qs[2]]
      val[j] = dens.c.value[mask[1], i, j]
      dens_group = vcat(dens_group, repeat([j], inner=sum(mask[1])))
    end # for
    for (index, chain) in enumerate(dens.c.chains)
      idxs = findall(==(chain), dens_group)
      @series begin
        label --> string(index)
        title --> dens.c.names[i]
        seriestype := :density
        [val...;][idxs]
      end # series
    end # for
  end # for
  primary := false
end # recipe

"""
    f(mean::Mean)

--- INTERNAL ---
Recipe for mean plots
"""
@recipe function f(mean::Mean)
  xguide --> "Iteration"
  yguide --> "Mean"
  layout --> (1, length(mean.indeces))

  nrows, nvars, nchains = size(mean.c.value)
  val = cummean(mean.c.value)
  x = repeat(collect(mean.c.range), outer=[nchains])
  for (index, i) in enumerate(mean.indeces)
    subplot := index
    y = vec(val[:, i, :])
    mean_group = repeat(mean.c.chains, inner=[length(mean.c.range)])
    for (index, chain) in enumerate(mean.c.chains)
      idxs = findall(==(chain), mean_group)
      @series begin
        label --> string(index)
        title --> mean.c.names[i]
        seriestype := :line
        x[idxs], y[idxs]
      end # series
    end # for
  end # for
  primary := false
end # recipe

"""
    mixeddensityplot(c::AbstractChains,, indeces::Vector{Int64};
                     barbounds::Tuple{Real, Real}=(0, Inf), args...):Plots.Plot

--- INTERNAL ---
Helper function that creates a barplot for each discrete and a density for each
continuous variable.
"""
function mixeddensityplot(c::AbstractChains, indeces::Vector{Int64};
                          barbounds::Tuple{Real, Real}=(0, Inf), args...)
  plots = Array{Plots.Plot}(undef, length(indeces))
  ilength = length(indeces)
  discrete_temp = MCPhylo.indiscretesupport(c, barbounds)
  discrete = Bool[]
  for index in indeces
    try
      push!(discrete, discrete_temp[index])
    catch BoundsError
    end # try / catch
  end
  for i in 1:length(discrete)
    if discrete[i] == true
      plots[i] = bar_int(c, [indeces[i]]; args...)
    else
      plots[i] = Plots.plot(c, [indeces[i]]; ptype=:density, size=(500, 300), args...)
    end # if / else
  end # for
  p = Plots.plot(plots..., layout=(1, ilength), size=(ilength * 500, 300))
  return p
end # function

"""
    f(trace::Trace)

--- INTERNAL ---
Recipe for trace plots
"""
@recipe function f(trace::Trace)
  xguide --> "Iteration"
  yguide --> "Value"
  layout --> (1, length(trace.indeces))
  legendtitle --> "Chain"
  widen --> false

  nrows, nvars, nchains = size(trace.c.value)
  x = repeat(collect(trace.c.range), outer=[nchains])
  for (index, i) in enumerate(trace.indeces)
    subplot := index
    y = vec(trace.c.value[:, i, :])
    trace_group = repeat(trace.c.chains, inner=[length(trace.c.range)])
    for (index, chain) in enumerate(trace.c.chains)
      idxs = findall(==(chain), trace_group)
      @series begin
        label --> string(index)
        title --> trace.c.names[i]
        seriestype := :line
        x[idxs], y[idxs]
      end # series
    end # for
  end # for
  primary := false
end # recipe

#=
This is a slightly modified version of the GroupedBar recipe found in the
StatsPlots package:
  https://github.com/JuliaPlots/StatsPlots.jl/blob/master/src/bar.jl
=#
@userplot GroupBar
recipetype(::Val{:groupbar}, args...) = GroupBar(args)
Plots.group_as_matrix(g::GroupBar) = true
grouped_xy(x::AbstractVector, y::AbstractArray) = x, y
grouped_xy(y::AbstractArray) = 1:size(y,1), y
"""
    f(g::GroupBar; spacing=0)

--- INTERNAL ---
Recipe for a grouped bar plot.
"""
@recipe function f(g::GroupBar; spacing=0)
    xguide --> "Value"
    yguide --> "Density"
    grid --> :dash
    gridalpha --> 0.5
    legend --> false
    legendtitle --> "Chain"
    legendtitlefonthalign := :left
    margin --> 7mm
    x, y = grouped_xy(g.args...)

    nr, nc = size(y)
    isstack = pop!(plotattributes, :bar_position, :dodge) == :stack
    isylog = pop!(plotattributes, :yscale, :identity) ∈ (:log10, :log)
    ymax = maximum(isstack ? mapslices(sum, y, dims=2) : y)
    ylims --> (0.0, ymax)

    # extract xnums and set default bar width.
    # might need to set xticks as well
    xnums = if eltype(x) <: Number
        xdiff = length(x) > 1 ? mean(diff(x)) : 1
        bar_width --> 0.8 * xdiff
        x
    else
        bar_width --> 0.8
        ux = unique(x)
        xnums = (1:length(ux)) .- 0.5
        xticks --> (xnums, ux)
        xnums
    end
    @assert length(xnums) == nr

    # compute the x centers.  for dodge, make a matrix for each column
    x = if isstack
        x
    else
        bws = plotattributes[:bar_width] / nc
        bar_width := bws * clamp(1 - spacing, 0, 1)
        xmat = zeros(nr,nc)
        for r=1:nr
            bw = Plots._cycle(bws, r)
            farleft = xnums[r] - 0.5 * (bw * nc)
            for c=1:nc
                xmat[r,c] = farleft + 0.5bw + (c-1)*bw
            end
        end
        xmat
    end

    fill_bottom = if isylog
        if isfinite(the_ylims[1])
            min(minimum(y) / 100, the_ylims[1])
        else
            minimum(y) / 100
        end
    else
        0
    end
    # compute fillrange
    y, fr = isstack ? groupedbar_fillrange(y) : (y, get(plotattributes, :fillrange, [fill_bottom]))
    if isylog
        replace!( fr, 0 => fill_bottom )
    end
    fillrange := fr

    seriestype := :bar
    x, y
end

"""
    groupedbar_fillrange(y)

--- INTERNAL ---
Helper function for the groupbar plot
"""
function groupedbar_fillrange(y)
    nr, nc = size(y)
    # bar series fills from y[nr, nc] to fr[nr, nc], y .>= fr
    fr = zeros(nr, nc)
    y = copy(y)
    y[.!isfinite.(y)] .= 0
    for r = 1:nr
        y_neg = 0
        # upper & lower bounds for positive bar
        y_pos = sum([e for e in y[r, :] if e > 0])
        # division subtract towards 0
        for c = 1:nc
            el = y[r, c]
            if el >= 0
                y[r, c] = y_pos
                y_pos -= el
                fr[r, c] = y_pos
            else
                fr[r, c] = y_neg
                y_neg += el
                y[r, c] = y_neg
            end
        end
    end
    y, fr
end


text:
#--
# Threads.nthreads() displays number of available threads
# if ==1, start julia with JULIA_NUM_THREADS=4 julia
#--

#--
# If no graphics (X-windows) required, then ENV["GKSwstype"]="nul"
#--

using FDModBiotCyl
#
using CPUTime
using ProgressMeter
using Plots
using Printf
#

#--Suppressing graphics
#ENV["GKSwstype"]="nul"

@__DIR__
@__FILE__

#================
Drawing a model
================#
function drawmodel(Field_Var,nz,dz,nr,dr,src_index,index_allrec,nrec,LPML_r,LPML_z)
   println("Drawing a model...")
   zvec=[0:dz:(nz-1)*dz]
   rvec=[0:dr:(nr-1)*dr]
   R=(nr-1)*dr
#   plt1=heatmap(rvec,zvec,Vp,yflip=true,xlabel="R (m)",ylabel="Z (m)",ratio=1)
   plt1=heatmap(rvec,zvec,Field_Var,yflip=true,xlabel="R (m)",ylabel="Z (m)")
   plot!([(src_index[2]-1)*dr],[(src_index[1]-1)*dz],label="",markershape=:circle,markersize=5,markerstrokewidth=0,seriescolor=:red)
   plot!((index_allrec[:,2]-ones(nrec))*dr,(index_allrec[:,1]-ones(nrec))*dz,label="",markershape=:cross,markersize=5,markerstrokewidth=0,seriescolor=:red)
   plot!([0;R],ones(2)*(LPML_z-1)*dz,label="",color=:black)
   plot!([0;R],ones(2)*((nz-LPML_z+1)*dz),label="",color=:black)
   plot!(ones(2)*((nr-LPML_r+1)*dr),[0;(nz-1)*dz],label="",color=:black)
   display(plot(plt1))
   tmp_filename=@sprintf "tmp_FDmodel.png"
   tmp_filename_full=string(@__DIR__,"/",tmp_filename)
   png(tmp_filename_full)
   println("--------")
   println("Model config figure is saved in:", tmp_filename_full)
   println("--------")
end

#================
Drawing a snapshot
================#
# Four snapshots(2,2)
function drawsnap_2by2(Data1,Data2,Data3,Data4,
         nz,dz,nr,dr,src_index,index_allrec_tii,nrec,LPML_r,LPML_z,
         title_str1,title_str2,title_str3,title_str4)
   println("Drawing a snapshot...")
   zvec=[0:dz:(nz-1)*dz]
   rvec=[0:dr:(nr-1)*dr]
   R=(nr-1)*dr
   plt1=heatmap(rvec,zvec,Data1,yflip=true,xlabel="R (m)",ylabel="Z (m)",legend=:none)
   heatmap!(title=title_str1)
   plot!([0;R],ones(2)*(LPML_z-1)*dz,label="",color=:black)
   plot!([0;R],ones(2)*((nz-LPML_z+1)*dz),label="",color=:black)
   plot!(ones(2)*((nr-LPML_r+1)*dr),[0;(nz-1)*dz],label="",color=:black)

   plt2=heatmap(rvec,zvec,Data2,yflip=true,xlabel="R (m)",ylabel="Z (m)",ylims=(50,60),legend=:none)
   heatmap!(title=title_str2)
   plot!([0;R],ones(2)*(LPML_z-1)*dz,label="",color=:black)
   plot!([0;R],ones(2)*((nz-LPML_z+1)*dz),label="",color=:black)
   plot!(ones(2)*((nr-LPML_r+1)*dr),[0;(nz-1)*dz],label="",color=:black)

   plt3=heatmap(rvec,zvec,Data3,yflip=true,xlabel="R (m)",ylabel="Z (m)",ylims=(50,60),legend=:none)
   heatmap!(title=title_str3)
   plot!([0;R],ones(2)*(LPML_z-1)*dz,label="",color=:black)
   plot!([0;R],ones(2)*((nz-LPML_z+1)*dz),label="",color=:black)
   plot!(ones(2)*((nr-LPML_r+1)*dr),[0;(nz-1)*dz],label="",color=:black)

   plt4=heatmap(rvec,zvec,Data4,yflip=true,xlabel="R (m)",ylabel="Z (m)",ylims=(50,60),legend=:none)
   heatmap!(title=title_str4)
   plot!([0;R],ones(2)*(LPML_z-1)*dz,label="",color=:black)
   plot!([0;R],ones(2)*((nz-LPML_z+1)*dz),label="",color=:black)
   plot!(ones(2)*((nr-LPML_r+1)*dr),[0;(nz-1)*dz],label="",color=:black)

   display(plot(plt1,plt2,plt3,plt4,layout=(2,2)))

end

#==========================================
Creating a model for a borehole embedded in
 a poroelastic three-layer medium
==========================================#
function makemodel_3L()
   #Necessary input matrices to the main loop function are:
   #Rho: Bulk density
   #Rhof: Fluid density
   #M,C,H: Poroelastic moduli
   #G: Formation shear modulus
   #D1,D2: Poroelastic moduli in Ou's formulation
   #Flag_AC,Flag_E: flag specifying acoustic region or elastic region
   #Other input parameters to the main loop function are:
   #ir_wall: grid number in the r direction at which a borehole wall starts (single value)
   #
   # This function creates the above input parameters assuming a borehole with a constant radius
   # embedded in a three-layer medium.
   # The medium is a thin poroelastic layer sandwiched between two elastic layers

   #Model size
   nr=301 #samples
   nz=701 #samples
   dr=0.01 #meter
   dz=0.2 #meter

   #==============================================
   Initializing grain/frame parameters
   ==============================================#
   #First creating grain/frame parameters and then converting them
   #to necessary poroelastic parameters for our FD
   #----solid phase-----
   Km=zeros(nz,nr) #frame
   Ks=zeros(nz,nr) #grain
   G=zeros(nz,nr) #bulk
   Rho=zeros(nz,nr) #bulk (weighted average of Rhos and Rhof)
   Rhos=zeros(nz,nr) #grain
   #---fluid phase----
   Kf=zeros(nz,nr)
   Rhof=zeros(nz,nr)
   Kappa0=zeros(nz,nr) #static permeability (m^2, 1D=9.869E-13m^2)
   Eta=zeros(nz,nr) #dynamic viscosity Pa.s=kg/(m.s)
   #---other parameters---
   Phi=zeros(nz,nr)
   Tot=zeros(nz,nr) #Tortuosity factor (see Sidler 2014)
   #End initializing poroelastic parameters

   #==============================================
   Fluid parameters
   ==============================================#
   #Assuming that fluid properties are identical everywhere
   Kf[:,:]=ones(nz,nr)*2.25*10^9
   Rhof[:,:]=ones(nz,nr)*1000
   Kappa0[:,:]=ones(nz,nr)*1*9.869*10^(-13) #m^2
   Eta[:,:]=ones(nz,nr)*1.0*10^(-3) #water 1E-3 Pa.s

   #==============================
    Background poroelastic medium
   ==============================#
   Phi[:,:]=ones(nz,nr)*0.3
   #Tortuosity factor from Sidler
   #Tot=0.5*(ones(nz,nr)+1 ./Phi) #Empirical formula (see Sidler)
   #Tortuosity factor from Ou
   m_Ou=8.0 #
   alpha_inf_Ou=3.0 #tortuosity
   Tot=ones(nz,nr)*(1+2/m_Ou)*alpha_inf_Ou

   #-Solid parameters from dry_Vp, dry_Vs, grain density, and grain bulk modulus
   #Parameters giving sat_Vp=5000, sat_Vs=3000, bulk density=2500
   Vp_dry=5170.0
   Vs_dry=3198.0
   Rhos[:,:]=ones(nz,nr)*3143.0
   Ks[:,:]=ones(nz,nr)*10*1E10
   #

   Rho_dry=(ones(nz,nr)-Phi).*Rhos
   G[:,:]=Vs_dry^2*Rho_dry[:,:]
   Km[:,:]=Vp_dry^2*Rho_dry[:,:]-4/3*G[:,:]
   Rho=(ones(nz,nr)-Phi).*Rhos+Phi.*Rhof

   #========================================
   Inclusion of a borehole (acoustic media)
   ========================================#
   Flag_AC=zeros(nz,nr)
   ir_wall=6 #A3-A5
   for iz=1:nz
       Flag_AC[iz,1:ir_wall]=ones(ir_wall)
       Rho[iz,1:ir_wall]=ones(ir_wall)*Rhof[1,1] #Rho=Rhof
       G[iz,1:ir_wall]=ones(ir_wall)*0.0 #G=0
       Phi[iz,1:ir_wall]=ones(ir_wall) #Phi=1
       Km[iz,1:ir_wall]=ones(ir_wall)*0.0 #Km=0 (G=0,phi=0->M=C=H=Kf, Ou's)
#       Tot[iz,1:ir]=Phi[iz,1:ir] #Tot=Phi (-> D1=rhof, Ou's)
#       Tot[iz,1:ir]=Phi[iz,1:ir] #Tot=Phi (-> D1=1, Guan's)
       Eta[iz,1:ir_wall]=ones(ir_wall)*0.0 #eta=0 (-> D2=0, Ou's)
   end
   #

   #==========================
   Inclusion of Top and Bottom
   Elastic Layers
   ==========================#
   #3Layer model (Elastic-PoroElastic-Elastic)
   #boundary @iz=271, L~1m or k=5 samples
   #
   Flag_E=zeros(nz,nr)
   Vp_elastic=5000.0 #Or, Specify K_elastic
   Vs_elastic=3000.0
   Rho_elastic=2500.0
   #

   for iz=1:270 #Upper Elastic medium

       G_elastic=Vs_elastic^2*Rho_elastic
       K_elastic=Vp_elastic^2*Rho_elastic-4/3*G_elastic #Activate here if you have specified Vp_elastic

       Flag_E[iz,ir_wall+1:end]=ones(nr-ir_wall)
       Rho[iz,ir_wall+1:end]=ones(nr-ir_wall)*Rho_elastic #Rho=Rho_elastic
       Rhof[iz,ir_wall+1:end]=ones(nr-ir_wall)*Rho_elastic #Rhof=Rho_elastic
       G[iz,ir_wall+1:end]=ones(nr-ir_wall)*G_elastic #G=Gs (Ou's)
       #--
       Phi[iz,ir_wall+1:end]=ones(nr-ir_wall)*0.0 #Phi=0
       Km[iz,ir_wall+1:end]=ones(nr-ir_wall)*K_elastic #Km
       Ks[iz,ir_wall+1:end]=ones(nr-ir_wall)*K_elastic #Ks=Km (-> M=Inf, Ou's)
       #--
       Kappa0[iz,ir_wall+1:end]=ones(nr-ir_wall)*0.0 #k0=0 (-> D1=D2=Inf, Ou's)
   end

   for iz=276:nz #Lower Elastic medium

       G_elastic=Vs_elastic^2*Rho_elastic
       K_elastic=Vp_elastic^2*Rho_elastic-4/3*G_elastic #Activate here if you have specified Vp_elastic

       Flag_E[iz,ir_wall+1:end]=ones(nr-ir_wall)
       Rho[iz,ir_wall+1:end]=ones(nr-ir_wall)*Rho_elastic #Rho=Rho_elastic
       Rhof[iz,ir_wall+1:end]=ones(nr-ir_wall)*Rho_elastic #Rhof=Rho_elastic
       G[iz,ir_wall+1:end]=ones(nr-ir_wall)*G_elastic #G=Gs (Ou's)
       #--
       Phi[iz,ir_wall+1:end]=ones(nr-ir_wall)*0.0 #Phi=0
       Km[iz,ir_wall+1:end]=ones(nr-ir_wall)*K_elastic #Km
       Ks[iz,ir_wall+1:end]=ones(nr-ir_wall)*K_elastic #Ks=Km (-> M=Inf, Ou's)
       #--
       Kappa0[iz,ir_wall+1:end]=ones(nr-ir_wall)*0.0 #k0=0 (-> D1=D2=Inf, Ou's)
   end

   #=====================================================
   Converting grain/frame parameters into poroelastic parameters
   =====================================================#
   #Rho, Rhof, M, C, H, G, D1, D2
   #Rho as is
   #Rhof as is
   Alpha=ones(nz,nr)-Km./Ks
   M=( (Alpha-Phi)./Ks+Phi./Kf ).^(-1)
   C=M.*Alpha
   H=Km+4/3*G+M.*Alpha.^2
   #G as is
   D1=Tot.*Rhof./Phi #(1+2/m)*alpha_inf*rhof/phi
   D2=Eta./Kappa0

   #Correction of H and C for elastic media
   for ir=1:nr
      for iz=1:nz
         if(Flag_E[iz,ir]==1)
            H[iz,ir]=Ks[iz,ir]+4/3*G[iz,ir]
            C[iz,ir]=Ks[iz,ir] #Ou's
         end

         if(Flag_AC[iz,ir]==1)
#            D1[iz,ir]=1.0 #Guan's
            D1[iz,ir]=1000. #Ou's D1=rhof
            D2[iz,ir]=0. #Ou's D2=0
         end
      end
   end

   return nr,nz,dr,dz,Rho,Rhof,M,C,H,G,D1,D2,Flag_AC,Flag_E,ir_wall
end

#---function to create receiver geometry for a plane wave (VSP) experiment
function make_receivers()
   #-------------------------
   #  Receiver positions
   #-------------------------
   # Center of the borehole
   nrec1=nz
   recgeom1=zeros(nrec1,2) #(z,r)
   for irec=1:nrec1
      recgeom1[irec,1]=dz*(irec-1)
      recgeom1[irec,2]=0.0
   end

   #Receivers without interpolation
   rec_vr1,rec_vz1,rec_tii1,index_allrec_vr1,index_allrec_vz1,index_allrec_tii1=init_receiver(recgeom1,nrec1,dr,dz,nr,nz,nt)

   # Borehole wall
   # Exact boundrary is at (ir_wall-1)*dr+dr/2.
   # The boundary is the same location as vr and trz.
   # dr/2 inside of solid media: vz, tii, vphi, tpz
   #
   recgeom2=zeros(nrec1,2) #(z,r)
   #homogeneous R version
   for irec=1:nrec1
      recgeom2[irec,1]=dz*(irec-1)
      recgeom2[irec,2]=(ir_wall-1)*dr+dr/2
   end

   #inhomogeneous R version (assuming Vp[1:ir_wall@z]=1500)
#   for irec=1:nrec1
#      ir_wall_now=maximum(findall(x->x==1500,Vp[irec,:])) #assuming nrec=nz
#      recgeom2[irec,1]=dz*(irec-1)
#      recgeom2[irec,2]=(ir_wall_now-1)*dr+dr/2 #this is a position of wall (good for vr)
#   end

   #--Receivers without interpolation
   rec_vr,rec_vz,rec_tii,index_allrec_vr2,index_allrec_vz_dummy,index_allrec_tii_dummy=init_receiver(recgeom2,nrec1,dr,dz,nr,nz,nt)
   # repeat for tii
   for irec=1:nrec1
      recgeom2[irec,1]=dz*(irec-1)
      recgeom2[irec,2]=(ir_wall-1)*dr+dr
   end
#   for irec=1:nrec1
#      ir_wall_now=maximum(findall(x->x==1500,Vp[irec,:])) #assuming nrec=nz
#      recgeom2[irec,1]=dz*(irec-1) #this is a position good for tii (same as Vp)
#      recgeom2[irec,2]=(ir_wall_now-1)*dr+dr #this is a position of wall+dr/2 (good for vz and tii)
#   end
   #--Receivers without interpolation
#   rec_vr,rec_vz,rec_tii,index_allrec_vr_dummy,index_allrec_vz2,index_allrec_tii2=init_receiver(recgeom2,nrec1,dr,dz,nr,nz,nt)
   rec_vr,rec_vz,rec_tii,index_allrec_vr_dummy,index_allrec_vz_dummy,index_allrec_tii2=init_receiver(recgeom2,nrec1,dr,dz,nr,nz,nt)

   # repeat for vz (be careful where to define "solid phase". Vp is defined at tii. vr and vz is located in between.)
   Vp_E=(H./Rho).^0.5
   Vs_E=(G./Rho).^0.5
   for irec=1:nrec1-1
      ir_wall_up=maximum(findall(x->x==1500,Vp_E[irec,:])) #assuming nrec=nz
      ir_wall_down=maximum(findall(x->x==1500,Vp_E[irec+1,:])) #assuming nrec=nz
      if (ir_wall_up==ir_wall_down)
         recgeom2[irec,1]=dz*(irec-1)+dz/2 #this is a position good for vz
         recgeom2[irec,2]=(ir_wall_up-1)*dr+dr #this is a position of wall+dr/2 (good for vz and tii)
      elseif (ir_wall_up>ir_wall_down) #radius gets smaller
         recgeom2[irec,1]=dz*(irec-1)+dz/2 #this is a position good for vz
         recgeom2[irec,2]=(ir_wall_up-1)*dr+dr #
      elseif (ir_wall_up<ir_wall_down) #radius gets larger
         recgeom2[irec,1]=dz*(irec-1)+dz/2 #this is a position good for vz
         recgeom2[irec,2]=(ir_wall_down-1)*dr+dr #
      else
         error("ERROR in RECGEOM!")
      end

   end
   #--Receivers without interpolation
   rec_vr,rec_vz,rec_tii,index_allrec_vr_dummy,index_allrec_vz2,index_allrec_tii_dummy=init_receiver(recgeom2,nrec1,dr,dz,nr,nz,nt)

   #--initializing arrays with total nrec
   nrec=nrec1*2
   rec_vr=zeros(nt,nrec)
   rec_vz=zeros(nt,nrec)
   rec_tii=zeros(nt,nrec)
   #---merge indexes
   index_allrec_vr=zeros(nrec,2) #z,r
   index_allrec_vz=zeros(nrec,2) #z,r
   index_allrec_tii=zeros(nrec,2) #z,r
   index_allrec_vr[1:nrec1,:]=index_allrec_vr1[1:nrec1,:]
   index_allrec_vr[nrec1+1:2*nrec1,:]=index_allrec_vr2[1:nrec1,:]
   index_allrec_vz[1:nrec1,:]=index_allrec_vz1[1:nrec1,:]
   index_allrec_vz[nrec1+1:2*nrec1,:]=index_allrec_vz2[1:nrec1,:]
   index_allrec_tii[1:nrec1,:]=index_allrec_tii1[1:nrec1,:]
   index_allrec_tii[nrec1+1:2*nrec1,:]=index_allrec_tii2[1:nrec1,:]

   index_allrec_vr=map(Int,index_allrec_vr)
   index_allrec_vz=map(Int,index_allrec_vz)
   index_allrec_tii=map(Int,index_allrec_tii)

   #--Receivers vfr and pf----
   nrec_PE=nr
   recgeom_PE=zeros(nrec_PE,2) #(z,r)
   for irec=1:nrec_PE
         recgeom_PE[irec,1]=(273-1)*dz
#      recgeom_PE[irec,1]=(320-1)*dz
      recgeom_PE[irec,2]=dr*(irec-1)
   end
   #--Receivers without interpolation
   rec_vfr_PE,rec_vfz_PE,rec_pf_PE,index_allrec_vfr_PE,index_allrec_vfz_PE,index_allrec_pf_PE=init_receiver(recgeom_PE,nrec_PE,dr,dz,nr,nz,nt)
   rec_tii_PE=zeros(size(rec_pf_PE)) #for mean normal stress

return  nrec,rec_vr,rec_vz,rec_tii,index_allrec_vr,index_allrec_vz,index_allrec_tii,
      nrec_PE,rec_vfr_PE,rec_pf_PE,rec_tii_PE,index_allrec_vfr_PE,index_allrec_pf_PE

end

#----MAIN loop function----
function main_loop!(nr,nz,dr,dz,Rho,Rhof,M,C,H,G,D1,D2,dt,nt,T,
   vr,vz,trr,tpp,tzz,trz,
   vfr,vfz,pf,
   ir_wall,
   Flag_AC,Flag_E,
   nrec,index_allrec_vr,index_allrec_vz,index_allrec_tii,rec_vr,rec_vz,rec_tii,
   nrec_PE,index_allrec_vfr_PE,index_allrec_pf_PE,rec_vfr_PE,rec_pf_PE,rec_tii_PE,
   LPML_r,LPML_z,PML_Wr,PML_Wz,PML_IWr,PML_Wr2,PML_Wz2,PML_IWr2,
   snapshots_vr,snapshots_vz,snapshots_pf,nsnap,itvec_snap,nskip)

#println("start")

Flag_vf_zero=get_Flag_vf_zero(Flag_AC,Flag_E,nr,nz)

#initializing matrices of PML
Prz_T,Pzz_T,Rz_T,Srz_T,PzzPE_T,RzPE_T,
Prz_B,Pzz_B,Rz_B,Srz_B,PzzPE_B,RzPE_B,
Prr_R,Qrp_R,Pzr_R,Qzp_R,Rr_R,Rp_R,Rrz_R,PrrPE_R,RrPE_R,RpPE_R,
Prr_TR,Qrp_TR,Prz_TR,Pzr_TR,Qzp_TR,Pzz_TR,Rr_TR,Rp_TR,Rz_TR,Rrz_TR,Srz_TR,PrrPE_TR,PzzPE_TR,RrPE_TR,RpPE_TR,RzPE_TR,
Prr_BR,Qrp_BR,Prz_BR,Pzr_BR,Qzp_BR,Pzz_BR,Rr_BR,Rp_BR,Rz_BR,Rrz_BR,Srz_BR,PrrPE_BR,PzzPE_BR,RrPE_BR,RpPE_BR,RzPE_BR=
init_PML_Por(LPML_r,LPML_z,nr,nz)

memT_vr,memT_vz,memT_trr,memT_tpp,memT_tzz,memT_trz,memT_vfr,memT_vfz,memT_pf,
memB_vr,memB_vz,memB_trr,memB_tpp,memB_tzz,memB_trz,memB_vfr,memB_vfz,memB_pf,
memR_vr,memR_vz,memR_trr,memR_tpp,memR_tzz,memR_trz,memR_vfr,memR_vfz,memR_pf=
init_memory_variables_Por(LPML_r,LPML_z,nr,nz)
#return

#initializing following matrices (required for BC at the borehole wall)
vr_old=zeros(nz,nr)
vz_old=zeros(nz,nr)
vfr_old=zeros(nz,nr)
vfz_old=zeros(nz,nr)

#to check the Skempton coefficient
mean_tii=zeros(nz,nr)

# Making sure RightBC for stress (zero values are trp and trz)
ApplyBCRight_stress!(vr,vz, #Use with flag_zero=1 (see ApplyBCRight_stress1D01)
    trz,
    G,nr,nz,dr,dz,dt)

@showprogress for ii=1:nt

#----Time at (ii-1)*dt---(updating velocities)-----

#--PML: save velocity at previous step
PML_save_vel!(memT_vr,memT_vz,memT_vfr,memT_vfz,
              memB_vr,memB_vz,memB_vfr,memB_vfz,
              memR_vr,memR_vz,memR_vfr,memR_vfz,
              vr,vz,vfr,vfz,nr,nz,LPML_z,LPML_r)
#println("update vel")

#Copying the current velocity values as _old in order to apply Acoustic-Poroelastic BC later
mycopy_mat(vr,vr_old,nz,nr)
mycopy_mat(vz,vz_old,nz,nr)
mycopy_mat(vfr,vfr_old,nz,nr)
mycopy_mat(vfz,vfz_old,nz,nr)

# Main velocity update!
update_velocity_1st_Por!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,Rho,Rhof,D1,D2,Flag_vf_zero,nr,nz,dr,dz,dt,LPML_z,LPML_r)

#PML: update velocity
PML_update_vel!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
              Rho,Rhof,D1,D2,Flag_vf_zero,nr,nz,dr,dz,dt,
              Prz_T,Pzz_T,PzzPE_T,
              Prz_B,Pzz_B,PzzPE_B,
              Prr_R,Qrp_R,Pzr_R,Qzp_R,Rr_R,Rp_R,Rrz_R,
              PrrPE_R,RrPE_R,RpPE_R,
              Prz_TR,Pzz_TR,PzzPE_TR,
              Prr_TR,Qrp_TR,Pzr_TR,Qzp_TR,Rr_TR,Rp_TR,Rrz_TR,
              PrrPE_TR,RrPE_TR,RpPE_TR,
              Prz_BR,Pzz_BR,PzzPE_BR,
              Prr_BR,Qrp_BR,Pzr_BR,Qzp_BR,Rr_BR,Rp_BR,Rrz_BR,
              PrrPE_BR,RrPE_BR,RpPE_BR,
              LPML_z,LPML_r)

# Making sure RightBC for velocity (zero values are vr)
ApplyBCRight_vel!(vr,vz,
      trr,tpp,tzz,trz,
      vfr,vfz,pf,
      Rho,Rhof,D1,D2,Flag_vf_zero,nr,nz,dr,dz,dt)

#---velocity B.Cs-----
#Periodic Left Edge
ApplyBCLeft_vel!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
                          Rho,Rhof,D1,D2,Flag_vf_zero,
                          nr,nz,dr,dz,dt,
                          Prz_T,Pzz_T,PzzPE_T,
                          Prz_B,Pzz_B,PzzPE_B,
                          LPML_z)

#Boundary conditions at the borehole wall (Fluid and PoroElastic media): replacing velocity values
update_vr_vfr_1st_vertical(vr,trr,tpp,tzz,trz,vfr,pf,
    vr_old,vfr_old,
    Rho,Rhof,D1,D2,nr,nz,dr,dz,dt,LPML_z,
    Prz_T,Prz_B,
    ir_wall,Flag_AC,Flag_E)

#PML: update memory variables for stress (Rx and Sxx) using velocity at two time steps
PML_update_memRS!(Rz_T,Srz_T,RzPE_T,
                  Rz_B,Srz_B,RzPE_B,
                  Rr_R,Rp_R,Rrz_R,RrPE_R,RpPE_R,
                  Rr_TR,Rp_TR,Rz_TR,Rrz_TR,Srz_TR,RrPE_TR,RpPE_TR,RzPE_TR,
                  Rr_BR,Rp_BR,Rz_BR,Rrz_BR,Srz_BR,RrPE_BR,RpPE_BR,RzPE_BR,
                  memT_vr,memT_vz,memT_vfr,memT_vfz,
                  memB_vr,memB_vz,memB_vfr,memB_vfz,
                  memR_vr,memR_vz,memR_vfr,memR_vfz,
                  vr,vz,vfr,vfz,nr,nz,dr,dz,dt,LPML_z,LPML_r,PML_Wz,PML_Wr,PML_IWr,PML_Wz2,PML_Wr2,PML_IWr2)

#--src injection (velocity)
#srcamp=src_func[ii]
#srcapply!(vz,src_index,src_dn,srcamp)
#srcapply!(vr,src_index,src_dn,srcamp)

#----Time at (ii-1)*dt+dt/2----(updating stress)---

#--PML: save stress at previous step
PML_save_stress!(memT_trr,memT_tpp,memT_tzz,memT_trz,memT_pf,
                 memB_trr,memB_tpp,memB_tzz,memB_trz,memB_pf,
                 memR_trr,memR_tpp,memR_tzz,memR_trz,memR_pf,
                 trr,tpp,tzz,trz,pf,nr,nz,LPML_z,LPML_r)

#println("update stress")

# Main stress update!
update_stress_1st_Por!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
     M,C,H,G,nr,nz,dr,dz,dt,LPML_z,LPML_r)
#error()

#PML: stress update
PML_update_stress!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
                  M,C,H,G,nr,nz,dr,dz,dt,
                  Rz_T,Srz_T,RzPE_T,
                  Rz_B,Srz_B,RzPE_B,
                  Rr_R,Rp_R,Rrz_R,RrPE_R,RpPE_R,
                  Rr_TR,Rp_TR,Rz_TR,Rrz_TR,Srz_TR,RrPE_TR,RpPE_TR,RzPE_TR,
                  Rr_BR,Rp_BR,Rz_BR,Rrz_BR,Srz_BR,RrPE_BR,RpPE_BR,RzPE_BR,
                  LPML_z,LPML_r)

# Making sure RightBC for stress (zero values are trp and trz)
ApplyBCRight_stress!(vr,vz, #Use with flag_zero=1 (see ApplyBCRight_stress1D01)
     trz,
     G,nr,nz,dr,dz,dt)

#---stress B.Cs
#println("update stress Left")
ApplyBCLeft_stress!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
                    M,C,H,G,nr,nz,dr,dz,dt,
                    Rz_T,Srz_T,RzPE_T,
                    Rz_B,Srz_B,RzPE_B,
                    LPML_z)

#--src injection (stress:monopole/dipole src)
#srcamp=-src_func[ii]
#srcapply!(pf,src_index,src_dn,-srcamp)
#srcapply!(tzz,src_index,src_dn,srcamp)
#srcapply!(trr,src_index,src_dn,srcamp)
#srcapply!(tpp,src_index,src_dn,srcamp)

#Additional conditions when Flag_Acoustic/Flag_Elastic
ApplyBC_stress_AcousticMedia!(trr,tpp,tzz,trz,pf,Flag_AC,nr,nz) #when Flag_AC==1, then set pf=-(trr+tpp+tzz)/3 and trr=tpp=tzz=-pf
ApplyBC_stress_ElasticMedia_Ou!(pf,Flag_E,nr,nz) #when Flag_E==1, then set pf=0

#--PML: update memory variables for velocity (Pxx and Qxx) using stress at two time steps
PML_update_memPQ!(Prz_T,Pzz_T,PzzPE_T,
                  Prz_B,Pzz_B,PzzPE_B,
                  Prr_R,Qrp_R,Pzr_R,Qzp_R,PrrPE_R,
                  Prr_TR,Qrp_TR,Prz_TR,Pzr_TR,Qzp_TR,Pzz_TR,PrrPE_TR,PzzPE_TR,
                  Prr_BR,Qrp_BR,Prz_BR,Pzr_BR,Qzp_BR,Pzz_BR,PrrPE_BR,PzzPE_BR,
                  memT_trr,memT_tpp,memT_tzz,memT_trz,memT_pf,
                  memB_trr,memB_tpp,memB_tzz,memB_trz,memB_pf,
                  memR_trr,memR_tpp,memR_tzz,memR_trz,memR_pf,
                  trr,tpp,tzz,trz,pf,nr,nz,dr,dz,dt,LPML_z,LPML_r,PML_Wz,PML_Wr,PML_IWr,PML_Wz2,PML_Wr2,PML_IWr2)

#-----End of updating field----

#println("receiver")
#-Receiver field (extraction)
getRecData_from_index!(vr,rec_vr,index_allrec_vr,nrec,ii)
getRecData_from_index!(vz,rec_vz,index_allrec_vz,nrec,ii)
getRecData_from_index!(trr,rec_tii,index_allrec_tii,nrec,ii)

#Additional receiver array along porous layer
getRecData_from_index!(vfr,rec_vfr_PE,index_allrec_vfr_PE,nrec_PE,ii)
getRecData_from_index!(pf,rec_pf_PE,index_allrec_pf_PE,nrec_PE,ii)
mean_tii=(trr+tpp+tzz)/3
getRecData_from_index!(mean_tii,rec_tii_PE,index_allrec_pf_PE,nrec_PE,ii)

#--Snapshots
#println("check snap")
check_snap=findall(x -> x==ii,itvec_snap)
if (length(check_snap)!=0)
   println("ii=",ii)
   cnt_snap=Int(check_snap[1])
   get_snapshots!(snapshots_vz,cnt_snap,vz,nr,nz)
   get_snapshots!(snapshots_vr,cnt_snap,vr,nr,nz)
   get_snapshots!(snapshots_pf,cnt_snap,pf,nr,nz)

   drawsnap_2by2(snapshots_vz[:,:,cnt_snap],snapshots_pf[:,:,cnt_snap],
                 snapshots_vr[:,:,cnt_snap],snapshots_vz[:,:,cnt_snap],nz,dz,nr,dr,
                 src_index,index_allrec_tii,nrec,LPML_r,LPML_z,
                 "Vz","Pf","Vr","Vz")

end

end #i Time

end #function
#------------------------------

#-----ENTRY POINT HERE----
CPUtic()
start=time()

#======================
General variables
======================#

#time samples
dt=0.125E-5
nt=16001
T=(nt-1)*dt
tvec=range(0.0,T,length=nt) # range object (no memory allocation)
tvec=collect(tvec) # a vector

#PML thickness in samples
LPML_r=51
LPML_z=61

#======================
Creating a model
======================#
nr,nz,dr,dz,Rho,Rhof,M,C,H,G,D1,D2,Flag_AC,Flag_E,ir_wall=makemodel_3L()
#error()

#======================
Stability check
======================#
Vmax=maximum((H./Rho).^(0.5)) #for PML
check_stability01(dt,dr,Vmax,0)
check_stability01(dt,dz,Vmax,0)

#======================
Creating a src wavelet
======================#
f0=200 #src Freq
delay=1/f0*1.0
src_func=myricker2(tvec,f0,delay,2) #when using 2nd derivative Gaussian (Ricker wavelet with a negative sign)
tmp_maxamp=maximum(map(abs,src_func))
src_func=src_func/tmp_maxamp
#Note: Peng_solution() assumes unit vz. Therefore, src_func above will be interpreted as vz.
#The following scaling of 1/(rho*Vp) makes the incident wave unit amplitude Ricker wavelet for tzz.
src_func=src_func/(2500*5000)

display(plot(tvec,src_func[:],title="src function"))
tmpfilename_src=string(@__DIR__,"/src_function.png")
png(tmpfilename_src)
println("Source signature figure saved in ",tmpfilename_src)

#================================
Src depth where a plane P wave
starts to propagate downward
================================#
srcdepth=28. #meter

#==============================
Initializing field variables
==============================#
# Initial conditions of a plane P wave
vr,vz,trr,tpp,tzz,trz,
vfr,vfz,pf,
vz_srcBC,vr_srcBC,
trr_srcBC,tpp_srcBC,tzz_srcBC,trz_srcBC,
vz_init,vr_init,trr_init,tpp_init,tzz_init,
trz_init,
src_iz,noffset,
f0=initialize_planewave(nr,nz,dr,dz,ir_wall,Rho,H,G,nt,tvec,dt,src_func,srcdepth,f0,delay)

# Initializing PML field variables
LPML_r,LPML_z,PML_Wr,PML_Wz,
PML_IWr,PML_Wr2,PML_Wz2,PML_IWr2=init_PML_profile(LPML_r,LPML_z,Vmax,dr,dz,nr,f0)
#PML_check(LPML_r,LPML_z,Vmax,dr,dz,f0)

#error()

#==============================
Receiver geometry (plane wave experiments)
==============================#
nrec,rec_vr,rec_vz,rec_tii,index_allrec_vr,index_allrec_vz,index_allrec_tii,
nrec_PE,rec_vfr_PE,rec_pf_PE,rec_tii_PE,index_allrec_vfr_PE,index_allrec_pf_PE=make_receivers()

#==============================
Check model
==============================#
src_index=[1 1]; #dummy
drawmodel(vz_init,nz,dz,nr,dr,src_index,index_allrec_tii,nrec,LPML_r,LPML_z)
#error()

#==============================
Snapshot settings
==============================#
nskip=100 #change here to adjust time sampling (1~nt)
snapshots_vz,nsnap,itvec_snap=init_snapshots(nt,nz,nr,nskip)
snapshots_vr,nsnap,itvec_snap=init_snapshots(nt,nz,nr,nskip)
snapshots_pf,nsnap,itvec_snap=init_snapshots(nt,nz,nr,nskip)

#error()

#==============================
Start main FD Loop
==============================#

main_loop!(nr,nz,dr,dz,Rho,Rhof,M,C,H,G,D1,D2,dt,nt,T,
   vr_init,vz_init,trr_init,tpp_init,tzz_init,trz_init,
   vfr,vfz,pf,
   ir_wall,
   Flag_AC,Flag_E,
   nrec,index_allrec_vr,index_allrec_vz,index_allrec_tii,rec_vr,rec_vz,rec_tii,
   nrec_PE,index_allrec_vfr_PE,index_allrec_pf_PE,rec_vfr_PE,rec_pf_PE,rec_tii_PE,
   LPML_r,LPML_z,PML_Wr,PML_Wz,PML_IWr,PML_Wr2,PML_Wz2,PML_IWr2,
   snapshots_vr,snapshots_vz,snapshots_pf,nsnap,itvec_snap,nskip)

CPUtoc()
println("elapsed real time: ", round(time() - start;digits=3)," seconds")

println("FD finished without a problem.")

println("Plotting receiver responses...")
FD_z0=272*dz #the reference depth of the porous layer
recdepth=(float(index_allrec_tii[1:nz,1])-ones(nz,1))*dz #receiver depth

offset=recdepth-FD_z0*ones(size(recdepth))
plt1=heatmap(offset[:],tvec[1:10:end],-rec_tii[1:10:end,1:nz],xlabel="Depth (m)",ylabel="time (s)",title="pressure")
heatmap!(yflip=true)
display(plot(plt1))
plt2=plot(tvec[:],-rec_tii[:,223],xlabel="Depth (m)",ylabel="time (s)",title="-10.0m",label="",ylims=(-0.03,0.03))
plt3=plot(tvec[:],-rec_tii[:,323],xlabel="Depth (m)",ylabel="time (s)",title="+10.0m",label="",ylims=(-0.03,0.03))
display(plot(plt2,plt3,layout=(2,1)))

rvec=range(0.0,(nr-1)*dr,length=nr) # range object (no memory allocation)
rvec=collect(rvec) # a vector
plt4=plot(rvec[:],rec_pf_PE[4001,:],xlabel="r (m)",ylabel="pressure (Pa)",title="t=0.005 s",label="Pf",ylims=(-0.08,-0.06))
plot!(rvec[:],-rec_tii_PE[4001,:],xlabel="r (m)",ylabel="pressure (Pa)",title="t=0.005 s",label="-tii/3",ylims=(-0.08,-0.06))


text:
module FitBellCurves

using Optim
using ..Utils: gaussian, lorentzian, voigtian, flatten_arrays, array_arrays
using ..FTFSpectrum: _objective_function, MeanAbs

abstract type FitShapeCurves end
struct FitCurveModel{T0, T1} <: FitShapeCurves where {T0<:Symbol, T1<:Float64}
    model::T0
    xdata::Vector{T1}
    ydata::Vector{T1}
    seed::Array{Array{T1,1},1}
    lb::Array{Array{T1,1},1}
    ub::Array{Array{T1,1},1}
    solution
    ymodel::Vector{T1}
    optParams::Array{Array{T1,1},1}
    sigma::Array{T1,1}
end

export FitCurveModel, fit_curve_model

function fit_curve_model(
    model::T0,
    xdata::AbstractVector{T1},
    ydata::AbstractVector{T1},
    seed::Array{Array{T1,1},1};
    lb::Array{Array{T1,1},1}=-2.0.*array_arrays((abs.(flatten_arrays(seed)).+1.0), seed),
    ub::Array{Array{T1,1},1}=2.0.*array_arrays((abs.(flatten_arrays(seed)).+1.0), seed),
    σ::Array{T1,1}=ones(length(ydata)),
    alg=SAMIN(),
    options=Optim.Options(
            g_abstol=1e-8,
            g_reltol=1e-8,
            iterations=10^5,
            show_trace=false,
            store_trace=true,
    ),
) where {T0<:Symbol, T1<:Real}
    isequal(length(xdata), length(ydata)) || throw("Input xdata and ydata must have the same length.")
    oftype = MeanAbs()
    sol = optimize(p -> _fit_mse(p, model, seed, xdata, ydata, σ, oftype),
                   flatten_arrays(lb), flatten_arrays(ub), flatten_arrays(seed), alg, options)
    xopt = array_arrays(sol.minimizer, seed)
    sol_ = FitCurveModel(
                model,
                float.(xdata),
                float.(ydata),
                float.(seed),
                float.(lb),
                float.(ub),
                sol,
                eval(model)(xdata, xopt),
                xopt,
                σ,
    )
    return sol_
end

function _fit_mse(p, model, seed, xdata, ydata, σ, oftype)
    ŷ = eval(model)(xdata, array_arrays(p, seed))
    mse = _objective_function(oftype, vec(ŷ), ydata; σ=σ)
    return mse
end

end # module


text:
function demo_naive_update()
    # Update PCCF after each observed change
    # just by recalculating mean and standard deviation values
    #
    observed_chps_locations = [1, 11+0.9, 21-0.5, 31+1.5, 41-1.2, 51+1.25, 61, 71-1.0]
    #observed_chps_locations = [1, 10, 20, 30, 40, 50, 60, 70]
    n_chp = length(observed_chps_locations)
    mu_0 = 10.0
    sigma_0 = 1.2
    N = 100
    M = zeros(N, n_chp) # matrix to store updated PCCFs
    pccf_init = Pccf.pccf(N, mu_0, sigma_0)
    M[:,1] = pccf_init
    for i in 1:length(observed_chps_locations)
        if i > 1
            subdat = observed_chps_locations[1:i]
            mu = mean(diff(subdat))
            if length(subdat)>2
                sigma = std(diff(subdat))
            else
                sigma = sigma_0
            end
            print("new mu=",mu," new sigma=", sigma, "\n")
            M[:,i] = Pccf.pccf(N, mu, sigma)
        end
    end
   writedlm("D:/MyTemp/tmp_exp_out/upd_pccf.txt", M, ",")
end

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
module Sum_of_Powers

export Solution, err
export search, Best, SubSet, SlidingWindow, MabeyBest
export write_file

import Base.string

include("Solution.jl")
include("OnePositions.jl")
include("track_best.jl")
include("search.jl")
include("write_file.jl")

end  #module


text:
const KB = 1024
const MB = KB*KB
const BLKSZ = 128 * MB
const RECSZ = 100
const RECSPERBLK = round(Int, BLKSZ/RECSZ)

const DFSUSER = ENV["USER"]
const UGI = UserGroupInformation(DFSUSER)
const DFSHOST = "localhost"
const DFSPORT = 9000

const INDIR = "/sortinp"
const OUTDIR = "/sortout"
const SAMPDIR = "/sortsamp"

function cleanup()
    println("cleaning up last run...")
    dfs = HDFSClient(DFSHOST, DFSPORT, UGI)
    for d in (INDIR, OUTDIR, SAMPDIR)
        exists(dfs, d) && rm(dfs, d, true)
    end
    nothing
end

function startcluster(clean::Bool=false)
    if clean
        cleanup()
        println("creating empty folder...")
        dfs = HDFSClient(DFSHOST, DFSPORT, UGI)
        for d in (INDIR, OUTDIR, SAMPDIR)
            mkdir(dfs, d)
        end
    end

    println("creating app master...")
    yarncm = YarnManager(yarnhost="localhost", rmport=8032, schedport=8030, launch_timeout=60)
    nnodes = nodecount(yarncm.clnt)

    println("adding $nnodes nodes")
    appenv = Dict("JULIA_PKGDIR" => Pkg.dir())
    addprocs(yarncm; np=nnodes, env=appenv)
    if nworkers() != nnodes
        stopcluster()
        error("Error bringing up all nodes")
    end
    yarncm
end

function stopcluster(yarncm::YarnManager)
    rmprocs(workers())
    Elly.disconnect(yarncm)
end

function splitpoints(len::Int, n::Int)
    out = Int[]
    i = 0
    nrem = n
    while i < len
        s = ceil(Int, ((len - i) / nrem))
        nrem -= 1
        push!(out, i+=s)
   end
   out[1:(n-1)]
end


text:
using Test
using SparseIR
using Random
using MultiFloats

include("__conftest.jl")
include("_util.jl")

@testset verbose = true "SparseIR.jl" begin
    include("gauss.jl")
    include("kernel.jl")
    include("poly.jl")
    include("sve.jl")
    include("basis_set.jl")
    include("sampling.jl")
    include("matsubara.jl")
    include("augment.jl")
    include("bessels.jl")
    include("composite.jl")
    include("spr.jl")
    include("_linalg.jl")
    # Works only if Mathematica and MathLink.jl are available.
    # include("_bessels.jl")
end

nothing


text:
module GenerativePreTrain

using Flux

using ..Transformers: Abstract3DTensor
using ..Basic

include("./gpt.jl")
include("./load_pretrain.jl")

end


text:
include("newtypes.jl")

# using DataFrames
# importall Base

########################
### TwoWayDict class ###
########################

function bijective_invert{K,V}(d::Dict{K,V})
    result =  Dict{V,K}()
    for (k, v) in d
        result[v] = k
    end
    result
end

type TwoWayDict{K,V}
    first  :: Dict{K, V}
    second :: Dict{V, K}
end

TwoWayDict{K,V}(first::Dict{K,V}) = TwoWayDict(first, bijective_invert(first))
TwoWayDict(pairs::Vararg) = TwoWayDict(Dict(pairs))
getindex{K,V}(d::TwoWayDict{K,V}, k::K) = d.first[k]
getindex{K,V}(d::TwoWayDict{K,V}, v::V) = d.second[v]

#################
### Key class ###
#################

immutable Key
  root    :: ModInt{12}
  ismajor :: Bool
end

const key_dict = TwoWayDict(
    Key(0,  true) => "C",
    Key(1,  true) => "Db",
    Key(2,  true) => "D",
    Key(3,  true) => "Eb",
    Key(4,  true) => "E",
    Key(5,  true) => "F",
    Key(6,  true) => "Gb",
    Key(7,  true) => "G",
    Key(8,  true) => "Ab",
    Key(9,  true) => "A",
    Key(10, true) => "Bb",
    Key(11, true) => "B",
    Key(0, false) => "c",
    Key(1, false) => "c#",
    Key(2, false) => "d",
    Key(3, false) => "eb",
    Key(4, false) => "e",
    Key(5, false) => "f",
    Key(6, false) => "f#",
    Key(7, false) => "g",
    Key(8, false) => "g#",
    Key(9, false) => "a",
    Key(10,false) => "bb",
    Key(11,false) => "b")

Key(root::ModInt{12}, mode::AbstractString) =
    Key(root, mode == "major" ? true : false)
Key(root::Int, mode::AbstractString) =
    Key(ModInt{12}(root), mode == "major" ? true : false)
Key(s::AbstractString) = key_dict[s]

show(io::IO, k::Key) = print(io, key_dict[k])
root(key::Key) = key.root
ismajor(key::Key) = key.ismajor
isminor(k::Key) = !ismajor(k)
mode(k::Key)    = ismajor(k) ? "major" : "minor"
change_mode(k::Key) = Key(root(k), !ismajor(k))

function uses_flats(k::Key)
    if ismajor(k)
        root(k) * 7 > ModInt{12}(5)
    else
        (root(k) + 3) * 7 > ModInt{12}(5)
    end
end

+(k::Key, i) = Key(root(k) + i, ismajor(k))
-(k::Key, i) = Key(root(k) - i, ismajor(k))

###########################
### Scalae degree class ###
###########################

immutable ScaleDeg
    deg :: ModInt{7}
end

const roman_dict = TwoWayDict(
    "I"  => ModInt{7}(0),
    "II" => ModInt{7}(1),
    "III"=> ModInt{7}(2),
    "IV" => ModInt{7}(3),
    "V"  => ModInt{7}(4),
    "VI" => ModInt{7}(5),
    "VII"=> ModInt{7}(6)
)

ScaleDeg(str::String) = ScaleDeg(roman_dict[str])

==(s::ScaleDeg, str::String) = numeral(s) == str
==(str::String, s::ScaleDeg) = s == str

numeral(s::ScaleDeg) = roman_dict[s.deg]
show(io::IO, s::ScaleDeg) = print(io, numeral(s))

+(s::ScaleDeg, i::Int) = ScaleDeg(s.deg + i)
-(s::ScaleDeg, i::Int) = ScaleDeg(s.deg - i)

const major_pattern = (0,2,4,5,7,9,11)
const minor_pattern = (0,2,3,5,7,8,10)

function root(s::ScaleDeg, k::Key) :: ModInt{12}
    pattern = ismajor(k) ? major_pattern : minor_pattern
    pattern[s.deg] + root(k)
end

const natural_tone_dict = TwoWayDict(
    "C" => ModInt{12}(0),
    "D" => ModInt{12}(2),
    "E" => ModInt{12}(4),
    "F" => ModInt{12}(5),
    "G" => ModInt{12}(7),
    "A" => ModInt{12}(9),
    "B" => ModInt{12}(11)
)

function pitch_class(spelled_pitch_class::AbstractString)
    m = match(r"([A-G])(#*)(b*)", spelled_pitch_class)
    natural_tone_dict[String(m[1])] + length(m[2]) - length(m[3])
end

function spelled_root(s::ScaleDeg, k::Key) :: String
    pc = root(s, k) # pitch class of the root
    if pc in major_pattern # then it is the pitch class of a natural tone
        natural_tone_dict[pc]
    elseif uses_flats(k)
        natural_tone_dict[pc+1] * "b"
    else # key uses sharps
        natural_tone_dict[pc-1] * "#"
    end
end

const major_modes = ("major","minor","minor","major","major","minor","minor")
const minor_modes = ("minor","minor","major","minor","minor","major","major")
mode(s::ScaleDeg, k::Key) = (ismajor(k) ? major_modes : minor_modes)[s.deg]

modulate(s::ScaleDeg, k::Key) = Key(root(s, k), mode(s, k))

const major_triad_forms = ("maj","min","min","maj","maj","min","dim")
const minor_triad_forms = ("min","hdim","maj","min","min","maj","maj")
chord_form(s::ScaleDeg, k::Key) = (ismajor(k) ? major_triad_forms : minor_triad_forms)[s.deg]

terminal(s::ScaleDeg, k::Key, form::String=chord_form(s)) = string(Int(root(s, k)), form)

######################
### JazzCategory class ###
######################

abstract type AbstractJazzCategory end

immutable JazzTerminal <: AbstractJazzCategory
    symbol :: String
    root   :: ModInt{12}
    form   :: String
end
==(t1::JazzTerminal, t2::JazzTerminal) = t1.root == t2.root && t1.form == t2.form
hash(t::JazzTerminal, h::UInt) = hash(hash("JazzTerminal", hash(t.root, hash(t.form))), h)

JazzTerminal(s::ScaleDeg, k::Key, form::String=chord_form(s, k)) =
    JazzTerminal(spelled_root(s, k) * form, root(s, k), form)

function JazzTerminal(str::AbstractString)
    if ismatch(r"^(\d+)(.*)", str)
        m = match(r"^(\d*)(.*)", str)
        JazzTerminal(str, parse(m[1]), m[2])
    elseif ismatch(r"^([A-G][b#]*)(.*)", str)
        m = match(r"^([A-G][b#]*)(.*)", str)
        pc = pitch_class(m[1])
        form = m[2]
        JazzTerminal(str, pc, form)
    end
end
show(io::IO, t::JazzTerminal) = print(io, t.symbol)

function make_terminals(scaledegs, key)
   lst = split(scaledegs)
   [JazzTerminal(ScaleDeg(String(s)),Key(key)) for s in lst]
end

immutable JazzCategory <: AbstractJazzCategory
  scaledeg   :: ScaleDeg
  key        :: Key
end
JazzCategory(str::String, k::Key) = JazzCategory(ScaleDeg(str), k)
JazzCategory(deg::ScaleDeg, str::String) = JazzCategory(deg, Key(str))
JazzCategory(str1::String, str2::String) = JazzCategory(ScaleDeg(str1), Key(str2))

show(io::IO, cat::JazzCategory) =
    print(io, scaledeg(cat), "_{", key(cat), "}")

scaledeg(cat::JazzCategory) = cat.scaledeg
key(cat::JazzCategory) = cat.key

keyroot(cat::JazzCategory) = root(key(cat))
keymode(cat::JazzCategory) = mode(key(cat))
root(cat::JazzCategory) = root(scaledeg(cat), key(cat))
spelled_root(cat::JazzCategory) = spelled_root(scaledeg(cat), key(cat))
modulate(cat::JazzCategory) = modulate(scaledeg(cat), key(cat))
change_mode(cat::JazzCategory) = JazzCategory(scaledeg(cat), change_mode(key(cat)))

JazzTerminal(cat::JazzCategory) = JazzTerminal(scaledeg(cat), key(cat))
JazzTerminal(cat::JazzCategory, form::AbstractString) = JazzTerminal(scaledeg(cat), key(cat), form)

# diatonic transposition
Base.transpose(cat::JazzCategory, i) = JazzCategory(scaledeg(cat) + i, key(cat))

const chord_form_dict = Dict(
    ":maj7" => "^7",
    "7" => "7",
    ":min7" => "m7",
    "h7" => "%7",
    "7#9" => "7",
    "7b9" => "7",
    ":maj7#5" => "^7",
    "" => "^",
    "6" => "6",
    "9" => "7",
    ":min6" => "m6",
    "o7" => "o7",
    "+" => "+",
    "7#5" => "7",
    ":maj7#11" => "^7",
    "69" => "6",
    ":min7b5" => "%7",
    "9sus" => "sus",
    "7b9#5" => "7",
    "7b13" => "7",
    "7alt" => "7",
    "7sus" => "sus",
    ":min" => "m",
    "7#9b5" => "7",
    ":min9" => "m",
    "7#11" => "7",
    "13" => "7",
    "h" => "%7",
    "^9" => "^7",
    ":min:maj7" => "m^7",
    "11" => "7",
    ":min11" => "m7",
    "o" => "o7",
    "9#5" => "7",
    "7b9sus" => "sus",
    "7b5" => "7",
    "add9" => "7",
    "7b9b5" => "7",
    "^" => "^",
    "9#11" => "7",
    "maj7#11" => "^7",
    "maj7#5" => "^7",
    "6;" => "6",
    ":maj7;" => "^7",
    ":min6;" => "m6",
    "7#9#5" => "7",
    "13b9" => "7",
    ":minb6" => "m7",
    "7#11;" => "7",
    "13sus" => "sus",
    ":maj7#11;" => "^7",
    ":min69" => "m6",
    "maj7" => "^7",
    "^9#11" => "^7",
    ";" => "^",
    "7;" => "7",
    "13#11" => "7",
    "7b9#9" => "7",
    ":maj9" => "^7",
    "13#11;" => "7",
    "7#9#11" => "7",
    ":maj" => "^",
    "7susadd3" => "sus",
    ":min;" => "m",
    "7b9b13" => "7",
    "7b9#11" => "7",
    "7#9;" => "7",
    ":min:maj7;" => "m^7",
    ":min9;" => "m7",
    ":7" => "7",
    "769" => "7",
    ":7sus" => "sus",
    "dim7" => "o7",
    "sus" => "sus",
    "7b13sus" => "sus",
    "9;" => "7",
    "9#11;" => "7",
    ":maj7#9#11;" => "^7",
    "7#9#5;" => "7",
    "69;" => "6",
    ":maj7#9#11" => "^7",
    "9b5" => "7",
    "^;" => "^"
)

function read_iReal_songs(dir="../iRealPro/iRb_thru/")
    files = readdir(dir)

    pieces_with_meter_change = [
        "conferenceofthebirds.jazz",
        "homecoming.jazz",
        "howmyheartsings.jazz",
        "imeanyou.jazz",
        "isayalittleprayerforyou.jazz",
        "joshua.jazz",
        "lookoflove.jazz",
        "midnightattheoasis.jazz",
        "tellmeabedtimestory.jazz",
        "walktall.jazz"
    ]

    pieces_with_strange_meter = [ # non-strange meters are "4/4", "3/4", and "6/4"
        "litha.jazz",             # 6/8 meter
        "looktotherainbow.jazz",  # 3/2. meter
        "takefive.jazz",          # 5/4
        "thebalance.jazz"         # 5/4
    ]

    iReal_songs = DataFrame(
        file = String[],
        title = String[],
        composers = Vector{String}[],
        year = Int[],
        form = String[],
        meter = String[],
        key = String[],
        chords = DataFrame[]
    )

    for (i, file) in enumerate(setdiff(files, [pieces_with_meter_change; pieces_with_strange_meter]))
        if ismatch(r"^\.", file) continue end # skip hidden files

        lines = readlines(joinpath(dir, file))

        title = ""
        composers = String[]
        year = 0
        form = ""
        meter = ""
        keystring = ""

        bar = 1
        beat = 1.0
        chords = DataFrame(
            bar=Int[],
            beat=Float64[],
            duration=Float64[],
            spelledroot=String[],
            pcroot=Int[],
            chordform=String[],
            slash=Union{Void, String}[],
            alternative=Union{Void, String}[]
        )

        for line in lines
            if ismatch(r"[!*]", line) # meta data line
                if ismatch(r"^!!!OTL", line)
                    @assert title == ""
                    title = match(r"^!!!OTL: (.*)", line)[1]
                elseif ismatch(r"^!!!COM", line)
                    push!(composers, match(r"^!!!COM\d?:?\t? ?(.*)", line)[1])
                elseif ismatch(r"^!!!ODT", line)
                    @assert year == 0
                    year = parse(match(r"^!!!ODT: (.*)", line)[1])
                elseif ismatch(r"^\*>\[", line)
                    @assert form == ""
                    form = match(r"^\*>(.*)", line)[1]
                elseif ismatch(r"^\*M", line)
                    # @assert meter == ""
                    meter = match(r"^\*M(.*)", line)[1]
                elseif ismatch(r"^\*(.*):", line)
                    # @assert keystring == ""
                    keystring = match(r"^\*(.*):", line)[1]
                end
            elseif ismatch(r"^=", line) # new bar
                bar += 1
                beat = 1.0
            elseif ismatch(r"^\dr", line) # rest line
                duration = 4/parse(match(r"^(\d)r", line)[1])
                beat += duration
            elseif ismatch(r"^\d", line) # note line
                @assert ismatch(r"(\d+)(\.*)([A-G])([-#]?)([^/]*)([A-G][-#]?)?(\(.*\))?", line)
                m = match(r"(\d+)(\.*)([A-G])([-#]?)([^/()]*)([A-G][-#]?)?(\(.*\))?", line)
                duration, punktuation, chordroot, accidential, chordform, slash, alternative =
                    4/parse(m[1]), m[2], m[3], replace(m[4], "-", "b"), m[5], m[6], m[7]
                for p in punktuation
                    duration += duration / 2
                end
                pcroot = natural_tone_dict[String(chordroot)] +
                    (contains(accidential, "#") ? length(accidential) : -length(accidential))
                spelledroot = chordroot*accidential
                push!(chords,
                    (bar, beat, duration, spelledroot, pcroot.val, chord_form_dict[chordform], slash, alternative)
                )

                beat += duration
            end
        end
        push!(iReal_songs, (file, title, composers, year, form, meter, keystring, chords))
    end
    iReal_songs
end

iReal_songs = read_iReal_songs();
@show number_chord_symbols = size(vcat(iReal_songs[:chords]...))

function terminals(iReal_chord_df)
    Vector(
        JazzTerminal.(
            iReal_chord_df[:spelledroot] .* iReal_chord_df[:chordform],
            iReal_chord_df[:pcroot],
            iReal_chord_df[:chordform]
        )
    )
end

@show terminals(iReal_songs[4, :chords])[1:13]

const major_keys = [Key(r,"major") for r in 0:11]
const minor_keys = [Key(r,"minor") for r in 0:11]
const scale_degrees = [ScaleDeg(s) for s in 0:6]

categories(scale_deg::String) = categories([scale_deg])

function categories(scale_degrees::Vector{String}=["I", "II", "III", "IV", "V", "VI", "VII"]; without="")
    [
        JazzCategory(s, k)
        for s in ScaleDeg.(setdiff(scale_degrees, [without]))
        for k in vcat(major_keys, minor_keys)
    ]
end

function init_jazz_grammar(dependent_components=cat -> (scaledeg(cat), keymode(cat)))
    prolongation = CFRule(categories(), :prolongation) do c
        [c, c]
    end

    diatonic_prep = CFRule(categories(without="IV"), :diatonic_prep) do c
        [transpose(c, 4), c]
    end
    dominant_prep = CFRule(categories(without="I"), :dominant_prep) do c
        [JazzCategory("V", modulate(c)), c]
    end
    plagal_prep = CFRule(categories("I"), :plagal_prep) do c
        [JazzCategory("IV", key(c)), c]
    end

    modulation = CFRule(categories(without="I"), :modulation) do c
        [JazzCategory("I", modulate(c))]
    end
    mode_change = CFRule(categories("I"), :mode_change) do c
        [change_mode(c)]
    end

    diatonic_subs = CFRule(categories(["I", "II", "V"]), :diatonic_subs) do c
        if ismajor(key(c))
            if scaledeg(c) == "I"
                [JazzCategory("VI", key(c))]
            elseif scaledeg(c) == "II"
                [JazzCategory("IV", key(c))]
            else
                [JazzCategory("VII", key(c))]
            end
        else
            if scaledeg(c) == "I"
                [JazzCategory("III", key(c))]
            elseif scaledeg(c) == "II"
                [JazzCategory("IV", key(c))]
            else
                [JazzCategory("VII", key(c))]
            end
        end
    end

    for i in (3, 6, 9)
        eval(quote
            $(Symbol("octa_dom_subs", i)) = CFRule(categories("V"), Symbol("octa_dom_subs", $i)) do c
                [JazzCategory(scaledeg(c), key(c) + $i)]
            end
        end)
    end

    category_rules = [
        prolongation,
        diatonic_prep,
        dominant_prep,
        plagal_prep,
        modulation,
        mode_change,
        diatonic_subs,
        octa_dom_subs3,
        octa_dom_subs6,
        octa_dom_subs9
    ];

    pitchclass(c::JazzCategory) = root(modulate(c))

    chord_form_dict_ = Dict(
        "major" => Dict(
            "I"   => ["^", "^7", "6"],
            "II"  => ["m", "m7"],
            "III" => ["m", "m7"],
            "IV"  => ["^", "^7", "6"],
            "V"   => ["^", "7", "+", "sus"],
            "VI"  => ["m", "m7"],
            "VII" => ["o7", "%7"]
        ),
        "minor" => Dict(
            "I"   => ["m", "m7", "m^7"],
            "II"  => ["%7"],
            "III" => ["^", "^7", "6"],
            "IV"  => ["m", "m7"],
            "V"   => ["m", "m7", "^", "7", "+", "sus"],
            "VI"  => ["^", "^7", "6"],
            "VII" => ["^", "7", "+", "sus"]
        )
    )

    terminal_rules = push!(
        [
            CFRule([JazzCategory(s, Key(r, mode)) for r in 0:11], Symbol(form, "_termination")) do c
                [JazzTerminal(spelled_root(c) * form, pitchclass(c), form)]
            end
            for mode in ("major", "minor")
            for s in ("I", "II", "III", "IV", "V", "VI", "VII")
            for form in chord_form_dict_[mode][s]
        ],
        CFRule(JazzCategory.("V", major_keys), Symbol("m6_termination")) do c
            [JazzTerminal(spelled_root(transpose(c, 4)) * "m6", pitchclass(transpose(c, 4)), "m6")]
        end
    )

    CFGrammar(category_rules, terminal_rules, categories("I"), dependent_components)
end

iReal_songs[1:10, [:file, :title, :composers, :year, :form, :meter, :key]]

# grammar = init_jazz_grammar()
# chords = terminals(iReal_songs[4, :chords])[1:13]
# @time forest = run_chartparser(chords, grammar, cyclic=false); iscomplete(forest)
# @time forest = run_chartparser(chords, grammar, cyclic=true); iscomplete(forest)
# println(sample_tree(forest))
# iscomplete(forest)
# score(forest)
# category.(heads(forest))
# println(sample_tree(forest))
# @time category.(heads(run_chartparser(chords, grammar)))
#
# terminals(iReal_songs[4, :chords])


text:
export LazyDiscretePost,
       ApproximatingDiscretePost

import LazySets.use_precise_ρ

"""
    LazyDiscretePost <: DiscretePost

Textbook implementation of a discrete post operator, but with lazy intersections.

### Fields

- `options` -- an `Options` structure that holds the algorithm-specific options

### Algorithm

The algorithm is based on [Flowpipe-Guard Intersection for Reachability
Computations with Support Functions](http://spaceex.imag.fr/sites/default/files/frehser_adhs2012.pdf).
"""
struct LazyDiscretePost <: DiscretePost
    options::Options

    function LazyDiscretePost(𝑂::Options)
        𝑂copy = copy(𝑂)
        # TODO: pass 𝑂 directly?
        check_aliases_and_add_default_value!(𝑂.dict, 𝑂copy.dict, [:check_invariant_intersection], false)
        check_aliases_and_add_default_value!(𝑂.dict, 𝑂copy.dict, [:overapproximation], Hyperrectangle)
        check_aliases_and_add_default_value!(𝑂.dict, 𝑂copy.dict, [:lazy_R⋂I], false)
        check_aliases_and_add_default_value!(𝑂.dict, 𝑂copy.dict, [:lazy_R⋂G], true)
        check_aliases_and_add_default_value!(𝑂.dict, 𝑂copy.dict, [:lazy_A⌜R⋂G⌟], "invertible")
        check_aliases_and_add_default_value!(𝑂.dict, 𝑂copy.dict, [:lazy_A⌜R⋂G⌟⋂I], true)
        check_aliases_and_add_default_value!(𝑂.dict, 𝑂copy.dict, [:combine_invariant_guard], 𝑂copy[:lazy_R⋂I])

        if 𝑂copy[:combine_invariant_guard] && !𝑂copy[:lazy_R⋂I]
            throw(ArgumentError("option :combine_invariant_guard only makes " *
                                "sense in combination with option :lazy_R⋂I"))
        end

        return new(𝑂copy)
    end
end

# convenience constructor from pairs of symbols
LazyDiscretePost(𝑂::Pair{Symbol,<:Any}...) = LazyDiscretePost(Options(Dict{Symbol,Any}(𝑂)))

# default options for the LazyDiscretePost discrete post operator
LazyDiscretePost() = LazyDiscretePost(Options())

"""
    ApproximatingDiscretePost()

Textbook implementation of a discrete post operator, but with lazy intersections
followed by an overapproximation. This is a particular case of the
`LazyDiscretePost`.
"""
function ApproximatingDiscretePost()
    return LazyDiscretePost(:check_invariant_intersection=>false,
                            :overapproximation=>Hyperrectangle,
                            :lazy_R⋂I=>false,
                            :lazy_R⋂G=>false,
                            :lazy_A⌜R⋂G⌟⋂I=>false)
end

function ApproximatingDiscretePost(𝑂::Options)
    𝑂_default = Options(:lazy_R⋂I=>false,
                        :lazy_R⋂G=>false,
                        :lazy_A⌜R⋂G⌟⋂I=>false)
    merge!(𝑂_default, 𝑂)
    LazyDiscretePost(𝑂_default)
end

init(𝒫::LazyDiscretePost, 𝒮::AbstractSystem, 𝑂::Options) = init!(𝒫, 𝒮, copy(𝑂))

# TODO: use 𝑂 only?
function init!(𝒫::LazyDiscretePost, 𝒮::AbstractSystem, 𝑂::Options)
    𝑂[:n] = statedim(𝒮, 1)

    # solver-specific options (adds default values for unspecified options)
    𝑂out = validate_solver_options_and_add_default_values!(𝑂)

    return 𝑂out
end

function tube⋂inv(𝒫::LazyDiscretePost,
                  reach_tube::Vector{<:AbstractReachSet{<:LazySet{N}}},
                  invariant,
                  start_interval
                 ) where {N}

    dirs = 𝒫.options[:overapproximation]

    Rsets = Vector{AbstractReachSet{<:LazySet{N}}}()
    @inbounds for reach_set in reach_tube
        R⋂I = Intersection(set(reach_set), invariant)
        if 𝒫.options[:check_invariant_intersection] && isempty(R⋂I)
            break
        end
        if !𝒫.options[:lazy_R⋂I]
            R⋂I = overapproximate(R⋂I, dirs)
        end
        push!(Rsets,
              substitute(reach_set, set=R⋂I,
                         time_start=time_start(reach_set) + start_interval[1],
                         time_end=time_end(reach_set) + start_interval[2]))
    end

    return Rsets
end

function post(𝒫::LazyDiscretePost,
              HS::HybridSystem,
              waiting_list::Vector{Tuple{Int, <:AbstractReachSet{LazySet{N}}, Int}},
              passed_list,
              source_loc_id,
              tube⋂inv,
              jumps,
              options
             ) where {N}
    jumps += 1
    oa = 𝒫.options[:overapproximation]
    source_invariant = HS.modes[source_loc_id].X
    inv_isa_Hrep, inv_isa_H_polytope = get_Hrep_info(source_invariant)

    for trans in out_transitions(HS, source_loc_id)
        info("Considering transition: $trans")
        target_loc_id = target(HS, trans)
        target_loc = HS.modes[target(HS, trans)]
        target_invariant = target_loc.X
        constrained_map = resetmap(HS, trans)
        guard = stateset(constrained_map)

        guard_isa_Hrep, guard_isa_H_polytope = get_Hrep_info(guard)
        combine_constraints = 𝒫.options[:combine_invariant_guard] &&
                              inv_isa_Hrep && guard_isa_Hrep
        if combine_constraints # combine the constraints of invariant and guard
            T = inv_isa_H_polytope || guard_isa_H_polytope ? HPolytope : HPolyhedron
            # TODO: remove redundant constraints => use intersection(..)
            invariant_guard = T([constraints_list(source_invariant);
                                 constraints_list(guard)])
        end

        # perform jumps
        post_jump = Vector{ReachSet{LazySet{N}}}()
        sizehint!(post_jump, length(tube⋂inv))
        for reach_set in tube⋂inv
            # check intersection with guard
            if combine_constraints
                R⋂G = Intersection(set(reach_set).X, invariant_guard)
            else
                R⋂G = Intersection(set(reach_set), guard)
            end
            if isempty(R⋂G)
                continue
            end
            if !𝒫.options[:lazy_R⋂G]
                R⋂G = overapproximate(R⋂G, oa)
            end

            # apply assignment
            A⌜R⋂G⌟ = apply_assignment(𝒫, constrained_map, R⋂G)
            if 𝒫.options[:lazy_A⌜R⋂G⌟] == "always" ||
                    (𝒫.options[:lazy_A⌜R⋂G⌟] == "invertible" &&
                     !isinvertible(constrained_map))
                A⌜R⋂G⌟ = overapproximate(A⌜R⋂G⌟, oa)
            end

            # intersect with target invariant
            A⌜R⋂G⌟⋂I = Intersection(target_invariant, A⌜R⋂G⌟)
            if isempty(A⌜R⋂G⌟⋂I)
                continue
            end
            if !𝒫.options[:lazy_A⌜R⋂G⌟⋂I]
                A⌜R⋂G⌟⋂I = overapproximate(A⌜R⋂G⌟⋂I, oa)
            end

            # store result
            push!(post_jump, ReachSet{LazySet{N}}(A⌜R⋂G⌟⋂I,
                                                  time_start(reach_set),
                                                  time_end(reach_set)))
        end

        postprocess(𝒫, HS, post_jump, options, waiting_list, passed_list,
            target_loc_id, jumps)
    end
end

function get_Hrep_info(set::LazySet)
    return (false, false)
end

function get_Hrep_info(set::AbstractPolytope)
    return (true, true)
end

function get_Hrep_info(set::AbstractPolyhedron)
    return (true, false)
end

# --- line-search policies ---

# usually do not use line search
function use_precise_ρ(𝒫::LazyDiscretePost,
                       cap::Intersection{N})::Bool where N<:Real
    return false
end

# use line search for the outermost level, which is a LinearMap
function use_precise_ρ(𝒫::LazyDiscretePost,
                       cap::Intersection{N, <:LinearMap{N}}
                       )::Bool where N<:Real
    return true
end