text:
# This file contains a script to convert `*.ipynb` files
# to Markdown files to further transfer to `Documenter`

@info "ConvertScript: converting notebooks to markdown files."

# Get path to the notebooks folder
notebooks_folder_name = "notebooks"
notebooks_folder = normpath(joinpath(@__DIR__, "..", notebooks_folder_name))

# Get paths to notebooks and their (base)names
names = String[]
notebooks = String[]
preambles = String[]
for (root, dirs, files) in walkdir(notebooks_folder)
    for file in files
        if endswith(file, ".ipynb")
            name = replace(basename(file), ".ipynb" => "")
            push!(notebooks, joinpath(root, file))
            push!(names, name)
            contains(name, "KIC") && push!(preambles, name)
        end
    end
end

# Get path to the `make.jl` script
make = joinpath(@__DIR__, "make.jl")

# Determine the Markdown lines sets
md_lines_preamble = Dict{Int, String}(
    [
        1 => """
        Let's add a couple of packages. If you try this out, know that they take up a
        pretty decent amount of memory.
        """,

        2 => """
        Now let's get the data from the
        [Kepler Archive](https://archive.stsci.edu/kepler/data_search/search.php).
        To make it easier, let's use the Python's [`kplr`](https://github.com/dfm/kplr)
        package:
        """,

        3 => """
        It's useful to look at this data, but before that let's make the plots a little
        prettier. It's never harmful.
        """,

        4 => """
        Okay, now all the power of TeX is under our control. Let's make that plot then:
        """,

        5 => """
        These sharp leaps are associated with the inconsistency of observational periods.
        Let's pick one of the quarters:
        """,

        6 => """
        That will be enough. Let's look at the time series:
        """,

        7 => """
        Gaussian processes are cool, of course, but the models they will use will require
        initial values. Some of these have already been determined before (mean and
        variance), but we will also need an estimate of the period. Let's go the classic
        way and make this estimate by determining the peak position on the Lomb-Scargle
        periodogram.
        """,

        8 => """
        Here's another look at it, just with a logarithmic scale:
        """,

        9 => """
        Yeah, close the plot. It's just good practice. This ends the preamble to this
        object, and we can now proceed with tests of different models for Gaussian
        processes. Click on any of the model IDs in the sidebar to continue. If you're not
        sure what they mean, check out the [`IDs`](@ref IDs) page.
        """,
    ]
)

md_lines_kernel = Dict{Int, String}(
    [
        1 => """
        Let's execute the preamble in the current scope. Text output will go here as a
        brief reminder of the input data.
        """,

        2 => """
        Let's define a function for calculating the negative log marginal likelihood and an
        auxiliary function for unpacking the tuple of parameters. See the [`IDs`](@ref IDs)
        page on the sidebar for decrypting model IDs.
        """,

        3 => """
        Initial values can be quite important in the optimization process, so it is
        advisable to set them well. Fortunately, a small manual selection with creating the
        plot can help in this. Let's try to choose such parameters that a realization
        of the Gaussian process is similar to the original time series.
        """,

        4 => """
        Let's optimize the negative log marginal likelihood function. This process
        can take quite a lot of time, so the output of the optimizer is very large and
        therefore placed under the spoiler. Don't mind the failure, by the way. It's
        because of the prohibition of the BFGS solver from getting out of the pits.
        """,

        5 => """
        Now let's create a plot similar to the previous one using a set of parameters
        obtained during the optimization process.
        """,

        6 => """
        We can now compare this realization of the Gaussian process with the original
        time series for similarity of characteristics. This is for a test of matching our
        subjective opinion of similarity with the objective result of the solver.
        """,

        7 => """
        This ends the notebook about this model. There are others on the sidebar.
        """,
    ]
)

# Go to the output directory
generated_folder = joinpath(@__DIR__, "src", "generated")
output_folder = joinpath(generated_folder, notebooks_folder_name)
!isdir(generated_folder) && mkdir(generated_folder)
!isdir(output_folder) && mkdir(output_folder)
cd(output_folder)

for (index, notebook) in enumerate(notebooks)

    # Get the name of the notebook
    name = names[index]

    # Initialize an empty Markdown lines set
    md_lines = Dict{Int, String}()

    # A flag showing whether current notebook is about a kernel
    about_kernel = false

    # Determine the preamble
    for preamble in preambles
        if name == preamble

            # Create a working directory
            !isdir(name) && mkdir(name)
            cd(name)

            # Choose a Markdown lines set
            md_lines = md_lines_preamble

            break

        elseif count(preamble, notebook) == 1

            # Create a working directory
            !isdir(preamble) && mkdir(preamble)
            cd(preamble)
            !isdir(name) && mkdir(name)
            cd(name)

            # Choose a Markdown lines set
            md_lines = md_lines_kernel

            about_kernel = true

            break

        end
    end

    # Construct the name of a Markdown file
    md = "$name.md"

    # Convert the notebook to a Markdown file
    output_cmd_part = ["--output-dir", ".", "--output", "$md", "--log-level", "WARN"]
    run(`jupyter nbconvert --to markdown $notebook $output_cmd_part`)

    # Get the content of the generated Markdown file
    lines = readlines("$md")

    # Initialize auxiliary variables
    inside_julia_block = false
    inside_output_block = false
    inside_interrupt_exception = false
    last_index = size(lines, 1)

    # Initialize the code block counter
    code_block = 0

    # Put the text output in the text blocks (and more!)
    for (index, line) in enumerate(lines)

        # Condition of being inside an interrupt exception block
        if inside_interrupt_exception && startswith(line, "    ")

            lines[index] = ""

        # Condition of entering a code block
        elseif startswith(line, "```julia") && !inside_julia_block

            code_block += 1

            # Condition of exiting an output block (code block)
            if inside_output_block

                # Condition of placing the end of a spoiler (before a block of code)
                if about_kernel && code_block == 5
                    lines[index - 1] = """
                    ```
                    ```@raw html
                    </details><br>
                    ```
                    """
                else
                    lines[index - 1] = "```\n"
                end

                inside_output_block = false

            end

            # Insert a Markdown line if it is defined for the current cell
            lines[index] = get(md_lines, code_block, "") * '\n' * lines[index]

            inside_interrupt_exception = false
            inside_julia_block = true

        # Condition of exiting a code block
        elseif startswith(line, "```") && inside_julia_block

            inside_julia_block = false

        # Condition of entering a block of interrupt exception
        elseif !inside_interrupt_exception && inside_output_block #=
            =# && line == "    InterruptException:"

            # Condition of placing the end of a spoiler (before an interrupt exception)
            if about_kernel && code_block == 4
                lines[index - 1] = """
                ```
                ```@raw html

                ```
                """
            else
                lines[index - 1] = "```\n"
            end

            lines[index] = ""

            inside_output_block = false
            inside_interrupt_exception = true

        # Condition of exiting an output block (Markdown line)
        elseif inside_output_block && !isempty(line) && !startswith(line, "    ")

            # Condition of placing the end of a spoiler (before a Markdown line)
            if about_kernel && code_block == 4
                lines[index - 1] = """
                ```
                ```@raw html

                ```
                """
            else
                lines[index - 1] = "```\n"
            end

            inside_output_block = false

        # Condition of entering an output block
        elseif !inside_output_block && !inside_julia_block && !inside_interrupt_exception #=
            =# && startswith(line, "    ")

            # Condition of placing the start of a spoiler
            if about_kernel && code_block == 4
                lines[index - 1] = """
                ```@raw html
                <details>
                <summary>Optimizer output</summary>
                ```
                ```text
                """
            else
                lines[index - 1] = "\n```text\n"
            end

            inside_output_block = true

        # Condition of getting the last line of the last block of output
        elseif inside_output_block && index == last_index

            lines[index] = "\n```"

        end

    end

    # Construct a meta block to hide the `Edit on GitHub` URL
    meta_block = join(
        [
            "```@meta",
            "EditURL = \"nothing\"",
            "```\n\n",
        ],
        '\n',
    )

    # Write the meta block
    open("$md", "w") do io
        lines[1] = meta_block * lines[1]
        print(io, join(lines, '\n'))
    end

    # Go back to the output folder
    cd(output_folder)

end

# Go back to the current directory
cd(@__DIR__)

max_stars_repo_name: paveloom-c/GPKernels

max_issues_repo_name: paveloom-c/GPKernels

max_forks_repo_name: paveloom-c/GPKernels


text:
module CovarianceRealism

using   SIMDPirates, SLEEFPirates, SLEEFwrap,
        LoopVectorization,
        StaticArrays,# BandedMatrices,
        Random,
        Distributions,
        Base.Cartesian,
        SpecialFunctions,
        UnsafeArrays,
        Statistics,
        LinearAlgebra,
        ScatteredArrays,
        VectorizationBase,
        KernelDensityDistributionEsimates, KernelDensity,
        RandomNumbers, VectorizedRNG,
        DifferentiableObjects,
        ForwardDiff,
        PaddedMatrices,
        Base.Threads
        # Gadfly # plotting

export  process_BPP!,
        run_sample!,
        thread_sample!,
        WorkDataChains,
        MCMCResult, MCMCResults,
        WeightedSamples,
        UniformSamples,
        ResizableMatrix,
        InverseWishart, CholInvWishart, RevCholWishart,
        MahalanobisDistances, PercentileMatch,
        sample_distances!,
        KDE,
        PercentileMatch!,
        mixture_fit,
        pc2dfoster_RIC,
        InvCholCovar,
        sample_Pc!

using PaddedMatrices: MutableFixedSizeVector,
    MutableFixedSizeMatrix

# RandomNumbers.jl exports an AbstractRNG{T} that is a subset of AbstractRNG.
# We want our AbstractRNGs to refer to the more general Random.AbstractRNG.
import Random: AbstractRNG

set_zero_subnormals(true)

struct PCG_Scalar_and_Vector <: AbstractRNG
    scalar::RandomNumbers.PCG.PCGStateUnique{UInt64,Val{:RXS_M_XS},UInt64}
    vector::VectorizedRNG.PCG{2}
end
@inline Random.rand(pcg::PCG_Scalar_and_Vector) = rand(pcg.scalar)
@inline function Random.rand(pcg::PCG_Scalar_and_Vector, T::Union{Float16,Float32,Float64})
    rand(pcg.scalar, T)
end
@inline function Random.rand(pcg::PCG_Scalar_and_Vector, A::AbstractArray)
    rand(pcg.scalar, A)
end
@inline function Random.randn(pcg::PCG_Scalar_and_Vector)
    randn(pcg.scalar)
end
@inline function Random.randn(pcg::PCG_Scalar_and_Vector, T::Union{Float16,Float32,Float64})
    randn(pcg.scalar, T)
end
@inline function Random.randexp(pcg::PCG_Scalar_and_Vector)
    randexp(pcg.scalar)
end
@inline function Random.randexp(pcg::PCG_Scalar_and_Vector, T::Union{Float16,Float32,Float64})
    randexp(pcg.scalar, T)
end
@inline function Random.rand(pcg::PCG_Scalar_and_Vector, ::Type{NTuple{N,Core.VecElement{T}}}) where {N,T}
    rand(pcg.vector, NTuple{N,Core.VecElement{T}})
end
@inline function Random.randn(pcg::PCG_Scalar_and_Vector, ::Type{NTuple{N,Core.VecElement{T}}}) where {N,T}
    randn(pcg.vector, NTuple{N,Core.VecElement{T}})
end
@inline function Random.randexp(pcg::PCG_Scalar_and_Vector, ::Type{NTuple{N,Core.VecElement{T}}}) where {N,T}
    randexp(pcg.vector, NTuple{N,Core.VecElement{T}})
end
@inline function Random.rand(pcg::PCG_Scalar_and_Vector, ::Type{VectorizationBase.SVec{N,T}}) where {N,T}
    SVec(rand(pcg.vector, NTuple{N,Core.VecElement{T}}))
end
@inline function Random.randn(pcg::PCG_Scalar_and_Vector, ::Type{VectorizationBase.SVec{N,T}}) where {N,T}
    SVec(randn(pcg.vector, NTuple{N,Core.VecElement{T}}))
end
@inline function Random.randexp(pcg::PCG_Scalar_and_Vector, ::Type{VectorizationBase.SVec{N,T}}) where {N,T}
    SVec(randexp(pcg.vector, NTuple{N,Core.VecElement{T}}))
end

const PCG_Scalar = Union{PCG_Scalar_and_Vector, RandomNumbers.PCG.PCGStateUnique}
const PCG_Vector = Union{PCG_Scalar_and_Vector, VectorizedRNG.PCG}

const GLOBAL_PCG = PCG_Scalar_and_Vector(
    RandomNumbers.PCG.PCGStateUnique(PCG.PCG_RXS_M_XS),
    VectorizedRNG.PCG{2}()
)

include("utilities.jl")
include("distance_samples.jl")
include("ar_process.jl")
# include("kernel_density_estimate.jl")
# include("gaussian_process.jl")

include("mahalanobis_distances.jl")
include("misc_linear_algebra.jl") # old code, needs updating.

include("wisharts.jl")
include("groups.jl")
include("rng.jl")

include("mcmc_result.jl")
include("working_data.jl")

include("process_inputs.jl")
include("gibbs.jl")

include("simplex.jl")
include("mixture.jl")
include("percentile_matching.jl")

# include("gaussian_process.jl")
include("RIC2ECI.jl")
include("pc_2d_foster.jl")

end # module

ext: jl


text:
################################## Platform ###########################################
"""
    struct Platform{M,C}
      models::Dict{Symbol,M}     # dictionary storing Models
      scenarios::Dict{Symbol,C} # dictionary storing Scenarios
    end

The main storage representing a modeling platform.
Typically HetaSimulator works with one platform object which can include several models and scenarios.

Usually a `Platform` is created based on Heta formatted files using [`load_platform`]{@ref}.

To get the platform content use methods: `models(platform)`, `scenarios(platform).
"""
struct Platform{M,C}
  models::OrderedDict{Symbol,M}
  scenarios::OrderedDict{Symbol,C}
end

models(p::Platform) = p.models
scenarios(p::Platform) = p.scenarios

function Base.show(io::IO, mime::MIME"text/plain", p::Platform)
  models_names = join(keys(p.models), ", ")
  scn_names = join(keys(p.scenarios), ", ")

  measurements_count = 0
  for (x, y) in scenarios(p)
    measurements_count += length(measurements(y))
  end

  println(io, "Platform with $(length(models(p))) model(s), $(length(scenarios(p))) scenario(s), $measurements_count measurement(s)")
  println(io, "   Models: $models_names")
  println(io, "   Scenarios: $scn_names")
end

################################## Model ###########################################

abstract type AbstractModel end

"""
    struct Model{IF,OF,EV,SG,C,EA} <: AbstractModel
      init_func::IF
      ode_func::OF
      events::EV
      saving_generator::SG
      records_output::AbstractVector{Pair{Symbol,Bool}}
      constants::C
      events_active::EA
    end

Structure storing core properties of ODE model.
This represent the content of one namespace from a Heta platform.

To get list of model content use methods: constants(model), records(model), switchers(model).

To get the default model options use methods: `parameters(model)`,
`events_active(model)`, `events_save(model)`, `observables(model)`.
These options can be rewritten by a [`Scenario`]{@ref}.
"""
struct Model{IF,OF,EV,SG,C,EA} <: AbstractModel
  init_func::IF
  ode_func::OF
  events::EV # IDEA: use (:TimeEvent, ...) instead of TimeEvent(...)
  saving_generator::SG
  records_output::AbstractVector{Pair{Symbol,Bool}}
  constants::C # LArray{Float64,1,Array{Float64,1},(:a, :b)}
  events_active::EA
end

constants(m::Model) = [keys(m.constants)...]
records(m::Model) = first.(m.records_output)
switchers(m::Model) = [keys(m.events)...]

#parameters(m::Model) = collect(Pair{Symbol, Real}, pairs(m.constants))
events_active(m::Model) = collect(Pair{Symbol, Bool}, pairs(m.events_active))
events_save(m::Model) = [first(x) => (true,true) for x in pairs(m.events)]
observables(m::Model) = begin # observables
  only_true = filter((p) -> last(p), m.records_output)
  first.(only_true)
end

# auxilary function to display first n components of vector
function print_lim(x::Union{Vector, Tuple}, n::Int)
  first_n = ["$y" for y in first(x, n)]
  if length(x) > n
    push!(first_n, "...")
  end
  return join(first_n, ", ")
end
function print_lim(::Nothing, n::Int)
  return "-"
end

function Base.show(io::IO, mime::MIME"text/plain", m::AbstractModel)
  const_str = print_lim(constants(m), 10)
  record_str = print_lim(records(m), 10)
  switchers_str = print_lim(switchers(m), 10)

  println(io, "Model contains $(length(m.constants)) constant(s), $(length(m.records_output)) record(s), $(length(m.events)) switcher(s).")
  println(io, "   Constants: $const_str")
  println(io, "   Records: $record_str")
  println(io, "   Switchers (events): $switchers_str")
end

################################## Params ###########################################

struct Params{C,S}
  constants::C
  static::S
end

################################## Measurement ###########################################

# abstract type AbstractMeasurement end # not used
abstract type AbstractMeasurementPoint end

struct NormalMeasurementPoint{M,SD} <: AbstractMeasurementPoint
  t::Float64
  val::Float64
  scope::Symbol
  μ::M
  σ::SD
end

NormalMeasurementPoint(t,val,scope::Missing,μ,σ) = NormalMeasurementPoint(t,val,:ode_,μ,σ)

struct LogNormalMeasurementPoint{M,SD} <: AbstractMeasurementPoint
  t::Float64
  val::Float64
  scope::Symbol
  μ::M
  σ::SD
end

LogNormalMeasurementPoint(t,val,scope::Missing,μ,σ) = LogNormalMeasurementPoint(t,val,:ode_,μ,σ)

const MeasurementVector{P} = AbstractVector{P} where P<:AbstractMeasurementPoint

################################## Scenario ###########################################
abstract type AbstractScenario end

"""
  struct Scenario{F,P,M} <: AbstractScenario
    init_func::F
    prob::P
    measurements::M
  end

  Type representing simulation conditions, i.e. model variant with updated constants and outputs.

  To get the internal properties use methods: `saveat(scenario)`, `tspan(scenario)`, `parameters(scenario)`, `measurements(scenario)`
"""
struct Scenario{F,P,M} <: AbstractScenario
  init_func::F
  prob::P
  measurements::M
end

saveat(scn::Scenario) = scn.prob.kwargs[:callback].discrete_callbacks[1].affect!.saveat_cache
tspan(scn::Scenario) = scn.prob.tspan
parameters(scn::Scenario) = scn.prob.p.constants
measurements(scn::Scenario) = scn.measurements
# events_active(scn::Scenario)
# events_save(scn::Scenario)
# observables(scn::Scenario)

function Base.show(io::IO, mime::MIME"text/plain", scn::Scenario)
  saveat_str = print_lim(saveat(scn), 10)
  parameters_str = print_lim(
    keys(parameters(scn)),
    10
    )
  measurements_count = length(measurements(scn))

  if length(saveat(scn)) == 0
    time_points_str = "for tspan=$(tspan(scn))"
  else
    time_points_str = "for saveat=$(saveat_str)"
  end

  println(io, "Scenario $time_points_str")
  println(io, "   Time range (tspan): $(tspan(scn))")
  println(io, "   Time points to save the solution at (saveat): $(saveat_str)")
  println(io, "   Parameters: $(parameters_str)")
  println(io, "   Number of measurement points: $(measurements_count)")
end

################################## SimResult ###########################################
abstract type AbstractResult end

struct SavedValues{uType,tType,scopeType}
  u::uType
  t::tType
  scope::scopeType
end

function clear_savings(sv::SavedValues)
  !isempty(sv.u) && deleteat!(sv.u,1:length(sv.u))
  !isempty(sv.t) && deleteat!(sv.t,1:length(sv.t))
  !isnothing(sv.scope) && !isempty(sv.scope) && deleteat!(sv.scope,1:length(sv.scope))
  return nothing
end

struct Simulation{V,scopeType,P}
  vals::V
  scope::scopeType
  params::P
  status::Symbol
end

# copy fix is tmp needed not to rewrite SavedValues with new simulation
Simulation(sv::SavedValues, params, status::Symbol) = Simulation(
  DiffEqBase.SciMLBase.DiffEqArray(copy(sv.u),copy(sv.t)),
  sv.scope,
  params,
  status
)

status(s::Simulation) = s.status
times(s::Simulation) = s.vals.t
vals(s::Simulation) = s.vals.u
parameters(s::Simulation) = s.params

@inline Base.length(S::Simulation) = length(S.vals.t)

# tmp fix to support https://github.com/SciML/SciMLBase.jl/blob/ccaba96f4d7e29e9980cd4cd7270086fc5e542d6/src/ensemble/ensemble_analysis.jl#L56
function Base.getproperty(s::Simulation, sym::Symbol)
  if sym == :t
    return s.vals.t
  else
    return getfield(s, sym)
  end
end

"""
    struct SimResult{S, C<:Scenario} <: AbstractResult
      sim::S
      scenario::C
    end

Structure storing results from [`sim`]{@ref} method applied for one [`Scenario`]{@ref}.

To get the content use methods: `status(results)`, `times(results)`, `vals(results)`, `parameters(results)`, observables(results).

The results can be transformed using `DataFrame` method or visualized using `plot` method.
"""
struct SimResult{S, C<:Scenario} <: AbstractResult
  sim::S
  scenario::C
end

status(sr::SimResult) = status(sr.sim)
times(sr::SimResult) = times(sr.sim)
vals(sr::SimResult) = vals(sr.sim)
parameters(sr::SimResult) = parameters(sr.sim)
measurements(sr::SimResult) = sr.scenario.measurements

@inline Base.length(sr::SimResult) = length(sr.sim)

function Base.show(io::IO, mime::MIME"text/plain", sr::SimResult)
  dim2 = length(keys(sr.sim[1])) # number of observables
  dimentions_str = "$(length(sr))x$dim2"
  times_str = print_lim(times(sr), 10)
  outputs_str = print_lim(observables(sr), 10)
  parameters_str = print_lim(keys(parameters(sr)), 10)

  println(io, "$dimentions_str SimResult with status :$(status(sr)).")
  println(io, "    Solution status: $(status(sr))")
  println(io, "    Time points (times): $(times_str)")
  println(io, "    Observables (outputs): $(outputs_str)")
  println(io, "    Parameters: $(parameters_str)")
end

function Base.show(io::IO, mime::MIME"text/plain", srp::Pair{Symbol, S}, short::Bool = false) where S<:SimResult
  sr = last(srp)
  dim2 = length(keys(sr.sim[1])) # number of observables
  dimentions_str = "$(length(sr))x$dim2"

  short || println(io, "Pair{Symbol, SimResult}")
  println(io, "    :$(first(srp)) => $dimentions_str SimResult with status :$(status(sr)).")
end

function Base.show(io::IO, mime::MIME"text/plain", vector::Vector{Pair{Symbol, S}}) where S<:SimResult
  println(io, "$(length(vector))-element Vector{Pair{Symbol, SimResult}}")

  for x in vector
    show(io, mime, x, true)
  end
end

function Base.getindex(vector::Vector{Pair{Symbol, S}}, id::Symbol) where S<:SimResult
  ind = findfirst((x) -> first(x)===id, vector)
  if ind === nothing
    throw("Index :$id is not found.")
  else
    return vector[ind]
  end
end

################################## Monte-Carlo Simulation ##############################

"""
    struct MCResult{S,C} <: AbstractResult
      sim::S
      saveat::Bool
      scenario::C
    end

Structure storing results of [`mc`]{@ref} method applied for one `Scenario`.

To convert into tabular format, use `DataFrame` method.

The base visulaization can be done with `plot` method.
"""
struct MCResult{S,C} <: AbstractResult
  sim::S
  saveat::Bool
  scenario::C
  # converged
  # elapsed_time
end

@inline Base.length(mcr::MCResult) = length(mcr.sim)
parameters(mcr::MCResult) = [parameters(mcr[i]) for i in 1:length(mcr)]
vals(mcr::MCResult) = [vals(mcr[i]) for i in 1:length(mcr)]
status_summary(mcr::MCResult) = counter([s.status for s in mcr.sim])

function Base.show(io::IO, mime::MIME"text/plain", mcr::MCResult, short::Bool = false)
  # dimentions
  dim0 = length(mcr)
  if mcr.saveat
    dim1 = length(times(mcr.sim[1]))
  else
    dim1 = "?"
  end
  dim2 = length(mcr.sim[1]) # test only first Simulation
  dimentions_str = "$(dim0)x$(dim1)x$(dim2)"

  # other
  status_pairs = [":$(first(x)) x $(last(x))" for x in status_summary(mcr)]
  status = join(status_pairs, ", ")
  outputs_str = print_lim(observables(mcr.sim[1]), 10)
  parameters_1 = parameters(mcr.sim[1])
  parameters_str = isnothing(parameters_1) ? "-" : print_lim(keys(parameters_1), 10)

  println(io, "$(dimentions_str) MCResult with status $(status)" )
  if !short
  println(io, "    Solution status: $(status)")
  println(io, "    Observables (outputs): $(outputs_str)")
  println(io, "    Parameters: $(parameters_str)")
  end
end

function Base.show(io::IO, mime::MIME"text/plain", mcrp::Pair{Symbol, S}, short=false) where S<:MCResult
  short || println(io, "Pair{Symbol, MCResult}")
  print(io, "    :$(first(mcrp)) => ")
  Base.show(io, mime, last(mcrp), true)
end

function Base.show(io::IO, mime::MIME"text/plain", vector::Vector{Pair{Symbol, S}}) where S<:MCResult
  println(io, "$(length(vector))-element Vector{Pair{Symbol, MCResult}}")

  for x in vector
    show(io, mime, x, true)
  end
end

function Base.getindex(V::Vector{Pair{Symbol, S}}, id::Symbol) where S<:MCResult
  ind = findfirst((x) -> first(x)===id, V)
  if ind === nothing
    throw("Index :$id is not found.")
  else
    return V[ind]
  end
end

################################## Events ##############################################
abstract type AbstractEvent end

struct TimeEvent{F1,F2} <: AbstractEvent
  condition_func::F1
  affect_func::F2
  atStart::Bool
end

struct CEvent{F1,F2} <: AbstractEvent
  condition_func::F1
  affect_func::F2
  atStart::Bool
end

struct StopEvent{F1} <: AbstractEvent
  condition_func::F1
  atStart::Bool
end

################################## Fitting ###########################################
"""
    struct FitResult{L<:Real, I}
      obj::L
      optim::Vector{Pair{Symbol,Float64}}
      status::Symbol
      numevals::I
    end

Result of [`fit`](@ref).

Use `optim` method to get optimal values.

Use `object` method to get optimal objective function.

Use `status` to get status.

The optimal parameters can be saved in heta file, see [`save_as_heta`](@ref) method.
"""
struct FitResult{L<:Real, I}
  obj::L
  optim::Vector{Pair{Symbol,Float64}}
  status::Symbol
  numevals::I
end

function Base.show(io::IO, mime::MIME"text/plain", fr::FitResult)

  println(io, "FitResult with status :$(fr.status)")
  println(io, "   Status: $(fr.status)")
  println(io, "   Optimal values: $(fr.optim)")
  println(io, "   OF value: $(fr.obj)")
  println(io, "   OF count: $(fr.numevals)")
end

optim(f::FitResult) = f.optim
obj(f::FitResult) = f.obj
status(f::FitResult) = f.status


text:
module TestPkg

include("Module.jl")

greet() = print("Hello World")

export greet

end # module


text:
using Atom, Test, Atom.JSON, Logging, Atom.CSTParser, Pkg

joinpath′(files...) = Atom.fullpath(joinpath(files...))

atomjldir = joinpath′(@__DIR__, "..")
atomsrcdir = joinpath′(atomjldir, "src")
atomjlfile = joinpath′(atomsrcdir, "Atom.jl")
webiofile = joinpath′(atomsrcdir, "display", "webio.jl")
traceurfile = joinpath′(atomsrcdir, "profiler", "traceur.jl")

# files in `Atom` module except files in submodules
# TODO:
# currently both Revise-like and CSTPraser-based module traverse fails
# to detect lazily loaded files even if they are actually loaded.
atommodfiles = let
    files = []
    debuggerdir = joinpath′(atomsrcdir, "debugger")
    profilerdir = joinpath′(atomsrcdir, "profiler")

    for (d, ds, fs) in walkdir(atomsrcdir)
        # NOTE: update directories below when you create an new submodule
        # the 2 files below are in Atom module
        if d == debuggerdir
            push!(files, joinpath′(d, "debugger.jl"))
            continue
        end
        if d == profilerdir
            push!(files, joinpath′(d, "profiler.jl"))
            # push!(files, joinpath′(d, "traceur.jl"))
            continue
        end

        for f in fs
            f == "webio.jl" && continue
            # .jl check is needed for travis, who creates hoge.cov files
            endswith(f, ".jl") && push!(files, joinpath′(d, f))
        end
    end

    files
end

# mock modules
fixturedir = joinpath′(@__DIR__, "fixtures")
junkpath = joinpath′(fixturedir, "Junk.jl")
subjunkspath = joinpath′(fixturedir, "SubJunks.jl")
include(junkpath)

@testset "Atom.jl" begin
    include("utils.jl")
    include("misc.jl")
    include("display.jl")
    include("static/static.jl")
    include("environments.jl")
    include("modules.jl")

    include("eval.jl")
    include("outline.jl")
    include("completions/completions.jl")
    include("goto.jl")
    include("datatip.jl")
    include("workspace.jl")
    include("docs.jl")
    include("formatter.jl")
end


text:
# This file is a part of SIS3316Digitizers.jl, licensed under the MIT License (MIT).

function multiple_tries(thunk, n_tries::Integer, error_types::NTuple{N,Type}, action_description::String) where N
    for i in 1:n_tries
        # i > 1 && @info "$action_description, retry no $(i - 1)"
        try
            return thunk()
        catch err
            if !(typeof(err) in error_types)
                @error "$action_description failed with error $err, does not qualify for retry"
                rethrow()
            elseif i < n_tries
                @info "$action_description failed with error $err, retrying $((n_tries - i)) more times"
            else
                @error "$action_description failed with error $err, giving up after $n_tries tries"
                rethrow()
            end
        end
    end
    @assert false "Shouldn't have reached this part of code"
end
export multiple_tries

const MemAddr = UInt32
const RegVal = UInt32

struct RegisterReadOperation{A<:Unsigned,V<:Unsigned}
    addr::A
    val_ch::Channel{V}
end

struct RegisterWriteOperation{A<:Unsigned,V<:Unsigned}
    addr::A
    op::BitWriteOperation{V}
    ack_ch::Channel{Nothing}
end

const RegisterOperation{A<:Unsigned,V<:Unsigned} = Union{RegisterReadOperation{A,V}, RegisterWriteOperation{A,V}}

const OpsChannel = Channel{RegisterOperation{MemAddr, RegVal}}

mutable struct SIS3316Memory{GW<:AbstractVMEGateway}
    gw::GW
    ops_queue::OpsChannel
    fifo_lock::ReentrantLock
end

export SIS3316Memory

function SIS3316Memory(gw::AbstractVMEGateway)
    ops_queue = start_rwtask(gw)
    fifo_lock = ReentrantLock()
    SIS3316Memory(gw, ops_queue, fifo_lock)
end

Base.show(io::IO, mem::SIS3316Memory) = print(io, "SIS3316Memory($(repr(mem.gw)))")

Base.deepcopy(mem::SIS3316Memory) = SIS3316Memory(deepcopy(mem.gw))

@inline Base.Broadcast.broadcastable(mem::SIS3316Memory) = Ref(mem)

Base.isopen(mem::SIS3316Memory) = isopen(mem.ops_queue)

function Base.close(mem::SIS3316Memory)
    if isopen(mem.ops_queue)
        # @info "Closing $mem" # logging causes trouble in finalizers
        close(mem.ops_queue)
        # close(mem) # Could possibly cause race condition with rwtask shutdown
        # @info "Closed $mem" # logging causes trouble in finalizers
    end
end

function Base.getindex(mem::SIS3316Memory, addr::Integer)
    conv_addr = convert(MemAddr, addr)
    val_ch = Channel{RegVal}(1)
    push!(mem.ops_queue, RegisterReadOperation(conv_addr, val_ch))
    x = take!(val_ch)
    close(val_ch)
    x
end

function Base.getindex(mem::SIS3316Memory, ref::Union{Register,BitSelRef})
    x = mem[getaddress(ref)]
    getval(x, getlayout(ref))
end

function Base.setindex!(mem::SIS3316Memory, value::Integer, addr::Integer)
    conv_addr = convert(MemAddr, addr)
    conv_value = convert(RegVal, value)
    setindex!(mem, BitWriteOperation(conv_value, ~zero(conv_value), false), conv_addr)
end

function Base.setindex!(mem::SIS3316Memory, value::BitWriteOperation{RegVal}, addr::Integer)
    conv_addr = convert(MemAddr, addr)
    ack_ch = Channel{Nothing}(1)
    push!(mem.ops_queue, RegisterWriteOperation(conv_addr, value, ack_ch))
    take!(ack_ch)
    close(ack_ch)
    mem
end

function Base.setindex!(mem::SIS3316Memory, value, ref::Union{Register,BitSelRef})
    op = setval(BitWriteOperation{UInt32}(getaccessmode(ref)), getlayout(ref), value)
    setindex!(mem, op, getaddress(ref))
end

function start_rwtask(gw::AbstractVMEGateway)
    OpsChannel(100; taskref = nothing, spawn = true) do ops_queue
        write_ops = Vector{RegisterWriteOperation{MemAddr,RegVal}}()
        read_ops = Vector{RegisterReadOperation{MemAddr,RegVal}}()

        store_op(op::RegisterWriteOperation) = push!(write_ops, op)
        store_op(op::RegisterReadOperation) = push!(read_ops, op)

        try
            @debug "SIS3316MemState r/w task started"
            while true
                # @debug "Waiting for operations"
                wait(ops_queue)

                while isready(ops_queue)
                    op = store_op(take!(ops_queue))
                end

                # Give other tasks a chance to push more operations
                yield()

                # Try to maximize ops merging, don't process while more ops incoming
                if !isready(ops_queue)
                    execute_mem_ops!(gw, write_ops, read_ops)
                    @assert isempty(write_ops)
                    @assert isempty(read_ops)
                end
            end
        catch err
            if err isa InvalidStateException && !isopen(ops_queue)
                @debug "SIS3316MemState r/w task closed and terminating"
            else
                @error err
                rethrow()
            end
        finally
            close(gw)
        end
    end
end

function execute_mem_ops!(
    gw::AbstractVMEGateway,
    write_ops::Vector{RegisterWriteOperation{MemAddr,RegVal}},
    read_ops::Vector{RegisterReadOperation{MemAddr,RegVal}}
)
    try
        # @debug "Operations to execute" write_ops read_ops

        if !isempty(write_ops)
            let
                # @debug "Executing write operations"
                combined_writes = IdDict{MemAddr,BitWriteOperation{RegVal}}()
                for reg_write in write_ops
                    if haskey(combined_writes, reg_write.addr)
                        combined_writes[reg_write.addr] = merge(combined_writes[reg_write.addr], reg_write.op)
                    else
                        combined_writes[reg_write.addr] = reg_write.op
                    end
                end
                @debug "To write:" combined_writes

                write_addrs = sort(collect(keys(combined_writes)))
                wrtops = getindex.(Ref(combined_writes), write_addrs)
                write_values = map(x -> x.value, wrtops)
                with_mask_idxs = findall(is_masked, wrtops)
                with_mask_addrs = write_addrs[with_mask_idxs]
                with_mask_ops = wrtops[with_mask_idxs]
                currvals = Vector{RegVal}(undef, length(with_mask_addrs))
                if !isempty(with_mask_addrs)
                    read_registers!(gw, with_mask_addrs, currvals)
                end
                # @debug "Read current values $(repr(currvals)) from addresses $(repr(with_mask_addrs)) before writing"
                write_values[with_mask_idxs] .= bset.(currvals, with_mask_ops)
                write_registers!(gw, write_addrs, write_values)

                while !isempty(write_ops)
                    op = first(write_ops)
                    result = nothing
                    push!(op.ack_ch, result)
                    popfirst!(write_ops)
                end
            end
        end

        if !isempty(read_ops)
            let
                # @debug "Executing read operations"
                combined_reads = sort(collect(Set(reg_read.addr for reg_read in read_ops)))
                @debug "To read:" combined_reads

                read_addrs = sort(collect(Set(reg_read.addr for reg_read in read_ops)))
                read_values = Vector{RegVal}(undef, length(read_addrs))
                read_registers!(gw, read_addrs, read_values)
                read_results = IdDict{MemAddr,RegVal}()
                for (a, v) in zip(read_addrs, read_values)
                    read_results[a] = v
                end

                # @debug "Read results" read_addrs read_values read_results
                while !isempty(read_ops)
                    op = first(read_ops)
                    result = read_results[op.addr]
                    @debug "Pushing to result channel: $(repr(result))"
                    push!(op.val_ch, result)
                    popfirst!(read_ops)
                end
            end
        end
    catch err
        if err isa EOFError
            @warn "UDP Timeout"
        else
            @warn err
        end

        # Cancel all pending operations:

        for op in write_ops
            close(op.ack_ch)
        end
        empty!(write_ops)

        for op in read_ops
            close(op.val_ch)
        end
        empty!(read_ops)
    finally
        @assert isempty(write_ops)
        @assert isempty(read_ops)
    end
end


text:
@testset "SplitPath" begin
    p, q = synthetic_gradient(SynthSphere(), radius = 50)
    Z = convert_gradient(SplitPath(), p, q)
    @test Z[38,38] ≈ 0.0
    @test Z[1,1] ≈ 0.0
    @test Z[76,76] ≈ 46.69574348683868
end


text:
#
# Definitions to move to NDTensors
#

storage(T::Tensor) = store(T)


text:
"""
    stfr(M0, T1, T2, Δω, κ, Tfree, Tg, α, β, ϕ, TE, [spoil])

Compute the steady-state signal generated from an STFR scan of the given tissue
parameters with the given scan parameters, assuming instantaneous RF pulses and
ideal spoiling.

# Arguments
- `M0::Real`: Equilibrium magnetization
- `T1::Real`: Spin-lattice recovery time constant (ms)
- `T2::Real`: Spin-spin recovery time constant (ms)
- `Δω::Real`: Off-resonance frequency (rad/s)
- `κ::Real`: Variation in tip angles
- `Tfree::Real`: Free precession time (ms)
- `Tg::Real`: Spoiler gradient time (ms)
- `α::Real`: Tip-down angle (rad)
- `β::Real`: Tip-up angle (rad)
- `ϕ::Real`: Phase of tip-up RF pulse (rad)
- `TE::Real`: Echo time (ms)
- `spoil::Val = Val(true)`: Whether to include RF spoiling

# Return
- `signal::Complex`: Signal generated from the STFR scan
"""
function stfr(
    M0::Real,
    T1::Real,
    T2::Real,
    Δω::Real,
    κ::Real,
    Tfree::Real,
    Tg::Real,
    α::Real,
    β::Real,
    ϕ::Real,
    TE::Real,
    ::Val{spoil} = Val(true)
) where {spoil}

    if spoil
        return   stfr_spoil(M0, T1, T2, Δω, κ, Tfree, Tg, α, β, ϕ, TE)
    else
        return stfr_unspoil(M0, T1, T2, Δω, κ, Tfree, Tg, α, β, ϕ, TE)
    end

end

"""
    stfr(M0, ff, T1f, T1s, T2f, T2s, Δωf, Δω, κ, Tfree, Tg, α, β, ϕ, TE, [spoil])

Compute the steady-state signal generated from an STFR scan of the given
two-compartment tissue parameters with the given scan parameters, assuming
instantaneous RF pulses, ideal spoiling, and no exchange between the two
compartments.

# Arguments
- `M0::Real`: Equilibrium magnetization
- `ff::Real`: Fraction of tissue that is fast compartment
- `T1f::Real`: Fast compartment spin-lattice recovery time constant (ms)
- `T1s::Real`: Slow compartment spin-lattice recovery time constant (ms)
- `T2f::Real`: Fast compartment spin-spin recovery time constant (ms)
- `T2s::Real`: Slow compartment spin-spin recovery time constant (ms)
- `Δωf::Real`: Fast compartment off-resonance frequency (added to Δω) (rad/s)
- `Δω::Real`: Bulk off-resonance frequency (rad/s)
- `κ::Real`: Variation in tip angles
- `Tfree::Real`: Free precession time (ms)
- `Tg::Real`: Spoiler gradient time (ms)
- `α::Real`: Tip-down angle (rad)
- `β::Real`: Tip-up angle (rad)
- `ϕ::Real`: Phase of tip-up RF pulse (rad)
- `TE::Real`: Echo time (ms)
- `spoil::Val = Val(true)`: Whether to include RF spoiling

# Return
- `signal::Complex`: Signal generated from the two-compartment STFR scan
"""
function stfr(
    M0::Real,
    ff::Real,
    T1f::Real,
    T1s::Real,
    T2f::Real,
    T2s::Real,
    Δωf::Real,
    Δω::Real,
    κ::Real,
    Tfree::Real,
    Tg::Real,
    α::Real,
    β::Real,
    ϕ::Real,
    TE::Real,
    ::Val{spoil} = Val(true)
) where {spoil}

    if spoil
        return   ff  * stfr_spoil(M0, T1f, T2f, Δω + Δωf,
                                  κ, Tfree, Tg, α, β, ϕ, TE) +
            (1 - ff) * stfr_spoil(M0, T1s, T2s, Δω,
                                  κ, Tfree, Tg, α, β, ϕ, TE)
    else
        return ff  * stfr_unspoil(M0, T1f, T2f, Δω + Δωf,
                                  κ, Tfree, Tg, α, β, ϕ, TE) +
          (1 - ff) * stfr_unspoil(M0, T1s, T2s, Δω,
                                  κ, Tfree, Tg, α, β, ϕ, TE)
    end

end

"Compute the STFR signal using the spoiled STFR signal model."
function stfr_spoil(M0, T1, T2, Δω, κ, Tfree, Tg, α, β, ϕ, TE)

    κα = κ * α
    κβ = κ * β
    θfree = Δω * Tfree/1000

    n = exp(-Tg/T1) * (1 - exp(-Tfree/T1)) * sin(κα) * cos(κβ) +
        (1 - exp(-Tg/T1)) * sin(κα)

    d = 1 - exp(-Tg/T1 - Tfree/T2) * sin(κα) * sin(κβ) * cos(θfree - ϕ) -
            exp(-Tg/T1 - Tfree/T1) * cos(κα) * cos(κβ)

    return (M0 * n / d) * exp(-TE/T2 - im*Δω*TE/1000)

end

"Compute the STFR signal using the unspoiled STFR signal model."
function stfr_unspoil(M0, T1, T2, Δω, κ, Tfree, Tg, α, β, ϕ, TE)

    κα = κ * α
    κβ = κ * β
    θfree = Δω * Tfree/1000
    Ef1 = exp(-Tfree/T1)
    Ef2 = exp(-Tfree/T2)
    Eg1 = exp(-Tg/T1)
    Eg2 = exp(-Tg/T2)

    a = -im * Eg2 * (Ef2 * (-1 + Eg1 + (-1 + Ef1) * Eg1 * cos(κβ)) *
        cos(θfree - ϕ) * sin(κα) + (Ef1 * (-1 + Eg1) + (-1 + Ef1) * cos(κα)) *
        sin(κβ) + im * Ef2 * (-1 + Eg1 + (-1 + Ef1) * Eg1 * cos(κβ)) * sin(κα) *
        sin(θfree - ϕ))

    b = Eg2 * (Ef2 * ((-1 + Ef1) * Eg1 + (-1 + Eg1) * cos(κβ)) * cos(θfree - ϕ) *
        sin(κα) - (-1 + Ef1 + Ef1 * (-1 + Eg1) * cos(κα)) * sin(κβ) + im * Ef2 *
        ((-1 + Ef1) * Eg1 + (-1 + Eg1) * cos(κβ)) * sin(κα) * sin(θfree - ϕ))

    c = im * ((-1 + Eg1 + (-1 + Ef1) * Eg1 * cos(κβ)) * sin(κα) + Ef2 * Eg2^2 *
        (Ef1 * (-1 + Eg1) + (-1 + Ef1) * cos(κα)) * sin(κβ) * (cos(θfree - ϕ) +
        im * sin(θfree - ϕ)))

    d = Eg2 * (-Ef2 * (-1 + Ef1 * Eg1) * (1 + cos(κα) * cos(κβ)) * cos(θfree - ϕ)
        + (Ef1 - Ef2^2 * Eg1) * sin(κα) * sin(κβ))

    e = Ef2 * (-1 + Ef1 * Eg1) * Eg2 * (cos(κα) + cos(κβ)) * sin(θfree - ϕ)

    f = -1 + Ef1 * Ef2^2 * Eg1 * Eg2^2 + (Ef1 * Eg1 - Ef2^2 * Eg2^2) * cos(κα) *
        cos(κβ) + Ef2 * (Eg1 - Ef1 * Eg2^2) * cos(θfree - ϕ) * sin(κα) * sin(κβ)

    s = sqrt(f^2 - d^2 - e^2)

    tmp1 = -c / s
    tmp2 = (a * d + b * e) / (d^2 + e^2)
    tmp3 = (-f - s) / s

    return M0 * (tmp1 - tmp2 * tmp3) * exp(-TE/T2 - im*Δω*TE/1000)

end


text:
"""
    NonlinearSolvers

A set of solvers for systems of non-linear equations
"""
module NonlinearSolvers

using DocStringExtensions
using ForwardDiff

export solve!

const FTypes = Union{AbstractFloat, AbstractArray}

"""
    AbstractNonlinearSolverMethod

A super type for non-linear systems
"""
abstract type AbstractNonlinearSolverMethod{FTypes} end

"""
    method_args(method::AbstractNonlinearSolverMethod)

Return tuple of positional args for `AbstractNonlinearSolverMethod`.
"""
function method_args end

include("tolerance_types.jl")
include("solution_types.jl")

include("newton_method.jl")
include("newton_method_ad.jl")

# Main entry point: Dispatch to specific method
"""
    solve!(
        method::AbstractNonlinearSolverMethod{FT},
        soltype::SolutionType = CompactSolution(),
        tol::Union{Nothing, AbstractTolerance} = nothing,
        maxiters::Union{Nothing, Int} = 10_000,
    )

Solve the non-linear system given
 - `method` the numerical method
 - `soltype` the solution type (`CompactSolution` or `VerboseSolution`)
 - `tol` the stopping tolerance
 - `maxiters` the maximum number of iterations to perform
"""
function solve!(
    method::AbstractNonlinearSolverMethod{FT},
    soltype::SolutionType = CompactSolution(),
    tol::Union{Nothing, AbstractTolerance} = nothing,
    maxiters::Union{Nothing, Int} = 10_000,
) where {FT <: FTypes, F <: Function}
    if tol === nothing
        tol = ResidualTolerance{eltype(FT)}(sqrt(eps(FT)))
    end
    return solve!(method, method_args(method)..., soltype, tol, maxiters)
end

end


text:
######################################################################
# smc_multi.jl
#
# Like smc.jl, but for multi-agent trajectories
#
# In theory, this shouldn't be so different from single version
######################################################################

using Distances
export smc_trajectory_s
export smc_trajectory_st

function smc_trajectory(em::ErgodicManager, vtm::VTM; verbose::Bool=true, umax::Float64=1.0)

	num_agents = length(vtm)
	N = vtm[1].N
	dt = vtm[1].h
	#ud = VVF(tm.N)
	#xd = VVF(tm.N+1)
	xds = VVVF(num_agents)
	uds = VVVF(num_agents)
	for j = 1:num_agents
		xds[j] = VVF(N+1)
		uds[j] = VVF(N)
		xds[j][1] = deepcopy(vtm[j].x0)
	end

	ck = zeros(em.K+1, em.K+1)
	for n = 1:N

		# update ck
		update_ck(em, xds, n, ck)

		# compute B
		for j = 1:num_agents
			Bx, By = compute_Bj(em, xds[j], n, dt, ck, num_agents)

			# normalize
			den = sqrt(Bx*Bx + By*By)
			ux = -umax * Bx / den
			uy = -umax * By / den

			uds[j][n] = [ux, uy]
			xds[j][n+1] = integrate(vtm[j], xds[j][n], uds[j][n])
		end
	end

	return vvvf2vvf(xds, uds)
end

function update_ck(em::ErgodicManager, xds::VVVF, n::Int, ck::MF)
	num_agents = length(xds)
	Lx = em.domain.lengths[1]
	Ly = em.domain.lengths[2]
	for k1 = 0:em.K
		for k2 = 0:em.K
			hk = em.hk[k1+1, k2+1]
			ck[k1+1, k2+1] *= (n-1)
			for j = 1:num_agents
				dxn = xds[j][n][1] - x_min(em)
				dyn = xds[j][n][2] - y_min(em)

				cx = cos(k1*pi*dxn/Lx)
				cy = cos(k2*pi*dyn/Ly)

				# compute ck
				ck[k1+1, k2+1] += cx*cy/(hk*num_agents)
			end
			ck[k1+1, k2+1] /= n
		end
	end
end

# TODO: this can be done a lot quicker
function compute_Bj(em::ErgodicManager, xd::VVF,n::Int,h::Float64,ck::MF, num_agents::Int)
	Lx = em.domain.lengths[1]
	Ly = em.domain.lengths[2]
	dxn = xd[n][1] - x_min(em)
	dyn = xd[n][2] - y_min(em)
	Bx = By = 0.0
	for k1 = 0:em.K
		for k2 = 0:em.K
			hk = em.hk[k1+1, k2+1]
			cx = cos(k1*pi*dxn/Lx)
			cy = cos(k2*pi*dyn/Ly)

			# compute ck
			#fk = ( (n-1)*ck[k1+1,k2+1] + cx*cy/hk ) / n
			#ck[k1+1, k2+1] = fk

			# t * Lambda_k * (c_k - phi_k)
			de = ck[k1+1, k2+1] - em.phik[k1+1,k2+1]
			LS = em.Lambda[k1+1,k2+1] * h*n * num_agents * de

			# multiplying by gradient
			Bx += -LS * k1*pi * sin(k1*pi*dxn/Lx) * cy / (hk*Lx)
			By += -LS * k2*pi * cx * sin(k2*pi*dyn/Ly) / (hk*Ly)
		end
	end
	return Bx, By
end

######################################################################
# This is a modified version of smc_multi.  It allows for
# entirely separate local agent trajectories - smc_multi has been kept
# for comparison purposes.
######################################################################

function smc_trajectory_s(em::ErgodicManager, vtm::VTM, threshold::Float64; verbose::Bool=true, umax::Float64=1.0)

	num_agents = length(vtm)
	N = vtm[1].N
	dt = vtm[1].h
	#ud = VVF(tm.N)
	#xd = VVF(tm.N+1)
	xds = VVVF(num_agents)
	uds = VVVF(num_agents)
	for j = 1:num_agents
		xds[j] = VVF(N+1)
		uds[j] = VVF(N)
		xds[j][1] = deepcopy(vtm[j].x0)
	end

	#ck = zeros(em.K+1, em.K+1)
	for n = 1:N

		# update ck
		update_ck_s(em, xds, n, vtm, threshold)

		# compute B
		for j = 1:num_agents
			Bx, By = compute_Bj_s(em, xds[j], n, dt, vtm[j].ck, num_agents)

			# normalize
			den = sqrt(Bx*Bx + By*By)
			ux = -umax * Bx / den
			uy = -umax * By / den

			uds[j][n] = [ux, uy]
			xds[j][n+1] = integrate(vtm[j], xds[j][n], uds[j][n])
			#println(n)
			vtm[j].xy = xds[j][n]
		end
	end

	return vvvf2vvf(xds, uds)
end

function update_ck_s(em::ErgodicManager, xds::VVVF, n::Int, vtm::VTM, threshold::Float64)
	num_agents = length(xds)
	Lx = em.domain.lengths[1]
	Ly = em.domain.lengths[2]
	local_agents = 0

	for agent = 1: num_agents
		for k1 = 0:em.K
			for k2 = 0:em.K
				hk = em.hk[k1+1, k2+1]
				vtm[agent].ck[k1+1, k2+1] *= (n-1)
				for j = 1:num_agents
					#Check for local agents
					local_agents = 0
					for i = 1:num_agents
						eu = euclidean(vtm[agent].xy, vtm[i].xy)
						if eu < threshold
							local_agents += 1
						end
					end
					#if statement goes here to check for agent distance
					euc = euclidean(vtm[agent].xy, vtm[j].xy)
					#print(vtm[agent].x0, vtm[j].x0)
					#println(euc)
					if euc < threshold
						dxn = xds[j][n][1] - x_min(em)
						dyn = xds[j][n][2] - y_min(em)

						cx = cos(k1*pi*dxn/Lx)
						cy = cos(k2*pi*dyn/Ly)

						# compute ck
						vtm[agent].ck[k1+1, k2+1] += cx*cy/(hk*local_agents)
					end
				end
				vtm[agent].ck[k1+1, k2+1] /= n
			end
		end
	end
end

# TODO: this can be done a lot quicker
function compute_Bj_s(em::ErgodicManager, xd::VVF,n::Int,h::Float64,ck::MF, num_agents::Int)
	Lx = em.domain.lengths[1]
	Ly = em.domain.lengths[2]
	dxn = xd[n][1] - x_min(em)
	dyn = xd[n][2] - y_min(em)
	Bx = By = 0.0
	local_agents = 0

	for k1 = 0:em.K
		for k2 = 0:em.K
			hk = em.hk[k1+1, k2+1]
			cx = cos(k1*pi*dxn/Lx)
			cy = cos(k2*pi*dyn/Ly)

			# compute ck
			#fk = ( (n-1)*ck[k1+1,k2+1] + cx*cy/hk ) / n
			#ck[k1+1, k2+1] = fk

			# t * Lambda_k * (c_k - phi_k)
			de = ck[k1+1, k2+1] - em.phik[k1+1,k2+1]
			LS = em.Lambda[k1+1,k2+1] * h*n * num_agents * de

			# multiplying by gradient
			Bx += -LS * k1*pi * sin(k1*pi*dxn/Lx) * cy / (hk*Lx)
			By += -LS * k2*pi * cx * sin(k2*pi*dyn/Ly) / (hk*Ly)
		end
	end
	return Bx, By
end

# For plotting time-elapsed GIF's (R2T Domain)
function smc_trajectory_st(em::ErgodicManagerR2T, vtm::VTM, threshold::Float64; verbose::Bool=true, umax::Float64=1.0)

	num_agents = length(vtm)
	N = vtm[1].N
	dt = vtm[1].h
	#ud = VVF(tm.N)
	#xd = VVF(tm.N+1)
	xds = VVVF(num_agents)
	uds = VVVF(num_agents)
	for j = 1:num_agents
		xds[j] = VVF(N+1)
		uds[j] = VVF(N)
		xds[j][1] = deepcopy(vtm[j].x0)
	end

	#ck = zeros(em.K+1, em.K+1)
	for n = 1:N

		# update ck
		update_ck_st(em, xds, n, vtm, threshold)

		# compute B
		for j = 1:num_agents
			Bx, By = compute_Bj_st(em, xds[j], n, dt, vtm, vtm[j].ckt, num_agents)

			# normalize
			den = sqrt(Bx*Bx + By*By)
			ux = -umax * Bx / den
			uy = -umax * By / den

			uds[j][n] = [ux, uy]
			xds[j][n+1] = integrate(vtm[j], xds[j][n], uds[j][n])
			#println(n)
			vtm[j].xy = xds[j][n]
		end
	end

	return vvvf2vvf(xds, uds)
end

function update_ck_st(em::ErgodicManagerR2T, xds::VVVF, n::Int, vtm::VTM, threshold::Float64)
	num_agents = length(xds)
	Lx = em.domain.lengths[1]
	Ly = em.domain.lengths[2]
	t = n*vtm[1].h
	T = vtm[1].N*vtm[1].h
	local_agents = 0
	# for i = 1: num_agents
	# 	euc = euclidean(vtm[agent].xy, vtm[j].xy)
	# 	if euc < threshold
	# 		local_agents += 1
	# 	end
	# end

	for agent = 1: num_agents
		for k1 = 0:em.K
			for k2 = 0:em.K
				for k3 = 0: em.K
					hk = em.hk[k1+1, k2+1, k3+1]
					vtm[agent].ckt[k1+1, k2+1, k3+1] *= (n-1)
					for j = 1:num_agents
						#This is for averaging local neighbours
						local_agents = 0
						for i = 1:num_agents
							eu = euclidean(vtm[agent].xy, vtm[i].xy)
							if eu < threshold
								local_agents += 1
							end
						end
						#if statement goes here to check for agent distance
						euc = euclidean(vtm[agent].xy, vtm[j].xy)
						if euc < threshold
							dxn = xds[j][n][1] - x_min(em)
							dyn = xds[j][n][2] - y_min(em)

							cx = cos(k1*pi*dxn/Lx)
							cy = cos(k2*pi*dyn/Ly)
							ct = cos(k3*pi*t/T)

							# compute ck
							vtm[agent].ckt[k1+1, k2+1, k3+1] += cx*cy*ct/(hk*local_agents)
						end
					end
					vtm[agent].ckt[k1+1, k2+1, k3+1] /= n
				end
			end
		end
	end
end

# TODO: this can be done a lot quicker
function compute_Bj_st(em::ErgodicManagerR2T, xd::VVF,n::Int,h::Float64, vtm::VTM,ckt::Array{Float64,3}, num_agents::Int)
	Lx = em.domain.lengths[1]
	Ly = em.domain.lengths[2]
	dxn = xd[n][1] - x_min(em)
	dyn = xd[n][2] - y_min(em)
	Bx = By = 0.0
	t = n*vtm[1].h
	T = vtm[1].N*vtm[1].h
	local_agents = 0

	for k1 = 0:em.K
		for k2 = 0:em.K
			for k3 = 0:em.K

				hk = em.hk[k1+1, k2+1, k3+1]
				cx = cos(k1*pi*dxn/Lx)
				cy = cos(k2*pi*dyn/Ly)
				ct = cos(k3*pi*t/T)
				# compute ck
				#fk = ( (n-1)*ck[k1+1,k2+1] + cx*cy/hk ) / n
				#ck[k1+1, k2+1] = fk

				# t * Lambda_k * (c_k - phi_k)
				de = ckt[k1+1, k2+1, k3+1] - em.phik[k1+1,k2+1, k3+1]
				LS = em.Lambda[k1+1,k2+1, k3+1] * h*n * num_agents * de

				# multiplying by gradient
				Bx += -LS * k1*pi * sin(k1*pi*dxn/Lx) * cy  * ct/ (hk*Lx)
				By += -LS * k2*pi * cx * sin(k2*pi*dyn/Ly) * ct / (hk*Ly)

			end
		end
	end
	return Bx, By
end


text:
using OAuth, HTTP, JSON

# Set basic variables for authentication
oauth_consumer_key = "3APw4LnxlOgKzPoaa7EFwbOXZ"
oauth_consumer_secret = "oP9GrVUow2ivMsijrO0yxIhVhA2T4iltjeKPBaFbCbskeL1VN5"

oauth_token = "1207655773783769088-Ar8GUF6tkNe7tvZtvxzCgZ4IYivZWV"
oauth_token_secret = "TBqmc5bBzbrAeaNEpu1AfLqcbpF5IhKht9DF91vBgESKg"

# All the API URLs to call for various GET and POST requests
retweet_url = "https://api.twitter.com/1.1/statuses/retweet/:id.json"
get_following = "https://api.twitter.com/1.1/statuses/home_timeline.json"

# Number of tweets to get every 5 minutes
# I'm using just below Twitter's limit
# News is fast, and the bot should keep pace.
number_of_tweets = 14

# Words that determine if a tweet is considered news
trigger_words = ["now", "new", "latest", "see", "update", "feature", "check", "announce", "announcing", "announced", "today", "unveil", "introducing", "release"]

# GET and POST functions
get_oauth(endpoint, options) = oauth_request_resource(endpoint, "GET", options, oauth_consumer_key, oauth_consumer_secret, oauth_token, oauth_token_secret)
tweet_oauth(endpoint, options) = oauth_request_resource(endpoint, "POST", options, oauth_consumer_key, oauth_consumer_secret, oauth_token, oauth_token_secret)

# Our main function
function tweet()

    # GET the 10 most recent tweest in the last 15 minutes
    tweets = get_oauth(get_following, Dict("count" => "$number_of_tweets"))

    # Convert HTTP request to JSON and process it
    jsonCode = String(tweets.body)
    data = JSON.parse(jsonCode)

    # List of tweets
    acquired_tweets = []

    # Check if each tweet is considered news, and append its id to the list of tweets if it is news
    for i in 1:length(data)
        id = data[i]["id_str"]
        for trigger_word in trigger_words
            if occursin(trigger_word, lowercase(data[i]["text"]))
                push!(acquired_tweets, id)
            end
        end
    end

    # For every new tweet, retweet the tweet using its id.
    for j in 1:length(acquired_tweets)
    # If a 403 error is returned, ignore it and move on to the next tweet
    # This happens when we try to tweet a tweet that has already been retweeted by us.
        try
            id = acquired_tweets[j]
            options = Dict("id" => "$id")
            retweet_url = "https://api.twitter.com/1.1/statuses/retweet/$id.json"
            mytweet = tweet_oauth(retweet_url, options)
        catch
            continue
        end
    end
end

tweet()


text:
using Test
using CastepReader
using Unitful
using UnitfulAtomic
using FFTW

const cr=CastepReader
const testdata = joinpath(pkgdir(cr), "test/testdata")

gaas = joinpath(testdata, "GaAs")

# Test file name functions
@test cr.suffix("a.b") == "b"
@test cr.stem("x/a.b") == "x/a"
@test cr.filename("c/x/a.b") == "a.b"
@test cr.stem(cr.filename("c/x/a.b")) == "a"

@testset "Bands" begin
    kpts, kw, kidx, eigen, ne, efermi, cell_file = cr.read_bands_castep_raw(joinpath(gaas, "GaAs.bands"))
    @test size(kpts) == (3, 10)
    @test size(kw)[1] == 10
    @test size(eigen) == (1, 10, 40)
    @test kpts[:, 8] == [0., 0., 0.]
    @test kw[8] ≈ 0.008
    @test kidx[8] == 8
    @test ne == 28
    @test efermi[1] == uconvert(u"eV", 0.113556u"Eh_au")
    @test cell_file[1] == uconvert(u"Å", 5.432963u"a0_au")

    # Read the optical matrix
    ns, nk, nb = size(eigen)
    bin = cr.read_ome_bin(joinpath(gaas, "GaAs.ome_bin"), nb, nk, ns)
    @test size(bin) == (nb, nb, 3, nk, ns)
    @test bin[1] ≈ (0.00022207179206183356 - 0.0im) * cr.HaBohr

    # TODO - add test for cst_ome
    cst = cr.read_cst_ome(joinpath(gaas, "GaAs.cst_ome"), nb, nk, ns)

    @test size(cst) == (nb, nb, 3, nk, ns)
    @test cst[1] ≈ (0.0004183035272252899 - 0.0im) * cr.HaBohr2perAng

end

# Test reading the cell file
@testset "read cell" begin
    cell_file = cr.read_cell_castep(joinpath(gaas, "GaAs.cell"))
    @test :symmetry_ops in keys(cell_file.blocks)
    @test begin
        :kpoints_mp_spacing in keys(cell_file.kws) && cell_file.kws[:kpoints_mp_spacing] == "0.07"
    end
    cellmat = cr.read_lattice(cell_file)
    @test cellmat[1] == 2.875u"Å"

    symm, disp = cr.read_symm_disp(cell_file)
    @test begin
        symm[1, 1, 2] == 0. && symm[2, 1, 2] == -1. && all(disp .== 0.)
    end
end

@testset "castep_bin" begin
    bin = joinpath(gaas, "GaAs.castep_bin")
    output = read_castep_check(bin)
    @test :FORCES in keys(output)
    @test output[:VERSION_NUMBER] ≈ 20.1100006103
    @test output[:NUM_IONS]  == 2
    @test output[:NUM_SPECIES]  == 2
    @test output[:MAX_IONS_IN_SPECIES]  == 1

    @test size(output[:FORCES]) == (3,1, 2)
end

@testset "wave" begin
    checkfile = joinpath(gaas, "GaAs.check")
    output = read_castep_check(checkfile)
    wavef = CastepReader.WaveFunction(output)
    ngx = 25
    @test wavef.ngx == ngx
    out = CastepReader.sliceband(wavef, 1)
    @test size(out, 1) == ngx
    ifft!(wavef)
    @test wavef.real_space == true
    fft!(wavef)
    @test wavef.real_space == false
    # Charge density - with out interpolation
    ifft!(wavef)
    out = CastepReader.chargedensity(wavef, ones(wavef.nbands, wavef.nkpts, wavef.nspins), ones(wavef.nkpts))
    @assert out.ngx == ngx

    # For test updating wave coefficients
    @testset "update" begin
        dst = joinpath(mktempdir(), "GaAs.check")
        cp(checkfile, dst)
        data = read_castep_check(dst)
        getcoeff = data -> data[:WAVEFUNCTION_DATA].coeffs
        getcoeff(data)[1011] = 0.
        # Update the data
        CastepReader.update_wavefunction(dst, getcoeff(data))
        new_data = read_castep_check(dst)
        @test getcoeff(new_data)[1011] == 0.

        # Cannot update castep_bin files - should throw assertion error
        binfile = joinpath(gaas, "GaAs.castep_bin")
        dst = joinpath(mktempdir(), "GaAs.castep_bin")
        cp(binfile, dst)
        @test_throws AssertionError  CastepReader.update_wavefunction(dst, getcoeff(data))

        # Update with given wavefunction type
        wf = CastepReader.WaveFunction(data)
        dst = joinpath(mktempdir(), "GaAs.check")
        CastepReader.update_wavefunction(checkfile, dst, wf, data)
        new_data = read_castep_check(dst)
        @test getcoeff(new_data)[1011] == 0.
    end
end

@testset "pdos_bin" begin
    bin = joinpath(gaas, "GaAs.pdos_bin")
    output = read_pdos_bin(bin)
    @test size(output.kpoints)[2] == 10
    @test size(output.pdos_weights) == (20, 40, 10, 1)
    @test maximum(output.species) == 2

    m_mask = cr.m_channel_mask(output)
    @test m_mask[1] == Orbitals.dz2
    @test m_mask[6] == Orbitals.s
    @test m_mask[10] == Orbitals.s

    indices = cr.pdos_index_by_site(output)
    @test indices[1][Orbitals.dz2] == [1]
    @test indices[1][Orbitals.s] == [6, 10]
end

@testset "om" begin
    seed = joinpath(gaas, "GaAs")
    @testset "ome_bin" begin
        output = cr.read_ome_bin(seed)
        @test size(output)[4] == 10
    end

    @testset "cst_ome" begin
        output = cr.read_cst_ome(seed)
        @test size(output)[4] == 10
    end

    output = cr.read_om_castep(seed)
    @test size(output)[4] == 10

end

@testset "writers" begin
    @testset "chgcar" begin
        buffer = IOBuffer()
        CastepReader.write_chgcar(buffer, [rand(10, 10, 10)])
        seek(buffer, 0)
        content = take!(buffer)
        @test endswith(String(content), "\n")
    end
end