text:
###########∇VolumeParts
#this struct/function computes the Jacobian
#vec∇(Θ::VolumeParts, ∇Θ) = reduce(vcat, (p->vec(∇Θ[getproperty(Θ,p)])).(Θ.Πelements))

function ∇VolumeParts!(loss, X::NamedTuple, y, Θ::AbstractVolumeParts)
  function ∇Θ!(out::AbstractVector, Θ)
    ∇Θ = gradient(()->loss(X,y), Flux.params(Θ))
    out .= vec∇(Θ, ∇Θ)
    return Θ
  end

  return ∇Θ!
end

Flux.@functor VolumePartsLSq
#Flux.trainable(a::VolumePartsLSq) = (a.A, a.V₀,)
Flux.params(Θ::VolumePartsLSq) = Flux.params(Θ.A, Θ.V₀,)

Flux.@functor VolumePartsIter
Flux.params(Θ::VolumePartsIter) = Flux.params(Θ.G,)

#WARNING- OPT_INDEX was originally in Capacity.jl- best to leave it here if possible, but it may
#need to be defined earlier
OPT_INDEX = Dict{Symbol,Any}(
  :adam=>Flux.ADAM, #47s (9.7->3.7)
  :radam=>Flux.RADAM, #fail
  :rmsprop=>Flux.RMSProp, #44s (8.4->2.7)
  :nadam=>Flux.NADAM, #70s (9.7->3.6)
  :adamw=>Flux.ADAMW,# 43s, (9.1->3.5)
  :momentum=>Flux.Momentum, #fail
  :nesterov=>Flux.Nesterov,#fail
  :descent=>Flux.Descent, #fail
  #:adadelta=>Flux.ADADelta, #fail #WARNING: uses non-standard learning parameterization
  :adamax=>Flux.AdaMax, #51s (8.8->3.1)
  :adagrad=>Flux.ADAGrad, #44s (8.1->-1.9(!))
  :amsgrad=>Flux.AMSGrad) #42s->(9.6->-0.3)

function test∇VolumeParts!(loss, X::NamedTuple, y, Θ::VolumePartsLSq)
  @info "test ∇VolumeParts!"
  Π = Flux.params(Θ) #formatted parameters
  @info "length: $(length(Π))"
  gs = gradient(()->loss(X,y),Π)

  rms(v) = sqrt(mean(v.^2))
  vrms(m) = rms(vec(m))

  @info gs[Θ.A]

  @info "test gradrms: Θ.A-$(vrms(gs[Θ.A])); Θ.V₀-$(vrms(gs[Θ.V₀]))"

  ∇Θ! = ∇VolumeParts!(loss, X,y, Θ)

  #out .= reduce(vcat, [∇Θ[Θ.A], vec(∇Θ[Θ.V₀])])#vec∇(Θ, ∇Θ)
  out = similar(Θ.Π)
  ∇Θ!(out, Θ)
  NvecA::Int = length(vec(Θ.A))
  NvecV₀::Int = length(vec(Θ.V₀))
  @info "test ∇Θ!: Θ.A-$(rms(out[1:NvecA]));" *
    " Θ.V₀-$(vrms(out[(NvecA+1):(NvecA + NvecV₀)]))"
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
abstract type Perturbation end

"""
Use to model the trip of an AC Branch in in the system. Accepts lines or transformer lines
"""
mutable struct BranchTrip <: Perturbation
    time::Float64
    branch_name::String
end

function ybus_update!(
    ybus::SparseMatrixCSC{ComplexF64, Int},
    b::PSY.ACBranch,
    num_bus::Dict{Int, Int},
    mult::Float64 = 1.0,
)
    arc = PSY.get_arc(b)
    bus_from_no = num_bus[PSY.get_number(arc.from)]
    bus_to_no = num_bus[arc.to.number]

    Y_l = (1 / (PSY.get_r(b) + PSY.get_x(b) * 1im))
    Y11 = Y_l + (1im * PSY.get_b(b).from)
    Y22 = Y_l + (1im * PSY.get_b(b).to)
    ybus[bus_from_no, bus_from_no] += mult * Y11
    Y12 = -Y_l
    ybus[bus_from_no, bus_to_no] += mult * Y12
    #Y21 = Y12
    ybus[bus_to_no, bus_from_no] += mult * Y12
    Y22 = Y_l + (1im * PSY.get_b(b).to)
    ybus[bus_to_no, bus_to_no] += mult * Y22
    return
end

function ybus_update!(
    ybus::SparseMatrixCSC{ComplexF64, Int},
    b::PSY.TapTransformer,
    num_bus::Dict{Int, Int},
    mult::Float64 = 1.0,
)
    arc = PSY.get_arc(b)
    bus_from_no = num_bus[PSY.get_number(arc.from)]
    bus_to_no = num_bus[arc.to.number]

    Y_t = 1 / (PSY.get_r(b) + PSY.get_x(b) * 1im)
    c = 1 / PSY.get_tap(b)
    b = PSY.get_primary_shunt(b)

    Y11 = (Y_t * c^2)
    ybus[bus_from_no, bus_from_no] += mult * Y11
    Y12 = (-Y_t * c)
    ybus[bus_from_no, bus_to_no] += mult * Y12
    #Y21 = Y12
    ybus[bus_to_no, bus_from_no] += mult * Y12
    Y22 = Y_t
    ybus[bus_to_no, bus_to_no] += mult * (Y22 + (1im * b))

    return
end

function ybus_update!(
    ybus::SparseMatrixCSC{ComplexF64, Int},
    b::PSY.PhaseShiftingTransformer,
    num_bus::Dict{Int, Int},
    mult::Float64 = 1.0,
)
    arc = PSY.get_arc(b)
    bus_from_no = num_bus[PSY.get_number(arc.from)]
    bus_to_no = num_bus[arc.to.number]

    Y_t = 1 / (PSY.get_r(b) + PSY.get_x(b) * 1im)
    tap = (PSY.get_tap(b) * exp(PSY.get_α(b) * 1im))
    c_tap = (PSY.get_tap(b) * exp(-1 * PSY.get_α(b) * 1im))
    b = PSY.get_primary_shunt(b)

    Y11 = (Y_t / abs(tap)^2)
    ybus[bus_from_no, bus_from_no] += mult * Y11
    Y12 = (-Y_t / c_tap)
    ybus[bus_from_no, bus_to_no] += mult * Y12
    Y21 = (-Y_t / tap)
    ybus[bus_to_no, bus_from_no] += mult * Y21
    Y22 = Y_t
    ybus[bus_to_no, bus_to_no] += mult * (Y22 + (1im * b))
    return
end

function ybus_update!(
    ybus::SparseMatrixCSC{ComplexF64, Int},
    fa::PSY.FixedAdmittance,
    num_bus::Dict{Int, Int},
    mult::Float64 = 1.0,
)
    bus = PSY.get_bus(fa)
    bus_ix = num_bus[PSY.get_number(bus)]
    ybus[bus_ix, bus_ix] += mult * fa.Y
    return
end

function ybus_update!(integrator_params, branch::PSY.ACBranch, mult::Float64)
    ybus_update!(integrator_params.Ybus, branch, integrator_params.lookup, mult)
    return
end

function get_affect(sys::PSY.System, pert::BranchTrip)
    branch = first(
        PSY.get_components(PSY.ACBranch, sys, x -> PSY.get_name(x) == pert.branch_name),
    )
    return (integrator) -> begin
        ybus_update!(integrator.p, branch, -1.0)
    end
    return
end

# Pending implementation
#=
mutable struct BusTrip <: Perturbation
    time::Float64
    bus::PSY.Bus
end

function get_affect(::PSY.System, pert::BusTrip)
    return (integrator) -> begin
        0.0
    end
end

mutable struct ThreePhaseFault <: Perturbation
    time::Float64
    bus::PSY.Bus
    impedance::Complex
end

function get_affect(::PSY.System, pert::ThreePhaseFault)
    return (integrator) -> begin
        0.0
    end
end
=#

"""
Use to model a change in the network by switching the underlying Ybus in the simulation
"""
mutable struct NetworkSwitch <: Perturbation
    time::Float64
    Ybus::SparseMatrixCSC{Complex{Float64}, Int}
end

function get_affect(::PSY.System, pert::NetworkSwitch)
    return (integrator) -> begin
        # TODO: This code can be more performant using SparseMatrix methods
        for (i, v) in enumerate(pert.Ybus)
            integrator.p.Ybus[i] = v
        end
    end
end

"""
Use to model control reference changes in devices of the model
"""
mutable struct ControlReferenceChange <: Perturbation
    time::Float64
    device::PSY.DynamicInjection
    signal_index::Int
    ref_value::Float64
end

function get_affect(system::PSY.System, pert::ControlReferenceChange)
    device = PSY.get_component(typeof(pert.device), system, PSY.get_name(pert.device))
    pert.device = device
    return (integrator) -> begin
        control_ref = PSY.get_ext(device)[CONTROL_REFS]
        return control_ref[pert.signal_index] = pert.ref_value
    end
end

ext: jl


text:
using ExprRules
using AbstractTrees

using BenchmarkTools
using Random

include("cas_theory.jl")
include("cas_simplify.jl")

grammar = @grammar begin
    Real = x | y | z | a | b | c     # symbol
    Real = -Real
    Real = Real * Real | Real + Real | Real - Real | Real / Real           # julia expression
    Real = _(Base.rand(1.0:100.0))  # special syntax, eval argument of _() at derivation time
    Real = sin(Real) | cos(Real)        # multiple rules on a single line
    Real = |(1:100)                 # same as Real = 4 | 5 | 6
end

Random.seed!(rand(UInt))

for i ∈ 1:100
    rulenode = rand(RuleNode, grammar, :Real, 10)
    println(get_executable(rulenode, grammar))
end


text:
import ImageView
using Colors, Gtk.ShortNames
using TestImages

grid, frames, c = ImageView.canvasgrid((2,2))
showall(Window("canvasgrid", 800, 600) |> grid)
ImageView.imshow(c[1,1], testimage("lighthouse"))
ImageView.imshow(c[1,2], testimage("mountainstream"))
ImageView.imshow(c[2,1], testimage("moonsurface"))
ImageView.imshow(c[2,2], testimage("mandrill"))


text:
"""
Gray Code Generator
===================

## References

- [knuth]: The Art of Computer Programming, Volume 4 Fascicle 2A, 7.2.1.1 Generating all n-tuples.
- [manber]: Introduction to Algorithms: A Creative Approach, 2.9 Gray Codes.
- [wiki]: https://en.wikipedia.org/wiki/Gray_code#Constructing_an_n-bit_Gray_code
"""
module GrayCode

# From [wiki]:
#
# "Once a binary number appears in $G_n$ it appears in the same position in all longer lists; so it
# makes sense to talk about the reflective Gray code value of a number."
#
# "These characteristics suggest a simple and fast method of translating a binary value into the
# corresponding Gray code. Each bit is inverted if the next higher bit of the input value is set to
# one."
#
# From [knuth]:
#
# "[R]eversing the order of Gray binary code is equivalent to complementing the fist bit."
#
# - - - -
#
# We have
#
#     g(i) = b(i) ⊕ b(i + 1)
#
# So g(i) of the MSbit is b(i) as b(i + 1) is 0. And for converting from Gray code to binary
# (c.f. https://math.stackexchange.com/a/1829662):
#
#     b(i) = g(i) ⊕ b(i + 1)
#
# or,
#
#     b(i) = g(i) ⊕ g(i + 1) ⊕ g(i + 2) ⊕ ...

#===== Gray code <-> Binary =====#

# b -> g  (n-th is 1-based)
"Computes the `n`-th reflective gray code."
function gray_nth(n::Integer)::Vector{Int8}
    @assert n > 0
    n -= 1
    g = xor(n, n >>> 1)
    digits(Int8, g, base=2)
end

# g -> b
function gray_to_binary(g::Vector{Int8})::Vector{Int8}
    b = copy(g)
    for i in (length(g) - 1):-1:1
        b[i] = xor(g[i], b[i + 1])
    end
    b
end

#===== Constructs Gray code by recursion =====#

# A slow implementation, just for proof-of-concept. See [manber].
function gray_codes_of_len(len::Integer)::Vector{Vector{Int8}}
    @assert iseven(len) "Gray code of odd length does not exist."
    @assert len > 0 "Length should be greater than 0."
    gray_codes_of_len_r(len)
end

function gray_codes_of_len_r(len::Integer)::Vector{Vector{Int8}}
    if len == 0
        return [[]]
    elseif len == 2
        g =
        [
            Int8[0],
            Int8[1],
        ]
        return g
    end

    if isodd(len)
        _2k = len - 1
        g_2k = gray_codes_of_len_r(_2k)
        if ispow2(_2k)
            g =
            [
                [[0; g] for g in g_2k]...,
                [1; g_2k[end]],
            ]
            return g
        else
            # There are some unused strings of length `ceil(log2(k * 2))`, one of which is connected
            # to one of the used strings.

            # Oops, this construction is useful for the proof but cannot be elegantly expressed as
            # code.
            len_pow2 = nextpow(2, _2k)
            g_pow2 = gray_codes_of_len_r(len_pow2)  # slow, fvck
            used = Set(g_2k)
            for (s, t) in [(i, mod1(i + 1, len_pow2)) for i in 1:len_pow2]
                if g_pow2[s] in used && !(g_pow2[t] in used)
                    index_2k = findfirst(x -> x == g_pow2[s], g_2k)
                    g =
                    [
                        g_2k[(index_2k + 1):end]
                        g_2k[1:index_2k]
                        [g_pow2[t]]
                    ]
                    return g
                end
            end
            error("Should not reach here.")
        end
    end

    half_len = div(len, 2)
    g_half = gray_codes_of_len_r(half_len)
    g_half_rev = reverse(g_half)

    [
        [[0; g] for g in g_half];
        [[1; g] for g in g_half_rev];
    ]
end

#===== Gray code iterating =====#

# c.f. [knuth] Algorithm G.
mutable struct GrayCodeIter
    n::Int
    code::Union{UInt128, Nothing}
    parity::Int8

    GrayCodeIter(n::Int) = new(n, 0, 0)
end

function Base.iterate(iter::GrayCodeIter, state=iter)
    if iter.code == nothing
        nothing
    else
        curr_code = iter.code
        iter.parity ⊻= 1
        if iter.parity == 1
            iter.code ⊻= 1
        elseif iter.code == 1 << (iter.n - 1)
            iter.code = nothing
        else
            iter.code ⊻= curr_code ⊻ (curr_code - 1) + 1
        end
        (curr_code, iter)
    end
end

Base.length(iter::GrayCodeIter) = UInt128(2) ^ iter.n

# TODO: [knuth] Fig. 14

#===== Unit Tests =====#
using Test
using Random: seed!

# Helpers for testing
function diff_one(a::Vector{T}, b::Vector{T})::Bool where T<:Integer
    len_diff = length(a) - length(b)
    if len_diff == 0
        return sum(xor.(a, b)) == 1
    elseif len_diff == 1
        # does nothing
    elseif len_diff == -1
        a, b = b, a
    else
        return false
    end

    a[end] == 1 && @view(a[1:(end - 1)]) == b
end

function binary_array_to_value(ba::Vector{Int8})::UInt64
    @assert length(ba) <= 64
    reduce((acc, x) -> (acc << 1) + x, @view ba[end:-1:1]; init=UInt64(0))
end

function check_grayness(gray_codes::Vector{Vector{Int8}})
    len = length(gray_codes)
    for i in 1:len
        @test diff_one(gray_codes[i], gray_codes[mod1(i + 1, len)])
    end
    @test length(unique(map(binary_array_to_value, gray_codes))) == len
end

@testset "N-th Gray code" begin
    seed!()
    N = 2^16
    for _ in 1:100
        n = rand(1:N)
        n_minus1 = mod1(n - 1, N)
        n_plus1 = mod1(n + 1, N)
        g = gray_nth(n)
        g_minus1 = gray_nth(n_minus1)
        g_plus1 = gray_nth(n_plus1)
        @test diff_one(g_minus1, g)
        @test diff_one(g, g_plus1)

        b = gray_to_binary(g)
        nb = binary_array_to_value(b) + 1
        @test n == nb
    end
end

@testset "Gray code of given length" begin
    for len in 2:2:64
        gray_codes = gray_codes_of_len(len)
        check_grayness(gray_codes)
    end
end

@testset "Gray code - Chinese ring puzzle" begin
    for n in 2:8
        gc_gen = GrayCodeIter(n)
        gray_codes = [digits(Int8, g, base=2, pad=n) for g in gc_gen]
        check_grayness(gray_codes)
    end
end
end  # module

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
function distance(x,y,nx,ny,x0,x1,y0,y1,vx,vy)

## set parameter
paramsx0 = x0;
paramsy0 = y0;

paramsx1 = x1;
paramsy1 = y1;

paramsnx = nx;
paramsny = ny;

paramsdx = (paramsx1-paramsx0)/paramsnx;
paramsdy = (paramsy1-paramsy0)/paramsny;

if y >= y1
    y = y1-1e-8;
end
if x >= x1
    x = x1-1e-8;
end

savex = x;
savey = y;

## traverse

i,j,valid = getij(x, y, paramsx0,paramsy0,paramsx1,paramsy1,paramsdx,paramsdy,paramsnx,paramsny);
I =[]; J = []; L = [];
while (valid)
    I = [I,i] ;
    J = [J,j];
    lasti = i;
    lastj = j;
    northy, eastx, southy, westx = getwallpositions(i, j, paramsx0,paramsy0,paramsx1,paramsy1,paramsdx,paramsdy,paramsnx,paramsny);

    if vx > 0
        if vy > 0
            nexty = southy;
            nextx = eastx;
                else #vy<0
            nexty = northy;
            nextx = eastx;
        end
    else #vx<0
        if vy > 0
            nexty = southy;
            nextx = westx;
    else #vy<0
            nexty = northy;
            nextx = westx;
        end
    end

    # get scaled distance to next walls
    disttowallx = (nextx-x)/vx;
    disttowally = (nexty-y)/vy;

    if disttowallx == 0
        disttowallx = paramsdx/abs(vx);
    end
    if disttowally == 0
        disttowally = paramsdy/abs(vy);
    end

    if disttowallx < disttowally
        l = disttowallx;
    else
        l = disttowally;
    end
    l = l + 1e-8;

    x = x +l*vx;
    y = y +l*vy;
    L = [L,l];

    i, j, valid = getij(x, y, paramsx0,paramsy0,paramsx1,paramsy1,paramsdx,paramsdy,paramsnx,paramsny);
    if valid
        if abs(i-lasti)+abs(j-lastj)==1
        else
            warning("this should not happen!")
            exit()
        end
    end
end

xresult = x;
yresult = y;

return I,J,L,x,y;

end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    SingleBridgeOptimizer{BT<:AbstractBridge}(model::MOI.ModelLike)

Return `AbstractBridgeOptimizer` that always bridges any objective function
supported by the bridge `BT`.

This is in contrast with the [`MathOptInterface.Bridges.LazyBridgeOptimizer`](@ref),
which only bridges the objective function if it is supported by the bridge `BT`
and unsupported by `model`.

## Example

```jldoctest obj_singlebridgeoptimizer; setup=:(using MathOptInterface; const MOI = MathOptInterface)
julia> struct MyNewBridge{T} <: MOI.Bridges.Objective.AbstractBridge end

julia> bridge = MOI.Bridges.Objective.SingleBridgeOptimizer{MyNewBridge{Float64}}(
           MOI.Utilities.Model{Float64}(),
       )
MOIB.Objective.SingleBridgeOptimizer{MyNewBridge{Float64}, MOIU.Model{Float64}}
with 0 objective bridges
with inner model MOIU.Model{Float64}
```

## Implementation notes

All bridges should simplify the creation of `SingleBridgeOptimizer`s by defining
a constant that wraps the bridge in a `SingleBridgeOptimizer`.
```jldoctest obj_singlebridgeoptimizer
julia> const MyNewBridgeModel{T,OT<:MOI.ModelLike} =
           MOI.Bridges.Objective.SingleBridgeOptimizer{MyNewBridge{T},OT};
```
This enables users to create bridged models as follows:
```jldoctest obj_singlebridgeoptimizer
julia> MyNewBridgeModel{Float64}(MOI.Utilities.Model{Float64}())
MOIB.Objective.SingleBridgeOptimizer{MyNewBridge{Float64}, MOIU.Model{Float64}}
with 0 objective bridges
with inner model MOIU.Model{Float64}
```
"""
mutable struct SingleBridgeOptimizer{BT<:AbstractBridge,OT<:MOI.ModelLike} <:
               MOI.Bridges.AbstractBridgeOptimizer
    model::OT
    map::Map # `MOI.ObjectiveFunction` -> objective bridge
end

function SingleBridgeOptimizer{BT}(model::OT) where {BT,OT<:MOI.ModelLike}
    return SingleBridgeOptimizer{BT,OT}(model, Map())
end

bridges(::MOI.Bridges.AbstractBridgeOptimizer) = EmptyMap()

bridges(bridge::SingleBridgeOptimizer) = bridge.map

MOI.Bridges.supports_constraint_bridges(::SingleBridgeOptimizer) = false

function MOI.Bridges.is_bridged(
    ::SingleBridgeOptimizer,
    ::Type{<:MOI.AbstractSet},
)
    return false
end

function MOI.Bridges.is_bridged(
    ::SingleBridgeOptimizer,
    ::Type{<:MOI.AbstractFunction},
    ::Type{<:MOI.AbstractSet},
)
    return false
end

MOI.Bridges.is_bridged(::SingleBridgeOptimizer, ::MOI.VariableIndex) = false

# We need the next two methods to remove method ambiguities.

function MOI.Bridges.is_bridged(
    ::SingleBridgeOptimizer,
    ::MOI.ConstraintIndex{MOI.VariableIndex},
)
    return false
end

function MOI.Bridges.is_bridged(
    ::SingleBridgeOptimizer,
    ::MOI.ConstraintIndex{MOI.VectorOfVariables},
)
    return false
end

function MOI.Bridges.supports_bridging_objective_function(
    ::SingleBridgeOptimizer{BT},
    F::Type{<:MOI.AbstractScalarFunction},
) where {BT}
    return supports_objective_function(BT, F)
end

function MOI.Bridges.is_bridged(
    bridge::SingleBridgeOptimizer,
    F::Type{<:MOI.AbstractScalarFunction},
)
    return MOI.Bridges.supports_bridging_objective_function(bridge, F)
end

function MOI.Bridges.bridge_type(
    ::SingleBridgeOptimizer{BT},
    ::Type{<:MOI.AbstractScalarFunction},
) where {BT}
    return BT
end

MOI.Bridges.recursive_model(bridge::SingleBridgeOptimizer) = bridge.model


text:
#####
##### Abstract Particle and Particle Group
#####

abstract type AbstractParticle end

position(p::AbstractParticle) = p.position

abstract type AbstractParticleGroup end

position(g::AbstractParticleGroup) = g.positions

Base.iterate(g::AbstractParticleGroup, state=1) = state > length(g) ? nothing : (g[state], state+1)
Base.length(g::AbstractParticleGroup) = length(g.positions)
Base.size(g::AbstractParticleGroup) = length(g)

Base.firstindex(g::AbstractParticleGroup) = 1
Base.lastindex(g::AbstractParticleGroup) = length(g)

#####
##### Nucleus and Nuclei
#####

struct Nucleus <: AbstractParticle
    position::Vec3{Float64}
    charge::Float64
    number::Int
end

function Nucleus(
    position::AbstractVector{Float64},
    number::Int,
    )
    Nucleus(position, Float64(number), number)
end

function Nucleus(
    p::AbstractParticle,
    )
    Nucleus(p.position, p.number)
end

struct Nuclei <: AbstractParticleGroup
    positions::Vector{Vec3{Float64}}
    charges::Vector{Float64}
    numbers::Vector{Int}
end

function Nuclei(
    positions::Vector{Vec3{Float64}},
    numbers::Vector{Int}
    )
    return Nuclei(positions, Float64.(numbers), numbers)
end

Nuclei() = Nuclei(Vector{Vec3{Float64}}[], Vector{Float64}[], Vector{Int}[])
Nuclei(file::String) = Nuclei(parse_geom(file)...)

function Base.push!(g::Nuclei, p::Nucleus)
    push!(g.positions, p.position)
    push!(g.charges, p.charge)
    push!(g.numbers, p.number)
end

Base.getindex(g::Nuclei, i::Int) = Nucleus(g.positions[i], g.charges[i], g.numbers[i])

#####
##### Point Charge and Point Charges
#####

struct PointCharge <: AbstractParticle
    position::Vec3{Float64}
    charge::Float64
end

function PointCharge(
    p::AbstractParticle,
    )
    return PointCharge(p.position, p.charge)
end

function PointCharge(
    p::AbstractParticle,
    charge::Float64
    )
    return PointCharge(p.position, charge)
end

struct PointCharges <: AbstractParticleGroup
    positions::Vector{Vec3{Float64}}
    charges::Vector{Float64}
end

PointCharges(g::AbstractParticleGroup) = PointCharges(g.positions, g.charges)
PointCharges() = PointCharges(Vector{Vec3{Float64}}[], Vector{Float64}[])
PointCharges(file::String) = PointCharges(parse_pc(file)...)

function Base.push!(g::PointCharges, p::PointCharge)
    push!(g.positions, p.position)
    push!(g.charges, p.charge)
end

Base.getindex(g::PointCharges, i::Int) = PointCharge(g.positions[i], g.charges[i])

#####
##### Gaussian Charge and Gaussian Charges
#####

struct GaussianCharge <: AbstractParticle
    position::Vec3{Float64}
    charge::Float64
    sigma::Float64
end

function GaussianCharge(
    p::AbstractParticle,
    sigma::Float64,
    )
    return GaussianCharge(p.position, p.charge, sigma)
end

struct GaussianCharges <: AbstractParticleGroup
    positions::Vector{Vec3{Float64}}
    charges::Vector{Float64}
    sigmas::Vector{Float64}
end

GaussianCharges(g::AbstractParticleGroup, sigmas::Vector{Float64}) = GaussianCharges(g.positions, g.charges, sigmas)
GaussianCharges() = GaussianCharges(Vector{Vec3{Float64}}[], Vector{Float64}[], Vector{Float64}[])

function Base.push!(g::GaussianCharges, p::GaussianCharge)
    push!(g.positions, p.position)
    push!(g.charges, p.charge)
    push!(g.sigmas, p.sigma)
end

Base.getindex(g::GaussianCharges, i::Int) = GaussianCharge(g.positions[i], g.charges[i], g.sigmas[i])

#####
##### Point and Points
#####

struct Point <: AbstractParticle
    position::Vec3{Float64}
end

struct Points <: AbstractParticleGroup
    positions::Vector{Vec3{Float64}}
end

Points(g::AbstractParticleGroup) = Points(g.positions)
Points() = Points(Vector{Vec3{Float64}}[])

function Base.push!(g::Points, p::Point)
    push!(g.positions, p.position)
end

Base.getindex(g::Points, i::Int) = Point(g.positions[i])

#####
##### Typed Atom and Typed Atoms
#####

struct TypedAtom <: AbstractParticle
    position::Vec3{Float64}
    charge::Float64
    type::Int
end

struct TypedAtoms <: AbstractParticleGroup
    positions::Vector{Vec3{Float64}}
    charges::Vector{Float64}
    types::Vector{Int}
end

TypedAtoms() = TypedAtom(Vector{Vec3{Float64}}[], Vector{Float64}[], Vector{Int}[])

function Base.push!(g::TypedAtoms, p::TypedAtom)
    push!(g.positions, p.position)
    push!(g.charges, p.charge)
    push!(g.types, p.type)
end

Base.getindex(g::TypedAtoms, i::Int) = TypedAtom(g.positions[i], g.charges[i], g.types[i])


text:
using Base.Test
# morlet1d.jl
import WaveletScattering: gauss

# gauss
@test_approx_eq gauss(0.0, 1.0) 1.0
for ω in 1.0:10.0
    for den in logspace(0, 3, 4)
        g = gauss(ω, den)
        @test g >= 0.0
        @test_approx_eq g gauss(-ω, den)
        @test_approx_eq sqrt(-log(g) * den) ω
    end
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct StaticReserveNonSpinning <: ReserveNonSpinning
        name::String
        available::Bool
        time_frame::Float64
        requirement::Float64
        ext::Dict{String, Any}
        operation_cost::Union{Nothing, TwoPartCost}
        internal::InfrastructureSystemsInternal
    end

Data Structure for a non-spinning reserve product for system simulations.

# Arguments
- `name::String`
- `available::Bool`
- `time_frame::Float64`: the relative saturation time_frame, validation range: `(0, nothing)`, action if invalid: `error`
- `requirement::Float64`: the static value of required reserves in system p.u., validation range: `(0, nothing)`, action if invalid: `error`
- `ext::Dict{String, Any}`
- `operation_cost::Union{Nothing, TwoPartCost}`: Cost for providing reserves  [`TwoPartCost`](@ref)
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct StaticReserveNonSpinning <: ReserveNonSpinning
    name::String
    available::Bool
    "the relative saturation time_frame"
    time_frame::Float64
    "the static value of required reserves in system p.u."
    requirement::Float64
    ext::Dict{String, Any}
    "Cost for providing reserves  [`TwoPartCost`](@ref)"
    operation_cost::Union{Nothing, TwoPartCost}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function StaticReserveNonSpinning(name, available, time_frame, requirement, ext=Dict{String, Any}(), operation_cost=nothing, )
    StaticReserveNonSpinning(name, available, time_frame, requirement, ext, operation_cost, InfrastructureSystemsInternal(), )
end

function StaticReserveNonSpinning(; name, available, time_frame, requirement, ext=Dict{String, Any}(), operation_cost=nothing, internal=InfrastructureSystemsInternal(), )
    StaticReserveNonSpinning(name, available, time_frame, requirement, ext, operation_cost, internal, )
end

# Constructor for demo purposes; non-functional.
function StaticReserveNonSpinning(::Nothing)
    StaticReserveNonSpinning(;
        name="init",
        available=false,
        time_frame=0.0,
        requirement=0.0,
        ext=Dict{String, Any}(),
        operation_cost=TwoPartCost(nothing),
    )
end

InfrastructureSystems.get_name(value::StaticReserveNonSpinning) = value.name
"""Get [`StaticReserveNonSpinning`](@ref) `available`."""
get_available(value::StaticReserveNonSpinning) = value.available
"""Get [`StaticReserveNonSpinning`](@ref) `time_frame`."""
get_time_frame(value::StaticReserveNonSpinning) = value.time_frame
"""Get [`StaticReserveNonSpinning`](@ref) `requirement`."""
get_requirement(value::StaticReserveNonSpinning) = value.requirement
"""Get [`StaticReserveNonSpinning`](@ref) `ext`."""
get_ext(value::StaticReserveNonSpinning) = value.ext
"""Get [`StaticReserveNonSpinning`](@ref) `operation_cost`."""
get_operation_cost(value::StaticReserveNonSpinning) = value.operation_cost
"""Get [`StaticReserveNonSpinning`](@ref) `internal`."""
get_internal(value::StaticReserveNonSpinning) = value.internal

InfrastructureSystems.set_name!(value::StaticReserveNonSpinning, val) = value.name = val
"""Set [`StaticReserveNonSpinning`](@ref) `available`."""
set_available!(value::StaticReserveNonSpinning, val) = value.available = val
"""Set [`StaticReserveNonSpinning`](@ref) `time_frame`."""
set_time_frame!(value::StaticReserveNonSpinning, val) = value.time_frame = val
"""Set [`StaticReserveNonSpinning`](@ref) `requirement`."""
set_requirement!(value::StaticReserveNonSpinning, val) = value.requirement = val
"""Set [`StaticReserveNonSpinning`](@ref) `ext`."""
set_ext!(value::StaticReserveNonSpinning, val) = value.ext = val
"""Set [`StaticReserveNonSpinning`](@ref) `operation_cost`."""
set_operation_cost!(value::StaticReserveNonSpinning, val) = value.operation_cost = val


text:
using DiffEqFlux, OrdinaryDiffEq, Test

u0 = Float32[2.; 0.]
datasize = 30
tspan = (0.0f0,1.5f0)

function trueODEfunc(du,u,p,t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end
t = range(tspan[1],tspan[2],length=datasize)
prob = ODEProblem(trueODEfunc,u0,tspan)
ode_data = Array(solve(prob,Tsit5(),saveat=t))

model = Chain(x -> x.^3,
             Dense(2,50,tanh),
             Dense(50,2))
neuralde = NeuralODE(model,tspan,Rodas5(),saveat=t,reltol=1e-7,abstol=1e-9)

function predict_n_ode()
  neuralde(u0)
end
loss_n_ode() = sum(abs2,ode_data .- predict_n_ode())

data = Iterators.repeated((), 10)
opt = ADAM(0.1)
cb = function () #callback function to observe training
  display(loss_n_ode())
end

# Display the ODE with the initial parameter values.
cb()

neuralde = NeuralODE(model,tspan,Rodas5(),saveat=t,reltol=1e-7,abstol=1e-9)
ps = Flux.params(neuralde)
loss1 = loss_n_ode()
Flux.train!(loss_n_ode, ps, data, opt, cb = cb)
loss2 = loss_n_ode()
@test loss2 < loss1

neuralde = NeuralODE(model,tspan,KenCarp4(),saveat=t,reltol=1e-7,abstol=1e-9)
ps = Flux.params(neuralde)
loss1 = loss_n_ode()
Flux.train!(loss_n_ode, ps, data, opt, cb = cb)
loss2 = loss_n_ode()
@test loss2 < loss1

neuralde = NeuralODE(model,tspan,RadauIIA5(),saveat=t,reltol=1e-7,abstol=1e-9)
ps = Flux.params(neuralde)
loss1 = loss_n_ode()
Flux.train!(loss_n_ode, ps, data, opt, cb = cb)
loss2 = loss_n_ode()
@test loss2 < loss1


text:
using Revise

using Pkg
Pkg.activate(".")
Pkg.instantiate()

using RadMod2D

using GLMakie
# using CairoMakie

include("./plot2D.jl")
include("./tests.jl")

# test_mesh2D()
# test_raycast2D()
# test_view2D_blocking_2elem() # not working
# test_view2D_blocking_bf_vs_tiles()
# test_view2D_vf()
# test_view2D_tiles(n = 30)
test_therm2D()

# include("./test/run_test.jl")


text:
using DocSeeker
using Test

include("dynamicsearch.jl")
include("finddocs.jl")


text:
include("../src/SpecialGraphs.jl")
include("runtests.jl")


text:
struct Midpoints{T <: AbstractFloat} <: AbstractVector{T}
    data::LinRange{T}
    function Midpoints{T}(n) where {T}
        if n == 0
            return new{T}(LinRange{T}(0,0,0))
        else
            return new{T}(LinRange{T}(-1,1,2n+1)[2:2:end-1])
        end
    end
end
Midpoints(n) = Midpoints{Float64}(n)
Base.size(x::Midpoints) = size(x.data)
Base.getindex(x::Midpoints, i::Int) = x.data[i]


text:
type Group
	dim::Dimension
	name::ASCIIString
	reduction::ASCIIString
end
Base.write(io::IO, group::Group) = print(io, "var ", group.name, " = ", group.dim.name, ".group().", group.reduction, ";")

"""
  reduce_count

A reduction for a Group that is simply reduceCount()
"""
reduce_count(dim::Dimension) = Group(dim, string(dim.name)*"_count", "reduceCount()")

"""
  reduce_sum

A reduction for a Group that sums the values.
"""
reduce_sum(dim::Dimension) = Group(dim, string(dim.name)*"_sum", @sprintf("reduceSum(function(d){ return d.%s; })", dim.name))

"""
  reduce_sum

 A master reduction which sums values from all provided columns and tallies a count.
 Useful for making more complex charts like bubble charts.
"""
function reduce_master(dim::Dimension, columns::Vector{Symbol})
  reduction_str = IOBuffer()
  write(reduction_str, "reduce(function (p, v) {
  ++p.DCCount;
")
  for col in columns
    write(reduction_str, "  p.", col, "_sum += v.", col, ";\n")
  end
  write(reduction_str,"  return p;
},
function (p, v) {
  --p.DCCount;
")
  for col in columns
    write(reduction_str, "  p.", col, "_sum -= v.", col, ";\n")
  end
  write(reduction_str, "  return p;
},
function () {
  return {
    DCCount: 0,
")
  for col in columns
    write(reduction_str, "    ", col, "_sum: 0,\n")
  end
  write(reduction_str, "  };
})")
  # Assign a random number to master reduction so it doesn't get overwritten
  Group(dim, string(dim.name, "_master", rand(0:999)), takebuf_string(reduction_str))
end

"""
  infer_group

Infer construction of a group based on the array datatype.
"""
infer_group{I<:Integer}(arr::AbstractDataArray{I}, dim::Dimension) = reduce_sum(dim)
infer_group{F<:AbstractFloat}(arr::AbstractDataArray{F}, dim::Dimension) = reduce_sum(dim)
infer_group{S<:AbstractString}(arr::AbstractDataArray{S}, dim::Dimension) = reduce_count(dim)


text:
"""
EcoTrialStructure is a package for representing behavioral/physiological data from economic choice experiments.

You can learn more about each of the following from its own docstring, e.g., `?CellsTrial`:

- `CellsTrial`: Store cellular responses for a single trial.
- `EventTiming`: Specify the timing of events during a trial.
- `FrameSeq`: A sequence of `nframes` frames starting nearest to a specified time.
- `OfferType`: Encode the offer quantities and side (left/right).
- `TrialResult`: Encode the `OfferType` and the animal's choice.
- `isforced`: Return `true` for a forced-choice trial.
- `iswrong`: Return `true` if the animal chose incorrectly on a forced-choice trial.
- `madechoice`: Return `true` if the animal made a choice in the trial.
- `ncells`: Return the number of cells in a `CellsTrial`
- `parsemat`: Parse the trial data in `filename`.
- `positive_cells`: Extract the vector of `idx_positive_cells` from the .mat file.
"""
module EcoTrialStructure

using Unitful
using Unitful: ms, s
using IntervalSets
using OffsetArrays
using OrderedCollections
using MAT

export CellsTrial, FrameSeq, OfferType, TrialResult, EventTiming
export isdeferred, isforced, iswrong, madechoice, ncells, .., ms, s
export parsemat, positive_cells

include("types.jl")
include("utils.jl")
include("matlab.jl")

Base.@deprecate_binding TrialType OfferType

end


text:
# Descriptor of a cartesian grid

"""
    struct CartesianDescriptor{D,T,F<:Function}
      origin::Point{D,T}
      sizes::NTuple{D,T}
      partition::NTuple{D,Int}
      map::F
    end

Struct that stores the data defining a Cartesian grid.
"""
struct CartesianDescriptor{D,T,F<:Function} <: GridapType
  origin::Point{D,T}
  sizes::NTuple{D,T}
  partition::NTuple{D,Int}
  map::F
  isperiodic::NTuple{D,Bool}
  @doc """
      CartesianDescriptor(
        origin::Point{D},
        sizes::NTuple{D},
        partition;
        map::Function=identity,
        isperiodic::NTuple{D,Bool}=tfill(false,Val{D})) where D

  `partition` is a 1D indexable collection of arbitrary type.
  """
  function CartesianDescriptor(
    origin::Point{D},
    sizes::NTuple{D},
    partition;
    map::Function=identity,
    isperiodic::NTuple{D,Bool}=tfill(false,Val{D}())) where D

    for i in 1:D
      if isperiodic[i]
        @assert (partition[i] > 2) "A minimum of 3 elements is required in any
          periodic direction"
      end
    end

    T = eltype(sizes)
    F = typeof(map)
    new{D,T,F}(origin,sizes,Tuple(partition),map,isperiodic)
  end
end

"""
    CartesianDescriptor(
      domain,
      partition;
      map::Function=identity,
      isperiodic::NTuple{D,Bool}=tfill(false,Val{D}))

`domain` and `partition` are 1D indexable collections of arbitrary type.
"""
function CartesianDescriptor(
  domain,
  partition;
  map::Function=identity,
  isperiodic::NTuple=tfill(false,Val{length(partition)}()))

  D = length(partition)
  limits = [(domain[2*d-1],domain[2*d]) for d in 1:D]
  sizes = Tuple([(limits[d][2]-limits[d][1])/partition[d] for d in 1:D])
  origin = Point([ limits[d][1] for d in 1:D]...)
  CartesianDescriptor(origin,sizes,partition;map=map,isperiodic=isperiodic)
end

"""
    CartesianDescriptor(
      pmin::Point{D},
      pmax::Point{D},
      partition;
      map::Function=identity,
      isperiodic::NTuple{D,Bool}=tfill(false,Val{D})) where D

`partition` is a 1D indexable collection of arbitrary type.
"""
function CartesianDescriptor(
  pmin::Point{D},
  pmax::Point{D},
  partition;
  map::Function=identity,
  isperiodic::NTuple{D,Bool}=tfill(false,Val{D}())) where D

  T = eltype(pmin)
  domain = zeros(T,2*D)
  for d in 1:D
    domain[2*(d-1)+1] = pmin[d]
    domain[2*(d-1)+2] = pmax[d]
  end
  CartesianDescriptor(domain,partition;map=map,isperiodic=isperiodic)
end

# Deprecated signatures for backwards compatibility

@deprecate CartesianDescriptor(origin::Point{D}, sizes::NTuple{D}, partition, map::Function) where D CartesianDescriptor(origin, sizes, partition; map=map)

@deprecate CartesianDescriptor(domain,partition,map::Function) CartesianDescriptor(domain,partition;map=map)

@deprecate CartesianDescriptor(pmin::Point{D}, pmax::Point{D}, partition, map::Function) where D CartesianDescriptor(pmin, pmax, partition; map=map)

# Coordinates

struct CartesianCoordinates{D,T,F} <: AbstractArray{Point{D,T},D}
  data::CartesianDescriptor{D,T,F}
end

Base.size(a::CartesianCoordinates) = a.data.partition .+ 1

Base.IndexStyle(::Type{<:CartesianCoordinates}) = IndexCartesian()

function Base.getindex(a::CartesianCoordinates{D,T}, I::Vararg{Int,D}) where {D,T}
  p = zero(mutable(Point{D,T}))
  x0 = a.data.origin
  dx = a.data.sizes
  @inbounds for d in 1:D
    p[d] =  x0[d] + (I[d]-1)*dx[d]
  end
  a.data.map(Point(p))
end

# Cell nodes

struct CartesianCellNodes{D} <: AbstractArray{Vector{Int},D}
  partition::NTuple{D,Int}
  function CartesianCellNodes(partition)
    D = length(partition)
    new{D}(partition)
  end
end

Base.size(a::CartesianCellNodes) = a.partition

Base.IndexStyle(::Type{<:CartesianCellNodes}) = IndexCartesian()

function array_cache(a::CartesianCellNodes{D}) where D
  zeros(Int,2^D)
end

@inline function getindex!(cache,a::CartesianCellNodes,i::Integer)
  cis = CartesianIndices(size(a))
  ci = cis[i]
  getindex!(cache,a,Tuple(ci)...)
end

@inline function getindex!(v,a::CartesianCellNodes{1},i::Integer)
  v[1] = i
  v[2] = i+1
  v
end

@inline function getindex!(v,a::CartesianCellNodes{D},i::Integer...) where D
  nodes = LinearIndices(size(a).+1)
  lnodes = CartesianIndices(tfill(2,Val{D}()))
  j = i .- 1
  for (i,lnode) in enumerate(lnodes)
    k = Tuple(lnode) .+ j
    v[i] = nodes[ k... ]
  end
  v
end

function Base.getindex(a::CartesianCellNodes,i::Integer...)
  cache = array_cache(a)
  getindex!(cache,a,i...)
end

# Cartesian Grid

"""
    struct CartesianGrid{D,T,F} <: Grid{D,D}
      # private fields
    end
"""
struct CartesianGrid{D,T,F} <: Grid{D,D}
  node_coords::CartesianCoordinates{D,T,F}
  cell_nodes::CartesianCellNodes{D}
  cell_type::Fill{Int8,1,Tuple{Base.OneTo{Int}}}
  @doc """
      CartesianGrid(desc::CartesianDescriptor)
  """
  function CartesianGrid(desc::CartesianDescriptor{D,T,F}) where {D,T,F}
    node_coords = CartesianCoordinates(desc)
    cell_nodes = CartesianCellNodes(desc.partition)
    cell_type = Fill(Int8(1),length(cell_nodes))
    new{D,T,F}(node_coords,cell_nodes,cell_type)
  end
end

OrientationStyle(::Type{<:CartesianGrid}) = Val{true}()

"""
    get_cartesian_descriptor(grid::CartesianGrid)

Get the descriptor of the Cartesian grid
"""
function get_cartesian_descriptor(a::CartesianGrid)
  a.node_coords.data
end

get_node_coordinates(g::CartesianGrid) = g.node_coords

get_cell_type(g::CartesianGrid) = g.cell_type

get_cell_nodes(g::CartesianGrid) = g.cell_nodes

function get_reffes(g::CartesianGrid{D}) where D
  p = Polytope(tfill(HEX_AXIS,Val{D}()))
  order = 1
  reffe = LagrangianRefFE(Float64,p,order)
  [reffe,]
end

get_reffes(g::CartesianGrid{1}) = [SEG2,]

get_reffes(g::CartesianGrid{2}) = [QUAD4,]

get_reffes(g::CartesianGrid{3}) = [HEX8,]

"""
    CartesianGrid(args...;kwargs...)

Same args needed to construct a `CartesianDescriptor`
"""
function CartesianGrid(args...;kwargs...)
  desc = CartesianDescriptor(args...;kwargs...)
  CartesianGrid(desc)
end

# Cell map

struct CartesianMap{D,T,L} <: AbstractArray{AffineMap{D,T,L},D}
  data::CartesianDescriptor{D,T,typeof(identity)}
  function CartesianMap(des::CartesianDescriptor{D,T}) where {D,T}
    L = D*D
    new{D,T,L}(des)
  end
end

Base.size(a::CartesianMap) = a.data.partition

Base.IndexStyle(::Type{<:CartesianMap}) = IndexCartesian()

function Base.getindex(a::CartesianMap{D,T},I::Vararg{Int,D}) where {D,T}
  p = zero(mutable(Point{D,T}))
  x0 = a.data.origin
  dx = a.data.sizes
  @inbounds for d in 1:D
    p[d] =  x0[d] + (I[d]-1)*dx[d]
  end
  origin =  Point(p)
  jacobian = diagonal_tensor(VectorValue(dx))
  AffineMap(jacobian,origin)
end

function field_array_gradient(a::CartesianMap)
  dx = a.data.sizes
  jacobian = diagonal_tensor(VectorValue(dx))
  j = AffineMapGrad(jacobian)
  Fill(j,length(a))
end

function field_array_gradient(a::Reindexed{T,N,A}) where {T,N,A<:CartesianMap}
  g = field_array_gradient(a.i_to_v)
  reindex(g,a.j_to_i)
end

function get_cell_map(grid::CartesianGrid{D,T,typeof(identity)} where {D,T})
  CartesianMap(grid.node_coords.data)
end

# Cartesian grid topology with periodic BC

function _cartesian_grid_topology_with_periodic_bcs(grid::UnstructuredGrid,
  isperiodic::NTuple,
  partition)

  cell_to_vertices, vertex_to_node =
    _generate_cell_to_vertices_from_grid(grid, isperiodic, partition)
  _generate_grid_topology_from_grid(grid,cell_to_vertices,vertex_to_node)
end

function _generate_cell_to_vertices_from_grid(grid::UnstructuredGrid,
  isperiodic::NTuple, partition)

  if is_first_order(grid)
    nodes = get_cell_nodes(grid)
    cell_to_vertices = copy(nodes)

    nnodes = num_nodes(grid)
    num_nodes_x_dir = [partition[i]+1 for i in 1:length(partition)]
    point_to_isperiodic, slave_point_to_point, slave_point_to_master_point =
      _generate_slave_to_master_point(num_nodes_x_dir,isperiodic, nnodes)

    vertex_to_point = findall( .! point_to_isperiodic)
    point_to_vertex = fill(-1,length(point_to_isperiodic))
    point_to_vertex[vertex_to_point] = 1:length(vertex_to_point)
    point_to_vertex[slave_point_to_point] = point_to_vertex[slave_point_to_master_point]

    cell_to_vertices = Table(LocalToGlobalArray(nodes,point_to_vertex))

    vertex_to_node = vertex_to_point
    node_to_vertex = point_to_vertex
  else
    @notimplemented
  end
  (cell_to_vertices,vertex_to_node, node_to_vertex)
end

function _generate_slave_to_master_point(num_nodes_x_dir::Vector{Int},
  isperiodic::NTuple, num_nodes::Int)

  point_to_isperiodic = fill(false,num_nodes)

  slave_ijk_bounds = Array{Any,1}(undef,length(isperiodic))
  master_ijk_bounds = Array{Any,1}(undef,length(isperiodic))

  linear_indices = LinearIndices(Tuple(num_nodes_x_dir))
  periodic_dirs = findall(x->x==true, isperiodic)
  for periodic_dir in periodic_dirs
    for i in 1:length(isperiodic)
      if i == periodic_dir
        slave_ijk_bounds[i] = num_nodes_x_dir[i]
      else
        slave_ijk_bounds[i] = 1:num_nodes_x_dir[i]
      end
    end
    slave_ijk_set = CartesianIndices(Tuple(slave_ijk_bounds))
    point_to_isperiodic[linear_indices[slave_ijk_set]] .= true
  end

  slave_point_to_point = findall( point_to_isperiodic)
  slave_point_to_master_point = Array{Int,1}(undef,length(slave_point_to_point))

  cartesian_indices = CartesianIndices(Tuple(num_nodes_x_dir))
  for (i,point) in enumerate(slave_point_to_point)
    ijk = collect(cartesian_indices[point].I)
    for i in periodic_dirs
      if ijk[i] == num_nodes_x_dir[i]
        ijk[i] = 1
      end
    end

    master_point_ijk = CartesianIndex(Tuple(ijk))
    slave_point_to_master_point[i] = linear_indices[master_point_ijk]
  end

  point_to_isperiodic, slave_point_to_point, slave_point_to_master_point
end


text:
const _SCI_FUN_POOL = Dict{Int, Symbol}()
let
    empty!(_SCI_FUN_POOL)
    ds = 1:10
    for d in ds
        srtf = "%0.$(d)e"
        fun = Symbol("_sci", d)
        @eval begin
            $fun(n) = @sprintf($srtf, n)
        end
        _SCI_FUN_POOL[d] = fun
    end
end

function sci_str(n; d::Int = 1)
    d = clamp(d, 1, 10)
    scifun = getproperty(@__MODULE__, _SCI_FUN_POOL[d])
    scifun(n)
end


text:
mutable struct SparseMatrix{T} <: AbstractMatrix{T}
	obj::Ptr{Cvoid}
end


text:
module SymArrays

export SymArray, contract, contract!, SymArr_ifsym, symgrp_size, symgrps, nsymgrps, storage_type

include("helpers.jl")

include("symarray.jl")

include("contractions.jl")

include("cuda_contractions.jl")

end # module


text:
mutable struct LinkedQueue{T} <: AbstractQueue{T}
    front::LinkNode{T}
    rear::LinkNode{T}
    size::Int
    function LinkedQueue{T}() where T
        front = rear = LinkNode{T}()
        size = 0
        new(front, rear, size)
    end
end

isempty(q::LinkedQueue) = q.size == 0
length(q::LinkedQueue) = q.size
frontvalue(q::LinkedQueue) = isempty(q) ? throw(ArgumentError("queue must be non-empty")) : q.front.next.data

function enqueue!(q::LinkedQueue, item)
    q.rear = q.rear.next = LinkNode{eltype(q)}(item)
    q.size += 1
end

function dequeue!(q::LinkedQueue)
    (q.size == 0) && throw(ArgumentError("queue must be non-empty"))
    item = q.front.next.data
    if q.front.next === q.rear
        q.rear = q.front
    end
    q.front.next = q.front.next.next
    q.size -= 1
    return item
end


text:
include("bike-attributes.jl")

mutable struct Bike

    bike_attributes :: BikeAttributes
    x_position :: Int
    y_position :: Int

end


text:
"""
    UnderlyingMDP(m::POMDP)

Transform `POMDP` `m` into an `MDP` where the states are fully observed.

    UnderlyingMDP(m::MDP)

Return `m`
"""
struct UnderlyingMDP{P <: POMDP, S, A} <: MDP{S, A}
    pomdp::P
end

function UnderlyingMDP(pomdp::POMDP{S, A, O}) where {S,A,O}
    P = typeof(pomdp)
    return UnderlyingMDP{P,S,A}(pomdp)
end

UnderlyingMDP(m::MDP) = m

POMDPs.transition(mdp::UnderlyingMDP{P, S, A}, s::S, a::A) where {P,S,A}= transition(mdp.pomdp, s, a)
POMDPs.initialstate_distribution(mdp::UnderlyingMDP) = initialstate_distribution(mdp.pomdp)
POMDPs.initialstate(mdp::UnderlyingMDP, rng::AbstractRNG) = initialstate(mdp.pomdp, rng)
POMDPs.states(mdp::UnderlyingMDP) = states(mdp.pomdp)
POMDPs.actions(mdp::UnderlyingMDP) = actions(mdp.pomdp)
POMDPs.reward(mdp::UnderlyingMDP{P, S, A}, s::S, a::A) where {P,S,A} = reward(mdp.pomdp, s, a)
POMDPs.reward(mdp::UnderlyingMDP{P, S, A}, s::S, a::A, sp::S) where {P,S,A} = reward(mdp.pomdp, s, a, sp)
POMDPs.isterminal(mdp ::UnderlyingMDP{P, S, A}, s::S) where {P,S,A} = isterminal(mdp.pomdp, s)
POMDPs.discount(mdp::UnderlyingMDP) = discount(mdp.pomdp)
POMDPs.stateindex(mdp::UnderlyingMDP{P, S, A}, s::S) where {P,S,A} = stateindex(mdp.pomdp, s)
POMDPs.stateindex(mdp::UnderlyingMDP{P, Int, A}, s::Int) where {P,A} = stateindex(mdp.pomdp, s) # fix ambiguity with src/convenience
POMDPs.stateindex(mdp::UnderlyingMDP{P, Bool, A}, s::Bool) where {P,A} = stateindex(mdp.pomdp, s)
POMDPs.actionindex(mdp::UnderlyingMDP{P, S, A}, a::A) where {P,S,A} = actionindex(mdp.pomdp, a)
POMDPs.actionindex(mdp::UnderlyingMDP{P,S, Int}, a::Int) where {P,S} = actionindex(mdp.pomdp, a)
POMDPs.actionindex(mdp::UnderlyingMDP{P,S, Bool}, a::Bool) where {P,S} = actionindex(mdp.pomdp, a)

POMDPs.gen(d::DDNOut, mdp::UnderlyingMDP, s, a, rng) = gen(d, mdp.pomdp, s, a, rng)
POMDPs.gen(d::DDNNode, mdp::UnderlyingMDP, s, a, rng) = gen(d, mdp.pomdp, s, a, rng)
POMDPs.gen(mdp::UnderlyingMDP, s, a, rng) = gen(m.pomdp, s, a, rng)

# deprecated in POMDPs v0.8
POMDPs.n_actions(mdp::UnderlyingMDP) = n_actions(mdp.pomdp)
POMDPs.n_states(mdp::UnderlyingMDP) = n_states(mdp.pomdp)
POMDPs.generate_s(mdp::UnderlyingMDP, s, a, rng::AbstractRNG) = generate_s(mdp.pomdp, s, a, rng)
POMDPs.generate_sr(mdp::UnderlyingMDP, s, a, rng::AbstractRNG) = generate_sr(mdp.pomdp, s, a, rng)


text:
n = 12345352345634
my_string = "Hello, I am Srujan"

# Integer to string
str = string(n)
println("string of Integer = $str")

# String to Integer
int = parse(Int64, str)
println("Integer of the string = $int")

# Reverse a string
reversed = reverse(string(n))
println("Reversed string = $reversed")

# Split a string
split_string = split(my_string, " ")
println("Split string = $split_string")

# Join into a single string
joined_string = join(split_string, "-")
println("Joined string = $joined_string")


text:
using OrdinaryDiffEq, StaticArrays, BenchmarkTools, Sundials, LSODA

println("Lorenz")
# Lorenz
function lorenz_static(u,p,t)
    @inbounds begin
        dx = p[1]*(u[2]-u[1])
        dy = u[1]*(p[2]-u[3]) - u[2]
        dz = u[1]*u[2] - p[3]*u[3]
    end
    SA[dx,dy,dz]
end
u0 = SA[1.0,0.0,0.0]
p  = SA[10.0,28.0,8/3]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz_static,u0,tspan,p)
@btime solve(prob,Tsit5(),saveat=0.1,reltol=1.0e-8,abstol=1.0e-8) # 1.340 ms (27 allocations: 64.55 KiB)
@btime solve(prob,Vern8(),saveat=0.1,reltol=1.0e-8,abstol=1.0e-8); # 718.000 μs (28 allocations: 71.52 KiB)

println("\nRober")
# rober
function rober(du,u,p,t)
    y₁,y₂,y₃ = u
    k₁,k₂,k₃ = p
    @inbounds begin
        du[1] = -k₁*y₁+k₃*y₂*y₃
        du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
        du[3] =  k₂*y₂^2
    end
    nothing
end
prob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e5),[0.04,3e7,1e4])

@btime solve(prob,Rodas5(),reltol=1.0e-8,abstol=1.0e-8, saveat = 1000) # 227.200 μs (1288 allocations: 52.56 KiB)
@btime solve(prob,TRBDF2(),reltol=1.0e-8,abstol=1.0e-8, saveat = 1000) # 2.202 ms (3424 allocations: 194.41 KiB)
@btime solve(prob,CVODE_BDF(),reltol=1.0e-8,abstol=1.0e-8, saveat = 1000) # 627.500 μs (7183 allocations: 214.33 KiB)
@btime solve(prob,lsoda(),reltol=1.0e-8,abstol=1.0e-8, saveat = 1000) # 254.400 μs (2370 allocations: 127.09 KiB)

# rober
function rober_static(u,p,t)
    y₁,y₂,y₃ = u
    k₁,k₂,k₃ = p
    du1 = -k₁*y₁+k₃*y₂*y₃
    du2 =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
    du3 =  k₂*y₂^2
    SA[du1,du2,du3]
end
prob = ODEProblem{false}(rober_static,SA[1.0,0.0,0.0],(0.0,1e5),SA[0.04,3e7,1e4])
@btime solve(prob,Rodas5(),reltol=1.0e-8,abstol=1.0e-8, saveat = 1000) # 42.300 μs (29 allocations: 9.75 KiB)

function rober_jac(u,p,t)
  y₁,y₂,y₃ = u
  k₁,k₂,k₃ = p
  J11 = k₁ * -1
  J21 = k₁
  J31 = 0
  J12 = y₃ * k₃
  J22 = y₂ * k₂ * -2 + y₃ * k₃ * -1
  J32 = y₂ * 2 * k₂
  J13 = k₃ * y₂
  J23 = k₃ * y₂ * -1
  J33 = 0
  SA[J11 J12 J13
     J21 J22 J23
     J31 J32 J33]
end

ff = ODEFunction(rober_static, jac=rober_jac)
prob2 = ODEProblem{false}(ff,SA[1.0,0.0,0.0],(0.0,1e5),SA[0.04,3e7,1e4])
@btime solve(prob2,Rodas5(), reltol=1.0e-8, abstol=1.0e-8, saveat = 1000); # 42.000 μs (30 allocations: 9.70 KiB)

max_stars_repo_name: Nicholaswogan/NumbaLSODA

max_issues_repo_name: Nicholaswogan/NumbaLSODA

max_forks_repo_name: Nicholaswogan/NumbaLSODA


text:
module MusicologyPlots

using DigitalMusicology
using VegaLite
using DataFrames

export pianoroll, plotpolygrams

function pianoroll(notes; color=nothing, width=800, height=400)
    notedf = DataFrame(pitch=map(note -> convert(Int, pitch(note)), notes),
                       onset=map(onset, notes),
                       offset=map(offset, notes))
    maxp = maximum(notedf[:pitch])
    minp = minimum(notedf[:pitch])
    if color==nothing
        notedf |>
            @vlplot(mark={typ=:rect,clip=true},
                    x = :onset,
                    x2 = :offset,
                    y = {
                        field=:pitch,
                        typ="ordinal",
                        scale={domain=collect(maxp:-1:minp)}
                    },
                    tooltip=:pitch,
                    selection={grid={typ=:interval, bind=:scales}},
                    width=width, height=height)
    else
        notedf |>
            @vlplot(mark={typ=:rect,clip=true},
                    x = :onset,
                    x2 = :offset,
                    y = {
                        field=:pitch,
                        typ="ordinal",
                        scale={domain=collect(maxp:-1:minp)}
                    },
                    color={value=color},
                    tooltip=:pitch,
                    selection={grid={typ=:interval, bind=:scales}},
                    width=wigth, height=height)
    end
end

flatten(v) = vcat(v...)

function plotpolygrams(notes, polys;
                       margin=1, width=800, height=400,
                       start=nothing, stop=nothing)
    #lower = minimum(onset(pgram[1][1]) for pgram in polygrams)
    #upper = maximum(maximum(map(offset, pgram[end])) for pgram in polygrams)
    #inrange(note) = onset(note) > lower-margin && offset(note) < upper+margin
    #bg = pianoroll(notes, color="#888888")
    flattened = map(flatten, polys)
    polynames = map(1:length(flattened)) do i
        fill(string("poly ", i), length(flattened[i]))
    end
    polynotes = flatten(flattened)
    allnotes = vcat(notes, polynotes)
    polydf = DataFrame(pitch=map(note -> convert(Int, pitch(note)), allnotes),
                       onset=map(onset, allnotes),
                       offset=map(offset, allnotes),
                       name=vcat(fill("notes", length(notes)), flatten(polynames)))

    maxp = maximum(polydf[:pitch])
    minp = minimum(polydf[:pitch])
    if start == nothing
        start = minimum(map(onset,polynotes)) - margin
    end
    if stop == nothing
        stop = maximum(map(onset,polynotes)) + margin
    end
    dom = [start, stop]

    polydf |>
        @vlplot(mark={typ=:rect,clip=true},
                x={field=:onset, scale={domain=dom}},
                x2=:offset,
                y={field=:pitch, typ="ordinal", scale={domain=collect(maxp:-1:minp)}},
                color={field=:name, scale={scheme="tableau20"}}, tooltip=:pitch,
                #selection={grid={typ=:interval, bind=:scales}},
                width=width, height=height)
end

include("widgets.jl")

end # module


text:
module TopologyPreprocessing
    # MatrixOrganization.jl
    export matrix_poling,
            subsample_matrix,
            add_random_patch

    # MatrixProcessing.jl
    export shift_to_non_negative,
            normalize_to_01,
            diagonal_symmetrize,
            group_distances,
            generate_indices,
            reduce_arrs_to_min_len,
            increase_arrs_to_max_len,
            get_ordered_matrix,
            group_distances,
            generate_indices,
            arr_to_vec,
            cartesianInd_to_vec,
            sort_indices_by_values,
            set_values!

    # BettiCurves.jl
    export get_bettis,
            normalise_bettis,
            get_vectorized_bettis,
            plot_bettis,
            get_bettis_color_palete

    # BarCodes.jl
    export get_barcodes,
            plot_barcodes,
            plot_barcodes!,
            get_birth_death_ratio,
            get_barcode_lifetime,
            get_barcode_max_lifetime,
            boxplot_birth_death,
            boxplot_lifetime,
            get_barcode_max_db_ratios

    include("MatrixOrganization.jl")
    include("MatrixProcessing.jl")
    include("BettiCurves.jl")
    include("Barcodes.jl")

end # module


text:
@testset "types.jl" begin

    using LinearAlgebra

    p1 = Point(80.0, 190.0)
    p2 = Point(40.0, -10.0)

    @test p1.ϕ ≈ 80.0 atol = 0.1
    @test p1.λ ≈ 190.0 atol = 0.1
    @test (p1-p2).ϕ ≈ 40.0 atol = 0.1
    @test (p1-p2).λ ≈ 200.0 atol = 0.1

    @test rad2deg(deg2rad(p1)).ϕ ≈ 80.0 atol = 0.1
    @test rad2deg(deg2rad(p1)).λ ≈ 190.0 atol = 0.1

    ls = Arc(p1, p2)
    @test ls.point₁.ϕ ≈ 80.0 atol = 0.1
    @test ls.point₂.λ ≈ -10.0 atol = 0.1

    l = Line(p1, 90.0)
    @test l.point.ϕ ≈ 80.0 atol = 0.1
    @test l.azimuth ≈ 90.0 atol = 0.1

    p3 = Point(60.0, 30.0)
    lss1 = Arcs([p1, p2, p3])

    @test lss1.points[1].ϕ ≈ 80.0 atol = 0.1
    @test lss1.points[3].λ ≈ 30.0 atol = 0.1

    @test_throws BoundsError lss2 = Arcs([p1])

    pg = Polygon(Point(60.0, 31.0), [p1, p2, p3])

    @test pg.points[1].ϕ ≈ 80.0 atol = 0.1
    @test length(pg.points) == 4
    @test pg.points[4].λ ≈ 190.0 atol = 0.1

    @test_throws BoundsError pg2 = Polygon(Point(60.0, 31.0), [p1, p2])

    pa1 = Point(50.0, 5.0)
    pa2 = Point(60.0, 5.0)
    pa3 = Point(60.0, 15.0)
    pa4 = Point(50.0, 15.0)

    pg_simple = Polygon(Point(55.0, 6.0), [pa1, pa2, pa3, pa4])
    pg_complex = Polygon(Point(55.0, 6.0), [pa1, pa2, pa4, pa3])

    @test issimple(pg_simple)
    @test iscomplex(pg_complex)

    #concave polygon
    pc1 = Point(64.0, 27.0)
    pc2 = Point(47.0, 0.0)
    pc3 = Point(42.0, 25.0)
    pc4 = intermediate_point(pc2,pc3,0.33)
    pc5 = Point(52.0, 16.0)
    pc6 = intermediate_point(pc2,pc3,0.67)
    pc7 = Point(55.0, -7.0)
    pc8 = intermediate_point(pc2,pc1,0.33)
    pc9 = Point(62.0, 4.0)
    pc10 = intermediate_point(pc2,pc1,0.67)
    pc11 = Point(56.0, 21.0)

    polygon_concave = Polygon(pc11,
    [pc1,pc3,pc6,pc5,pc4,pc2,pc7,pc8,pc9,pc10])

    @test boundingbox(pg_simple)[1][2] ≈ 60.0945 atol = 0.0001
    arcs_simple = Arcs(pg_simple.points)
    @test boundingbox(arcs_simple)[1][2] ≈ 60.0945 atol = 0.0001

    @test isconcave(pg_simple) == false
    @test isconcave(polygon_concave)
    @test isconvex(pg_simple)
    @test isconvex(polygon_concave) == false
    @test isconcave(pg_complex) == false
    @test isconvex(pg_complex) == false

    #concave polygon
    pc1 = Point(64.0, 27.0)
    pc2 = Point(42.0, 25.0)
    pc6 = Point(47.0, 0.0)
    pc3 = intermediate_point(pc2,pc6,0.33)
    # pc3 = Point(46.0, 17.0)
    pc4 = Point(52.0, 16.0)
    pc5 = intermediate_point(pc2,pc6,0.67)
    # pc5 = Point(47.0, 9.0)

    pc7 = Point(55.0, -7.0)
    pc8 = intermediate_point(pc6,pc1,0.33)
    pc9 = Point(62.0, 4.0)
    pc10 = intermediate_point(pc6,pc1,0.67)
    pc11 = Point(56.0, 21.0)

    polygon_concave = Polygon(pc11,
    [pc1,pc2,pc3,pc4,pc5,pc6,pc7,pc8,pc9,pc10])

    polygon_convexhull = Polygon(pc11,
    [pc1,pc2,pc6,pc7,pc9])

    @test length(convexhull(polygon_concave).points) == 6
    @test area(polygon_convexhull) ≈ area(convexhull(polygon_concave))
    @test area(polygon_convexhull) ≈ area(convexhull(polygon_convexhull))

    # pa1 = Point(50.0, 5.0)
    # pa2 = Point(60.0, 5.0)
    # pa3 = Point(60.0, 15.0)
    # pa4 = Point(50.0, 15.0)

    # pg_simple = Polygon(Point(55.0, 6.0), [pa1, pa2, pa3, pa4])
    @test area(pg_simple) ≈ area(convexhull(pg_simple))
end


text:
module VonMisesMixedPlasticity
    mktempdir() do dir
        cd(dir) do
            include(joinpath(@__DIR__, "../docs/src/literate/von_mises_plasticity_mixed_hardening.jl"))
        end
    end
end


text:
f0 = Complex(0.0, 0.0)
f1 = Complex(1.0, 0.0)
f1i = Complex(0.0, 1.0)
f11 = Complex(1.0, 1.0)
f2h = Complex(2.0, 0.5)

d0 = Complex(Double64(0), Double64(0))
d1 = Complex(Double64(1), Double64(0))
d1i = Complex(Double64(0), Double64(1))
d11 = Complex(Double64(1), Double64(1))
d2h = Complex(Double64(2), Double64(0.5))

@testset "Complex" begin
  @test isapprox(log(f11), Complex{Float64}(log(d11)))
  @test isapprox(exp(f11), Complex{Float64}(exp(d11)))

  @test isapprox(sin(f11), Complex{Float64}(sin(d11)))
  @test isapprox(cos(f11), Complex{Float64}(cos(d11)))
  @test isapprox(tan(f11), Complex{Float64}(tan(d11)))
  @test isapprox(csc(f11), Complex{Float64}(csc(d11)))
  @test isapprox(sec(f11), Complex{Float64}(sec(d11)))
  @test isapprox(cot(f11), Complex{Float64}(cot(d11)))

  @test isapprox(sinh(f11), Complex{Float64}(sinh(d11)))
  @test isapprox(cosh(f11), Complex{Float64}(cosh(d11)))
  @test isapprox(tanh(f11), Complex{Float64}(tanh(d11)))
  @test isapprox(csch(f11), Complex{Float64}(csch(d11)))
  @test isapprox(sech(f11), Complex{Float64}(sech(d11)))
  @test isapprox(coth(f11), Complex{Float64}(coth(d11)))

  @test isapprox(asin(f11), Complex{Float64}(asin(d11)))
  @test isapprox(acos(f11), Complex{Float64}(acos(d11)))
  @test isapprox(atan(f2h), Complex{Float64}(atan(d2h)))
  @test isapprox(acsc(f11), Complex{Float64}(acsc(d11)))
  @test isapprox(asec(f11), Complex{Float64}(asec(d11)))
  @test isapprox(acot(f2h), Complex{Float64}(acot(d2h)))

  @test isapprox(asinh(f11), Complex{Float64}(asinh(d11)))
  @test isapprox(acosh(f11), Complex{Float64}(acosh(d11)))
  @test isapprox(atanh(f11), Complex{Float64}(atanh(d11)))
  @test isapprox(acsch(f11), Complex{Float64}(acsch(d11)))
  @test isapprox(asech(f11), Complex{Float64}(asech(d11)))
  @test isapprox(acoth(f11), Complex{Float64}(acoth(d11)))
end

@testset "Complex Accuracy" begin
  x = ComplexDF64(sqrt(df64"2"), cbrt(df64"3"))
  y = Complex{BigFloat}(x)

  @test abs(sin(y) - sin(x)) < 1.0e-31
  @test abs(cos(y) - cos(x)) < 1.0e-31
  @test abs(tan(y) - tan(x)) < 1.0e-31
  @test abs(csc(y) - csc(x)) < 1.0e-31
  @test abs(sec(y) - sec(x)) < 1.0e-31
  @test abs(tan(y) - tan(x)) < 1.0e-31

  @test abs(sinh(y) - sinh(x)) < 1.0e-31
  @test abs(cosh(y) - cosh(x)) < 1.0e-31
  @test abs(tanh(y) - tanh(x)) < 1.0e-31
  @test abs(csch(y) - csch(x)) < 1.0e-31
  @test abs(sech(y) - sech(x)) < 1.0e-31
  @test abs(tanh(y) - tanh(x)) < 1.0e-31

  @test abs(asin(y) - asin(x)) < 1.0e-31
  @test abs(acos(y) - acos(x)) < 1.0e-31
  @test abs(atan(y) - atan(x)) < 1.0e-31
  @test abs(acsc(y) - acsc(x)) < 1.0e-31
  @test abs(asec(y) - asec(x)) < 1.0e-31
  @test abs(atan(y) - atan(x)) < 1.0e-31

  @test abs(asinh(y) - asinh(x)) < 1.0e-31
  @test abs(acosh(y) - acosh(x)) < 1.0e-31
  @test abs(atanh(y) - atanh(x)) < 1.0e-31
  @test abs(acsch(y) - acsch(x)) < 1.0e-31
  @test abs(asech(y) - asech(x)) < 1.0e-31
  @test abs(atanh(y) - atanh(x)) < 1.0e-31

end

@testset "Complex exp(huge)" begin
    x = Double64(-1.168869267851442e79) + Double64(4.140780141319954e87)*im
    @test ComplexDF64(exp(x)) == exp(ComplexDF64(x))
end

@testset "Complex f(Inf)" begin
   @test isapprox(exp(ComplexDF64(Inf)), exp(ComplexF64(Inf)))
   @test isapprox(exp(ComplexDF64(-Inf)), exp(ComplexF64(-Inf)))
   @test isapprox(log(ComplexDF64(Inf)), log(ComplexF64(Inf)))
   @test isapprox(log(ComplexDF64(-Inf)), log(ComplexF64(-Inf)))

   @test isapprox(sinh(ComplexDF64(Inf)), sinh(ComplexF64(Inf)))
   @test isapprox(sinh(ComplexDF64(-Inf)), sinh(ComplexF64(-Inf)))
   @test isapprox(cosh(ComplexDF64(Inf)), cosh(ComplexF64(Inf)))
   @test isapprox(cosh(ComplexDF64(-Inf)), cosh(ComplexF64(-Inf)))
   @test isapprox(tanh(ComplexDF64(Inf)), tanh(ComplexF64(Inf)))
   @test isapprox(tanh(ComplexDF64(-Inf)), tanh(ComplexF64(-Inf)))
   @test isapprox(csch(ComplexDF64(Inf)), csch(ComplexF64(Inf)))
   @test isapprox(csch(ComplexDF64(-Inf)), csch(ComplexF64(-Inf)))
   @test isapprox(sech(ComplexDF64(Inf)), sech(ComplexF64(Inf)))
   @test isapprox(sech(ComplexDF64(-Inf)), sech(ComplexF64(-Inf)))
   @test isapprox(coth(ComplexDF64(Inf)), coth(ComplexF64(Inf)))
   @test isapprox(coth(ComplexDF64(-Inf)), coth(ComplexF64(-Inf)))

   @test isapprox(asin(ComplexDF64(Inf)), asin(ComplexF64(Inf)))
   @test isapprox(asin(ComplexDF64(-Inf)), asin(ComplexF64(-Inf)))
   @test isapprox(acos(ComplexDF64(Inf)), acos(ComplexF64(Inf)))
   @test isapprox(acos(ComplexDF64(-Inf)), acos(ComplexF64(-Inf)))
   @test isapprox(atan(ComplexDF64(Inf)), atan(ComplexF64(Inf)))
   @test isapprox(atan(ComplexDF64(-Inf)), atan(ComplexF64(-Inf)))
   @test isapprox(acsc(ComplexDF64(Inf)), acsc(ComplexF64(Inf)))
   @test isapprox(acsc(ComplexDF64(-Inf)), acsc(ComplexF64(-Inf)))
   @test isapprox(asec(ComplexDF64(Inf)), asec(ComplexF64(Inf)))
   @test isapprox(asec(ComplexDF64(-Inf)), asec(ComplexF64(-Inf)))
   @test isapprox(acot(ComplexDF64(Inf)), acot(ComplexF64(Inf)))
   @test isapprox(acot(ComplexDF64(-Inf)), acot(ComplexF64(-Inf)))

   @test isapprox(asinh(ComplexDF64(Inf)), asinh(ComplexF64(Inf)))
   @test isapprox(asinh(ComplexDF64(-Inf)), asinh(ComplexF64(-Inf)))
   @test isapprox(acosh(ComplexDF64(Inf)), acosh(ComplexF64(Inf)))
   @test isapprox(acosh(ComplexDF64(-Inf)), acosh(ComplexF64(-Inf)))
   @test isapprox(atanh(ComplexDF64(Inf)), atanh(ComplexF64(Inf)))
   @test isapprox(atanh(ComplexDF64(-Inf)), atanh(ComplexF64(-Inf)))
   @test isapprox(acsch(ComplexDF64(Inf)), acsch(ComplexF64(Inf)))
   @test isapprox(acsch(ComplexDF64(-Inf)), acsch(ComplexF64(-Inf)))
  # @test isapprox(asech(ComplexDF64(Inf)), asech(ComplexF64(Inf)))     needs fix in Julia v1
  # @test isapprox(asech(ComplexDF64(-Inf)), asech(ComplexF64(-Inf)))
   @test isapprox(acoth(ComplexDF64(Inf)), acoth(ComplexF64(Inf)))
   @test isapprox(acoth(ComplexDF64(-Inf)), acoth(ComplexF64(-Inf)))
end


text:
using PyPlot #: quiver, gcf, figure, imshow, subplots, title, xlabel, ylabel, gca
using Printf, Colors, Images

# set plot size defaults
SMALL = 8
MEDIUM = 10
LARGE = 12

rc("font", size=SMALL)          # controls default text sizes
rc("axes", titlesize=LARGE)     # fontsize of the axes title
rc("axes", labelsize=SMALL)     # fontsize of the x and y labels
rc("xtick", labelsize=SMALL)    # fontsize of the tick labels
rc("ytick", labelsize=SMALL)    # fontsize of the tick labels
rc("legend", fontsize=SMALL)    # legend fontsize
rc("figure", titlesize=LARGE)   # fontsize of the figure title
rc("figure", figsize=(6,6))
rc("figure", dpi=400)
rc("savefig", dpi=400)
rc("image", origin="lower")
# rcdefaults()

"""
    viewimg(img)

View image as grayscale in Juno.
"""
viewimg(img) = colorview(Gray, img)

"""
    addpoints(inds; ret::Symbol=:figure, labels=[])

Add points to a plot at the locations `inds` optionaly provide `labels`.
"""
function addpoints(inds::Array{CartesianIndex{2},1}; ret::Symbol=:figure, labels=[])
    pos_x = [ind[1] for ind in inds]
    pos_y = [ind[2] for ind in inds]

    ax = PyPlot.scatter(pos_y, pos_x, marker = :x)

    for (k, label) in enumerate(labels)
        annotate(label, (pos_y[k], pos_x[k]))
    end

    if ret == :figure
        return gcf()
    elseif ret == :pyobject
        return ax
    end
end

function addpoints(locations::Matrix; ret::Symbol=:figure, labels=[])
    pos_x = map(x -> locations[x, 1], 1:size(locations, 1))
    pos_y = map(x -> locations[x, 2], 1:size(locations, 1))
    ax = PyPlot.scatter(pos_y, pos_x, marker = :x)

    for (k, label) in enumerate(labels)
        annotate(label, (pos_y[k], pos_x[k]))
    end

    if ret == :figure
        return gcf()
    elseif ret == :pyobject
        return ax
    end
end

"""
    imgoverlay_v2(orange_img::Image, blue_img::Image)

Show an overlay of the images `orange_img` and `blue_img`, each in the specified colors, to see the differences between them.

See also: [`imgshow`](@ref), [`imgoverlay`](@ref)
"""
function imgoverlay_v2(orange_img, blue_img; flipped = false)

    normed_orange_img = orange_img.*(1/extrema(orange_img)[2])
    normed_blue_img = blue_img.*(1/extrema(blue_img)[2])

    if flipped
        normed_orange_img = reverse(normed_orange_img, dims=1)
        normed_blue_img = reverse(normed_blue_img, dims=1)
    end

    color1 = [1 0.7 0.5]
    color2 = [0 0.3 0.5]

    a = normed_orange_img
    b = color1[2].*normed_orange_img .+ color2[2].*normed_blue_img
    c = color1[3].*normed_orange_img .+ color2[3].*normed_blue_img

    imm = collect(colorview(RGB, a, b, c))
    return imm
end

"""
    imgoverlay(img1, img2; fig=nothing, figtitle::String="Image overlay", ret::Symbol=:figure)

Show an overlay of the images `img1` and `img2`, each in a different colors, to better see the differences between them.

See also: [`imgshow`](@ref), [`imgoverlay_v2`](@ref)
"""
function imgoverlay(img1, img2; fig=nothing, figtitle::String="Image overlay", ret::Symbol=:figure, origin_bot=false)

    if fig == nothing
        fig, ax = subplots(dpi = 400)
    else
        fig = fig
        ax = gca()
    end

    img1 = Float64.(img1)
    img2 = Float64.(img2)
    imshow(img1, cmap = :Blues_r, alpha = 0.5)
    imshow(img2, cmap = :Oranges_r, alpha = 0.5)

    if origin_bot == true
        ax.invert_yaxis()
    end

    title(figtitle)

    if ret == :figure
        return gcf()
    elseif ret == :pyobject
        return ax
    end
end

"""
    imgshow(img; <keyword arguments>)

Return a figure with image `img`. Plot origin is in the botom left.

# Keyword Arguments
- `fig=nothing`: add a figure to plot in. By defaults creates a blank new figure.
- `figtitle::String="Image"`: add title to the figure.
- `ret::Symbol=:figure`: set return object, by default returns Figure, other options: :pyobject returns a PyObject. (Using figure makes Juno directly plot.)
- `origin_bot::Bool=false`: set the origin to the lower left corner, which is not typical for images (they will appear flipped).

See also: [`showflow`](@ref), [`imgshowflow`](@ref), [`warp_imgshowflow`](@ref)
"""
function imgshow(img;
                 fig=nothing,
                 figtitle::String="Image",
                 ret::Symbol=:figure,
                 origin_bot::Bool=false)

    if fig == nothing
        fig, ax = subplots(dpi = 400)
    else
        fig = fig
        ax = gca()
    end

    img = Float64.(img)

    imshow(img, cmap = :gray);
    # ax.set_ylim(0, size(img)[1]-1);
    # ax.set_xlim(0, size(img)[2]-1);
    if origin_bot == true
        ax.invert_yaxis()
    end

    xlabel("x")
    ylabel("y")
    title(figtitle)

    if ret == :figure
        return gcf()
    elseif ret == :pyobject
        return ax
    end
end

"""
    showflow(flow::Flow; <keyword arguments>)

Return a figure with the displacement field `flow` by default skipping some vectors to make it easy to read.

# Keyword Arguments
- `flow::Flow`: the vector flow to be plotted.
- `disp_type::Symbol=:full` : display mode, either `:full` -> display all the vectors, or `:sparse` -> displays only the vectors above a threshold, or `:auto` which decide for you based on the data.
- `skip_count=nothing`: the number of vectors to skip between each displayed vector. By default set so that the output is ``20 × 20`` vectors.
- `fig=nothing`: add a figure to plot in. By defaults creates a blank new figure.
- `mag::Real=1`: magnify the plotted vectors.
- `key::Bool=true`: add key with maximum vector length.
- `figtitle::String="Flow"`: add title to the figure.
- `ret::Symbol=:figure`: set return object, by default returns Figure, other options: :pyobject returns a PyObject. (Using figure makes Juno directly plot.)

See also: [`imgshow`](@ref), [`imgshowflow`](@ref), [`warp_imgshowflow`](@ref)
"""
function showflow(flow::Flow; disp_type::Symbol=:auto, skip_count=nothing, fig=nothing, mag::Real=1, key::Bool=true, figtitle::String="Flow", ret::Symbol=:figure)

    # set defaults
    vecs_in_one_dim = 20
    if skip_count == nothing
        skip_count = ceil(Int64, maximum(size(flow))/vecs_in_one_dim)-1
    end

    # prepare data
    uv_flow = zeros(size(flow)..., 2)
    uv_flow[:, :, 1] = real(flow)
    uv_flow[:, :, 2] = imag(flow)

    n = skip_count+1
    # if n == 1
    #     n = 2
    # end

    # meshgrid(x, y) = (repeat(x, outer=length(y)), repeat(y, inner=length(x)))
    # The x coordinates of the arrow locations
    X = [1:size(flow, 2);]

    # The y coordinates of the arrow locations
    Y = [1:size(flow, 1);]

    min_threshold = 1e-10

    if disp_type == :auto
        if count(x -> vec_len(x) > min_threshold, flow)/length(flow) <= 0.1
            disp_type = :sparse
        else
            disp_type = :full
        end
    end

    if disp_type == :full
        trimmed_real = fill(NaN, size(uv_flow[:,:,1]))
        trimmed_real[1:n:end, 1:n:end] .= uv_flow[1:n:end, 1:n:end, 1]
        trimmed_imag = fill(NaN, size(uv_flow[:,:,2]))
        trimmed_imag[1:n:end, 1:n:end] .= uv_flow[1:n:end, 1:n:end, 2]
    elseif disp_type == :sparse
        cond = ((abs.(real(flow)) .> min_threshold) .| (abs.(imag(flow)) .> min_threshold))
        trimmed_real = fill(NaN, size(uv_flow[:,:,1]))
        trimmed_real[cond] .= uv_flow[:, :, 1][cond]
        trimmed_imag = fill(NaN, size(uv_flow[:,:,2]))
        trimmed_imag[cond] .= uv_flow[:, :, 2][cond]
        if all(isnan, trimmed_real) || all(isnan, trimmed_imag)
            trimmed_imag = zeros(size(trimmed_imag))
            trimmed_real = zeros(size(trimmed_real))
        end
    end

    maxi = maximum(filter(!isnan, abs.([trimmed_imag trimmed_real])))

    scaling = maxi/(n * mag)
    # color = atan.(trimmed_real, trimmed_imag)
    # color = (color .- minimum(color)) ./ maximum(color)

    if fig == nothing
        fig, ax = subplots(dpi = 400)
    else
        fig = fig
        ax = gca()
    end

    q = ax.quiver(X, Y, trimmed_real, trimmed_imag, atan.(trimmed_real, trimmed_imag),
            scale_units = "xy", scale = scaling)

    # angles='uv' sets the angle of vector by atan2(u,v), angles='xy' draws the vector from (x,y) to (x+u, y+v)
    ax.set_aspect(1.)
    xlabel("x\n\n")
    ylabel("y")
    title(figtitle)

    # ax.invert_yaxis()

    if key == true
        # subplots_adjust(bottom=0.1)
        max_mag = max_displacement(flow)
        label_text = "length = " * @sprintf("%0.2f", max_mag)
        ax.quiverkey(q, X=0.75, Y = 0.2, U = max_mag, coordinates="inches",
        label=label_text, labelpos = "E")
    end

    if ret == :figure
        return gcf()
    elseif ret == :pyobject
        return q
    end
end

"""
    imgshowflow(img, flow; <keyword arguments>)

Return a figure with an image `img` and displacement field `flow`.

# Keyword Arguments
- `img`: the image to be plotted.
- `flow::Flow`: the vector flow to be plotted.
- `disp_type::Symbol=:full` : display mode, either `:full` -> display all the vectors, or `:sparse` -> displays only the vectors above a threshold.
- `skip_count=nothing`: the number of vectors to skip between each displayed vector. By default set so that the output is ``20 × 20`` vectors.
- `fig=nothing`: add a figure to plot in. By defaults creates a blank new figure.
- `mag::Real=1`: magnify the plotted vectors.
- `key::Bool=true`: add key with maximum vector length.
- `figtitle::String="Iamge with Flow"`: add title to the figure.
- `ret::Symbol=:figure`: set return object, by default returns Figure, other options: :pyobject returns a PyObject. (Using figure makes Juno directly plot.)
- `origin_bot::Bool=true`: set the origin to the lower left corner, which is not typical for images (they will appear flipped).

See also: [`imgshow`](@ref), [`showflow`](@ref), [`warp_imgshowflow`](@ref)

"""
function imgshowflow(img,
                     flow::Flow;
                     disp_type::Symbol=:full,
                     skip_count=nothing, fig=nothing,
                     mag::Real=1,
                     key::Bool=true,
                     figtitle::String="Image with Flow",
                     ret::Symbol=:figure,
                     origin_bot::Bool=true)

    fig = imgshow(img, fig=fig, figtitle="", ret=:figure, origin_bot=origin_bot)
    showflow(flow, disp_type=disp_type, skip_count=skip_count, fig=fig, mag=mag, key=key, figtitle=figtitle, ret=ret);
end

"""
    warp_imgshowflow(img, flow; <keyword arguments>)

Return a figure with an image and a displacement field, where the image is warped by the displacement field.

# Keyword Arguments
- `img`: the image to be plotted.
- `flow::Flow`: the vector flow to be plotted.
- `disp_type::Symbol=:full` : display mode, either `:full` -> display all the vectors, or `:sparse` -> displays only the vectors above a threshold.
- `skip_count=nothing`: the number of vectors to skip between each displayed vector. By default set so that the output is ``20 × 20`` vectors.
- `fig=nothing`: add a figure to plot in. By defaults creates a blank new figure.
- `mag::Real=1`: magnify the plotted vectors.
- `key::Bool=true`: add key with maximum vector length.
- `figtitle::String="Iamge with Flow"`: add title to the figure.
- `ret::Symbol=:figure`: set return object, by default returns Figure, other options: :pyobject returns a PyObject. (Using figure makes Juno directly plot.)
- `origin_bot::Bool=true`: set the origin to the lower left corner, which is not typical for images (they will appear flipped).

See also: [`imgshow`](@ref), [`showflow`](@ref), [`imgshowflow`](@ref)
"""
function warp_imgshowflow(img,
                          flow::Flow;
                          disp_type::Symbol=:full,
                          skip_count=nothing,
                          fig=nothing, mag::Real=1,
                          key::Bool=true,
                          figtitle::String="Warped Image with Flow",
                          ret::Symbol=:figure,
                          origin_bot::Bool=true)

    imgw = warp_img(img, -real(flow), -imag(flow))
    imgshowflow(imgw, flow, disp_type=disp_type, skip_count=skip_count, fig=fig, mag=mag, key=key, figtitle=figtitle, ret=ret, origin_bot=origin_bot);
end


text:
"""
Function which calculates the covariance matrix of a cohort with respect to the relevant labels (for the first matrix calculation)
Arguments:
mycohort::Discohort -> contains all information of one cohort
wantedlabels::Array{ASCIIString,1} -> Array which contains the wantedlabels to give information about which covariances need to be calculated
"""
function calc_covarmat(mycohort::Discohort, wantedlabels::Array{ASCIIString,1})
	ncor = length(mycohort.y) - 1
	labelmap1 = labelmap(mycohort, wantedlabels)
	covmat = triu(ones(length(wantedlabels),length(wantedlabels)))
	@simd for i = 1 : length(labelmap1) - 1
		@simd for j = i + 1 : length(labelmap1)
			@inbounds covmat[i,j] = sum(mycohort.x[:, labelmap1[i]] .* mycohort.x[:, labelmap1[j]]) / ncor
		end
	end
	return Covarmat(covmat, wantedlabels)
end
"""
Function which calculates the expansions of the covariance matrix of a cohort with respect to the relevant labels
Arguments:
mycohort::Discohort -> contains all information of one cohort
wantedlabels::Array{ASCIIString,1} -> Array which contains the wantedlabels to give information about which covariances need to be calculated
usedlabels::Array{ASCIIString,1} -> Array which contains the labels which were already used to give information about which covariances need to be calculated
"""
function calc_covarmat(mycohort::Discohort, wantedlabels::Array{ASCIIString,1}, usedlabels::Array{ASCIIString,1})
	ncor = length(mycohort.y) - 1
	labelmap1 = labelmap(mycohort, wantedlabels)
	covmat = triu(ones(length(wantedlabels), length(wantedlabels)))
	covmat2 = zeros(length(usedlabels), length(wantedlabels))
	labelmap2 = labelmap(mycohort, usedlabels)
	@simd for i = 1 : max(length(labelmap1), length(labelmap2))
		@simd for j = 1 : max(length(labelmap1), length(labelmap2))
			if(i <= length(usedlabels) && j <= length(wantedlabels))
				@inbounds covmat2[i,j] = sum(mycohort.x[:,labelmap1[j]] .* mycohort.x[:,labelmap2[i]]) / ncor
			end
			if(j >= i+1 && i <= length(wantedlabels) && j <= length(wantedlabels))
				@inbounds covmat[i,j]  = sum(mycohort.x[:,labelmap1[i]] .* mycohort.x[:,labelmap1[j]]) / ncor
			end
		end
	end
	return Covarmat(vcat(covmat2, covmat), vcat(usedlabels, wantedlabels))
end
"""
Function for the calculation of the univariate estimators
Arguments:
cohort::Discohort -> contains all information about one cohort
"""
function calc_unibeta(cohort::Discohort)
	n = length(cohort.y)
	unibeta::Array{Float64} = zeros(size(cohort.x)[2])
	for i in 1:size(cohort.x)[2]
		@simd for j in 1:n
			@inbounds unibeta[i] += cohort.x[j,i]*cohort.y[j]
		end
		unibeta[i] /= n-1
	end
	return Unibeta(unibeta,ones(2500),cohort.labels)
end

max_stars_repo_name: danielazoeller/DistributedBoosting

max_issues_repo_name: danielazoeller/DistributedBoosting

max_forks_repo_name: danielazoeller/DistributedBoosting


text:
using Turing, Test, AbstractMCMC, DynamicPPL, Random, Turing.RandomMeasures, Libtask

@model infiniteGMM(x) = begin
    # Hyper-parameters, i.e. concentration parameter and parameters of H.
    α = 1.0
    μ0 = 0.0
    σ0 = 1.0

    # Define random measure, e.g. Dirichlet process.
    rpm = DirichletProcess(α)

    # Define the base distribution, i.e. expected value of the Dirichlet process.
    H = Normal(μ0, σ0)

    # Latent assignment.
    z = tzeros(Int, length(x))

    # Locations of the infinitely many clusters.
    μ = tzeros(Float64, 0)

    for i in 1:length(x)

        # Number of clusters.
        K = maximum(z)
        nk = Vector{Int}(map(k -> sum(z .== k), 1:K))

        # Draw the latent assignment.
        z[i] ~ ChineseRestaurantProcess(rpm, nk)

        # Create a new cluster?
        if z[i] > K
            push!(μ, 0.0)

            # Draw location of new cluster.
            μ[z[i]] ~ H
        end

        # Draw observation.
        x[i] ~ Normal(μ[z[i]], 1.0)
    end
end

# Generate some test data.
Random.seed!(1)

data = vcat(randn(10), randn(10) .- 5, randn(10) .+ 10)
data .-= mean(data)
data /= std(data)

# MCMC sampling
Random.seed!(2)
iterations = 500
model_fun = infiniteGMM(data)

m = Turing.Core.TracedModel(model_fun, Sampler(SMC(50)), VarInfo())
f = m.evaluator[1]
args = m.evaluator[2:end]

t = TapedTask(f, args...)

t.tf(args...)

@show Libtask.result(t.tf)

max_stars_repo_name: TuringLang/Libtask

max_issues_repo_name: TuringLang/Libtask

max_forks_repo_name: TuringLang/Libtask


text:
using LightGraphs: AbstractGraph, adjacency_matrix

struct Graph
    neighbors::Vector{Set{Int}}
    active::BitSet
    peo::Vector{Int}
end

# # SIMPLE LINEAR-TIME ALGORITHMS TO TEST CHORDALITY OF GRAPHS, TEST ACYCLICITY OF HYPERGRAPHS, AND SELECTIVELY REDUCE ACYCLIC HYPERGRAPHS
# # ROBERT E. TARJAN AND MIHALIS YANNAKAKIS
# Algorithm 4.2 in VA
"""
    maximum_cardinality_search(A)

Compute the perfect elimination ordering for a chordal graph represented by
a symmetric sparse matrix `A`. This function can also be used with AbstractGraph
objects from `LightGraphs.jl`.
"""
function maximum_cardinality_search(A::SparseMatrixCSC)
    !issymmetric(A) && error(ArgumentError("Matrix must be symmetric"))
    n = size(A, 1)
    peo = zeros(Int, n)
    cache = zeros(Int, n)

    _maximum_cardinality_search!(peo, A; n = n, cache=cache)
    return peo
end
maximum_cardinality_search(G::AbstractGraph) = maximum_cardinality_search(adjacency_matrix(G))

function _maximum_cardinality_search!(peo, A; n=n, cache=cache)
    v = 1
    for i in n:-1:1
        max_val = 0
        for j in 1:n
            peo[j] != 0 && continue
            if cache[j] > max_val
                max_val = cache[j]
                v = j
            end
        end
        peo[v] = i
        @views cache[(rowvals(A)[nzrange(A, v)])] .+= 1
    end

    return peo
end

"""
    is_chordal(A)

Tests if the graph represented by symmetric sparse matrix `A` is chordal. This
function can also be used with AbstractGraph objects from `LightGraphs.jl`.

# References
* [Simple linear-time algorithms to test chordality of graphs, test acyclicity of hypergraphs, and selectively reduce acyclic hypergraphs](https://epubs.siam.org/doi/pdf/10.1137/0213035?casa_token=A22jkwrsrL0AAAAA:rx-G6F21ubTkMiJmRTH3IKqxmFTo_IVWDDZfJig5lsZxnQtNH2vUKWfZ3eZJKUv9CiKbIPt1VQs) by Robert Tarjan and Mihalis Yannakakis
"""
function is_chordal(A::SparseMatrixCSC; peo=nothing)
    if isnothing(peo)
        peo = maximum_cardinality_search(A)
    end
    i_peo = invperm(peo)

    n = length(peo)
    f = zeros(Int, n)
    ind = zeros(Int, n)

    for i in 1:n
        w = i_peo[i]
        f[w] = w
        ind[w] = i
        for v in rowvals(A)[nzrange(A, w)]
            peo[v] >= i && continue
            ind[v] = i
            if f[v] == v
                f[v] = w
            end
        end

        for v in rowvals(A)[nzrange(A, w)]
            peo[v] >= i && continue
            ind[f[v]] < i && return false
        end
    end
    return true
end

is_chordal(G::AbstractGraph; peo=nothing) = is_chordal(adjacency_matrix(G); peo=peo)


text:
struct PriorityQueue{K, V}
    values::Array{Pair{K,V}}
    PriorityQueue{K,V}() where {K,V} = new(Array{Pair{K,V}}(undef, 0))
end

function push!(pq::PriorityQueue{K,V}, pair::Pair{K,V}) where {K,V}
    Base.push!(pq.values, pair)
    index = length(pq.values)
    while index > 1
        parent_index = div(index - 1, 2) + 1
        (pq.values[parent_index].first >= pair.first) && break
        pq.values[index] = pq.values[parent_index]
        index = parent_index
    end
    pq.values[index] = pair
    return nothing
end

function pop_highest!(pq::PriorityQueue)
    isempty(pq.values) && return nothing
    max_value = highest(pq)
    value = pop!(pq.values)
    index = 1
    while (index * 2) <= length(pq.values)
        cind1, cind2 = index * 2, index * 2 + 1
        flag = (cind2 >= length(pq.values) || pq.values[cind1] > pq.values[cind2])
        bigger_index = flag ? cind1 : cind2
        pq.values[bigger_index] <= value && break
        pq.values[index] = pq.values[bigger_index]
        index = bigger_index
    end
    pq.values[index] = value
    return max_value
end

pop_lowest!(pq::PriorityQueue) = !isempty(pq.values) ? pop!(pq.values) : nothing
highest(pq::PriorityQueue) = !isempty(pq.values) ? pq.values[1] : nothing
lowest(pq::PriorityQueue) = !isempty(pq.values) ? pq.values[end] : nothing


text:
abstract type DataSet end

copy(ds::DS) where {DS<:DataSet} = DS(fields(ds)...)

hash(ds::DataSet, h::UInt64) = foldr(hash, (typeof(ds), fieldvalues(ds)...), init=h)

function show(io::IO, ds::DataSet)
    println(io, typeof(ds), ": ")
    ds_dict = OrderedDict(k => getproperty(ds,k) for k in propertynames(ds) if k!=Symbol("_super"))
    for line in split(sprint(show, MIME"text/plain"(), ds_dict, context = (:limit => true)), "\n")[2:end]
        println(io, line)
    end
end

# util for distributing a singleton global dataset to workers
"""
    set_distributed_dataset(ds, [storage])
    get_distributed_dataset()

Sometimes it's more performant to distribute a DataSet object to
parallel workers just once, and have them refer to it from the global
state, rather than having it get automatically but repeatedly sent as
part of closures. This provides that functionality. Use
`set_distributed_dataset(ds)` from the master process to set the
global DataSet and `get_distributed_dataset()` from any process to
retrieve it. Repeated calls will not resend `ds` if it hasn't changed
(based on `hash(ds)`) and if no new workers have been added since the
last send. The optional argument `storage` will also adapt the dataset
to a particular storage on the workers, and can be a symbol, e.g.
`:CuArray`, in the case that CUDA is not loaded on the master process.
"""
function set_distributed_dataset(ds, storage=nothing)
    h = hash((procs(), ds, storage))
    if h != _distributed_dataset_hash
        @everywhere @eval CMBLensing _distributed_dataset = adapt(eval($storage), $ds)
        global _distributed_dataset_hash = h
    end
    nothing
end
get_distributed_dataset() = _distributed_dataset
_distributed_dataset = nothing
_distributed_dataset_hash = nothing

# Stores variables needed to construct the posterior
@kwdef mutable struct BaseDataSet <: DataSet
    d                # data
    Cϕ               # ϕ covariance
    Cf               # unlensed field covariance
    Cf̃ = nothing     # lensed field covariance (not always needed)
    Cn               # noise covariance
    Cn̂ = Cn          # approximate noise covariance, diagonal in same basis as Cf
    M  = I           # user mask
    M̂  = M           # approximate user mask, diagonal in same basis as Cf
    B  = I           # beam and instrumental transfer functions
    B̂  = B           # approximate beam and instrumental transfer functions, diagonal in same basis as Cf
    D  = I           # mixing matrix for mixed parametrization
    G  = I           # reparametrization for ϕ
    L  = LenseFlow   # lensing operator, possibly cached for memory reuse
    Nϕ = nothing     # some estimate of the ϕ noise, used in several places for preconditioning
end

function subblock(ds::DS, block) where {DS<:DataSet}
    DS(map(collect(pairs(fields(ds)))) do (k,v)
        @match (k,v) begin
            ((:Cϕ || :G || :L), v)              => v
            (_, L::Union{Nothing,FuncOp,Real})  => L
            (_, L)                              => getindex(L,block)
        end
    end...)
end

function (ds::DataSet)(θ)
    DS = typeof(ds)
    DS(map(fieldvalues(ds)) do v
        (v isa Union{ParamDependentOp,DataSet}) ? v(θ) : v
    end...)
end
(ds::DataSet)(;θ...) = ds((;θ...))

adapt_structure(to, ds::DS) where {DS <: DataSet} = DS(adapt(to, fieldvalues(ds))...)

@doc doc"""
    resimulate(ds::DataSet; [f, ϕ, n, f̃, rng, seed])

Make a new DataSet with the data replaced by a simulation. Keyword
argument fields will be used instead of new simulations, if they are
provided.

Returns a named tuple of `(;ds, f, ϕ, n, f̃)`.
"""
resimulate(ds::DataSet; kwargs...) = resimulate!(copy(ds); kwargs...)

@doc doc"""
    resimulate!(ds::DataSet; [f, ϕ, n, f̃, rng, seed])

Replace the data in this DataSet in-place with a simulation. Keyword
argument fields will be used instead of new simulations, if they are
provided.

Returns a named tuple of `(;ds, f, ϕ, n, f̃)`.
"""
function resimulate!(
    ds::DataSet;
    f=nothing, ϕ=nothing, n=nothing, f̃=nothing,
    Nbatch=(isnothing(ds.d) ? nothing : ds.d.Nbatch),
    rng=global_rng_for(ds.d), seed=nothing
)

    @unpack M,B,L,Cϕ,Cf,Cn,d = ds()

    if isnothing(f̃)
        if isnothing(ϕ)
            ϕ = simulate(Cϕ; Nbatch, rng, seed)
        end
        if isnothing(f)
            f = simulate(Cf; Nbatch, rng, seed = (isnothing(seed) ? nothing : seed+1))
        end
        f̃ = L(ϕ)*f
    else
        f = ϕ = nothing
    end
    if isnothing(n)
        n = simulate(Cn; Nbatch, rng, seed = (isnothing(seed) ? nothing : seed+2))
    end

    ds.d = d = M*B*f̃ + n

    (;ds,f,ϕ,n,f̃,d)

end

@doc doc"""

    load_sim(;kwargs...)

The starting point for many typical sessions. Creates a `BaseDataSet`
object with some simulated data, returing the DataSet and simulated
truths, which can then be passed to other maximization / sampling
functions. E.g.:

```julia
@unpack f,ϕ,ds = load_sim(;
    θpix  = 2,
    Nside = 128,
    pol   = :P,
    T     = Float32
)
```

Keyword arguments:

* `θpix` — Angular resolution, in arcmin.
* `Nside` — Number of pixels in the map as an `(Ny,Nx)` tuple, or a
  single number for square maps.
* `pol` — One of `:I`, `:P`, or `:IP` to select intensity,
  polarization, or both.
* `T = Float32` — Precision, either `Float32` or `Float64`.
* `storage = Array` — Set to `CuArray` to use GPU.
* `Nbatch = 1` — Number of batches of data in this dataset.
* `μKarcminT = 3` — Noise level in temperature in μK-arcmin.
* `ℓknee = 100` — 1/f noise knee.
* `αknee = 3` — 1/f noise slope.
* `beamFWHM = 0` — Beam full-width-half-max in arcmin.
* `pixel_mask_kwargs = (;)` — NamedTuple of keyword arguments to
  pass to `make_mask` to create the pixel mask.
* `bandpass_mask = LowPass(3000)` — Operator which performs
  Fourier-space masking.
* `fiducial_θ = (;)` — NamedTuple of keyword arguments passed to
  `camb()` for the fiducial model.
* `seed = nothing` — Specific seed for the simulation.
* `L = LenseFlow` — Lensing operator.

Returns a named tuple of `(;f, f̃, ϕ, n, ds, Cℓ, proj)`.

"""
function load_sim(;

    # basic configuration
    θpix,
    Nside,
    pol,
    T = Float32,
    storage = Array,
    Nbatch = 1,

    # noise parameters, or set Cℓn or even Cn directly
    μKarcminT = 3,
    ℓknee = 100,
    αknee = 3,
    Cℓn = nothing,
    Cn = nothing,

    # beam parameters, or set B directly
    beamFWHM = 0,
    B = nothing, B̂ = nothing,

    # mask parameters, or set M directly
    pixel_mask_kwargs = nothing,
    bandpass_mask = LowPass(3000),
    M = nothing, M̂ = nothing,

    # theory
    Cℓ = nothing,
    fiducial_θ = (;),
    rfid = nothing,

    rng = global_rng_for(storage),
    seed = nothing,
    D = nothing,
    G = nothing,
    Nϕ_fac = 2,
    L = LenseFlow,

)

    # projection
    Ny, Nx = Nside .* (1,1)
    proj = ProjLambert(;Ny, Nx, θpix, T, storage)

    # the biggest ℓ on the 2D fourier grid
    ℓmax = round(Int,ceil(√2*proj.nyquist)+1)

    # CMB Cℓs
    if (rfid != nothing)
        @warn "`rfid` will be removed in a future version. Use `fiducial_θ=(r=...,)` instead."
        fiducial_θ = merge(fiducial_θ,(r=rfid,))
    end
    Aϕ₀ = T(get(fiducial_θ, :Aϕ, 1))
    fiducial_θ = Base.structdiff(fiducial_θ, NamedTuple{(:Aϕ,)}) # remove Aϕ key if present
    if (Cℓ == nothing)
        Cℓ = camb(;fiducial_θ..., ℓmax=ℓmax)
    else
        if !isempty(fiducial_θ)
            error("Can't pass both `Cℓ` and `fiducial_θ` parameters which affect `Cℓ`, choose one or the other.")
        elseif maximum(Cℓ.total.TT.ℓ) < ℓmax
            error("ℓmax of `Cℓ` argument should be higher than $ℓmax for this configuration.")
        end
    end
    r₀ = T(Cℓ.params.r)

    # noise Cℓs (these are non-debeamed, hence beamFWHM=0 below; the beam comes in via the B operator)
    if (Cℓn == nothing)
        Cℓn = noiseCℓs(μKarcminT=μKarcminT, beamFWHM=0, ℓknee=ℓknee, αknee=αknee, ℓmax=ℓmax)
    end

    # some things which depend on whether we chose :I, :P, or :IP
    pol = Symbol(pol)
    ks,F,F̂,nF = @match pol begin
        :I  => ((:TT,),            FlatMap,    FlatFourier,    1)
        :P  => ((:EE,:BB),         FlatQUMap,  FlatEBFourier,  2)
        :IP => ((:TT,:EE,:BB,:TE), FlatIQUMap, FlatIEBFourier, 3)
        _   => throw(ArgumentError("`pol` should be one of :I, :P, or :IP"))
    end

    # covariances
    Cϕ₀ = Cℓ_to_Cov(:I,  proj, (Cℓ.total.ϕϕ))
    Cfs = Cℓ_to_Cov(pol, proj, (Cℓ.unlensed_scalar[k] for k in ks)...)
    Cft = Cℓ_to_Cov(pol, proj, (Cℓ.tensor[k]          for k in ks)...)
    Cf̃  = Cℓ_to_Cov(pol, proj, (Cℓ.total[k]           for k in ks)...)
    Cn̂  = Cℓ_to_Cov(pol, proj, (Cℓn[k]                for k in ks)...)
    if (Cn == nothing); Cn = Cn̂; end
    Cf = ParamDependentOp((;r=r₀,   _...)->(Cfs + (T(r)/r₀)*Cft))
    Cϕ = ParamDependentOp((;Aϕ=Aϕ₀, _...)->(T(Aϕ) * Cϕ₀))

    # data mask
    if (M == nothing)
        Mfourier = Cℓ_to_Cov(pol, proj, ((k==:TE ? 0 : 1) * bandpass_mask.diag.Wℓ for k in ks)...; units=1)
        if (pixel_mask_kwargs != nothing)
            Mpix = adapt(storage, Diagonal(F(repeated(T.(make_mask(Nside,θpix; pixel_mask_kwargs...).Ix),nF)..., proj)))
        else
            Mpix = I
        end
        M = Mfourier * Mpix
        if (M̂ == nothing)
            M̂ = Mfourier
        end
    else
        if (M̂ == nothing)
            M̂ = M
        end
    end
    if (M̂ isa DiagOp{<:BandPass})
        M̂ = Diagonal(M̂ * one(diag(Cf)))
    end

    # beam
    if (B == nothing)
        B = Cℓ_to_Cov(pol, proj, ((k==:TE ? 0 : 1) * sqrt(beamCℓs(beamFWHM=beamFWHM)) for k=ks)..., units=1)
    end
    if (B̂ == nothing)
        B̂ = B
    end

    # creating lensing operator cache
    Lϕ = alloc_cache(L(diag(Cϕ)),diag(Cf))

    # put everything in DataSet
    ds = BaseDataSet(;d=nothing, Cn, Cn̂, Cf, Cf̃, Cϕ, M, M̂, B, B̂, D, L=Lϕ)

    # simulate data
    @unpack ds,f,f̃,ϕ,n = resimulate(ds; rng, seed)

    # with the DataSet created, we now more conveniently create the mixing matrices D and G
    ds.Nϕ = Nϕ = quadratic_estimate(ds).Nϕ / Nϕ_fac
    if (G == nothing)
        G₀ = sqrt(I + Nϕ * pinv(Cϕ()))
        ds.G = ParamDependentOp((;Aϕ=Aϕ₀, _...)->(pinv(G₀) * sqrt(I + 2 * Nϕ * pinv(Cϕ(Aϕ=Aϕ)))))
    end
    if (D == nothing)
        σ²len = T(deg2rad(5/60)^2)
        ds.D = ParamDependentOp(
            function (;r=r₀, _...)
                Cfr = Cf(;r=r)
                sqrt((Cfr + I*σ²len + 2*Cn̂) * pinv(Cfr))
            end,
        )
    end

    if Nbatch > 1
        ds.d *= batch(ones(Int,Nbatch))
        ds.L = alloc_cache(L(ϕ*batch(ones(Int,Nbatch))), ds.d)
    end

    return (;f, f̃, ϕ, n, ds, ds₀=ds(), Cℓ, proj)

end

@deprecate load_sim_dataset(args...; kwargs...) load_sim(args...; kwargs...)