text:
using Adapt

struct Field{X, Y, Z, A, D, G, T, B} <: AbstractDataField{X, Y, Z, A, G, T}
                   data :: D
           architecture :: A
                   grid :: G
    boundary_conditions :: B

    function Field{X, Y, Z}(data::D, arch::A, grid::G, bcs::B) where {X, Y, Z, D, A, G, B}
        T = eltype(grid)
        return new{X, Y, Z, A, D, G, T, B}(data, arch, grid, bcs)
    end
end

"""
    Field(X, Y, Z, [arch = CPU()], grid,
          [ bcs = FieldBoundaryConditions(grid, (X, Y, Z)),
           data = new_data(eltype(grid), arch, grid, (X, Y, Z))])

Construct a `Field` on `grid` with `data` on architecture `arch` with
boundary conditions `bcs`. Each of `(X, Y, Z)` is either `Center` or `Face` and determines
the field's location in `(x, y, z)`.

Example
=======

```jldoctest
julia> using Oceananigans, Oceananigans.Grids

julia> ω = Field(Face, Face, Center, CPU(), RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)))
Field located at (Face, Face, Center)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (3, 3, 3)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=1, Ny=1, Nz=1)
└── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=ZeroFlux, top=ZeroFlux)
```
"""
function Field(X, Y, Z,
               arch::AbstractArchitecture,
               grid::AbstractGrid,
               bcs = FieldBoundaryConditions(grid, (X, Y, Z)),
               data = new_data(eltype(grid), arch, grid, (X, Y, Z)))

    validate_field_data(X, Y, Z, data, grid)

    return Field{X, Y, Z}(data, arch, grid, bcs)
end

# Default CPU architecture
Field(X, Y, Z, grid::AbstractGrid, args...) = Field(X, Y, Z, CPU(), grid, args...)

# Canonical `similar` for AbstractField (doesn't transfer boundary conditions)
Base.similar(f::AbstractField{X, Y, Z, Arch}) where {X, Y, Z, Arch} = Field(X, Y, Z, Arch(), f.grid)

# Type "destantiation": convert Face() to Face and Center() to Center if needed.
destantiate(X) = typeof(X)
destantiate(X::DataType) = X

"""
    Field(L::Tuple, arch, grid, data, bcs)

Construct a `Field` at the location defined by the 3-tuple `L`,
whose elements are `Center` or `Face`.
"""
Field(L::Tuple, args...) = Field(destantiate.(L)..., args...)

#####
##### Special constructors for tracers and velocity fields
#####

"""
    CenterField([arch=CPU()], grid
                bcs = TracerBoundaryConditions(grid),
                data = new_data(eltype(grid), arch, grid, (Center, Center, Center)))

Return a `Field{Center, Center, Center}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function CenterField(arch::AbstractArchitecture,
                     grid,
                     bcs = TracerBoundaryConditions(grid),
                     data = new_data(eltype(grid), arch, grid, (Center, Center, Center)))

    return Field(Center, Center, Center, arch, grid, bcs, data)
end

"""
    XFaceField([arch=CPU()], grid
               bcs = UVelocityBoundaryConditions(grid),
               data = new_data(eltype(grid), arch, grid, (Face, Center, Center)))

Return a `Field{Face, Center, Center}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function XFaceField(arch::AbstractArchitecture,
                    grid,
                    bcs = UVelocityBoundaryConditions(grid),
                    data = new_data(eltype(grid), arch, grid, (Face, Center, Center)))

    return Field(Face, Center, Center, arch, grid, bcs, data)
end

"""
    YFaceField([arch=CPU()], grid,
               bcs = VVelocityBoundaryConditions(grid),
               data = new_data(eltype(grid), arch, grid, (Center, Face, Center)))

Return a `Field{Center, Face, Center}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function YFaceField(arch::AbstractArchitecture,
                    grid,
                    bcs = VVelocityBoundaryConditions(grid),
                    data = new_data(eltype(grid), arch, grid, (Center, Face, Center)))

    return Field(Center, Face, Center, arch, grid, bcs, data)
end

"""
    ZFaceField([arch=CPU()], grid
               bcs = WVelocityBoundaryConditions(grid),
               data = new_data(eltype(grid), arch, grid, (Center, Center, Face)))

Return a `Field{Center, Center, Face}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function ZFaceField(arch::AbstractArchitecture,
                    grid,
                    bcs = WVelocityBoundaryConditions(grid),
                    data = new_data(eltype(grid), arch, grid, (Center, Center, Face)))

    return Field(Center, Center, Face, arch, grid, bcs, data)
end

# Default CPU architectures...
CenterField(grid::AbstractGrid, args...) = CenterField(CPU(), grid, args...)
 XFaceField(grid::AbstractGrid, args...) =  XFaceField(CPU(), grid, args...)
 YFaceField(grid::AbstractGrid, args...) =  YFaceField(CPU(), grid, args...)
 ZFaceField(grid::AbstractGrid, args...) =  ZFaceField(CPU(), grid, args...)

Adapt.adapt_structure(to, field::Field) = Adapt.adapt(to, field.data)


text:
# Successive Projection Algorithm (SPA) for separable NMF
#
#   Reference: N. Gillis and S. A. Vavasis, "Fast and robust recursive
#   algorithms for separable nonnegative matrix factorization,"
#   IEEE Transactions on Pattern Analysis and Machine Intelligence,
#   vol. 36, no. 4, pp. 698-714, 2013.

mutable struct SPA{T}
    obj::Symbol   # objective :mse or :div

    function SPA{T}(;obj=:mse) where T
        obj == :mse || obj == :div || throw(ArgumentError("Invalid value for obj."))
        new{T}(obj)
    end
end

"""
separable_data(m,n,k)
Generate a (m x n) matrix X of nonnegative, separable data
with nonnegative rank k. The rows of X correspond to
observations and the columns of X correspond to features.
The separability condition implies that the columns of H
can be permuted to form a (k x k) diagonal block, and
the sum of the entries of each column of H is at most one.
Thus, the (scaled) columns of W appear in X.
"""
function separable_data(m, n, k)
    W = rand(m, k)

    # impose separability
    V = rand(k, n-k)
    V ./= sum(V, dims=1)
    H = [Matrix(I, k, k) V]
    # permute columns of H
    H = H[:, randperm(n)]

    return W, H
end

# initialization
function spa(X::Matrix{T}, k::Integer; nnls_alg::Tuple{Symbol, Symbol}=(:pivot, :cache)) where T

    # Normalize data so that columns of X sum to one
    R = X ./ sum(X, dims=1)

    # W = R[:,ai], where ai are the "anchor indices"
    # (ai forms the convex hull of columns in R)
    ai = Vector{Int}(undef, k)

    # Add columns of X that are furthest from span(W)
    for j = 1:k
        # Add column with the largest residual
        ai[j] = argmax(vec(sum(R.^2, dims=1)))

        # Project R onto the selected column
        p = R[:,ai[j]]         	# column we're projecting on
        R -= p*(p'*R) ./(p'*p) 	# new residual matrix
    end

    # Estimate W as the anchor columns of X
    W = X[:,ai]

    # Estimate H by non-negative least squares: minimize ||X - W*H||
    H = nonneg_lsq(W, X, alg=:fnnls)
    projectnn!(H)

    return W, H
end

# calculate statistics for result
function solve!(alg::SPA{T}, X, W, H) where T
    if alg.obj == :mse
        objv = convert(T, 0.5) * sqL2dist(X, W*H)
    elseif alg.obj == :div
        objv = gkldiv(X, W*H)
    else
        error("Invalid value for obj.")
    end
    return Result{T}(W, H, 0, true, objv)
end

ext: jl


text:
include("operators.jl")

### Terms
# - A Term is something that, given a state, returns a named tuple (E, hams) with an energy
#   and a list of RFO (for each kpoint).
# - Each term must overload
#     `ene_ops(t, ψ, occ; kwargs...)` -> (E::Real, ops::Vector{RealFourierOperator}).
# - Note that terms are allowed to hold on to references to ψ (eg Fock term),
#   so ψ should not mutated after ene_ops

# The Hamiltonian is defined as half of the gradient of the energy
# with respect to the density matrix sum_n fn |ψn><ψn|.
# In particular, dE/dψn = 2 fn |Hψn> (plus weighting for kpoint sampling)
abstract type Term end

### Builders are objects X that store the term parameters, and produce a
# XTerm <: Term when instantiated with a `basis`

"""
A term with a constant zero energy.
"""
struct NoopTerm <: Term
    basis::PlaneWaveBasis
end
function ene_ops(term::NoopTerm, ψ, occ; kwargs...)
    (E=zero(eltype(term.basis)), ops=[NoopOperator(term.basis, kpoint)
                                      for kpoint in term.basis.kpoints])
end

include("Hamiltonian.jl")

include("kinetic.jl")
include("local.jl")
include("nonlocal.jl")
include("hartree.jl")
include("power_nonlinearity.jl")
include("xc.jl")
include("ewald.jl")
include("psp_correction.jl")
include("entropy.jl")
include("magnetic.jl")

# breaks_symmetries on a term builder answers true if this term breaks
# the symmetries of the lattice/atoms (in which case kpoint reduction
# is invalid)
breaks_symmetries(term_builder::Any) = false
breaks_symmetries(term_builder::Magnetic) = true
breaks_symmetries(term_builder::ExternalFromReal) = true
breaks_symmetries(term_builder::ExternalFromFourier) = true

# forces computes either nothing or an array forces[el][at][α]
function forces(term::Term, ψ, occ; kwargs...)
    nothing  # by default, no force
end
@timing function forces(basis::PlaneWaveBasis, ψ, occ; kwargs...)
    if !any(iszero(kpt.coordinate) for kpt in basis.kpoints)
        @warn "Forces for shifted k-Grids not tested"
    end

    # TODO optimize allocs here
    T = eltype(basis)
    f = [zeros(Vec3{T}, length(positions)) for (type, positions) in basis.model.atoms]
    for term in basis.terms
        f_term = forces(term, ψ, occ; kwargs...)
        if !isnothing(f_term)
            f += f_term
        end
    end
    f
end
forces(scfres) = forces(scfres.basis, scfres.ψ, scfres.occupation, ρ=scfres.ρ)

"""
    compute_kernel(term::Term; kwargs...)

Computes a matrix representation of the full response kernel
(derivative of potential with respect to density) in real space.
"""
function compute_kernel(term::Term; kwargs...)
    false  # "strong zero": false + something = something
end
"""
    apply_kernel(term::Term, dρ; kwargs...)

Computes the potential response to a perturbation dρ in real space
"""
function apply_kernel(term::Term, dρ; kwargs...)
    nothing  # by default, no kernel
end

function compute_kernel(basis::PlaneWaveBasis{T}; kwargs...) where {T}
    k = zeros(T, prod(basis.fft_size), prod(basis.fft_size))
    for term in basis.terms
        k .+= compute_kernel(term; kwargs...)
    end
    k
end

function apply_kernel(basis, dρ; RPA=false, kwargs...)
    dV = RealFourierArray(basis)
    for term in basis.terms
        # Skip XC term if RPA is selected
        RPA && term isa TermXc && continue

        dV_term = apply_kernel(term, dρ; kwargs...)
        if !isnothing(dV_term)
            dV.real .+= dV_term.real
        end
    end
    dV
end


text:
__precompile__(true)

module Graft

using DataFrames
using ProgressMeter

import Base: deepcopy, ==, +, -, |>

# Export DataFrame stuff
export NA

# Package Wide Utilities
include("util.jl")

# Edge Iteration
include("edgeiter.jl")

# SparseMatrixCSC compatibilty
include("sparsematrix.jl")

# Vertex Labelling
include("labelmap.jl")

# Basic Graph Definition
include("graph.jl")

# Subgraph
include("subgraph.jl")

# Graph IO
include("graphio.jl")

# Querying
include("query.jl")

# Graph Algortithms
include("algorithms.jl")

# Operations
include("operations/export.jl")
include("operations/merging.jl")

end # module


text:
"""
    GridData <: AbstractArray

Simulation data specific to a single grid.
"""
abstract type GridData{S,R,T,N} <: StaticArray{S,T,N} end

function (::Type{G})(d::GridData{S,R,T,N}) where {G<:GridData,S,R,T,N}
    args = source(d), dest(d), mask(d), proc(d), opt(d), boundary(d), padval(d),
        sourcestatus(d), deststatus(d)
    G{S,R,T,N,map(typeof, args)...}(args...)
end
function ConstructionBase.constructorof(::Type{T}) where T<:GridData{S,R} where {S,R}
    T.name.wrapper{S,R}
end

GridDataOrReps = Union{GridData, Vector{<:GridData}}

# Array interface
# Base.size(d::GridData{<:Tuple{S}}) where {S} = (Y, X)
# Base.axes(d::GridData) = map(Base.OneTo, size(d))
# Base.eltype(d::GridData{<:Any,<:Any,T}) where T = T
# Base.firstindex(d::GridData) = 1
# Base.lastindex(d::GridData{<:Tuple{S}}) where {S} = Y * X

# Getters
mask(d::GridData) = d.mask
radius(d::GridData{<:Any,R}) where R = R
proc(d::GridData) = d.proc
opt(d::GridData) = d.opt
boundary(d::GridData) = d.boundary
padval(d::GridData) = d.padval
source(d::GridData) = d.source
dest(d::GridData) = d.dest
sourcestatus(d::GridData) = d.sourcestatus
deststatus(d::GridData) = d.deststatus

gridsize(d::GridData) = size(d)
gridsize(A::AbstractArray) = size(A)
gridsize(nt::NamedTuple) = gridsize(first(nt))
gridsize(nt::NamedTuple{(),Tuple{}}) = 0, 0
gridview(d::GridData) = sourceview(d)
sourceview(d::GridData) = view(parent(source(d)), map(a -> a .+ radius(d), axes(d))...)
destview(d::GridData) = view(parent(dest(d)), map(a -> a .+ radius(d), axes(d))...)
source_array_or_view(d::GridData) = source(d) isa OffsetArray ? sourceview(d) : source(d)
dest_array_or_view(d::GridData) = dest(d) isa OffsetArray ? destview(d) : dest(d)

function Base.copy!(grid::GridData{<:Any,R}, A::AbstractArray) where R
    pad_axes = map(ax -> ax .+ R, axes(A))
    copyto!(parent(source(grid)), CartesianIndices(pad_axes), A, CartesianIndices(A))
    _updatestatus!(grid)
    return grid
end
function Base.copy!(A::AbstractArray, grid::GridData{<:Any,R}) where R
    pad_axes = map(ax -> ax .+ R, axes(A))
    copyto!(A, CartesianIndices(A), parent(source(grid)), CartesianIndices(pad_axes))
    return A
end
function Base.copy!(A::AbstractDimArray{T,N}, grid::GridData{<:Any,R}) where {T,N,R}
    copy!(parent(A), grid)
    return A
end
function Base.copy!(grid::GridData{<:Any,R}, A::AbstractDimArray{T,N}) where {R,T,N}
    copy!(grid, parent(A))
    return grid
end
function Base.copy!(dst::GridData{<:Any,RD}, src::GridData{<:Any,RS}) where {RD,RS}
    dst_axes = map(s -> RD:s + RD, gridsize(dst))
    src_axes = map(s -> RS:s + RS, gridsize(src))
    copyto!(parent(source(dst)), CartesianIndices(dst_axes),
            parent(source(src)), CartesianIndices(src_axes)
    )
    return dst
end

@propagate_inbounds Base.getindex(d::GridData{s}, I...) where s = getindex(source(d), I...)
@propagate_inbounds function Base.getindex(d::GridData{s}, i1::Int, I::Int...) where s
    getindex(source(d), i1, I...)
end

# _swapsource => ReadableGridData
# Swap source and dest arrays of a grid
_swapsource(d::Tuple) = map(_swapsource, d)
function _swapsource(grid::GridData)
    src = grid.source
    dst = grid.dest
    @set! grid.dest = src
    @set! grid.source = dst
    srcstatus = grid.sourcestatus
    dststatus = grid.deststatus
    @set! grid.deststatus = srcstatus
    return @set grid.sourcestatus = dststatus
end

# _updatestatus!
# Initialise the block status array.
# This tracks whether anything has to be done in an area of the main array.
_updatestatus!(grid::GridData) = _updatestatus!(opt(grid), grid)
function _updatestatus!(opt::SparseOpt, grid)
    blocksize = 2 * radius(grid)
    src = parent(source(grid))
    for I in CartesianIndices(src)
        # Mark the status block (by default a non-zero value)
        if _isactive(src[I], opt)
            bi = _indtoblock.(Tuple(I), blocksize)
            @inbounds sourcestatus(grid)[bi...] = true
            @inbounds deststatus(grid)[bi...] = true
        end
    end
    return nothing
end
_updatestatus!(opt, grid) = nothing

# _addpadding => OffsetArray{T,N}
# Find the maximum radius required by all rules
# Add padding around the original init array, offset into the negative
# So that the first real cell is still 1, 1
function _addpadding(init::AbstractArray{T,N}, r, padval) where {T,N}
    h, w = size(init)
    paddedsize = h + 4r, w + 2r
    paddedindices = -r + 1:h + 3r, -r + 1:w + r
    sourceparent = fill(convert(eltype(init), padval), paddedsize)
    source = OffsetArray(sourceparent, paddedindices...)
    # Copy the init array to the middle section of the source array
    for j in 1:size(init, 2), i in 1:size(init, 1)
        @inbounds source[i, j] = init[i, j]
    end
    return source
end

# _build_status => (Matrix{Bool}, Matrix{Bool})
# Build block-status arrays
# We add an additional block that is never used so we can
# index into it in the block loop without checking
function _build_status(opt::SparseOpt, source, r)
    hoodsize = 2r + 1
    blocksize = 2r
    nblocs = _indtoblock.(size(source), blocksize) .+ 1
    sourcestatus = zeros(Bool, nblocs)
    deststatus = zeros(Bool, nblocs)
    sourcestatus, deststatus
end
_build_status(opt::PerformanceOpt, init, r) = nothing, nothing

"""
    ReadableGridData <: GridData

    ReadableGridData(grid::GridData)
    ReadableGridData{S,R}(init::AbstractArray, mask, opt, boundary, padval)

Simulation data and storage passed to rules for each timestep.

# Type parameters

- `Y`: number of rows
- `X`: number of columns
- `R`: grid padding radius
"""
struct ReadableGridData{
    S<:Tuple,R,T,N,Sc,D,M,P<:Processor,Op<:PerformanceOpt,Bo,PV,SSt,DSt
} <: GridData{S,R,T,N}
    source::Sc
    dest::D
    mask::M
    proc::P
    opt::Op
    boundary::Bo
    padval::PV
    sourcestatus::SSt
    deststatus::DSt
end
function ReadableGridData{S,R}(
    source::Sc, dest::D, mask::M, proc::P, opt::Op, boundary::Bo,
    padval::PV, sourcestatus::SSt, deststatus::DSt
) where {S,R,Sc<:AbstractArray{T,N},D,M,P,Op,Bo,PV,SSt,DSt} where {T,N}
    ReadableGridData{S,R,T,N,Sc,D,M,P,Op,Bo,PV,SSt,DSt}(
        source, dest, mask, proc, opt, boundary, padval, sourcestatus, deststatus
    )
end
# Generate simulation data to match a ruleset and init array.
@inline function ReadableGridData{S,R}(
    init::AbstractArray, mask, proc, opt, boundary, padval
) where {S,R}
    # We add one extra row and column of status blocks so
    # we dont have to worry about special casing the last block
    if R > 0
        source = _addpadding(init, R, padval)
        dest = _addpadding(init, R, padval)
    else
        if opt isa SparseOpt
            opt = NoOpt()
        end
        source = deepcopy(init)
        dest = deepcopy(init)
    end
    sourcestatus, deststatus = _build_status(opt, source, R)

    grid = ReadableGridData{S,R}(
        source, dest, mask, proc, opt, boundary, padval, sourcestatus, deststatus
    )
    _updatestatus!(grid)
    return grid
end

Base.parent(d::ReadableGridData{S,<:Any,T,N}) where {S,T,N} = SizedArray{S,T,N}(sourceview(d))

"""
    WritableGridData <: GridData

    WritableGridData(grid::GridData)

Passed to rules as write grids, and can be written to directly as an array,
or preferably using `add!` etc. All writes handle updates to SparseOpt()
and writing to the correct source/dest array.

Reads are _always from the source array_, as rules must not be sequential between
cells. This means using e.g. `+=` is not supported, instead use `add!`.
"""
struct WritableGridData{
    S<:Tuple,R,T,N,Sc<:AbstractArray{T,N},D<:AbstractArray{T,N},
    M,P<:Processor,Op<:PerformanceOpt,Bo,PV,SSt,DSt
} <: GridData{S,R,T,N}
    source::Sc
    dest::D
    mask::M
    proc::P
    opt::Op
    boundary::Bo
    padval::PV
    sourcestatus::SSt
    deststatus::DSt
end
function WritableGridData{S,R}(
    source::Sc, dest::D, mask::M, proc::P, opt::Op,
    boundary::Bo, padval::PV, sourcestatus::SSt, deststatus::DSt
) where {S,R,Sc<:AbstractArray{T,N},D<:AbstractArray{T,N},M,P,Op,Bo,PV,SSt,DSt} where {T,N}
    WritableGridData{S,R,T,N,Sc,D,M,P,Op,Bo,PV,SSt,DSt}(
        source, dest, mask, proc, opt, boundary, padval, sourcestatus, deststatus
    )
end

Base.parent(d::WritableGridData{S}) where {S} = SizedArray{S}(destview(d))

### UNSAFE / LOCKS required

# Base.setindex!
# This is not safe for general use.
# It can be used where only identical transformations of a cell
# can happen from any other cell, such as setting all 1s to 2.
@propagate_inbounds function Base.setindex!(d::WritableGridData, x, I...)
    _setindex!(d, proc(d), x, I...)
end
@propagate_inbounds function Base.setindex!(d::WritableGridData, x, i1::Int, I::Int...)
    _setindex!(d, proc(d), x, i1, I...)
end

@propagate_inbounds function _setindex!(d::WritableGridData, opt::SingleCPU, x, I...)
    @boundscheck checkbounds(dest(d), I...)
    @inbounds _setdeststatus!(d, x, I...)
    @inbounds dest(d)[I...] = x
end
@propagate_inbounds function _setindex!(d::WritableGridData, proc::ThreadedCPU, x, I...)
    # Dest status is not threadsafe, even if the
    # setindex itself is safe. So we LOCK
    lock(proc)
    @inbounds _setdeststatus!(d, x, I...)
    unlock(proc)
    @inbounds dest(d)[I...] = x
end

# _setdeststatus!
# Sets the status of the destination block that the current index is in.
# It can't turn of block status as the block is larger than the cell
# But should be used inside a LOCK
_setdeststatus!(d::WritableGridData{<:Any,R}, x, I...) where R =
    _setdeststatus!(d, opt(d), x, I...)
function _setdeststatus!(d::WritableGridData{<:Any,R}, opt::SparseOpt, x, I...) where R
    blockindex = _indtoblock.(I .+ R, 2R)
    @inbounds deststatus(d)[blockindex...] |= !(opt.f(x))
    return nothing
end
_setdeststatus!(d::WritableGridData{<:Any,R}, opt::PerformanceOpt, x, I...) where R = nothing


text:
module MultilaneDriving

using StaticArrays
using Distributions
using Parameters
using AutomotiveDrivingModels
using Base.Test

importall POMDPs

import Distributions: rand, pdf

export
    NoCrashMLDriver,
    MLDriverState,
    HighwayVehicle,
    HighwayDrivingContext,
    Highway,
    AccelSpeed,
    LongAccel,

    MLDPOMDP,
    MLDState,
    MLDObs,

    action_distribution,
    propagate

# package code goes here
include("cached.jl")
include("copula.jl")

include("roadway.jl")

include("vehicles.jl")
include("contexts.jl")
include("safety.jl")
include("idm.jl")
include("drivers.jl")

include("pomdp.jl")

end # module


text:
#-------------DMRG.jl----------------------------------------------------------#
#
# Purpose: A test of a simple DMRG algorithm, following the simple DMRG guide:
#              http://simple-dmrg.readthedocs.io/en/latest/
#
#   Notes: This is done in julia because of the difficultly in the DMRG
#          algorithm, itself.
#------------------------------------------------------------------------------#

# Compatibility necessary for later in the algorithm
using Compat

# Data structure for block and enlarged block for later
# immutable are types that are passed by copying rather than reference
# May change in the future, depending on how efficient we want to be
immutable Block
    length::Int
    basis_size::Int
    operator_dict::Dict{Symbol, AbstractMatrix{Float64}}
end

immutable EnlargedBlock
    length::Int
    basis_size::Int
    operator_dict::Dict{Symbol, AbstractMatrix{Float64}}
end

# Defining global variables
# Variable Declarations
# Specific to the XXZ chain
model_d = 2

# Single site S^z
Sz1 = [0.5 0; 0 -0.5]

# Single site S^+
Sp1 = [0 1; 0 0]

# Single site H, magnetic field
H1 = [0 0; 0 0]
#H1 = [1 1; 1 1]

# Creating initial block, conn is the connection operator for the inside edge of
# of block
initial_block = Block(1, model_d, @compat Dict{Symbol, AbstractMatrix{Float64}}(
    :H => H1,
    :conn_Sz => Sz1,
    :conn_Sp => Sp1,
))

# This is computer science magic -- Deprecated
# Basis_size must be the dimension of each operator matrix
isvalid(block::Union(Block, EnlargedBlock)) =
    # All checks to make sure all interior components are true
    all(op -> size(op) == (block.basis_size, block.basis_size),
                          values(block.operator_dict))

# Function to combine two sites together with kron product
function site_join(Sz1, Sp1, Sz2, Sp2)
    const J = 1.0
    const Jz = 1.0
    combined_site = (J/2) * (kron(Sp1, Sp2') + kron(Sp1', Sp2))
                    + Jz * kron(Sz1, Sz2)
    return(combined_site)
end

# Function to enlarge provided block
function enlarge_block(block::Block)
    mblock = block.basis_size
    Op = block.operator_dict

    # Creating operators for enlarged block
    enlarged_Op_dict = @compat Dict{Symbol, AbstractMatrix{Float64}}(
        :H => kron(Op[:H], speye(model_d)) + kron(speye(mblock), H1)
              + site_join(Op[:conn_Sz], Op[:conn_Sp], Sz1, Sp1),
        :conn_Sz => kron(speye(mblock), Sz1),
        :conn_Sp => kron(speye(mblock), Sp1),
    )

    return EnlargedBlock(block.length + 1, block.basis_size * model_d,
                         enlarged_Op_dict)
end

# Function to perform a single DMRG step
function DMRG_step(sys::Block, env::Block, m::Int)

    # Checking system and environment variables are valid
    @assert isvalid(sys)
    @assert isvalid(env)

    # Enlarging each block
    sys_1 = enlarge_block(sys)
    if sys === env
        env_1 = sys_1
    else
        env_1 = enlarge_block(env)
    end

    # Checking the new system and environment
    @assert isvalid(sys_1)
    @assert isvalid(env_1)

    # Constructing the full superblock hamiltonian
    m_sys_1 = sys_1.basis_size
    m_env_1 = env_1.basis_size
    sys_1_op = sys_1.operator_dict
    env_1_op = env_1.operator_dict
    super_ham = kron(sys_1_op[:H], speye(m_env_1))
                + kron(speye(m_sys_1), env_1_op[:H])
                + site_join(sys_1_op[:conn_Sz], sys_1_op[:conn_Sp],
                            env_1_op[:conn_Sz], env_1_op[:conn_Sp])

    # Modification to superblock hamiltonian so it is hermitian
    super_ham = (super_ham + super_ham') / 2
    (energy,), psi0 = eigs(super_ham, nev=1, which=:SR)

    # Manipulating columns and rows sys and environment
    psi0 = transpose(reshape(psi0, (env_1.basis_size, sys_1.basis_size)))
    rho = Hermitian(psi0 * psi0')

    # Diagonalization and sorting
    fact = eigfact(rho)
    evals, evecs = fact[:values], fact[:vectors]
    permutation = sortperm(evals, rev=true)

    # Build TF matrix from m most significant eigenvectors
    min_m = min(length(evals),m)
    indices = permutation[1:min_m]
    TFmatrix = evecs[:, indices]

    truncation_error = 1 - sum(evals[indices])
    println("truncation error: ", truncation_error)

    # Rotate and truncate
    new_op_dict = Dict{Symbol, AbstractMatrix{Float64}}()
    for (name, op) in sys_1.operator_dict
        new_op_dict[name] = rotate_and_truncate(op, TFmatrix)
    end

    newblock = Block(sys_1.length, min_m, new_op_dict)

    return newblock, energy

end

# transforms the operator to the new (possibly truncated) basis with a
# provided transformation matrix
function rotate_and_truncate(Op, TFmatrix)
    return TFmatrix' * (Op * TFmatrix)
end

# function for the infinite system algorithm
function infinite_sys(L::Int, m::Int)

    # opening outputfile
    out = open("out.dat", "w+")

    block = initial_block
    i = 0
    while 2 * block.length < L
        i = block.length * 2 + 2
        println("L = ", i)
        block, energy = DMRG_step(block, block, m)
        EoL = energy / (block.length * 2)
        println("E/L = ", EoL)
        print(out, i, '\t', EoL, '\n')
        #showcompact(block.operator_dict[:conn_Sz])
        println(size(block.operator_dict[:conn_Sp]))
    end
end

infinite_sys(100, 20)

println("Everything's golden!")

max_issues_repo_name: shiffman/simuleios

max_forks_repo_name: shiffman/simuleios


text:
module WeavePynb
# cf: https://github.com/fredrikekre/Literate.jl/tree/master/src

## Goal here
## have a weave that:
## markdownToHTML: creates web page with questions
## markdownToPynb: creates notebooks, evaluated
## markdownToLaTeXQ: creates latex with questions handled so they can be picked up
## markdownToLaTeX: creates latex files from markdown (XXX??)

using JSON, Mustache, LaTeXStrings, Plots
using Markdown
using Dates
using Base64
using Random

include("evalit.jl")
include("pandoc.jl")
include("questions.jl")
include("formatting.jl")
include("markdown-additions.jl")

include("bootstrap.jl")
include("mmdTomd.jl")
include("markdownToPynb.jl")    # notebook for questions,
include("markdownToHTML.jl")   # for making webpages
include("markdownToRmd.jl")   # for making Rmd files (to run through ...)
include("markdownToLaTeXQ.jl")  # for CSI questions, not of more general usage
include("markdownToLaTeX.jl")

export mmd_to_md
export markdownToPynb
export markdownToHTML
export markdownToRmd
export markdownToLaTeXQ
export markdownToLaTeX

export Verbatim, Invisible, Outputonly, ImageFile, HTMLonly
export alert, warning, note
export example, popup, table
export gif_to_data

end # module


text:
module ConvergenceMHDTest

using Gridap
using Test
using Polynomials: fit

u(x) = VectorValue(x[1]^4+x[2]^4,-x[1]^2-x[2]^3,one(x[1])+sin(x[2]^4))
p(x) = sin(x[1]+x[2]+x[3])
j(x) = VectorValue(x[1]^4+x[2],-x[1]-cos(x[2]),one(x[1])+x[3]^4)
φ(x) = sin(x[1]+x[2]+x[3])
B = VectorValue(1.0,1.0,1.0)

∇u(x) = ∇(u)(x)
Δu(x) = Δ(u)(x)
∇p(x) = ∇(p)(x)
∇φ(x) = ∇(φ)(x)

f_u(x) = (∇u(x)')⋅u(x) - Δu(x) + ∇p(x) - j(x)×B
f_p(x) = (∇⋅u)(x)
f_j(x) = j(x) + ∇φ(x) - u(x)×B
f_φ(x) = (∇⋅j)(x)

g_u(x) = u(x)
g_j(x) = j(x)

eu_l2 = Vector{Float64}()
eu_h1 = Vector{Float64}()
ep_l2 = Vector{Float64}()
ej_l2 = Vector{Float64}()
ej_hdiv = Vector{Float64}()
eφ_l2 = Vector{Float64}()

order = 2
nxs = 2:5
domain = (-0.5,0.5,-0.5,0.5,-0.5,0.5)

for n=nxs
  partition = (n,n,n)
  model = CartesianDiscreteModel(domain,partition)

  labels = get_face_labeling(model)
  add_tag_from_tags!(labels,"dirichlet_u",collect(1:24))
  add_tag_from_tags!(labels,"dirichlet_j",collect(1:24))
  add_tag_from_tags!(labels,"neumann_u",collect(25:26))
  add_tag_from_tags!(labels,"neumann_j",collect(25:26))

  Vu = FESpace(model, ReferenceFE(lagrangian,VectorValue{3,Float64},order);
      conformity=:H1, dirichlet_tags="dirichlet_u")

  Vp = FESpace(model, ReferenceFE(lagrangian,Float64,order-1,space=:P);
      conformity=:L2)

  Vj = FESpace(model, ReferenceFE(raviart_thomas,Float64,order-1);
      conformity=:Hdiv, dirichlet_tags="dirichlet_j")

  Vφ = FESpace(model, ReferenceFE(lagrangian,Float64,order-1,space=:Q);
      conformity=:L2)

  U = TrialFESpace(Vu,g_u)
  P = TrialFESpace(Vp)
  J = TrialFESpace(Vj,g_j)
  Φ = TrialFESpace(Vφ)

  Y = MultiFieldFESpace([Vu, Vp, Vj, Vφ])
  X = MultiFieldFESpace([U, P, J, Φ])

  trian = Triangulation(model)
  degree = 2*(order)
  dΩ = Measure(trian,degree)

  uk = interpolate(u,U)
  function a(X,Y)
    u  , p  , j  , φ   = X
    v_u, v_p, v_j, v_φ = Y

    ∫((∇(u)'⋅uk)⋅v_u + inner(∇(u),∇(v_u)) - p*(∇⋅v_u) - (j×B)⋅v_u +
      (∇⋅u)*v_p +
      j⋅v_j - φ*(∇⋅v_j) - (u×B)⋅v_j +
      (∇⋅j)*v_φ )*dΩ
  end

  function l(Y)
    v_u, v_p, v_j, v_φ = Y

    ∫( v_u⋅f_u + v_p*f_p + v_j⋅f_j + v_φ*f_φ )*dΩ
  end

  btrian_u = BoundaryTriangulation(model,tags=["neumann_u"])
  dΓ_u = Measure(btrian_u,degree)
  nb_u = get_normal_vector(btrian_u)

  function l_Γ_u(Y)
    v_u, v_p, v_j, v_φ = Y

    ∫( v_u⋅(nb_u⋅∇u) - (nb_u⋅v_u)*p)*dΓ_u
  end

  btrian_j = BoundaryTriangulation(model,tags=["neumann_j"])
  dΓ_j = Measure(btrian_j,degree)
  nb_j = get_normal_vector(btrian_j)

  function l_Γ_j(Y)
    v_u, v_p, v_j, v_φ = Y

    ∫( -(v_j⋅nb_j)*φ)*dΓ_j
  end

  b(Y) = l(Y) + l_Γ_u(Y) + l_Γ_j(Y)

  op  = AffineFEOperator(a,b,X,Y)

  xh = solve(op)
  uh, ph, jh, φh = xh

  eu = uh - u
  ep = ph - p
  ej = jh - j
  eφ = φh - φ

  l2(v) = sqrt(sum(∫(v⋅v)*dΩ ))
  h1(v) = sqrt(sum(∫(v⋅v + inner(∇(v),∇(v)) )*dΩ ))
  hdiv(v) = sqrt(sum(∫(v⋅v + inner((∇⋅v),(∇⋅v)) )*dΩ ))

  append!(eu_l2, l2(eu) )
  append!(eu_h1, h1(eu) )
  append!(ep_l2, l2(ep) )
  append!(ej_l2, l2(ej) )
  append!(ej_hdiv, hdiv(ej) )
  append!(eφ_l2, l2(eφ) )

end

slope_eu_l2 = fit([log(x) for x in nxs], [log(y) for y in eu_l2], 1)[end]
slope_eu_h1 = fit([log(x) for x in nxs], [log(y) for y in eu_h1], 1)[end]
slope_ep_l2 = fit([log(x) for x in nxs], [log(y) for y in ep_l2], 1)[end]
slope_ej_l2 = fit([log(x) for x in nxs], [log(y) for y in ej_l2], 1)[end]
slope_ej_hdiv = fit([log(x) for x in nxs], [log(y) for y in ej_hdiv], 1)[end]
slope_eφ_l2 = fit([log(x) for x in nxs], [log(y) for y in eφ_l2], 1)[end]

@test abs(slope_eu_l2 + order + 1.0) < 2e-1
@test abs(slope_eu_h1 + order) < 2e-1
@test abs(slope_ep_l2 + order) < 2e-1
@test abs(slope_ej_l2 + order + 1.0) < 2e-1
@test abs(slope_ej_hdiv + order) < 2e-1
@test abs(slope_eφ_l2 + order) < 2e-1

end #module


text:
function test_encode_decode()
    @test HPack.huffman_encode_bytes(HPack.bytearr("www.example.com")) ==
        [0xf1; 0xe3; 0xc2; 0xe5; 0xf2; 0x3a; 0x6b; 0xa0; 0xab; 0x90; 0xf4; 0xff]

    @test HPack.huffman_decode_bytes(HPack.huffman_encode_bytes(HPack.bytearr("www.example.com"))) == convert(Vector{UInt8}, b"www.example.com")

    @test HPack.decode(HPack.DynamicTable(), IOBuffer([0x82; 0x84])) == [(":method", "GET"), (":path", "/")]

    ## Below are examples copied from Appendix C

    ### C.2.1

    headers = [("custom-key", "custom-header")]
    headers_raw =
        [0x00; 0x0a; 0x63; 0x75; 0x73; 0x74; 0x6f; 0x6d; 0x2d; 0x6b; 0x65; 0x79;
         0x0d; 0x63; 0x75; 0x73; 0x74; 0x6f; 0x6d; 0x2d; 0x68; 0x65; 0x61; 0x64;
         0x65; 0x72]

    @test HPack.encode(HPack.DynamicTable(), headers; huffman=false) == headers_raw

    @test HPack.decode(HPack.DynamicTable(), IOBuffer(headers_raw)) == headers

    ### C.2.2

    headers = [(":path", "/sample/path")]
    headers_raw =
        [0x04; 0x0c; 0x2f; 0x73; 0x61; 0x6d; 0x70; 0x6c; 0x65; 0x2f; 0x70; 0x61;
         0x74; 0x68]

    @test HPack.decode(HPack.DynamicTable(), IOBuffer(HPack.encode(HPack.DynamicTable(), headers; huffman=true))) == headers

    @test HPack.decode(HPack.DynamicTable(), IOBuffer(headers_raw)) == headers
end

function test_client_server()
    ### Server and client examples
    request_headers = [(":method", "GET"),
                       ( ":path", "/"),
                       ( ":scheme", "http"),
                       ( ":authority", "127.0.0.1:9000"),
                       ( "accept", "*/*"),
                       ( "accept-encoding", "gzip, deflate"),
                       ( "user-agent", "HTTP2.jl")]

    response_headers = [(":status", "404"),
                        ("server", "nghttpd nghttpd2/1.10.0"),
                        ("date", "Thu, 02 Jun 2016 19:00:13 GMT"),
                        ("content-type", "text/html; charset=UTF-8")]

    client_dt = HPack.DynamicTable()
    server_dt = HPack.DynamicTable()

    request_headers_raw = HPack.encode(client_dt, request_headers; huffman=true)
    @test HPack.decode(server_dt, IOBuffer(request_headers_raw)) == request_headers

    request_headers_raw = HPack.encode(client_dt, request_headers; huffman=false)
    @test HPack.decode(server_dt, IOBuffer(request_headers_raw)) == request_headers

    response_headers_raw = HPack.encode(server_dt, response_headers; huffman=true)
    @test HPack.decode(client_dt, IOBuffer(response_headers_raw)) == response_headers
end

function test_nghttp_response()
    #### Recorded nghttp server response
    response_headers_raw = [0x8d,0x76,0x90,0xaa,0x69,0xd2,0x9a,0xe4,0x52,0xa9,0xa7,0x4a,0x6b,0x13,0x01,0x5c,0x20,0x5c,0x1f,0x61,0x96,0xdf,0x3d,0xbf,0x4a,0x00,0x4a,0x65,0xb6,0xa5,0x04,0x00,0xb8,0xa0,0x5f,0xb8,0x26,0xee,0x32,0xda,0x98,0xb4,0x6f,0x5f,0x92,0x49,0x7c,0xa5,0x89,0xd3,0x4d,0x1f,0x6a,0x12,0x71,0xd8,0x82,0xa6,0x0e,0x1b,0xf0,0xac,0xf7]

    client_dt = HPack.DynamicTable()
    #HPack.encode(client_dt, request_headers; huffman=false)
    response_headers = HPack.decode(client_dt, IOBuffer(response_headers_raw))
    println("Headers")
    println("======================")
    for header in response_headers
        println(header)
    end
end

test_encode_decode()
test_client_server()
test_nghttp_response()


text:
# https://www.hackerrank.com/challenges/30-conditional-statements/problem

n = parse(Int, readline())

function check_weird(v::Int)
    r_weird = "Weird"
    r_not_weird = "Not Weird"
    is_even = v % 2 == 0
    if is_even
        if 2 ≤ v ≤ 5
            return r_not_weird
        elseif 6 ≤ v ≤ 20
            return r_weird
        else
            return r_not_weird
        end
    else
        return r_weird
    end
end

println(check_weird(n))

max_stars_repo_path: Practice/Tutorials/30 Days of Code/Day 03 Intro to Conditional Statements.jl

max_stars_repo_name: mtaufanr/hackerrank-code-julia

max_issues_repo_name: mtaufanr/hackerrank-code-julia

max_forks_repo_name: mtaufanr/hackerrank-code-julia


text:
###### This file is part of the MendelImpute.jl package.
###### It contains relevant code to decompress a vector of `HaplotypeMosaicPair`
###### to a genotype matrix

"""
    convert_compressed(t<:Real, phaseinfo::String, reffile::String)

Converts `phaseinfo` into a phased genotype matrix of type `t` using the full
reference haplotype panel `H`

# Inputs
- `t`: Type of matrix. If `bool`, genotypes are converted to a `BitMatrix`
- `phaseinfo`: Vector of `HaplotypeMosaicPair`s stored in `.jlso` format
- `reffile`: The complete (uncompressed) haplotype reference file

# Output
+ `X1`: allele 1 of the phased genotype. Each column is a sample. `X = X1 + X2`.
+ `X2`: allele 2 of the phased genotype. Each column is a sample. `X = X1 + X2`.
+ `phase`: the original data structure after phasing and imputation.
+ `sampleID`: The ID's of each imputed person.
+ `H`: the complete reference haplotype panel. Columns of `H` are haplotypes.
"""
function convert_compressed(
    t::Type{T},
    phaseinfo::AbstractString,
    reffile::AbstractString
    ) where T <: Real
    endswith(phaseinfo, ".jlso") || error("phaseinfo does not end with '.jlso'")
    H = convert_ht(Bool, reffile, trans=true, msg="importing reference data...")

    phase_and_sampleID = JLSO.load(phaseinfo)
    phase = phase_and_sampleID[:ph]
    sampleID = phase_and_sampleID[:sampleID]

    X1, X2 = convert_compressed(t, phase, H)
    return X1, X2, phase, sampleID, H
end

"""
    convert_compressed(t<:Real, phaseinfo::Vector{HaplotypeMosaicPair}, H::AbstractMatrix)

Columns of `H` are haplotypes.
"""
function convert_compressed(
    t::Type{T},
    phaseinfo::Vector{HaplotypeMosaicPair},
    H::AbstractMatrix
    ) where T <: Real

    people = length(phaseinfo)
    snps = phaseinfo[1].strand1.length
    M = (t == Bool ? BitArray{2} : Matrix{t})

    X1 = M(undef, snps, people)
    X2 = M(undef, snps, people)
    impute!(X1, X2, H, phaseinfo)

    return X1, X2
end


text:
const CuPosition{T} = NamedTuple{(:x, :y, :z), NTuple{3, T}}

struct CuPositionVector{T} <: AbstractArray{T, 1}
    dim::Int
    vec::CuDeviceVector{T, AS.Global}
end

@inline Base.size(v::CuPositionVector{T}) where T = (v.dim,)

@inline Base.getindex(
    v::CuPositionVector{T},
    i::Int
) where T = (x = v.vec[i], y = v.vec[v.dim + i], z = v.vec[2 * v.dim + i])

@inline Base.setindex!(
    v::CuPositionVector{T},
     ::Any,
     ::Int
) where T = error("setindex! not defined for ", typeof(v))

const CuTriangle{T} = NamedTuple{
    (:v1, :v2, :v3, :normal, :distorig, :area),
    Tuple{CuPosition{T}, CuPosition{T}, CuPosition{T}, CuPosition{T}, T, T}
}

struct CuTriangleVector{T} <: AbstractArray{T, 1}
    dim::Int
    vec::CuDeviceVector{T, AS.Global}
end

@inline Base.size(v::CuTriangleVector{T}) where T = (v.dim,)

@inline function Base.getindex(
    v::CuTriangleVector{T},
    i::Int
) where T
    i = (i - 1) * 14 + 1
    (
        v1 = (x = v.vec[i], y = v.vec[i + 1], z = v.vec[i + 2]),
        v2 = (x = v.vec[i + 3], y = v.vec[i + 4], z = v.vec[i + 5]),
        v3 = (x = v.vec[i + 6], y = v.vec[i + 7], z = v.vec[i + 8]),
        normal = (x = v.vec[i + 9], y = v.vec[i + 10], z = v.vec[i + 11]),
        distorig = v.vec[i + 12],
        area = v.vec[i + 13]
    )
end

@inline Base.setindex!(
    v::CuTriangleVector{T},
     ::Any,
     ::Int
) where T = error("setindex! not defined for ", typeof(v))


text:
@testset "testpattern" begin
    tp = testpattern()
    @test typeof(tp) <: Matrix
    @test eltype(tp) <: RGBA
    @test size(tp) === (300,400)
    tp2 = augment(NoOp())
    @test tp == tp2
end

@testset "safe_rand" begin
    num = @inferred AugmentorBase.safe_rand()
    @test 0 <= num <= 1
    @test typeof(num) <: Float64
    num = @inferred AugmentorBase.safe_rand(2)
    @test all(0 .<= num .<= 1)
    @test typeof(num) <: Vector{Float64}
end

@testset "maybe_copy" begin
    A = [1 2 3; 4 5 6; 7 8 9]
    Ao = OffsetArray(A, (-2,-1))
    @test @inferred(AugmentorBase.maybe_copy(A)) === A
    @test @inferred(AugmentorBase.maybe_copy(Ao)) === Ao
    let Ast = @SMatrix [1 2 3; 4 5 6; 7 8 9]
        @test @inferred(AugmentorBase.maybe_copy(Ast)) === Ast
    end
    let v = view(A, 2:3, 1:2)
        @test @inferred(AugmentorBase.maybe_copy(v)) == A[2:3, 1:2]
        @test typeof(AugmentorBase.maybe_copy(v)) <: Array
    end
    let v = view(OffsetArray(A, (-2,-1)), 0:1, 0:1)
        @test @inferred(AugmentorBase.maybe_copy(v)) == A[2:3, 1:2]
        @test typeof(AugmentorBase.maybe_copy(v)) <: Array
    end
    let v = view(A, IdentityRange(2:3), IdentityRange(1:2))
        @test @inferred(AugmentorBase.maybe_copy(v)) == OffsetArray(A[2:3,1:2], (1, 0))
        @test typeof(AugmentorBase.maybe_copy(v)) <: OffsetArray
    end

    let Ar = reshape(view(A,:,:),1,3,3)
        @test @inferred(AugmentorBase.maybe_copy(Ar)) == reshape(A,1,3,3)
        @test typeof(AugmentorBase.maybe_copy(Ar)) <: Array{Int,3}
    end
    let Ar = reshape(view(Ao,:,:),1,3,3)
        @test @inferred(AugmentorBase.maybe_copy(Ar)) == reshape(A,1,3,3)
        @test typeof(AugmentorBase.maybe_copy(Ar)) <: Array{Int,3}
    end
end

@testset "plain_array" begin
    A = [1 2 3; 4 5 6; 7 8 9]
    @test @inferred(AugmentorBase.plain_array(A)) === A
    @test @inferred(AugmentorBase.plain_array(OffsetArray(A, (-2,-1)))) === A
    let As = sparse(A)
       @test @inferred(AugmentorBase.plain_array(As)) == A
       @test typeof(AugmentorBase.plain_array(As)) <: Array
       Ar = reshape(As, 3, 3, 1)
       @test typeof(Ar) <: Base.ReshapedArray
       @test @inferred(AugmentorBase.plain_array(Ar)) == reshape(A,3,3,1)
       @test typeof(AugmentorBase.plain_array(Ar)) <: Array
    end
    let Ast = @SMatrix [1 2 3; 4 5 6; 7 8 9]
        @test @inferred(AugmentorBase.plain_array(Ast)) === Ast
    end
    let v = view(A, 2:3, 1:2)
        @test @inferred(AugmentorBase.plain_array(v)) == A[2:3, 1:2]
        @test typeof(AugmentorBase.plain_array(v)) <: Array
    end
    let v = view(OffsetArray(A, (-2,-1)), 0:1, 0:1)
        @test @inferred(AugmentorBase.plain_array(v)) == A[2:3, 1:2]
        @test typeof(AugmentorBase.plain_array(v)) <: Array
    end
    let v = view(A, IdentityRange(2:3), IdentityRange(1:2))
        @test @inferred(AugmentorBase.plain_array(v)) == A[2:3, 1:2]
        @test typeof(AugmentorBase.plain_array(v)) <: Array
    end
    let Ar = reshape(view(A,:,:),1,3,3)
        @test @inferred(AugmentorBase.plain_array(Ar)) == reshape(A,1,3,3)
        @test typeof(AugmentorBase.plain_array(Ar)) <: Array{Int,3}
    end
end

@testset "plain_axes" begin
    A = [1 2 3; 4 5 6; 7 8 9]
    @test @inferred(AugmentorBase.plain_axes(A)) === A
    @test @inferred(AugmentorBase.plain_axes(OffsetArray(A, (-2,-1)))) === A
    let v = view(A, 2:3, 1:2)
        @test @inferred(AugmentorBase.plain_axes(v)) === v
    end
    let v = view(OffsetArray(A, (-2,-1)), 0:1, 0:1)
        @test @inferred(AugmentorBase.plain_axes(v)) === v
    end
    let v = view(A, IdentityRange(2:3), IdentityRange(1:2))
        @test @inferred(AugmentorBase.plain_axes(v)) === view(A, 2:3, 1:2)
    end

    let Ar = reshape(view(A,:,:),1,3,3)
        @test @inferred(AugmentorBase.plain_axes(Ar)) === Ar
    end
end

@testset "match_idx" begin
    A = [1 2 3; 4 5 6; 7 8 9]
    @test @inferred(AugmentorBase.match_idx(A, axes(A))) === A
    let img = @inferred AugmentorBase.match_idx(A, (2:4, 2:4))
        @test axes(img) === Base.IdentityUnitRange.((2:4, 2:4))
        @test typeof(img) <: OffsetArray
    end
    let B = view(A,1:3,1:3)
        @test @inferred(AugmentorBase.match_idx(B, axes(B))) === B
    end
    let B = view(A,1:3,1:3)
        img = @inferred(AugmentorBase.match_idx(B, B.indices))
        @test axes(img) === Base.IdentityUnitRange.((1:3, 1:3))
        @test typeof(img) <: OffsetArray
    end
    let img = @inferred AugmentorBase.match_idx(view(A,1:3,1:3), (2:4,2:4))
        @test axes(img) === Base.IdentityUnitRange.((2:4, 2:4))
        @test typeof(img) <: OffsetArray
    end

end

@testset "direct_view" begin
    A = [1 2 3; 4 5 6; 7 8 9]
    Av = view(A, IdentityRange(2:3), IdentityRange(1:2))
    As = view(A, 3:-1:2, 1:1:2)
    @test_throws MethodError AugmentorBase.indirect_axes((), ())
    @test_throws MethodError AugmentorBase.direct_axes((), ())
    @test_throws MethodError AugmentorBase.direct_view(A, ())
    @test @inferred(AugmentorBase.direct_view(A, (2:3,1:2))) === Av
    @test @inferred(AugmentorBase.indirect_view(A, (2:3,1:2))) === Av
    @test @inferred(AugmentorBase.direct_view(A, (IdentityRange(2:3),IdentityRange(1:2)))) === Av
    @test @inferred(AugmentorBase.indirect_view(A, (IdentityRange(2:3),IdentityRange(1:2)))) === Av
    @test @inferred(AugmentorBase.direct_view(A, (3:-1:2,1:1:2))) === As
    @test @inferred(AugmentorBase.indirect_view(A, (3:-1:2,1:1:2))) === As
    @test @inferred(AugmentorBase.direct_view(Av, (3:3,1:2))) === view(Av,3:3,1:2)
    @test @inferred(AugmentorBase.indirect_view(Av, (2:2,1:2))) === view(Av,3:3,1:2)
    @test @inferred(AugmentorBase.direct_view(Av, (IdentityRange(3:3),IdentityRange(1:2)))) === view(Av,3:3,1:2)
    @test @inferred(AugmentorBase.indirect_view(Av, (IdentityRange(2:2),IdentityRange(1:2)))) === view(Av,3:3,1:2)
    @test_throws MethodError AugmentorBase.direct_view(Av, (3:-1:2,2:-1:1))
    @test @inferred(AugmentorBase.indirect_view(Av, (2:-1:1,2:-1:1))) === view(A,3:-1:2,2:-1:1)
end

@testset "vectorize" begin
    @test @inferred(AugmentorBase.vectorize(2)) === 2:2
    @test @inferred(AugmentorBase.vectorize(2.3)) === 2.3:2.3
    @test @inferred(AugmentorBase.vectorize(3:4)) === 3:4
    @test @inferred(AugmentorBase.vectorize(3:1:4)) === 3:1:4
    @test @inferred(AugmentorBase.vectorize(3.0:4)) === 3.0:4
    @test @inferred(AugmentorBase.vectorize(3.0:1:4)) === 3.0:1:4
    @test @inferred(AugmentorBase.vectorize(Base.OneTo(4))) === Base.OneTo(4)
end

@testset "round_if_float" begin
    @test @inferred(AugmentorBase.round_if_float(3,2)) === 3
    @test @inferred(AugmentorBase.round_if_float(3.1111,2)) === 3.11
    @test @inferred(AugmentorBase.round_if_float((3,3.1111),2)) === (3,3.11)
end

@testset "unionrange" begin
    @test_throws MethodError AugmentorBase.unionrange(1:1:2, 4:5)
    @test @inferred(AugmentorBase.unionrange(1:5, 2:6)) === 1:6
    @test @inferred(AugmentorBase.unionrange(2:6, 1:5)) === 1:6
    @test @inferred(AugmentorBase.unionrange(1:2, 4:5)) === 1:5
    @test @inferred(AugmentorBase.unionrange(Base.OneTo(2), 4:5)) === 1:5
    @test @inferred(AugmentorBase.unionrange(1:6, 2:3)) === 1:6
    @test @inferred(AugmentorBase.unionrange(2:3, 1:6)) === 1:6
end

@testset "_2dborder!" begin
    A = rand(2, 5, 6)
    @test @inferred(AugmentorBase._2dborder!(A, 0.)) === A
    s = 0.
    ndim, h, w = size(A)
    for i = 1:h, j = 1:w
        if i == 1 || j == 1 || i == h || j == w
            for d = 1:ndim
                s += abs(A[d, i, j])
            end
        end
    end
    @test s == 0.
end


text:
export subsample

############################################################################################################################
# Base multi source abstract type
abstract type judiMultiSourceVector{T} <: DenseVector{T} end

mutable struct judiMultiSourceException <: Exception
    msg :: String
end

display(ms::judiMultiSourceVector) = println("$(string(typeof(ms))) with $(ms.nsrc) sources")
show(io::IO, ms::judiMultiSourceVector) = print(io, "$(string(typeof(ms))) with $(ms.nsrc) sources")
showarg(io::IO, ms::judiMultiSourceVector, toplevel) = print(io, "$(string(typeof(ms))) with $(ms.nsrc) sources")
summary(io::IO, ms::judiMultiSourceVector) = print(io, string(typeof(ms)))
show(io::IO, ::MIME{Symbol("text/plain")}, ms::judiMultiSourceVector) = println(io, "$(typeof(ms)) with $(ms.nsrc) sources")

# Size and type
eltype(::judiMultiSourceVector{T}) where T = T
size(jv::judiMultiSourceVector) = (jv.nsrc,)
length(jv::judiMultiSourceVector{T}) where {T} = sum([length(jv.data[i]) for i=1:jv.nsrc])

# Comparison
isequal(ms1::judiMultiSourceVector, ms2::judiMultiSourceVector) = ms1 == ms2
==(ms1::judiMultiSourceVector, ms2::judiMultiSourceVector) = all(getfield(ms1, s) == getfield(ms2, s) for s in fieldnames(typeof(ms1)))

isapprox(x::judiMultiSourceVector, y::judiMultiSourceVector; rtol::AbstractFloat=sqrt(eps()), atol::AbstractFloat=0.0) =
    all(isapprox(getfield(x, f), getfield(y, f); rtol=rtol, atol=atol) for f in fieldnames(typeof(x)))

check_compat(ms::Vararg{judiMultiSourceVector, N}) where N = true
check_compat(x::Number, ms::judiMultiSourceVector) = true
check_compat(ms::judiMultiSourceVector, x::Number) = true

# Copy
copyto!(ms::judiMultiSourceVector, a::Vector{Array{T, N}}) where {T, N} = copyto!(ms.data, a)
copyto!(a::Vector{Array{T, N}}, ms::judiMultiSourceVector) where {T, N} = copyto!(a, ms.data)
copy(ms::judiMultiSourceVector{T}) where {T} = begin y = zero(T, ms); y.data = deepcopy(ms.data); y end
deepcopy(ms::judiMultiSourceVector{T}) where {T} = copy(ms)
unsafe_convert(::Type{Ptr{T}}, msv::judiMultiSourceVector{T}) where {T} = unsafe_convert(Ptr{T}, msv.data)

# indexing
IndexStyle(::Type{<:judiMultiSourceVector}) = IndexLinear()
setindex!(ms::judiMultiSourceVector{T}, v::Array{T, N}, i::Integer) where {T, N} = begin ms.data[i] = v; nothing end
getindex(ms::judiMultiSourceVector{T}, i::Integer) where {T} = i > ms.nsrc ? 0 : ms[i:i]
getindex(ms::judiMultiSourceVector{T}, ::Colon) where {T} = vec(ms)
firstindex(ms::judiMultiSourceVector) = 1
lastindex(ms::judiMultiSourceVector) = ms.nsrc
iterate(S::judiMultiSourceVector, state::Integer=1) = state > S.nsrc ? nothing : (S[state], state+1)
# Backward compat subsample
subsample(ms::judiMultiSourceVector, i) = getindex(ms, i)

zero(::Type{T}, x::judiMultiSourceVector) where T = throw(judiMultiSourceException("$(typeof(x)) does not implement zero copy zero(::Type{T}, x)"))
similar(x::judiMultiSourceVector{T}) where T = zero(T, x)
similar(x::judiMultiSourceVector, nsrc::Integer) = nsrc < x.nsrc ? zero(eltype(ET), x; nsrc=nsrc) : zero(eltype(ET), x)
similar(x::judiMultiSourceVector, ::Type{ET}) where ET = zero(eltype(ET), x)
similar(x::judiMultiSourceVector, ::Type{ET}, dims::AbstractUnitRange) where ET = similar(x, ET)
similar(x::judiMultiSourceVector, ::Type{ET}, nsrc::Integer) where ET = nsrc <= x.nsrc ? zero(eltype(ET), x; nsrc=nsrc) : similar(x, ET)

jo_convert(::Type{Array{T, N}}, v::judiMultiSourceVector, B::Bool) where {T, N} = jo_convert(T, v, B)

fill!(x::judiMultiSourceVector, val) = fill!.(x.data, val)
sum(x::judiMultiSourceVector) = sum(sum(x.data))

isfinite(v::judiMultiSourceVector) = all(all(isfinite.(v.data[i])) for i=1:v.nsrc)

conj(A::judiMultiSourceVector{vDT}) where vDT = A
transpose(A::judiMultiSourceVector{vDT}) where vDT = A
adjoint(A::judiMultiSourceVector{vDT}) where vDT = A

maximum(a::judiMultiSourceVector{avDT}) where avDT = max([maximum(a.data[i]) for i=1:a.nsrc]...)
minimum(a::judiMultiSourceVector{avDT}) where avDT = min([minimum(a.data[i]) for i=1:a.nsrc]...)

vec(x::judiMultiSourceVector) = vcat(vec.(x.data)...)

time_sampling(ms::judiMultiSourceVector) = [1 for i=1:ms.nsrc]
############################################################################################################################
# Linear algebra `*`
(msv::judiMultiSourceVector{T})(x::Vector{T}) where {T<:AbstractFloat} = x
(msv::judiMultiSourceVector{T})(x::judiMultiSourceVector{T}) where {T} = x
(msv::judiMultiSourceVector{mT})(x::Vector{T}) where {mT, T<:Array} = begin y = deepcopy(msv); y.data .= x; return y end

function *(J::Union{Matrix{vDT}, joAbstractLinearOperator}, x::judiMultiSourceVector{vDT}) where vDT
    outvec = try J.fop(x) catch; J*vec(x) end
    outdata = try reshape(outvec, size(x.data[1]), x.nsrc) catch; outvec end
    return x(outdata)
end

function *(J::joCoreBlock, x::judiMultiSourceVector{vDT}) where vDT
    outvec = vcat([J.fop[i]*x for i=1:J.l]...)
    outdata = try reshape(outvec, size(x.data[1]), J.l*x.nsrc); catch; outvec end
    return x(outdata)
end

# Propagation input
make_input(ms::judiMultiSourceVector) = throw(judiMultiSourceException("$(typeof(ms)) must implement `make_input(ms, dt)` for propagation"))
make_input(a::Array) = a

############################################################################################################################
# Linear algebra norm/abs/cat...
function norm(a::judiMultiSourceVector{T}, p::Real=2) where T
    if p == Inf
        return maximum([norm(a.data[i], p) for i=1:a.nsrc])
    end
    x = 0.f0
    dt = time_sampling(a)
    for j=1:a.nsrc
        x += dt[j] * norm(a.data[j], p)^p
    end
    return T(x^(1.f0/p))
end

# inner product
function dot(a::judiMultiSourceVector{T}, b::judiMultiSourceVector{T}) where T
	# Dot product for data containers
	size(a) == size(b) || throw(judiMultiSourceException("dimension mismatch: $(size(a)) != $(size(b))"))
	dotprod = 0f0
    dt = time_sampling(a)
	for j=1:a.nsrc
		dotprod += dt[j] * dot(a.data[j], b.data[j])
	end
	return T(dotprod)
end

# abs
function abs(a::judiMultiSourceVector{T}) where T
	b = deepcopy(a)
	for j=1:a.nsrc
		b.data[j] = abs.(a.data[j])
	end
	return b
end

# vcat
vcat(a::Array{<:judiMultiSourceVector, 1}) = vcat(a...)

function vcat(ai::Vararg{<:judiMultiSourceVector{T}, N}) where {T, N}
    N == 1 && (return ai[1])
    N > 2 && (return vcat(ai[1], vcat(ai[2:end]...)))
    a, b = ai
    res = deepcopy(a)
    push!(res, b)
    return res
end

############################################################################################################################
# Diff/cumsum
function cumsum(x::judiMultiSourceVector;dims=1)
    y = deepcopy(x)
    cumsum!(y, x; dims=dims)
    return y
end

function cumsum!(y::judiMultiSourceVector, x::judiMultiSourceVector;dims=1)
    dims == 1 || dims == 2 || throw(judiMultiSourceException("Dimension $(dims) is out of range for a 2D array"))
    h = dims == 1 ? time_sampling(x) : [1f0 for i=1:x.nsrc] # receiver dimension is non-dimensional
    for i = 1:x.nsrc
        cumsum!(y.data[i], x.data[i], dims=dims)
        lmul!(h[i], y.data[i])
    end
    return y
end

function diff(x::judiMultiSourceVector; dims=1)
    # note that this is not the same as default diff in julia, the first entry stays in the diff result
    dims == 1 || dims == 2 || throw(judiMultiSourceException("Dimension $(dims) is out of range for a 2D array"))
    y = 1f0 * x        # receiver dimension is non-dimensional
    h = dims == 1 ? time_sampling(x) : [1f0 for i=1:x.nsrc]
    for i = 1:x.nsrc
        n = size(y.data[i],dims)
        selectdim(y.data[i], dims, 2:n) .-= selectdim(y.data[i], dims, 1:n-1)
        lmul!(1 / h[i], y.data[i])
    end
    return y
end

############################################################################################################################
# Type conversions
tof32(x::Number) = [Float32(x)]
tof32(x::Array{T, N}) where {N, T<:Real} = T==Float32 ? x : Float32.(x)
tof32(x::Array{Array{T, N}, 1}) where {N, T<:Real} = T==Float32 ? x : tof32.(x)
tof32(x::Array{Any, 1}) = try Float32.(x) catch e tof32.(x) end
tof32(x::StepRangeLen) = tof32.(x)
tof32(x::Array{StepRangeLen}) = tof32.(x)


text:
using Catlab.CategoricalAlgebra.CSets
using Catlab.Theories
using Catlab.Theories: adom, attr, attrtype, attr, adom, acodom
using Catlab.Present

using AutoHashEquals
using PermutationGroups

# Color assigned to each elem of each compoennt
const CDict = Dict{Symbol, Vector{Int}}

"""    pseudo_cset(g::ACSet{CD,AD})::Tuple{CSet, Vector{Vector{Any}}} where {CD, AD}

To compute automorphisms of Attributed CSets, we create a pseudo CSet which has
additional components for each data type.

This is inefficient for attributes which have a total order on them
(e.g. integers/strings) since we solve for a canonical permutation of the
attributes. Future work could address this by initializing the coloring with
the 'correct' canonical order.
"""
function pseudo_cset(g::StructACSet{S}
                    )::Tuple{StructACSet, Dict{Symbol,Vector{Any}}} where {S}
  tabs, arrs = collect(ob(S)), collect(hom(S))
  src, tgt = dom(S), codom(S)
  dtabs, darrs = collect(attrtype(S)), collect(attr(S))
  dsrc, dtgt = adom(S), acodom(S)

  # Create copy of schema (+ an extra component for each datatype)
  pres = Presentation(FreeSchema)
  xobs = [Ob(FreeSchema,t) for t in vcat(tabs,dtabs)]
  xobsdic = Dict([t=>Ob(FreeSchema,t) for t in vcat(tabs,dtabs)])
  n = length(tabs)
  for x in xobs add_generator!(pres, x) end

  for (arr, s, t) in zip(map(collect, [arrs, src, tgt])...)
    add_generator!(pres, Hom(arr, xobsdic[s], xobsdic[t]))
  end

  # Add an arrow for each attr and store all possible Data values
  attrvals = Dict([t=>Set() for t in dtabs])
  for (arr, s, t) in zip(darrs, dsrc, dtgt)
    add_generator!(pres, Hom(arr, xobsdic[s], xobsdic[t]))
    union!(attrvals[t], Set(g[arr]))
  end

  # Use Julia ordering to give each value an index
  attrvals = Dict([k=>sort(collect(x)) for (k,x) in collect(attrvals)])

  # Create and populate pseudo-cset
  name = Symbol("Pseudo_$(typeof(g).name.name)")
  expr = struct_acset(name, StructACSet, pres, index=vcat(arrs,darrs))
  eval(expr)
  csettype = eval(name)
  res = Base.invokelatest(eval(name))

  # Copy the original non-attr data
  for t in tabs
    add_parts!(res, t, nparts(g,t))
  end
  for a in arrs
    set_subpart!(res, a, g[a])
  end

  # initialize pseudo data components
  for t in dtabs
    add_parts!(res, t, length(attrvals[t]))
  end

  # Replace data value with an index for each attribute
  for (a,t) in zip(darrs, dtabs)
    fks = [findfirst(==(v), attrvals[t]) for v in g[a]]
    set_subpart!(res, a, fks)
  end

  return res, attrvals
end

"""  pseudo_cset_inv(g::CSet, orig::ACSet{CD,AD}, attrvals::Vector{Vector{Any}})::ACSet{CD,AD} where {CD,AD}

Inverse of pseudo_cset. Requires mapping (generated by `pseudo_cset`) of indices
for each Data to the actual data values.
"""
function pseudo_cset_inv(g::StructACSet,
                         orig::StructACSet{S},
                         attrvals::Dict{Symbol,Vector{Any}}
                        )::StructACSet{S} where {S}
  orig = deepcopy(orig)
  arrs = hom(S)
  darrs, dtabs = attr(S), attrtype(S)
  for arr in arrs
    set_subpart!(orig, arr, g[arr])
  end
  for (darr,tgt) in zip(darrs, dtabs)
    set_subpart!(orig, darr, attrvals[tgt][g[darr]])
  end
  return orig
end

# The maximum color of an empty color list is 0
function max0(x::Vector{Int})::Int
  return isempty(x) ? 0 : maximum(x)
end

function check_auto(x::CDict)::Bool
  return all(map(Base.isperm, values(x)))
end

"""    apply_automorphism(c::CSet{S},d::CDict)::CSet{S} where {S}

Apply a coloring to a Cset to get an isomorphic cset
"""
function apply_automorphism(c::StructACSet{S}, d::CDict)::StructACSet{S} where {S}
  check_auto(d) || error("received coloring that is not an automorphism: $d")
  new = deepcopy(c)
  tabs, arrs, srcs, tgts = ob(S), hom(S), dom(S), codom(S)
  for (arr, src, tgt) in zip(arrs,srcs,tgts)
    set_subpart!(new, d[src],arr, d[tgt][c[arr]])
  end
  return new
end

"""    canonical_iso(g::StructACSet)::StructACSet

Lexicographic minimum of all automorphisms
"""
function canonical_iso(g::StructACSet{S})::StructACSet{S} where {S}
  isos = sort([apply_automorphism(g, Dict(a)) for a in autos(g)], by=string)
  return isempty(isos) ? g : isos[1]
end

"""
Compute automorphisms for the pseudo-cset, but then substitute in
the actual attribute values before evaluating the lexicographic order
"""
function canonical_iso(g::StructACSet{S})::StructACSet{S} where {S}
  p, avals = pseudo_cset(g)
  isos = sort([pseudo_cset_inv(apply_automorphism(p, Dict(a)), g, avals)
               for a in autos(p)], by=string)
  return isempty(isos) ? g : isos[1]
end

"""    canonical_hash(g::ACSet)::UInt64

Hash of canonical isomorphism.
"""
function canonical_hash(g::StructACSet)::UInt64
  return hash(string(canonical_iso(g)))
end

"""Data for an individual component (each vector corresponds to its elements)
1.) how many of each color (for each in-arrow) targets this point
2.) what color this point targets (for each out arrow)

This could be extended to add extra automorphism-invariant properties.
"""
@auto_hash_equals struct CDataPoint
  indata::Vector{Vector{Int}}
  outdata::Vector{Int}
end

# Data required to color a CSet (each element of each component)
const CData = Dict{Symbol, Vector{CDataPoint}}

"""
Computes colors for a CSet, distinguishing nodes by their immediate
connectivity. It is not sufficient to compute the automorphism group, but it is
a good starting point.

This does not generalize to ACSets. We cannot naively throw the attributes as
raw data into the color data. It will make indistinguishable elements (e.g. two
elements that map to different data but otherwise can be permuted) as
distinguishable.
"""
function compute_color_data(g::StructACSet{S}, color::CDict)::CData where {S}
  tabs, arrs, srcs, tgts = ob(S), hom(S), dom(S), codom(S)
  ntab = eachindex(tabs)
  res = CData()
  a_in  = [[a_ind for (a_ind, t_ind) in enumerate(tgts) if t_ind==i]
           for i in ntab]
  a_out = [[a_ind for (a_ind, s_ind) in enumerate(srcs) if s_ind==i]
           for i in ntab]

  for (tgt_ind, tgt) in enumerate(tabs)
    subres = Vector{Vector{Int}}[]  # for each table
    for arr_ind in a_in[tgt_ind]
      src = tabs[srcs[arr_ind]]
      color_src = color[src]
      n_color_src = max0(color_src)
      preimg = g.indices[arrs[arr_ind]]
      subsubres = Vector{Int}[]  # for particular elem in tgt
      for tgt_elem in 1:nparts(g, tgt)
        precolor_elem = color_src[preimg[tgt_elem]]
        n_precolor = [count(==(x), precolor_elem) for x in 1:n_color_src]
        push!(subsubres, n_precolor)
      end
      push!(subres, subsubres)
    end

    # Also compute per-element data for table `tgt` (now, regard as a src)
    outgoing_arrows = a_out[tgt_ind]
    out_subres = Vector{Int}[color[tabs[tgts[oga]]][g[arrs[oga]]]
                             for oga in outgoing_arrows]

    # Combine the two pieces of data for each elmeent in tgt, store in res
    res[tgt] = [CDataPoint([ssr[i] for ssr in subres],
                           [osr[i] for osr in out_subres])
          for i in 1:nparts(g,tgt)]
  end
  return res
end

"""
Construct permutation σ⁻¹ such that σσ⁻¹=id
"""
function invert_perm(x::CDict)::CDict
  return Dict([k=>Base.invperm(v) for (k, v) in collect(x)])
end

"""
Compose permutations
"""
function compose_perms(x::CDict, y::CDict)::CDict
  function compose_comp(xs::Vector{Int},ys::Vector{Int})::Vector{Int}
    return [ys[xs[i]] for i in eachindex(xs)]
  end
  return Dict([k=>compose_comp(v1,y[k]) for (k, v1) in collect(x)])
end

"""    color_refine(g::CSet{CD}; init_color::Union{Nothing,CDict}=nothing)::CDict where {CD}

Iterative color refinement based on the number (and identity) of incoming and
outgoing arrows.
"""
function color_refine(g::StructACSet{S};
                      init_color::Union{Nothing,CDict}=nothing
                     )::CDict where {S}
  # Default: uniform coloring
  new_color = (init_color === nothing
              ? CDict([k => ones(Int, nparts(g, k)) for k in ob(S)])
              : init_color)

  prev_n, curr_n, iter = 0, 1, 0
  while prev_n != curr_n
    iter += 1
    prev_color = new_color
    # All that matters about newdata's type is that it is hashable
    newdata = compute_color_data(g, prev_color)
    # Distinguish by both color AND newly computed color data
    new_datahash = Dict{Symbol, Vector{UInt}}(
      [k=>map(hash, zip(prev_color[k],v)) for (k, v) in collect(newdata)])
    # Identify set of new colors for each component
    hashes = Dict{Symbol, Vector{UInt}}(
      [k=>sort(collect(Set(v))) for (k, v) in new_datahash])
    # Assign new colors by hash value of color+newdata
    new_color = CDict([
      k=>[findfirst(==(new_datahash[k][i]), hashes[k])
          for i in 1:nparts(g, k)]
      for (k, v) in new_datahash])
    prev_n = sum(map(max0, values(prev_color)))
    curr_n = sum(map(max0, values(new_color)))
  end
  return new_color
end

"""
Find index at which two vectors diverge (used in `search_tree`)
"""
function common(v1::Vector{T}, v2::Vector{T})::Int where {T}
  for (i, (x, y)) in enumerate(zip(v1, v2))
    if x != y
      return i-1
    end
  end
  return i
end

"""Convert sorted permutation into a single permutation"""
function to_perm(d::CDict, csum::Vector{Int})::Perm
  res = vcat([[vi + offset for vi in v]
              for ((_, v), offset) in zip(sort(collect(d)), csum)]...)
  return Perm(res)
end

"""Convert single permutation back into sorted permutation"""
function from_perm(p::Perm, syms::Vector{Symbol}, lens::Vector{Int},
                   csum::Vector{Int})::CDict
  res = CDict()
  for (k, l, off) in zip(syms, lens, csum)
    res[k] = [p[i+off]-off for i in 1:l]
  end
  return res
end

"""Enumerate the elements of a permutation group from its generators"""
function all_perms(perm_gens::Vector{CDict})::Set{CDict}
  syms = sort(collect(keys(perm_gens[1])))
  lens = map(length, [perm_gens[1][x] for x in syms])
  csum = vcat([0], cumsum(lens)[1:end-1])
  _,_,Cs = schreier_sims([to_perm(g, csum) for g in perm_gens])
  result = Set{CDict}()

  """
  Recursive function based on (sub)chain C and partial product r
  Algorithm from Fig II.1 of Alexander Hulpke's notes on Computational
  Group Theory (Spring 2010).
  """
  function enum(i, r)
    leaf = i == length(Cs)
    C = Cs[i]
    for d ∈ C.orb
      xr = C[d]*r
      leaf ? push!(result, from_perm(xr, syms, lens, csum)) : enum(i+1, xr)
    end
  end

  enum(1, Perm(1:sum(lens)))

  return result
end

"""
DFS tree of colorings, with edges being choices in how to break symmetry
Goal is to acquire all leaf nodes.

Algorithm from "McKay’s Canonical Graph Labeling Algorithm" by Hartke and
Radcliffe (2009).

Note, there is another use of automorphisms not yet implemented:
 -"Automorphisms discovered during the search can also be used to prune the
   search tree in another way. Let d be a node being re-visited by the depth
   first search. Let Γ be the group generated by the automorphisms discovered
   thus far, and let Φ be the subgroup of Γ that fixes d. Suppose that b and c
   are children of d where some element of Φ maps b to c. If T(G,b) has already
   been examined, then, as above, there is no need to examine T(G,c). Hence
   T(G,c) can be pruned from the tree."
"""
function search_tree!(g::StructACSet{S}, res::Vector{CDict},
                     coloring::CDict,
                     split_seq::Vector{Int},
                     tree::Dict{Vector{Int},CDict},
                     perms::Set{Vector{Int}},
                     skip::Set{Vector{Int}}
                    )::Nothing where {S}

  tree[split_seq] = coloring # add the current color to the tree

  # To reduce branching factor, split on the SMALLEST nontrivial partition
  colors_by_size = []
  for (k, v) in coloring
    for color in 1:max0(v)
      n_c = count(==(color), v)
      if n_c > 1
        # Store which table and which color
        push!(colors_by_size, n_c => (k, color))
      end
    end
  end

  if isempty(colors_by_size) # We found a leaf!
    # Construct automorphisms between leaves
    # to possibly prune the search tree. See Figure 4
    tau_inv = invert_perm(coloring)
    for p in perms
      pii = tree[p]
      auto = compose_perms(pii,tau_inv)
      i = common(p, split_seq)
      a = tree[p[1:i]]
      if compose_perms(auto, a) == a
        b = tree[p[1:i+1]]
        c_location = split_seq[1:i+1]
        c = tree[c_location]
        if compose_perms(auto, b) == c
          push!(skip, c_location)
        end
      end
    end
    # Add permutation to the list of results
    push!(perms, split_seq)
    push!(res, coloring)
  else
    sort!(colors_by_size)
    split_tab, split_color = colors_by_size[1][2]
    colors = coloring[split_tab]
    split_inds = findall(==(split_color), colors)
    for split_ind in split_inds
      if  !(split_seq in skip)
        new_coloring = deepcopy(coloring)
        new_seq = vcat(split_seq, [split_ind])
        new_coloring[split_tab][split_ind] = max0(colors) + 1
        refined = color_refine(g; init_color=new_coloring)
        search_tree!(g, res, refined, new_seq, tree, perms, skip)
      end
    end
  end
  return nothing
end

"""    autos(g::CSet)::Vector{CDict}

Compute the automorphisms of a CSet
"""
function autos(g::StructACSet)::Set{CDict}
  res = CDict[]
  search_tree!(g, res, color_refine(g), Int[],
             Dict{Vector{Int},CDict}(),
             Set{Vector{Int}}(),
             Set{Vector{Int}}())
  return all_perms(res)
end


text:
function init_preconditioner(d, q, elType=typeof(1.0))
    D = d * (q + 1)
    P = Diagonal(ones(elType, D))
    PI = Diagonal(ones(elType, D))
    return P, PI
end

function make_preconditioners!(cache::GaussianODEFilterCache, dt)
    @unpack P, PI, d, q = cache
    make_preconditioner!(P, dt, d, q)
    make_preconditioner_inv!(PI, dt, d, q)
    return nothing
end
function make_preconditioners!(post::GaussianODEFilterPosterior, dt)
    @unpack P, PI, d, q = post
    make_preconditioner!(P, dt, d, q)
    make_preconditioner_inv!(PI, dt, d, q)
    return nothing
end

@fastmath @inbounds function make_preconditioner!(P, h, d, q)
    val = factorial(q) / h^(q + 1 / 2)
    for j in 0:q
        @simd for i in 0:d-1
            P[j+i*(q+1)+1, j+i*(q+1)+1] = val
        end
        val /= (q - j) / h
    end
    return P
end

@fastmath @inbounds function make_preconditioner_inv!(PI, h, d, q)
    val = h^(q + 1 / 2) / factorial(q)
    for j in 0:q
        @simd for i in 0:d-1
            PI[j+i*(q+1)+1, j+i*(q+1)+1] = val
        end
        val *= (q - j) / h
    end
    return PI
end


text:
"""
    AbstractZero <: AbstractDifferential

Supertype for zero-like differentials—i.e., differentials that act like zero when
added or multiplied to other values.
If an AD system encounters a propagator that takes as input only subtypes of `AbstractZero`,
then it can stop performing AD operations.
All propagators are linear functions, and thus the final result will be zero.

All `AbstractZero` subtypes are singleton types.
There are two of them: [`Zero()`](@ref) and [`DoesNotExist()`](@ref).
"""
abstract type AbstractZero <: AbstractDifferential end
Base.iszero(::AbstractZero) = true

Base.iterate(x::AbstractZero) = (x, nothing)
Base.iterate(::AbstractZero, ::Any) = nothing

Base.Broadcast.broadcastable(x::AbstractZero) = Ref(x)
Base.Broadcast.broadcasted(::Type{T}) where T<:AbstractZero = T()

# Linear operators
Base.adjoint(z::AbstractZero) = z
Base.transpose(z::AbstractZero) = z
Base.:/(z::AbstractZero, ::Any) = z

Base.convert(::Type{T}, x::AbstractZero) where T <: Number = zero(T)

"""
    Zero() <: AbstractZero

The additive identity for differentials.
This is basically the same as `0`.
A derivative of `Zero()` does not propagate through the primal function.
"""
struct Zero <: AbstractZero end

extern(x::Zero) = false  # false is a strong 0. E.g. `false * NaN = 0.0`

Base.eltype(::Type{Zero}) = Zero

Base.zero(::AbstractDifferential) = Zero()
Base.zero(::Type{<:AbstractDifferential}) = Zero()

"""
    DoesNotExist() <: AbstractZero

This differential indicates that the derivative does not exist.
It is the differential for primal types that are not differentiable,
such as integers or booleans (when they are not being used to represent
floating-point values).
The only valid way to perturb such values is to not change them at all.
As a consequence, `DoesNotExist` is functionally identical to `Zero()`,
but it provides additional semantic information.

Adding this differential to a primal is generally wrong: gradient-based
methods cannot be used to optimize over discrete variables.
An optimization package making use of this might want to check for such a case.

!!! note:
    This does not indicate that the derivative is not implemented,
    but rather that mathematically it is not defined.

This mostly shows up as the derivative with respect to dimension, index, or size
arguments.
```
    function rrule(fill, x, len::Int)
        y = fill(x, len)
        fill_pullback(ȳ) = (NO_FIELDS, @thunk(sum(Ȳ)), DoesNotExist())
        return y, fill_pullback
    end
```
"""
struct DoesNotExist <: AbstractZero end

function extern(x::DoesNotExist)
    throw(ArgumentError("Derivative does not exit. Cannot be converted to an external type."))
end


text:
export getLinearInterMatrixTranspose

function findValid(js::Array{Int,2},n::Array{Int})
	nv  =  size(js,1)
	val = BitArray(nv)
	for k=1:nv
	 val[k] =  (js[k,1]>0) && (js[k,1] <=n[1]) && (js[k,2] > 0) && (js[k,2]<=n[2])
	end
	return val
end

"""
function [A,dATc] = getLinearInterMatrixTranspose(domain,m,x,T)

builds transpose of linear interpolation matrix A,
                           s.t. A(x)'*T(:) = linearInter(T,domain,x)

Input:
 domain    - description of computational domain
 m        - discretization size of image
 x        - interpolation points
 T        - interpolation coefficients
 varargin - optional parameters

Output:
 A     - linear interpolation matrix
 dAT   - derivative of (A'*T) w.r.t. x

"""
function getLinearInterMatrixTranspose(domain::Array{Float64},n::Array{Int},x::Array{Float64};doDerivative::Bool=false,Img::Vector{Float64}=zeros(0))

h   = vec(domain[2:2:end]-domain[1:2:end])./vec(n)
dim = round(Int,length(domain)/2)
np  = round(Int,length(x)/dim)
x   = reshape(copy(x),np,dim)

# Convert x and y to the coordinate system 1:m, 1:n
for i=1:dim
    x[:,i] = (x[:,i]-domain[2*i-1])/h[i] + .5
end

Valid(j) = trues(size((x[:,j])));#; (0<x(:,j) & x(:,j)<m(j)+1);      # determine indices of valid points

if dim==1
	valid = find( Valid(1) )
elseif dim===2
	valid = find( Valid(1) .& Valid(2) )
elseif dim==3
	valid = find( Valid(1) .& Valid(2) .& Valid(3) )
end

if isempty(valid)
    return spzeros(np,prod(n))
end;

P     = round.(Int,floor.(x)); x = x-P;    # split x into integer/remainder
p(j)  = P[valid,j]
xi(j) = x[:,j]

if dim==1
	error("nyi")
elseif dim==2
	ind  = valid

	# get indices of neighbouring points
	jk    = [p(1)   p(2)  ]
	j1k   = [p(1)+1 p(2)  ]
	jk1   = [p(1)   p(2)+1]
	j1k1  = [p(1)+1 p(2)+1]

	# compute weights for interpolation
	s1 = (1-xi(1)).*(1-xi(2))
	s2 = (xi(1)).*(1-xi(2))
	s3 = (1-xi(1)).*(xi(2))
	s4 =  (xi(1)).*(xi(2))
	ii = [ind  ind  ind  ind]
	jj = [jk ; j1k; jk1; j1k1]

	# points that are outside the domain
	# valid = (js::Array{Int,2}) -> find((js[:,1].>0)&(js[:,1].<=n[1])&(js[:,2].>0)&(js[:,2].<=n[2]))
	# valid = (js::Array{Int,2}) -> (js[:,1].>0)&(js[:,1].<=n[1])&(js[:,2].>0)&(js[:,2].<=n[2])
    valid = (ij::Array{Int,2}) -> findValid(ij,n)
    s2ind = (jj::Array{Int,2}) -> jj[:,1]+(jj[:,2]-1)*n[1]

        jkt = valid(jk); ijkt = ind[jkt]
		A1 = sparse(s2ind(jk[jkt,:])     ,ijkt   ,s1[ijkt]  ,prod(n),np);
        A2 = sparse(s2ind(j1k[valid(j1k),:])   ,ind[valid(j1k)]  ,s2[ind[valid(j1k)]] ,prod(n),np);
        A3 = sparse(s2ind(jk1[valid(jk1),:])   ,ind[valid(jk1)]  ,s3[ind[valid(jk1)]] ,prod(n),np);
        A4 = sparse(s2ind(j1k1[valid(j1k1),:]) ,ind[valid(j1k1)] ,s4[ind[valid(j1k1)]],prod(n),np);

        A =  A1 + A2 + A3 + A4;
        if doDerivative
			if isempty(Img); error("need to provide direction"); end

            dv1dxi = zeros(np); dv1deta = zeros(np)
            dv1dxi[ind] = -(1-xi(2)); dv1deta[ind] = -(1-xi(1))

            dv2dxi = zeros(np); dv2deta = zeros(np);
            dv2dxi[ind] = (1-xi(2)); dv2deta[ind] = -(xi(1));

            dv3dxi = zeros(np); dv3deta = zeros(np);
            dv3dxi[ind] = (-xi(2)); dv3deta[ind] = 1-xi(1);

            dv4dxi = zeros(np); dv4deta = zeros(np);
            dv4dxi[ind] = (xi(2)); dv4deta[ind] = xi(1);
            # built sdiags before creating function handle
            dv1dxi = sdiag(dv1dxi);
            dv2dxi = sdiag(dv2dxi);
            dv3dxi = sdiag(dv3dxi);
            dv4dxi = sdiag(dv4dxi);

			A1 = sparse(s2ind(jk[valid(jk),:])     ,ind[valid(jk)]   ,Img[ind[valid(jk)]]  ,prod(n),np);
	        A2 = sparse(s2ind(j1k[valid(j1k),:])   ,ind[valid(j1k)]  ,Img[ind[valid(j1k)]] ,prod(n),np);
	        A3 = sparse(s2ind(jk1[valid(jk1),:])   ,ind[valid(jk1)]  ,Img[ind[valid(jk1)]] ,prod(n),np);
	        A4 = sparse(s2ind(j1k1[valid(j1k1),:]) ,ind[valid(j1k1)] ,Img[ind[valid(j1k1)]],prod(n),np);

			Dxi =  A1*dv1dxi + A2*dv2dxi +  A3*dv3dxi + A4*dv4dxi

            # built sdiags before creating function handle
            dv1deta = sdiag(dv1deta);
            dv2deta = sdiag(dv2deta);
            dv3deta = sdiag(dv3deta);
            dv4deta = sdiag(dv4deta);
            Deta =  A1*dv1deta + A2*dv2deta + A3*dv3deta + A4*dv4deta
            dAT =  [Dxi/h[1]  Deta/h[2]];
		end
   elseif dim==3
	error("NYI")
	end
if !doDerivative
	return A
else
	return A, dAT;
end
end


text:
# Find and add MDP's path to Julia search path
cdir = pwd(); ii = findlast("MDP", cdir);
if ii===nothing
    # MDPpath = /path/to/MDP;
    display("MDP's path is not found. Please uncomment the above line and set the path to MDP.");
else
    MDPpath = cdir[1:ii[end]] * "/";
    include(MDPpath * "Installation/setpath.jl");
end

# External packages
using Revise, DelimitedFiles, SymPy

# MDP packages
using Preprocessing, Gencode, Postprocessing

app = Preprocessing.initializeapp(cdir,MDPpath);
app.platform = "cpu";    # app.platform = "gpu" => MDP runs on Nvidia GPU
app.unitstyle = "metal"; # unit system
app.pbc = [1 1 1];       # periodic boundary conditions
nlattice = 4;            # note that # of atoms = 2 * nlattice^3
app.appname = "snapTA06A" * string(nlattice);

#  neighbor list
app.neighskin = 1.0;
app.neighcell = 0;       # 0 -> O(N^2) algorithm, 1 -> Cell-linked list algorithm to form neighbor list
app.neighmax = 100;      # maximum number of neighbors allowed
app.neighevery = 1;      # perform neighbor rebuild list for "every" iterations
app.neighdelay = 0;      # delay neighbor rebuild
app.neighcheck = 1;      # 0 -> neighbor rebuild by delay and every settings,
                         # 1 -> neighbor rebuild occurs if some atom has moved more than half the skin distance
app.neighpair = 0;       # 0 -> full neighbor list, 1 -> half neighbor list for pair potentials

# atom
app.atommasses = reshape([180.88],1,1);
app.atomcharges = reshape([0.0],1,1);

# lattice, region, domain
app.lattice = Preprocessing.setlattice("bcc", 3.316);
app.region = Preprocessing.setregion([nlattice nlattice nlattice]);

# nonparametric potential descriptors
app.descriptor = "snap";    # (snap, shp)
app.snaprcutfac = 4.67637;
app.snaptwojmax = 6;
app.snaprfac0 = 0.99363;
app.snaprmin0 = 0;
app.snapbzeroflag = 0;
app.snapquadraticflag = 0;
app.snapswitchflag = 1;
app.snapchemflag = 0;
app.snapbnormflag = 0;
app.snapwselfallflag = 0;
app.snapnelem = 1;    # # of elements
app.snapncoeff = 31;  # # of coeffients per element
app.snapelemradius = reshape([0.5],1,1); # radius per element
app.snapelemweight = reshape([1.0],1,1);   # weight per element
app.snapcoefffile = "snapTA06A.coeff";
app.rcutml = app.snaprcutfac;       # cut-off radius for nonparametric potential

# parametric potential descriptors
app.potentialfile = "ZBLpotential.jl";   # empirical potential file
include(app.potentialfile);  # include the potential file
zblcutinner = 4;
zblcutouter = 4.8;
zblz = 73;
app.pot2a = reshape([1],1,1);     # list of nonbonded pair potentials
app.mu2a = [zblcutinner zblcutouter zblz];    # parameters for ZBL potential
app.rcut2a = reshape([zblcutouter],1,1);  # cut-off radius

# time integration and outputs
app.ensemblemode = "nve";
app.time = 0.0;             # initial time
app.dt = 0.5e-3;            # time step size
app.ntimesteps = 100;       # # time steps
app.globalfreq = 10;

temp = 300;  seed0 = 4928459; distflag = 0; sumflag = 0; loopflag = 2; momflag = 1; rotflag = 0;
app.createvelocity = [temp seed0 distflag sumflag loopflag momflag rotflag];

# generate input files and run simulation
app, config = mdp(app);

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
module  Utils
    using GLMakie
    using Dates
    using Statistics
    using StatsBase
    using EnergyStatistics

    struct Region
        latitude_min
        latitude_max
        longitude_min
        longitude_max
    end

    struct Point
        latitude
        longitude
    end

    function get_area(loc, width, height)
        offset_x = width / 2
        offset_y = height / 2
        return Region(
            round(loc[1] - offset_x, digits = 6),
            round(loc[1] + offset_x, digits = 6),
            round(loc[2] - offset_y, digits = 6),
            round(loc[2] + offset_y, digits = 6),
        )
    end

    function intersection_over_union(boxA::Region, boxB::Region)
        # determine the (x, y)-coordinates of the intersection rectangle
        xA = max(boxA.latitude_min, boxB.latitude_min)
        yA = max(boxA.longitude_min, boxB.longitude_min)
        xB = min(boxA.latitude_max, boxB.latitude_max)
        yB = min(boxA.longitude_max, boxB.longitude_max)

        # compute the area of intersection rectangle
        interArea = max(0, (xB - xA)) * max(0, (yB - yA))
        if interArea == 0
            return 0
        end

        # compute the area of both the prediction and ground-truth
        # rectangles
        boxAArea = abs((boxA.latitude_max - boxA.latitude_min) * (boxA.longitude_max - boxA.longitude_min))
        boxBArea = abs((boxB.latitude_max - boxB.latitude_min) * (boxB.longitude_max - boxB.longitude_min))

        # compute the intersection over union by taking the intersection
        # area and dividing it by the sum of prediction + ground-truth
        # areas - the interesection area
        iou = interArea / (boxAArea + boxBArea - interArea)

        # return the intersection over union value
        return iou
    end

    #my_pseudolog10(x) = sign(x) * exp(abs(x) + 1)

    function set_colormap!(hm)
        if hm.colorrange[][1] >= 0
            #hm.colorrange[] = (0, 1) .* maximum(abs, hm.colorrange[])
            hm.colormap = cgrad(:amp) #, scale = my_pseudolog10)
        else
            if hm.colorrange[][2] > 0
                hm.colorrange[] = (-1, 1) .* maximum(abs, hm.colorrange[])
                hm.colormap = :balance
            else
                #hm.colorrange[] = (-1, 0) .* maximum(abs, hm.colorrange[])
                hm.colormap = cgrad(:ice) #, scale = my_pseudolog10)
            end
        end

        #println(hm.colorrange[][1], ", ", hm.colorrange[][2])
    end

    function create_heatmap(X, Y, width, height, colName; windDir = false)
        data1 = Array{Float32}(undef, width, height)
        data2 = Array{Float32}(undef, width, height)
        data3 = Array{Float32}(undef, width, height)

        k = 1
        for j::Int in 1:height
            for i::Int in 1:width
                data1[i, j] = cor(X[!, "$(colName)$(k)"], Y[!, "Irradiance"]) * 100
                data2[i, j] = corspearman(X[!, "$(colName)$(k)"], Y[!, "Irradiance"]) * 100
                data3[i, j] = dcor(X[!, "$(colName)$(k)"], Y[!, "Irradiance"]) * 100
                k = k + 1
            end
        end

        fig = Figure(resolution = (1720, 600))

        ax1, hm1 = heatmap(fig[1, 1], data1)
        ax1.title = "$(colName) - Pearson correlation"
        ax1.xlabel = "Longitude"
        ax1.ylabel = "Latitude"
        set_colormap!(hm1)
        Colorbar(fig[1, 2], hm1, label="Percent [%]")
        ax2, hm2 = heatmap(fig[1, 3], data2)
        ax2.title = "$(colName) - Spearman correlation"
        ax2.xlabel = "Longitude"
        ax2.ylabel = "Latitude"
        set_colormap!(hm2)
        Colorbar(fig[1, 4], hm2, label="Percent [%]")
        ax3, hm3 = heatmap(fig[1, 5], data3)
        ax3.title = "$(colName) - Distance correlation"
        ax3.xlabel = "Longitude"
        ax3.ylabel = "Latitude"
        set_colormap!(hm3)
        Colorbar(fig[1, 6], hm3, label="Percent [%]")

        if windDir
            xs = 1:width
            ys = 1:height
            dir = mean(Matrix(X[!, ["WindDirection$(k)" for k::Int in 1:width * height]]), dims=1)
            Wy = -sind.(mod.(-dir .+ 90, 360))
            Wx = -cosd.(mod.(-dir .+ 90, 360))

            arrows!(fig[1, 1], xs, ys, Wx, Wy, arrowsize = 15, lengthscale = 0.2, linewidth = 3)
            arrows!(fig[1, 3], xs, ys, Wx, Wy, arrowsize = 15, lengthscale = 0.2, linewidth = 3)
            arrows!(fig[1, 5], xs, ys, Wx, Wy, arrowsize = 15, lengthscale = 0.2, linewidth = 3)
        end

        save("imgs/$(colName).png", fig)
    end
end


text:
using HDF5,JLD,KUparser,KUnet

info("Load data")
@date dev=load("acl11.dev.jld4")
@date tst=load("acl11.tst.jld4")
@date trn=load("acl11.trn.jld4")
@assert trn["deprel"] == dev["deprel"] == tst["deprel"]
@assert trn["postag"] == dev["postag"] == tst["postag"]
nd = length(dev["deprel"])
pt = ArcEager13
ft = Flist.zn11single
@date (pdev,xdev,ydev) = oparse(pt, dev["corpus"], nd, ft)
@date (ptst,xtst,ytst) = oparse(pt, tst["corpus"], nd, ft)
@date (ptrn,xtrn,ytrn) = oparse(pt, trn["corpus"], nd, ft)
nrows = max(maximum(xtrn.rowval), maximum(xtst.rowval), maximum(xdev.rowval))
xtrn.m = xtst.m = xdev.m = nrows
ytrn = sparse(ytrn); ydev = sparse(ydev); ytst = sparse(ytst)
@save "zn11oparse1.jld" xtrn ytrn xdev ydev xtst ytst
# @save "zn11oparse1_hash.jld" KUparser.SFhash

# @date (pdev,xdev,ydev) = oparse(pt, dev["corpus"][1:100], nd, ft)
# @save "foo.zn11orig.jld" xdev ydev

:ok

# p0 = pt(1,ndeps)
# s1 = corpus[1]

# ft = Array[["s0w"], ["s0p"], ["s0w", "s0p"]]
# @date (p1,x1,y1)=oparse(pt, s1, ndeps, ft)

# SFarray = Array(Any, length(KUparser.SFhash))
# for k in keys(KUparser.SFhash)
#     v = KUparser.SFhash[k]
#     SFarray[v] = k
# end

# for j=1:size(x1,2)
#     println(SFarray[find(x1[:,j])])
# end


text:
# Have to do the funky push! on arrays due to
# the deprecated [] auto-concatenation still
# enabled in Julia 0.4

include("types.jl")

function _wrap(w)

  if ismatch(r"^Array", string(typeof(w)))
    for k in 1:length(w)
      w[k] = _wrap(w[k])
    end
    p = []
    push!(p, 2)
    push!(p, w)
    return p
  end

  if ismatch(r"^Dict", string(typeof(w)))
    for k in Base.keys(w)
      w[k] = _wrap(w[k])
    end
    return w
  end

  w
end

function wrap(o)
  if typeof(o) != ReqlObject && typeof(o) != ReqlSequence
    return o
  end

  _wrap(o.value)
  o
end

function convert_type(from, to)
  if typeof(from) == to
    return from
  end
  to(-1, from)
end

function convert_arr_type(from, to)
  new_arr = []

  for i in 1:length(from)
    if typeof(from[i]) == to
      push!(new_arr, from[i])
    else
      push!(new_arr, to(-1, from[i]))
    end
  end
  new_arr
end

macro reql_zero(op_code::Int, name::Symbol, RET)
  quote
    function $(esc(name))()
      $(RET)($op_code, [])
    end
  end
end

macro reql_one(op_code::Int, name::Symbol, T1, RET)

  quote
    function $(esc(name))(arg1::$(reql_type(eval(T1))))

      arg1 = convert_type(arg1, $T1)

      local retval = []
      push!(retval, wrap(arg1))

      $(RET)($op_code, retval)
    end
  end
end

macro reql_onearr(op_code::Int, name::Symbol, T1, RET)
  quote
    function $(esc(name))(arg1::$(reql_type(eval(T1)))...)

      arg1 = convert_arr_type(arg1, $T1)

      local retval = []
      for i in arg1
        push!(retval, wrap(i))
      end

      $(RET)($op_code, retval)
    end
  end
end

macro reql_one_two(op_code::Int, name::Symbol, T1, T2, RET)
  quote
    function $(esc(name))(arg1::$(reql_type(eval(T1))), arg2::$(reql_type(eval(T2))))

      arg1 = convert_type(arg1, $T1)
      arg2 = convert_type(arg2, $T2)

      local retval = []
      push!(retval, wrap(arg1))
      push!(retval, wrap(arg2))

      $(RET)($op_code, retval)
    end
  end
end

macro reql_one_twoarr(op_code::Int, name::Symbol, T1, T2, RET)
  quote
    function $(esc(name))(arg1::$(reql_type(eval(T1))), arg2::$(reql_arr_type(eval(T2))))

      arg1 = convert_type(arg1, $T1)
      arg2 = convert_arr_type(arg2, $T2)

      local retval = []
      push!(retval, wrap(arg1))
      for i in arg2
        push!(retval, wrap(i))
      end

      $(RET)($op_code, retval)
    end
  end
end

macro reql_one_two_three(op_code::Int, name::Symbol, T1, T2, T3, RET)
  quote
    function $(esc(name))(arg1::$(reql_type(eval(T1))), arg2::$(reql_type(eval(T2))), arg3::$(reql_type(eval(T3))))

      arg1 = convert_type(arg1, $T1)
      arg2 = convert_type(arg2, $T2)
      arg3 = convert_type(arg3, $T3)

      local retval = []
      push!(retval, wrap(arg1))
      push!(retval, wrap(arg2))
      push!(retval, wrap(arg3))

      $(RET)($op_code, retval)
    end
  end
end

max_stars_repo_name: mrcsparker/julia-rethinkdb

max_issues_repo_name: mrcsparker/julia-rethinkdb

max_forks_repo_name: mrcsparker/julia-rethinkdb


text:
#!/usr/bin/julia

# Copyright (C) 2021, 2022 Matheus Fernandes Bigolin <mfrdrbigolin@disroot.org>
# SPDX-License-Identifier: MIT

"Day Nine, Smoke Basin."

#= At first I left this day uncompleted because I lacked the graph theory
 = knowledge to solve part two; but after learning the Dijkstra’s algorithm for
 = day fifteen, I managed to solve this problem in a similar way.
 =#

include("Utils.jl")

using .Utils

"Calculate the sum of the risk levels of all low points in `heightmap`."
function solve1(heightmap :: Matrix{Int})
    risklevel_count = 0

    for ind in CartesianIndices(heightmap)
        neighbors = find_neighbors(heightmap, ind.I)

        if heightmap[ind] < minimum([heightmap[n...] for n in neighbors])
            risklevel_count += heightmap[ind] + 1
        end
    end

    return risklevel_count
end

"Calculate the product of the size of the three largest basins in `heightmap`."
function solve2(heightmap :: Matrix{Int})
    # Mark all locations with height nine as visited.
    unvisited = BitArray{2}([row .!= 9 for row in heightmap])

    current = (1, 1)
    next = Array{Tuple{Int, Int}}([])

    basins = Array{Int}([])
    current_basin = 0

    while count(unvisited) != 0
        unvisited[current...] = false
        current_basin += 1

        neighbors = find_neighbors(heightmap, current)

        filter!(x -> unvisited[x...] && x ∉ next, neighbors)
        push!(next, neighbors...)

        if isempty(next)
            push!(basins, current_basin)
            current_basin = 0

            if count(unvisited) != 0
                current = findfirst(unvisited).I
            end
        else
            current = pop!(next)
        end
    end

    return prod(sort(basins)[end - 2:end])
end

function main()
    usage_and_exit(length(ARGS) != 1)

    data = collect.(readlines(first(ARGS)))
    heightmap = hcat([parse.(Int, l) for l in data]...)

    println(solve1(heightmap))
    println(solve2(heightmap))
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
function load_data(dataset::AbstractString)
    if dataset == "facebook_data"
        return load_facebook_data()
    end
end

function load_facebook_data()
    pathname = joinpath(dirname(dirname(@__FILE__)),"data")
    fname = joinpath(pathname, "facebook_data.txt")
    facebook_adj_matrix = fill(0, (4039, 4039))
    open(fname) do file
        for ln in eachline(file)
            pieces = split(ln,' ', keepempty=false)
            facebook_adj_matrix[(parse(Int64,pieces[1]) + 1), (parse(Int64,pieces[2]) + 1)] = 1
        end
    end

    facebook_adj_matrix = sparse(facebook_adj_matrix)
    facebook_adj_matrix = facebook_adj_matrix + facebook_adj_matrix'
    facebook_adj_matrix = min.(facebook_adj_matrix, 1)
    facebook_adj_matrix = MatrixNetwork(facebook_adj_matrix)

    return facebook_adj_matrix
end


text:
function mdl_tg_ode!(
    device_states,
    output_ode,
    device::PSY.DynamicGenerator{M, S, A, PSY.TGFixed, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS}

    #Update inner vars
    P_ref = PSY.get_ext(device)[CONTROL_REFS][P_ref_index]
    get_inner_vars(device)[τm_var] = P_ref * PSY.get_efficiency(PSY.get_prime_mover(device))

    return
end

function mdl_tg_ode!(
    device_states,
    output_ode,
    device::PSY.DynamicGenerator{M, S, A, PSY.TGTypeI, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS}

    #Obtain references
    ω_ref = PSY.get_ext(device)[CONTROL_REFS][ω_ref_index]
    P_ref = PSY.get_ext(device)[CONTROL_REFS][P_ref_index]

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(device, PSY.TGTypeI)

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    x_g1 = internal_states[1]
    x_g2 = internal_states[2]
    x_g3 = internal_states[3]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(device, PSY.TGTypeI)
    ω = @view device_states[external_ix]

    #Get Parameters
    tg = PSY.get_prime_mover(device)
    R = PSY.get_R(tg)
    Ts = PSY.get_Ts(tg)
    Tc = PSY.get_Tc(tg)
    T3 = PSY.get_T3(tg)
    T4 = PSY.get_T4(tg)
    T5 = PSY.get_T5(tg)
    P_min = PSY.get_P_min(tg)
    P_max = PSY.get_P_max(tg)

    #Compute auxiliary parameters
    P_in = P_ref + (1.0 / R) * (ω_ref - ω[1])

    #Set anti-windup for P_in. #TODO in callbacks
    #if P_in > P_max
    #    P_in = P_max
    #elseif P_in < P_min
    #    P_in = P_min
    #end

    #Compute 3 States TG ODE:
    output_ode[local_ix[1]] = (1.0 / Ts) * (P_in - x_g1)
    output_ode[local_ix[2]] = (1.0 / Tc) * ((1.0 - T3 / Tc) * x_g1 - x_g2)
    output_ode[local_ix[3]] =
        (1.0 / T5) * ((1.0 - T4 / T5) * (x_g2 + (T3 / Tc) * x_g1) - x_g3)

    #Update mechanical torque
    get_inner_vars(device)[τm_var] = x_g3 + (T4 / T5) * (x_g2 + (T3 / Tc) * x_g1)

    return
end

function mdl_tg_ode!(
    device_states,
    output_ode,
    device::PSY.DynamicGenerator{M, S, A, PSY.TGTypeII, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS}

    #Obtain references
    ω_ref = PSY.get_ext(device)[CONTROL_REFS][ω_ref_index]
    P_ref = PSY.get_ext(device)[CONTROL_REFS][P_ref_index]

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(device, PSY.TGTypeII)

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    xg = internal_states[1]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(device, PSY.TGTypeII)
    ω = @view device_states[external_ix]

    #Get Parameters
    tg = PSY.get_prime_mover(device)
    R = PSY.get_R(tg)
    T1 = PSY.get_T1(tg)
    T2 = PSY.get_T2(tg)
    τ_min = PSY.get_τ_min(tg)
    τ_max = PSY.get_τ_max(tg)

    #Compute auxiliary parameters
    τ_m = (1.0 / R) * (T1 / T2) * (ω_ref - ω[1]) + P_ref / 1.0 + xg

    #Set anti-windup for τ_m. NOT WORKING
    #if τ_m > τ_max
    #    τ_m = τ_max
    #elseif τ_m < τ_min
    #    τ_m = τ_min
    #end

    #Compute 1 State TG ODE:
    output_ode[local_ix[1]] = (1.0 / T2) * ((1.0 / R) * (1 - T2 / T1) * (ω_ref - ω[1]) - xg)

    #Update mechanical torque
    get_inner_vars(device)[τm_var] = τ_m

    return
end


text:
# closed-loop-bank-model.jl
# SimLynx.jl Example Closed-Loop Monte Carlo Model

using SimLynx
SimLynx.greet()

using Distributions: Exponential, Uniform
using Random

const N_TELLERS = 2
const MEAN_INTERARRIVAL_TIME = 4.0
const MIN_SERVICE_TIME = 2.0
const MAX_SERVICE_TIME = 10.0

tellers = nothing

"Process to generate n customers arriving into the system."
@process generator(n::Integer) begin
    dist = Exponential(MEAN_INTERARRIVAL_TIME)
    for i = 1:n
        @schedule now customer(i)
        # We don't want to wait after the last customer, which could extend the
        # simulation time past the last customer leaving.
        if i < n
            wait(rand(dist))
        end
    end
end

"The ith customer into the system."
@process customer(i::Integer) begin
    dist = Uniform(MIN_SERVICE_TIME, MAX_SERVICE_TIME)
    @with_resource tellers begin
        work(rand(dist))
    end
end

"Run the simulation for n customers."
function run_simulation(n1::Integer, n2::Integer)
    @assert n1 > 0 && n2 > 0
    println("SimLynx Closed Loop Processing Example")
    println("Number of runs = $n1, number of customers = $n2")
    println("Single Queue, $N_TELLERS Teller Bank Model")
    @simulation begin
        avg_tellers = Variable{Float64}(data=:tally, history=true)
        max_tellers = Variable{Int64}(data=:tally, history=true)
        for i = 1:n1
            @simulation begin
                global tellers = Resource("tellers")
                @schedule at 0.0 generator(n2)
                start_simulation()
                avg_tellers.value = tellers.allocated.stats.mean
                max_tellers.value = tellers.allocated.stats.max
            end
        end
        print_stats(avg_tellers, title="Average Tellers Statistics")
        plot_history(avg_tellers, title="Average Tellers History")
        print_stats(max_tellers, title="Maximum Tellers Statistics")
        plot_history(max_tellers, title="Maximum Tellers History")
    end
end

run_simulation(10_000, 100)


text:
plothomingsearching(filter(kv -> haskey(last(kv).runs[1].metadata.setup, :displace_direction), trackdata), :displace_direction)

# TODO
# you have to change the nest to originalnest becuase in dispalced there is a originalnest but in transfer they might not be...
# basically, I need to make sure that these  functions work for any combination of experiments, wihch might mean that I need to rethink what should be missing and what should not, or maybe have different summary functions for diffeent experiments
# find more bad runs and say what's wrong with them

trackdata = JLSO.load(joinpath(coffeesource, "trackdata.jlso"))
badruns = [DateTime(2017,11,17,8,54), # the particle filter really fails there
           DateTime(2017,11,15,11,44) # the turningpoint is detected way too early
          ]
kill = Pair{String, Int}[]
for (experimentid, v) in data, (i, r) in enumerate(v.runs)
    if start(r.data[:track].data.video) ∈ badruns
        push!(kill, experimentid => i)
    end
end
for (experimentid,i) in kill
    deleteat!(trackdata[experimentid].runs, i)
end

for (k,v) in filter(kv -> haskey(last(kv).runs[1].metadata.setup, :transfer), trackdata)
    plothomingsearching(Dict(k => v), :transfer)
end

for (k,v) in filter(kv -> first(kv) == "nest#closed person#therese", trackdata)
    plothomingsearching(Dict(k => v), :nest_coverage)
end