text:
"""
    GD53.jl

This is a Julia port of a MATLAB implementation of batch and incremental
Generalized Dunn's Index 53 (GD53) Cluster Validity Index.

Authors:
MATLAB implementation: Leonardo Enzo Brito da Silva
Julia port: Sasha Petrenko <sap625@mst.edu>

REFERENCES
[1] A. Ibrahim, J. M. Keller, and J. C. Bezdek, "Evaluating Evolving Structure
in Streaming Data With Modified Dunn's Indices," IEEE Transactions on Emerging
Topics in Computational Intelligence, pp. 1–12, 2019.
[2] M. Moshtaghi, J. C. Bezdek, S. M. Erfani, C. Leckie, and J. Bailey,
"Online Cluster Validity Indices for Streaming Data," ArXiv e-prints, 2018,
arXiv:1801.02937v1 [stat.ML].
[3] M. Moshtaghi, J. C. Bezdek, S. M. Erfani, C. Leckie, J. Bailey, "Online
cluster validity indices for performance monitoring of streaming data clustering,"
Int. J. Intell. Syst., pp. 1–23, 2018.
[4] J. C. Dunn, "A fuzzy relative of the ISODATA process and its use in detecting
compact well-separated clusters," J. Cybern., vol. 3, no. 3 , pp. 32–57, 1973.
[5] J. C. Bezdek and N. R. Pal, "Some new indexes of cluster validity," IEEE
Trans. Syst., Man, and Cybern., vol. 28, no. 3, pp. 301–315, Jun. 1998.
"""

"""
    GD53

The stateful information of the Generalized Dunn's Index 53 (GD53) Cluster Validity Index.

# References
1. A. Ibrahim, J. M. Keller, and J. C. Bezdek, "Evaluating Evolving Structure in Streaming Data With Modified Dunn's Indices," IEEE Transactions on Emerging Topics in Computational Intelligence, pp. 1–12, 2019.
2. M. Moshtaghi, J. C. Bezdek, S. M. Erfani, C. Leckie, and J. Bailey, "Online Cluster Validity Indices for Streaming Data," ArXiv e-prints, 2018, arXiv:1801.02937v1 [stat.ML].
3. M. Moshtaghi, J. C. Bezdek, S. M. Erfani, C. Leckie, J. Bailey, "Online cluster validity indices for performance monitoring of streaming data clustering," Int. J. Intell. Syst., pp. 1–23, 2018.
4. J. C. Dunn, "A fuzzy relative of the ISODATA process and its use in detecting compact well-separated clusters," J. Cybern., vol. 3, no. 3 , pp. 32–57, 1973.
5. J. C. Bezdek and N. R. Pal, "Some new indexes of cluster validity," IEEE Trans. Syst., Man, and Cybern., vol. 28, no. 3, pp. 301–315, Jun. 1998.
"""
mutable struct GD53 <: CVI
    label_map::LabelMap
    dim::Integer
    n_samples::Integer
    mu_data::RealVector     # dim
    n::IntegerVector        # dim
    v::RealMatrix           # dim x n_clusters
    CP::RealVector          # dim
    G::RealMatrix           # dim x n_clusters
    D::RealMatrix           # n_clusters x n_clusters
    inter::Float
    intra::Float
    n_clusters::Integer
    criterion_value::Float
end # GD53 <: CVI

"""
    GD53()

Default constructor for the Generalized Dunn's Index 53 (GD53) Cluster Validity Index.
"""
function GD53()
    GD53(
        LabelMap(),                     # label_map
        0,                              # dim
        0,                              # n_samples
        Array{Float, 1}(undef, 0),      # mu_data
        Array{Integer, 1}(undef, 0),    # n
        Array{Float, 2}(undef, 0, 0),   # v
        Array{Float, 1}(undef, 0),      # CP
        Array{Float, 2}(undef, 0, 0),   # G
        Array{Float, 2}(undef, 0, 0),   # D
        0.0,                            # inter
        0.0,                            # intra
        0,                              # n_clusters
        0.0                             # criterion_value
    )
end # GD53()

"""
    setup!(cvi::GD53, sample::Vector{T}) where {T<:RealFP}
"""
function setup!(cvi::GD53, sample::Vector{T}) where {T<:RealFP}
    # Get the feature dimension
    cvi.dim = length(sample)
    # Initialize the augmenting 2-D arrays with the correct feature dimension
    # NOTE: R is emptied and calculated in evaluate!, so it is not defined here
    cvi.v = Array{T, 2}(undef, cvi.dim, 0)
    cvi.G = Array{T, 2}(undef, cvi.dim, 0)
end # setup!(cvi::GD53, sample::Vector{T}) where {T<:RealFP}

function param_inc!(cvi::GD53, sample::RealVector, label::Integer)
    # Get the internal label
    i_label = get_internal_label!(cvi.label_map, label)

    n_samples_new = cvi.n_samples + 1
    if isempty(cvi.mu_data)
        mu_data_new = sample
        setup!(cvi, sample)
    else
        mu_data_new = (1 - 1/n_samples_new) .* cvi.mu_data + (1/n_samples_new) .* sample
    end

    if i_label > cvi.n_clusters
        n_new = 1
        v_new = sample
        CP_new = 0.0
        G_new = zeros(cvi.dim)
        if cvi.n_clusters == 0
            D_new = zeros(1,1)
        else
            D_new = zeros(cvi.n_clusters + 1, cvi.n_clusters + 1)
            D_new[1:cvi.n_clusters, 1:cvi.n_clusters] = cvi.D
            d_column_new = zeros(cvi.n_clusters + 1)
            for jx = 1:cvi.n_clusters
                d_column_new[jx] = (CP_new + cvi.CP[jx]) / (n_new + cvi.n[jx])
            end
            D_new[:, i_label] = d_column_new
            D_new[i_label, :] = transpose(d_column_new)
        end
        # Update 1-D parameters with a push
        cvi.n_clusters += 1
        push!(cvi.CP, CP_new)
        push!(cvi.n, n_new)
        # Update 2-D parameters with appending and reassignment
        cvi.v = [cvi.v v_new]
        cvi.G = [cvi.G G_new]
        cvi.D = D_new
    else
        n_new = cvi.n[i_label] + 1
        v_new = (1 - 1/n_new) .* cvi.v[:, i_label] + (1/n_new) .* sample
        delta_v = cvi.v[:, i_label] - v_new
        diff_x_v = sample .- v_new
        CP_new = cvi.CP[i_label] + transpose(diff_x_v)*diff_x_v + cvi.n[i_label]*transpose(delta_v)*delta_v + 2*transpose(delta_v)*cvi.G[:, i_label]
        G_new = cvi.G[:, i_label] + diff_x_v + cvi.n[i_label].*delta_v
        d_column_new = zeros(cvi.n_clusters)
        for jx = 1:cvi.n_clusters
            # Skip the current i_label index
            if jx == i_label
                continue
            end
            # d_column_new[jx] = sum((v_new - cvi.v[:, jx]).^2)
            d_column_new[jx] = (CP_new + cvi.CP[jx])/(n_new + cvi.n[jx])
        end
        # Update parameters
        cvi.n[i_label] = n_new
        cvi.v[:, i_label] = v_new
        cvi.CP[i_label] = CP_new
        cvi.G[:, i_label] = G_new
        cvi.D[:, i_label] = d_column_new
        cvi.D[i_label, :] = transpose(d_column_new)
    end
    cvi.n_samples = n_samples_new
    cvi.mu_data = mu_data_new
end # param_inc!(cvi::GD53, sample::RealVector, label::Integer)

function param_batch!(cvi::GD53, data::RealMatrix, labels::IntegerVector)
    cvi.dim, cvi.n_samples = size(data)
    # Take the average across all samples, but cast to 1-D vector
    cvi.mu_data = mean(data, dims=2)[:]
    # u = findfirst.(isequal.(unique(labels)), [labels])
    u = unique(labels)
    cvi.n_clusters = length(u)
    cvi.n = zeros(Integer, cvi.n_clusters)
    cvi.v = zeros(cvi.dim, cvi.n_clusters)
    cvi.CP = zeros(cvi.n_clusters)
    cvi.D = zeros(cvi.n_clusters, cvi.n_clusters)
    for ix = 1:cvi.n_clusters
        subset = data[:, findall(x->x==u[ix], labels)]
        cvi.n[ix] = size(subset, 2)
        cvi.v[1:cvi.dim, ix] = mean(subset, dims=2)
        diff_x_v = subset - cvi.v[:, ix] * ones(1, cvi.n[ix])
        cvi.CP[ix] = sum(diff_x_v.^2)
    end
    for ix = 1 : (cvi.n_clusters - 1)
        for jx = ix + 1 : cvi.n_clusters
            # cvi.D[jx, ix] = sum((cvi.v[:, ix] - cvi.v[:, jx]).^2)
            cvi.D[jx, ix] = (cvi.CP[ix] + cvi.CP[jx]) / (cvi.n[ix] + cvi.n[jx])
        end
    end
    cvi.D = cvi.D + transpose(cvi.D)
end # param_batch!(cvi::GD53, data::RealMatrix, labels::IntegerVector)

function evaluate!(cvi::GD53)
    cvi.intra = 2*maximum(cvi.CP ./ cvi.n)
    if cvi.n_clusters > 1
        # Between-group measure of separation/isolation
        cvi.inter = minimum(cvi.D[triu(ones(Bool, cvi.n_clusters, cvi.n_clusters), 1)])
        # GD53 index value
        cvi.criterion_value = cvi.inter/cvi.intra
    else
        cvi.criterion_value = 0.0
    end
end # evaluate(cvi::GD53)


text:
"""
    couple_models(aero::Peters, stru::LiftingLineSection)

Create an aerostructural model using a using the unsteady aerodynamic model
defined by Peters et al. and a lifting line section model.  The existence of this
coupling allows [`Peters`](@ref) to be used with [`LiftingLine`](@ref).  This
model introduces the freestream air density ``\\rho`` as an additional parameter.
"""
couple_models(aero::Peters, stru::LiftingLineSection) = (aero, stru)

# --- Traits --- #

number_of_additional_parameters(::Type{<:Peters}, ::Type{LiftingLineSection}) = 1

coupling_inplaceness(::Type{<:Peters}, ::Type{LiftingLineSection}) = OutOfPlace()

coupling_rate_jacobian_type(::Type{<:Peters}, ::Type{LiftingLineSection}) = Linear()
coupling_state_jacobian_type(::Type{<:Peters}, ::Type{LiftingLineSection}) = Nonlinear()
coupling_parameter_jacobian_type(::Type{<:Peters}, ::Type{LiftingLineSection}) = Nonlinear()
coupling_time_gradient_type(::Type{<:Peters}, ::Type{LiftingLineSection}) = Zeros()

# --- Methods --- #

function get_coupling_inputs(aero::Peters{N,TF,SV,SA}, stru::LiftingLineSection,
    dx, x, p, t) where {N,TF,SV,SA}

    # extract rate variables
    dλ = dx[SVector{N}(1:N)]
    dvx, dvy, dvz, dωx, dωy, dωz = dx[SVector{6}(N+1:N+6)]

    # extract state variables
    λ = x[SVector{N}(1:N)]
    vx, vy, vz, ωx, ωy, ωz = x[SVector{6}(N+1:N+6)]

    # extract parameters
    a, b, a0, α0, cd0, cm0, ρ = p

    # extract model constants
    bbar = aero.b

    # freestream velocity components
    u, v, ω = liftingline_velocities(vx, vz, ωy)
    udot, vdot, ωdot = liftingline_accelerations(dvx, dvz, dωy)

    # aerodynamic loads
    Na, Aa, Ma = peters_loads(a, b, ρ, a0, α0, cd0, cm0, bbar, u, v, ω, vdot, ωdot, λ)

    # loads per unit span
    f = SVector(Aa, 0, Na)
    m = SVector(0, Ma, 0)

    # return inputs
    return SVector(u, ω, vdot, ωdot, f..., m...)
end

# --- Performance Overloads --- #

#TODO: State rate, state, and parameter jacobians

# --- Convenience Methods --- #

function set_additional_parameters!(padd, aero::Peters, stru::LiftingLineSection; rho)

    padd[1] = rho

    return padd
end

function separate_additional_parameters(aero::Peters, stru::LiftingLineSection, padd)

    return (rho = padd[1],)
end

ext: jl


text:
export AbstractAgent

"""
An agent is a functional object, which takes in an observation and returns an action. Agents must also implement the `update!(agent::AbstractAgent, obs_action::Pair)` method to indicate how to update the internal state of the agent.
"""
abstract type AbstractAgent end


text:
function make_gui_block(cur_block, cur_symbols)
  cur_func = gensym()

  cur_lambda = Expr(
    :(->),
    Expr(:tuple, cur_symbols...),
    cur_block
  )

  quote
    $(cur_func) = $(cur_lambda)

    cur_observer = SimplePlots.Observables.Observable{Any}(
      $(cur_func)(
        $(map(
          cur_symbol -> :( observe( $(cur_symbol) )[] ),
          cur_symbols
        )...)
      )
    )

    map!(
      $(cur_func), cur_observer, $(map(
        cur_symbol -> :( observe( $(cur_symbol) ) ),
        cur_symbols
      )...)
    )

    cur_observer
  end
end


text:
using NipponColors, Test

plot_color()
dict_color(; mode = :hex)
dict_color(; mode = :rgb)


text:
using WindMixing
using BSON
using OceanParameterizations
using Flux
using OrdinaryDiffEq, DiffEqSensitivity

train_files = ["strong_wind"]

PATH = pwd()
OUTPUT_PATH = joinpath(PATH, "Output")
FILE_PATH = joinpath("D:\\University Matters\\Massachusetts Institute of Technology\\CLiMA Project\\OceanParameterizations.jl\\training_output", "test.jld2")

𝒟train = data(train_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=true)

uw_NN_model = BSON.load(joinpath(PATH, "NDEs", "uw_NDE_SWNH_100.bson"))[:neural_network]
vw_NN_model = BSON.load(joinpath(PATH, "NDEs", "vw_NDE_SWNH_100.bson"))[:neural_network]
wT_NN_model = BSON.load(joinpath(PATH, "NDEs", "wT_NDE_SWNH_100.bson"))[:neural_network]

train_epochs = [100]
train_tranges = [1:5:100]
opts = [ADAM(0.01)]

write_metadata_NDE_training(FILE_PATH, train_files, train_epochs, train_tranges, opts, uw_NN_model, vw_NN_model, wT_NN_model)

train_NDE_convective_adjustment(uw_NN_model, vw_NN_model, wT_NN_model, 𝒟train, 1:5:100, ROCK4(), [ADAM(0.01)], 100, OUTPUT_PATH, "weights_SWNH_convective_adjustment_100")


text:
using Distances, Statistics
using MultivariateStats
using PyPlot
using WordTokenizers
using TextAnalysis
using DelimitedFiles
using Embeddings
using Parameters: @with_kw
@with_kw mutable struct Args
    path::String = "data/pg345.txt" # path to data
    string::String = "my favorite food is strawberry ice cream" # Max length of the doc to be considered
end
const embtable = load_embeddings(GloVe{:en},1) # or load_embeddings(FastText_Text) or ...
#Function to return the index of the word in the embedding (returns 0 if the word is not found)
const get_word_index = Dict(word=>ii for (ii,word) in enumerate(embtable.vocab))
embeddings = embtable.embeddings
vocab = embtable.vocab
vec_size, vocab_size = size(embeddings)
println("Loaded embeddings, each word is represented by a vector with $vec_size features. The vocab size is $vocab_size")
vec_idx(s) = findfirst(x -> x == s, vocab)
function vec(s)
    if vec_idx(s)!= nothing
        embeddings[:, vec_idx(s)]
    end
end
# To compare the similarity of document vectors
cosine(x,y) = 1 - cosine_dist(x, y)
function closest(v, n=20)
    list = [(x,cosine(embeddings'[x,:], v)) for x in 1:size(embeddings)[2]]
    topn_idx = sort(list, by = x -> x[2], rev = true)[1:n]
    return [vocab[a] for (a,_) in topn_idx]
end
# Function to get a mean sentence vector
function sentvec(s,sentences)
    local arr = []
    for w in split(sentences[s])
        if vec(w) != nothing
            push!(arr, vec(w))
        end
    end
    if length(arr) == 0
        ones(Float32, (50,1))*999
    else
        mean(arr)
    end
end
function closest_sent(input_str, n=20)
    mean_vec_input = mean([vec(w) for w in split(input_str)])
    list = [(x,cosine(mean_vec_input, sentvec(x))) for x in 1:length(sentences)]
    topn_idx=sort(list, by = x -> x[2], rev = true)[1:n]
    return [sentences[a] for (a,_) in topn_idx]
end
function match_docs(; kws...)
    args = Args(; kws...)
    txt = open(args.path) do file
        read(file, String)
    end
    println("Loaded $(args.path), length=$(length(txt)) characters")
    txt = replace(txt, r"\n|\r|_|,|—" => " ")
    txt = replace(txt, r"[\"""*();!]" => "")
    sd = StringDocument(txt)
    prepare!(sd, strip_whitespace)
    global sentences = split_sentences(sd.text)
    i = 1
    for s in 1:length(sentences)
        if length(split(sentences[s]))>3
            sentences[i]=lowercase(replace(sentences[s], "."=>""))
            i += 1
        end
    end
    save_vector(sentences)
    return closest_sent(args.string)
end
function save_vector(sentences)
    i = length(sentences) + 1
    sent_vecs = []
    for s in 1:length(sentences)
        i == 1 ? sent_vecs=sentvec(s) : push!(sent_vecs,sentvec(s))
    end
    writedlm( "data/sent_vec.csv",  sent_vecs, ',')
    writedlm( "data/sentences.csv",  sentences, ',')
end
function closest_sent_pretrained(pretrained_arr, input_str, n=20)
    mean_vec_input = mean([vec(w) for w in split(input_str)])
    list = [(x,cosine(mean_vec_input, pretrained_arr[x,:])) for x in 1:length(sentences)]
    topn_idx = sort(list, by = x -> x[2], rev=true)[1:n]
    return [sentences[a] for (a,_) in topn_idx]
end
function load_and_match(sentences_path, sent_vecs_path, string)
    sentences = readdlm("sentences.csv", '!', String, header=false)
    sent_vecs = readdlm("sent_vec.csv", ',', Float32, header=false)
    return closest_sent_pretrained(sent_vecs, string)
end

max_stars_repo_name: akaysh/Rango

max_issues_repo_name: akaysh/Rango

max_forks_repo_name: akaysh/Rango


text:
"""
    vtk_mesh_file(mesh::GenericMesh[, Ω::Domain], name::String)

Creates a `VTK` file (.vtu) with name `name` containing the mesh information.
It is possible to export only a `Domain` (i.e. only a part of the mesh).

Note that all the heavy lifting is done by the package `WriteVTK`.
We refer to its documentation for more information.

Warning: the output is not an actual file (on disk).
To save it, simply write:

    vtk_mesh_file(mesh, "my_mesh") |> vtk_save

To add some data (scalar or vector values at the mesh nodes or mesh cells)
for visualization purposes, you can do for instance:

    vtkfile = vtk_mesh_file(mesh, name)
    vtkfile["my_point_data", VTKPointData()] = pdata
    vtkfile["my_cell_data", VTKCellData()] = cdata
    vtk_save(vtkfile)

It is possible also to export a partition `Ωs::Vector{Domain}` using
`Multiblock` files (.vtm), for instance like so

    vtmfile = vtk_multiblock(name)
    for (Ω, pdata) in zip(Ωs, pdatas)
        vtkfile = vtk_mesh_file(mesh, Ω, name)
        vtkfile["my_point_data", VTKPointData()] = pdata
    end
    vtk_save(vtmfile)

To save a sequence of solutions (time steps, iterations), simply append
the number of the element to the file name.
Paraview will recognize the sequence automatically.
"""
function vtk_mesh_file(mesh::GenericMesh, name::String)
    points = _vtk_points(mesh)
    cells = _vtk_cells(mesh)
    return vtk_grid(name * ".vtu", points, cells)
end
function vtk_mesh_file(mesh::GenericMesh, Ω::Domain, name::String)
    points = _vtk_points(mesh)
    cells = _vtk_cells(mesh, Ω)
    return vtk_grid(name * ".vtu", points, cells)
end

"""
    _vtk_points(mesh::GenericMesh)

Creates the matrix of nodes in the format required by `WriteVTK`.
"""
function _vtk_points(mesh::GenericMesh)
    vtx = zeros(Float64, 3, length(nodes(mesh)))
    for (i, nd) in enumerate(nodes(mesh))
        vtx[1:ambient_dimension(mesh),i] = nd
    end
    return vtx
end

"""
    _vtk_cells(mesh::GenericMesh, E::DataType)
    _vtk_cells(mesh::GenericMesh, Ω::Domain)

Creates the vector of all elements contained in the mesh in the format
required by `WriteVTK` for a particular element type `E<:AbstractElement`
or a `Domain` instance.
"""
function _vtk_cells end

function _vtk_cells(tags, E::DataType)
    vtk_cell_type, ind = etype_to_vtk_cell_type[E]
    return [MeshCell(vtk_cell_type, tags[ind, i]) for i in 1:size(tags, 2)]
end
function _vtk_cells(mesh::GenericMesh, E::DataType)
    tags = elements(mesh)[E]
    return _vtk_cells(tags, E)
end
function _vtk_cells(mesh::GenericMesh, Ω::Domain)
    cells = MeshCell[]
    # Loop on `ElementaryEntity`
    for ω in Ω
        # Loop on `AbstractElement`
        for (E, ind) in ent2tags(mesh)[ω]
            # Subset corresponding to the `ElementaryEntity`
            tags = elements(mesh)[E][:,ind]
            append!(cells,  _vtk_cells(tags, E))
        end
    end
    return cells
end
function _vtk_cells(mesh::GenericMesh)
    cells = MeshCell[]
    # Loop on `AbstractElement`
    for E in etypes(mesh)
        # Export only the cells of the largest geometrical dimension
        if geometric_dimension(E) == ambient_dimension(mesh)
            append!(cells,  _vtk_cells(mesh, E))
        end
    end
    return cells
end

"""
    const etype_to_vtk_cell_type

OrderedDictionary mapping internal element types to a tuple containing:
    - the corresponding `WriteVTK` cell types (following the convention
    chosen by `VTK`, see below);
    - the indices in the `elements` column that defines the element.
    This is because we want to support the export of more than just the flat
    elements available in the `VTK` specification, hence which may require
    a conversion of some sort.

See VTK specification [Fig. 2]
http://www.vtk.org/VTK/img/file-formats.pdf

VTK_VERTEX (=1)
VTK_POLY_VERTEX (=2)
VTK_LINE (=3)
VTK_POLY_LINE (=4)
VTK_TRIANGLE (=5)
VTK_TRIANGLE_STRIP (=6)
VTK_POLYGON (=7)
VTK_PIXEL (=8)
VTK_QUAD (=9)
VTK_TETRA (=10)
VTK_VOXEL (=11)
VTK_HEXAHEDRON (=12)
VTK_WEDGE (=13)
VTK_PYRAMID (=14)
"""
const etype_to_vtk_cell_type = OrderedDict(
    SVector{3,Float64} => (VTKCellTypes.VTK_VERTEX, collect(1:1)),
    LagrangeLine{2,SVector{3,Float64}} => (VTKCellTypes.VTK_LINE, collect(1:2)),
    LagrangeTriangle{3,SVector{2,Float64}} => (VTKCellTypes.VTK_TRIANGLE, collect(1:3)),
    LagrangeTriangle{3,SVector{3,Float64}} => (VTKCellTypes.VTK_TRIANGLE, collect(1:3)),
    LagrangeTetrahedron{4,SVector{3,Float64}} => (VTKCellTypes.VTK_TETRA, collect(1:4)),
)

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
using Oceananigans
using Oceananigans.Grids
using JLD2
import Oceananigans.Fields: interpolate
using Oceananigans.Fields: fractional_indices, _interpolate

"""
    interpolate(field, LX, LY, LZ, grid, x, y, z)
Interpolate `field` to the physical point `(x, y, z)` using trilinear interpolation. The location of
the field is specified with `(LX, LY, LZ)` and the field is defined on `grid`.
Note that this is a lower-level `interpolate` method defined for use in CPU/GPU kernels.
"""
@inline function interpolate(field, LX, LY, LZ, grid, x, y, z)
    i, j, k = fractional_indices(x, y, z, (LX, LY, LZ), grid)

    # We use mod and trunc as CUDA.modf is not defined.
    # For why we use Base.unsafe_trunc instead of trunc see:
    # https://github.com/CliMA/Oceananigans.jl/issues/828
    # https://github.com/CliMA/Oceananigans.jl/pull/997
    ξ, i = mod(i, 1), Base.unsafe_trunc(Int, i)
    η, j = mod(j, 1), Base.unsafe_trunc(Int, j)
    ζ, k = mod(k, 1), Base.unsafe_trunc(Int, k)

    return _interpolate(field, ξ, η, ζ, i+1, j+1, k+1)
end

function interpolate(field::AbstractArray, LX, LY, LZ, grid, x::AbstractArray, y::AbstractArray, z::AbstractArray)
    newgrid = zeros((length(x), length(y), length(z)))
    for k in eachindex(z), j in eachindex(y), i in eachindex(x)
        newgrid[i,j,k] = interpolate(field, LX, LY, LZ, grid, x[i], y[j], z[k])
    end
    return newgrid
end

function getgrid(file)
    coarsetobe = jldopen(file)
    grid = coarsetobe["grid"]
    close(coarsetobe)
    return grid
end

function interpolate(states::Array{Array{S, 3},1}, grid, newgrid) where {S}
    newx, newy, newz = nodes((Cell, Cell, Cell), new_grid, reshape=true)
    interpolatedstates = [interpolate(states, LX, LY, LZ, grid, newx, newy, newz) for state in states]
    return interpolatedstates
end
##
#=
# Example
files = [pwd() * "/Channel_1_checkpoint_iteration404905.jld2",
         pwd() * "/Channel_2_checkpoint_iteration728850.jld2",
         pwd() * "/Channel_3_checkpoint_iteration1087211.jld2",
         pwd() * "/Channel_4_checkpoint_iteration6160221.jld2",
         pwd() * "/Channel_8_checkpoint_iteration3164301.jld2",
         pwd() * "/Channel_16_checkpoint_iteration6317902.jld2",
         pwd() * "/Channel_24_checkpoint_iteration8612852.jld2",
         pwd() * "/Channel_32_checkpoint_iteration1272125.jld2"
]

file = files[end]
coarsetobe = jldopen(file)
grid = coarsetobe["grid"]
u = coarsetobe["velocities"]["u"]
LX, LY, LZ = typeof.(u["location"])
array = u["data"]
# newpoint = interpolate(array, LX, LY, LZ, grid, 0.0, 0.0, 00.0)
i, j, k = fractional_indices(0, 0, -2000, (LX, LY, LZ), grid)
##

topo = (Periodic, Bounded, Bounded)
domain = (x=(0, 1e6), y=(0, 1e6), z=(-3000, 0))
new_grid = RegularCartesianGrid(topology=topo, size=(150, 150, 12); domain...)
newx, newy, newz = nodes((Cell, Cell, Cell), new_grid, reshape=true)

field = u["data"]
newgrid = interpolate(field, LX, LY, LZ, grid, newx, newy, newz)

# newgridvals = interpolate.(Ref(array), Ref(LX), Ref(LY), Ref(LZ), Ref(grid),newx, newy, newx)

scene = volume(0..1, 0..1, 0..1, newgrid)
volume!(scene, 1..2, 0..1, 0..1, field)
=#

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
using PyPlot

s = im
z = 0.5
n = 20
k = 1:n
w = z.^k ./ k.^s

figure(1)
plot( real(w), imag(w), "o-")
axis("equal")
xlim([-0.5, 0.5])


text:
#=----------------------------------------------------------------------------

Perceptually Uniform Colour Maps

Copyright (c) 2015 Peter Kovesi
pk@peterkovesi.com

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

The Software is provided "as is", without warranty of any kind.

PK October 2015

----------------------------------------------------------------------------=#

"""
**PerceptualColourMaps**

Perceptually Uniform Colour maps

Peter Kovesi

[peterkovesi.com](http://peterkovesi.com)

*Colour Map Generation and Application*

* cmap - Library of perceptually uniform colour maps.
* equalisecolourmap - Equalise colour contrast over a colour map.
* linearrgbmap - Linear rgb colourmap from black to a specified colour.
* applycolourmap - Applies colourmap to a single channel image to obtain an RGB result.
* applycycliccolourmap - Applies a cyclic colour map to an image of angular data.
* applydivergingcolourmap - Applies a diverging colour map to an image.
* ternaryimage - Perceptualy uniform ternary image from 3 bands of data.
* relief - Generates a relief shaded image.
* viewlabspace - Visualisation of Lab colour space.

*Images for testing colour maps*

* sineramp - Generates sine on a ramp colour map test image.
* circlesineramp - Generates a test image for evaluating cyclic colour maps.

*Utilities*

* histtruncate - Truncates ends of an image histogram.
* normalise - Normalises image values to 0-1, or to desired mean and variance.
* srgb2lab - Convert RGB colour map or RGB image to Lab.
* lab2srgb - Convert Lab colour map or Lab image to RGB.

*Reference:*

[Good Colour Maps: How to Design Them.  arXiv:1509.03700 [cs.GR] 2015.](http://arxiv.org/abs/1509.03700)

"""
module PerceptualColourMaps

using Printf, LinearAlgebra, Statistics
using Images

include("cmap.jl")
include("applycolourmaps.jl")
include("viewlabspace.jl")
include("relief.jl")
include("utilities.jl")

end  # module


text:
macro Σein!_str(s::AbstractString)
    tensor_lhs_sym, tensors_rhs_sym, ex_rhs_lhs, ex_rhs_rhs = parser_strmac(s)
    ex_rhs =
        macroexpand(TensorOperations, :(TensorOperations.@tensor $ex_rhs_lhs = $ex_rhs_rhs))
    ex_lhs = Expr(:tuple, tensor_lhs_sym, tensors_rhs_sym...)
    return :($ex_lhs -> $ex_rhs)
end

macro Σein_str(s::AbstractString)
    tensor_lhs_sym, tensors_rhs_sym, ex_rhs_lhs, ex_rhs_rhs = parser_strmac(s)
    ex_rhs = macroexpand(
        TensorOperations,
        :(TensorOperations.@tensor $ex_rhs_lhs := $ex_rhs_rhs),
    )
    ex_lhs = Expr(:tuple, tensors_rhs_sym...)
    return :($ex_lhs -> $ex_rhs)
end

function parser_strmac(s::AbstractString)
    tensors_rhs = parse_tensors_rhs(s)
    indices_lhs = parse_indices_lhs(s)
    pairs_index = parse_pairs_index(s)
    tree_contract = parse_tree_contract(s)
    for ((kind, klab), v) in pairs_index
        replace!(tensors_rhs[klab].indices, (kind, klab) => v)
    end
    indices_all = union([x.indices for x in tensors_rhs]...)
    indices_sym = Dict(x => gensym() for x in indices_all)
    tensors_rhs_sym = [gensym() for _ = 1:length(tensors_rhs)]
    tensors_rhs_ex = Vector{Expr}()
    for (i, a) in enumerate(tensors_rhs)
        ex = Expr(:ref, tensors_rhs_sym[i], [indices_sym[i] for i in a.indices]...)
        ex = a.isconj ? Expr(:call, :conj, ex) : ex
        push!(tensors_rhs_ex, ex)
    end
    tensor_lhs_sym = gensym()
    ex_rhs_lhs = Expr(:ref, tensor_lhs_sym, [indices_sym[i] for i in indices_lhs]...)
    ex_rhs_rhs = if isnothing(tree_contract)
        Expr(:call, :(*), tensors_rhs_ex...)
    else
        evaltree(i, j) = Expr(
            :call,
            :(*),
            i isa Int ? tensors_rhs_ex[i] : evaltree(i...),
            j isa Int ? tensors_rhs_ex[j] : evaltree(j...),
        )
        evaltree(i, j, k...) = evaltree((i, j), k...)
        evaltree(tree_contract...)
    end
    return tensor_lhs_sym, tensors_rhs_sym, ex_rhs_lhs, ex_rhs_rhs
end


text:
module Impero

using LinearAlgebra, SymbolicUtils

include("abstract_core.jl")
include("symbolic_utils_hooks.jl")
include("tree_plotting.jl")

end # module


text:
module ReactiveMPDistributionsTest

using Test
using ReactiveMP

@testset "Distributions" begin

end

end


text:
# ---
# title: circsine
# id: demo_circsine
# cover: assets/circsine.png
# author: Peter Kovesi
# date: 2018-10-26
# ---

using Images
using ImagePhaseCongruency

## Circular features at a phase congruent angle of pi/4 and
## an amplitude decay exponent of 1.5
img = circsine(offset = pi/4, ampexponent = -1.5)
adjust_histogram(Gray.(img), LinearStretching())

# save cover image #src
isdir("assets") || mkdir("assets") #src
save(joinpath("assets", "circsine.png"), adjust_histogram(Gray.(img), LinearStretching())) #src

max_stars_repo_name: peterkovesi/ImagePhaseCongruency

max_issues_repo_name: peterkovesi/ImagePhaseCongruency

max_forks_repo_name: peterkovesi/ImagePhaseCongruency


text:
"""
Specify the folder for the Literate scripts, by default this is `scripts/`.
"""
function literate_folder(rp::String="")
    isempty(rp) && return PATHS[:literate]
    path = joinpath(PATHS[:folder], rp)
    !isdir(path) && error("Specified Literate path not found ($rp -- $path)")
    PATHS[:literate] = path
    return path
end

#
# Convenience functions to work with strings and substrings
#

"""
    str(s)

Returns the string corresponding to `s`: `s` itself if it is a string, or the
parent string if `s` is a substring. Do not confuse with `String(s::SubString)`
which casts `s` into its own object.

# Example

```julia-repl
julia> a = SubString("hello Fraknlin", 3:8);
julia> Franklin.str(a)
"hello Franklin"
julia> String(a)
"llo Fr"
```
"""
str(s::String)::String    = s
str(s::SubString)::String = s.string

"""
    subs(s, from, to)
    subs(s, from)
    subs(s, range)
    subs(s)

Convenience functions to take a substring of a string.

# Example
```julia-repl
julia> Franklin.subs("hello", 2:4)
"ell"
```
"""
subs(s::AS, from::Int, to::Int)::SubString    = SubString(s, from, to)
subs(s::AS, from::Int)::SubString             = subs(s, from, from)
subs(s::AS, range::UnitRange{Int})::SubString = SubString(s, range)
subs(s::AS) = SubString(s)

"""
$(SIGNATURES)

Given a substring `ss`, returns the position in the parent string where the substring starts.

# Example
```julia-repl
julia> ss = SubString("hello", 2:4); Franklin.from(ss)
2
```
"""
from(ss::SubString)::Int = nextind(str(ss), ss.offset)
from(s::String) = 1

"""
$(SIGNATURES)

Given a substring `ss`, returns the position in the parent string where the substring ends.

# Example
```julia-repl
julia> ss = SubString("hello", 2:4); Franklin.to(ss)
4
```
"""
to(ss::SubString)::Int = ss.offset + ss.ncodeunits
to(s::String) = lastindex(s)

"""
$(SIGNATURES)

Returns the string span of a regex match. Assumes there is no unicode in the match.

# Example
```julia-repl
julia> Franklin.matchrange(match(r"ell", "hello"))
2:4
```
"""
matchrange(m::RegexMatch)::UnitRange{Int} = m.offset .+ (0:(length(m.match)-1))

# Other convenience functions

"""
$(SIGNATURES)

Convenience function to display a time since `start`.
"""
function time_it_took(start::Float64)
    comp_time = time() - start
    mess = comp_time > 60 ? "$(round(comp_time/60;   digits=1))m" :
           comp_time >  1 ? "$(round(comp_time;      digits=1))s" :
                            "$(round(comp_time*1000; digits=1))ms"
    return "[done $(lpad(mess, 6))]"
end

```
$(SIGNATURES)

Nicer printing of processes.
```
function print_final(startmsg::AS, starttime::Float64)::Nothing
    tit = time_it_took(starttime)
    rprint("$startmsg$tit")
    print("\n")
end

"""
$(SIGNATURES)

Convenience function to denote a string as being in a math context in a recursive parsing
situation. These blocks will be processed as math blocks but without adding KaTeX elements to it
given that they are part of a larger context that already has KaTeX elements.
NOTE: this happens when resolving latex commands in a math environment. So for instance if you have
`\$\$ x \\in \\R \$\$` and `\\R` is defined as a command that does `\\mathbb{R}` well that would be
an embedded math environment. These environments are marked as such so that we don't add additional
KaTeX markers around them.
"""
mathenv(s::AS)::String = "_\$>_$(s)_\$<_"

"""
$(SIGNATURES)

Takes a string `s` and replace spaces by underscores so that that we can use it
for hyper-references. So for instance `"aa  bb"` will become `aa_bb`.
It also defensively removes any non-word character so for instance `"aa bb !"` will be `"aa_bb"`
"""
function refstring(s::AS)::String
    # remove html tags
    st = replace(s, r"<[a-z\/]+>"=>"")
    # remove non-word characters
    st = replace(st, r"&#[0-9]+;" => "")
    st = replace(st, r"[^\p{L}0-9_\-\s]" => "")
    # replace spaces by dashes
    st = replace(lowercase(strip(st)), r"\s+" => "_")
    # to avoid clashes with numbering of repeated headers, replace
    # double underscores by a single one (see convert_header function)
    st = replace(st, r"__" => "_")
    # in the unlikely event we don't have anything here, return the hash
    # of the original string
    isempty(st) && return string(hash(s))
    return st
end

"""
    context(parent, position)

Return an informative message of the context of a position and where the
position is, this is useful when throwing error messages.
"""
function context(par::AS, pos::Int)
    # context string
    lidx = lastindex(par)
    if pos > 20
        head = max(1, prevind(par, pos-20))
    else
        head = 1
    end
    if pos <= lidx-20
        tail = min(lidx, nextind(par, pos+20))
    else
        tail = lidx
    end
    prepend  = ifelse(head > 1, "...", "")
    postpend = ifelse(tail < lidx, "...", "")

    ctxt = prepend * subs(par, head, tail) * postpend

    # line number
    lines  = split(par, "\n", keepempty=false)
    nlines = length(lines)
    ranges = zeros(Int, nlines, 2)
    cs = 0
    for (i, l) in enumerate(lines[1:end-1])
        tmp = [nextind(par, cs), nextind(par, lastindex(l) + cs)]
        ranges[i, :] .= tmp
        cs = tmp[2]
    end
    ranges[end, :] = [nextind(par, cs), lidx]

    lno = findfirst(i -> ranges[i,1] <= pos <= ranges[i,2], 1:nlines)

    # Assemble to form a message
    mess = """
    Context:
    \t$(strip(ctxt)) (near line $lno)
    \t$(" "^(pos-head+length(prepend)))^---
    """
    return mess
end

"""
    rprint(s)

Print an overwriting line being aware of the display width.
"""
function rprint(s::AS)::Nothing
    dwidth  = displaysize(stdout)[2]

    padded_s = rpad(s, dwidth)
    trunc_s  = padded_s[1:prevind(padded_s, min(dwidth, lastindex(padded_s)))]
    padded_s = rpad(trunc_s, dwidth)

    print("\r$padded_s")
    return nothing
end

"""
    check_ping(ipaddr)

Try a single ping to an address `ipaddr`. There's a timeout of 1s.
"""
function check_ping(ipaddr)
    # portable count (creds https://docdave.science/writing-ping-in-julia/)
    opt = ifelse(Sys.iswindows(), "-n", "-c")
    return success(`ping $opt 1 -t 1 $ipaddr`)
end

"""
    invert_dict(dict)

Invert a dictionary i.e transform a=>[1,2],b=>[1] to 1=>[a,b], 2=>[a]
"""
function invert_dict(dict)
    fe = first(dict)
    TK = typeof(fe.first)
    TV = Vector{eltype(fe.second)}
    inv_dict = LittleDict{TK, Vector{eltype(dict).types[1]}}()
    for (key, val) in dict
        for nkey in val
            if haskey(inv_dict, nkey)
                push!(inv_dict[nkey], key)
            else
                inv_dict[nkey] = [key]
            end
        end
    end
    return inv_dict
end


text:
# Not using `@syms` in this file to be more flexible (allowing vargs) and avoid using `@doc`.
# Use macroexpand(Main, :(@syms f(x))) to see what @syms would generate.

"""
    P(x)
    P(x, y)
    P(vargs...)

A named variable to denote probability.
"""
const P = SymbolicUtils.Sym{SymbolicUtils.FnType{Tuple, Number}}(:P)

"""
    Do(x)

Do-operator.
Pearl denotes this operator with `do`.
However, Julia's parser already uses `do`, so `Do` seemed to be the most natural alternative.
"""
const Do = SymbolicUtils.Sym{(SymbolicUtils.FnType){Tuple{Number}, Number}}(:Do)

"""
    u ¦ v

`u` given `v`.
Using the broken bar (¦) instead of the vertical line (|), because the vertical line is already defined by Julia base.

Probably, the vertical line could be used by defining it for sets.
We can just use comma for P(A, B)
"""
const (¦) = SymbolicUtils.Sym{(SymbolicUtils.FnType){Tuple{Number, Number}, Number}}(:¦)


text:
# Autogenerated wrapper script for Qt5Tools_jll for armv7l-linux-gnueabihf-cxx03
export qdbus, qtdiag, qtpaths, qtplugininfo

using Glib_jll
using Qt5Base_jll
JLLWrappers.@generate_wrapper_header("Qt5Tools")
JLLWrappers.@declare_executable_product(qdbus)
JLLWrappers.@declare_executable_product(qtdiag)
JLLWrappers.@declare_executable_product(qtpaths)
JLLWrappers.@declare_executable_product(qtplugininfo)
function __init__()
    JLLWrappers.@generate_init_header(Glib_jll, Qt5Base_jll)
    JLLWrappers.@init_executable_product(
        qdbus,
        "bin/qdbus",
    )

    JLLWrappers.@init_executable_product(
        qtdiag,
        "bin/qtdiag",
    )

    JLLWrappers.@init_executable_product(
        qtpaths,
        "bin/qtpaths",
    )

    JLLWrappers.@init_executable_product(
        qtplugininfo,
        "bin/qtplugininfo",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=
This object represents a Telegram user or bot.

    Field       Type	   Description
    id	        Integer	   Unique identifier for this user or bot
    first_name	String	   User‘s or bot’s first name
    last_name	String	   Optional. User‘s or bot’s last name
    username	String	   Optional. User‘s or bot’s username
=#

type User
    __id::Int32
    __first_name::UTF8String
    __last_name::UTF8String
    __username::UTF8String
end

#=
This object represents a chat.

    Field	    Type      Description
    id	        Integer	  Unique identifier for this chat, not exceeding 1e13 by absolute value
    type	    String	  Type of chat, can be either “private”, “group”, “supergroup” or “channel”
    title	    String	  Optional. Title, for channels and group chats
    username	String	  Optional. Username, for private chats and channels if available
    first_name	String	  Optional. First name of the other party in a private chat
    last_name	String	  Optional. Last name of the other party in a private chat
=#

type Chat
    __id::Int32
    __type::UTF8String
    __title::UTF8String
    __username::UTF8String
    __first_name::UTF8String
    __last_name::UTF8String
end

#=
This object represents one special entity in a text message. For example, hashtags, usernames, URLs, etc.

    Field	  Type	       Description
    type	  String	   Type of the entity. One of mention (@username), hashtag, bot_command, url, email, bold (bold text), italic (italic text), code (monowidth string), pre (monowidth block), text_link (for clickable text URLs)
    offset	  Integer	   Offset in UTF-16 code units to the start of the entity
    length	  Integer	   Length of the entity in UTF-16 code units
    url	      String	   Optional. For “text_link” only, url that will be opened after user taps on the text
=#

type MessageEntity
    __type::UTF8String
    __offset::Int32
    __length::Int32
    __url::UTF8String
end

#=
This object represents an audio file to be treated as music by the Telegram clients.

    Field	      Type	    Description
    file_id	      String	Unique identifier for this file
    duration	  Integer	Duration of the audio in seconds as defined by sender
    performer	  String	Optional. Performer of the audio as defined by sender or by audio tags
    title	      String	Optional. Title of the audio as defined by sender or by audio tags
    mime_type	  String	Optional. MIME type of the file as defined by sender
    file_size	  Integer	Optional. File size
=#

type Audio
    __file_id::UTF8String
    __duration::Int32
    __performer::UTF8String
    __title::UTF8String
    __mime_type::UTF8String
    __file_size::Int32
end

#=
This object represents one size of a photo or a file / sticker thumbnail.

    Field	      Type	    Description
    file_id	      String	Unique identifier for this file
    width	      Integer	Photo width
    height	      Integer	Photo height
    file_size	  Integer	Optional. File size
=#

type PhotoSize
    __file_id::UTF8String
    __width::Int32
    __height::Int32
    __file_size::Int32
end

#=
This object represents a general file (as opposed to photos, voice messages and audio files).

    Field	      Type	       Description
    file_id	      String	   Unique file identifier
    thumb	      PhotoSize	   Optional. Document thumbnail as defined by sender
    file_name	  String	   Optional. Original filename as defined by sender
    mime_type	  String	   Optional. MIME type of the file as defined by sender
    file_size	  Integer	   Optional. File size
=#

type Document
    __file_id::UTF8String
    __thumb::PhotoSize
    __file_name::UTF8String
    __mime_type::UTF8String
    __file_size::Int32
end

#=
This object represents a sticker.

    Field	      Type	       Description
    file_id	      String	   Unique identifier for this file
    width	      Integer	   Sticker width
    height	      Integer	   Sticker height
    thumb	      PhotoSize	   Optional. Sticker thumbnail in .webp or .jpg format
    file_size	  Integer	   Optional. File size
=#

type Sticker
    __file_id::UTF8String
    __width::Int32
    __height::Int32
    __thumb::PhotoSize
    __file_size::Int32
end

#=
This object represents a video file.

    Field	      Type	       Description
    file_id	      String	   Unique identifier for this file
    width	      Integer	   Video width as defined by sender
    height	      Integer	   Video height as defined by sender
    duration	  Integer	   Duration of the video in seconds as defined by sender
    thumb	      PhotoSize	   Optional. Video thumbnail
    mime_type	  String	   Optional. Mime type of a file as defined by sender
    file_size	  Integer	   Optional. File size
=#

type Video
    __file_id::UTF8String
    __width::Int32
    __height::Int32
    __duration::Int32
    __thumb::PhotoSize
    __mime_type::UTF8String
    __file_size::Int32
end

#=
This object represents a voice note.

    Field	      Type	    Description
    file_id	      String	Unique identifier for this file
    duration	  Integer	Duration of the audio in seconds as defined by sender
    mime_type	  String	Optional. MIME type of the file as defined by sender
    file_size	  Integer	Optional. File size
=#

type Voice
    __file_id::UTF8String
    __duration::Int32
    __mime_type::UTF8String
    __file_size::Int32
end

#=
This object represents a phone contact.

    Field	          Type	       Description
    phone_number	  String       Contact's phone number
    first_name	      String	   Contact's first name
    last_name	      String	   Optional. Contact's last name
    user_id	          Integer	   Optional. Contact's user identifier in Telegram
=#

type Contact
    __phone_number::UTF8String
    __first_name::UTF8String
    __last_name::UTF8String
    __user_id::Int32
end

#=
This object represents a point on the map.

    Field	    Type	Description
    longitude	Float	Longitude as defined by sender
    latitude	Float	Latitude as defined by sender
=#

type Location
    __longitude::Float32
    __latitude::Float32
end

#=
This object represents a venue.

    Field	          Type	    Description
    location	      Location	Venue location
    title	          String	Name of the venue
    address	          String	Address of the venue
    foursquare_id	  String	Optional. Foursquare identifier of the venue
=#

type Venue
    __location::Location
    __title::UTF8String
    __address::UTF8String
    __foursquare_id::UTF8String
end

#=
This object represents a message.

    Field	                  Type	                   Description
    message_id	              Integer	               Unique message identifier
    from	                  User	                   Optional. Sender, can be empty for messages sent to channels
    date	                  Integer	               Date the message was sent in Unix time
    chat	                  Chat	                   Conversation the message belongs to
    forward_from	          User	                   Optional. For forwarded messages, sender of the original message
    forward_date	          Integer	               Optional. For forwarded messages, date the original message was sent in Unix time
    reply_to_message	      Message	               Optional. For replies, the original message. Note that the Message object in this field will not contain further reply_to_message fields even if it itself is a reply.
    text	                  String	               Optional. For text messages, the actual UTF-8 text of the message, 0-4096 characters.
    entities	              Array of MessageEntity   Optional. For text messages, special entities like usernames, URLs, bot commands, etc. that appear in the text
    audio	                  Audio	                   Optional. Message is an audio file, information about the file
    document	              Document	               Optional. Message is a general file, information about the file
    photo	                  Array of PhotoSize	   Optional. Message is a photo, available sizes of the photo
    sticker	                  Sticker	               Optional. Message is a sticker, information about the sticker
    video	                  Video	                   Optional. Message is a video, information about the video
    voice	                  Voice	                   Optional. Message is a voice message, information about the file
    caption	                  String	               Optional. Caption for the document, photo or video, 0-200 characters
    contact	                  Contact	               Optional. Message is a shared contact, information about the contact
    location	              Location	               Optional. Message is a shared location, information about the location
    venue	                  Venue	                   Optional. Message is a venue, information about the venue
    new_chat_member	          User	                   Optional. A new member was added to the group, information about them (this member may be the bot itself)
    left_chat_member	      User	                   Optional. A member was removed from the group, information about them (this member may be the bot itself)
    new_chat_title	          String	               Optional. A chat title was changed to this value
    new_chat_photo	          Array of PhotoSize	   Optional. A chat photo was change to this value
    delete_chat_photo	      True	                   Optional. Service message: the chat photo was deleted
    group_chat_created	      True	                   Optional. Service message: the group has been created
    supergroup_chat_created	  True	                   Optional. Service message: the supergroup has been created
    channel_chat_created	  True	                   Optional. Service message: the channel has been created
    migrate_to_chat_id	      Integer	               Optional. The group has been migrated to a supergroup with the specified identifier, not exceeding 1e13 by absolute value
    migrate_from_chat_id	  Integer	               Optional. The supergroup has been migrated from a group with the specified identifier, not exceeding 1e13 by absolute value
    pinned_message	          Message	               Optional. Specified message was pinned. Note that the Message object in this field will not contain further reply_to_message fields even if it is itself a reply.
=#

type Message
    __message_id::Int32
    __from::User
    __date::Int32
    __chat::Chat
    __forward_from::User
    __forward_date::Int32
    __replay_to_message::Message
    __text::UTF8String
    __entities::Array{ MessageEntity, 1 }
    __audio::Audio
    __document::Document
    __photo::Array{ PhotoSize, 1 }
    __sticker::Sticker
    __video::Video
    __voice::Voice
    __caption::UTF8String
    __contact::Contact
    __location::Location
    __venue::Venue
    __new_chat_member::User
    __left_chat_member::User
    __new_chat_title::UTF8String
    __new_chat_photo::Array{ PhotoSize, 1 }
    __delete_chat_photo::Bool=true
    __group_chat_created::Bool=true
    __channel_chat_created::Bool=true
    __migrate_to_chat_id::Int32
    __pinned_message::Message
end

#=
This object represents an incoming inline query. When the user sends an empty query, your bot could return some default or trending results.

    Field	    Type	     Description
    id	        String	     Unique identifier for this query
    from	    User	     Sender
    location	Location	 Optional. Sender location, only for bots that request user location
    query	    String	     Text of the query
    offset	    String	     Offset of the results to be returned, can be controlled by the bot
=#

type InlineQuery
    __id::UTF8String
    __from::User
    __location::Location
    __query::UTF8String
    __offset::UTF8String
end

#=
Represents a result of an inline query that was chosen by the user and sent to their chat partner.

    Field	               Type	      Description
    result_id	           String	  The unique identifier for the result that was chosen
    from	               User	      The user that chose the result
    location	           Location	  Optional. Sender location, only for bots that require user location
    inline_message_id	   String	  Optional. Identifier of the sent inline message. Available only if there is an inline keyboard attached to the message. Will be also received in callback queries and can be used to edit the message.
    query	               String	  The query that was used to obtain the result
=#

type ChosenInlineResult
    __result_id::UTF8String
    __from::User
    __location::Location
    __inline_message_id::UTF8String
    __query::UTF8String
end

#=
This object represents an incoming callback query from a callback button in an inline keyboard. If the button that originated the query was attached to a message sent by the bot, the field message will be presented. If the button was attached to a message sent via the bot (in inline mode), the field inline_message_id will be presented.

    Field	              Type	       Description
    id	                  String	   Unique identifier for this query
    from	              User	       Sender
    message	              Message	   Optional. Message with the callback button that originated the query. Note that message content and message date will not be available if the message is too old
    inline_message_id	  String	   Optional. Identifier of the message sent via the bot in inline mode, that originated the query
    data	              String	   Data associated with the callback button. Be aware that a bad client can send arbitrary data in this field
=#

type CallbackQuery
    __id::UTF8String
    __from::User
    __message::Message
    __inline_message_id::UTF8String
    __data::UTF8String
end

#=
This object represents an incoming update.
Only one of the optional parameters can be present in any given update.

    Field	              Type	               Description
    update_id	          Integer	           The update‘s unique identifier. Update identifiers start from a certain positive number and increase sequentially. This ID becomes especially handy if you’re using Webhooks, since it allows you to ignore repeated updates or to restore the correct update sequence, should they get out of order.
    message	              Message	           Optional. New incoming message of any kind — text, photo, sticker, etc.
    inline_query	      InlineQuery	       Optional. New incoming inline query
    chosen_inline_result  ChosenInlineResult   Optional. The result of an inline query that was chosen by a user and sent to their chat partner.
    callback_query	      CallbackQuery	       Optional. New incoming callback query
=#

type Update
    __update_id::Int32
    __message::Message
    __inline_query::InlineQuery
    __chosen_inline_result::ChosenInlineResult
    __callback_query::CallbackQuery
end

#=
This object represent a user's profile pictures.

    Field	        Type	                        Description
    total_count	    Integer	                        Total number of profile pictures the target user has
    photos	        Array of Array of PhotoSize	    Requested profile pictures (in up to 4 sizes each)
=#

type UserProfilePhotos
    __total_count::Int32
    __photos::Array{ Array{ PhotoSize, 1 }, 1 }
end

#=
This object represents a file ready to be downloaded. The file can be downloaded via the link https://api.telegram.org/file/bot<token>/<file_path>. It is guaranteed that the link will be valid for at least 1 hour. When the link expires, a new one can be requested by calling getFile.
Maximum file size to download is 20 MB

    Field	      Type	    Description
    file_id	      String	Unique identifier for this file
    file_size	  Integer	Optional. File size, if known
    file_path	  String	Optional. File path. Use https://api.telegram.org/file/bot<token>/<file_path> to get the file.
=#

type File
    __file_id::UTF8String
    __file_size::Int32
    __file_path::UTF8String
end

#=
This object represents one button of the reply keyboard. For simple text buttons String can be used instead of this object to specify text of the button. Optional fields are mutually exclusive.

    Field	              Type	    Description
    text	              String	Text of the button. If none of the optional fields are used, it will be sent to the bot as a message when the button is pressed
    request_contact	      Boolean	Optional. If True, the user's phone number will be sent as a contact when the button is pressed. Available in private chats only
    request_location	  Boolean	Optional. If True, the user's current location will be sent when the button is pressed. Available in private chats only
=#

type KeyboardButton
    __text::UTF8String
    __requested_contact::Bool
    __requested_location::Bool
end

#=
This object represents a custom keyboard with reply options (see Introduction to bots for details and examples).

    Field	              Type	                               Description
    keyboard           	  Array of Array of KeyboardButton	   Array of button rows, each represented by an Array of KeyboardButton objects
    resize_keyboard	      Boolean	                           Optional. Requests clients to resize the keyboard vertically for optimal fit (e.g., make the keyboard smaller if there are just two rows of buttons). Defaults to false, in which case the custom keyboard is always of the same height as the app's standard keyboard.
    one_time_keyboard	  Boolean	                           Optional. Requests clients to hide the keyboard as soon as it's been used. The keyboard will still be available, but clients will automatically display the usual letter-keyboard in the chat – the user can press a special button in the input field to see the custom keyboard again. Defaults to false.
    selective	          Boolean	                           Optional. Use this parameter if you want to show the keyboard to specific users only. Targets: 1) users that are @mentioned in the text of the Message object; 2) if the bot's message is a reply (has reply_to_message_id), sender of the original message.

    Example: A user requests to change the bot‘s language, bot replies to the request with a keyboard to select the new language. Other users in the group don’t see the keyboard.
=#

type ReplyKeyboardMarkup
    __keyboard::Array{ Array{ KeyboardButton, 1 }, 1 }
    __resize_keyboard::Bool
    __one_time_keyboard::Bool
    __selective::Bool
end

#=
Upon receiving a message with this object, Telegram clients will hide the current custom keyboard and display the default letter-keyboard. By default, custom keyboards are displayed until a new keyboard is sent by a bot. An exception is made for one-time keyboards that are hidden immediately after the user presses a button (see ReplyKeyboardMarkup).

    Field	          Type	       Description
    hide_keyboard	  True	       Requests clients to hide the custom keyboard
    selective	      Boolean	   Optional. Use this parameter if you want to hide keyboard for specific users only. Targets: 1) users that are @mentioned in the text of the Message object; 2) if the bot's message is a reply (has reply_to_message_id), sender of the original message.

    Example: A user votes in a poll, bot returns confirmation message in reply to the vote and hides keyboard for that user, while still showing the keyboard with poll options to users who haven't voted yet.
=#

type ReplyKeyboardHide
    __hide_keyboard::Bool=true
    __selective::Bool
end

#=
This object represents one button of an inline keyboard. You must use exactly one of the optional fields.

    Field	              Type	    Description
    text	              String	Label text on the button
    url	                  String	Optional. HTTP url to be opened when button is pressed
    callback_data	      String	Optional. Data to be sent in a callback query to the bot when button is pressed, 1-64 bytes
    switch_inline_query	  String	Optional. If set, pressing the button will prompt the user to select one of their chats, open that chat and insert the bot‘s username and the specified inline query in the input field. Can be empty, in which case just the bot’s username will be inserted.

Note: This offers an easy way for users to start using your bot in inline mode when they are currently in a private chat with it. Especially useful when combined with switch_pm… actions – in this case the user will be automatically returned to the chat they switched from, skipping the chat selection screen.
Note: This will only work in Telegram versions released after 9 April, 2016. Older clients will display unsupported message.
=#

type InlineKeyboardButton
    __text::UTF8String
    __url::UTF8String
    __callback_data::UTF8String
    __switch_inline_query::UTF8String
end

#=
This object represents an inline keyboard that appears right next to the message it belongs to.
Warning: Inline keyboards are currently being tested and are only available in one-on-one chats (i.e., user-bot or user-user in the case of inline bots).

    Field	          Type	                                    Description
    inline_keyboard	  Array of Array of InlineKeyboardButton	Array of button rows, each represented by an Array of InlineKeyboardButton objects

Note: This will only work in Telegram versions released after 9 April, 2016. Older clients will display unsupported message
=#

type InlineKeyboardMarkup
    __inline_keyboard::Array{ Array{ InlineKeyboardButton, 1 }, 1 }
end

#=
Upon receiving a message with this object, Telegram clients will display a reply interface to the user (act as if the user has selected the bot‘s message and tapped ’Reply'). This can be extremely useful if you want to create user-friendly step-by-step interfaces without having to sacrifice privacy mode.

    Field	       Type	      Description
    force_reply	   True	      Shows reply interface to the user, as if they manually selected the bot‘s message and tapped ’Reply'
    selective	   Boolean    Optional. Use this parameter if you want to force reply from specific users only. Targets: 1) users that are @mentioned in the text of the Message object; 2) if the bot's message is a reply (has reply_to_message_id), sender of the original message.

Example: A poll bot for groups runs in privacy mode (only receives commands, replies to its messages and mentions). There could be two ways to create a new poll:
Explain the user how to send a command with parameters (e.g. /newpoll question answer1 answer2). May be appealing for hardcore users but lacks modern day polish.
Guide the user through a step-by-step process. ‘Please send me your question’, ‘Cool, now let’s add the first answer option‘, ’Great. Keep adding answer options, then send /done when you‘re ready’.
The last option is definitely more attractive. And if you use ForceReply in your bot‘s questions, it will receive the user’s answers even if it only receives replies, commands and mentions — without any extra work for the user.
=#

type ForceReply
    __force_reply::Bool=true
    __selective::Bool
end

max_stars_repo_name: Fazendaaa/Julia-Telegram-Bot

max_issues_repo_name: Fazendaaa/Julia-Telegram-Bot

max_forks_repo_name: Fazendaaa/Julia-Telegram-Bot


text:
abstract AbstractWeighting

immutable EqualWeighting <: AbstractWeighting end
weight_frequencies(::EqualWeighting, ξs) = ones(ξs)

immutable LoudnessWeighting <: AbstractWeighting
    samplerate::Int
end

function weight_frequencies(weighting::LoudnessWeighting, ξs)
    freqs = weighting.samplerate * ξs
    freqs2 = freqs .* freqs
    numerator = 12200.0 * freqs2
    sqden1 = (freqs2 + 20.6 * 20.6)
    sqden2 = sqrt.((freqs2 + 107.7 * 107.7) .* (freqs2 + 737.9 * 737.9))
    sqden3 = (freqs2 + 12200.0 * 12200.0)
    sqdenominator = sqden1 .* sqden2 .* sqden3
    denominator = sqrt.(sqdenominator)
    return numerator ./ denominator
end


text:
using Roofline

@noinline function experiment(A)
    acc = zero(eltype(A))
    @simd for i in eachindex(A)
        @inbounds acc += A[i]
    end
    acc
end

function setup(N)
    data = rand(Float32, N)
    return (data,)
end

bench = Roofline.benchmark(experiment, setup, 2^27)


text:
"""
VAR1SimPs(A,epsilon,T,x0=0.0)

Calculate impulse response function of a VAR(1) system
x(t) = A * x(t-1) +  epsilon(t), where x(t) is nx1

# Input
- `A::Matrix`:                  nxn VAR(1) matrix, see above
- `epsilon::Number or Vector`:  n-vector of shocks in inital period, or Txn matrix with shocks in all periods
- `T::Number`:                  scalar, last period to calculate for
- `x0::Number or Vector`:       n-vector with starting values, optional

# Output
- `xM::Matrix`:               Txn matrix, impulse response function

Paul.Soderlind@unisg.ch, to Julia Nov 2015

"""
function VAR1SimPs(A,epsilon,T,x0=0.0)

  n = size(A,1)

  isa(x0,Number)         && (x0      = fill(x0,n))           #if scalar
  isa(epsilon,Number)    && (epsilon = fill(epsilon,n))
  (length(epsilon) == n) && (epsilon = vcat(vec(epsilon)',zeros(T-1,n)))

  xM      = fill(NaN,(T,n))                        #to put results in
  xM[1,:] = A*vec(x0) + epsilon[1,:]
  for t = 2:T                                      #loop over time periods
    xM[t,:] = A*xM[t-1,:] + epsilon[t,:]
  end

  return xM

end
#-----------------------------------------------------------------------

max_stars_repo_name: PaulSoderlind/PaulSoderlindCode

max_issues_repo_name: PaulSoderlind/PaulSoderlindCode

max_forks_repo_name: PaulSoderlind/PaulSoderlindCode


text:
using DataStructures

function getNbors(M::Matrix{Int}, pos::CartesianIndex{2}, expand::Int) :: Vector{CartesianIndex{2}}
  s = size(M) .* expand
  nbors :: Vector{CartesianIndex{2}} = []
  (pos[1] != 1   ) && push!(nbors, pos + CartesianIndex(-1,0))
  (pos[1] != s[1]) && push!(nbors, pos + CartesianIndex(1,0))
  (pos[2] != 1   ) && push!(nbors, pos + CartesianIndex(0,-1))
  (pos[2] != s[2]) && push!(nbors, pos + CartesianIndex(0,1))
  return nbors
end

getNbors(M::Matrix{Int}, pos::CartesianIndex{2}) = getNbors(M,pos,1)

function getWeights(M::Matrix{Int}, pos::CartesianIndex{2}) :: Int
  s = size(M)
  return ((M[(pos[1]-1)%s[1]+1, (pos[2]-1)%s[2]+1]
          + ((pos[1]-1)÷s[1]) + ((pos[2]-1)÷s[2])) - 1) % 9 +1
end

# Checks if neighbor belongs in the queue:
#  - If neighbor has not been visited, adds it to the queue
#  - If neighbor has been visited but with bigger weight, checks if it is still in the queue
#    and either adds it (if not present) or updates weight (if present).
function addNbor!(q::PriorityQueue{CartesianIndex{2},Int}, n::CartesianIndex{2}, w::Int, W::Matrix{Int})
  (W[n] > w) && begin (W[n] = w) ; (in(n,keys(q)) ? q[n] = w : enqueue!(q,n,w)) end
end

# Dijkstra-based algorithm to find the path with less risk.
function bestPath(M::Matrix{Int}, pstart::CartesianIndex, pend::CartesianIndex, expand::Int) :: Int
  W :: Matrix{Int} = 5*sum(M)*ones(size(M) .* expand) # Weight matrix
  q = PriorityQueue{CartesianIndex{2}, Int}()
  enqueue!(q, pstart, 0)
  W[pstart] = 0

  while !isempty(q)
    # Retrieve next position
    pos, weight = dequeue_pair!(q)

    # Have we arrived?
    pos == pend && return weight

    # Otherwise, enqueue new path options
    nbors = getNbors(M,pos,expand)
    map(n -> addNbor!(q,n,weight+getWeights(M,n),W), nbors)
  end
  return -1
end

function readInput() :: Matrix{Int}
  f = open("./inputs/day15.txt")
  data :: Matrix{Int} = zeros(100,100)

  i = 1
  for line in eachline(f)
    data[i,:] = map(x -> parse(Int,x), collect(line))
    i += 1
  end
  return data
end

M = readInput()

expand = 5
res = bestPath(M, CartesianIndex(1,1), CartesianIndex(100*expand,100*expand), expand)
println(res)


text:
using BenchmarkTools

using LinearAlgebra

using Krylov
using LinearOperators
using MatrixMarket

include("../test/get_div_grad.jl")
include("../test/test_utils.jl")
include("fetch_matrices.jl")

# we don't want to use matrixdepot(matrix, :read) to read in matrices because of
# https://github.com/JuliaMatrices/MatrixDepot.jl/issues/26

const SUITE = BenchmarkGroup()

SUITE["CG"] = BenchmarkGroup(["CG", "SPD"])

for N in [32, 64, 128]
  A = get_div_grad(N, N, N)
  n = size(A, 1)
  b = ones(n)
  op = preallocated_LinearOperator(A)
  M = nonallocating_opEye(n)
  rtol = 1.0e-6
  SUITE["CG"]["DivGrad N=$N"] = @benchmarkable cg($op, $b, M=$M, atol=0.0, rtol=$rtol, itmax=$n)
end

SUITE["CG"]["UFL-small"] = BenchmarkGroup()
for matrix in spd_small
  name = basename(matrix)
  A = MatrixMarket.mmread(joinpath(matrix_path, "..", "data", "uf", matrix, "$(name).mtx"))
  n = size(A, 1)
  b = ones(n)
  op = preallocated_LinearOperator(A)
  M = nonallocating_opEye(n)
  rtol = 1.0e-6
  SUITE["CG"]["UFL-small"][matrix] = @benchmarkable cg($op, $b, M=$M, atol=0.0, rtol=$rtol, itmax=$n)
end

SUITE["CG"]["UFL-medium"] = BenchmarkGroup()
for matrix in spd_med
  name = basename(matrix)
  A = MatrixMarket.mmread(joinpath(matrix_path, "..", "data", "uf", matrix, "$(name).mtx"))
  n = size(A, 1)
  b = ones(n)
  op = preallocated_LinearOperator(A)
  M = nonallocating_opEye(n)
  rtol = 1.0e-6
  SUITE["CG"]["UFL-medium"][matrix] = @benchmarkable cg($op, $b, M=$M, atol=0.0, rtol=$rtol, itmax=$n)
end

SUITE["CG"]["UFL-large"] = BenchmarkGroup()
for matrix in spd_large
  name = basename(matrix)
  A = MatrixMarket.mmread(joinpath(matrix_path, "..", "data", "uf", matrix, "$(name).mtx"))
  n = size(A, 1)
  b = ones(n)
  op = preallocated_LinearOperator(A)
  M = nonallocating_opEye(n)
  rtol = 1.0e-6
  SUITE["CG"]["UFL-large"][matrix] = @benchmarkable cg($op, $b, M=$M, atol=0.0, rtol=$rtol, itmax=$n)
end


text:
mview(c, inds, ::Val{1}) = view(c, :, inds[1], :)
mview(c, inds, ::Val{2}) = view(c, :, inds[1], inds[2], :)
mview(c, inds, ::Val{3}) = view(c, :, inds[1], inds[2], inds[3], :)
mview(c, inds, ::Val{4}) = view(c, :, inds[1], inds[2], inds[3], inds[4], :)

tensor_contraction(
    A,
    B,
    ::Val{1},
) = @tullio C[o, a, b] := A[i, a, o] * B[i, a, b]
tensor_contraction(A, B, ::Val{2}) =
    @tullio C[o, a₁, a₂, b] := A[i, a₁, a₂, o] * B[i, a₁, a₂, b]
tensor_contraction(A, B, ::Val{3}) =
    @tullio C[o, a₁, a₂, a₃, b] := A[i, a₁, a₂, a₃, o] * B[i, a₁, a₂, a₃, b]
tensor_contraction(A, B, ::Val{4}) = @tullio C[o, a₁, a₂, a₃, a₄, b] :=
    A[i, a₁, a₂, a₃, a₄, o] * B[i, a₁, a₂, a₃, a₄, b]

sparse_mean(w, c, ::Val{1}) = @tullio μ[o, a, b] := w[i, o] * c[i, a, b]
sparse_mean(
    w,
    c,
    ::Val{2},
) = @tullio μ[o, a₁, a₂, b] := w[i, o] * c[i, a₁, a₂, b]
sparse_mean(w, c, ::Val{3}) =
    @tullio μ[o, a₁, a₂, a₃, b] := w[i, o] * c[i, a₁, a₂, a₃, b]
sparse_mean(w, c, ::Val{4}) =
    @tullio μ[o, a₁, a₂, a₃, a₄, b] := w[i, o] * c[i, a₁, a₂, a₃, a₄, b]

sparse_covariance(
    w,
    c,
    ::Val{1},
) = @tullio μ[o, a, r, b] := w[i, r, o] * c[i, a, b]
sparse_covariance(w, c, ::Val{2}) =
    @tullio μ[o, a₁, a₂, r, b] := w[i, r, o] * c[i, a₁, a₂, b]
sparse_covariance(w, c, ::Val{3}) =
    @tullio μ[o, a₁, a₂, a₃, r, b] := w[i, r, o] * c[i, a₁, a₂, a₃, b]
sparse_covariance(w, c, ::Val{4}) =
    @tullio μ[o, a₁, a₂, a₃, a₄, r, b] := w[i, r, o] * c[i, a₁, a₂, a₃, a₄, b]


text:
# Copyright 2021 The CLVR Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# This file contains gerenal helper functions.

"""Read results CSV file into a Results object."""
function csv_to_results(filepath)

    csv_reader = CSV.File(filepath)

    results = Results()
    for row in csv_reader
        logresult!(results, row.iterations, row.times, row.fvaluegaps, row.metricLPs)
    end
    return results
end

max_stars_repo_name: ericlincc/Efficient-GLP

max_issues_repo_name: ericlincc/Efficient-GLP

max_forks_repo_name: ericlincc/Efficient-GLP


text:
using ValueHistories
using Base.Test

tests = [
    "tst_history.jl"
    "tst_mvhistory.jl"
]

perf = [
    "bm_history.jl"
]

for t in tests
    @testset "[->] $t" begin
        include(t)
    end
end

for p in perf
    @testset "[->] $p" begin
        include(p)
    end
end


text:
module AXMethods

# Dependencies
using DataFrames
using Distributions
using LinearAlgebra
using Random
using Statistics

# Module Export
export AXEstimator, AXLeastSquares, AX2StageLeastSquares
export coef, predict, infer

# Module Content
include("G:\\My Drive\\PhD\\AXMethods\\src\\FUN\\AXEstimator.jl")
include("G:\\My Drive\\PhD\\AXMethods\\src\\FUN\\AXLeastSquares.jl")
include("G:\\My Drive\\PhD\\AXMethods\\src\\FUN\\AX2StageLeastSquares.jl")

end


text:
@testset "Simulation: Particle (3D)" begin
    s = ContactImplicitMPC.get_simulation("particle", "flat_3D_lc", "flat_lc")
    model = s.model

    # time
    h = 0.01
    T = 100

    ## DROP
    # initial conditions
    q0 = @SVector [0.0, 0.0, 1.0]
    q1 = @SVector [0.0, 0.0, 1.0]

    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8, κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = false))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test all(isapprox.(sim.traj.q[end], 0.0, atol = 1.0e-6))

    ## DROP (warmstart)
    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8,
            κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(
            warmstart = true,
            z_warmstart = 0.001,
            κ_warmstart = 0.001))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test all(isapprox.(sim.traj.q[end], 0.0, atol = 1.0e-6))

    ## SLIDE
    # initial conditions
    v1 = @SVector [1.0, 2.0, 0.0]
    q1 = @SVector [0.0, 0.0, 1.0]
    q0 = q1 - h * v1

    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8, κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = false))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test isapprox.(sim.traj.q[end][3], 0.0, atol = 1.0e-6)
    @test all(isapprox.((sim.traj.q[end] - sim.traj.q[end-1]) ./ h, 0.0, atol = 1.0e-6))
    # @show (sim.traj.q[end] - sim.traj.q[end-1]) ./ h

    ## DROP (warmstart)
    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8,
            κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(
            warmstart = true,
            z_warmstart = 0.001,
            κ_warmstart = 0.001))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test isapprox.(sim.traj.q[end][3], 0.0, atol = 1.0e-6)
    @test all(isapprox.((sim.traj.q[end] - sim.traj.q[end-1]) ./ h, 0.0, atol = 1.0e-6))
end

@testset "Simulation: Particle (3D), Code-Gen Solver" begin
    s = ContactImplicitMPC.get_simulation("particle", "flat_3D_lc", "flat_lc")
    model = s.model

    # time
    h = 0.01
    T = 100

    ## DROP
    # initial conditions
    q0 = @SVector [0.0, 0.0, 1.0]
    q1 = @SVector [0.0, 0.0, 1.0]

    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8, κ_tol = 1.0e-8, solver = :lu_solver),
        sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = false))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test all(isapprox.(sim.traj.q[end], 0.0, atol = 1.0e-6))

    ## SLIDE
    # initial conditions
    v1 = @SVector [1.0, 2.0, 0.0]
    q1 = @SVector [0.0, 0.0, 1.0]
    q0 = q1 - h * v1

    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8, κ_tol = 1.0e-8, solver = :lu_solver),
        sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = false))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test isapprox.(sim.traj.q[end][3], 0.0, atol = 1.0e-6)
    @test all(isapprox.((sim.traj.q[end] - sim.traj.q[end-1]) ./ h, 0.0, atol = 1.0e-6))
end

@testset "Simulation: Particle (3D), Approximate Residual Jacobian" begin

    s = ContactImplicitMPC.get_simulation("particle", "flat_3D_lc", "flat_lc")
    model = s.model

    dir_dyn = joinpath(module_dir(), "src/dynamics/particle")
    dir_sim = joinpath(module_dir(), "src/simulation/particle")
    model2 = deepcopy(ContactImplicitMPC.particle)
    env2 = environment_3D_flat()
    s2 = Simulation(model2, env2)

    path_base = joinpath(dir_dyn, "dynamics/base.jld2")
    path_dyn = joinpath(dir_dyn, "dynamics/dynamics.jld2")
    path_res = joinpath(dir_sim, "quadratic_approx/residual.jld2")
    path_jac = joinpath(dir_sim, "quadratic_approx/sparse_jacobians.jld2")

    instantiate_base!(s2.model, path_base)

    expr_dyn = generate_dynamics_expressions(s2.model, derivs = true)
    save_expressions(expr_dyn, path_dyn, overwrite=true)
    instantiate_dynamics!(s2.model, path_dyn, derivs = true)

    expr_res, rz_sp, rθ_sp = generate_residual_expressions(s2.model, s2.env, jacobians = :approx)
    save_expressions(expr_res, path_res, overwrite=true)
    @save path_jac rz_sp rθ_sp
    instantiate_residual!(s2, path_res, path_jac, jacobians = :approx)

    r0 = zeros(num_var(s.model, s.env))
    s.res.r!(r0, 5.0 * ones(num_var(s.model, s.env)), 0.1 * ones(num_data(model)), 1.0)
    s.res.rz!(s.rz, 5.0 * ones(num_var(s.model, s.env)), 0.1 * ones(num_data(model)))
    s.res.rθ!(s.rθ, ones(num_var(s.model, s.env)), ones(num_data(model)))

    r1 = zeros(num_var(s2.model, s2.env))
    s2.res.r!(r1, 5.0 * ones(num_var(s2.model, s2.env)), 0.1 * ones(num_data(s2.model)), 1.0)
    s2.res.rz!(s2.rz, 5.0 * ones(num_var(s2.model, s2.env)), 0.1 * ones(num_data(s2.model)))
    s2.res.rθ!(s2.rθ, ones(num_var(s2.model, s2.env)), ones(num_data(s2.model)))

    @test norm(r0 - r1) < 1.0e-8
    @test norm(s.rz - s2.rz) < 1.0e-8
    @test norm(s.rθ - s2.rθ) < 1.0e-8

    # time
    h = 0.01
    T = 100

    ## DROP
    # initial conditions
    q0 = @SVector [0.0, 0.0, 1.0]
    q1 = @SVector [0.0, 0.0, 1.0]

    # simulator
    sim = ContactImplicitMPC.simulator(s2, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8, κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = false))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test all(isapprox.(sim.traj.q[end], 0.0, atol = 1.0e-6))
end

@testset "Simulation: Particle (3D) in Quadratic Bowl" begin
    s = get_simulation("particle", "quadratic_bowl_3D_lc", "quadratic")
    s.model.μ_world = 0.1

    # time
    h = 0.01
    T = 1000

    # initial conditions
    q1 = @SVector [1.0, 0.5, 2.0]
    q0 = @SVector [1.1, 0.5, 2.0]

    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
    	ip_opts = ContactImplicitMPC.InteriorPointOptions(
    		r_tol = 1.0e-6, κ_tol = 1.0e-6,
    		diff_sol = true),
    	sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = true))

    # simulate
    @time status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test abs(sim.traj.q[end][1]) < 0.05
    @test abs(sim.traj.q[end][2]) < 0.05
    @test abs(sim.traj.q[end][3]) < 0.001

    ## Drop
    # initial conditions
    q1 = @SVector [1.0, 0.5, 2.0]
    q0 = @SVector [1.0, 0.5, 2.0]

    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8, κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = false))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test abs(sim.traj.q[end][1]) < 0.05
    @test abs(sim.traj.q[end][2]) < 0.05
    @test abs(sim.traj.q[end][3]) < 0.001

    ## DROP (warmstart)
    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8,
            κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(
            warmstart = true,
            z_warmstart = 0.001,
            κ_warmstart = 0.001))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test abs(sim.traj.q[end][1]) < 0.05
    @test abs(sim.traj.q[end][2]) < 0.05
    @test abs(sim.traj.q[end][3]) < 0.001
end

@testset "Simulation: Particle (3D) Quadratic Bowl, Approximate Residual Jacobian" begin
    dir_dyn = joinpath(module_dir(), "src/dynamics/particle")
    dir_sim = joinpath(module_dir(), "src/simulation/particle")
    model2 = deepcopy(ContactImplicitMPC.particle)
    model2.μ_world = 0.1
    env2 = deepcopy(ContactImplicitMPC.quadratic_bowl_3D_lc)
    s2 = Simulation(model2, env2)

    path_base = joinpath(dir_dyn, "dynamics/base.jld2")
    path_dyn = joinpath(dir_dyn, "dynamics/dynamics.jld2")
    path_res = joinpath(dir_sim, "quadratic_approx/residual.jld2")
    path_jac = joinpath(dir_sim, "quadratic_approx/sparse_jacobians.jld2")

    instantiate_base!(s2.model, path_base)

    expr_dyn = generate_dynamics_expressions(s2.model, derivs = true)
    save_expressions(expr_dyn, path_dyn, overwrite=true)
    instantiate_dynamics!(s2.model, path_dyn, derivs = true)

    expr_res, rz_sp, rθ_sp = generate_residual_expressions(s2.model, s2.env, jacobians = :approx)
    save_expressions(expr_res, path_res, overwrite=true)
    @save path_jac rz_sp rθ_sp
    instantiate_residual!(s2, path_res, path_jac, jacobians = :approx)

    # time
    h = 0.01
    T = 1000

    ## DROP
    # initial conditions
    q1 = @SVector [1.0, 0.5, 2.0]
    q0 = @SVector [1.1, 0.5, 2.0]

    # simulator
    sim = ContactImplicitMPC.simulator(s2, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8, κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = false))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test abs(sim.traj.q[end][1]) < 0.05
    @test abs(sim.traj.q[end][2]) < 0.05
    @test abs(sim.traj.q[end][3]) < 0.001
end

@testset "Simulation: Particle (2D)" begin
    s = get_simulation("particle_2D", "flat_2D_lc", "flat_lc")

    # time
    h = 0.01
    T = 100

    ## DROP
    # initial conditions
    q0 = @SVector [0.0, 1.0]
    q1 = @SVector [0.0, 1.0]

    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8, κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = false))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test all(isapprox.(sim.traj.q[end], 0.0, atol = 1.0e-6))

    ## DROP (warmstart)
    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8,
            κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(
            warmstart = true,
            z_warmstart = 0.001,
            κ_warmstart = 0.001))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test all(isapprox.(sim.traj.q[end], 0.0, atol = 1.0e-6))

    ## SLIDE
    # initial conditions
    v1 = @SVector [1.0, 0.0]
    q1 = @SVector [0.0, 1.0]
    q0 = q1 - h * v1

    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8, κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = false))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test isapprox.(sim.traj.q[end][2], 0.0, atol = 1.0e-6)
    @test all(isapprox.((sim.traj.q[end] - sim.traj.q[end-1]) ./ h, 0.0, atol = 1.0e-6))

    ## DROP (warmstart)
    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8,
            κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(
            warmstart = true,
            z_warmstart = 0.001,
            κ_warmstart = 0.001))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test isapprox.(sim.traj.q[end][2], 0.0, atol = 1.0e-6)
    @test all(isapprox.((sim.traj.q[end] - sim.traj.q[end-1]) ./ h, 0.0, atol = 1.0e-6))
end

@testset "Simulation: Particle (2D + nonlinear cone)" begin
    s = get_simulation("particle_2D", "flat_2D_nc", "flat_nc")

    # time
    h = 0.01
    T = 100

    ## DROP
    # initial conditions
    q0 = @SVector [0.0, 1.0]
    q1 = @SVector [0.0, 1.0]

    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8, κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = false))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test all(isapprox.(sim.traj.q[end], 0.0, atol = 1.0e-6))

    ## DROP (warmstart)
    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8,
            κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(
            warmstart = true,
            z_warmstart = 0.001,
            κ_warmstart = 0.001))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test all(isapprox.(sim.traj.q[end], 0.0, atol = 1.0e-6))

    ## SLIDE
    # initial conditions
    v1 = @SVector [1.0, 0.0]
    q1 = @SVector [0.0, 1.0]
    q0 = q1 - h * v1

    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8, κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = false))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test isapprox.(sim.traj.q[end][2], 0.0, atol = 1.0e-6)
    @test all(isapprox.((sim.traj.q[end] - sim.traj.q[end-1]) ./ h, 0.0, atol = 1.0e-6))

    ## DROP (warmstart)
    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8,
            κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(
            warmstart = true,
            z_warmstart = 0.001,
            κ_warmstart = 0.001))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test isapprox.(sim.traj.q[end][2], 0.0, atol = 1.0e-6)
    @test all(isapprox.((sim.traj.q[end] - sim.traj.q[end-1]) ./ h, 0.0, atol = 1.0e-6))
end

@testset "Simulation: Particle (2D) On Slope" begin
    s = get_simulation("particle_2D", "slope1_2D_lc", "slope")
    s.model.μ_world = 0.1

    # time
    h = 0.01
    T = 100

    ## DROP
    # initial conditions
    q0 = @SVector [0.0, 1.0]
    q1 = @SVector [0.0, 1.0]

    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8, κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = false))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test sim.traj.q[end][1] < 0.0
    @test sim.traj.q[end][2] < 0.0
end

@testset "Simulation: Particle (3D nonlinear cone)" begin
    s = ContactImplicitMPC.get_simulation("particle", "flat_3D_nc", "flat_nc")
    model = s.model

    # time
    h = 0.01
    T = 100

    ## DROP
    # initial conditions
    q0 = @SVector [0.0, 0.0, 1.0]
    q1 = @SVector [0.0, 0.0, 1.0]

    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8, κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = false))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test all(isapprox.(sim.traj.q[end], 0.0, atol = 1.0e-6))

    ## DROP (warmstart)
    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8,
            κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(
            warmstart = true,
            z_warmstart = 0.001,
            κ_warmstart = 0.001))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test all(isapprox.(sim.traj.q[end], 0.0, atol = 1.0e-6))

    ## SLIDE
    # initial conditions
    v1 = @SVector [1.0, 2.0, 0.0]
    q1 = @SVector [0.0, 0.0, 1.0]
    q0 = q1 - h * v1

    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8, κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = false))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test isapprox.(sim.traj.q[end][3], 0.0, atol = 1.0e-6)
    @test all(isapprox.((sim.traj.q[end] - sim.traj.q[end-1]) ./ h, 0.0, atol = 1.0e-6))

    ## DROP (warmstart)
    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
            r_tol = 1.0e-8,
            κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(
            warmstart = true,
            z_warmstart = 0.001,
            κ_warmstart = 0.001))

    # simulate
    status = ContactImplicitMPC.simulate!(sim)
    @test status
    @test isapprox.(sim.traj.q[end][3], 0.0, atol = 1.0e-6)
    @test all(isapprox.((sim.traj.q[end] - sim.traj.q[end-1]) ./ h, 0.0, atol = 1.0e-6))
end


text:
export @subham
using MacroTools
using MacroTools: postwalk, @match

merge_placeholder_expr(::Tuple{}, holders::NTuple{N, Symbol}) where N = holders

function merge_placeholder_expr(a::NTuple{N1, Symbol}, b::NTuple{N2, Symbol}) where {N1, N2}
    (a..., Tuple(each for each in b if !(each in a))...)
end

function is_hamiltonian_expr(expr)
    @match expr begin
        # Z[i]
        # A[i, j, k]
        _[__] => true

        # Z[i]Z[j]
        # A[i, j, k]B[i, j]
        *(args__) => all(map(is_hamiltonian_expr, args))
        ⊗(args__) => all(map(is_hamiltonian_expr, args))

        # Expr(:call, [])
        _ => false
    end
end

function decompose(expr)
    is_hamiltonian_expr(expr) || throw(Meta.ParseError("Invalid Syntax, expect hamiltonian expr, got $expr"))

    ref_list = Expr[]
    postwalk(expr) do x
        Meta.isexpr(x, :ref) && push!(ref_list, x)
        x
    end

    ref_list
end

abstract type HamiltonianExpr end

struct SubHamExpr{PT <: Tuple, VT <: Union{Tuple, Symbol}} <: HamiltonianExpr
    placeholders::PT
    val::VT
end

function SubHamExpr(ex::Expr) # A[i]
    @capture(ex, name_[vars__]) && return SubHamExpr(Tuple(vars), name)

    vals = Tuple(SubHamExpr(each) for each in decompose(ex))

    placeholders = ()
    for each in vals
        placeholders = merge_placeholder_expr(placeholders, each.placeholders)
    end
    SubHamExpr(placeholders, vals)
end

macro subham_lower(ex)
    SubHamExpr(ex)
end

function Base.show(io::IO, ex::SubHamExpr{<:Tuple, Symbol})
    print(io, ex.val, "[", join(ex.placeholders, ", "), "]")
end

function Base.show(io::IO, ex::SubHamExpr{<:Tuple, <:Tuple})
    for each in ex.val
        show(io, each)
    end
end

struct SubHam{PT <: Tuple, VT}
    placeholders::PT
    val::VT
    name::Symbol

    SubHam(ph::PT, val::VT) where {PT, VT} = new{PT, VT}(ph, val)
    SubHam(ph::PT, val::VT, name::Symbol) where {PT, VT} = new{PT, VT}(ph, val, name)
end

Base.eltype(x::SubHam{<:Tuple, <:LocalOperator}) = eltype(x.val)
Base.eltype(x::SubHam{<:Tuple, <:Tuple}) = promote_type((eltype(each) for each in x.val)...)

function make_subham(ex::SubHamExpr{<:Tuple, Symbol})
    placeholders = Expr(:tuple, (esc(each) for each in ex.placeholders)...)
    quote
        SubHam($(placeholders), $(esc(ex.val)), $(QuoteNode(ex.val)))
    end
end

function make_subham(ex::SubHamExpr{<:Tuple, <:Tuple})
    placeholders = Expr(:tuple, (esc(each) for each in ex.placeholders)...)
    vals = Expr(:tuple, (make_subham(each) for each in ex.val)...)
    quote
        SubHam($(placeholders), $(vals))
    end
end

macro subham(expr)
    ex = SubHamExpr(expr)
    make_subham(ex)
end

function _plain(h::SubHam{<:Tuple, <:LocalOperator})
    index_str = join(h.placeholders, ", ")
    "$(h.name)[$index_str]"
end

function _plain(h::SubHam{<:Tuple, <:Tuple})
    join(_plain(each) for each in h.val)
end

Base.show(io::IO, h::SubHam) = print(io, _plain(h))

import LaTeXStrings: latexstring

function _latexstring(h::SubHam{<:Tuple, <:LocalOperator})
    index_str = join(h.placeholders, ", ")
    "$(h.name)_{$index_str}"
end

function _latexstring(h::SubHam{<:Tuple, <:Tuple})
    join(_latexstring(each) for each in h.val)
end

latexstring(h::SubHam) = latexstring(_latexstring(h))

Base.show(io::IO, ::MIME"text/latex", h::SubHam) = print(io, latexstring(h))

nplaceholders(h::SubHam) = length(h.placeholders)
placeholder_type(h::SubHam{PT}) where PT = PT

function merge_placeholder(a::Tuple, b::Tuple)
    (a..., Tuple(each for each in b if !(each in a))...)
end

import Base: *, ==

# NOTE: we only check if placeholders are the same type
==(lhs::SubHam, rhs::SubHam) = false

function ==(lhs::SubHam{PT, <:Tuple}, rhs::SubHam{PT, <:Tuple}) where {PT <: Tuple}
    all(map(==, lhs.val, rhs.val))
end

function ==(lhs::SubHam{PT, <:LocalOperator}, rhs::SubHam{PT, <:LocalOperator}) where {PT <: Tuple}
    lhs.val == rhs.val
end

function *(lhs::SubHam, rhs::SubHam)
    placeholders = merge_placeholder(lhs.placeholders, rhs.placeholders)
    SubHam(placeholders, (lhs, rhs))
end

# error
function *(lhs::SubHam{NTuple{2, Neighbor{K}}}, rhs::SubHam) where K
    error("Neighbor correlation can only be applied on 2-local Hamiltonian")
end

# indexing
# single site
Base.getindex(h::SubHam{<:Tuple, <:LocalOperator}, inds...) = getindex(h.val, inds...)

struct IndexTable{N, PT <: Tuple}
    syms::PT
    lhs_inds::NTuple{N, Int}
    rhs_inds::NTuple{N, Int}
end

IndexTable(h::SubHam, lhs, rhs) = IndexTable(h.placeholders, lhs, rhs)

function pick(table::IndexTable, inds::NTuple{K}) where K
    indices = findall(x->in(x, inds), table.syms)
    table.lhs_inds[indices], table.rhs_inds[indices]
end

# TODO: unroll the loop with generated & Cassette
# NOTE: this should directly generate the following, since we have almost everything
# in type.
# A[(i1, i2, ..., iN), (j1, j2, ..., jN)] := Z[i1, j1] * Z[i2, j2] * C[i3, i4, j3, j4] ...
function Base.getindex(h::SubHam{<:Tuple, <:Tuple}, lhs::NTuple{N, Int}, rhs::NTuple{N, Int}) where N
    table = IndexTable(h, lhs, rhs)
    elem = one(eltype(h))
    for each in h.val
        each_lhs, each_rhs = pick(table, each.placeholders)
        elem *= getindex(each, each_lhs, each_rhs)
    end
    elem
end

max_stars_repo_name: Roger-luo/QuHamiltonian

max_issues_repo_name: Roger-luo/QuHamiltonian

max_forks_repo_name: Roger-luo/QuHamiltonian


text:
module PeeEmTee

using HDF5
using DocStringExtensions
using StatsBase
using Distributions
using Optim
using RecipesBase
using DelimitedFiles
import Base
export high_voltages, WaveSet, calculate_charges, ChargeDist, pmtresp,
       pre_fit, gauss, norm_gauss, pmtresp_fit, pmtresp_uap, plot, plot!,
       calculate_transit_times, bin_data, subtract_baseline, simulate_charges,
       baseline_max, PMTRespFit, PMTRespUapFit, make_qfunc, PreFit, mean_waveform

include("tools.jl")
include("fit.jl")
include("plot.jl")

end # module


text:
## Main loop to extract in e2e the DBSCAN parameters.
## The list of votable is selected directly from the directory

using DataFrames, Query
using CSV, Glob, Dates
using Statistics, Random
using Printf

rootdir =  ENV["GAIA_ROOT"]

push!(LOAD_PATH,"$rootdir/run/src")
using GaiaClustering

## directory
wdir    = "$rootdir/e2e_products"
votdir  = "$wdir/votable"
plotdir = "$wdir/plotSelect"
ocdir   = "$wdir/oc"

## load a liist of votable and update the file if done
## add results
##

function readlist_votable(filelist::String)
    vot = DataFrames.copy(CSV.read(filelist))
    return(vot)
end

function getdata(filevot)
    voname = filevot

    data       = read_votable(voname)
    df         = filter_data(data, [0,2100])
    dfcart     = add_cartesian(df)
    blck       = [[1,2,3],[4,5], [6,7,8]]
    wghtblck   = [4.0,5.0,1.0]
    norm       = "identity"

    dfcartnorm , scale8 = normalization_PerBlock(dfcart, blck, wghtblck , norm, false)
    return(df, dfcart , dfcartnorm)
end

## to check if done and record
## check and updt if votname analyzed. If not done return false
function updt_votcompleted(fileres, votname , cycletot=1, flag= 0 , onlycheck=true)
    let
        if  onlycheck
            if !isfile(fileres)
                return(0, false)
            else
                res = DataFrames.copy(CSV.read(fileres, delim=";"))
                if votname in res.votname
                    x = @from i in res begin
                        @where i.votname == votname
                        @select i
                        @collect DataFrame
                    end
                    return(x, true)
                else
                    return(0, false)
                end
            end
            ## UPDTE
        else
            if !isfile(fileres)
                res = DataFrame(votname=votname, cycle=cycletot, flag=flag)
                CSV.write(fileres,res,delim=';')
                println("## $fileres created...")
                return(res,true)
            else
                res = DataFrames.copy(CSV.read(fileres, delim=";"))
                newrow = DataFrame(votname=votname,cycle=cycletot, flag=flag)
                append!(res,newrow)
                CSV.write(fileres,res,delim=';')
                return(res,true)
            end
        end
    end
end
## check for blacklist
##
function read_blacklist(blackname)
    if isfile(blackname)
        df= DataFrames.copy(CSV.read(blackname, delim=";"))
        blacklist= df.votname
    else
        blacklist= [""]
    end

    return(blacklist)
end

#######################
## Main loop
##
function main(filelist, metafile)
    let
        wd= pwd() ; nfile= size(filelist)[1] ; totalTime= 0.
        println("## Starting main loop using optimization with cycles...")
        println("## It can be very long but will be resumed to the last reduced file.")
        println("## Working directory: $wd")
        println("## $nfile files to analyze...")

        m= read_params(metafile, false)
        m.plotdir= plotdir
        m.ocdir= ocdir
        fileres= "$(m.prefile).done.csv"

        # read a possible votname blacklist
        blackname= "blacklist.csv"
        blacklist= read_blacklist(blackname)
        println("## blacklist read...")

        for i in 1:nfile
            votname = filelist[i]
            res, votfound= updt_votcompleted(fileres, votname, 0, 0, true)

            ## test blacklist
            if votname in blacklist
                println("##")
                println("## $votname in blacklist...")
                println("##")
                votfound= true
            end

            if !votfound
                tstart= now()
                println("###########################")
                println("## Starting with $votname")
                println("## Starting at $tstart")

                df , dfcart , dfcartnorm = getdata(votdir*"/"*votname)
                m.votname= votname
                cycle, flag= cycle_extraction(df, dfcart, m)
                res,  votfound= updt_votcompleted(fileres, votname , cycle, flag, false)

                tend= now()
                println("## Ending at $tend")
                println("## number of cycle: $cycle , flag:$flag ")
                println("##########################")
                println("##")

                duration= Dates.value(tend-tstart) / (1000*3600)
                totalTime += duration
                meanTime= totalTime / i
                ETA= meanTime * (nfile-i) / 24
                nleft= nfile-i
                ETAstr= @sprintf("%3.3f", ETA) ; durationstr= @sprintf("%3.3f", duration)
                @printf("## %s \n",specialstr("Duration: $durationstr hours","YELLOW"))
                @printf("## %s \n",specialstr("ETA: $ETAstr days","YELLOW"))
                @printf("## %s \n",specialstr("Votable done: $votname","YELLOW"))
                @printf("## %s \n",specialstr("Files analyzed: $i","YELLOW"))
                @printf("## %s \n",specialstr("Files to go: $nleft","YELLOW"))
                println("##\n##")

            end
        end
    end
    println("## Main loop done.")
end
###############################################################################
header_extract()

cd(votdir)
votlist= glob("*.vot")
cd(wdir)

main(votlist,"configAll.ext")

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock