text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1PodTemplate <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    template::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PodTemplateSpec } # spec name: template

    function IoK8sApiCoreV1PodTemplate(;apiVersion=nothing, kind=nothing, metadata=nothing, template=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PodTemplate, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiCoreV1PodTemplate, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiCoreV1PodTemplate, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiCoreV1PodTemplate, Symbol("template"), template)
        setfield!(o, Symbol("template"), template)
        o
    end
end # type IoK8sApiCoreV1PodTemplate

const _property_map_IoK8sApiCoreV1PodTemplate = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("template")=>Symbol("template"))
const _property_types_IoK8sApiCoreV1PodTemplate = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("template")=>"IoK8sApiCoreV1PodTemplateSpec")
Base.propertynames(::Type{ IoK8sApiCoreV1PodTemplate }) = collect(keys(_property_map_IoK8sApiCoreV1PodTemplate))
Swagger.property_type(::Type{ IoK8sApiCoreV1PodTemplate }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PodTemplate[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PodTemplate }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PodTemplate[property_name]

function check_required(o::IoK8sApiCoreV1PodTemplate)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1PodTemplate }, name::Symbol, val)
end


text:
using Zygote: @adjoint
import Flux: trainable

mutable struct State{T, S}
    L::T
    p::S
    u::S
    ww::T
    wr::T
end

State(N::Int, R::Int, B::Int) = State(
    zeros(Float32, N, N, B),
    zeros(Float32, N, B),
    zeros(Float32, N, B),
    zeros(Float32, N, 1, B),
    zeros(Float32, N, R, B)
    )

mutable struct MemoryAccess{R, T, S}
    M::T
    state::State{T, S}
    inputmaps::R
end

MemoryAccess(inputsize, N, W, R, B; init=Flux.glorot_normal) =
    MemoryAccess(zeros(Float32, N, W, B), State(N, R, B), inputmappings(inputsize, R, W;init=init))

trainable(ma::MemoryAccess) = ma.inputmaps

function (ma::MemoryAccess)(inputs)
    p_prev, u_prev, ww_prev, wr_prev = ma.state.p, ma.state.u, ma.state.ww, ma.state.wr
    inputs = split_ξ(inputs, ma.inputmaps)
    u = usage(u_prev, ww_prev, wr_prev, inputs.f)
    ww = writeweights(ma.M, inputs, u)
    ma.M = eraseandadd(ma.M, ww, inputs.e, inputs.v)
    update_state_after_write!(ma.state, ww, u)
    wr = readweights(ma.M, inputs, ma.state.L, wr_prev)
    update_state_after_read!(ma.state, wr)
    readvectors = readmem(ma.M, wr)
    readvectors
end

function readmem(M::AbstractArray{T, 3}, wr::AbstractArray{T, 3}) where T
    batched_mul(PermutedDimsArray(M, (2, 1, 3)), wr)
end

"""
    readweights(M, inputs, L, prev_wr)

# Arguments
- `M`: (N, W, B) memory
- `inputs`: Named Tuple of controller parameters
- `L`: (N, N, B) link matrix
- `prev_wr`: (N, R, B) last iterations read weights
"""
function readweights(M, inputs, L, prev_wr)
    k, β, readmode= inputs.kr, inputs.βr, inputs.readmode
    cr = contentaddress(k, M, β)
    b = backwardweight(L, prev_wr)
    f = forwardweight(L, prev_wr)
    wr = readweight(b, cr, f, readmode)
end

"""
    writeweights(M, inputs, usage)

Fuzzy write to memory. Location is based on either content similarity or row usage.

"""
function writeweights(M, inputs, usage)
    k, β, ga, gw = inputs.kw, inputs.βw, inputs.ga, inputs.gw
    cw = contentaddress(k, M, β)
    a = allocationweighting(usage)
    ww = writeweight(cw, a, gw, ga)
end

function update_state_after_write!(state, ww, usage)
    state.u = usage
    state.ww = ww
    updatelinkmatrix!(state.L, state.p, state.ww)
    state.p = precedenceweight(state.p, state.ww)
end

@adjoint update_state_after_write!(state::State, ww, usage) =
    update_state_after_write!(state, ww, usage), _ -> nothing

function update_state_after_read!(state, wr)
    state.wr = wr
end

@adjoint update_state_after_read!(state, wr) =
    update_state_after_read!(state, wr), _ -> nothing

function eraseandadd(M::AbstractArray{T, 3}, ww::AbstractArray{T, 3}, e::AbstractArray{T,2}, a::AbstractArray{T, 2}) where T
    e = reshape(e, size(e, 1), 1, size(e, 2))
    a = reshape(a, size(a, 1), 1, size(a, 2))
    erasematrix = batched_mul(ww, PermutedDimsArray(e, (2, 1, 3)))
    addmatrix = batched_mul(ww, PermutedDimsArray(a, (2, 1, 3)))
    @cast newM[i, j, k] := M[i, j, k] * (1-erasematrix[i, j, k]) + addmatrix[i, j, k]
end

ext: jl


text:
"""
    hubbard_dispersion(k)
Dispersion relation for [`HubbardMom1D`](@ref). Returns `-2cos(k)`.

See also [`continuum_dispersion`](@ref).
"""
hubbard_dispersion(k) = -2cos(k)

"""
    continuum_dispersion(k)
Dispersion relation for [`HubbardMom1D`](@ref). Returns `k^2`.

See also [`hubbard_dispersion`](@ref).
"""
continuum_dispersion(k) = k^2

"""
    HubbardMom1D(address; u=1.0, t=1.0, dispersion=hubbard_dispersion)

Implements a one-dimensional Bose Hubbard chain in momentum space.

```math
\\hat{H} =  \\sum_{k} ϵ_k n_k + \\frac{u}{M}\\sum_{kpqr} a^†_{r} a^†_{q} a_p a_k δ_{r+q,p+k}
```

# Arguments

* `address`: the starting address, defines number of particles and sites.
* `u`: the interaction parameter.
* `t`: the hopping strength.
* `dispersion`: defines ``ϵ_k =``` t*dispersion(k)`
    - [`hubbard_dispersion`](@ref): ``ϵ_k = -2t \\cos(k)``
    - [`continuum_dispersion`](@ref): ``ϵ_k = tk^2``

# See also

* [`HubbardReal1D`](@ref)
* [`ExtendedHubbardReal1D`](@ref)
"""
struct HubbardMom1D{TT,M,AD<:AbstractFockAddress,U,T} <: AbstractHamiltonian{TT}
    add::AD # default starting address, should have N particles and M modes
    ks::SVector{M,TT} # values for k
    kes::SVector{M,TT} # values for kinetic energy
end

function HubbardMom1D(
    add::Union{BoseFS,FermiFS2C};
    u=1.0, t=1.0, dispersion = hubbard_dispersion,
)
    M = num_modes(add)
    U, T = promote(float(u), float(t))
    step = 2π/M
    if isodd(M)
        start = -π*(1+1/M) + step
    else
        start = -π + step
    end
    kr = range(start; step = step, length = M)
    ks = SVector{M}(kr)
    # kes = SVector{M}(-2T*cos.(kr))
    kes = SVector{M}(T .* dispersion.(kr))
    return HubbardMom1D{typeof(U),M,typeof(add),U,T}(add, ks, kes)
end

function Base.show(io::IO, h::HubbardMom1D)
    print(io, "HubbardMom1D($(h.add); u=$(h.u), t=$(h.t))")
end

function starting_address(h::HubbardMom1D)
    return h.add
end

LOStructure(::Type{<:HubbardMom1D{<:Real}}) = IsHermitian()

Base.getproperty(h::HubbardMom1D, s::Symbol) = getproperty(h, Val(s))
Base.getproperty(h::HubbardMom1D, ::Val{:ks}) = getfield(h, :ks)
Base.getproperty(h::HubbardMom1D, ::Val{:kes}) = getfield(h, :kes)
Base.getproperty(h::HubbardMom1D, ::Val{:add}) = getfield(h, :add)
Base.getproperty(h::HubbardMom1D{<:Any,<:Any,<:Any,U}, ::Val{:u}) where {U} = U
Base.getproperty(h::HubbardMom1D{<:Any,<:Any,<:Any,<:Any,T}, ::Val{:t}) where {T} = T

ks(h::HubbardMom1D) = getfield(h, :ks)

"""
    num_singly_doubly_occupied_sites(address)

Returns the number of singly and doubly occupied sites for a bosonic bit string address.

# Example

```jldoctest
julia> Hamiltonians.num_singly_doubly_occupied_sites(BoseFS{3,3}((1, 1, 1)))
(3, 0)
julia> Hamiltonians.num_singly_doubly_occupied_sites(BoseFS{3,3}((2, 0, 1)))
(2, 1)
```
"""
function num_singly_doubly_occupied_sites(b::BoseFS)
    singlies = 0
    doublies = 0
    for (n, _, _) in occupied_modes(b)
        singlies += 1
        doublies += n > 1
    end
    return singlies, doublies
end

function num_singly_doubly_occupied_sites(onrep::AbstractArray)
    # this one is faster by about a factor of 2 if you already have the onrep
    # returns number of singly and doubly occupied sites
    singlies = 0
    doublies = 0
    for n in onrep
        singlies += n > 0
        doublies += n > 1
    end
    return singlies, doublies
end

# standard interface function
function num_offdiagonals(ham::HubbardMom1D, add::BoseFS)
    singlies, doublies = num_singly_doubly_occupied_sites(add)
    return num_offdiagonals(ham, add, singlies, doublies)
end

# 4-argument version
@inline function num_offdiagonals(ham::HubbardMom1D, add::BoseFS, singlies, doublies)
    M = num_modes(ham)
    return singlies * (singlies - 1) * (M - 2) + doublies * (M - 1)
end
@inline function num_offdiagonals(ham::HubbardMom1D, add::FermiFS2C{N1,N2}) where {N1,N2}
    M = num_modes(ham)
    return N1 * N2 * (M - 1)
end

"""
    momentum_transfer_diagonal(H, map::OccupiedModeMap)

Compute diagonal interaction energy term.

# Example

```jldoctest
julia> a = BoseFS{6,5}((1,2,3,0,0))
BoseFS{6,5}((1, 2, 3, 0, 0))

julia> H = HubbardMom1D(a);

julia> Hamiltonians.momentum_transfer_diagonal(H, OccupiedModeMap(a))
5.2
```
"""
@inline function momentum_transfer_diagonal(
    h::HubbardMom1D{<:Any,M,<:BoseFS}, map
) where {M}
    return h.u / 2M * momentum_transfer_diagonal(map)
end
@inline function momentum_transfer_diagonal(
    h::HubbardMom1D{<:Any,M,<:FermiFS2C}, map_a, map_b
) where {M}
    return h.u / 2M * momentum_transfer_diagonal(map_a, map_b)
end

@inline function diagonal_element(h::HubbardMom1D, add::BoseFS)
    map = OccupiedModeMap(add)
    return dot(h.kes, map) + momentum_transfer_diagonal(h, map)
end
@inline function diagonal_element(h::HubbardMom1D, add::FermiFS2C)
    map_a = OccupiedModeMap(add.components[1])
    map_b = OccupiedModeMap(add.components[2])
    return dot(h.kes, map_a) + dot(h.kes, map_b) +
        momentum_transfer_diagonal(h, map_a, map_b)
end

@inline function get_offdiagonal(
    ham::HubbardMom1D{<:Any,M,A}, add::A, chosen, map=OccupiedModeMap(add)
) where {M,A<:BoseFS}
    add, onproduct = momentum_transfer_excitation(add, chosen, map)
    return add, ham.u/(2*M)*onproduct
end
@inline function get_offdiagonal(
    ham::HubbardMom1D{<:Any,M,A}, add::A, chosen,
    map_a=OccupiedModeMap(add.components[1]), map_b=OccupiedModeMap(add.components[2])
) where {M,A<:FermiFS2C}
    add_a, add_b = add.components
    new_add_a, new_add_b, onproduct = momentum_transfer_excitation(
        add_a, add_b, chosen, map_a, map_b
    )
    return CompositeFS(new_add_a, new_add_b), ham.u/M * onproduct
end

###
### offdiagonals
###
"""
    OffdiagonalsBoseMom1D

Specialized [`AbstractOffdiagonals`](@ref) that keeps track of singly and doubly occupied
sites in current address.
"""
struct OffdiagonalsBoseMom1D{
    A<:BoseFS,T,H<:AbstractHamiltonian{T},O<:OccupiedModeMap
} <: AbstractOffdiagonals{A,T}
    hamiltonian::H
    address::A
    length::Int
    map::O
end

function offdiagonals(h::HubbardMom1D, a::BoseFS)
    map = OccupiedModeMap(a)
    singlies = length(map)
    doublies = count(i -> i.occnum ≥ 2, map)
    num = num_offdiagonals(h, a, singlies, doublies)
    return OffdiagonalsBoseMom1D(h, a, num, map)
end

function Base.getindex(s::OffdiagonalsBoseMom1D{A,T}, i)::Tuple{A,T} where {A,T}
    @boundscheck begin
        1 ≤ i ≤ s.length || throw(BoundsError(s, i))
    end
    new_address, matrix_element = get_offdiagonal(s.hamiltonian, s.address, i, s.map)
    return (new_address, matrix_element)
end

Base.size(s::OffdiagonalsBoseMom1D) = (s.length,)

struct OffdiagonalsFermiMom1D2C{
    F<:FermiFS2C,T,H<:AbstractHamiltonian{T},O1,O2
} <: AbstractOffdiagonals{F,T}
    hamiltonian::H
    address::F
    length::Int
    map_a::O1
    map_b::O2
end

function offdiagonals(h::HubbardMom1D, f::FermiFS2C)
    comp_a, comp_b = f.components
    map_a = OccupiedModeMap(comp_a)
    map_b = OccupiedModeMap(comp_b)
    num = num_offdiagonals(h, f)
    return OffdiagonalsFermiMom1D2C(h, f, num, map_a, map_b)
end

Base.size(s::OffdiagonalsFermiMom1D2C) = (s.length,)

function Base.getindex(s::OffdiagonalsFermiMom1D2C{A,T}, i)::Tuple{A,T} where {A,T}
    @boundscheck begin
        i ≤ i ≤ s.length || throw(BoundsError(s, i))
    end
    new_address, matrix_element = get_offdiagonal(
        s.hamiltonian, s.address, i, s.map_a, s.map_b
    )
    return (new_address, matrix_element)
end

###
### momentum
###
struct MomentumMom1D{T,H<:AbstractHamiltonian{T}} <: AbstractHamiltonian{T}
    ham::H
end
LOStructure(::Type{MomentumMom1D{H,T}}) where {H,T <: Real} = IsHermitian()
num_offdiagonals(ham::MomentumMom1D, add) = 0
diagonal_element(mom::MomentumMom1D, add) = mod1(onr(add)⋅ks(mom.ham) + π, 2π) - π # fold into (-π, π]

momentum(ham::HubbardMom1D) = MomentumMom1D(ham)


text:
module ClutteredEnvPathOpt
using Pipe
import LightGraphs

export LabeledGraph, find_feg_separator_lt, find_separator_lt, find_separator_fcs, find_separator_fcs_best, pp_expell, find_biclique_cover, find_feg_separator_lt_best
include("types.jl")
include("separator.jl")
include("biclique.jl")
include("obstacles.jl")
end # module


text:
using CauseMap

para, didi = load_example_data("ParaDidi")

E_vals     = 2:6  # range to test of system dimensionality
tau_s_vals = 1:1 # range for lag length for manifold reconstruction
tau_p_vals = 0:5  # range to test for time lag of causal effect

# makesingleplot(para, didi ,
#                                         4, 1, 0,
#                                        "x", "y")

makeoptimizationplots(repeat(para, outer=[1]),
                                    repeat(didi, outer=[1]),
                                        2:2, 1:1, 1:1,
                                        "Para.", "Didi."; plot=false
                                    )

maxTSlen = 6
TSlens = 1:maxTSlen
outarray = fill(NaN, (maxTSlen, 2))
for TSlen in 2:2
    tic()
    makeoptimizationplots(repeat(para, outer=[TSlen]),
                                        repeat(didi, outer=[TSlen]),
                                            E_vals, tau_s_vals, tau_p_vals,
                                            "Para.", "Didi."; plot=true
                                        )
    time = toc()
    outarray[TSlen, :] = [TSlen * length(para), time]
end
writedlm("TSresults.txt", outarray, '\t')


text:
import Dates
import HTTP
import JSON

function compare_http_date_header(header_value::String, timestamp_request_completed::Dates.DateTime) :: Nothing
    header_value_timestamp::Dates.DateTime = Dates.DateTime(split(header_value, " UTC")[1], "e, d u Y H:M:S")
    @test header_value_timestamp <= timestamp_request_completed
    nothing
end

function compare_http_header(headers::Array, key::String, value::String) :: Nothing
    @test header_get_value(headers::Array, key::String)==value
    nothing
end

function header_get_value(headers::Array, key::String) :: String
    for item in headers
        if item[1]==key
            return item[2]
        end
    end
end

function parse_and_test_request(r::HTTP.Messages.Response, status::Int64, headers::Vector{Pair{String,String}}, content_length::Int64, is_json::Bool, is_get::Bool, body::Union{Dict, String}="")
    timestamp_request_completed::Dates.DateTime = Dates.now(Dates.UTC)

    @test r.status==status

    for (key, value) in headers
        compare_http_header(r.headers, key, value)
    end
    compare_http_date_header(header_get_value(r.headers, "Date"), timestamp_request_completed)
    # TODO: fix Windows HTML file longer due to line final char
    if !Sys.iswindows()
        compare_http_header(r.headers, "Content-Length", string(content_length))
    end

    if is_json
        header_get_value(r.headers, "Content-Type") == "application/json"
        @test JSON.json(body) == String(r.body)
    elseif is_get
        header_get_value(r.headers, "Content-Type") == "text/html; charset=UTF-8"
        @test occursin("Content ...", String(r.body))
    else
        length(r.body) == 0
    end

    nothing
end


text:
# Modules

handle("editor.julia.module.update") do editor, data
  file = data["path"] == nothing ? "This file" : basename(data["path"])
  mod = data["path"] == nothing ? "Main" : Jewel.filemodule(data["path"])
  mod == "" && (mod = "Main")
  if Jewel.getthing(mod) == nothing
    notify("$file's module, `$mod`, isn't loaded yet.", class = "error")
    mod = "Main"
  end
  raise(editor, "julia.set-module", mod)
  raise(global_client, "julia.set-modules", @d(:modules => [string(m) for m in Jewel.allchildren(Main)]))
end

# Browser

objs(mod) =
  @>> mod begin
    (m->names(m, true))
    map(n-> c(n, Jewel.getthing(mod, [n])))
    filter(p->!isa(p[2], Module) && p[2] ≠ nothing && !ismatch(r"#", string(p[1])))
    map(p-> c(p[1], applydisplayinline(p[2])))
    map(p-> c(p[1], stringmime(bestmime(p[2]), p[2])))
    xs->sort(xs, by = first)
  end

handle("browser.get-objects") do browser, data
  mod = data["module"] == nothing ? Main : Jewel.getthing(data["module"])
  raise(browser, "update", @d(:objs => objs(mod)))
end


text:
# Tensor type
# Holds all the needed information for general calculations
# and info needed for backprop
mutable struct Tensor{T,N} <: AbstractArray{T,N}

    data::Array{T,N}

    req_grad::Bool
    grad::Array{Float32, N}

    # Attributes to keep track of computation
    # Operation by which the Tensor was created
    op::Type
    # Other elements (Tensors, Arrays) whose computation created this Tensor
    prev::Tuple

    function Tensor(data::Array{T,N},
                    req_grad::Bool=true; # Tensors should keep track of computation by default
                    grad::Array{Float64,N}=zeros(size(data)),
                    op=Type{},
                    prev=()) where {T,N}
        new{T,N}(data, req_grad, grad, op, prev)
    end
end

# Tensor should behave like an AbstractArray with keeping track of computations and the ability to backpropagate
# implement AbstractArray interface
import Base: size, getindex, setindex!

Base.IndexStyle(::Tensor) = IndexLinear()

Base.size(tensor::Tensor) = size(tensor.data)

@inline function Base.getindex(tensor::Tensor, i::Int)
    getindex(tensor.data, i)
end

@inline function Base.setindex!(tensor::Tensor, v, i::Int)
    setindex!(tensor.data, v, i)
end

include("ops.jl")
include("grad.jl")
include("broadcast.jl")

function backward!(t::Tensor, result::Bool=true)
    # End condition
    # No previous nodes, backward computation ends here
    if t.prev == ()
        return
    end

    # Implicit grad creation
    if result
        t.grad = ones(size(t))
    end

    calculategrad!(t.op, t.prev[1], t.prev[2], t)

    # traverse computation graph in reverse order
    backward!(t.prev[1], false)
    backward!(t.prev[2], false)
end

# Gradient of this node should not/cannot be determined
function backward!(t::Any, result::Bool) end

max_stars_repo_name: leonvol/magrad

max_issues_repo_name: leonvol/magrad

max_forks_repo_name: leonvol/magrad


text:
__precompile__()

module GuidingCenterOrbits

#import standard libraries
using Printf
using LinearAlgebra
using Statistics

using Equilibrium
using OrdinaryDiffEq
using ForwardDiff
using Optim
using Base.Iterators
using Interpolations
using StaticArrays
using Contour
#using Plots

const e0 = 1.60217733e-19 # Coulombs / Joules
const mu0 = 4*pi*1e-7 # N/A^2
const c0 = 2.99792458e8 # m/s

const mass_u = 1.6605402e-27 # kg
const e_amu = 5.48579909070e-4 # amu
const H1_amu = 1.007276466879 # amu
const H2_amu = 2.0141017778 # amu
const H3_amu = 3.01550071632 # amu
const He3_amu = 3.01602931914 # amu
const He4_amu = 4.00150617913 # amu
const B5_amu = 10.81 # amu
const C6_amu = 12.011 # amu

include("particles.jl")
export AbstractParticle, Particle
export Electron, Proton, Deuteron, Triton, Alpha
export GCParticle, EPRParticle
export GCElectron, GCProton, GCDeuteron, GCTriton, GCAlpha

include("coordinates.jl")
export AbstractOrbitCoordinate, EPRCoordinate, HamiltonianCoordinate

include("callbacks.jl")
export standard_callback

include("orbit.jl")
export Orbit, OrbitPath, integrate, get_orbit, gc_velocity

include("fullorbit.jl")
export FullOrbitPath, get_full_orbit

include("projection.jl")
export OrbitProjection, orbit_projection, in_orbit

include("utils.jl")
export get_pitch, hits_wall_path, hits_wall, get_kinetic_energy
export lorentz_factor, cyclotron_frequency, cyclotron_period, normalize
export perpendicular_vectors, gyro_step, larmor_radius, gcde_check

using ForwardDiff
using ForwardDiff: Dual, partials, value

include("jacobian.jl")
export get_jacobian, transform

end


text:
const TA0_DN = UInt8[
	0x30, 0x1C, 0x31, 0x0B, 0x30, 0x09, 0x06, 0x03, 0x55, 0x04, 0x06, 0x13,
	0x02, 0x43, 0x41, 0x31, 0x0D, 0x30, 0x0B, 0x06, 0x03, 0x55, 0x04, 0x03,
	0x13, 0x04, 0x52, 0x6F, 0x6F, 0x74
]

const TA0_RSA_N = UInt8[
	0xB6, 0xD9, 0x34, 0xD4, 0x50, 0xFD, 0xB3, 0xAF, 0x7A, 0x73, 0xF1, 0xCE,
	0x38, 0xBF, 0x5D, 0x6F, 0x45, 0xE1, 0xFD, 0x4E, 0xB1, 0x98, 0xC6, 0x60,
	0x83, 0x26, 0xD2, 0x17, 0xD1, 0xC5, 0xB7, 0x9A, 0xA3, 0xC1, 0xDE, 0x63,
	0x39, 0x97, 0x9C, 0xF0, 0x5E, 0x5C, 0xC8, 0x1C, 0x17, 0xB9, 0x88, 0x19,
	0x6D, 0xF0, 0xB6, 0x2E, 0x30, 0x50, 0xA1, 0x54, 0x6E, 0x93, 0xC0, 0xDB,
	0xCF, 0x30, 0xCB, 0x9F, 0x1E, 0x27, 0x79, 0xF1, 0xC3, 0x99, 0x52, 0x35,
	0xAA, 0x3D, 0xB6, 0xDF, 0xB0, 0xAD, 0x7C, 0xCB, 0x49, 0xCD, 0xC0, 0xED,
	0xE7, 0x66, 0x10, 0x2A, 0xE9, 0xCE, 0x28, 0x1F, 0x21, 0x50, 0xFA, 0x77,
	0x4C, 0x2D, 0xDA, 0xEF, 0x3C, 0x58, 0xEB, 0x4E, 0xBF, 0xCE, 0xE9, 0xFB,
	0x1A, 0xDA, 0xA3, 0x83, 0xA3, 0xCD, 0xA3, 0xCA, 0x93, 0x80, 0xDC, 0xDA,
	0xF3, 0x17, 0xCC, 0x7A, 0xAB, 0x33, 0x80, 0x9C, 0xB2, 0xD4, 0x7F, 0x46,
	0x3F, 0xC5, 0x3C, 0xDC, 0x61, 0x94, 0xB7, 0x27, 0x29, 0x6E, 0x2A, 0xBC,
	0x5B, 0x09, 0x36, 0xD4, 0xC6, 0x3B, 0x0D, 0xEB, 0xBE, 0xCE, 0xDB, 0x1D,
	0x1C, 0xBC, 0x10, 0x6A, 0x71, 0x71, 0xB3, 0xF2, 0xCA, 0x28, 0x9A, 0x77,
	0xF2, 0x8A, 0xEC, 0x42, 0xEF, 0xB1, 0x4A, 0x8E, 0xE2, 0xF2, 0x1A, 0x32,
	0x2A, 0xCD, 0xC0, 0xA6, 0x46, 0x2C, 0x9A, 0xC2, 0x85, 0x37, 0x91, 0x7F,
	0x46, 0xA1, 0x93, 0x81, 0xA1, 0x74, 0x66, 0xDF, 0xBA, 0xB3, 0x39, 0x20,
	0x91, 0x93, 0xFA, 0x1D, 0xA1, 0xA8, 0x85, 0xE7, 0xE4, 0xF9, 0x07, 0xF6,
	0x10, 0xF6, 0xA8, 0x27, 0x01, 0xB6, 0x7F, 0x12, 0xC3, 0x40, 0xC3, 0xC9,
	0xE2, 0xB0, 0xAB, 0x49, 0x18, 0x3A, 0x64, 0xB6, 0x59, 0xB7, 0x95, 0xB5,
	0x96, 0x36, 0xDF, 0x22, 0x69, 0xAA, 0x72, 0x6A, 0x54, 0x4E, 0x27, 0x29,
	0xA3, 0x0E, 0x97, 0x15
]

const TA0_RSA_E = UInt8[
	0x01, 0x00, 0x01
]

const TA1_DN = UInt8[
	0x30, 0x1C, 0x31, 0x0B, 0x30, 0x09, 0x06, 0x03, 0x55, 0x04, 0x06, 0x13,
	0x02, 0x43, 0x41, 0x31, 0x0D, 0x30, 0x0B, 0x06, 0x03, 0x55, 0x04, 0x03,
	0x13, 0x04, 0x52, 0x6F, 0x6F, 0x74
]

const TA1_EC_Q = UInt8[
	0x04, 0x71, 0x74, 0xBA, 0xAB, 0xB9, 0x30, 0x2E, 0x81, 0xD5, 0xE5, 0x57,
	0xF9, 0xF3, 0x20, 0x68, 0x0C, 0x9C, 0xF9, 0x64, 0xDB, 0xB4, 0x20, 0x0D,
	0x6D, 0xEA, 0x40, 0xD0, 0x4A, 0x6E, 0x42, 0xFD, 0xB6, 0x9A, 0x68, 0x25,
	0x44, 0xF6, 0xDF, 0x7B, 0xC4, 0xFC, 0xDE, 0xDD, 0x7B, 0xBB, 0xC5, 0xDB,
	0x7C, 0x76, 0x3F, 0x41, 0x66, 0x40, 0x6E, 0xDB, 0xA7, 0x87, 0xC2, 0xE5,
	0xD8, 0xC5, 0xF3, 0x7F, 0x8D
]

struct br_x500_name
    data::Vector{UInt8}
    len::Csize_t
end

abstract type PublicKey end

struct br_rsa_public_key <: PublicKey
    n::Vector{UInt8}
    nlen::Csize_t
    e::Vector{UInt8}
    elen::Csize_t
end

struct br_ec_public_key <: PublicKey
    curve::Cint
    q::Vector{UInt8}
    qlen::Csize_t
end

struct br_x509_pkey{T <: PublicKey}
    key_type::UInt8
    key::T
end

struct br_x509_trust_anchor{T <: PublicKey}
    dn::br_x500_name
    flags::Cuint
    pkey::br_x509_pkey{T}
end

const BR_X509_TA_CA = Cuint(0x0001)
const BR_KEYTYPE_RSA = 0x01
const BR_KEYTYPE_EC = 0x02
const BR_EC_secp256r1 = Cint(23)

const TAs = br_x509_trust_anchor[
	br_x509_trust_anchor(
		br_x500_name(TA0_DN, sizeof(TA0_DN)),
		BR_X509_TA_CA,
		br_x509_pkey(
			BR_KEYTYPE_RSA,
			br_rsa_public_key(
				TA0_RSA_N, sizeof(TA0_RSA_N),
				TA0_RSA_E, sizeof(TA0_RSA_E),
            )
		)
	),
	br_x509_trust_anchor(
		br_x500_name(TA1_DN, sizeof(TA1_DN)),
		BR_X509_TA_CA,
		br_x509_pkey(
			BR_KEYTYPE_EC,
			br_ec_public_key(
				BR_EC_secp256r1,
                TA1_EC_Q, sizeof(TA1_EC_Q)
            )
        )
    )
]

const TAs_NUM = 2


text:
@testset "Video encode/decode accuracy (read, encode, read, compare)" begin
    file = joinpath(videodir, "annie_oakley.ogg")
    imgstack_rgb = VideoIO.openvideo(collect, file)
    imgstack_gray = convert.(Array{Gray{N0f8}}, imgstack_rgb)
    @testset "Lossless Grayscale encoding" begin
        encoder_options = (color_range=2, crf=0, preset="medium")
        VideoIO.save(tempvidpath, imgstack_gray, codec_name = "libx264", encoder_options = encoder_options)
        imgstack_gray_copy = VideoIO.openvideo(collect, tempvidpath, target_format = VideoIO.AV_PIX_FMT_GRAY8)
        @test eltype(eltype(imgstack_gray)) == eltype(eltype(imgstack_gray_copy))
        @test length(imgstack_gray) == length(imgstack_gray_copy)
        @test size(imgstack_gray[1]) == size(imgstack_gray_copy[1])
        @test !any(.!(imgstack_gray .== imgstack_gray_copy))
    end

    @testset "Lossless RGB encoding" begin
        encoder_options = (color_range=2, crf=0, preset="medium")
        codec_name="libx264rgb"
        VideoIO.save(tempvidpath, imgstack_rgb, codec_name = codec_name, encoder_options = encoder_options)
        imgstack_rgb_copy = VideoIO.openvideo(collect, tempvidpath)
        @test eltype(imgstack_rgb) == eltype(imgstack_rgb_copy)
        @test length(imgstack_rgb) == length(imgstack_rgb_copy)
        @test size(imgstack_rgb[1]) == size(imgstack_rgb_copy[1])
        @test !any(.!(imgstack_rgb .== imgstack_rgb_copy))
    end

    @testset "UInt8 accuracy during read & lossless encode" begin
        # Test that reading truth video has one of each UInt8 value pixels (16x16 frames = 256 pixels)
        raw_img = VideoIO.openvideo(read,
                                    joinpath(testdir, "references", "precisiontest_gray_truth.mp4"),
                                    target_format = VideoIO.AV_PIX_FMT_GRAY8)
        frame_truth = collect(rawview(channelview(raw_img)))
        h_truth = fit(Histogram, frame_truth[:], 0:256)
        @test h_truth.weights == fill(1,256) #Test that reading is precise

        # Test that encoding new test video has one of each UInt8 value pixels (16x16 frames = 256 pixels)
        img = Array{UInt8}(undef,16,16)
        for i in 1:256
            img[i] = UInt8(i-1)
        end
        imgstack = []
        for i=1:24
            push!(imgstack,img)
        end

        encoder_options = (color_range=2, crf=0, preset="medium")
        VideoIO.save(tempvidpath, imgstack, encoder_options = encoder_options)
        f = VideoIO.openvideo(tempvidpath, target_format = VideoIO.AV_PIX_FMT_GRAY8)
        try
            frame_test = collect(rawview(channelview(read(f))))
            h_test = fit(Histogram, frame_test[:], 0:256)
            @test h_test.weights == fill(1,256) #Test that encoding is precise (if above passes)
            @test VideoIO.counttotalframes(f) == 24
        finally
            close(f)
        end
    end

    @testset "Correct frame order when reading & encoding" begin
        @testset "Frame order when reading ground truth video" begin
            # Test that reading a video with frame-incremental pixel values is read in in-order
            f = VideoIO.openvideo(joinpath(testdir,"references","ordertest_gray_truth.mp4"),target_format=VideoIO.AV_PIX_FMT_GRAY8)
            frame_ids_truth = []
            while !eof(f)
                img = collect(rawview(channelview(read(f))))
                push!(frame_ids_truth,img[1,1])
            end
            @test frame_ids_truth == collect(0:255) #Test that reading is in correct frame order
            @test VideoIO.counttotalframes(f) == 256
        end
        @testset "Frame order when encoding, then reading video" begin
            # Test that writing and reading a video with frame-incremental pixel values is read in in-order
            imgstack = []
            img = Array{UInt8}(undef,16,16)
            for i in 0:255
                push!(imgstack,fill(UInt8(i),(16,16)))
            end

            encoder_options = (color_range=2, crf=0, preset="medium")
            VideoIO.save(tempvidpath, imgstack, encoder_options = encoder_options)
            f = VideoIO.openvideo(tempvidpath, target_format = VideoIO.AV_PIX_FMT_GRAY8)
            try
                frame_ids_test = []
                while !eof(f)
                    img = collect(rawview(channelview(read(f))))
                    push!(frame_ids_test,img[1,1])
                end
                @test frame_ids_test == collect(0:255) #Test that reading is in correct frame order
                @test VideoIO.counttotalframes(f) == 256
            finally
                close(f)
            end
        end
    end
end


text:
function sphereFunc(x::AbstractArray)
    return sum(x.^2)
end

function rosenbrockfunc(x::AbstractVector)
    sum = 0.0
    for i in 1:length(x) - 1
        sum += 100*(x[i + 1] - x[i]^2)^2 + (x[i] - 1)^2
    end
    return sum
end

function rastriginfunc(x::AbstractVector)
    d = length(x)

    sum = 0.0
    for i in 1:d
        sum += (x[i]^2 - 10*cos(2*π*x[i]))
    end
    return 10*d + sum
end

function ackleyfunc(x::AbstractVector)
    a = 20
    b = 0.2
    c = 2*π
    d = length(x)

    sum1 = 0.0
    sum2 = 0.0
    for i in 1:d
        sum1 += x[i]^2
        sum2 += cos(c*x[i])
    end
    return -a*exp(-b*sqrt(sum1/d))  - exp(sum2/d) + a + exp(1)
end


text:
# AbstractTable interface

"""
    `fields(tbl::Table)`

Obtain an ordered list of fields (column names) of `tbl`.

Notes: `fields(tbl)` is dual to `index(tbl)` in the sense that

fields(tbl)[(index(tbl)[fld]] == fld
index(tbl)[(fields(tbl)[i]] == i
"""
AbstractTables.fields(tbl::Table) = tbl.fields

AbstractTables.eltypes(tbl::Table) = map(eltype, columns(tbl))
AbstractTables.eltypes(tbl::Table, fields::Symbol...) =
    map(eltype, [ tbl[field] for field in fields ])

"""
    `nrow(tbl::Table)`

Obtain the number of rows contained in `tbl`.
"""
AbstractTables.nrow(tbl::Table) = ncol(tbl) > 0 ? length(tbl.columns[1]) : 0

# Column indexable interface

"""
    `columns(tbl::Table)`

Obtain the columns of `tbl`.
"""
AbstractTables.columns(tbl::Table) = tbl.columns

"""
    `index(tbl::Table)`

Obtain the index of `tbl`.
"""
AbstractTables.index(tbl::Table) = tbl.index

# Other

"""
    `copy(tbl)`

Return a copy of `tbl`.

Notes: applies `copy` to each column and inserts the copies into an `empty`
`Table`.
"""
function Base.copy(tbl::Table)
    new_tbl = Table()
    for (fld, col) in eachcol(tbl)
        new_tbl[fld] = copy(col)
    end
    return new_tbl
end

function Base.isequal(tbl1::Table, tbl2::Table)
    isequal(ncol(tbl1), ncol(tbl2)) || return false
    # isequal(tbl1.allowsnulls, tbl2.allowsnulls) || return false
    isequal(tbl1.hasnulls, tbl2.hasnulls) || return false
    for ((fld1, col1), (fld2, col2)) in zip(eachcol(tbl1), eachcol(tbl2))
        isequal(fld1, fld2) || return false
        isequal(col1, col2) || return false
    end
    return true
end

function Base.hash(tbl::Table)
    # h = hash(tbl.allowsnulls) + 1
    h = hash(tbl.hasnulls) + 1
    for (i, (fld, col)) in enumerate(eachcol(tbl))
        h = hash(i, h)
        h = hash(fld, h)
        h = hash(tbl[fld], h)
    end
    return @compat UInt(h)
end


text:
using mdfReader
using Test

@testset "mdfReader.jl" begin
    @test mdfReader.open("testfiles/test.mdf")
end


text:
# Setup:

isbinder() = "jovyan" in split(pwd(), "/")

const REPO = "https://github.com/ablaom/MachineLearningInJulia2020"
using Pkg
uuid = Pkg.TOML.parsefile(joinpath(DIR,"Project.toml"))["uuid"]
uuid == "4764ce03-6504-4302-ab9a-b32cdba420f3" ||
    error("It appears this file is not in the same directory as other "*
          "files it needs (and in particular, the intended "*
          "Project.toml file). "*
          "A complete tutorial can be obtained by cloning the"*
          "MachineLearningInJulia repository from $REPO. ")
if !isbinder()
    Pkg.activate(DIR)
    Pkg.instantiate()
    using CategoricalArrays
    import MLJLinearModels
    import DataFrames
    import CSV
    import DecisionTree
    using MLJ
    import MLJLinearModels
    import MultivariateStats
    import MLJFlux
    import Plots
else
    @info "Skipping package instantiation as binder notebook. "
end
@info "Done loading"


text:
# Add our own pairwise function to be able to apply it on vectors

pairwise(d::PreMetric, X::AbstractVector, Y::AbstractVector) = broadcast(d, X, permutedims(Y))

pairwise(d::PreMetric, X::AbstractVector) = pairwise(d, X, X)

function pairwise!(
    out::AbstractMatrix,
    d::PreMetric,
    X::AbstractVector,
    Y::AbstractVector,
)
    broadcast!(d, out, X, Y')
end

pairwise!(out::AbstractMatrix, d::PreMetric, X::AbstractVector) = pairwise!(out, d, X, X)

function pairwise(d::PreMetric, x::AbstractVector{<:Real})
    return Distances.pairwise(d, reshape(x, :, 1); dims = 1)
end

function pairwise(
    d::PreMetric,
    x::AbstractVector{<:Real},
    y::AbstractVector{<:Real},
)
    return Distances.pairwise(d, reshape(x, :, 1), reshape(y, :, 1); dims = 1)
end

function pairwise!(out::AbstractMatrix, d::PreMetric, x::AbstractVector{<:Real})
    return Distances.pairwise!(out, d, reshape(x, :, 1); dims = 1)
end

function pairwise!(
    out::AbstractMatrix,
    d::PreMetric,
    x::AbstractVector{<:Real},
    y::AbstractVector{<:Real},
)
    return Distances.pairwise!(out, d, reshape(x, :, 1), reshape(y, :, 1); dims=1)
end


text:
module BinaryTreeTassellation

using Luxor
using Random

export Tile, TileLeaf, TileNode
export base, height, rescale, rescale!, size
export paint_tree

include("Tile.jl")
include("palettes.jl")
include("draw.jl")
include("tree.jl")

end


text:
using DataFrames
using CSV

data = readdlm("55.4_roykenes_ptq.txt")

date = [DateTime(data[i,1], data[i,2], data[i,3], data[i,4]) for i in 1:size(data, 1)]

prec = data[:, 5:14]
tair = data[:, 15:24]
qobs = data[:, end]

qobs = qobs * (3 * 3600 * 1000) / (50.09 * 1e6)

qobs = round(qobs, 2)

df_prec = DataFrame(hcat(date, prec))
df_tair = DataFrame(hcat(date, tair))
df_qobs = DataFrame(hcat(date, qobs))

CSV.write("Prec.txt", df_prec, delim = ';', header = false, dateformat = "yyyy-mm-dd HH:MM")
CSV.write("Tair.txt", df_tair, delim = ';', header = false, dateformat = "yyyy-mm-dd HH:MM")
CSV.write("Q_obs.txt", df_qobs, delim = ';', header = false, dateformat = "yyyy-mm-dd HH:MM")

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
mutable struct CircularBuffer{T}

    function CircularBuffer{T}(capacity::Integer) where {T}

    end
end

function Base.push!(cb::CircularBuffer, item; overwrite::Bool=false)

end

function Base.popfirst!(cb::CircularBuffer)

end

function Base.empty!(cb::CircularBuffer)

end


text:
using LinearAlgebra, Plots, Statistics, TaylorSeries, Printf
# if !(pwd() in LOAD_PATH) push!(LOAD_PATH, pwd()) end
# using TaylorBigF

path1=realpath(dirname(@__FILE__)*"/..")
include(string(path1,"/src/TaylorBigF.jl"))

# 3.6 Solution of Ordinary Differential Equations

p=50;setprecision(p); setrounding(BigFloat, RoundUp);
L=BigFloat("1"); n=BigInt(100); iN=convert(Int64,n);
dx=BigFloat("1")*L/(n-1); x=0:dx:(n-1)*dx;
I=zeros(BigFloat, n, n);for i=1:n I[i,i]=BigFloat("1") end

V=TaylorBigF.calc_Vandermonde(x,n,n+BigInt(4))
# Vi=calc_Vandermonde(xi,n-1,n+BigInt(4))
dV=TaylorBigF.calc_Vandermonde_dot(x,n,n+BigInt(4))
ddV=TaylorBigF.calc_Vandermonde_dot_dot(x,n,n+BigInt(4))
dddV=TaylorBigF.calc_Vandermonde_dot_dot_dot(x,n,n+BigInt(4))
ddddV=TaylorBigF.calc_Vandermonde_dot_dot_dot_dot(x,n,n+BigInt(4))

xi=(x.+dx/BigFloat("2"))[1:end-1];
Vi=TaylorBigF.calc_Vandermonde(xi,n-1,n+BigInt(4));

q=zeros(BigFloat,n)
# q.-=1e2sin.(exp.(3x))
# q.+=convert(Array{BigFloat},(rand(n).-1/2).*10)
# plot(q)
Eq1=V[1,:]
q1=BigFloat("0")
Eq2=dV[1,:]
q2=BigFloat("0")
Eq3=V[end,:]
q3=BigFloat("1")
Eq4=dV[end,:]
q4=BigFloat("0")

VV=[ddddV;Eq1';Eq2';Eq3';Eq4']
q_all=[q;q1;q2;q3;q4]
a=VV[5:end,:]\q_all[5:end]
# a=convert(Array{Float64},VV[5:end,:])\convert(Array{Float64},q_all[5:end])
# plot(a,xlabel="n",ylabel="a",color=:black,width=2,
#         label=string("p=",p),linestyle=:dash)
# plot!(a,xlabel="n",ylabel="a",color=:black,width=2,
#         label=string("p=",p),linestyle=:dashdotdot)
# savefig(string("ode_beam.pdf"))

# r_root=Array{Float64}(undef,0); for i=1:Int64(n-1) push!(r_root,abs(a[i])^(1/(n-1))) end; r_rootF=convert(Array{Float64},r_root);
# plot(r_rootF)
sol1=V*a
scatter(convert(Array{Float64},x),sol1)

exact1=-2x.^3+3x.^2
plot!(convert(Array{Float64},x),exact1)

error1=exact1-sol1
# plot(error1)
@sprintf("error max=%1.3e",maximum(abs.(error1)))

exacti1=-2xi.^3+3xi.^2
soli1=Vi*a
errori1=exacti1-soli1
@sprintf("error i max=%1.3e",maximum(abs.(errori1)))
plot!(convert(Array{Float64},xi),soli1)

M=ddV*a
plot(convert(Array{Float64},x),M)

precision(BigFloat)
m1=string(M[1])


text:
#"""
#For now, we have to define our own indexing functions for performance.
#because of this issue:
#https://github.com/JuliaLang/julia/issues/21957
#"""
#@inline function Base.getindex(S::SharedArray, i::Real)
#	getindex(S.s, i)
#end
#@inline function Base.setindex!(S::SharedArray, x, i::Real)
#	setindex!(S.s, x, i)
#end
#@inline function Base.getindex(S::SharedArray, i::Real...)
#	getindex(S.s, i...)
#end
#@inline function Base.setindex!(S::SharedArray, x, i::Real...)
#	setindex!(S.s, x, i...)
#end

function setindex_sharedarray_worker(ba::BigArray{D,T},
                                sharedBuffer::OffsetArray{T,N,SharedArray{T,N}},
                                chunkGlobalRange::CartesianIndices{N},
                                rangeInBuffer::CartesianIndices{N}) where {D,T,N}
    C = get_encoding(ba)
    try
        @inbounds block = sharedBuffer[rangeInBuffer]
        @inbounds ba.kvStore[ cartesian_range2string(chunkGlobalRange) ] = encode( block, C)
    catch err
        println("catch an error while saving in BigArray: $err with type $(typeof(err))")
        rethrow()
    end
    nothing
end

"""
    put array in RAM to a BigArray
this version uses channel to control the number of asynchronized request
"""
function setindex_sharedarray!( ba::BigArray{D,T}, buf::Array{T,N},
                       idxes::Union{UnitRange, Int, Colon} ... ) where {D,T,N}
    idxes = colon2unit_range(buf, idxes)
    sharedBuffer = OffsetArray(SharedArray(buf), idxes...)
    # check alignment
    @assert all(map((x,y,z)->mod(x.start - 1 - y, z),
                    idxes, ba.offset.I, ba.chunkSize).==0)
                    "the start of index should align with BigArray chunk size"
    t1 = time()
    baIter = ChunkIterator(idxes, ba.chunkSize; offset=ba.offset)
    @sync @distributed for (blockID, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer) in baIter
        chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer =
            adjust_volume_boundary(ba, chunkGlobalRange, globalRange,
                                   rangeInChunk, rangeInBuffer)
        setindex_sharedarray_worker(ba, sharedBuffer, chunkGlobalRange, rangeInBuffer)
        #@async remote_do(setindex_remote_worker, WORKER_POOL,
        #                           ba, sharedBuffer, chunkGlobalRange, rangeInBuffer)
    end
    elapsed = time() - t1 # sec
    println("saving speed: $(sizeof(sharedBuffer)/1024/1024/elapsed) MB/s")
end

function getindex_sharedarray_worker!(ba::BigArray{D,T}, jobs::RemoteChannel,
                                sharedBuffer::OffsetArray{T,N,SharedArray{T,N}}) where {D,T,N}
    baRange = CartesianIndices(ba)
    blockId, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer = take!(jobs)
    if any(map((x,y)->x>y, first(globalRange).I, last(baRange).I)) ||
            any(map((x,y)->x<y, last(globalRange).I, first(baRange).I))
        # @warn("out of volume range, keep it as zeros")
        return
    end
    chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer =
            adjust_volume_boundary(ba, chunkGlobalRange, globalRange,
                                   rangeInChunk, rangeInBuffer)
    # finalize to avoid memory leak, see
    # https://discourse.julialang.org/t/understanding-distributed-memory-garbage-collection/8726/2
    #finalize(jobs)
    chunkSize = (last(chunkGlobalRange) - first(chunkGlobalRange) + one(CartesianIndex{N})).I
    #println("processing block in global range: $(cartesian_range2string(globalRange))")
    try
        data = ba.kvStore[ cartesian_range2string(chunkGlobalRange) ]
        chk = Codings.decode(data, get_encoding(ba))
        chk = reshape(reinterpret(T, chk), chunkSize)
        @inbounds sharedBuffer[globalRange] = chk[rangeInChunk]
    catch err
        if isa(err, KeyError) && ba.fillMissing
            println("no such key in file system: $(err), will fill this block as zeros")
            return
        else
            println("catch an error while getindex: $err with type of $(typeof(err))")
            rethrow()
        end
    end
    nothing
end

function getindex_sharedarray(ba::BigArray, idxes::Union{UnitRange, Int}...)
    t1 = time()
    sz = map(length, idxes)
    # it seems that the default value is automatically set to zero
    sharedBuffer = OffsetArray(SharedArray{T}(sz), idxes...)

    channelSize = cld( nworkers(), 2 )
    jobs    = RemoteChannel(()->Channel{Tuple}( channelSize ));

    baIter = ChunkIterator(idxes, ba.chunkSize; offset=ba.offset)

    @sync begin
        @async begin
            for iter in baIter
                put!(jobs, iter)
            end
            close(jobs)
        end
        # control the number of concurrent requests here
        for iter in baIter
            @async remote_do(getindex_sharedarray_worker!, WORKER_POOL, ba, jobs, sharedBuffer)
        end
    end
    ret = OffsetArray(sdata(sharedBuffer |> parent), axes(sharedBuffer))
    elapsed = time() - t1 # seconds
    println("cutout speed: $(sizeof(ret)/1024/1024/elapsed) MB/s")
    ret
end


text:
using Documenter, Koma

makedocs(
    modules = [Koma],
    sitename = "Koma.jl: General MRI simulation framework",
    authors = "Carlos Castillo Passi",
    pages = [
        "Home" => "index.md",
        "API" => "API.md"
    ]
)

deploydocs(
    repo = "github.com/cncastillo/Koma.jl.git",
    # versions = ["stable" => "v^", "v#", "dev" => "dev"],
)

#julia --color=yes make.jl


text:
function sum(x, y)
  x + y
end

println(sum(10, 30))

# like js
sum2(x, y) = x + y
println(sum2(10, 20))

# Can use unicode chars
足す(x, y) = x + y
println(足す(1, 3))

function p(v)
    println(v)
    # return nothing
    return nothing
end
p("aaa")

function optionalArgs(n::Int64 = 1)::Int64
    n * 10
end
println(optionalArgs(10))

function varargs(n::Int64, other...) # tuple
    map(x -> x * n, other)
end
println(varargs(3, 1, 3, 2)) # -> (3, 9, 6)
println(varargs(10)) # -> ()

function tup(a::Int64, b::Int64)
    a + b, a * b # return tuple
end
(add, mult) = tup(3, 5)
println(add)
println(mult)

# anonymous function
sum3 = (a::Int64, b::Int64) -> a + b
println(sum3(1, 3))

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
"""
Tableau for the 2-stage stochastic LobattoIIIA-IIIB-IIID method
  Tableau for the 2-stage stochastic LobattoIIIA-IIIB-IIID method
  (based on the deterministic LobattoIIIA-IIIB-IIID due to L. Jay)
  It satisfies the conditions for convergence of order 1.0 for one Wiener process,
  but it doesn't satisfy the conditions for Lagrange-d'Alembert integrators
"""
function TableauStochasticLobattoIIIABD2()

    TableauSISPRK(:StochasticLobattoIIIABD2, CoefficientsLobattoIIIA(2), CoefficientsLobattoIIIA(2),
                                         CoefficientsLobattoIIIB(2), CoefficientsLobattoIIID(2),
                                         CoefficientsLobattoIIIB(2), CoefficientsLobattoIIID(2))
end

"""
Tableau for the 2-stage modified stochastic LobattoIIIA-IIIB method
  Tableau for the 2-stage modified stochastic LobattoIIIA-IIIB method
  Satisfies the conditions for Lagrange-d'Alembert integrators
  and the conditions for convergence of order 1.0 for one Wiener process
"""
function TableauModifiedStochasticStoermerVerlet(c::Number=0.0)

    @assert c ≥ 0.0
    @assert c ≤ 1.0

    a_drift2 = [[c 0.0]
                [c 0.0]]

    b_drift2 = [c, 1-c]

    c_drift2 = [c, c]

    TableauSISPRK(:StochasticModifiedStormerVerlet,
                  CoefficientsLobattoIIIA(2), CoefficientsLobattoIIIA(2),
                  CoefficientsLobattoIIIB(2), CoefficientsRK(typeof(c), :cTableau, 1, a_drift2, b_drift2, c_drift2),
                  CoefficientsLobattoIIIB(2), CoefficientsLobattoIIIB(2))
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 82abdf71-6f38-4957-8bf7-044d36008bbd
using FFTW, SampledSignals, Plots

# ╔═╡ 02e74af4-b5ac-11eb-3971-e79cdb113be1
md"# Discrete Summation fomulas
Discrete summation fomulas can be used to create waveforms with
harmonics at regular intervals.
Discrete summation formulae can be parameterized and combined to make a variety of spectra but here we will look at basic every harmonic (saw-like) and an odd harmonic (square-like) tones"

# ╔═╡ f1e8126d-02c0-42c7-9b84-f84cfc17b391
begin
    fs = 8000
	freqs = SampleBuf(Array(LinRange(0, 2*pi*220, fs)), fs)
	t = Array(0:1/100:1) * 2 * pi
end

# ╔═╡ 65e90c08-5ca0-4ea6-a9d5-2eb2a5a85c81
# A naive sin sum
function sine_sum(a, t, b, n)
	sum(a^i*sin(t + i*b) for i=0:n)
end

# ╔═╡ 94bf7a3f-752c-4130-b908-c8c14a32e247
# An equivalent closed form calculation
function dsf(a, t, b, n)
	"""
	a: ratio between adjacent harmonics
	t: base frequency
	b: frequency offset between harmonics
	   (harmonics of t + mb for 0 <= m <= n)
	n: number of harmonics
	"""
	out = sin(t) - a*sin(t-b)
	out -= a^(n+1)*(sin(t + (n+1)*b) - a*sin(t + n*b))
	out /= 1 + a*a - 2*a*cos(b)
end

# ╔═╡ ea9f553a-b624-46d0-a3db-df1ec5a1dc7b
begin
	# Demonstration that they are the same, try random values
	a, t2, b, n = .234, .34534, .234, 34
	sine_sum(a, t2, b, n), dsf(a, t2, b, n)
end

# ╔═╡ 09eec3c8-55c6-4b69-8c6a-21f64d261713
function normalize(tone)
	m = minimum(tone)
	M = maximum(tone)
	2*(tone.-m)/(M-m) .-1
end

# ╔═╡ 4e5cc60c-5d67-4cd3-afa0-6ab938e22f65
function plotfreqs(tone, title)
	freqs = fft(tone) |> fftshift
	domain = FFTW.fftfreq(length(tone), fs) |> fftshift
    plot(domain, 10*log10.(abs.(freqs)/fs), title=title,
		xlimit=(0, fs/2), ylimit=(-50, 0))
end

# ╔═╡ ed2e86f6-17cd-4ae7-821a-ab935111aca6
function dsfsaw(t)
	dsf(.8, t, t, 16)/3
end

# ╔═╡ c377d413-6dbc-4175-b418-0f71bf2b29aa
dsfsaw_tone = normalize(dsfsaw.(freqs))

# ╔═╡ 989d247f-84d8-46b9-ba18-abaaec50a94d
plotfreqs(dsfsaw_tone, "DSF saw Spectrum")

# ╔═╡ e6930c81-3007-4e35-a695-a8b13f32b0e9
md"DSF can be used to create very bright sounds like a saw or square without aliasing. It also can mimic the sharp cutoff above a limit that a filter would produce. However the amplitude of the harmonics is geometric $(a^n)$ rather than harmonic $(1/n)$ so an exact bandlimited saw or square can't be replicated"

# ╔═╡ Cell order:
# ╟─02e74af4-b5ac-11eb-3971-e79cdb113be1
# ╠═82abdf71-6f38-4957-8bf7-044d36008bbd
# ╠═f1e8126d-02c0-42c7-9b84-f84cfc17b391
# ╠═65e90c08-5ca0-4ea6-a9d5-2eb2a5a85c81
# ╠═94bf7a3f-752c-4130-b908-c8c14a32e247
# ╠═ea9f553a-b624-46d0-a3db-df1ec5a1dc7b
# ╟─09eec3c8-55c6-4b69-8c6a-21f64d261713
# ╟─4e5cc60c-5d67-4cd3-afa0-6ab938e22f65
# ╠═ed2e86f6-17cd-4ae7-821a-ab935111aca6
# ╠═c377d413-6dbc-4175-b418-0f71bf2b29aa
# ╠═989d247f-84d8-46b9-ba18-abaaec50a94d
# ╟─e6930c81-3007-4e35-a695-a8b13f32b0e9

max_stars_repo_name: jnha/music

max_issues_repo_name: jnha/music

max_forks_repo_name: jnha/music


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApimachineryPkgApisMetaV1Initializers <: SwaggerModel
    pending::Any # spec type: Union{ Nothing, Vector{IoK8sApimachineryPkgApisMetaV1Initializer} } # spec name: pending
    result::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Status } # spec name: result

    function IoK8sApimachineryPkgApisMetaV1Initializers(;pending=nothing, result=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1Initializers, Symbol("pending"), pending)
        setfield!(o, Symbol("pending"), pending)
        validate_property(IoK8sApimachineryPkgApisMetaV1Initializers, Symbol("result"), result)
        setfield!(o, Symbol("result"), result)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1Initializers

const _property_map_IoK8sApimachineryPkgApisMetaV1Initializers = Dict{Symbol,Symbol}(Symbol("pending")=>Symbol("pending"), Symbol("result")=>Symbol("result"))
const _property_types_IoK8sApimachineryPkgApisMetaV1Initializers = Dict{Symbol,String}(Symbol("pending")=>"Vector{IoK8sApimachineryPkgApisMetaV1Initializer}", Symbol("result")=>"IoK8sApimachineryPkgApisMetaV1Status")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1Initializers }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1Initializers))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1Initializers }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1Initializers[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1Initializers }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1Initializers[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1Initializers)
    (getproperty(o, Symbol("pending")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgApisMetaV1Initializers }, name::Symbol, val)
end


text:
abstract type AbstractMillModel end

const MillFunction = Union{Flux.Dense, Flux.Chain, Function}

include("arraymodel.jl")
include("bagmodel.jl")
include("productmodel.jl")
include("missingmodel.jl")
include("lazymodel.jl")

import HierarchicalUtils: encode, stringify

reflectinmodel(x, db=d->Flux.Dense(d, 10), da=d->SegmentedMean(d); b = Dict(), a = Dict()) = _reflectinmodel(x, db, da, b, a, "")[1]

function _reflectinmodel(x::AbstractBagNode, db, da, b, a, s)
    im, d = _reflectinmodel(x.data, db, da, b, a, s * encode(1, 1))
    c = stringify(s)
    agg = c in keys(a) ? a[c](d) : da(d)
    bm, d = _reflectinmodel(BagModel(im, agg)(x), db, da, b, a, s)
    BagModel(im, agg, bm), d
end

function _reflectinmodel(x::AbstractProductNode, db, da, b, a, s)
    n = length(x.data)
    mm = [_reflectinmodel(xx, db, da, b, a, s * encode(i, n)) for (i, xx) in enumerate(x.data)]
    im = tuple([i[1] for i in mm]...)
    tm, d = _reflectinmodel(ProductModel(im)(x), db, da, b, a, s)
    ProductModel(im, tm), d
end

function _reflectinmodel(x::ProductNode{T,C}, db, da, b, a, s) where {T<:NamedTuple, C}
    n = length(x.data)
    ks = keys(x.data)
    ms = (;[k => _reflectinmodel(x.data[k], db, da, b, a, s * encode(i, n))[1] for (i, k) in enumerate(ks)]...)
    tm, d = _reflectinmodel(ProductModel(ms)(x), db, da, b, a, s)
    ProductModel(ms, tm), d
end

function _reflectinmodel(x::ArrayNode, db, da, b, a, s)
    c = stringify(s)
    t = c in keys(b) ? b[c](size(x.data, 1)) : db(size(x.data, 1))
    m = ArrayModel(t)
    m, size(m(x).data, 1)
end


text:
arg_dict = Dict{Any, Any}(
# specify the physics and SBP operator
"physics" => "Euler",  # specify physics to run
"operator_type" => "SBPDiagonalE",  # specify SBP operator
"dimensions" => 2,  # this is a two dimensional mesh

# specify temporal and spatial discretization
"run_type" => 5,  # use steady Newton's method
"jac_type" => 3,  # Petsc linear solver
"jac_method" => 2, # use complex step to compute Jacobian, or explicit
                   # Jacobian computation if possible
"real_time" => false,  # this is an steady problem
"order" => 1,  # p = 1 operator
"use_DG" => true,  # use discontinuous galerkin solver
"volume_integral_type" => 2,  # entropy-stable volume integral
"Volume_flux_name" => "IRFlux",  # flux function for volume integrals
"Flux_name" => "HLLFlux",  # numerical flux function for face integrals,
                           # HLL Riemann solver
# specify the problem itself
"IC_name" => "ICIsentropicVortex",  # initial condtiion
"numBC" => 1,  # number of boundary conditions
"BC1" => [0, 1, 2, 3],  # geometric edges to apply the BC to
"BC1_name" => "isentropicVortexBC",   # name of boundary condition

# specify mesh
"smb_name" => "../../meshes/square_vortex_p4_1_.smb",
"dmg_name" => ".null",

# misc options
"res_abstol" => 1e-10,
"res_reltol" => 1e-10,
"itermax" => 20,
"use_lps" => true,
"do_postproc" => true,  # calculate error at end of run
"exact_soln_func" => "ICIsentropicVortex",  # use this function for the exact soltuion (to calculate error)
)


text:
# Run tests on functions in src/by_reference_doc/markets.jl

@testset verbose = true "GET-request endpoints for markets" begin

    @test ~isempty(Spotify.markets_get()[1])

end


text:
"""
$(TYPEDEF)

Pose3 is currently a Euler angle mechanization of three Euclidean translations and three Circular rotation.

Future:
------
- Work in progress on AMP3D for proper non-Euler angle on-manifold operations.
"""
struct Pose3 <: IncrementalInference.InferenceVariable
  dims::Int
  labels::Vector{String}
  manifolds::Tuple{Symbol,Symbol,Symbol,Symbol,Symbol,Symbol}
  # TODO the AMP upgrade is aimed at resolving 3D to Quat/SE3/SP3 -- current Euler angles will be replaced
  Pose3() = new(6, String[], (:Euclid,:Euclid,:Euclid,:Circular,:Circular,:Circular) )
end


text:
################################################################################

module Strategy

@enum Choice begin
	Military
	Corruption
end

function Base.string(s::Choice)
	s == Military && return "military"
	s == Corruption && return "corruption"
	return nothing
end

function parse(s::AbstractString)
	for strat in instances(Choice)
		s == string(strat) && return strat
	end
	return nothing
end

end

function StrategyChoice(s::AbstractString)
	strat = Strategy.parse(s)
	isnothing(strat) && error("\"$(s)\" is not a valid strategy.")
	return strat
end

################################################################################

module Die

@enum Face begin
	Character
	Army
	Muster
	ArmyMuster
	Event
	# Eye # Not really used
	# WillOfTheWest # Not used
end

struct FaceInfo
	str::String
	c::Char
	article::String
end

const FaceInfos = Dict(
	Character => FaceInfo("Character Die",'C',"a"),
	Army => FaceInfo("Army Die",'A',"an"),
	Muster => FaceInfo("Muster Die",'M',"a"),
	ArmyMuster => FaceInfo("Muster/Army Die",'H',"a"),
	Event => FaceInfo("Event Die",'P',"an"),
	# Eye => FaceInfo("Eye",'E',"an"),
	# WillOfTheWest => FaceInfo("Will of the West",'W',"a"),
	)

article(f::Face) = FaceInfos[f].article

Base.string(f::Face) = FaceInfos[f].str
char(f::Face) = FaceInfos[f].c

function parse(c::Char)
	c = uppercase(c)
	for f in instances(Face)
		c == char(f) && return f
	end
	return nothing
end

function parse(s::AbstractString)
	s = filter(!isspace, collect(s)) # remove any spaces
	fs = parse.(s)
	(isempty(s) || any(isnothing.(fs))) && return nothing
	return fs
end

end

function DieFace(c)
	f = Die.parse(c)
	isnothing(f) && error("\'$(c)\' is not a valid die.")
	return f
end

max_stars_repo_name: mvmorin/queller-bot

max_issues_repo_name: mvmorin/queller-bot

max_forks_repo_name: mvmorin/queller-bot


text:
module SymmetryCriteria

using CrystallographyBase:
    CrystalSystem,
    Cubic,
    Hexagonal,
    Tetragonal,
    Trigonal,
    Orthorhombic,
    Monoclinic,
    Triclinic
using ..LinearElasticity: StiffnessMatrix, ComplianceMatrix

export issystem, whichsystem, isisotropic

function symmetry_criteria(::Cubic, x::Union{StiffnessMatrix,ComplianceMatrix})
    return (
        all(iszero, x[1:3, 4:6]),
        all(iszero, x[4, 5:6]),
        iszero(x[5, 6]),
        all(!iszero, (x[1, 1], x[4, 4], x[1, 3])),
        x[1, 1] == x[2, 2] == x[3, 3],
        x[4, 4] == x[5, 5] == x[6, 6],
        x[1, 2] == x[1, 3] == x[2, 3],
    )
end
function symmetry_criteria(::Hexagonal, c::StiffnessMatrix)
    return (
        all(iszero, c[1:3, 4:6]),
        all(iszero, c[4, 5:6]),
        iszero(c[5, 6]),
        all(!iszero, (c[1, 1], c[3, 3], c[4, 4], c[1, 2], c[1, 3], c[6, 6])),
        c[1, 1] == c[2, 2],
        c[4, 4] == c[5, 5],
        c[1, 3] == c[2, 3],
        2c[6, 6] == c[1, 1] - c[1, 2],
    )
end
function symmetry_criteria(::Hexagonal, s::ComplianceMatrix)
    return (
        all(iszero, s[1:3, 4:6]),
        all(iszero, s[4, 5:6]),
        iszero(s[5, 6]),
        all(!iszero, (s[1, 1], s[3, 3], s[4, 4], s[1, 2], s[1, 3], s[6, 6])),
        s[1, 1] == s[2, 2],
        s[4, 4] == s[5, 5],
        s[1, 3] == s[2, 3],
        s[6, 6] == 2(s[1, 1] - s[1, 2]),
    )
end
function symmetry_criteria(::Tetragonal, x::Union{StiffnessMatrix,ComplianceMatrix})
    return (
        all(iszero, x[1:3, 4:5]),
        all(iszero, x[4, 5:6]),
        all(iszero, (x[3, 6], x[5, 6])),
        all(!iszero, (x[1, 1], x[1, 2], x[1, 3], x[3, 3], x[4, 4], x[6, 6])),
        x[1, 1] == x[2, 2],
        x[4, 4] == x[5, 5],
        x[1, 3] == x[2, 3],
        if all(iszero, (x[2, 6], x[1, 6]))  # Tetragonal (I) class, 4mm, -42m, 422, 4/mmm
            true
        else  # Tetragonal (II) class, 4, -4, 4/m
            x[1, 6] == -x[2, 6]
        end,
    )
end
function symmetry_criteria(::Trigonal, c::StiffnessMatrix)
    return (
        all(iszero, c[1:3, 6]),
        all(iszero, c[3, 4:5]),
        iszero(c[4, 5]),
        all(!iszero, c[1, 1:4]),
        all(!iszero, c[2, 2:3]),
        all(!iszero, (c[3, 3], c[4, 4], c[6, 6])),
        c[1, 1] == c[2, 2],
        c[4, 4] == c[5, 5],
        c[1, 3] == c[2, 3],
        2c[6, 6] == c[1, 1] - c[1, 2],
        c[1, 4] == -c[2, 4] == c[5, 6],
        if iszero(c[1, 5])  # Rhombohedral (I) class, 32, -3m, 3m
            all(iszero, (c[2, 5], c[4, 6]))
        else  # Rhombohedral (II) class, 3, -3
            -c[1, 5] == c[2, 5] == c[4, 6]
        end,
    )
end
function symmetry_criteria(::Trigonal, s::ComplianceMatrix)
    return (
        all(iszero, s[1:3, 6]),
        all(iszero, s[3, 4:5]),
        iszero(s[4, 5]),
        all(!iszero, s[1, 1:4]),
        all(!iszero, s[2, 2:3]),
        all(!iszero, (s[3, 3], s[4, 4], s[6, 6])),
        s[1, 1] == s[2, 2],
        s[4, 4] == s[5, 5],
        s[1, 3] == s[2, 3],
        s[6, 6] == 2(s[1, 1] - s[1, 2]),
        2s[1, 4] == -2s[2, 4] == s[5, 6],
        if iszero(s[1, 5])  # Rhombohedral (I) class, 32, -3m, 3m
            all(iszero, (s[2, 5], s[4, 6]))
        else  # Rhombohedral (II) class, 3, -3
            -2s[1, 5] == 2s[2, 5] == s[4, 6]
        end,
    )
end
function symmetry_criteria(::Orthorhombic, x::Union{StiffnessMatrix,ComplianceMatrix})
    return (
        all(iszero, x[1:3, 4:6]),
        all(iszero, x[4, 5:6]),
        iszero(x[5, 6]),
        all(!iszero, x[1:2, 1:3]),
        all(!iszero, (x[3, 3], x[4, 4], x[5, 5], x[6, 6])),
    )
end
function symmetry_criteria(::Monoclinic, x::Union{StiffnessMatrix,ComplianceMatrix})
    return (
        all(iszero, x[1:3, 4]),
        iszero(x[5, 6]),
        all(!iszero, x[1:2, 1:3]),
        all(!iszero, (x[3, 3], x[4, 4], x[5, 5], x[6, 6])),
        if iszero(x[4, 5])  # Diad // x2, standard orientation
            all(iszero, x[1:3, 6]) && all(!iszero, x[1:3, 5]) && !iszero(x[4, 6])
        else
            all(!iszero, x[1:3, 6]) && all(iszero, x[1:3, 5]) && iszero(x[4, 6])
        end,
    )
end
symmetry_criteria(::Triclinic, x::Union{StiffnessMatrix,ComplianceMatrix}) =
    all(!iszero, x.data.data)

issystem(x::Union{StiffnessMatrix,ComplianceMatrix}, system::CrystalSystem) =
    all(symmetry_criteria(system, x))

function whichsystem(x::Union{StiffnessMatrix,ComplianceMatrix})
    for system in (
        Cubic(),
        Hexagonal(),
        Tetragonal(),
        Trigonal(),
        Orthorhombic(),
        Monoclinic(),
        Triclinic(),
    )
        if issystem(x, system)
            return system
        end
    end
end

function isisotropic(c::StiffnessMatrix)
    return all((
        all(iszero, c[1:3, 4:6]),
        all(iszero, c[4, 5:6]),
        iszero(c[5, 6]),
        all(!iszero, (c[1, 1], c[4, 4], c[1, 3])),
        c[1, 1] == c[2, 2] == c[3, 3],
        2c[4, 4] == 2c[5, 5] == 2c[6, 6] == c[1, 1] - c[1, 2],
        c[1, 2] == c[1, 3] == c[2, 3],
    ))
end
function isisotropic(s::ComplianceMatrix)
    return all((
        all(iszero, s[1:3, 4:6]),
        all(iszero, s[4, 5:6]),
        iszero(s[5, 6]),
        all(!iszero, (s[1, 1], s[4, 4], s[1, 3])),
        s[1, 1] == s[2, 2] == s[3, 3],
        s[4, 4] == s[5, 5] == s[6, 6] == 2(s[1, 1] - s[1, 2]),
        s[1, 2] == s[1, 3] == s[2, 3],
    ))
end

end


text:
@testset "110.balanced-binary-tree.jl" begin
    @test is_balanced_binary_tree(TreeNode{Int}([3, 9, 20, nothing, nothing, 15, 7])) ==
        true
    @test is_balanced_binary_tree(TreeNode{Int}([1, 2, 2, 3, 3, nothing, nothing, 4, 4])) ==
        false
    @test is_balanced_binary_tree(nothing) == true
end


text:
using EnvironmentRestore
using Test

@testset "EnvironmentRestore.jl" begin
    # Write your tests here.
end


text:
########################################################################################
# Post-Processing: Smoothing and uncertainty calibration
########################################################################################
function smooth_all!(integ)

    @unpack x, t, diffusions = integ.sol
    @unpack A, Q, Precond = integ.cache
    # x_pred is just used as a cache here

    for i in length(x)-1:-1:2
        dt = t[i+1] - t[i]
        if iszero(dt)
            copy!(x[i], x[i+1])
            continue
        end

        P = Precond(dt)
        PI = inv(P)

        Qh = apply_diffusion(Q, diffusions[i])

        x[i] = P * x[i]
        smooth!(x[i], P*x[i+1], A, Qh, integ)
        any(isnan.(x[i].μ)) && error("NaN mean after smoothing")
        any(isnan.(x[i].Σ)) && error("NaN cov after smoothing")
        x[i] = PI * x[i]
    end
end

function smooth!(x_curr, x_next, Ah, Qh, integ)
    # x_curr is the state at time t_n (filter estimate) that we want to smooth
    # x_next is the state at time t_{n+1}, already smoothed, which we use for smoothing
    @unpack d, q = integ.cache
    @unpack x_tmp = integ.cache

    # @info "smooth!" x_curr.Σ x_next.Σ Ah Qh PI
    if all((Qh) .< eps(eltype(Qh)))
        @warn "smooth: Qh is really small! The system is basically deterministic, so we just \"predict backwards\"."
        return inv(Ah) * x_next
    end

    # Prediction: t -> t+1
    predict!(x_tmp, x_curr, Ah, Qh)

    # Smoothing
    P_p = x_tmp.Σ
    P_p_inv = inv(P_p)
    G = x_curr.Σ * Ah' * P_p_inv
    x_curr.μ .+= G * (x_next.μ .- x_tmp.μ)

    # Joseph-Form:
    K_tilde = x_curr.Σ * Ah' * P_p_inv
    P_s = (
        X_A_Xt(x_curr.Σ, (I - K_tilde*Ah))
        + X_A_Xt(Qh, K_tilde)
        + X_A_Xt(x_next.Σ, G)
    )
    copy!(x_curr.Σ, P_s)

    assert_nonnegative_diagonal(x_curr.Σ)
    # PDMat(Symmetric(x_curr.Σ))

    return nothing
end


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

using EncodedArrays
using Test

using BitOperations

@testset "varlen_io" begin
    function test_encdec(T_read::Type{<:Integer}, f_read::Function, f_write::Function, x::Integer, encoded::AbstractVector{UInt8})
        buf_out = IOBuffer()
        @inferred f_write(buf_out, x)
        enc_data = take!(buf_out)
        @test enc_data == encoded

        buf_in = IOBuffer(enc_data)
        dec_x = @inferred f_read(buf_in, T_read)
        @test typeof(dec_x) == T_read
        @test dec_x == x
        @test eof(buf_in)
    end

    @testset "read_varlen, write_varlen" begin
        f_read = EncodedArrays.read_varlen
        f_write = EncodedArrays.write_varlen
        test_encdec(UInt64, f_read, f_write, UInt64(0x00), [0x00])
        test_encdec(UInt64, f_read, f_write, UInt64(0x7f), [0x7f])
        test_encdec(UInt64, f_read, f_write, UInt64(0x80), [0x80, 0x01])
        test_encdec(UInt64, f_read, f_write, UInt64(0x4c6b5f94759), [0xd9, 0x8e, 0xe5, 0xaf, 0xeb, 0x98, 0x01])
        @test_throws ErrorException test_encdec(UInt32, f_read, f_write, UInt64(0x4c6b5f94759), [0xd9, 0x8e, 0xe5, 0xaf, 0xeb, 0x98, 0x01])
    end

    @testset "read_autozz_varlen, write_autozz_varlen" begin
        f_read = EncodedArrays.read_autozz_varlen
        f_write = EncodedArrays.write_autozz_varlen
        test_encdec(UInt64, f_read, f_write, UInt64(0x4c6b5f94759), [0xd9, 0x8e, 0xe5, 0xaf, 0xeb, 0x98, 0x01])
        test_encdec(Int64, f_read, f_write, zigzagdec(UInt64(0x4c6b5f94759)), [0xd9, 0x8e, 0xe5, 0xaf, 0xeb, 0x98, 0x01])
    end
end # testset


text:
using ShockTube
using Test

@testset "ShockTube.jl" begin
    st = shockcalc(Species("N2"), Species("Ar"), 1.8)
    @test isapprox(ustrip(pressure(st.driver)), 1.426894182009389e6, rtol=2e-3)
    @test isapprox(ustrip(st.u2), 300.0825207636464, rtol=2e-3)
    @test isapprox(ustrip(density(st.shocked)), 3.3867027987446874, rtol=2e-3)
    @test isapprox(ustrip(soundspeed(st.reflected)), 544.5317089952606, rtol=5e-3)
end


text:
using Test, MPI

using ClimateMachine
using ClimateMachine.MPIStateArrays
using ClimateMachine.VariableTemplates: @vars, varsindex
using StaticArrays

ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()
const mpicomm = MPI.COMM_WORLD

const V = @vars begin
    a::Float32
    b::SVector{3, Float32}
    c::SMatrix{3, 8, Float32}
    d::Float32
end

@testset "MPIStateArray varsindex" begin
    Q = MPIStateArray{Float32, V}(mpicomm, ArrayType, 4, 29, 8)
    @test Q.a === view(MPIStateArrays.realview(Q), :, 1:1, :)
    @test Q.b === view(MPIStateArrays.realview(Q), :, 2:4, :)
    @test Q.c === view(MPIStateArrays.realview(Q), :, 5:28, :)
    @test Q.d === view(MPIStateArrays.realview(Q), :, 29:29, :)

    P = similar(Q)
    @test P.a === view(MPIStateArrays.realview(P), :, 1:1, :)
    @test P.b === view(MPIStateArrays.realview(P), :, 2:4, :)
    @test P.c === view(MPIStateArrays.realview(P), :, 5:28, :)
    @test P.d === view(MPIStateArrays.realview(P), :, 29:29, :)

    A = MPIStateArray{Float32}(mpicomm, ArrayType, 4, 29, 8)
    @test_throws ErrorException A.a
end


text:
# Anagrams are two or more words that composed of the same characters but in a different order

function is_anagram(s1::AbstractString, s2::AbstractString)
    # Disable case sensitivity
    s1 = lowercase(s1)
    s2 = lowercase(s2)

    # Similar and different-length strings are not anagrams
    if length(s1) != length(s2) || s1 == s2
        return false
    end

    # Calculate count of every character in the first string
    chr_count = Dict()
    for c in s1
        chr_count[c] = get(chr_count, c, 0) + 1
    end

    # Reduce the count by every character from the second string
    for c in s2
        t = get(chr_count, c, 0) - 1
        if t < 0
            # Got character that not exist in the first string
            return false
        else
            chr_count[c] = t
        end
    end

    # Check all counts to be zeroes
    return all(i->(i==0), values(chr_count))
end

function detect_anagrams(subject::AbstractString, candidates::AbstractArray)
    result = []
    for candidate = candidates
        if is_anagram(subject, candidate)
            push!(result, candidate)
        end
    end
    result
end


text:
"""
mySieve(y, x, control=nothing; basis="Bernstein", K=3, knots = nothing)

A sieve regression implementation.
"""
struct mySieve <: myEstimator
    β::Array{Float64}  # coefficients
    y::Array{Float64}  # response
    X::Array{Float64}  # matrix of regressors
    basis::String # type of polynominal basis
    K::Int # degree of sieve
	knots # knots for splines

    # Define constructor function
    function mySieve(y, x, control=nothing;
            basis="Bernstein", K=3, knots = nothing)
		# Data parameters
		N = length(y)

        # Define matrix of regressors
        ## Calculate from basis
        X = get_basis(x, basis, K, knots)

        ## Check whether add. variables are included
        if !isnothing(control)
            X = hcat(X, control)
        end

        # Estimate sieve regression
        fit_sieve = myLS(y, X)
        β = fit_sieve.β

        # Organize and return output
        new(β, y, X, basis, K, knots)
    end #MYSIEVE
end #MYSIEVE

# Methods ======================================================================

"""
predict(fit::mySieve)

A method to calculate predictions of a mySieve object.
"""
function predict(fit::mySieve; x=nothing, control=nothing)
    # Return fitted values
    if isnothing(x)
		# Calculate and return in fitted values
        return fit.X * fit.β
    else
		# Build new matrix of regressors
		X = get_basis(x, fit.basis, fit.K, fit.knots)
		## Check whether add. variables are included
        if !isnothing(control)
            X = hcat(X, control)
        end
		# Calculate and return fitted values
        return X * fit.β
    end
end #PREDICT.MYSIEVE

# Helper Functions =============================================================
"""
get_basis(x, basis, K, knots)

Function to construct matrix of regressors from a basis.
"""
function get_basis(x, basis, K, knots)
    # Check input
    if basis == "CSplines" || basis == "LSplines" || basis == "CSplines2"
        # Check whether sufficient knots are provided
		if isnothing(knots)
			#println("Warning: Knots defined as quantiles of x.")
			q_knots = collect(1:(K))./(K+1)
			knots = quantile(x[:], q_knots)
		else
			K = length(knots)
		end
    end
    # Comput basis terms
	if basis == "Bernstein"
        X = [binomial(K,k)*(x.^k).*(1 .-x).^(K-k) for k in 0:K] # omit comb.
		X = reduce(hcat,X)
    elseif basis == "Monomial"
        X = [x.^k for k in 0:K]
		X = reduce(hcat,X)
   elseif basis == "CSplines"
		# Calculate splines
		X = [(x.>knots[k]) for k in 1:K] # not very efficient...
		# Add constant
		X = hcat(ones(length(x)), reduce(hcat,X))
    elseif basis == "CSplines2"
		# Calculate splines
		X = [(knots[k+1].>x.>=knots[k]) for k in 1:(K-1)] # not very efficient...
		# Add outer terms
		X = hcat((x.<knots[1]), reduce(hcat,X), (knots[K].<=x))
	elseif basis == "LSplines"
		# Calculate splines
		X = [(x.>knots[k]).*(x.-knots[k]) for k in 1:K]
		# Add constant and linear term
		X = hcat(ones(length(x)), x, reduce(hcat,X))
    end
	# Return matrix of regressors
	return X
end #GET_BASIS


text:
module CircularList

import Base: insert!, delete!, length, size, eltype, iterate, show

export circularlist, length, size, current, previous, next,
    insert!, delete!, shift!, forward!, backward!, jump!,
    eltype, iterate, show, head, tail

"""
Doubly linked list implementation
"""
mutable struct Node{T}
    data::Union{T, Nothing}
    prev::Union{Node{T}, Nothing}
    next::Union{Node{T}, Nothing}
end

"""
List is used to hold a pre-allocated vector of nodes.
"""
mutable struct List{T}
    nodes::Vector{Node{T}}      # preallocated array of nodes
    current::Node{T}            # current "head" or the circular list
    length::Int                 # number of active elements
    last::Int                   # last index of the nodes array
    capacity::Int               # size of the nodes array
end

"Create a circular list with the specified data element."
function circularlist(data::T; capacity = 100) where T
    nodes = [Node{T}(nothing, nothing, nothing) for _ in 1:capacity]
    n = nodes[1]
    n.data = data
    n.prev = n
    n.next = n
    return List(nodes, n, 1, 1, capacity)
end

"Create a circular list from any vector"
function circularlist(vec::AbstractVector{T}; kw...) where T
    CL = circularlist(vec[1]; kw...)
    for i in 2:length(vec)
        insert!(CL, vec[i])
    end
    forward!(CL)  # move one more step to the original head
    return CL
end

"Returns the length of the circular list"
length(CL::List) = CL.length

"Allocates a new uninitialized node in the circular list"
function allocate!(CL::List, T::DataType)
    if CL.last == CL.capacity   # exceeded capacity...auto resize.
        newcapacity = CL.capacity * 2
        additional  = newcapacity - CL.capacity
        CL.nodes = vcat(CL.nodes, [Node{T}(nothing, nothing, nothing) for _ in 1:additional])
        CL.capacity = newcapacity
    end
    CL.length += 1
    CL.last += 1
    return CL.nodes[CL.last]
end

"Insert a new node after the current node and return the new node."
function insert!(CL::List, data)
    cl = CL.current

    n = allocate!(CL, typeof(data))  # make a new node and arrange prev/next pointers
    n.data = data
    n.prev = cl
    n.next = cl.next

    cl.next = n         # fix prev node's next pointer
    n.next.prev = n     # fix next node's prev pointer

    CL.current = n     # move pointer to newly inserted node
    return CL
end

"""
Delete current node and return the previous node.

_Warning_: removed nodes are not reclaimed from memory for simplicity reasons
"""
function delete!(CL::List)
    length(CL) == 1 && error("cannot remove last item in circular list")
    cl = CL.current
    cl.prev.next = cl.next   # fix prev node's next pointer
    cl.next.prev = cl.prev   # fix next node's prev pointer
    CL.current = cl.prev    # reset List's current pointer to prev
    CL.length -= 1
    return CL
end

"""
Shift the current pointer forward or backward.
"""
function shift!(CL::List, steps::Int, direction::Symbol)
    for i in 1:steps
        if direction == :forward
            CL.current = CL.current.next
        elseif direction == :backward
            CL.current = CL.current.prev
        else
            error("Wrong direction: $direction")
        end
    end
    return CL
end

"Shift the current pointer forward."
forward!(CL::List) = shift!(CL, 1, :forward)

"Shift the current pointer backward."
backward!(CL::List) = shift!(CL, 1, :backward)

"Return the current node."
current(CL::List) = CL.current

"Return the previous node."
previous(CL::List) = CL.current.prev

"Return the ext node."
next(CL::List) = CL.current.next

"Return the head of the list (current node)."
head(CL::List) = CL.current

"Return the tail of the list (last node)"
tail(CL::List) = CL.current.prev

jump!(CL::List, n::Node) = CL.current = n

"Iteration protocol implementation."
function iterate(CL::List, (el, i) = (CL.current, 1))
    i > CL.length && return nothing
    return (el.data, (el.next, i + 1))
end

"Return the element type of the list"
eltype(CL::List) = typeof(CL.current.data)

"Return the size of the list."
size(CL::List) = (CL.length, )

"Show list."
function show(io::IO, CL::List{T}) where T
    print(io, "CircularList.List(")
    i = 1
    for x in CL
        show(io, x)
        i += 1
        i <= length(CL) && print(io, ",")
    end
    print(io, ")")
end

"Show node"
function show(io::IO, node::Node{T}) where T
    print(io, "CircularList.Node(")
    show(io, node.data)
    print(io, ")")
end

end # module


text:
"""
    ObjectiveFunction

Represents an objective function optimized by algorithms.
Subtypes of `ObjectiveFunction` must implement at least
* `compute_value(::ObjectiveFunction, x)` for primal value evaluation
* `compute_gradient(::ObjectiveFunction, x)` for gradient evaluation.
and optionally `compute_value_gradient(::ObjectiveFunction, x)` returning the (primal, gradient) pair.
`compute_gradient` may always use the same storage and return a reference to it.
"""
abstract type ObjectiveFunction end

"""
    compute_value(f::ObjectiveFunction, x; [kwargs...])

Computes the objective `f` at `x`.
"""
function compute_value end

"""
    compute_gradient(f::ObjectiveFunction, x; [kwargs...])

Computes the gradient of `f` at `x`. May return a reference to an internal storage.
"""
function compute_gradient end

"""
    compute_value_gradient(f::ObjectiveFunction, x; [kwargs...])

Computes in one call the pair `(value, gradient)` evaluated at `x`.
By default, calls `compute_value` and `compute_gradient` with keywords `kwargs`
passed down to both.
"""
compute_value_gradient(f::ObjectiveFunction, x; kwargs...) =
    (compute_value(f, x; kwargs...), compute_gradient(f, x; kwargs...))

"""
    SimpleFunctionObjective{F,G,S}

An objective function built from separate primal objective `f(x)` and
in-place gradient function `grad!(storage, x)`.
It keeps an internal storage of type `s` used to evaluate the gradient in-place.
"""
struct SimpleFunctionObjective{F,G,S} <: ObjectiveFunction
    f::F
    grad!::G
    storage::S
end

compute_value(f::SimpleFunctionObjective, x) = f.f(x)
function compute_gradient(f::SimpleFunctionObjective, x)
    f.grad!(f.storage, x)
    return f.storage
end

"""
    StochasticObjective{F, G, XT, S}(f::F, grad!::G, xs::XT, storage::S)

Represents a composite function evaluated with stochastic gradient.
`f(θ, x)` evaluates the loss for a single data point `x` and parameter `θ`.
`grad!(storage, θ, x)` adds to storage the partial gradient with respect to data point `x` at parameter `θ`.
`xs` must be an indexable iterable (`Vector{Vector{Float64}}` for instance).
Functions using a `StochasticObjective` have optional keyword arguments `rng`, `batch_size`
and `full_evaluation` controlling whether the function should be evaluated over all data points.

Note: `grad!` must **not** reset the storage to 0 before adding to it.
"""
struct StochasticObjective{F,G,XT,S} <: ObjectiveFunction
    f::F
    grad!::G
    xs::XT
    storage::S
end

function compute_value(
    f::StochasticObjective,
    θ;
    batch_size::Integer=length(f.xs),
    rng=Random.GLOBAL_RNG,
    full_evaluation=false,
)
    (batch_size, rand_indices) = if full_evaluation
        (length(f.xs), eachindex(f.xs))
    else
        (batch_size, rand(rng, eachindex(f.xs), batch_size))
    end
    return sum(f.f(θ, f.xs[idx]) for idx in rand_indices) / batch_size
end

function compute_gradient(
    f::StochasticObjective,
    θ;
    batch_size::Integer=length(f.xs) ÷ 10 + 1,
    rng=Random.GLOBAL_RNG,
    full_evaluation=false,
)
    (batch_size, rand_indices) = _random_indices(f, batch_size, full_evaluation; rng=rng)

    f.storage .= 0
    for idx in rand_indices
        f.grad!(f.storage, θ, f.xs[idx])
    end
    f.storage ./= batch_size
    return f.storage
end

function compute_value_gradient(
    f::StochasticObjective,
    θ;
    batch_size::Integer=length(f.xs) ÷ 10 + 1,
    rng=Random.GLOBAL_RNG,
    full_evaluation=false,
)
    (batch_size, rand_indices) = _random_indices(f, batch_size, full_evaluation; rng=rng)
    # map operation, for each index, computes value and gradient
    f_val = zero(eltype(θ))
    f.storage .= 0
    for idx in rand_indices
        @inbounds x = f.xs[idx]
        f_val += f.f(θ, x)
        f.grad!(f.storage, θ, x)
    end
    f.storage ./= batch_size
    f_val /= batch_size
    return (f_val, f.storage)
end

function _random_indices(f::StochasticObjective, batch_size::Integer, full_evaluation::Bool; rng=Random.GLOBAL_RNG)
    if full_evaluation
        return (length(f.xs), eachindex(f.xs))
    end
    return (batch_size, rand(rng, eachindex(f.xs), batch_size))
end


text:
module timecorrelations

export correlation, spectrum, correlation2spectrum

using QuantumOpticsBase
using ..timeevolution, ..steadystate

using FFTW

"""
    timecorrelations.correlation([tspan, ]rho0, H, J, op1, op2; <keyword arguments>)

Calculate two time correlation values ``⟨A(t)B(0)⟩``.

The calculation is done by multiplying the initial density operator
with ``B`` performing a time evolution according to a master equation
and then calculating the expectation value ``\\mathrm{Tr} \\{A ρ\\}``

Without the `tspan` argument the points in time are chosen automatically from
the ode solver and the final time is determined by the steady state termination
criterion specified in [`steadystate.master`](@ref).

# Arguments
* `tspan`: Points of time at which the correlation should be calculated.
* `rho0`: Initial density operator.
* `H`: Operator specifying the Hamiltonian.
* `J`: Vector of jump operators.
* `op1`: Operator at time `t`.
* `op2`: Operator at time `t=0`.
* `rates=ones(N)`: Vector or matrix specifying the coefficients (decay rates)
        for the jump operators.
* `Jdagger=dagger.(J)`: Vector containing the hermitian conjugates of the jump
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function correlation(tspan::Vector, rho0::Operator{B,B}, H::AbstractOperator{B,B}, J::Vector,
                     op1::AbstractOperator{B,B}, op2::AbstractOperator{B,B};
                     rates::Union{Vector, Matrix, Nothing}=nothing,
                     Jdagger::Vector=dagger.(J),
                     kwargs...) where B<:Basis
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    function fout(t, rho)
        expect(op1, rho)
    end
    t,u = timeevolution.master(tspan_, op2*rho0, H, J; rates=rates, Jdagger=Jdagger,
                        fout=fout, kwargs...)
    u
end

function correlation(rho0::Operator{B,B}, H::AbstractOperator{B,B}, J::Vector,
                     op1::AbstractOperator{B,B}, op2::AbstractOperator{B,B};
                     tol=1e-4,
                     rates::Union{Vector, Matrix, Nothing}=nothing,
                     Jdagger::Vector=dagger.(J),
                     kwargs...) where B<:Basis
    op2rho0 = op2*rho0
    exp1 = expect(op1, op2rho0)
    function fout(t, rho)
        expect(op1, rho)
    end
    t,u = steadystate.master(H, J; rho0=op2rho0, tol=tol, fout=fout,
                       rates=rates, Jdagger=Jdagger, save_everystep=true,kwargs...)
end

"""
    timecorrelations.spectrum([omega_samplepoints,] H, J, op; <keyword arguments>)

Calculate spectrum as Fourier transform of a correlation function

This is done with the Wiener-Khinchin theorem

```math
S(ω, t) = 2\\Re\\left\\{\\int_0^{∞} dτ e^{-iωτ}⟨A^†(t+τ)A(t)⟩\\right\\}
```

The argument `omega_samplepoints` gives the list of frequencies where ``S(ω)``
is caclulated. A corresponding list of times is calculated internally by means
of a inverse discrete frequency fourier transform. If not given, the
steady-state is computed before calculating the auto-correlation function.

Without the `omega_samplepoints` arguments the frequencies are chosen
automatically.

# Arguments
* `omega_samplepoints`: List of frequency points at which the spectrum
        is calculated.
* `H`: Operator specifying the Hamiltonian.
* `J`: Vector of jump operators.
* `op`: Operator for which the auto-correlation function is calculated.
* `rho0`: Initial density operator.
* `tol=1e-4`: Tracedistance used as termination criterion.
* `rates=ones(N)`: Vector or matrix specifying the coefficients for the
        jump operators.
* `Jdagger=dagger.(J)`: Vector containing the hermitian conjugates of the
        jump operators. If they are not given they are calculated automatically.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function spectrum(omega_samplepoints,
                H::AbstractOperator{B,B}, J::Vector, op::AbstractOperator{B,B};
                rho0::Operator{B,B}=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))),
                tol=1e-4,
                rho_ss::Operator{B,B}=steadystate.master(H, J; tol=tol, rho0=rho0)[end][end],
                kwargs...) where B<:Basis
    domega = minimum(diff(omega_samplepoints))
    dt = 2*pi/abs(omega_samplepoints[end] - omega_samplepoints[1])
    T = 2*pi/domega
    tspan = [0.:dt:T;]
    exp_values = correlation(tspan, rho_ss, H, J, dagger(op), op, kwargs...)
    S = 2dt.*fftshift(real(fft(exp_values)))
    return omega_samplepoints, S
end

function spectrum(H::AbstractOperator{B,B}, J::Vector, op::AbstractOperator{B,B};
                rho0::Operator{B,B}=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))),
                tol=1e-4,
                rho_ss::Operator{B,B}=steadystate.master(H, J; tol=tol)[end][end],
                kwargs...) where B<:Basis
    tspan, exp_values = correlation(rho_ss, H, J, dagger(op), op, tol=tol, kwargs...)
    dtmin = minimum(diff(tspan))
    T = tspan[end] - tspan[1]
    tspan = Float64[0.:dtmin:T;]
    n = length(tspan)
    omega = mod(n, 2) == 0 ? [-n/2:n/2-1;] : [-(n-1)/2:(n-1)/2;]
    omega .*= 2pi/T
    return spectrum(omega, H, J, op; tol=tol, rho_ss=rho_ss, kwargs...)
end

"""
    timecorrelations.correlation2spectrum(tspan, corr; normalize_spec)

Calculate spectrum as Fourier transform of a correlation function with a given correlation function.

# Arguments
* `tspan`: List of time points corresponding to the correlation function.
* `corr`: Correlation function of which the Fourier transform is to be calculated.
* `normalize_spec`: Specify if spectrum should be normalized to its maximum.
"""
function correlation2spectrum(tspan, corr::Vector{T}; normalize_spec::Bool=false) where T <: Number
  n = length(tspan)
  if length(corr) != n
    ArgumentError("tspan and corr must be of same length!")
  end

  dt = tspan[2] - tspan[1]
  for i=2:length(tspan)-1
    if !isapprox(tspan[i+1] - tspan[i], dt)
      throw(ArgumentError("tspan must be equidistant!"))
    end
  end

  tmax = tspan[end] - tspan[1]
  omega = mod(n, 2) == 0 ? [-n/2:n/2-1;] : [-(n-1)/2:(n-1)/2;]
  omega .*= 2pi/tmax
  spec = 2dt.*fftshift(real(fft(corr)))

  omega, normalize_spec ? spec./maximum(spec) : spec
end

end # module


text:
using Plots

using ForwardDiff
eta_d = 1.0
eta_e = 0.5

pot(eta, e) = 1/eta * log(exp(eta*e) + exp(-eta*e))

e(y, y0) = abs(y0 - y)

yd = 0.5
ye = 1.0

potd(y) = pot(eta_d, e(y, yd))
pote(y) = pot(eta_e, e(y, ye))

pot2(y) = (pot(eta_d, e(y, yd)) + pot(eta_e, e(y, ye)))/2

# plot(potd, label="d")
# plot!(pote, label="e")
# plot!(pot2, label="sum")

nabla_pot2(y) = ForwardDiff.derivative(pot2, y)

y = -2.0:0.01:2.0
p2 = pot2.(y)
pd = potd.(y)
pe = pote.(y)

dp2 = [nabla_pot2(_y) for _y in y]

plot(y, p2, label="p2")
plot!(y, dp2,label="dp2")
plot!(y, pd, label="pd")
plot!(y, pe, label="pe")


text:
struct SChar{C} <: AbstractChar
    function SChar{C}() where C
        !(C isa Char) && throw(ArgumentError("SChar only supports static Char storage, got $(typeof(C))"))
        return new{C}()
    end
end
SChar(c::Char) = SChar{c}()
Base.values(::SChar{C}) where C = C::Char

Base.isvalid(::SChar{C}) where C = isvalid(C::Char)
Base.show(io::IO, x::SChar{C}) where C = print(io, C::Char)

struct SSymbol{S}
    function SSymbol{S}() where S
        !(S isa Symbol) && throw(ArgumentError("SSymbol only supports static Symbol storage, got $(typeof(S))"))
        return new{S}()
    end
end

SSymbol(s::Symbol) = SSymbol{s}()
Base.values(::SSymbol{S}) where S = S::Symbol

Base.show(io::IO, x::SSymbol{S}) where S = print(io, S::Symbol)

const AbstractSymbol = Union{SSymbol{<:Any},Symbol}


text:
arg_dict = Dict{String, Any}(
"physics" => "Advection",
"var4" => 1,
"smb_name" => "a string"
)


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# Inner produce of x and y
function inner( x, y )
    s = zero(eltype(x))
    @simd for i=1:length(x)
        @inbounds s += x[i]*y[i]
    end
    s
end

function flog_inner( m, x, y )
    s = zero(eltype(x))
    for j=1:m
        s += inner(x,y)
    end
    s
end

for t in [Float32,Float64]
    n = 1000
    x = rand(t,n)
    y = rand(t,n)
    bits = 8*sizeof(t)
    @timeit(flog_inner(100,x,y), "inner_$bits", "SIMD inner product for type $t", "SIMD")
end


text:
immutable PoissBinAliasSampler <: Sampleable{Univariate,Discrete}
    table::AliasTable
end

function PoissBinAliasSampler(p::AbstractVector)
    n = length(p)
    pv = poissonbinomial_pdf_fft(p)
    alias = Array(Int, n+1)
    StatsBase.make_alias_table!(pv, 1.0, pv, alias)
    BinomialAliasSampler(AliasTable(pv, alias, RandIntSampler(n+1)))
end

function PoissBinAliasSampler(d::PoissonBinomial)
    n = length(d.p)
    alias = Array(Int, n+1)
    pv = Array(Float64, n+1)
    StatsBase.make_alias_table!(d.pmf, 1.0, pv, alias)
    BinomialAliasSampler(AliasTable(pv, alias, RandIntSampler(n+1)))
end

rand(s::PoissBinAliasSampler) = rand(s.table) - 1


text:
export RelaxVankaFaces,RelaxVankaFacesColor,getVankaFacesPreconditioner,cellColor,cellRBColor

# const libdir = "C:/Users/YAEL/Dropbox/Eran/JuliaCode/Active/Multigrid.jl/src/";

# const libdir = "C:/Users/etrieste/Dropbox/Eran/JuliaCode/Active/Multigrid.jl/src/";
# const lib = string(libdir,"Vanka.dll");

function loc2cs3D(loc::Array{Int64,1},n::Array{Int64,1})
@inbounds cs = loc[1] + (loc[2]-1)*n[1] + (loc[3]-1)*n[1]*n[2];
return cs;
end

export loc2cs
function loc2cs(loc::Array{Int64,1},n::Array{Int64,1})
if length(n)==2
	@inbounds cs = loc[1] + (loc[2]-1)*n[1];
else
	@inbounds cs = loc[1] + (loc[2]-1)*n[1] + (loc[3]-1)*n[1]*n[2];
end
return cs;
end

function getVankaVariablesOfCell(i::Array{Int64},n::Array{Int64},nf::Array{Int64},Idxs::Array{Int64},includePressure::Bool)
if includePressure
	if length(i)==2
		@inbounds t1 = i[1] + (i[2]-1)*(n[1]+1);
		@inbounds t2 = nf[1] + i[1] + (i[2]-1)*n[1];
		Idxs[1] = t1;
		Idxs[2] = t1+1;
		Idxs[3] = t2;
		@inbounds Idxs[4] = t2 + n[1];
		@inbounds Idxs[5] = nf[2] + t2;
	else
		@inbounds t1 = loc2cs3D(i,n + [1;0;0]);
		Idxs[1] = t1;
		Idxs[2] = t1+1;
		@inbounds t2 = nf[1] + loc2cs3D(i,n + [0;1;0]);
		Idxs[3] = t2;
		@inbounds Idxs[4] = t2+n[1];
		@inbounds t3 = nf[1] + nf[2] + loc2cs3D(i,n); # it should be for n + [0;0;1] but it is the same for n
		Idxs[5] = t3;
		@inbounds Idxs[6] = t3+n[1]*n[2];
		@inbounds Idxs[7] = nf[3] + t3;
	end
else
	if length(i)==2
		@inbounds t1 = i[1] + (i[2]-1)*(n[1]+1);
		@inbounds t2 = nf[1] + i[1] + (i[2]-1)*n[1];
		Idxs[1] = t1;
		Idxs[2] = t1+1;
		Idxs[3] = t2;
		@inbounds Idxs[4] = t2 + n[1];
	else
		@inbounds t1 = loc2cs3D(i,n + [1;0;0]);
		Idxs[1] = t1;
		Idxs[2] = t1+1;
		@inbounds t2 = nf[1] + loc2cs3D(i,n + [0;1;0]);
		Idxs[3] = t2;
		@inbounds Idxs[4] = t2+n[1];
		@inbounds t3 = nf[1] + nf[2] + loc2cs3D(i,n); # it should be for n + [0;0;1] but it is the same for n
		Idxs[5] = t3;
		@inbounds Idxs[6] = t3 + n[1]*n[2];
	end
end

# getVankaVariablesOfCell(long long *i,long long *n,long long *nf,long long *Idxs,includePressure, long long dim)
# Idxs2 = copy(Idxs);
# ccall((:getVankaVariablesOfCell,lib),Void,(Ptr{Int64},Ptr{Int64},Ptr{Int64},Ptr{Int64},Int64, Int64,),i,n,nf,Idxs2,convert(Int64,includePressure),length(n));
# if sum(abs(Idxs - Idxs2)) != 0
	# error("Fix getVankaVariablesOfCell in C");
# end
return Idxs;
end

export cs2loc
function cs2loc(cs_loc::Int64,n::Array{Int64,1})
	if length(n)==3
		@inbounds loc1 = mod(cs_loc-1,n[1])+1;
		@inbounds loc2 = div(mod(cs_loc-1,n[1]*n[2]),n[1]) + 1;
		@inbounds loc3 = div(cs_loc-1,n[1]*n[2])+1;
		return [loc1;loc2;loc3];
	else
		@inbounds loc1 = mod(cs_loc-1,n[1]) + 1;
		@inbounds loc2 = div(cs_loc-1,n[1])+1;
		return [loc1;loc2];
	end
end

function cellRBColor(i::Array{Int64})
return mod(sum(i),2)+1;
end

function cellColor(i::Array{Int64})
color = 0;
if length(i)==2
	if mod(i[1],2)==1
		color =  mod(i[2],2)==1 ? 1 : 2;
	else
		color =  mod(i[2],2)==1 ? 3 : 4;
	end
else
	if mod(i[1],2)==1
		if mod(i[2],2)==1
			color =  mod(i[3],2)==1 ? 1 : 2;
		else
			color =  mod(i[3],2)==1 ? 3 : 4;
		end
	else
		if mod(i[2],2)==1
			color =  mod(i[3],2)==1 ? 5 : 6;
		else
			color =  mod(i[3],2)==1 ? 7 : 8;
		end
 	end
end

# cellColor(long long *i,long long dim)
# ccolor = ccall((:cellColor,lib),Int16,(Ptr{Int64},Int64,),i,length(i));
# if color != ccolor
	# error("Fix cellColor in C");
# end
return color;
end

# n = 100;
# A = sprandn(n,n,10/n) + 1im*sprandn(n,n,10/n);
# AT = A';

# for k=1:20
	# I = sort(randperm(n)[1:10]);
	# Acc = full(AT[I,I])';
	# AccNew = getDenseBlockFromAT(AT,I);
	# println(vecnorm(Acc-AccNew))
# end

function getDenseBlockFromAT(AT::SparseMatrixCSC,Idxs::Array{Int64},Acc::Array)
	# Acc = zeros(eltype(AT),length(Idxs),length(Idxs));
	Acc[:] = 0.0;
	for t = 1:length(Idxs)
		ii = AT.colptr[Idxs[t]];
		jj = 1;
		while ii < AT.colptr[Idxs[t]+1] && jj <= length(Idxs)
			if AT.rowval[ii] ==  Idxs[jj]
				Acc[t,jj] = conj(AT.nzval[ii]);
				ii+=1;
				jj+=1;
			elseif AT.rowval[ii] >  Idxs[jj]
				jj+=1;
			else
				ii+=1;
			end
		end
	end
	return Acc;
end

const VankaBlockType = Complex64

function getVankaFacesPreconditioner(AT::SparseMatrixCSC,M::RegularMesh,w::Float64,includePressure::Bool)
nf = 0;
n = M.n;
blockSize = 0;
if length(M.n)==2
	nf = [prod(n + [1; 0]),prod(n + [0; 1])];
	blockSize = includePressure ? 5 : 4;
else
	nf = [prod(n + [1; 0; 0]),prod(n + [0; 1; 0]),prod(n + [0; 0; 1])];
	blockSize = includePressure ? 7 : 6;
end

LocalBlocks = zeros(VankaBlockType,blockSize*blockSize,prod(M.n));

Acc = zeros(eltype(AT),blockSize,blockSize);
Idx_i = zeros(Int64,blockSize);
for ii = 1:prod(M.n)
	i = cs2loc(ii,M.n);
	Idx_i = getVankaVariablesOfCell(i,n,nf,Idx_i,includePressure);
	## THESE LINES ARE REALLY SLOW
	# Acc1 = AT[Idx_i,Idx_i];
	# Acc1 = full(Acc1');
	Acc = getDenseBlockFromAT(AT,Idx_i,Acc)
	AccInv = convert(Array{VankaBlockType},(w.*inv(Acc))');
	LocalBlocks[:,ii] = AccInv[:];
end

return LocalBlocks;
end

# function RelaxVankaFacesColorParallel(AT::SparseMatrixCSC,r::ArrayTypes,x::ArrayTypes,b::ArrayTypes,D::Array{VankaBlockType,2},numit::Int64,numCores::Int64,M::RegularMesh,includePressure::Bool,numCores::Int64)
	# # gcc -O3 -fopenmp -shared -fpic -DBUILD_DLL Vanka.c -o Vanka.dll
	# # const spmatveclib  = abspath(joinpath(splitdir(Base.source_path())[1],"..","deps","builds","ParSpMatVec"))

# end

function RelaxVankaFacesColor(AT::SparseMatrixCSC,x::ArrayTypes,b::ArrayTypes,y::ArrayTypes,D::Array{VankaBlockType,2},numit::Int64,numCores::Int64,M::RegularMesh,includePressure::Bool)
# function RelaxVankaFaces(AT::SparseMatrixCSC,M::RegularMesh)
	n = M.n;
	dim = M.dim;
	blockSize = 0;
	nf = 0;
	if dim==2
		nf = [prod(n + [1; 0]),prod(n + [0; 1])];
		blockSize = includePressure ? 5 : 4;
	else
		# Face sizes
		nf = [prod(n + [1; 0; 0]),prod(n + [0; 1; 0]),prod(n + [0; 0; 1])];
		blockSize = includePressure ? 7 : 6;
	end
	Idxs = zeros(Int64,blockSize);
	parallel = false;
	if parallel==false
		# y_t = copy(x);
		for k=1:numit
			for color = 1:(2^dim)
				y[:] = x;
				for i = 1:prod(n)
					i_vec = cs2loc(i,n);
					if cellColor(i_vec)==color
						Idxs = getVankaVariablesOfCell(i_vec,n,nf,Idxs,includePressure);
						r = computeResidualAtIdx(AT,b,x,Idxs);
						x[Idxs] = x[Idxs] + (reshape(D[:,i],blockSize,blockSize)'*r);
						# void updateSolution(float complex *mat, double complex *x, double complex *r, int n,long long* Idxs){
						# ccall((:updateSolution,lib),Void,(Ptr{VankaBlockType},Ptr{Complex128},Ptr{Complex128},Int16,Ptr{Int64},)
									# ,D[:,i],y ,r,convert(Int16,blockSize),Idxs);
						# if norm(x[Idxs] - y[Idxs]) > 1e-14
						# error("Fix updateSolution in C: ",norm(x[Idxs] - y[Idxs]));
						# end

					end
				end
			end
		end
	else
		# y[:] = 0.0;
		# void RelaxVankaFacesColor(long long *rowptr , double complex *valA ,long long *colA,
							  # long long *n,long long *nf,long long dim,double complex *x,double complex *b,
							  # float complex *D,long long numit,long long includePressure,long long numCores){
		ccall((:RelaxVankaFacesColor,lib),Void,(Ptr{Int64},Ptr{Complex128},Ptr{Int64},Ptr{Int64}, Ptr{Int64}, Int64,Ptr{Complex128},Ptr{Complex128},Ptr{Complex128},Ptr{VankaBlockType},Int64, Int64, Int64,Int64,),
			AT.colptr,AT.nzval,AT.rowval,n, nf, dim,x,b,y,D,numit,convert(Int64,includePressure),length(x),numCores);
	end
	# if norm(x - y_t) > 1e-14
		# error("Fix RelaxVankaFacesColor in C: ",norm(x - y_t));
	# end
	# println("Diff: ",norm(x-y_t));
	return x;
end


text:
using UUIDs
const BM3DDenoise = Base.PkgId(UUID("95fb3b36-088a-43fb-bb1b-b1f34fadbd7d"), "BM3DDenoise")

# Rewrite from ImageIO.jl
function checked_import(pkgid)
    Base.root_module_exists(pkgid) && return Base.root_module(pkgid)
    # If not available, load the library or throw an error.
    Base.require(pkgid)
end

@doc raw"""
    BM3D(σ [, config=bm3d_config()])

The BM3D(sparse 3D transform-domain collaborative filtering) denoising algorithm.

# Arguments

* `σ::Float64` is the variance of the noise.

* `config::bm3d_config` is `BM3DDenoise.bm3d_config`.

# Examples

```julia
img = testimage("lena_color_256")

n = AdditiveWhiteGaussianNoise(0.1)
noisy_img = apply_noise(img, n)

# use default arguments
f_denoise = BM3D(0.1)
denoised_img = reduce_noise(noisy_img, f_denoise)
```

See also: [`reduce_noise`](@ref), [`reduce_noise!`](@ref)
"""
struct BM3D <: AbstractImageDenoiseAlgorithm
    """degree of filtering"""
    σ::Float64
    """bm3d_config"""
    config
    function BM3D(σ, config)
        σ > 0 || @warn "σ is supposed to be positive"
        new(σ, config)
    end
end
BM3D(σ) = BM3D(σ, Base.invokelatest(checked_import(BM3DDenoise).bm3d_config))

function (f::BM3D)(out::AbstractArray{T},
                           img::AbstractArray) where T
    axes(out) == axes(img) || ArgumentError("Images should have the same axes.")
    out .= T.(Base.invokelatest(checked_import(BM3DDenoise).bm3d, img, f.σ, f.config))
    return out
end


text:
## ClusterSpace methods
for f in [:levels,:ground_state]
    @eval $(f)(c::ClusterSpace{<:NLevelSpace}, args...) = $(f)(c.original_space, args...)
end

has_cluster(::HilbertSpace,args...) = false
has_cluster(::ClusterSpace,args...) = true
function has_cluster(h::ProductSpace)
    for space in h.spaces
        has_cluster(space) && return true
    end
    return false
end
has_cluster(h::ProductSpace,aon) = has_cluster(h.spaces[get_i(aon)])
has_cluster(op::QNumber,args...) = has_cluster(hilbert(op),args...)
has_cluster(avg::Average,args...) = has_cluster(undo_average(avg),args...)

# ClusterAon methods
Base.hash(c::T, h::UInt) where T<:ClusterAon = hash(T, hash(c.i, hash(c.j, h)))
Base.getindex(v::Vector{<:HilbertSpace}, c::ClusterAon) = v[c.i]
get_i(x::Integer) = x
get_j(x::Integer) = x
get_i(x::ClusterAon) = x.i
get_j(x::ClusterAon) = x.j
Base.length(::ClusterAon) = 1

extract_names(names::Vector, i::Int) = names[i]
extract_names(names::Vector, c::ClusterAon) = names[c.i][c.j]
function extract_names(names::Vector, v::Vector)
    [extract_names(names, v_) for v_ in v]
end

Base.isequal(c1::T,c2::T) where T<:ClusterAon = (c1.i==c2.i && c1.j==c2.j)
Base.isless(i::Int,c::ClusterAon) = isless(i,c.i)
Base.isless(c::ClusterAon,i::Int) = isless(c.i,i)
function Base.isless(c1::ClusterAon,c2::ClusterAon)
    if isless(c1.i, c2.i)
        return true
    elseif isequal(c1.i, c2.i)
        return isless(c1.j, c2.j)
    else
        return false
    end
end
Base.iterate(c::ClusterAon, state=1) = isone(state) ? (c,state+1) : nothing

function _cluster(h::ProductSpace, op::QSym, aon::Int)
    order = h.spaces[aon].order
    return _cluster(h, op, aon, order)
end
function _cluster(h::ClusterSpace, op::QSym, aon::Int)
    return _cluster(h, op, aon, h.order)
end
function _cluster(h, op, aon, order)
    ops = QSym[]
    for i=1:order
        name = Symbol(op.name, :_, i)
        aon_i = ClusterAon(aon[1],i)
        op_ = _remake_op(op, h, name, aon_i)
        push!(ops, op_)
    end
    return ops
end

_remake_op(op::Transition, h, name, aon) = Transition(h, name, op.i, op.j, aon)
_remake_op(op::Destroy, h, name, aon) = Destroy(h, name, aon)
_remake_op(op::Create, h, name, aon) = Create(h, name, aon)