text:
export stationdensity_selection

"""
stationdensity_selection(InputDict::Dict)
return stationlist which is density-averaged over space.
"""
function stationdensity_selection(InputDict::Dict)

    starttime   = InputDict["starttime"]
    endtime     = InputDict["endtime"]
    kmlfile     = InputDict["kmlfile"]
    studyarea   = InputDict["studyarea"]
    includenet  = InputDict["includenet"]
    priority    = InputDict["priority"]
    TotalNumofStation = InputDict["TotalNumofStation"]
    nx          = InputDict["nx"]
    ny          = InputDict["ny"]
    IsPlotFigure     = InputDict["IsPlotFigure"]

    # parse kml station file
    println("---Parse KML file---")
    xdoc    = parse_file(kmlfile);
    xroot   = LightXML.root(xdoc);
    ces     = collect(child_elements(xroot));

    itmp = 1
    stid = 0
    while true
        itmp
        if !isempty(ces[itmp]["name"])
            stid = itmp
            break;
        end
        itmp += 1
        if itmp > 100
            error("Station network is more than 100. abort.")
        end
    end

    s1          = ces[stid:end]

    stationnum  = length(s1)
    stationlist = []
    # 1. read station info
    for i = 1:stationnum

        stationname = content(s1[i]["name"][1])
        stationloc = content(s1[i]["Point"][1]["coordinates"][1])
        coord_lon = parse.(Float64, split(stationloc, ","))[1]
        coord_lat = parse.(Float64, split(stationloc, ","))[2]

        # find operating range
        txt = string(s1[i]["Style"][1]["BalloonStyle"][1]["text"])
        tl  =split(string(txt), ";")
        operatingtime =tl[findfirst(occursin.("Operating Range", tl))]
        startoptime =  DateTime(split(operatingtime)[3][1:10])
        endoptime   =  DateTime(split(operatingtime)[5][1:10])

        datacentertxt =tl[findfirst(occursin.("Data Center", tl))]

        if (startoptime <= starttime && endtime <= endoptime
                && any(occursin.(stationname[1:2], includenet)))
            # take only stations which is operated over study time.
            st = Dict("stationname" => stationname,
                    "coord_lon" => coord_lon,
                    "coord_lat" => coord_lat)
            push!(stationlist, st)
        end
    end

    # 2. plot location
    if IsPlotFigure

        pyplot()

        fodir = "./fig"
        if !ispath(fodir) mkpath(fodir); end

        px = [stationlist[x]["coord_lon"] for x in 1:length(stationlist)]
        py = [stationlist[x]["coord_lat"] for x in 1:length(stationlist)]
        text = [stationlist[x]["stationname"] for x in 1:length(stationlist)]

        Plots.plot(bg=:white)
        p1 = Plots.scatter!(px, py,
                    markershape = :dtriangle,
                    markersize = 12,
                    label="",
                    size=(1200,1200),
                    xlabel = "lon",
                    xlims = (studyarea[4], studyarea[2]),
                    ylabel = "lat",
                    ylims = (studyarea[3], studyarea[1]),
                    aspect_ratio=:equal)

        Plots.savefig(p1, fodir*"/"*InputDict["figname"])

    end

    # 3. compute station density
    println("---Compute station density---")

    dx = (studyarea[2] - studyarea[4])/nx
    dy = (studyarea[1] - studyarea[3])/ny
    NG = nx * ny
    grid = []
    for i=1:nx
        for j=1:ny
            xl = (i-1) * dx + studyarea[4]
            xr = (i) * dx + studyarea[4]
            yb = (j-1) * dy + studyarea[3]
            yt = (j) * dy + studyarea[3]
            xm = (xl+xr)/2
            ym = (yb+yt)/2
            gridtmp = Dict("xm"=>xm, "ym"=>ym, "xl"=>xl, "xr"=>xr, "yb"=>yb, "yt"=>yt)
            push!(grid, gridtmp)
        end
    end

    stationidlist = []
    for gid = 1:NG
        xl = grid[gid]["xl"]
        xr = grid[gid]["xr"]
        yb = grid[gid]["yb"]
        yt = grid[gid]["yt"]
        numstationpergrid=0
        stationid=Int[]
        for i = 1:length(stationlist)
            x1 = stationlist[i]["coord_lon"]
            y1 = stationlist[i]["coord_lat"]

            if insquare(x1, y1, xl, xr, yb, yt, edgeinclude=true)
                numstationpergrid += 1
                push!(stationid, i)
            end
        end
        grid[gid]["numstationpergrid"] = numstationpergrid
        push!(stationidlist, stationid)
    end

    # 4. reduce the station density
    #   - sort from higher num of station
    #   - remove one of pairs with closest distance and priority
    println("---Averaging station density---")

    r_grid= deepcopy(grid)
    r_stationlist = deepcopy(stationlist)
    r_stationidlist = deepcopy(stationidlist)

    removestations!(TotalNumofStation, r_grid, r_stationlist, r_stationidlist, includenet, priority)

    r_stationlist = [stationlist[x] for x in collect(Iterators.flatten(r_stationidlist))]

    if IsPlotFigure

        px = [r_stationlist[x]["coord_lon"] for x in 1:length(r_stationlist)]
        py = [r_stationlist[x]["coord_lat"] for x in 1:length(r_stationlist)]
        text = [r_stationlist[x]["stationname"] for x in 1:length(r_stationlist)]

        Plots.plot(bg=:white)
        p2 = Plots.scatter!(px, py,
                    markershape = :dtriangle,
                    markersize = 12,
                    label="",
                    size=(1200,1200),
                    xlabel = "lon",
                    xlims = (studyarea[4], studyarea[2]),
                    ylabel = "lat",
                    ylims = (studyarea[3], studyarea[1]),
                    aspect_ratio=:equal)

        Plots.savefig(p2, fodir*"/"*"densityaveraged_"*InputDict["figname"])
    end

    #===
    save for SeisDownload input file
    ===#

    fofiledir = "./files"
    if !ispath(fofiledir) mkpath(fofiledir); end
    stationlistfoname = fofiledir*"/input_stationlist.jld2"
    jldopen(stationlistfoname, "w") do file
        file["stationlist"] = stationlist;
    end

    #===
    # save stationlist for gmt
    ===#

    # gather network name
    networks = Dict()
    println("---Station networks---")
    for i = 1:length(stationlist)
        networkname =  split(stationlist[i]["stationname"], ".")[1]
        if !haskey(networks, networkname)
            println(networkname)
            #add this network to dict
            networks["$(networkname)"] = Array[]
        end
    end

    for i = 1:length(stationlist)
        stationname = stationlist[i]["stationname"]
        coord_lon = stationlist[i]["coord_lon"]
        coord_lat = stationlist[i]["coord_lat"]

        # Here, elevation is not included because no elevetion in IRIS kml file.
        # However, you can find elevation when analysing dv/v as CorrData has the elevation of stations.

        for key = keys(networks)
            if occursin(stationname[1:3], "$(key).", )
                push!(networks["$(key)"], [coord_lon, coord_lat])
            end
        end
    end

    #save as jld2
    jldopen(fofiledir*"/gmt_stationcoords.jld2", "w") do file
        for key = keys(networks)
            file["$(key)"] = networks["$(key)"]
        end
    end

    return nothing

end

"""
insquare(x1, y1, xl, xr, yb, yt; edgeinclude::Bool=true)
return if x1, y1 is in square
"""
function insquare(x1, y1, xl, xr, yb, yt; edgeinclude::Bool=true)
    if edgeinclude
        if x1 >= xl && x1 <= xr && y1 >= yb && y1 <= yt
            return true
        else
            return false
        end
    else
        if x1 > xl && x1 < xr && y1 > yb && y1 < yt
            return true
        else
            return false
        end
    end
end

"""
removestations!(TotalNumofStation::Int, r_grid::Array{Any,1}, r_stationlist::Array{Any,1},
     r_stationidlist::Array{Any,1}, includenet::Array{String,1}, priority::Array{String,1})
remove stations with averaging density in grid
"""
function removestations!(TotalNumofStation::Int, r_grid::Array{Any,1}, r_stationlist::Array{Any,1},
     r_stationidlist::Array{Any,1}, includenet::Array{String,1}, priority::Array{String,1})

     NG = length(r_grid)

     r_numstlistsort = [r_grid[x]["numstationpergrid"] for x in 1:NG]
     r_sortid = sortperm(r_numstlistsort, rev=true)

     icount = 1
     while true

         if sum(r_numstlistsort) <= TotalNumofStation
             break
         elseif sum(r_numstlistsort) == 0
             error("All station is removed from list. abort.")
         elseif icount > 10000
             error("iteration number exceeds 10000. too many station removal. abort.")
         end

         r_gid = first(r_sortid) # remove from this grid
         st = r_stationidlist[r_gid]
         # compute distance with eachother
         rdist = 1e12
         r_ii = -1
         r_jj = -1
         r_netii = ""
         r_netjj = ""

         for i = 1:length(st)-1
             for j = i+1:length(st)

                 ii = st[i]
                 jj = st[j]

                 id1 = r_stationlist[ii]["stationname"]
                 id2 = r_stationlist[jj]["stationname"]
                 net1 =  string(split(id1, ".")[1])
                 net2 =  string(split(id2, ".")[1])
                 x1 = r_stationlist[ii]["coord_lon"]
                 y1 = r_stationlist[ii]["coord_lat"]
                 x2 = r_stationlist[jj]["coord_lon"]
                 y2 = r_stationlist[jj]["coord_lat"]
                 trial_dist = norm([(x2-x1), (y2-y1)])

                 # decide if we remove this
                 if  @isdefined includenet
                     if  any(occursin.(includenet, net1)) || any(occursin.(includenet, net2))
                         continue;
                     end
                 end

                 if trial_dist < rdist
                     #update potential removal station
                     r_ii = ii
                     r_jj = jj
                     r_netii = net1
                     r_netjj = net2
                     rdist = trial_dist
                 end
             end
         end

         if r_ii == -1 || r_jj == -1
             error("no station is chosen. abort.")
         end

         #evaluate removing either r_ii or r_jj
         if  !@isdefined priority
             priority = []
         end

         ri1 = findfirst(x -> x==r_netii, priority)
         ri2 = findfirst(x -> x==r_netjj, priority)

         if isnothing(ri1) || isnothing(ri2)
             removedid = r_ii #no priority
         elseif ri1 >= ri2
             removedid = r_ii
         else
             removedid = r_jj
         end

         # pop out removedid from list
         r_grid[r_gid]["numstationpergrid"] = r_numstlistsort[r_gid] - 1
         filter!(e -> e != removedid, r_stationidlist[r_gid])

         r_numstlistsort = [r_grid[x]["numstationpergrid"] for x in 1:NG]
         r_sortid = sortperm(r_numstlistsort, rev=true)

         icount += 1
     end
end


text:
# ---
# title: 996. Number of Squareful Arrays
# id: problem996
# author: Indigo
# date: 2021-05-31
# difficulty: Hard
# categories: Math, Backtracking, Graph
# link: <https://leetcode.com/problems/number-of-squareful-arrays/description/>
# hidden: true
# ---
#
# Given an array `A` of non-negative integers, the array is _squareful_ if for
# every pair of adjacent elements, their sum is a perfect square.
#
# Return the number of permutations of A that are squareful.  Two permutations
# `A1` and `A2` differ if and only if there is some index `i` such that `A1[i]
# != A2[i]`.
#
#
#
# **Example 1:**
#
#
#
#     Input: [1,17,8]
#     Output: 2
#     Explanation:
#     [1,8,17] and [17,8,1] are the valid permutations.
#
#
# **Example 2:**
#
#
#
#     Input: [2,2,2]
#     Output: 1
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 12`
#   2. `0 <= A[i] <= 1e9`
#
#
## @lc code=start
using LeetCode

function can_adj(u, v)
    sq = isqrt(u + v)
    return sq * sq == u + v
end

function num_squareful_perms(nums::Vector{Int})
    ctr = counter(nums)
    adj_set = Set{Tuple{Int, Int}}()
    for i1 in 1:length(nums), i2 in i1:length(nums)
        n1, n2 = nums[i1], nums[i2]
        if can_adj(n1, n2)
            push!(adj_set, (n1, n2))
            push!(adj_set, (n2, n1))
        end
    end
    function gen_perms(ctr, lst = nothing)
        sum(ctr) == 0 && return 1
        res = 0
        for k in keys(ctr)
            (ctr[k] == 0 || !isnothing(lst) && (lst, k) ∉ adj_set) && continue
            ctr[k] -= 1
            res += gen_perms(ctr, k)
            ctr[k] += 1
        end
        res
    end
    gen_perms(ctr)
end
## @lc code=end

ext: jl


text:
################
###Photoreceptor
################
abstract AbstractPhotoReceptor

numtimesteps(r::AbstractPhotoReceptor) = r.numsteps
numvilli(r::AbstractPhotoReceptor) = r.numvilli
numphotons(r::AbstractPhotoReceptor) = r.numphotons

photoreceptor(s::Symbol,photons::Vector,nvilli::Integer,args...) = _photoreceptor(Val{s},photons,nvilli,args...)

###############################################
###Photons stored in a matrix (dense or sparse)
###############################################

type PhotoReceptorMatrix{T<:Integer,A<:AbstractMatrix} <: AbstractPhotoReceptor
    microvilli::A
    numsteps::Int
    numvilli::Int
    numphotons::Int
    PhotoReceptorMatrix(mv::AbstractMatrix{T},ns::Int,nv::Int,np::Int) = new(mv,ns,nv,np)
end

### Dense Matrix factory functions

function _photoreceptor{T<:Integer}(::Type{Val{:full}},photons::Vector,nvilli::Integer,::Type{T})
    nsteps = length(photons)
    nphotons = sum(photons)
    m = zeros(T,nsteps,nvilli)
    for i = 1:nsteps
        r = rand(1:nvilli,photons[i])
        for j in r
            @inbounds m[i,j] += one(T)
        end
    end
    PhotoReceptorMatrix{T,typeof(m)}(m,nsteps,nvilli,nphotons)
end

### Sparse Matrix factory functions (row version)

@inline function _photontimeindex{I<:Integer}(photons::Vector,::Type{I})
    e = cumsum(photons)
    photontime = Vector{I}(e[end])
    for i=one(I):I(length(photons))
        @simd for j=e[i]-photons[i]+1:e[i]
            @inbounds photontime[j] = i
        end
    end
    photontime
end

@inline _randvilliindex{I<:Integer}(nvilli::Int,nphotons::Int,::Type{I}) = rand(one(I):I(nvilli),nphotons)

function _photoreceptor{T<:Integer,I<:Integer}(::Type{Val{:row}},photons::Vector,nvilli::Integer,::Type{T},::Type{I})
    pti = _photontimeindex(photons,I)
    nphotons = length(pti)
    nsteps = length(photons)
    m = sparse(_randvilliindex(nvilli,nphotons,I),pti,ones(T,nphotons),nvilli,nsteps)
    PhotoReceptorMatrix{T,typeof(m)}(m,nsteps,nvilli,nphotons)
end

@inline microvilliwithphotons{I<:Integer,M<:DenseMatrix}(r::PhotoReceptorMatrix{I,M},timepoint::Integer) = find(r.microvilli[timepoint,:])
@inline microvilliwithphotons{I<:Integer,S<:SparseMatrixCSC}(r::PhotoReceptorMatrix{I,S},timepoint::Integer) = view(rowvals(r.microvilli),nzrange(r.microvilli,timepoint))

################################
###Photons stored in cell matrix
################################

type PhotoReceptorCell{T<:Integer} <: AbstractPhotoReceptor
    microvilli::Vector{Vector{T}}
    numsteps::Int
    numvilli::Int
    numphotons::Int
    PhotoReceptorCell(mv::Vector{Vector{T}},ns::Int,nv::Int,np::Int) = new(mv,ns,nv,np)
end

function _photoreceptor{T<:Integer}(::Type{Val{:cell}},photons::Vector,nvilli::Integer,::Type{T})
    PhotoReceptorCell{T}([sort!(rand!(Vector{T}(p),one(T):T(nvilli))) for p in photons],length(photons),nvilli,sum(photons))
end

@inline microvilliwithphotons(r::PhotoReceptorCell,timepoint::Integer) = r.microvilli[timepoint]

function show(io::IO,r::AbstractPhotoReceptor)
    println(io,"PhotoReceptor with $(r.numvilli) microvilli and $(r.numsteps) timesteps")
    println(io,"  Additional fields: :microvilli")
    nothing
end


text:
export bin, bin2

"""
    bin(arr, binning)

Bins a `arr` by the factors given in `binning`.

For efficient 2 binning, see [`bin2`](@ref bin2).

# Examples
## Simple Usage
```jldoctest
julia> bin([1 2 3 4; 5 6 7 8], (2, 2))
1×2 Matrix{Int64}:
 14  22

julia> bin(ones((4, 4)), (2, 4))
2×1 Matrix{Float64}:
 8.0
 8.0

julia> bin([1,2,3,4,5], (2,))
2-element Vector{Int64}:
 3
 7
```

## You can leave out trailing 1s
```jldoctest
julia> bin(ones((2,2,2)), (2,))
1×2×2 Array{Float64, 3}:
[:, :, 1] =
 2.0  2.0

[:, :, 2] =
 2.0  2.0

julia> bin(ones((2,2,2)), (2,)) == bin(ones((2,2,2)), (2,1,1))
true
```
"""
function bin(arr::AbstractArray{T, N}, binning::NTuple{M, Int}) where {T, N, M}
    # an tuple indicating the binning in each dimnension
    # fill with 1s of binning only specifies a few dimensions
    binning = (binning..., ntuple(x -> 1, Val(N - M))...)

    # output array
    # largest size which is possible according to binning
    out = similar(arr, size(arr) .÷ binning)

    out_indices = CartesianIndices(out)
    z = zero(T)
    neighbours = CartesianIndices(binning)
    # iterate over the indices of the out put array
    #@inbounds and @simd make it faster
    @inbounds @simd for inds in out_indices
        # now we start to sum within binning
        s = z
        # sum all neighbours
        base_inds = binning .* Tuple(inds) .+ 1
        for neighbour in neighbours
            pos = base_inds .- Tuple(neighbour)
            s += arr[pos...]
        end
        # assign the sum
        out[inds] = s
    end
    return out
end

# was slightly slower with function barrier
# Function barrier.
# """
# function _bin_loop!(out::AbstractArray{T, N},
#                     arr::AbstractArray{T, N},
#                     binning::NTuple{N, Int}) where {T, N}
## end

"""
    bin2(arr)

Essentially calls `bin(arr, (2,2))`.
However, works only for `Union{AbstractArray{T, 2}, AbstractArray{T, 3}} where T`
since we use specialized methods of Tullio for that.

# Examples
```jldoctest
julia> x = [1 2 3 4; 5 6 7 8]
2×4 Matrix{Int64}:
 1  2  3  4
 5  6  7  8

julia> bin2(x)
1×2 Matrix{Int64}:
 14  22

julia> bin2(x) ≈ bin(x, (2,2))
true
```
"""
function bin2(arr::AbstractArray{T, 2}) where T
    @tullio out[i, j] := arr[2i, 2j] +  arr[2i - 1, 2j] + arr[2i - 1, 2j - 1] + arr[2i, 2j - 1]
end

function bin2(arr::AbstractArray{T, 3}) where T
    @tullio out[i, j, k] := arr[2i, 2j, k] +  arr[2i - 1, 2j, k] + arr[2i - 1, 2j - 1, k] + arr[2i, 2j - 1, k]
end

max_stars_repo_name: JuliaMicroscopy/MicroscopyTools

max_issues_repo_name: JuliaMicroscopy/MicroscopyTools


text:
struct EquispacedValueGrid{BC<:AbstractBC, T} <: AbstractGrid{BC, T}
  N::Int
  L::Float64
  values::Vector{T}
  bc::BC
end
function EquispacedValueGrid(N::Int, L::Real, ::Type{BC}=PeriodicGridBC,
   ::Type{T}=Float64) where {BC<:AbstractBC, T<:Number}
  return EquispacedValueGrid{BC,T}(N, L, zeros(T, N), BC(0.0, L))
end
Base.length(g::EquispacedValueGrid) = g.N
Base.ndims(::Type{EquispacedValueGrid{BC, T}}) where {BC<:AbstractBC, T} = 1

(g::EquispacedValueGrid)(x::AbstractFloat) = basis(g, cell(x, g))(x) * g.values[cell(x, g)]

cellwidth(g::EquispacedValueGrid) = g.L / g.N
cell(x, g::EquispacedValueGrid) = Int(floor(x / cellwidth(g))) + 1
cellcentre(g::EquispacedValueGrid, i::Integer) = (i - 0.5) * cellwidth(g)
cellcentres(g::EquispacedValueGrid) = (cellcentre(g, i) for i ∈ 1:g.N)
cellindices(ab, g::EquispacedValueGrid) = cellindices(ab..., g)
numberofunknowns(g::EquispacedValueGrid) = g.N
function cellindices(a, b, g::EquispacedValueGrid)
  return Int(floor(a / cellwidth(g))):Int(ceil(b / cellwidth(g)))
end
basis(g, i) = BasisFunction(TopHatShape(cellwidth(g)), cellcentre(g, i), g.values[i])
bases(g) = [basis(g, i) for i ∈ 1:g.N]

function Base.intersect(x::BasisFunction{S, T1}, g::EquispacedValueGrid{BC, T2}
    ) where {S<:AbstractShape, BC, T1, T2}
  #return enumerate(bases(gf)) # test against this
  indlo = Int(floor(lower(x) / cellwidth(g) - 0.5))
  indhi = Int(ceil(upper(x) / cellwidth(g) + 0.5))
  return ((mod1(i, g.N), basis(g, mod1(i, g.N))) for i ∈ indlo:indhi)
end

function Base.isapprox(a::T, b::T, atol=0, rtol=sqrt(eps())) where {T<:EquispacedValueGrid}
  a.N == b.N && return false
  a.L == b.L && return false
  return isapprox(a.values, b.values, atol=atol, rtol=rtol)
end

# AbstractArray interface
Base.size(g::EquispacedValueGrid) = (g.N,)
# EquispacedValueGrid
Base.getindex(g::EquispacedValueGrid, i) = g.values[i]
Base.setindex!(g::EquispacedValueGrid, v, i) = (g.values[i] = v)
Base.iterate(g::EquispacedValueGrid) = iterate(g.values)
Base.iterate(g::EquispacedValueGrid, state) = iterate(g.values, state)

zero!(g::EquispacedValueGrid) = (g.values .*= false)
#weight(g::EquispacedValueGrid, i) = g.values[i] * width(basis(g, i))

function deposit!(evg::EquispacedValueGrid{BC}, particle) where {BC<:AbstractBC}
  # loop over all items in evg that particle overlaps with
  bc = BC(0.0, evg.L)
  qw = charge(particle) * weight(particle)
  for (index, item) ∈ intersect(basis(particle), evg)
    # multiply by width because all basis functions are normalised
    amountincell = integral(item, basis(particle), bc) * width(item)
    evg[index] += amountincell * qw
  end
  return evg
end

function antideposit(evg::EquispacedValueGrid{BC}, particle) where {BC<:AbstractBC}
  # loop over all items in evg that particle overlaps with
  bc = BC(0.0, evg.L)
  amount = 0.0
  for (index, item) ∈ intersect(basis(particle), evg)
    # multiply by width because all basis functions are normalised
    amount += integral(item, basis(particle), bc) * evg.values[index]
  end
  return amount
end


text:
"""
    MutationOp
Perform mutation operations
"""
abstract type MutationOp end
abstract type MutationState <: MutationOp end

function reset_in!(s::MutationOp) end
function reset_out!(s::MutationOp) end
function reset!(s::MutationOp) end

"""
    nin(m::MutationState)
    nin(v::AbstractMutationVertex)

Return the size of the input data to the vertex

Computation may fail if input does not have size nin.
"""
function nin end

"""
    nout(m::MutationState)
    nout(v::AbstractMutationVertex)

Return the size of the output data from the vertex.
"""
function nout end

"""
    nin_org(m::MutationState)
    nin_org(v::AbstractMutationVertex)

Return the size of the input data to the vertex before any mutation was performed.

Computation may fail if input does not have size nin.
Computation may also fail if vertex has MutationState which does not track changes.
"""
function nin_org end

"""
    nout_org(m::MutationState)
    nout_org(v::AbstractMutationVertex)

Return the size of the output data from the vertex before any mutation was performed.
Computation may fail if vertex has MutationState which does not track changes.
"""
function nout_org end

"""
    Δnin(o::MutationOp, Δ)
    Δnin(v::AbstractMutationVertex, Δ)

Change input size by Δ. New size is nin + Δ

Might propagate to other vertices when performed on an AbstractMutationVertex
depending on vertex type
"""
function Δnin end

"""
    Δnout(o::MutationOp, Δ)
    Δnout(v::AbstractMutationVertex, Δ)

Change output size by Δ. New size is nout + Δ

Might propagate to other vertices when performed on an AbstractMutationVertex
depending on vertex type
"""
function Δnout end

"""
    in_inds(m::MutationState)

Return selected input indices.
"""
function in_inds end

"""
    out_inds(m::MutationState)

Return selected output indices.
"""
function out_inds end

"""
    NoOp
Does not perform any operations
"""
struct NoOp <: MutationOp end
clone(op::NoOp) = op
Δnin(s::NoOp, Δ...) = Δ
Δnout(s::NoOp, Δ) = Δ

"""
    IoSize
Size for input and output of a computation
"""
mutable struct IoSize{I, O} <: MutationState
    nin::I
    nout::O
end
IoSize(size::Integer) = IoSize(size, size)
IoSize(in::Integer, out::Integer) = IoSize([in], out)
clone(s::IoSize) = IoSize(copy(nin(s)), nout(s))

nin(s::IoSize) = s.nin
nout(s::IoSize) = s.nout

in_inds(s::IoSize) = [1:insize for insize in s.nin]
out_inds(s::IoSize) = 1:s.nout

Maybe{T} = Union{T, Missing}
Δnin(s::IoSize, Δ::Maybe{Integer}...) = s.nin .+= replace(collect(Δ), missing => 0)
Δnout(s::IoSize, Δ::Integer) = s.nout += Δ
function Δnin(s::IoSize, Δ::Maybe{AbstractArray{<:Integer,1}}...)
    for (i,Δi) in enumerate(Δ)
        if !ismissing(Δi)
            s.nin[i] = length(Δi)
        end
    end
end
Δnout(s::IoSize, Δ::AbstractArray{<:Integer,1}) = s.nout = length(Δ)

"""
    IoIndices
Indexes to retain for input and output of a computation

What those indices mean is up to the computation
"""
mutable struct IoIndices{I,O} <: MutationState
    in::I
    out::O
end
IoIndices(in::Integer, out::Integer) = IoIndices([in], out)
IoIndices(in::AbstractArray{<:Integer}, out::Integer) = IoIndices(collect.(range.(1, in, step=1)), collect(1:out))
IoIndices(s::MutationState) = IoIndices(nin(s), nout(s))
clone(s::IoIndices) = IoIndices(deepcopy(s.in), copy(s.out))

nin(s::IoIndices) = length.(s.in)
nout(s::IoIndices) = length(s.out)

in_inds(s::IoIndices) = s.in
out_inds(s::IoIndices) = s.out

Δnin(s::IoIndices, Δ::Maybe{AbstractArray{<:Integer,1}}...) = Δnin(s::IoIndices, map(i -> ismissing(Δ[i]) ? s.in[i] : Δ[i], eachindex(Δ))...)
Δnin(s::IoIndices, Δ::AbstractArray{<:Integer,1}...) = s.in = collect(deepcopy(Δ))
Δnout(s::IoIndices, Δ::AbstractArray{<:Integer,1}) = s.out = copy(Δ)
function reset_in!(s::IoIndices)
    foreach(s.in) do ini
        ini[1:end] = 1:length(ini)
    end
end
function reset_out!(s::IoIndices)
    s.out[1:end] = 1:length(s.out)
end

"""
    IoChange

Size for input and output of a computation stored as a change towards an original size.

Also maintains indices separately.

Useful for tracking change through several mutations before making a selection of which parameters to keep.
"""
mutable struct IoChange{S,I,DI,DO} <: MutationState
    size::S
    indices::I
    inΔ::DI
    outΔ::DO
end

IoChange(insize, outsize) =  IoChange(IoSize(insize, outsize), IoIndices(insize, outsize))
function IoChange(size::MutationState, inds::MutationState)
    nin(size) == nin(inds) || error("Insizes differ: $(nin(size)) vs $(nin(inds))")
    nout(size) == nout(inds) || error("Outsizes differ: $(nout(size)) vs $(nout(inds))")
    IoChange(size, inds, zeros(eltype(nin(size)), length(nin(size))), 0)
end
clone(s::IoChange) = IoChange(clone(s.size), clone(s.indices), copy(s.inΔ), s.outΔ)

nin(s::IoChange) = nin_org(s) .+ s.inΔ
nout(s::IoChange) = nout(s.size) + s.outΔ

in_inds(s::IoChange) = trunc_or_pad.(in_inds(s.indices), nin(s))
out_inds(s::IoChange) = trunc_or_pad(out_inds(s.indices), nout(s))

nin_org(s::IoChange) = nin(s.size)
nout_org(s::IoChange) = nout(s.size)

function trunc_or_pad(vec, size)
    res = -ones(eltype(vec), size)
    inds = 1:min(size, length(vec))
    res[inds] = vec[inds]
    return res
end

Δnin(s::IoChange, Δ::Maybe{Integer}...) = s.inΔ .+= replace(collect(Δ), missing => 0)
Δnout(s::IoChange, Δ::Integer) = s.outΔ += Δ
function Δnin(s::IoChange, Δ::Maybe{AbstractArray{<:Integer,1}}...)
    Δnin(s.indices, Δ...)

    for (i, Δi) in enumerate(Δ)
        if !ismissing(Δi)
            s.inΔ[i] = length(Δi) - nin(s.size)[i]
        end
    end
end

function Δnout(s::IoChange, Δ::AbstractArray{<:Integer,1})
    Δnout(s.indices, Δ)
    s.outΔ = nout(s.indices) - nout(s.size)
end

function reset_in!(s::IoChange)
    Δnin(s.size, s.inΔ...)
    s.inΔ[1:end] .= 0
    Δnin(s.indices, collect.(StepRange.(1,1, nin(s)))...)
    reset_in!(s.size)
    reset_in!(s.indices)
end

function reset_out!(s::IoChange)
    Δnout(s.size, s.outΔ)
    s.outΔ = 0
    Δnout(s.indices, collect(1:nout(s)))
    reset_out!(s.size)
    reset_out!(s.indices)
end


text:
using NLSolvers, Test

solve(NLEProblem(x->x^2), 0.001, Secant(;factor=0.99, shift=1e-1))
solve(NLEProblem(x->sin(x)), 0.001, Secant(;factor=0.99, shift=1e-1))
solve(NLEProblem(x->cos(x)), 0.001, Secant(;factor=0.99, shift=1e-1))

# @test_throws ArgumentError solve(NLEProblem(x->x^2), [0.001,3], Secant(;factor=0.99, shift=1e-8))


text:
using JuMP
using GLPK
using Graphs
using GraphPlot
using GraphIO

function maximum_independent_set(graph; show_model = false)
    model = Model(GLPK.Optimizer)
    @variable(model, x[1:Graphs.nv(graph)], Bin)
    @objective(model, Max, sum(x))
    for e in Graphs.edges(graph)
        # To access the first and second vertec from each edge in
        # from Graphs.edges(graph) you must call the src and dst
        # functions, respectively.
        @constraint(model, x[Graphs.src(e)] + x[Graphs.dst(e)] <= 1)
    end
    JuMP.optimize!(model)
    if show_model
        println(model)
    end
    return [xi for xi in x if JuMP.value(xi) == 1.0]
end

# Read in the Petersen graph edge list text file
graph = Graphs.loadgraph("PetersenGraph.txt", "graph_key", GraphIO.EdgeListFormat())
# Convert the directed graph to a simple graph
graph = Graphs.SimpleGraph(graph)

# Solve the maximimum independent set problem on the Petersen graph
max_independent_set = maximum_independent_set(graph, show_model = true)

println(max_independent_set)


text:
function create_partition(lengths::NTuple{N,Int}) where N
    num = 0
    partition = UnitRange{Int}[]
    len = [0; cumsum(collect(lengths))]
    partition = ntuple(i->(1:lengths[i]) .+ len[i], length(lengths))
    return partition
end

function create_partition(lengths::NTuple{N,Int}, names::NTuple{N,Symbol}) where {N,T}
    partition = Tuple(create_partition(lengths))
    named_part = NamedTuple{names}(partition)
    return named_part
    # return named_part::NamedTuple{M,NTuple{N,UnitRange{Int}}}
end

function create_partition2(len1::NTuple{N1,Int},len2::NTuple{N2,Int}) where {N1,N2}
    n1 = length(len1)
    n2 = length(len2)
    num1 = [0; cumsum(collect(len1))]
    num2 = [0; cumsum(collect(len2))]

    inds = [(i,j) for i = 1:n1 for j=1:n2]
    inds = Tuple(1:n1*n2)
    function get_inds(x)
        i = cart1(x,n1,n2)
        j = cart2(x,n1,n2)
        (1:len1[i]) .+ num1[i], (1:len2[j]) .+ num2[j]
    end
    partition = map(get_inds, ntuple(i->i,n1*n2))
    return partition
end
create_partition2(lengths::NTuple{N,Int}) where N = create_partition2(lengths,lengths)

function cart1(i,n,m)
    ((i-1) ÷ m)+1
end

function cart2(i,n,m)
    v = i % m
    v == 0 ? m : v
end

function create_partition2(len1::NTuple{N1,Int},len2::NTuple{N2,Int},
        names::Val{NAMES}) where {N1,N2,NAMES}
    n1 = length(len1)
    n2 = length(len2)
    partition = create_partition2(len1,len2)
    named_part = create_nt(names,partition)
end

function create_nt(::Val{names},part) where {names}
    NamedTuple{names}(part)
end

function combine_names(names1,names2)
    n1 = length(names1)
    n2 = length(names2)
    function get_inds(x)
        i = cart1(x,n1,n2)
        j = cart2(x,n1,n2)
        Symbol(string(names1[i])*string(names2[j]))
    end
    partition = map(get_inds, ntuple(i->i,n1*n2))
end


text:
#############################################################################
#Copyright (c) 2016-2019 Charles Le Losq
#
#The MIT License (MIT)
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the #Software without restriction, including without limitation the rights to use, copy, #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, #and to permit persons to whom the Software is furnished to do so, subject to the #following conditions:
#
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#
#############################################################################

"""
	tlcorrection(data::Array{Float64},temp::Float64,wave::Float64;correction="long",normalisation="area",density=2210.0)

Inputs
------
	data: Array{Float64}
		input spectrum with x and y in first and second columns respectively
	temp: Float64
		temperature in °C
	wave: Float64
		wavenumber at which the spectrum was acquirred in nm

Options
-------
	correction: String
		equation used for the correction. Choose between "long", "galeener", or "hehlen". Default = "long".
	normalisation: String
		indicate if you want to normalise your signal or not. Choose between "intensity", "area", or "no". Default = "area".
	density: Float64
		density of the studied material in kg m-3, to be used with the "hehlen" equation. Default = 2210.0 (density of silica).

Outputs
-------

	x: Array{Float64}
		x values
	long: Array{Float64}
		corrected y values
	eselong: Array{Float64}
		errors calculated as sqrt(y) on raw data and propagated after the correction.

Notes
-------

This correction uses the formula reported in Galeener and Sen (1978), Mysen et al. (1982), Brooker et al. (1988) and Hehlen et al. (2010).

The "galeener" equation is the exact one reported in Galeener and Sen (1978), which is a modification from Shuker and Gammon (1970) for accounting of (vo - v)^4 dependence of the Raman intensity. See also Brooker et al. (1988) for further discussion.

The "long" equation is that of Galeener and Sen (1978) corrected by a vo^3 coefficient for removing the cubic meter dimension of the equation of "galeener". This equation has been used in Mysen et al. (1982), Neuville and Mysen (1996) and Le Losq et al. (2012).

The "hehlen" equation is that reported in Hehlen et al. (2010). It actually originates before this publication (Brooker et al. 1988). It uses a different correction that avoid crushing the signal below 500 cm-1. THerefore, it has the advantage of keeping intact the Boson peak signal in glasses.
"""
function tlcorrection(data::Array{Float64},temp::Float64,wave::Float64;correction="long",normalisation="area",density=2210.0)

    h::Float64 = 6.626070040e-34   # J S    Plank constant from NIST
	hb::Float64 = 1.054571800e-34 # J S    Reduced Plank constant from NIST
    k::Float64 = 1.38064852e-23      # J K-1    Boltzman constant from NIST
    c::Float64 = 299792458.0         # M S-1    Speed of light from NIST
    nu0::Float64 = 1.0 ./wave .*1.0e9     # nu0 laser is in M-1 (wave is in nm)
    T::Float64 = temp + 273.15    # K temperature
	# density is in KG M-3

    x::Array{Float64} = data[:,1]
    y::Array{Float64} = data[:,2]

    # Calculate the error on data as sqrt(y). If y <= 0, then error = abs(y).
    ese::Array{Float64} = sqrt.(abs.(y))./abs.(y) # relative errors

	# get the Raman shift in m-1
	nu::Array{Float64} = 100.0.*x # cm-1 -> m-1 Raman shift

    # then we proceed to the correction
	if correction == "long"
		# Formula used in Mysen et al. (1982), Neuville and Mysen (1996) and Le Losq et al. (2012) (corrected for using the Planck constant in the last reference)
		# It is that reported in Brooker et al. (1988) with the addition of a scaling nu0^3 coefficient for adimentionality
    	frequency::Array{Float64} = nu0.^3.0.*nu./((nu0.-nu).^4) # frequency correction; dimensionless
    	boltzman::Array{Float64} = 1.0 .- exp.(-h.*c.*nu./(k.*T)) # temperature correction with Boltzman distribution; dimensionless
    	ycorr::Array{Float64} = y.*frequency.*boltzman; # correction

	elseif correction == "galeener"
		# This uses the formula reported in Galeener and Sen (1978) and Brooker et al. (1988); it uses the Bose-Einstein / Boltzman distribution
		# Formula from  without the scaling vo^3 coefficient reported in Mysen et al. (1982), Neuville and Mysen (1996) and Le Losq et al. (2012)
    	frequency = nu./((nu0.-nu).^4) # frequency correction; M^3
    	boltzman = 1.0 .- exp.(-h.*c.*nu./(k.*T)) # temperature correction with Boltzman distribution; dimensionless
    	ycorr = y.*frequency.*boltzman; # correction

	elseif correction =="hehlen"
		# this uses the formula reported in Hehlen et al. 2010
    	frequency = 1.0./(nu0.^3.0.*density) # frequency + density correction; M/KG
    	boltzman = 1.0 .- exp.(-h.*c.*nu./(k.*T)) # dimensionless
    	ycorr = nu.*y.*frequency.*boltzman; # correction
	else
		error("Not implemented, choose between long, galeener or hehlen.")
	end

	if normalisation == "area"
    	ycorr = ycorr./trapz(x,ycorr) # area normalisation
	elseif normalisation == "intensity"
		ycorr= ycorr./maximum(ycorr) # max. intensity normalisation
	elseif normalisation == "no"
		# nothing will happen
	else
		error("Set the optional normalisation parameter to area, intensity or no.")
	end

    esecorr::Array{Float64} = ese.*ycorr # error calculation

    return x, ycorr, esecorr
end


text:
################################################################################
#
#   BAND STRUCTURE CALCULATION
#
################################################################################
"""
    getBandstructure(
                unitcell::U,
                path::path
             [; resolution::Int64=-1,
                enforce_hermitian::Bool=false ]
            )

Calculates the band struture of a `Unitcell` object along some path given by a `Path` object and returns a `Bandstructure` object.
"""

function getBands(
                unitcell :: U,
                path :: RPA,
                bondInteractionMatrix:: HB,
                segment_resolution::Array{Int64, 1},
                kappa::Float64;
                enforce_hermitian::Bool=false,
                epsilon_degenerate::Float64=1e-6
            ) :: Bandstructure{RPA} where {
                 L,NB,S,B<:AbstractBond{L,NB},U<:AbstractUnitcell{S,B},HB<:AbstractBondHamiltonian{L,NS} where {
                 NS}, RPA<: AbstractReciprocalPath{P} where {
                 P <: AbstractReciprocalPoint{D} where {D}}}

   #get complete information of unitcell and couplings in a single object

    uc_bondInteractionMatrix = newHamiltonian(unitcell,bondInteractionMatrix)

    # build up the list of all bands of all segments
    segments_total    = Array{Array{Float64,1},1}[]

    # get the spin dimension

    spin_dimension = size(bondtermhopping(bondInteractionMatrix, bonds(unitcell)[rand(1:length(bonds(unitcell)))]))[1]

    # iterate over all path segments and push empty lists into the segments list
    for i in 1:length(segment_resolution)

        # build an empty band structure for this segment
        segment     = Array{Float64, 1}[]

        for b in 1:length(unitcell.sites)*spin_dimension
            push!(segment,     zeros(Float64, segment_resolution[i]))
           end

        # push the segment band structure into the complete segment list
        push!(segments_total,    segment)

    end

    # iterate over all segments
    for s in 1:length(segment_resolution)
        # get the grid in between two points
        resolution = segment_resolution[s]
        ss = (s+1)
        # get the local start and end point of the segment
        k1 = path.points[s].point
        k2 = path.points[ss].point
        kpoints  = getPointsOnLine(k1, k2, resolution)

        # calculate all energies
        for i in 1:resolution

            # get the current k
            k = kpoints[i]

            # get the interaction matrix for this k
            matrix = getMatrixAtK(uc_bondInteractionMatrix,k,kappa)
            # diagonalize the matrix
            eigenfactorization = eigen(matrix)
            eigenvalues  = (eigenfactorization.values)
            eigenvectors = eigenfactorization.vectors
            # save all the eigenvalues to their lists
            for b in 1:length(eigenvalues)
                segments_total[s][b][i] = eigenvalues[b]
            end

            ###########

        end
    end

    # generate a new band structure object
    bandstructure = newBandStructure(path, segments_total)

    # return the band structure
    return bandstructure
end
export getBands


text:
"""
Basic tests for the `Hsq` struct and associated function.
Tests all pass.
"""

include("../src/LDScoreJulia.jl")

using Main
using Main.LDScoreJulia
using Test

approx(x, y, eps) = abs(x - y) <= eps

function vector_approx(x, y, eps)
    if size(x) != size(y)
        return false
    end
    for (i, j) in enumerate(x)
        if !approx(x[i], y[i], eps)
            return false
        end
    end
    return true
end

@testset "Hsq struct" begin
    # column vectors d = 6 * 1
    y = [[0.2061] [0.2601] [4.3514] [6.1703] [5.0221] [2.418]]'
    x = [[6.1388] [6.4785] [8.88] [4.8064] [3.6219] [3.7244]]'
    w = [[6.1388] [6.4785] [8.88] [4.8064] [3.6219] [3.7244]]'
    N = [[592497.0] [591333.0] [696881.0] [617319.0] [687344.0] [685637.0]]'
    M = [[1173569.0]]
    n_blocks = 2
    intercept = nothing
    slow = false
    step1_ii = [[true] [true] [true] [true] [true] [true]]'
    old_weights = false

    hsq = LDScoreJulia.Hsq(
        y, x, w, N, M, n_blocks, intercept, slow, step1_ii, old_weights,
    )

    hsq.old_weights == false
    LDScoreJulia.ld_score_regression(
        hsq, y, x, w, N, M, n_blocks, intercept, slow, step1_ii, old_weights,
    )

    @test typeof(hsq) == LDScoreJulia.Hsq
    @test hsq.__null_intercept__ == 1
    @test hsq.y == y
    @test hsq.x == x
    @test hsq.w == w
    @test hsq.N == N
    @test hsq.M == M
    @test hsq.n_blocks == n_blocks
    @test hsq.intercept == intercept
    @test hsq.slow == slow
    @test hsq.step1_ii == step1_ii
end

@testset "Hsq functions" begin
    eps = 1.0e-6

    χ² = ones((4, 1)) .* 4
    ld = ones((4, 1))
    w_ld = ones((4, 1))
    N = 9 .* ones((4, 1))
    M = [[7]]
    n_blocks = 3
    intercept = 1
    slow = false
    twostep = nothing
    old_weights = false
    hsq = LDScoreJulia.Hsq(
        χ², ld, w_ld, N, M, n_blocks, intercept, slow, twostep, old_weights,
    )

    @testset "weights" begin
        x = [[0.09570313] [0.09570313] [0.09570313] [0.09570313]]'
        a = hsq_weights(ld, w_ld, N, M, 1)
        b = hsq_weights(ld, w_ld, N, M, 2)
        @test vector_approx(a, x, eps)
        @test vector_approx(a, b, eps)

        x = [[0.5] [0.5] [0.5] [0.5]]'
        a = hsq_weights(ld, w_ld, N, M, 0)
        b = hsq_weights(ld, w_ld, N, M, -1)
        @test vector_approx(a, x, eps)
        @test vector_approx(a, b, eps)
    end

    @testset "aggregate" begin
        χ² = ones((10, 1)) .* 3 ./ 2
        ld = ones((10, 1)) .* 100
        N = ones((10, 1)) .* 100000
        M = 1e7
        agg = aggregate(hsq, χ², ld, N, M)
        approx(agg, 0.5, eps)
        agg = aggregate(hsq, χ², ld, N, M; intercept=1.5)
        approx(agg, 0, eps)
    end
end


text:
mutable struct DMDALocalInfo{IT}
    dim::IT
    dof_per_node::IT
    stencil_width::IT
    global_size::NTuple{3, IT}
    local_start::NTuple{3, IT}
    local_size::NTuple{3, IT}
    ghosted_local_start::NTuple{3, IT}
    ghosted_local_size::NTuple{3, IT}
    boundary_type::NTuple{3, DMBoundaryType}
    stencil_type::DMDAStencilType
    ___padding___::NTuple{5, IT}
    DMDALocalInfo{IT}() where {IT} = new{IT}()
end

mutable struct DMDA{PetscLib} <: AbstractDM{PetscLib}
    ptr::CDM
    opts::Options{PetscLib}
    DMDA{PetscLib}(ptr, opts = Options(PetscLib)) where {PetscLib} =
        new{PetscLib}(ptr, opts)
end

"""
    empty(da::DMDA)

return an uninitialized `DMDA` struct.
"""
Base.empty(::DMDA{PetscLib}) where {PetscLib} = DMDA{PetscLib}(C_NULL)

"""
    DMDACreate1d(
        ::PetscLib
        comm::MPI.Comm,
        boundary_type::DMBoundaryType,
        global_dim,
        dof_per_node,
        stencil_width,
        points_per_proc::Union{Nothing, Vector{PetscInt}};
        dmsetfromoptions=true,
        dmsetup=true,
        options...
    )

Creates a 1-D distributed array with the options specified using keyword
arguments.

If keyword argument `dmsetfromoptions == true` then `setfromoptions!` called.
If keyword argument `dmsetup == true` then `setup!` is called.

# External Links
$(_doc_external("DMDA/DMDACreate1d"))
"""
function DMDACreate1d end

@for_petsc function DMDACreate1d(
    ::$UnionPetscLib,
    comm::MPI.Comm,
    boundary_type::DMBoundaryType,
    global_dim,
    dof_per_node,
    stencil_width,
    points_per_proc::Union{Nothing, Vector{$PetscInt}};
    dmsetfromoptions = true,
    dmsetup = true,
    options...,
)
    opts = Options($petsclib, options...)
    ref_points_per_proc = if isnothing(points_per_proc)
        C_NULL
    else
        @assert length(points_per_proc) == MPI.Comm_size(comm)
        points_per_proc
    end
    da = DMDA{$PetscLib}(C_NULL, opts)
    with(da.opts) do
        @chk ccall(
            (:DMDACreate1d, $petsc_library),
            PetscErrorCode,
            (
                MPI.MPI_Comm,
                DMBoundaryType,
                $PetscInt,
                $PetscInt,
                $PetscInt,
                Ptr{$PetscInt},
                Ptr{CDM},
            ),
            comm,
            boundary_type,
            global_dim,
            dof_per_node,
            stencil_width,
            ref_points_per_proc,
            da,
        )
    end
    dmsetfromoptions && setfromoptions!(da)
    dmsetup && setup!(da)
    # We can only let the garbage collect finalize when we do not need to
    # worry about MPI (since garbage collection is asyncronous)
    if comm == MPI.COMM_SELF
        finalizer(destroy, da)
    end
    return da
end

"""
    DMDACreate2d(
        ::PetscLib
        comm::MPI.Comm,
        boundary_type_x::DMBoundaryType,
        boundary_type_y::DMBoundaryType,
        stencil_type::DMDAStencilType,
        global_dim_x,
        global_dim_y,
        procs_x,
        procs_y,
        dof_per_node,
        stencil_width,
        points_per_proc_x::Union{Nothing, Vector{PetscInt}};
        points_per_proc_y::Union{Nothing, Vector{PetscInt}};
        dmsetfromoptions=true,
        dmsetup=true,
        options...
    )

Creates a 2-D distributed array with the options specified using keyword
arguments.

If keyword argument `dmsetfromoptions == true` then `setfromoptions!` called.
If keyword argument `dmsetup == true` then `setup!` is called.

# External Links
$(_doc_external("DMDA/DMDACreate2d"))
"""
function DMDACreate2d end

@for_petsc function DMDACreate2d(
    ::$UnionPetscLib,
    comm::MPI.Comm,
    boundary_type_x::DMBoundaryType,
    boundary_type_y::DMBoundaryType,
    stencil_type::DMDAStencilType,
    global_dim_x,
    global_dim_y,
    procs_x,
    procs_y,
    dof_per_node,
    stencil_width,
    points_per_proc_x::Union{Nothing, Vector{$PetscInt}},
    points_per_proc_y::Union{Nothing, Vector{$PetscInt}};
    dmsetfromoptions = true,
    dmsetup = true,
    options...,
)
    opts = Options($petsclib, options...)
    ref_points_per_proc_x = if isnothing(points_per_proc_x)
        C_NULL
    else
        @assert length(points_per_proc_x) == procs_x
        points_per_proc_x
    end
    ref_points_per_proc_y = if isnothing(points_per_proc_y)
        C_NULL
    else
        @assert length(points_per_proc_y) == procs_y
        points_per_proc_y
    end
    da = DMDA{$PetscLib}(C_NULL, opts)
    with(da.opts) do
        @chk ccall(
            (:DMDACreate2d, $petsc_library),
            PetscErrorCode,
            (
                MPI.MPI_Comm,
                DMBoundaryType,
                DMBoundaryType,
                DMDAStencilType,
                $PetscInt,
                $PetscInt,
                $PetscInt,
                $PetscInt,
                $PetscInt,
                $PetscInt,
                Ptr{$PetscInt},
                Ptr{$PetscInt},
                Ptr{CDM},
            ),
            comm,
            boundary_type_x,
            boundary_type_y,
            stencil_type,
            global_dim_x,
            global_dim_y,
            procs_x,
            procs_y,
            dof_per_node,
            stencil_width,
            ref_points_per_proc_x,
            ref_points_per_proc_y,
            da,
        )
    end
    dmsetfromoptions && setfromoptions!(da)
    dmsetup && setup!(da)
    # We can only let the garbage collect finalize when we do not need to
    # worry about MPI (since garbage collection is asyncronous)
    if comm == MPI.COMM_SELF
        finalizer(destroy, da)
    end
    return da
end

"""
    DMDACreate3d(
        ::PetscLib
        comm::MPI.Comm,
        boundary_type_x::DMBoundaryType,
        boundary_type_y::DMBoundaryType,
        boundary_type_z::DMBoundaryType,
        stencil_type::DMDAStencilType,
        global_dim_x,
        global_dim_y,
        global_dim_z,
        procs_x,
        procs_y,
        procs_z,
        global_dim_z,
        dof_per_node,
        stencil_width,
        points_per_proc_x::Union{Nothing, Vector{PetscInt}};
        points_per_proc_y::Union{Nothing, Vector{PetscInt}};
        points_per_proc_z::Union{Nothing, Vector{PetscInt}};
        dmsetfromoptions=true,
        dmsetup=true,
        options...
    )

Creates a 3-D distributed array with the options specified using keyword
arguments.

If keyword argument `dmsetfromoptions == true` then `setfromoptions!` called.
If keyword argument `dmsetup == true` then `setup!` is called.

# External Links
$(_doc_external("DMDA/DMDACreate3d"))
"""
function DMDACreate3d end

@for_petsc function DMDACreate3d(
    ::$UnionPetscLib,
    comm::MPI.Comm,
    boundary_type_x::DMBoundaryType,
    boundary_type_y::DMBoundaryType,
    boundary_type_z::DMBoundaryType,
    stencil_type::DMDAStencilType,
    global_dim_x,
    global_dim_y,
    global_dim_z,
    procs_x,
    procs_y,
    procs_z,
    dof_per_node,
    stencil_width,
    points_per_proc_x::Union{Nothing, Vector{$PetscInt}},
    points_per_proc_y::Union{Nothing, Vector{$PetscInt}},
    points_per_proc_z::Union{Nothing, Vector{$PetscInt}};
    dmsetfromoptions = true,
    dmsetup = true,
    options...,
)
    opts = Options($petsclib, options...)
    ref_points_per_proc_x = if isnothing(points_per_proc_x)
        C_NULL
    else
        @assert length(points_per_proc_x) == procs_x
        points_per_proc_x
    end
    ref_points_per_proc_y = if isnothing(points_per_proc_y)
        C_NULL
    else
        @assert length(points_per_proc_y) == procs_y
        points_per_proc_y
    end
    ref_points_per_proc_z = if isnothing(points_per_proc_z)
        C_NULL
    else
        @assert length(points_per_proc_z) == procs_z
        points_per_proc_z
    end
    da = DMDA{$PetscLib}(C_NULL, opts)
    with(da.opts) do
        @chk ccall(
            (:DMDACreate3d, $petsc_library),
            PetscErrorCode,
            (
                MPI.MPI_Comm,
                DMBoundaryType,
                DMBoundaryType,
                DMBoundaryType,
                DMDAStencilType,
                $PetscInt,
                $PetscInt,
                $PetscInt,
                $PetscInt,
                $PetscInt,
                $PetscInt,
                $PetscInt,
                $PetscInt,
                Ptr{$PetscInt},
                Ptr{$PetscInt},
                Ptr{$PetscInt},
                Ptr{CDM},
            ),
            comm,
            boundary_type_x,
            boundary_type_y,
            boundary_type_z,
            stencil_type,
            global_dim_x,
            global_dim_y,
            global_dim_z,
            procs_x,
            procs_y,
            procs_z,
            dof_per_node,
            stencil_width,
            ref_points_per_proc_x,
            ref_points_per_proc_y,
            ref_points_per_proc_z,
            da,
        )
    end
    dmsetfromoptions && setfromoptions!(da)
    dmsetup && setup!(da)
    # We can only let the garbage collect finalize when we do not need to
    # worry about MPI (since garbage collection is asyncronous)
    if comm == MPI.COMM_SELF
        finalizer(destroy, da)
    end
    return da
end

"""
    getinfo(da::DMDA)

Get the info associated with the distributed array `da`.

# External Links
$(_doc_external("DMDA/DMDAGetInfo"))
"""
getinfo(::DMDA)

@for_petsc function getinfo(da::DMDA{$PetscLib})
    dim = [$PetscInt(0)]
    glo_size = [$PetscInt(0), $PetscInt(0), $PetscInt(0)]
    procs_per_dim = [$PetscInt(0), $PetscInt(0), $PetscInt(0)]
    dof_per_node = [$PetscInt(0)]
    stencil_width = [$PetscInt(0)]
    boundary_type = [DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE]
    stencil_type = [DMDA_STENCIL_STAR]
    @chk ccall(
        (:DMDAGetInfo, $petsc_library),
        PetscErrorCode,
        (
            CDM,
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{DMBoundaryType},
            Ref{DMBoundaryType},
            Ref{DMBoundaryType},
            Ref{DMDAStencilType},
        ),
        da,
        dim,
        Ref(glo_size, 1),
        Ref(glo_size, 2),
        Ref(glo_size, 3),
        Ref(procs_per_dim, 1),
        Ref(procs_per_dim, 2),
        Ref(procs_per_dim, 3),
        dof_per_node,
        stencil_width,
        Ref(boundary_type, 1),
        Ref(boundary_type, 2),
        Ref(boundary_type, 3),
        stencil_type,
    )
    return (
        dim = dim[1],
        global_size = glo_size,
        procs_per_dim = procs_per_dim,
        dof_per_node = dof_per_node[1],
        boundary_type = boundary_type,
        stencil_width = stencil_width[1],
        stencil_type = stencil_type[1],
    )
end

"""
    getcorners(da::DMDA)

Returns a `NamedTuple` with the global indices (excluding ghost points) of the
`lower` and `upper` corners as well as the `size`.

# External Links
$(_doc_external("DMDA/DMDAGetCorners"))
"""
function getcorners(da::DMDA) end

@for_petsc function getcorners(da::DMDA{$PetscLib})
    info = DMDALocalInfo{$PetscInt}()
    corners = [$PetscInt(0), $PetscInt(0), $PetscInt(0)]
    local_size = [$PetscInt(0), $PetscInt(0), $PetscInt(0)]
    @chk ccall(
        (:DMDAGetCorners, $petsc_library),
        PetscErrorCode,
        (
            CDM,
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
        ),
        da,
        Ref(corners, 1),
        Ref(corners, 2),
        Ref(corners, 3),
        Ref(local_size, 1),
        Ref(local_size, 2),
        Ref(local_size, 3),
    )
    corners .+= 1
    return (
        lower = corners,
        upper = corners .+ local_size .- $PetscInt(1),
        size = local_size,
    )
end

"""
    getghostcorners(da::DMDA)

Returns a `NamedTuple` with the global indices (including ghost points) of the
`lower` and `upper` corners as well as the `size`.

# External Links
$(_doc_external("DMDA/DMDAGetGhostCorners"))
"""
function getghostcorners(da::DMDA) end

@for_petsc function getghostcorners(da::DMDA{$PetscLib})
    info = DMDALocalInfo{$PetscInt}()
    corners = [$PetscInt(0), $PetscInt(0), $PetscInt(0)]
    local_size = [$PetscInt(0), $PetscInt(0), $PetscInt(0)]
    @chk ccall(
        (:DMDAGetGhostCorners, $petsc_library),
        PetscErrorCode,
        (
            CDM,
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
            Ref{$PetscInt},
        ),
        da,
        Ref(corners, 1),
        Ref(corners, 2),
        Ref(corners, 3),
        Ref(local_size, 1),
        Ref(local_size, 2),
        Ref(local_size, 3),
    )
    corners .+= 1
    return (
        lower = corners,
        upper = corners .+ local_size .- 1,
        size = local_size,
    )
end

"""
    setuniformcoordinates!(
        da::DMDA
        xyzmin::NTuple{N, Real},
        xyzmax::NTuple{N, Real},
    ) where {N}

Set uniform coordinates for the `da` using the lower and upper corners defined
by the `NTuple`s `xyzmin` and `xyzmax`. If `N` is less than the dimension of the
`da` then the value of the trailing coordinates is set to `0`.

# External Links
$(_doc_external("DMDA/DMDASetUniformCoordinates"))
"""
function setuniformcoordinates! end

@for_petsc function setuniformcoordinates!(
    da::DMDA{$PetscLib},
    xyzmin::NTuple{N, Real},
    xyzmax::NTuple{N, Real},
) where {N}
    xmin = $PetscReal(xyzmin[1])
    xmax = $PetscReal(xyzmax[1])

    ymin = (N > 1) ? $PetscReal(xyzmin[2]) : $PetscReal(0)
    ymax = (N > 1) ? $PetscReal(xyzmax[2]) : $PetscReal(0)

    zmin = (N > 2) ? $PetscReal(xyzmin[3]) : $PetscReal(0)
    zmax = (N > 2) ? $PetscReal(xyzmax[3]) : $PetscReal(0)

    @chk ccall(
        (:DMDASetUniformCoordinates, $petsc_library),
        PetscErrorCode,
        (
            CDM,
            $PetscReal,
            $PetscReal,
            $PetscReal,
            $PetscReal,
            $PetscReal,
            $PetscReal,
        ),
        da,
        xmin,
        xmax,
        ymin,
        ymax,
        zmin,
        zmax,
    )
    return nothing
end


text:
struct NonlinearAdvection{𝒯} <: AbstractTerm end
struct LinearAdvection{𝒯} <: AbstractTerm end
struct VeryLinearAdvection{𝒯} <: AbstractTerm end

@inline calc_component!(flux, ::Nothing, _...) = nothing
@inline calc_component!(flux, ::AbstractTerm, _...) = nothing

@inline function calc_component!(flux, ::NonlinearAdvection{(:ρ, :ρu, :ρθ)}, state, aux, physics)
    ρ = state.ρ
    ρu = state.ρu
    ρθ = state.ρθ

    u = ρu / ρ

    flux.ρ += ρu
    flux.ρu += ρu ⊗ u
    flux.ρθ += ρθ * u

    nothing
end

@inline function calc_component!(flux, ::NonlinearAdvection{(:ρ, :ρu, :ρe)}, state, aux, physics)
    ρ = state.ρ
    ρu = state.ρu
    ρe = state.ρe
    ρq = state.ρq
    eos = physics.eos
    parameters = physics.parameters

    p = calc_pressure(eos, state, aux, parameters)
    u = ρu / ρ

    flux.ρ += ρu
    flux.ρu += ρu ⊗ u
    flux.ρe += (ρe + p) * u
    flux.ρq += ρq * u

    nothing
end

@inline function calc_component!(flux, ::LinearAdvection{(:ρ, :ρu, :ρe)}, state, aux, physics)
    ρu = state.ρu
    ρᵣ = aux.ref_state.ρ
    pᵣ = aux.ref_state.p
    ρeᵣ = aux.ref_state.ρe

    flux.ρ += ρu
    flux.ρe += (ρeᵣ + pᵣ) * ρu / ρᵣ

    nothing
end

@inline function calc_component!(flux, ::VeryLinearAdvection{(:ρ, :ρu, :ρe)}, state, aux, physics)
    # states
    ρ = state.ρ
    ρu = state.ρu
    ρe = state.ρe
    ρq = state.ρq

    # thermodynamics
    eos = physics.eos
    parameters = physics.parameters
    p = calc_very_linear_pressure(eos, state, aux, parameters)

    # Reference states
    ρᵣ = aux.ref_state.ρ
    ρuᵣ = aux.ref_state.ρu
    ρeᵣ = aux.ref_state.ρe
    ρqᵣ = aux.ref_state.ρq
    pᵣ = aux.ref_state.p

    # derived states
    u = ρu / ρᵣ - ρ * ρuᵣ / (ρᵣ^2)
    q = ρq / ρᵣ - ρ * ρqᵣ / (ρᵣ^2)
    e = ρe / ρᵣ - ρ * ρeᵣ / (ρᵣ^2)

    # derived reference states
    uᵣ = ρuᵣ / ρᵣ
    qᵣ = ρqᵣ / ρᵣ
    eᵣ = ρeᵣ / ρᵣ

    # can be simplified, but written this way to look like the VeryLinearKGVolumeFlux

    flux.ρ += ρᵣ * u + ρ * uᵣ # this is just ρu
    flux.ρu += p * I + ρᵣ .* (uᵣ .* u' + u .* uᵣ')
    flux.ρu += (ρ .* uᵣ) .* uᵣ'
    flux.ρe += (ρᵣ * eᵣ + pᵣ) * u
    flux.ρe += (ρᵣ * e + ρ * eᵣ + p) * uᵣ
    flux.ρq += ρᵣ * qᵣ * u + (ρᵣ * q + ρ * qᵣ) * uᵣ

    # flux.ρ  += ρu
    # flux.ρu += p * I
    # flux.ρe += (ρeᵣ + pᵣ) * ρu / ρᵣ

    nothing
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
using Test
using PotentialCalculation.UnitConversions

@testset "UnitConversions" begin
    e1 = 32.42
    e2 = 2.75E-4

    units = ["cm-1", "cm^-1", "ev", "eV", "kcal/mol", "kJ/mol", "kj/mol", "K"]

    for x in units
        t = energy_from(e1, x)
        @test t != e1
        @test e1 ≈ energy_to(t, x)
    end

    @test e2 ≈ change_energy_unit(e2, "eV", "eV")
    @test e2 == energy_from(e2, "hartree")
    @test e2 == energy_to(e2, "hartree")
end


text:
function distmesh2d(
    fd::Function,
    fh::Function,
    bbox,
    h0::Float64,
    pfix = [],
    ttol = 0.1,
    geps = 0.001 * h0,
    Fscale = 1.2,
    dptol = 0.001,
    deltat = 0.2,
    deps = sqrt(eps(Float64)) * h0,
)::Tuple{Array{Float64, 1}, Array{Float64, 1}}
    pold = Inf
    scaler = Scaler(bbox)
    v1, v2 = extractbox(bbox, h0, calculateh1(h0))
    x, y = meshgrid(v1, v2)
    shiftevenrows!(x, h0)
    p = buildinitialpoints(x,y,fd,scaler,geps)
    pfix = [vcat(scaledpoint(row, scaler)) for row in eachrow(pfix)]
    pfix = transpose(reshape(vcat(pfix...), 2, length(pfix)))
    while true
        del, vor, summ = deldir(p[:, 1], p[:, 2])
        trigs = triangles(del, summ)
        line_edges = validedges(trigs, del, scaler, fd, geps)
        pointstofvces = pointstoforces(line_edges, scaler, Fscale, pfix, fh)
        finalp, moveindex =
            finalpositions(pointstofvces, scaler, deltat, fd, geps, deps, h0)
        if moveindex < dptol
            x, y = plotedges(line_edges)
            break
        else
            p = finalp
        end
    end
    return [unscalex(x, scaler) for x in x], [unscaley(y, scaler) for y in y]
end

function buildinitialpoints(
    x::Array{Float64, 2},
    y::Array{Float64, 2},
    fd::Function,
    scaler::Scaler,
    geps::Float64
)::Array{Float64, 2}
    p = [x[:] y[:]]
    p = [vcat(row) for row in eachrow(p) if fd(row) < -geps]
    r0 = [1 ./ fh(row...) .^ 2 for row in p]
    r0max = maximum(r0 ./ maximum(r0))
    p = [vcat(scaledpoint(row, scaler)) for row in p if (rand(Float64, size(p)))[1] < r0max]
    p = transpose(reshape(vcat(p...), 2, length(p)))
    return p
end

function calculateh1(h0::Float64)::Float64
    return h0 * (sqrt(3) / 2)
end

function extractbox(
    bbox::Array{Float64,2},
    h0::Real,
    h1::Real
)::Tuple{Array{Float64, 1}, Array{Float64, 1}}
    v1 = bbox[1, 1]:h0:bbox[2, 1]
    v2 = bbox[1, 2]:h1:bbox[2, 2]
    return v1, v2
end

function shiftevenrows!(x::Array{Float64, 2}, h0::Float64)
    x[2:2:end, :] = x[2:2:end, :] .+ h0 / 2
end

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
function reduce_permuted(agg::OnlineStat, data, perm, rg)
    acc = deepcopy(agg)
    for i in rg
        fit!(acc, data[perm[i]])
    end
    return value(acc)
end

function reduce_permuted(agg, data, perm, rg)
    acc = data[perm[first(rg)]]
    for i in rg[2:end]
        acc = agg(acc, data[perm[i]])
    end
    return acc
end

function _reducer(args...; agg=Mean())
    key, data = StructArray(Base.front(args)), last(args)
    gp = GroupPerm(fast_sortable(key))
    perm = sortperm(gp)
    itr = (key[perm[first(rg)]] => reduce_permuted(agg, data, perm, rg) for rg in gp)
    keys, values = fieldarrays(StructArray(itr, unwrap = t -> t <: Tuple))
    return mapping(fieldarrays(keys)..., values)
end

"""
    reducer(args...; agg=Mean())

Reduce the last argument conditioned on the preceding ones using the online
statistic or binary function `agg`.
"""
const reducer = Analysis(_reducer)


text:
module PointsOfInterest

using AutoHashEquals, Base.Dates, IntervalSets

export POI

import Base: ∈, push!, delete!, empty!

@auto_hash_equals struct Video
    name::String
    time::DateTime
end

struct Metadata
    pois::Vector{String}
    files::Vector{Video}
end

const sides = (:left, :right)

@auto_hash_equals struct Point
    file::Int
    time::Millisecond

    function Point(md::Metadata, file::String, time::Millisecond)
        i = findfirst(x -> x.name == file, md.files)
        @assert i ≠ 0 "file not found in metadata"
        @assert time ≥ Millisecond(0) "negative times not allowed"
        new(i, time)
    end
end

Base.isless(p1::Point, p2::Point) = p1.file == p2.file ? p1.time < p2.time : p1.file < p2.file

@auto_hash_equals mutable struct POI
    name::Int
    interval::ClosedInterval{Point}
    comment::String

    function POI(md::Metadata, name::String, interval::ClosedInterval{Point}, comment::String)
        i = findfirst(md.pois, name)
        @assert i ≠ 0 "POI not found in metadata"
        @assert !isempty(interval) "interval is non-existent"
        new(i, interval, comment)
    end
end

POI(md::Metadata, name::String, start::Tuple{String, Millisecond}, stop::Tuple{String, Millisecond}, comment::String) = POI(md, name, Point(md, start...), Point(md, stop...), comment)

# POI(md::Metadata, name::String, label::String, start_file::String, start_time::Millisecond, stop_file::String, stop_time::Millisecond, comment::String) = POI(name, label, Point(md, start_file, start_time), Point(md, stop_file, stop_time), comment)

struct Log
    md::Metadata

    pois::Vector{POI}
end

Log(md::Metadata) = Log(md, POI[])

# findfirst

# findfirst(a::Log, x::Rep) = findfirst(a.reps, x)

# in

∈(x::POI, a::Log) = x ∈ a.pois

∈(x::Video, a::Log) = x ∈ a.files

# pushes

function push!(a::Log, x::Video)
    x ∈ a && return a
    i = findfirst(v -> v.time > x.time, a.md.files)
    insert!(a.md.files, i, x)
    for p in a.pois
        left = p.interval.left
        if left.file ≥ i
            right = p.interval.right
            p.interval = Point(left.file + 1, left.time)..Point(right.file + 1, right.time)
        end
    end
    return a
end

function push!(a::Log, x::POI)
    x ∉ a && push!(a.pois, x)
    return a
end

# deletes

function delete!(a::Log, x::POI)
    i = findfirst(a.pois, x)
    i ≠ 0 && deleteat!(a.pois, i)
    return a
end

function delete!(a::Log, x::Video)
    i = findfirst(a.md.files, x)
    i == 0 && return a
    deleteat!(a.md.files, i)
    filter!(p -> p.interval.left.file == i || p.interval.right.file == i, a.pois)
    for p in a.pois
        left = p.interval.left
        if left.file > i
            right = p.interval.right
            p.interval = Point(left.file - 1, left.time)..Point(right.file - 1, right.time)
        end
    end
    return a
end

# replace

function replace!(a::Log, o::POI, n::POI)
    o == n && return a
    n ∈ a && delete!(a, o)
    i = findfirst(a.pois, o)
    @assert i ≠ 0 "old POI not found"
    a.pois[i] = n
    return a
end

# empty

function empty!(a::Log, ::Type{POI})
    empty!(a.pois)
    return a
end

# pop

pop(md::Metadata, x::Point) = (md.files[x.file], x.time)

pop(md::Metadata, x::ClosedInterval{Point}) = (pop(md, x.left), pop(md, x.right))

pop(md::Metadata, x::POI) = (md.pois[x.name], pop(md, x.interval), x.comment)

# combine

function combine_mds(logs::Log...)
    pois = String[]
    videos = Video[]
    for log in logs
        union!(pois, log.md.pois)
        union!(videos, log.md.files)
    end
    sort!(videos, by = v -> v.time)
    return Metadata(pois, videos)
end

function combine(logs::Log...)
    md = combine_mds(logs)
    a = Log(md)
    for log in logs, p in log.pois
        push!(a, POI(md, pop(log.md, p)))
    end
    return a
end

end # module


text:
##################################################
# Useful general functions
##################################################
using StaticArrays # To have access to static arrays

using SpecialFunctions
using HypergeometricFunctions

##################################################
# Compute normalized variables
##################################################

"""
    _tr(r)

Computes the non-dimensional radius `x = r/b`.

# Arguments
- `r::Float64`: Radius.
"""
function _tr(r::Float64)
    return r/_b
end

"""
    _tE(E)

Computes the non-dimensional energy `tE = E/E0`.

# Arguments
- `E::Float64`: Energy.
"""
function _tE(E::Float64)
    return E/_E0
end

"""
    _tL(L)

Computes the non-dimensional angular momentum `tL = L/L0`.

# Arguments
- `L::Float64`: Angular momentum.
"""
function _tL(L::Float64)
    return L/_L0
end

##################################################
# Compute the normalized Plummer potential
##################################################

"""
    _tpsi(tr)

Computes the non-dimensional Plummer potential.

# Arguments
- `tr::Float64`: Non-dimensional radius.
"""
function _tpsi(tr::Float64)
    return 1/sqrt(1.0+tr^2)
end

"""
    _tdpsidr(tr)

Computes the derivative of the non-dimensional Plummer potential.

# Arguments
- `tr::Float64`: Non-dimensional radius.
"""
function _tdpsidr(tr::Float64)
    return -tr/(1.0+tr^2)^(3/2)
end

"""
    _td2psidr2(tr)

Computes the second derivative of the non-dimensional Plummer potential.

# Arguments
- `tr::Float64`: Non-dimensional radius.
"""
function _td2psidr2(tr::Float64)
    return -(1.0-tr^2)/(1.0+tr^2)^(5/2)
end

##################################################
# Compute the normalized effective Plummer potential
##################################################

"""
    _tpsiEff(tr,tL)

Computes the non-dimensional effective Plummer potential.

# Arguments
- `tr::Float64`: Non-dimensional radius.
- `tL::Float64`: Non-dimensional angular momentum.
"""
function _tpsiEff(tr::Float64, tL::Float64)
    return _tpsi(tr) - tL^2/(2*tr^2)
end

"""
    _tdpsiEffdr(tr,tL)

Computes the tr-derivative of the non-dimensional effective Plummer potential.

# Arguments
- `tr::Float64`: Non-dimensional radius.
- `tL::Float64`: Non-dimensional angular momentum.
"""
function _tdpsiEffdr(tr::Float64, tL::Float64)
    return _tdpsidr(tr) + tL^2/(tr^3)
end

"""
    _td2psiEffdr2(tr,tL)

Computes the second tr-derivative of the non-dimensional effective Plummer potential.

# Arguments
- `tr::Float64`: Non-dimensional radius.
- `tL::Float64`: Non-dimensional angular momentum.
"""
function _td2psiEffdr2(tr::Float64, tL::Float64)
    return _td2psidr2(tr) - 3* tL^2/(tr^4)
end

##################################################
# Compute the Plummer potential
##################################################

"""
    psi(r)

Computes the Plummer potential.

# Arguments
- `r::Float64`: Radius.
"""
function psi(r::Float64)
    return _E0*_tpsi(_tr(r))
end

"""
    dpsidr(r)

Computes the derivative of the Plummer potential.

# Arguments
- `r::Float64`: Radius.
"""
function dpsidr(r::Float64)
    return (_E0/_b)*_tdpsidr(_tr(r))
end

"""
    d2psidr2(r)

Computes the second derivative of the Plummer potential.

# Arguments
- `r::Float64`: Radius.
"""
function d2psidr2(r::Float64)
    return (_E0/_b^2)*_td2psidr2(_tr(r))
end

##################################################
# Compute the effective Plummer potential
##################################################

"""
    psiEff(r,L)

Computes the effective Plummer potential.

# Arguments
- `r::Float64`: Radius.
- `L::Float64`: Angular momentum.
"""
function psiEff(r::Float64, L::Float64)
    return _E0*_tpsiEff(_tr(r),_tL(L))
end

"""
    dpsiEffdr(r,L)

Computes the r-derivative of the effective Plummer potential.

# Arguments
- `r::Float64`: Radius.
- `L::Float64`: Angular momentum.
"""
function dpsiEffdr(r::Float64, L::Float64)
    return (_E0/_b)*_tdpsiEffdr(_tr(r),_tL(L))
end

"""
    d2psiEffdr2(r,L)

Computes the second r-derivative of the effective Plummer potential.

# Arguments
- `r::Float64`: Radius.
- `L::Float64`: Angular momentum.
"""
function d2psiEffdr2(r::Float64, L::Float64)
    return (_E0/_b^2)*_td2psiEffdr2(_tr(r),_tL(L))
end

##################################################
# Compute the Plummer density
##################################################

"""
    _trho(tr)

Computes the non-dimensional Plummer density.

# Arguments
- `tr::Float64`: Non-dimensional radius.
"""
function _trho(tr::Float64)
    return 1/(1+tr^2)^(5/2)
end

"""
    rho(r)

Computes the Plummer density.

# Arguments
- `r::Float64`: Radius.
"""
function rho(r::Float64)
    return _rho0*_trho(_tr(r))
end

##################################################
# Compute Lc(E) the angular momentum per unit mass of a circular orbit
##################################################

"""
    _tEta(s,tE)

Function used to compute the maximal, circular angular momentum of an orbit with non-dimensional energy `tE`.
Uses the radius parameter `s` defined by `s^2 = 1 + x^2` where `x = r/b`.

# Arguments
- `s::Float64`: Radius parameter.
- `tE::Float64`: Non-dimensional energy.
"""
function _tEta(s::Float64, tE::Float64)
    return (s^2-1)*(1/s - tE)
end

"""
    _sc(tE)

Compute the s-radius of the circular orbit with non-dimensional energy `tE`.
Uses the radius parameter `s` defined by `s^2 = 1 + x^2` where `x = r/b`.

# Remarks
- Uses an analytical expression (cubic root).

# Arguments
- `tE::Float64`: Non-dimensional energy.
"""
function _sc(tE::Float64)
    t1 = 1/(6*tE)
    t2 = (1+54*tE^2)/(216*tE^3)
    t3 = (1/(4*tE))*sqrt(1+1/(27*tE^2))
    return t1+cbrt(t2+t3)+cbrt(t2-t3)
end

"""
    _rc(tE)

Compute the radius of the circular orbit with non-dimensional energy `tE`.

# Remarks
- Uses an analytical expression (cubic root).

# Arguments
- `tE::Float64`: Non-dimensional energy.
"""
function _rc(tE::Float64)
    sc = _sc(tE)
    return _b*sqrt(sc^2-1)
end

"""
    Lc(tE)

Compute the angular momentum of the circular orbit with non-dimensional energy `tE`.

# Remarks
- Uses an analytical expression (cubic root).

# Arguments
- `tE::Float64`: Non-dimensional energy.
"""
function Lc(tE::Float64)
    if (tE == 1.0)
        return 0.0
    else
        return _L0*sqrt(2*_tEta(_sc(tE),tE))
    end
end

##################################################
# Convert (vr,vt) to (E,L) at a given radius r
##################################################

"""
    energy(r,vr,vt)

Compute the energy of a star at position `r` and velocity `(vr,vt)`.

# Arguments
- `r::Float64`: Radius.
- `vr::Float64`: Radial velocity.
- `vt::Float64`: Tangential velocity.
"""
function energy(r::Float64, vr::Float64, vt::Float64)
    return psi(r) + vr^2/(2.0) + vt^2/(2.0)
end

"""
    momentum(r,vr,vt)

Compute the angular momentum of a star at position `r` and velocity `(vr,vt)`.

# Arguments
- `r::Float64`: Radius.
- `vr::Float64`: Radial velocity.
- `vt::Float64`: Tangential velocity.
"""
function momentum(r::Float64, vr::Float64, vt::Float64)
    return r*vt
end

##################################################
# Convert (E,L) to (vr,vt) at a given radius r
# Convention vr >= 0
##################################################

"""
    radialVelocitySq(r,E,L)

Compute the square radial velocity of a star at position `r` of an orbit with parameters `(E,L)`.

# Arguments
- `r::Float64`: Radius.
- `E::Float64`: Energy.
- `L::Float64`: Angular momentum.
"""
function radialVelocitySq(r::Float64,E::Float64, L::Float64)
    return 2*(E-psiEff(r,L))
end

"""
    radialVelocity(r,E,L)

Compute the radial velocity of a star at position `r` of an orbit with parameters `(E,L)`.
Uses the convention `vr >= 0`.

# Arguments
- `r::Float64`: Radius.
- `E::Float64`: Energy.
- `L::Float64`: Angular momentum.
"""
function radialVelocity(r::Float64,E::Float64, L::Float64)
    return sqrt(2*abs((E-psiEff(r,L))))
end

"""
    tangentVelocity(r,E,L)

Compute the tangential velocity of a star at position `r` of an orbit with parameters `(E,L)`.

# Arguments
- `r::Float64`: Radius.
- `E::Float64`: Energy.
- `L::Float64`: Angular momentum.
"""
function tangentVelocity(r::Float64,E::Float64, L::Float64)
    return L/r
end

"""
    velocity(r,E,L)

Compute the velocity of a star at position `r` of an orbit with parameters `(E,L)`.

# Arguments
- `r::Float64`: Radius.
- `E::Float64`: Energy.
- `L::Float64`: Angular momentum.
"""
function velocity(r::Float64,E::Float64, L::Float64)
    return sqrt(2*abs((E-psi(r,L))))
end

##################################################
# Bisection algorithm to find zeroes of functions
# Here, used to compute Ec(L)
##################################################

"""
    bisection(fun, xl, xu, <keyword arguments>)

Algorithm which looks for zeroes of the functions `fun` given in argument.
It is a simple bisection algorithm, but it makes no allocations and is sufficiently fast.
It allows us not to have to use the Roots library that makes a lot of allocations.
The optional tolerances are set to the same as the ones found in Roots.jl, and are most likely overkill.

The zero is searched within the interval `[xl, xu]`, whose bounds are given as arguments.
"""
function bisection(fun, xl::Float64, xu::Float64, tolx::Float64=4.0*10^(-15), tolf::Float64=4.0*eps(Float64), iterMAX::Int64=50)
    if (xl > xu)
        xl, xu = xu, xl # Ordering the input arguments
    end
    #####
    fl, fu = fun(xl), fun(xu) # Evaluating the function on the bracket
    #####
    if (abs(fl) <= tolf) # We have already found a solution on the left bracket
        return xl # Returning the left bracket
    end
    #####
    if (abs(fu) <= tolf) # We have already found a solution on the right bracket
        return xu # Returning the right bracket
    end
    #####
    @assert fl*fu < 0.0 "bisection: NOT A BRACKET"
    #####
    iter = 0 # Counter for the iterations
    #####
    while true # Bisection loop
        #####
        xm = (xl+xu)*0.5 # Middle value
        #####
        if ((abs(xu-xl) <= tolx) || (iter > iterMAX)) # The considered bracket is smaller than the tolerance, or we have made too many iterations
            return xm # Returning the middle value
        end
        #####
        fm = fun(xm) # Value at the midpoint
        #####
        iter += 1 # Updating the counter of iterations
        #####
        if (abs(fm) <= tolf) # The middle value is below the threshold
            return xm # Returning the middle value
        end
        #####
        # Otherwise, we iterate the bisection
        if (fm*fl < 0.0) # One root can be found between the left point and the middle point
            xu, fu = xm, fm # The upper point becomes the midpoint
        else
            xl, fl = xm, fm # The lower point becomes the midpoint
        end
    end
end

"""
    _rc_L(L)

Compute the radius of the circular orbit with angular momentum `L`.

# Remarks
- Uses a bisection algorithm.

# Arguments
- `L::Float64`: Angular momentum.
"""
function _rc_L(L::Float64)
    tL = _tL(L)
    fct = (tr->_tdpsiEffdr(tr,tL))
    trmin = tL^(2/3) # derivative is positive
    trmax = trmin
    while (_tdpsiEffdr(trmax,tL) > 0) # find value of trmin such that derivative is negative
        trmax *= 2
    end

    troot = bisection(fct,trmin,trmax)

    return _b*troot
end

"""
    Ec(L)

Compute the energy of the circular orbit with angular momentum `L`.

# Remarks
- Uses a bisection algorithm.

# Arguments
- `L::Float64`: Angular momentum.
"""
function Ec(L::Float64)
    rc = _rc_L(L)
    return psiEff(rc,L)
end

max_stars_repo_name: KerwannTEP/CAT

max_issues_repo_name: KerwannTEP/CAT

max_forks_repo_name: KerwannTEP/CAT


text:
#=
Functions to reformat symbolic layer expressions for time steppers.
This is done in different ways depending on stepper and discretization
=#
function reformat_for_stepper(lhs, rhs, stepper)
    # rebuild the expressions depending on type of time stepper
    dt = symbols("dt");
    newlhs = [];
    newrhs = [];

    if typeof(rhs) <: Array && length(rhs) > 0 && typeof(rhs[1]) <: Array # recursively work on subarrays
        newlhs = copy(rhs);
        newrhs = copy(rhs);
        for vi=1:length(rhs)
            if length(lhs[1][vi]) > 0 # this dof has a time derivative term
                (newlhs[vi], newrhs[vi]) = reformat_for_stepper((lhs[1][vi], lhs[2][vi]), rhs[vi], stepper);
            else # no time derivative for this dof
                newlhs[vi] = lhs[2][vi];
                newrhs[vi] = rhs[vi];
            end
        end
    else
        # reformat depending on stepper type
        if stepper == EULER_IMPLICIT # lhs1 + dt*lhs2 = dt*rhs + lhs1
            for i=1:length(lhs[2])
                lhs[2][i] = lhs[2][i]*dt; # dt*lhs2
            end
            for i=1:length(rhs)
                rhs[i] = rhs[i]*dt; # dt*rhs
            end

            newlhs = copy(lhs[1]);
            append!(newlhs, lhs[2]); # lhs1 + dt*lhs2
            newrhs = copy(rhs);
            append!(newrhs, lhs[1]);# dt*rhs + lhs1

        elseif stepper == EULER_EXPLICIT || stepper == PECE # lhs1 = dt*rhs - dt*lhs2 + lhs1
            # # This version includes the lhs1 + dt*(  ) part
            # for i=1:length(lhs[2])
            #     lhs[2][i] = -lhs[2][i]*dt; # -dt*lhs2
            # end
            # for i=1:length(rhs)
            #     rhs[i] = rhs[i]*dt; # dt*rhs
            # end

            # newlhs = copy(lhs[1]);# lhs1
            # newrhs = copy(rhs);
            # append!(newrhs, lhs[2]);# dt*rhs - dt*lhs2
            # append!(newrhs, lhs[1]);# dt*rhs - dt*lhs2 + lhs1

            # This version doesn't include the lhs1 + dt*(  ) part
            for i=1:length(lhs[2])
                lhs[2][i] = -lhs[2][i]; # -lhs2
            end

            newlhs = copy(lhs[1]);# lhs1
            newrhs = copy(rhs);
            append!(newrhs, lhs[2]);# rhs - lhs2

        elseif stepper == CRANK_NICHOLSON # lhs1 + 0.5*dt*lhs2 = dt*rhs - 0.5*dt*lhs2 + lhs1
            lhs2l = copy(lhs[2]);
            lhs2r = copy(lhs[2]);
            for i=1:length(lhs[2])
                lhs2l[i] = Basic(0.5)*lhs2l[i]*dt; # 0.5*dt*lhs2
                lhs2r[i] = Basic(-0.5)*lhs2r[i]*dt; # -0.5*dt*lhs2
            end
            for i=1:length(rhs)
                rhs[i] = rhs[i]*dt; # dt*rhs
            end

            newlhs = copy(lhs[1]);
            append!(newlhs, lhs2l); # lhs1 + 0.5*dt*lhs2
            newrhs = copy(rhs);
            append!(newrhs, lhs[1]);# dt*rhs - 0.5*dt*lhs2 + lhs1
            append!(newrhs, lhs2r);

        elseif stepper == LSRK4 # (lhs1) : rhs - lhs2
            for i=1:length(lhs[2])
                lhs[2][i] = -lhs[2][i]; # -lhs2
            end
            # for i=1:length(rhs[1])
            #     rhs[1][i] = rhs[1][i]; # rhs
            # end

            newlhs = copy(lhs[1]);# lhs1
            newrhs = copy(rhs);
            append!(newrhs, lhs[2]);# rhs - lhs2

        elseif stepper == RK4 # (lhs1) : rhs - lhs2
            for i=1:length(lhs[2])
                lhs[2][i] = -lhs[2][i]; # -lhs2
            end
            # for i=1:length(rhs[1])
            #     rhs[1][i] = rhs[1][i]; # rhs
            # end

            newlhs = copy(lhs[1]);# lhs1
            newrhs = copy(rhs);
            append!(newrhs, lhs[2]);# rhs - lhs2

        end
    end

    return (newlhs, newrhs);
end

function reformat_for_stepper(lhs, rhs, face_lhs, face_rhs,stepper)
    # rebuild the expressions depending on type of time stepper
    dt = symbols("dt");
    newlhs = [];
    newrhs = [];
    newfacelhs = [];
    newfacerhs = [];

    if typeof(rhs) <: Array && length(rhs) > 0 && typeof(rhs[1]) <: Array # recursively work on subarrays
        newlhs = copy(rhs);
        newrhs = copy(rhs);
        newfacelhs = copy(face_rhs);
        newfacerhs = copy(face_rhs);
        for vi=1:length(rhs)
            if length(lhs[1][vi]) > 0 # this dof has a time derivative term
                (newlhs[vi], newrhs[vi], newfacelhs[vi], newfacerhs[vi]) = reformat_for_stepper((lhs[1][vi], lhs[2][vi]), rhs[vi], face_lhs[vi], face_rhs[vi], stepper);
            else # no time derivative for this dof
                newlhs[vi] = lhs[2][vi];
                newrhs[vi] = rhs[vi];
                newfacelhs[vi] = face_lhs[vi];
                newfacerhs[vi] = face_rhs[vi];
            end
        end
    else
        # reformat depending on stepper type
        if stepper == EULER_IMPLICIT # lhs1 + dt*lhs2 = dt*rhs + lhs1
            for i=1:length(lhs[2])
                lhs[2][i] = lhs[2][i]*dt; # dt*lhs2
            end
            for i=1:length(rhs)
                rhs[i] = rhs[i]*dt; # dt*rhs
            end
            for i=1:length(face_rhs)
                face_rhs[i] = face_rhs[i]*dt; # dt*facelhs
            end
            for i=1:length(face_lhs)
                face_lhs[i] = face_lhs[i]*dt; # dt*facerhs
            end

            newlhs = copy(lhs[1]);
            append!(newlhs, lhs[2]); # lhs1 + dt*lhs2
            newrhs = copy(rhs);
            append!(newrhs, lhs[1]);# dt*rhs + lhs1
            newfacerhs = copy(face_rhs);# dt*facerhs
            newfacelhs = copy(face_lhs);# dt*facelhs

        elseif stepper == EULER_EXPLICIT || stepper == PECE # lhs1 = lhs1 + dt*rhs - dt*lhs2
            # # This version includes the lhs1 + dt*(  ) part
            # for i=1:length(lhs[2])
            #     lhs[2][i] = -lhs[2][i]*dt; # -dt*lhs2
            # end
            # for i=1:length(rhs)
            #     rhs[i] = rhs[i]*dt; # dt*rhs
            # end
            # for i=1:length(face_rhs)
            #     face_rhs[i] = face_rhs[i]*dt; # dt*facerhs
            # end
            # for i=1:length(face_lhs)
            #     face_lhs[i] = -face_lhs[i]*dt; # -dt*facelhs
            # end

            # newlhs = copy(lhs[1]);# lhs1
            # newrhs = copy(rhs);
            # append!(newrhs, lhs[2]);# dt*rhs - dt*lhs2 + lhs1
            # #append!(newrhs, lhs[1]);
            # newfacerhs = copy(face_rhs);
            # append!(newfacerhs, face_lhs);# dt*facerhs - dt*facelhs
            # newfacelhs = [];

            # This version does not include the lhs1 + dt*(  ) part
            for i=1:length(lhs[2])
                lhs[2][i] = -lhs[2][i]; # -lhs2
            end
            for i=1:length(face_lhs)
                face_lhs[i] = -face_lhs[i]; # -facelhs
            end

            newlhs = copy(lhs[1]);# lhs1
            newrhs = copy(rhs);
            append!(newrhs, lhs[2]);# rhs - lhs2
            newfacerhs = copy(face_rhs);
            append!(newfacerhs, face_lhs);# facerhs - facelhs
            newfacelhs = [];

        elseif stepper == CRANK_NICHOLSON # lhs1 + 0.5*dt*lhs2 = dt*rhs - 0.5*dt*lhs2 + lhs1
            lhs2l = copy(lhs[2]);
            lhs2r = copy(lhs[2]);
            facelhs2l = copy(face_lhs);
            facelhs2r = copy(face_lhs);
            for i=1:length(lhs[2])
                lhs2l[i] = Basic(0.5)*lhs2l[i]*dt; # 0.5*dt*lhs2
                lhs2r[i] = Basic(-0.5)*lhs2r[i]*dt; # -0.5*dt*lhs2
            end
            for i=1:length(face_lhs)
                facelhs2l[i] = Basic(0.5)*facelhs2l[i]*dt; # 0.5*dt*lhs2
                facelhs2r[i] = Basic(-0.5)*facelhs2r[i]*dt; # -0.5*dt*lhs2
            end
            for i=1:length(rhs)
                rhs[i] = rhs[i]*dt; # dt*rhs
            end
            for i=1:length(face_rhs)
                face_rhs[i] = face_rhs[i]*dt; # dt*facelhs
            end

            newlhs = copy(lhs[1]);
            append!(newlhs, lhs2l); # lhs1 + 0.5*dt*lhs2
            newrhs = copy(rhs);
            append!(newrhs, lhs[1]);# dt*rhs - 0.5*dt*lhs2 + lhs1
            append!(newrhs, lhs2r);
            newfacelhs = facelhs2l;
            newfacerhs = face_rhs[1];
            append!(newfacerhs, facelhs2r);

        elseif stepper == LSRK4 # (lhs1) : rhs - lhs2
            for i=1:length(lhs[2])
                lhs[2][i] = -lhs[2][i]; # -lhs2
            end
            for i=1:length(face_lhs)
                face_lhs[i] = -face_lhs[i]; # -facelhs
            end

            newlhs = copy(lhs[1]);# lhs1
            newrhs = copy(rhs);
            append!(newrhs, lhs[2]);# rhs - lhs2
            newfacerhs = copy(face_rhs);
            append!(newfacerhs, face_lhs);# facerhs - facelhs
            newfacelhs = [];

        elseif stepper == RK4 # (lhs1) : rhs - lhs2
            for i=1:length(lhs[2])
                lhs[2][i] = -lhs[2][i]; # -lhs2
            end
            for i=1:length(face_lhs)
                face_lhs[i] = -face_lhs[i]; # -facelhs
            end

            newlhs = copy(lhs[1]);# lhs1
            newrhs = copy(rhs);
            append!(newrhs, lhs[2]);# rhs - lhs2
            newfacerhs = copy(face_rhs);
            append!(newfacerhs, face_lhs);# facerhs - facelhs
            newfacelhs = [];

        end
    end

    return (newlhs, newrhs, newfacelhs, newfacerhs);
end

# Special version for FV. Assumes a Dt(u) term that is not explicitly included.
function reformat_for_stepper_fv(flhs, frhs, slhs, srhs, stepper)
    (newflhs, newfrhs) = reformat_for_stepper_fv_flux(flhs, frhs, stepper);
    (newslhs, newsrhs) = reformat_for_stepper_fv_source(slhs, srhs, stepper);

    return (newflhs, newfrhs, newslhs, newsrhs);
end

# Special version for FV. Flux term only
function reformat_for_stepper_fv_flux(lhs, rhs, stepper)
    # rebuild the expressions depending on type of time stepper
    dt = symbols("dt");
    newlhs = [];
    newrhs = [];

    if typeof(rhs) <: Array && length(rhs) > 0 && typeof(rhs[1]) <: Array # recursively work on subarrays
        newlhs = copy(rhs); # copy just to set up the arrays right
        newrhs = copy(rhs);
        for vi=1:length(rhs);
            (newlhs[vi], newrhs[vi]) = reformat_for_stepper_fv_flux(lhs[vi], rhs[vi], stepper);
        end
    else
        # reformat depending on stepper type
        if stepper == EULER_EXPLICIT || stepper == RK4 || stepper == LSRK4 || stepper == PECE
            # du/dt = -lhs + rhs
            for i=1:length(lhs)
                lhs[i] = -lhs[i];
            end
            # for i=1:length(rhs)
            #     rhs[i] = -rhs[i]; # don't change sign for rhs, because already done by parser
            # end

            newlhs = []; # put everything in rhs for explicit steppers
            newrhs = copy(rhs);
            append!(newrhs, lhs);

        elseif stepper == EULER_IMPLICIT
            # du/dt + lhs = rhs  ->  u(+) + dt*lhs = u(-) + dt*rhs
            # Do the *dt part here, but adding u will be done by the solver
            for i=1:length(lhs)
                lhs[i] = lhs[i]*dt;
            end
            for i=1:length(rhs)
                rhs[i] = rhs[i]*dt;
            end

            newlhs = copy(lhs);
            newrhs = copy(rhs);

        elseif stepper == CRANK_NICHOLSON
            printerr("TODO: crank nicholson for FV. Choose something else, please.", fatal=true);
        end
    end

    return (newlhs, newrhs);
end

# Special version for FV. Assumes a Dt(u) term that is not explicitly included.
function reformat_for_stepper_fv_source(lhs, rhs, stepper)
    # rebuild the expressions depending on type of time stepper
    dt = symbols("dt");
    newlhs = [];
    newrhs = [];

    if typeof(rhs) <: Array && length(rhs) > 0 && typeof(rhs[1]) <: Array # recursively work on subarrays
        newlhs = copy(rhs); # copy just to set up the arrays right
        newrhs = copy(rhs);
        for vi=1:length(rhs);
            (newlhs[vi], newrhs[vi]) = reformat_for_stepper_fv_source(lhs[vi], rhs[vi], stepper);
        end
    else
        # reformat depending on stepper type
        if stepper == EULER_EXPLICIT || stepper == RK4 || stepper == LSRK4 || stepper == PECE
            # du/dt = -lhs + rhs
            # for i=1:length(lhs)
            #     lhs[i] = -lhs[i]; # don't change sign here because it is done by parser?
            # end
            for i=1:length(rhs)
                rhs[i] = -rhs[i]; # need to change sign here?
            end

            newlhs = []; # put everything in rhs for explicit steppers
            newrhs = copy(rhs);
            append!(newrhs, lhs);

        elseif stepper == EULER_IMPLICIT
            # du/dt + lhs = rhs  ->  u(+) + dt*lhs = u(-) + dt*rhs
            # Do the *dt part here, but adding u will be done by the solver
            for i=1:length(lhs)
                lhs[i] = -lhs[i]*dt;
            end
            for i=1:length(rhs)
                rhs[i] = -rhs[i]*dt;
            end

            newlhs = copy(lhs);
            newrhs = copy(rhs);

        elseif stepper == CRANK_NICHOLSON
            printerr("TODO: crank nicholson for FV. Choose something else, please.", fatal=true);
        end
    end

    return (newlhs, newrhs);
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
module Inflate

# package code goes here
export decompress, gunzip
include(Pkg.dir("Inflate","src","inflater.jl"))
end # module


text:
using Test
using QuantumControlBase.Shapes

@testset "flattop-blackman" begin
    shape(v) = flattop(v, t₀=10, T=20, t_rise=2, func=:blackman)
    @test shape(9.9) == 0
    @test shape(10) < 1e-14
    @test shape(20) < 1e-14
    @test shape(20.1) == 0
    @test shape(15) == 1

    tlist = collect(range(0, 30, step=1.0))
    shape_vec = flattop.(tlist, t₀=10, T=20, t_rise=2, func=:blackman)
    @test shape_vec[11] == shape(10)
end

@testset "flattop-sinsq" begin
    shape(v) = flattop(v, t₀=10, T=20, t_rise=2, func=:sinsq)
    @test shape(9.9) == 0
    @test shape(10) < 1e-14
    @test shape(20) < 1e-14
    @test shape(20.1) == 0
    @test shape(15) == 1

    tlist = collect(range(0, 30, step=1.0))
    shape_vec = flattop.(tlist, t₀=10, T=20, t_rise=2, func=:sinsq)
    @test shape_vec[11] == shape(10)
end