text:
# Format the code and notebooks in this repository.
# This can be useful to make things more consistent.

using Pkg

Pkg.activate(; temp=true)
package = PackageSpec(; name="JuliaFormatter", version="0.22")
Pkg.add(package)

using JuliaFormatter
format(".")

max_stars_repo_name: VEZY/PlantBiophysics-paper

max_issues_repo_name: VEZY/PlantBiophysics-paper

max_forks_repo_name: rikhuijzer/JuliaTutorialsTemplate


text:
using MCPhylo
using Test

start_tree = MCPhylo.parsing_newick_string("((A,(B,(C,(D,E)))),(F,(G,H)))")
descending_tree = MCPhylo.parsing_newick_string("(((((D,E),C),B),A),((G,H),F))")
ascending_tree = MCPhylo.parsing_newick_string("((F,(G,H)),(A,(B,(C,(D,E)))))")

start_newick = newick(start_tree)
desc_newick = newick(descending_tree)
asc_newick = newick(ascending_tree)

ladderized_asc_tree = MCPhylo.ladderize_tree(start_tree)
@test newick(ladderized_asc_tree) == asc_newick

ladderized_desc_tree = MCPhylo.ladderize_tree(start_tree,false)
@test newick(ladderized_desc_tree) == desc_newick

MCPhylo.ladderize_tree!(start_tree)
@test newick(start_tree) == asc_newick

start_tree = MCPhylo.parsing_newick_string("((A,(B,(C,(D,E)))),(F,(G,H)))")
MCPhylo.ladderize_tree!(start_tree, false)
@test newick(start_tree) == desc_newick

ext: jl


text:
module ParallelRun

using Base.Threads
using Base.Dates

import FITSIO
import JLD
import WCS

import ..Config
import ..Log
using ..Model
import ..SDSSIO
import ..Infer
import ..SDSSIO: RunCamcolField, IOStrategy, PlainFITSStrategy
import ..PSF
import ..SEP
import ..Coordinates: angular_separation, match_coordinates

import ..DeterministicVI: infer_source

include("joint_infer.jl")

abstract type ParallelismStrategy; end
struct ThreadsStrategy <: ParallelismStrategy; end

# In production mode, rather the development mode, always catch exceptions
const is_production_run = haskey(ENV, "CELESTE_PROD") && ENV["CELESTE_PROD"] != ""

#
# ------ bounding box ------
struct BoundingBox
    ramin::Float64
    ramax::Float64
    decmin::Float64
    decmax::Float64
end

function BoundingBox(ramin::String, ramax::String, decmin::String, decmax::String)
    BoundingBox(parse(Float64, ramin),
                parse(Float64, ramax),
                parse(Float64, decmin),
                parse(Float64, decmax))
end

# ------
# to time parts of Celeste
mutable struct InferTiming
    query_fids::Float64
    read_photoobj::Float64
    read_img::Float64
    preload_rcfs::Float64
    find_neigh::Float64
    load_wait::Float64
    proc_wait::Float64
    init_elbo::Float64
    opt_srcs::Float64
    num_srcs::Int64
    sched_ovh::Float64
    load_imba::Float64
    ga_get::Float64
    ga_put::Float64
    store_res::Float64
    write_results::Float64
    wait_done::Float64

    InferTiming() = new(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0, 0.0, 0.0,
                        0.0, 0.0, 0.0, 0.0, 0.0)
end

function add_timing!(i::InferTiming, j::InferTiming)
    i.query_fids = i.query_fids + j.query_fids
    i.read_photoobj = i.read_photoobj + j.read_photoobj
    i.read_img = i.read_img + j.read_img
    i.preload_rcfs = i.preload_rcfs + j.preload_rcfs
    i.find_neigh = i.find_neigh + j.find_neigh
    i.load_wait = i.load_wait + j.load_wait
    i.proc_wait = i.proc_wait + j.proc_wait
    i.init_elbo = i.init_elbo + j.init_elbo
    i.opt_srcs = i.opt_srcs + j.opt_srcs
    i.num_srcs = i.num_srcs + j.num_srcs
    i.sched_ovh = i.sched_ovh + j.sched_ovh
    i.load_imba = i.load_imba + j.load_imba
    i.ga_get = i.ga_get + j.ga_get
    i.ga_put = i.ga_put + j.ga_put
    i.store_res = i.store_res + j.store_res
    i.write_results = i.write_results + j.write_results
    i.wait_done = i.wait_done + j.wait_done
end

function puts_timing(i::InferTiming)
    i.num_srcs = max(1, i.num_srcs)
    Log.message("timing: query_fids=$(i.query_fids)")
    Log.message("timing: read_photoobj=$(i.read_photoobj)")
    Log.message("timing: read_img=$(i.read_img)")
    Log.message("timing: preload_rcfs=$(i.preload_rcfs)")
    Log.message("timing: find_neigh=$(i.find_neigh)")
    #Log.message("timing: load_wait=$(i.load_wait)")
    #Log.message("timing: proc_wait=$(i.proc_wait)")
    Log.message("timing: init_elbo=$(i.init_elbo)")
    Log.message("timing: opt_srcs=$(i.opt_srcs)")
    Log.message("timing: num_srcs=$(i.num_srcs)")
    #Log.message("timing: average opt_srcs=$(i.opt_srcs/i.num_srcs)")
    #Log.message("timing: sched_ovh=$(i.sched_ovh)")
    Log.message("timing: load_imba=$(i.load_imba)")
    #Log.message("timing: ga_get=$(i.ga_get)")
    #Log.message("timing: ga_put=$(i.ga_put)")
    #Log.message("timing: store_res=$(i.store_res)")
    Log.message("timing: write_results=$(i.write_results)")
    #Log.message("timing: wait_done=$(i.wait_done)")
end

function time_puts(elapsedtime, bytes, gctime, allocs)
    s = @sprintf("timing: total=%10.6f seconds", elapsedtime/1e9)
    if bytes != 0 || allocs != 0
        bytes, mb = Base.prettyprint_getunits(bytes, length(Base._mem_units),
                            Int64(1024))
        allocs, ma = Base.prettyprint_getunits(allocs, length(Base._cnt_units),
                            Int64(1000))
        if ma == 1
            s = string(s, @sprintf(" (%d%s allocation%s: ", allocs,
                                   Base._cnt_units[ma], allocs==1 ? "" : "s"))
        else
            s = string(s, @sprintf(" (%.2f%s allocations: ", allocs,
                                   Base._cnt_units[ma]))
        end
        if mb == 1
            s = string(s, @sprintf("%d %s%s", bytes,
                                   Base._mem_units[mb], bytes==1 ? "" : "s"))
        else
            s = string(s, @sprintf("%.3f %s", bytes, Base._mem_units[mb]))
        end
        if gctime > 0
            s = string(s, @sprintf(", %.2f%% gc time", 100*gctime/elapsedtime))
        end
        s = string(s, ")")
    elseif gctime > 0
        s = string(s, @sprintf(", %.2f%% gc time", 100*gctime/elapsedtime))
    end
    Log.message(s)
end

# ------
# initialization helpers

"""
    detect_sources(images)

Detect sources in a set of (possibly overlapping) `RawImage`s and combine
duplicates.

# Returns

- `catalog::Vector{CatalogEntry}`: Detected sources.
- `source_radii::Vector{Float64}`: Radius of circle containing all of each
  source's member pixels (degrees).

# Notes

Here, we're dealing with the `RawImage` type, which is SDSS-specific. The
pixel values are already sky subtracted (and calibrated), so we don't have to
background subtract. However, we still run a background analysis, just to
determine the rough image noise for the purposes of thresholding.

This could be slightly suboptimal in terms of
selecting faint sources: If the image noise varies significantly
across the image, the threshold might be too high or too low in some
places. However, we don't really trust the variable background RMS without
first masking sources. (We could add this.)
"""
function detect_sources(images::Vector{SDSSIO.RawImage})

    catalog = Vector{CatalogEntry}()
    source_radii = Vector{Float64}()

    for image in images

        # Run background, just to get background rms.
        bkg = SEP.Background(image.pixels; boxsize=(256, 256),
                             filtersize=(3, 3))
        sep_catalog = SEP.extract(image.pixels, 1.3; noise=SEP.global_rms(bkg))

        # convert pixel coordinates to world coordinates
        pixcoords = Array{Float64}(2, length(sep_catalog.x))
        for i in eachindex(sep_catalog.x)
            pixcoords[1, i] = sep_catalog.x[i]
            pixcoords[2, i] = sep_catalog.y[i]
        end
        worldcoords = WCS.pix_to_world(image.wcs, pixcoords)

        # Get angle offset between +RA axis and +x axis from the
        # image's WCS.  This assumes there is no skew, meaning the x
        # and y axes are perpindicular in world coordinates.
        cd = image.wcs[:cd]
        sgn = sign(det(cd))
        n_vs_y_rot = atan2(sgn * cd[1, 2],  sgn * cd[1, 1])  # angle of N CCW
                                                             # from +y axis
        x_vs_n_rot = -(n_vs_y_rot + pi/2)  # angle of +x CCW from N

        # convert sep_catalog entries to CatalogEntries
        im_catalog = Vector{CatalogEntry}(length(sep_catalog.npix))
        im_source_radii = Vector{Float64}(length(sep_catalog.npix))
        for i in eachindex(sep_catalog.npix)

            # For galaxy flux, use sum of all pixels above
            # threshold. Set other band fluxes to NaN, to indicate which
            # band the flux is measured in. We'll use this information below
            # when merging detections in different bands.
            gal_fluxes = fill(NaN, NUM_BANDS)
            gal_fluxes[image.b] = sep_catalog.flux[i]

            gal_ab = sep_catalog.a[i] / sep_catalog.b[i]

            # SEP angle is CCW from +x axis and in [-pi/2, pi/2].
            # Add offset of x axis from N to make angle CCW from N
            gal_angle = sep_catalog.theta[i] + x_vs_n_rot

            # A 2-d symmetric gaussian has CDF(r) =  1 - exp(-(r/sigma)^2/2)
            # The half-light radius is then r = sigma * sqrt(2 ln(2))
            sigma = sqrt(sep_catalog.a[i] * sep_catalog.b[i])
            gal_scale = sigma * sqrt(2. * log(2.))

            pos = worldcoords[:, i]
            im_catalog[i] = CatalogEntry(pos,
                                         false,  # is_star
                                         Float64[],  # will replace below
                                         gal_fluxes,
                                         0.5,  # gal_frac_dev
                                         gal_ab,
                                         gal_angle,
                                         gal_scale)

            # get object extent in degrees from bounding box in pixels
            xmin = sep_catalog.xmin[i]
            xmax = sep_catalog.xmax[i]
            ymin = sep_catalog.ymin[i]
            ymax = sep_catalog.ymax[i]
            corner_pixcoords = Float64[xmin xmin xmax xmax;
                                       ymin ymax ymin ymax]
            corners = WCS.pix_to_world(image.wcs, corner_pixcoords)
            im_source_radii[i] =
                maximum(angular_separation(pos[1], pos[2],
                                           corners[1, j], corners[2, j])
                        for j in 1:4)
        end

        # Combine the catalog for this image with the joined catalog
        if length(catalog) == 0
            catalog = im_catalog
            source_radii = im_source_radii
        else
            # for each detection in image, find nearest match in joined catalog
            idx, dist = match_coordinates([ce.pos[1] for ce in im_catalog],
                                          [ce.pos[2] for ce in im_catalog],
                                          [ce.pos[1] for ce in catalog],
                                          [ce.pos[2] for ce in catalog])

            for (i, ce) in enumerate(im_catalog)
                # if there is an existing object within 1 arcsec,
                # add it to the joined catalog
                if dist[i] < (1.0 / 3600.)
                    existing_ce = catalog[idx[i]]
                    if isnan(existing_ce.gal_fluxes[image.b])
                        existing_ce.gal_fluxes[image.b] =
                            ce.gal_fluxes[image.b]
                    end
                    if im_source_radii[i] > source_radii[idx[i]]
                        source_radii[idx[i]] = im_source_radii[i]
                    end

                # if not, add entry to joined catalog
                else
                    push!(catalog, ce)
                    push!(source_radii, im_source_radii[i])
                end
            end
        end
    end  # loop over images

    # clean up NaN flux entries in joined catalog,
    # copy gal_fluxes to star_fluxes
    for ce in catalog
        minflux = minimum(f for f in ce.gal_fluxes if !isnan(f))
        for i in eachindex(ce.gal_fluxes)
            if isnan(ce.gal_fluxes[i])
                ce.gal_fluxes[i] = minflux
            end
        end
        ce.star_fluxes = copy(ce.gal_fluxes)
    end

    return catalog, source_radii
end

"""
New version of infer_init() that uses detect_sources() instead of SDSS
photoObj catalog for initialization.
"""
function infer_init_new(rcfs::Vector{RunCamcolField},
                        strategy::SDSSIO.IOStrategy;
                        box=BoundingBox(-1000., 1000., -1000., 1000.),
                        timing=InferTiming())

    # Initialize variables to empty vectors in case try block fails
    catalog = CatalogEntry[]
    source_radii = Float64[]
    target_sources = Int[]
    images = Image[]

    try
        # Read in images for all RCFs
        tic()
        raw_images = SDSSIO.load_raw_images(strategy, rcfs)
        timing.read_img += toq()

        # detect sources on all raw images (before background added back)
        catalog, source_radii = detect_sources(raw_images)

        # Get indices of entries in the RA/Dec range of interest.
        # (Some images can have regions that are outside the box, so not
        # all sources are necessarily in the box.)
        entry_in_range = entry->((box.ramin < entry.pos[1] < box.ramax) &&
                                 (box.decmin < entry.pos[2] < box.decmax))
        target_sources = find(entry_in_range, catalog)

        Log.info("$(Time(now())): $(length(catalog)) primary sources, ",
                 "$(length(target_sources)) target sources in $(box.ramin), ",
                 "$(box.ramax), $(box.decmin), $(box.decmax)")

        # convert raw images to images
        try
            images = [convert(Image, raw_image) for raw_image in raw_images]
        catch exc
            Log.exception(exc)
            rethrow()
        end

    catch ex
        Log.exception(ex)
        empty!(target_sources)
        rethrow()
    end

    # build neighbor map based on source radii
    tic()
    neighbor_map = Vector{Int64}[Int64[] for s in target_sources]
    for ts in 1:length(target_sources)
        s = target_sources[ts]
        ce = catalog[s]
        ce_rad = source_radii[s]

        for s2 in 1:length(catalog)
            s2 == s && continue
            ce2 = catalog[s2]
            dist = angular_separation(ce.pos[1], ce.pos[2],
                                      ce2.pos[1], ce2.pos[2])

            if dist < source_radii[s] + source_radii[s2]
                push!(neighbor_map[ts], s2)
            end
        end
    end
    timing.find_neigh += toq()

    return catalog, target_sources, neighbor_map, images
end

"""
Given a list of RCFs, load the catalogs, determine the target sources,
load the images, and build the neighbor map.
"""
function infer_init(rcfs::Vector{RunCamcolField},
                    strategy::SDSSIO.IOStrategy;
                    box=BoundingBox(-1000., 1000., -1000., 1000.),
                    primary_initialization=true,
                    timing=InferTiming())
    catalog = Vector{CatalogEntry}()
    target_sources = Vector{Int}()
    neighbor_map = Vector{Vector{Int}}()
    images = Vector{Image}()

    # Read all primary objects in these fields.
    duplicate_policy = primary_initialization ? :primary : :first
    tic()
    timing.preload_rcfs += toq()

    tic()
    for rcf in rcfs
        this_cat = SDSSIO.read_photoobj_files(strategy, [rcf],
                                              duplicate_policy=duplicate_policy)
        append!(catalog, this_cat)
    end
    timing.read_photoobj += toq()

    # Get indices of entries in the RA/Dec range of interest.
    entry_in_range = entry->((box.ramin < entry.pos[1] < box.ramax) &&
                             (box.decmin < entry.pos[2] < box.decmax))
    target_sources = find(entry_in_range, catalog)

    Log.info("$(Time(now())): $(length(catalog)) primary sources, ",
             "$(length(target_sources)) target sources in $(box.ramin), ",
             "$(box.ramax), $(box.decmin), $(box.decmax)")

    # Load images and neighbor map for target sources
    if length(target_sources) > 0
        # Read in images for all (run, camcol, field).
        try
            tic()
            images = SDSSIO.load_field_images(strategy, rcfs)
            timing.read_img += toq()
        catch ex
            Log.exception(ex)
            empty!(target_sources)
        end

        tic()
        neighbor_map = Infer.find_neighbors(target_sources, catalog, images)
        timing.find_neigh += toq()
    end

    return catalog, target_sources, neighbor_map, images
end

# ------
# optimization result container

struct OptimizedSource
    init_ra::Float64
    init_dec::Float64
    vs::Vector{Float64}
end

function serialize(s::Base.AbstractSerializer, os::OptimizedSource)
    Base.serialize_type(s, typeof(os))
    write(s.io, os.init_ra)
    write(s.io, os.init_dec)
    for i = 1:length(Celeste.Model.ids)
        write(s.io, os.vs[i])
    end
end

function deserialize(s::Base.AbstractSerializer, ::Type{OptimizedSource})
    init_ra = read(s.io, Float64)::Float64
    init_dec = read(s.io, Float64)::Float64
    vs = zeros(Float64, length(Celeste.Model.ids))
    for i = 1:length(Celeste.Model.ids)
        vs[i] = read(s.io, Float64)
    end
    OptimizedSource(init_ra, init_dec, vs)
end

# ------
# optimization

"""
Optimize the `ts`th element of `target_sources`.
"""
function process_source(config::Config,
                        ts::Int,
                        catalog::Vector{CatalogEntry},
                        target_sources::Vector{Int},
                        neighbor_map::Vector{Vector{Int}},
                        images::Vector{Image})
    s = target_sources[ts]
    entry = catalog[s]
    neighbors = catalog[neighbor_map[ts]]

    tic()
    vs_opt = infer_source(config, images, neighbors, entry)
    Log.info("#$(ts) at ($(entry.pos[1]), $(entry.pos[2])): $(toq()) secs")
    return OptimizedSource(entry.pos[1], entry.pos[2], vs_opt)
end

"""
Use multiple threads to process each target source with the specified
callback and write the results to a file.
"""
function one_node_single_infer(config::Config,
                               catalog::Vector{CatalogEntry},
                               target_sources::Vector{Int},
                               neighbor_map::Vector{Vector{Int}},
                               images::Vector{Image};
                               infer_source_callback=infer_source,
                               timing=InferTiming())
    curr_source = 1
    last_source = length(target_sources)
    sources_lock = SpinLock()
    results = OptimizedSource[]
    results_lock = SpinLock()

    # iterate over sources
    function process_sources()
        tid = threadid()

        while true
            lock(sources_lock)
            ts = curr_source
            curr_source += 1
            unlock(sources_lock)
            if ts > last_source
                break
            end

            try
                result = process_source(config, ts, catalog, target_sources,
                                        neighbor_map, images)

                lock(results_lock)
                push!(results, result)
                unlock(results_lock)
            catch ex
                if is_production_run || nthreads() > 1
                    Log.exception(ex)
                else
                    rethrow(ex)
                end
            end
        end
    end

    tic()
    if nthreads() == 1
        process_sources()
    else
        ccall(:jl_threading_run, Void, (Any,), Core.svec(process_sources))
        ccall(:jl_threading_profile, Void, ())
    end
    timing.opt_srcs = toq()
    timing.num_srcs = length(target_sources)

    return results
end

# legacy wrapper
function one_node_single_infer(catalog::Vector{CatalogEntry},
                               target_sources::Vector{Int},
                               neighbor_map::Vector{Vector{Int}},
                               images::Vector{Image};
                               infer_source_callback=infer_source,
                               timing=InferTiming())
    one_node_single_infer(
        Config(),
        catalog,
        target_sources,
        neighbor_map,
        images,
        infer_source_callback=infer_source_callback,
        timing=timing,
    )
end

"""
Use multiple threads on one node to fit the Celeste model to sources in a given
bounding box.
"""
function one_node_infer(rcfs::Vector{RunCamcolField},
                        stagedir::String;
                        infer_callback=one_node_single_infer,
                        box=BoundingBox(-1000., 1000., -1000., 1000.),
                        primary_initialization=true,
                        timing=InferTiming())
    strategy = PlainFITSStrategy(stagedir)

    catalog, target_sources, neighbor_map, images =
        infer_init(rcfs,
                   strategy;
                   box=box,
                   primary_initialization=primary_initialization)

    Log.info("running with $(nthreads()) threads")

    # NB: All I/O happens above. The methods below don't touch disk.
    infer_callback(catalog, target_sources, neighbor_map, images)
end

"""
Query the SDSS database for all fields that overlap the given RA, Dec range.
"""
function get_overlapping_field_extents(query::BoundingBox, strategy::SDSSIO.IOStrategy)
    f = SDSSIO.readFITS(strategy, SDSSIO.FieldExtents())

    hdu = f[2]::FITSIO.TableHDU

    # read in the entire table.
    all_run = read(hdu, "run")::Vector{Int16}
    all_camcol = read(hdu, "camcol")::Vector{UInt8}
    all_field = read(hdu, "field")::Vector{Int16}
    all_ramin = read(hdu, "ramin")::Vector{Float64}
    all_ramax = read(hdu, "ramax")::Vector{Float64}
    all_decmin = read(hdu, "decmin")::Vector{Float64}
    all_decmax = read(hdu, "decmax")::Vector{Float64}

    close(f)

    ret = Tuple{RunCamcolField, BoundingBox}[]

    # The ramin, ramax, etc is a bit unintuitive because we're looking
    # for any overlap.
    for i in eachindex(all_ramin)
        if (all_ramax[i] > query.ramin && all_ramin[i] < query.ramax &&
                all_decmax[i] > query.decmin && all_decmin[i] < query.decmax)
            cur_box = BoundingBox(all_ramin[i], all_ramax[i],
                                  all_decmin[i], all_decmax[i])
            cur_fe = RunCamcolField(all_run[i], all_camcol[i], all_field[i])
            push!(ret, (cur_fe, cur_box))
        end
    end

    return ret
end
get_overlapping_field_extents(query::BoundingBox, stagedir::String) =
    get_overlapping_field_extents(query, SDSSIO.PlainFITSStrategy(stagedir))

"""
Like `get_overlapping_field_extents()`, but return a Vector of
(run, camcol, field) triplets.
"""
function get_overlapping_fields(query::BoundingBox, stagedir)
    fes = get_overlapping_field_extents(query, stagedir)
    [fe[1] for fe in fes]
end

"""
Store the contents of the `field_extents.fits` file, so it doesn't
have to be loaded repeatedly.
"""
struct FieldExtent
    run::Int16
    camcol::UInt8
    field::Int16
    ramin::Float64
    ramax::Float64
    decmin::Float64
    decmax::Float64
end

"""
Load `field_extents.fits` from `stagedir` and return a vector of
`FieldExtent`s.
"""
function load_field_extents(strategy)
    f = SDSSIO.readFITS(strategy, SDSSIO.FieldExtents())

    hdu = f[2]::FITSIO.TableHDU

    # read in the entire table.
    all_run = read(hdu, "run")::Vector{Int16}
    all_camcol = read(hdu, "camcol")::Vector{UInt8}
    all_field = read(hdu, "field")::Vector{Int16}
    all_ramin = read(hdu, "ramin")::Vector{Float64}
    all_ramax = read(hdu, "ramax")::Vector{Float64}
    all_decmin = read(hdu, "decmin")::Vector{Float64}
    all_decmax = read(hdu, "decmax")::Vector{Float64}

    close(f)

    fes = Vector{FieldExtent}()
    for i = 1:length(all_run)
        fe = FieldExtent(all_run[i], all_camcol[i], all_field[i],
                         all_ramin[i], all_ramax[i],
                         all_decmin[i], all_decmax[i])
        push!(fes, fe)
    end

    return fes
end

"""
Use the provided `FieldExtent`s to return a list of RCFs that overlap
the specified box.
"""
function get_overlapping_fields(query::BoundingBox, fes::Vector{FieldExtent})
    # The ramin, ramax, etc is a bit unintuitive because we're looking
    # for any overlap.
    rcfs = Vector{RunCamcolField}()
    for i in eachindex(fes)
        if (fes[i].ramax > query.ramin &&
                fes[i].ramin < query.ramax &&
                fes[i].decmax > query.decmin &&
                fes[i].decmin < query.decmax)
            push!(rcfs, RunCamcolField(fes[i].run, fes[i].camcol, fes[i].field))
        end
    end

    return rcfs
end

"""
Use the provided `FieldExtent`s to return a list of RCFs that overlap
all the provided boxes.
"""
function get_overlapping_fields(queries::Vector{Vector{BoundingBox}},
                                fes::Vector{FieldExtent})
    rcfs = Vector{RunCamcolField}()

    for boxlist in queries
        for query in boxlist
            # The ramin, ramax, etc is a bit unintuitive because we're looking
            # for any overlap.
            for i in eachindex(fes)
                if (fes[i].ramax > query.ramin &&
                        fes[i].ramin < query.ramax &&
                        fes[i].decmax > query.decmin &&
                        fes[i].decmin < query.decmax)
                    push!(rcfs, RunCamcolField(fes[i].run,
                                               fes[i].camcol,
                                               fes[i].field))
                end
            end
        end
        rcfs = unique(rcfs)
    end

    return rcfs
end

"""
Save provided results to a JLD file.
"""
function save_results(outdir, ramin, ramax, decmin, decmax, results)
    fname = @sprintf("%s/celeste-%.4f-%.4f-%.4f-%.4f-rank%d.jld",
                     outdir, ramin, ramax, decmin, decmax, grank())
    JLD.save(fname, "results", results)
end

save_results(outdir, box, results) =
    save_results(outdir, box.ramin, box.ramax, box.decmin, box.decmax, results)

"""
called from main entry point.
"""
function infer_box(strategy, box::BoundingBox, outdir::String)
    timing = InferTiming()

    Log.info("processing box $(box.ramin), $(box.ramax), $(box.decmin), ",
             "$(box.decmax) with $(nthreads()) threads")

    @time begin
        tic()
        # Get vector of (run, camcol, field) triplets overlapping this patch
        rcfs = get_overlapping_fields(box, strategy)
        timing.query_fids = toq()

        catalog, target_sources, neighbor_map, images =
            infer_init(rcfs,
                       strategy;
                       box=box,
                       primary_initialization=true,
                       timing=timing)

        #results = one_node_single_infer(catalog, target_sources,
        #                                neighbor_map, images; timing=timing)
        results = one_node_joint_infer(catalog, target_sources, neighbor_map,
                                       images; timing=timing)

        tic()
        save_results(outdir, box, results)
        timing.write_results = toq()
    end
    puts_timing(timing)
end

function infer_box(box::BoundingBox, stagedir::String, outdir::String)
    strategy = PlainFITSStrategy(stagedir)
    infer_box(strategy, box, outdir)
end

function do_infer_boxes(pstrategy::ThreadsStrategy,
                        all_rcfs::Vector{RunCamcolField},
                        all_rcf_nsrcs::Vector{Int16},
                        all_boxes::Vector{Vector{BoundingBox}},
                        all_boxes_rcf_idxs::Vector{Vector{Vector{Int32}}},
                        iostrategy::SDSSIO.IOStrategy,
                        prefetch::Bool,
                        outdir::String)
    for boxes in all_boxes
        for box in boxes
            infer_box(iostrategy, box, outdir)
        end
    end
end

"""
called from main entry point.
"""
function infer_boxes(pstrategy::ParallelismStrategy,
                     all_rcfs::Vector{RunCamcolField},
                     all_rcf_nsrcs::Vector{Int16},
                     all_boxes::Vector{Vector{BoundingBox}},
                     all_boxes_rcf_idxs::Vector{Vector{Vector{Int32}}},
                     iostrategy::SDSSIO.IOStrategy,
                     prefetch::Bool,
                     outdir::String)

    # Base.@time hack for distributed environment
    gc_stats = ()
    gc_diff_stats = ()
    gc_stats = Base.gc_num()
    elapsed_time = time_ns()

    do_infer_boxes(pstrategy, all_rcfs, all_rcf_nsrcs, all_boxes,
                   all_boxes_rcf_idxs, iostrategy, prefetch, outdir)

    # Base.@time hack for distributed environment
    elapsed_time = time_ns() - elapsed_time
    gc_diff_stats = Base.GC_Diff(Base.gc_num(), gc_stats)
    time_puts(elapsed_time, gc_diff_stats.allocd, gc_diff_stats.total_time,
              Base.gc_alloc_count(gc_diff_stats))
end

end


text:
module MyBigModule
include("files/big-module-part-1.jl")
include("files/big-module-part-2.jl")

export func2big

end #end of module

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples


text:
"""
Construct reactive power DeviceRangeConstraintInputs for specific types.
"""
function make_reactive_power_constraints_inputs(
    ::Type{T},
    ::Type{U},
    ::Type{V},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.Device, U <: AbstractDeviceFormulation, V <: PM.AbstractPowerModel}
    error(
        "make_reactive_power_constraints_inputs is not implemented for types $T / $U / $V",
    )
end

"""
Default implementation to add active_power constraints.

Users of this function must implement a method for
[`make_reactive_power_constraints_inputs`](@ref) for their specific types.
Users may also implement custom reactive_power_constraints! methods.
"""
function reactive_power_constraints!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, U},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.Device, U <: AbstractDeviceFormulation}
    use_parameters = model_has_parameters(optimization_container)
    use_forecasts = model_uses_forecasts(optimization_container)
    @assert !(use_parameters && !use_forecasts)
    inputs = make_reactive_power_constraints_inputs(
        T,
        U,
        PM.AbstractPowerModel,
        feedforward,
        use_parameters,
        use_forecasts,
    )
    device_range_constraints!(optimization_container, devices, model, feedforward, inputs)
end


text:
module TestTransducers
using Test

@testset "$file" for file in [
        "test_library.jl",
        "test_processes.jl",
        "test_ir.jl",
        ]
    include(file)
end

end  # module


text:
"""
Oja (non-Distributed) streaming PCA method

Author: Wayne Wang
Email: wayneyw@umich.edu
Last modified: 08/26/2021
"""

struct ADAGRAD end
struct ROBBINS_MONRO end

function oja(Xt::AbstractArray{<:Real}, V::AbstractArray{<:Real},
    α::AbstractArray{<:Real}; scheduling::AbstractString = "adagrad",
    stepsize_dim::Int = 1)
    """
    Oja's method
    """
    # initial setting
    if scheduling == "adagrad"
        scheduling = ADAGRAD()
    elseif scheduling == "robbins-monro"
        scheduling = ROBBINS_MONRO()
    end

    # update eigenvectors & learning rates
    V, α = ojaupdate(scheduling, Xt, V, α, stepsize_dim)
    # retraction
    V .= Array(qr!(V).Q)

    return V, α
end

function ojaupdate(scheduling::ROBBINS_MONRO, Xt::AbstractArray{<:Real}, V::AbstractArray{<:Real},
    α::AbstractArray{<:Real}, stepsize_dim::Int)
    """
    Oja × Robbins-Monro
    """
    V .+= α*Xt*(copy(Xt')*V)
    α = nothing
    return V, α
end

function ojaupdate(scheduling::ADAGRAD, Xt::AbstractArray{<:Real}, V::AbstractArray{<:Real},
    α::AbstractArray{<:Real}, stepsize_dim::Int)
    """
    Oja × Adagrad
    """
    grad = Xt*(copy(Xt')*V)
    if stepsize_dim == 1 # α ∈ k × 1
        α .= α .+ sum(abs2,grad,dims=1)[:]
        V .+= grad ./ reshape(sqrt.(α),(1,size(α,1)))
    elseif stepsize_dim == 2 # α ∈ d × k
        α .= α .+ grad .* grad
        V .+= grad ./ sqrt.(α)
    end
    return V, α
end

max_stars_repo_name: ayanbiswas/PRISM

max_issues_repo_name: ayanbiswas/PRISM

max_forks_repo_name: ayanbiswas/PRISM


text:
function writeoutput(s::AbstractSimulation)
    init = s.iteration[]
    #mycopy!(s.rhs,s.u)
    #setfourier!(s.rhs)
    #real!(s.rhs)
    write("u1.$init",s.u.rr.x)
    write("u2.$init",s.u.rr.y)
    write("u3.$init",s.u.rr.z)
    #dealias!(s.rhs)
    if haspassivescalar(s)
        #mycopy!(s.passivescalar.rhs.field.data, s.passivescalar.φ.field.data)
        #setfourier!(s.passivescalar.rhs)
        #real!(s.passivescalar.rhs)
        write("scalar.$init",s.passivescalar.φ.rr)
        #dealias!(s.passivescalar.ρrhs)
    end
    if hasdensity(s)
        #mycopy!(s.densitystratification.rhs.field.data,s.densitystratification.ρ.field.data)
        #setfourier!(s.densitystratification.rhs)
        #real!(s.densitystratification.rhs)
        write("rho.$init",s.densitystratification.ρ.rr)
        if hasdensityles(s)
            write("dflux1.$init",s.densitystratification.flux.rr.x)
            write("dflux2.$init",s.densitystratification.flux.rr.y)
            write("dflux3.$init",s.densitystratification.flux.rr.z)
        end
        #dealias!(s.densitystratification.ρrhs)
    end
    if hasles(s)
        t = s.lesmodel.tau
        isrealspace(t) || real!(t)
        write("T11.$init",t.rr.xx)
        write("T12.$init",t.rr.xy)
        write("T13.$init",t.rr.xz)
        write("T22.$init",t.rr.yy)
        write("T23.$init",t.rr.yz)
        if is_dynamic_les(s)
            write("cs.$init",s.lesmodel.c.rr)
            if is_dynP_les(s)
                write("cp.$init",s.lesmodel.cp.rr)
            end
        end
        #is_FakeSmagorinsky(s) || write("pr.$init",s.lesmodel.pr)
    end
    if typeof(s.forcing) <: RfForcing
        writedlm("R.$init",s.forcing.R)
        writedlm("Eh.$init",zip(s.forcing.avgK, s.forcing.Ef))
    end
    if typeof(s.forcing) <: AForcing
        writedlm("R.$init",s.forcing.R)
    end
    if typeof(s.forcing) <: NRfForcing
        writedlm("R.$init",s.forcing.R)
        writedlm("N.$init",s.forcing.N)
        writedlm("Eh.$init",zip(s.forcing.avgK, s.forcing.Ef))
    end

    if typeof(s.forcing) <: CNRfForcing
        writedlm("R.$init",s.forcing.R)
        writedlm("N.$init",s.forcing.N)
        writedlm("Eh.$init",zip(s.forcing.avgK, s.forcing.Ef))
    end

    if typeof(s.forcing) <: CRfForcing
        writedlm("R.$init",s.forcing.R)
        writedlm("Eh.$init",zip(s.forcing.avgK, s.forcing.Ef))
    end
end


text:
type Tokenizer
    dict::IdDict
    model
end

function Tokenizer()
    dict = IdDict(String["UNKNOWN", " ", "LF"])
    T = Float32
    embed = Embedding(T, 100, 10)
    ls = [Linear(T,70,70), Linear(T,70,4)]
    #=
    g = @graph begin
        chars = identity(:chars)
        x = Var(reshape(chars,1,length(chars)))
        x = embed(x)
        x = window2d(x,10,7,1,1,0,3)
        x = ls[1](x)
        x = relu(x)
        x = ls[2](x)
        x
    end
    =#
    Tokenizer(dict, g)
end

function (t::Tokenizer)(chars::Vector{Char})
    unk = t.dict["UNKNOWN"]
    x = map(chars) do c
        str = c == '\n' ? "LF" : string(c)
        get(t.dict, str, unk)
    end
    y = t.model(x).data
    tags = argmax(y,1)
    decode(chars, tags)
end
(t::Tokenizer)(str::String) = t(Vector{Char}(str))

function decode(chars::Vector{Char}, tags::Vector{Int})
    length(chars) == length(tags) || throw("Length unmatch")
    doc = Vector{String}[]
    tokens = String[]
    buffer = Char[]
    pos = 1
    function push_buffer!(i::Int)
        push!(buffer, chars[i])
    end
    function push_token!()
        length(buffer) == 0 && return
        form = String(buffer)
        #t = Token(pos=pos, form=form)
        push!(tokens, form)
        pos += length(buffer) + 1
        buffer = Char[]
    end
    function push_tokens!()
        length(tokens) == 0 && return
        push!(doc, tokens)
        tokens = String[]
    end

    for i = 1:length(chars)
        t = tags[i]
        if t == 1
            push_buffer!(i)
        elseif t == 2
            continue
        elseif t == 3
            push_buffer!(i)
            push_token!()
        elseif t == 4
            push_buffer!(i)
            push_token!()
            push_tokens!()
        else
            error("Invalid tag: $(t)")
        end
    end
    push_token!()
    push_tokens!()
    doc
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
using Documenter, Gadfly, Compose, Cairo

makedocs(
    modules = [Gadfly],
    clean = false,
    sitename = "Gadfly.jl",
    pages = Any[
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
        "Manual" => Any[
            "Plotting" => "man/plotting.md",
            "Compositing" => "man/compositing.md",
            "Backends" => "man/backends.md",
            "Themes" => "man/themes.md",
        ],
        "Gallery" => Any[
            "Geometries" => "gallery/geometries.md",
            "Guides" => "gallery/guides.md",
            "Statistics" => "gallery/statistics.md",
            "Coordinates" => "gallery/coordinates.md",
            "Scales" => "gallery/scales.md",
            "Shapes" => "gallery/shapes.md",
        ],
        "Library" => Any[
            "Gadfly" => "lib/gadfly.md",
            "Geometries" => "lib/geometries.md",
            "Guides" => "lib/guides.md",
            "Statistics" => "lib/statistics.md",
            "Coordinates" => "lib/coordinates.md",
            "Scales" => "lib/scales.md",
            "Shapes" => "lib/shapes.md",
        ],
        "Development" => Any[
            "Rendering Pipeline" => "dev/pipeline.md",
            "Regression Testing" => "dev/regression.md",
            "Relationship with Compose.jl" => "dev/compose.md",
        ]
    ]
)

deploydocs(
    repo   = "github.com/GiovineItalia/Gadfly.jl.git",
)


text:
module TestSort

using FactCheck
using DataCubes

facts("Sort tests") do
  @fact sort(@darr(a=[5 5 3;3 3 1],b=[1 3 4;9 8 7]), 1, :a) --> @darr(a=[3 3 1;5 5 3],b=[9 8 7;1 3 4])
  @fact sort(@darr(a=1.0*[5 5 3;3 3 1],b=[1 3 4;9 8 7]), 1, :a) --> @darr(a=1.0*[3 3 1;5 5 3],b=[9 8 7;1 3 4])
  @fact sort(@darr(a=[5 5 3;3 3 1],b=[1 3 4;9 8 7]), 2, :a) --> @darr(a=[3 5 5;1 3 3],b=[4 1 3;7 9 8])
  @fact sort(@darr(a=1.0*[5 5 3;3 3 1],b=[1 3 4;9 8 7]), 2, :a) --> @darr(a=1.0*[3 5 5;1 3 3],b=[4 1 3;7 9 8])
  @fact sort(@darr(a=[5 5 3;3 3 1],b=[3 1 4;8 9 7]), 2, :a,:b) --> @darr(a=[3 5 5;1 3 3],b=[4 1 3;7 9 8])
  @fact sort(@darr(a=[5 5 3;3 3 1],b=1.0*[3 1 4;8 9 7]), 2, :a,:b) --> @darr(a=[3 5 5;1 3 3],b=1.0*[4 1 3;7 9 8])
  @fact sort(@larr(a=[5 5 3;3 3 1],b=[1 3 4;9 8 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]]), 1, :a) --> @larr(a=[3 3 1;5 5 3],b=[9 8 7;1 3 4],axis1[k=[:b,:a]],axis2[r=[:X,:Y,:Z]])
  @fact sort(@larr(a=[5 5 3;3 3 1],b=[1 3 4;9 8 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]]), 2, :a) --> @larr(a=[3 5 5;1 3 3],b=[4 1 3;7 9 8],axis1[k=[:a,:b]],axis2[r=[:Z,:X,:Y]])
  @fact sort(@larr(a=1.0*[5 5 3;3 3 1],b=[1 3 4;9 8 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]]), 2, :a) --> @larr(a=1.0*[3 5 5;1 3 3],b=[4 1 3;7 9 8],axis1[k=[:a,:b]],axis2[r=[:Z,:X,:Y]])
  @fact sort(@larr(a=[5 5 3;3 3 1],b=[3 1 4;8 9 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]]), 2, :a,:b) --> @larr(a=[3 5 5;1 3 3],b=[4 1 3;7 9 8],axis1[k=[:a,:b]],axis2[r=[:Z,:Y,:X]])
  @fact sort(@larr(a=[5 5 3;3 3 1],b=[3 1 4;8 9 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]]), 2, :a,:b, a_rev=true) --> @larr(a=[5 5 3;3 3 1],b=[1 3 4;9 8 7],axis1[k=[:a,:b]],axis2[r=[:Y,:X,:Z]])
  @fact sort(@larr(a=[5 5 3;3 3 1],b=[3 1 4;8 9 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]]), 2, :a,:b, a_rev=true, b_rev=true) --> @larr(a=[5 5 3;3 3 1],b=[3 1 4;8 9 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]])
  @fact sort(@larr(a=[5 5 3;3 3 1],b=[3 1 4;8 9 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]]), 2, :b, b_coords=(1,:)) --> larr(a=[5 5 3;3 3 1],b=[1 3 4;9 8 7],axis1=darr(k=[:a,:b]),axis2=darr(r=[:Y,:X,:Z]))
  @fact sort(@larr(a=[5 5 3;3 3 1],b=[3 1 4;8 9 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]]), 2, :b, b_coords=(2,:)) --> larr(a=[3 5 5;1 3 3],b=[4 3 1;7 8 9],axis1=darr(k=[:a,:b]),axis2=darr(r=[:Z,:X,:Y]))

  @fact sort(@darr(a=[5.0 5.0 3.0;3.0 3.0 1.0],b=[1.0 3.0 4.0;9.0 8.0 7.0]), 1, :a) --> @darr(a=1.0*[3 3 1;5 5 3],b=1.0*[9 8 7;1 3 4])
  @fact sort(@darr(a=1.0*[5 5 3;3 3 1],b=1.0*[1 3 4;9 8 7]), 2, :a) --> @darr(a=1.0*[3 5 5;1 3 3],b=1.0*[4 1 3;7 9 8])
  @fact sort(@darr(a=1.0*[5 5 3;3 3 1],b=1.0*[3 1 4;8 9 7]), 2, :a,:b) --> @darr(a=1.0*[3 5 5;1 3 3],b=1.0*[4 1 3;7 9 8])
  @fact sort(@larr(a=1.0*[5 5 3;3 3 1],b=1.0*[1 3 4;9 8 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]]), 1, :a) --> @larr(a=1.0*[3 3 1;5 5 3],b=1.0*[9 8 7;1 3 4],axis1[k=[:b,:a]],axis2[r=[:X,:Y,:Z]])
  @fact sort(@larr(a=[5 5 3;3 3 1],b=1.0*[1 3 4;9 8 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]]), 2, :a) --> @larr(a=[3 5 5;1 3 3],b=1.0*[4 1 3;7 9 8],axis1[k=[:a,:b]],axis2[r=[:Z,:X,:Y]])
  @fact sort(@larr(a=1.0*[5 5 3;3 3 1],b=[3 1 4;8 9 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]]), 2, :a,:b) --> @larr(a=1.0*[3 5 5;1 3 3],b=[4 1 3;7 9 8],axis1[k=[:a,:b]],axis2[r=[:Z,:Y,:X]])
  @fact sort(@larr(a=1.0*[5 5 3;3 3 1],b=1.0*[3 1 4;8 9 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]]), 2, :a,:b, a_rev=true) --> @larr(a=1.0*[5 5 3;3 3 1],b=1.0*[1 3 4;9 8 7],axis1[k=[:a,:b]],axis2[r=[:Y,:X,:Z]])
  @fact sort(@larr(a=[5 5 3;3 3 1],b=1.0*[3 1 4;8 9 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]]), 2, :a,:b, a_rev=true, b_rev=true) --> @larr(a=[5 5 3;3 3 1],b=1.0*[3 1 4;8 9 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]])
  @fact sort(@larr(a=1.0*[5 5 3;3 3 1],b=1.0*[3 1 4;8 9 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]]), 2, :b, b_coords=(1,:)) --> larr(a=1.0*[5 5 3;3 3 1],b=1.0*[1 3 4;9 8 7],axis1=darr(k=[:a,:b]),axis2=darr(r=[:Y,:X,:Z]))
  @fact sort(@larr(a=1.0*[5 5 3;3 3 1],b=[3 1 4;8 9 7],axis1[k=[:a,:b]],axis2[r=[:X,:Y,:Z]]), 2, :b, b_coords=(2,:)) --> larr(a=1.0*[3 5 5;1 3 3],b=[4 3 1;7 8 9],axis1=darr(k=[:a,:b]),axis2=darr(r=[:Z,:X,:Y]))
  @fact sort(@larr(a=[5 NA 3;3 3 1],b=[3 1 4;8 9 7],axis1[k = [:a,:b]],axis2[r = [:X,:Y,:Z]]),2,:a,:b, a_rev=true) --> @larr(a=[5 3 NA;3 1 3],b=[3 4 1;8 7 9],axis1[k = [:a,:b]],axis2[r = [:X,:Z,:Y]])
  @fact sort(@larr(a=[5 NA 3;3 3 1],b=[3 1 4;8 9 7],axis1[k = [:a,:b]],axis2[r = [:X,:Y,:Z]]),2,:a,:b, a_rev=true,b_rev=true) --> @larr(a=[5 3 NA;3 1 3],b=[3 4 1;8 7 9],axis1[k = [:a,:b]],axis2[r = [:X,:Z,:Y]])
  @fact sort(@larr(a=[5 NA 3;3 3 1],b=[3 1 4;8 9 7],axis1[k = [:a,:b]],axis2[r = [:X,:Y,:Z]]),2,:a,:b) --> @larr(a=[NA 3 5;3 1 3],b=[1 4 3;9 7 8],axis1[k = [:a,:b]],axis2[r = [:Y,:Z,:X]])
  @fact sort(@larr(a=[5 NA 3;3 3 1],b=[3 1 4;8 9 7],axis1[k = [:a,:b]],axis2[r = [:X,:Y,:Z]]),2,:a,:b, a_rev=true,b_rev=true) --> @larr(a=[5 3 NA;3 1 3],b=[3 4 1;8 7 9],axis1[k = [:a,:b]],axis2[r = [:X,:Z,:Y]])
  @fact sort(@larr(a=[5 NA 3;3 3 1],b=[3 1 4;8 9 7],axis1[k = [:a,:b]],axis2[r = [:X,:Y,:Z]]),2,:a,:b) --> @larr(a=[NA 3 5;3 1 3],b=[1 4 3;9 7 8],axis1[k = [:a,:b]],axis2[r = [:Y,:Z,:X]])
  @fact sort(@larr(a=@enumeration([5 NA 3;3 3 1],[1,2,3,4,5]),b=[3 1 4;8 9 7],axis1[k = [:a,:b]],axis2[r = [:X,:Y,:Z]]),2,:a,:b) --> @larr(a=[NA 3 5;3 1 3],b=[1 4 3;9 7 8],axis1[k = [:a,:b]],axis2[r = [:Y,:Z,:X]])
  @fact sort(@larr(a=@enumeration([5 NA 3;3 3 1],[1,2,3,4,5]),b=enumeration([3 1 4;8 9 7],[10,9,8,7,6,5,4,3,2,1]),axis1[k = [:a,:b]],axis2[r = [:X,:Y,:Z]]),2,:a,:b,b_rev=true) --> @larr(a=[NA 3 5;3 1 3],b=[1 4 3;9 7 8],axis1[k = [:a,:b]],axis2[r = [:Y,:Z,:X]])
  @fact sort(@larr(a=[5.0 NA 3.0;3.0 3.0 1.0],b=[3.0 1.0 4.0;8.0 9.0 7.0],axis1[k = [:a,:b]],axis2[r = [:X,:Y,:Z]]),2,:a,:b, a_rev=true,b_rev=true) --> @larr(a=[5.0 3.0 NA;3.0 1.0 3.0],b=[3.0 4.0 1.0;8.0 7.0 9.0],axis1[k = [:a,:b]],axis2[r = [:X,:Z,:Y]])
  @fact sort(@larr(a=[5.0 NA 3.0;3.0 3.0 1.0],b=[3.0 1.0 4.0;8.0 9.0 7.0],axis1[k = [:a,:b]],axis2[r = [:X,:Y,:Z]]),2,:a,:b) --> @larr(a=[NA 3.0 5.0;3.0 1.0 3.0],b=[1.0 4.0 3.0;9.0 7.0 8.0],axis1[k = [:a,:b]],axis2[r = [:Y,:Z,:X]])
  @fact sort(larr(a=enumeration([1,7,4,3,5,3]),axis1=enumeration([1,12,13,14,15,16])),1,:a) --> larr(a=[1,7,4,3,3,5], axis1=[1,12,13,14,16,15])
  @fact sortperm(larr(a=enumeration([1,7,4,3,5,3]),axis1=enumeration([1,12,13,14,15,16])),1,:a) --> ([1,2,3,4,6,5],)
  @fact sortperm(larr(a=[1,7,4,3,5,3],axis1=enumeration([1,12,13,14,15,16])),1,:a) --> ([1,4,6,3,5,2],)
  for i in 1:100
    @fact (t=darr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5));sort(t)==sort(t,1)) --> true
    @fact (t=darr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5));sort(t)==sort(t,1)) --> true
    @fact (t=darr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5));sort(t)==sort(t,1,:a,:b,:c)) --> true
    @fact (t=darr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5));sort(t,2)==sort(t,2,:a,:b,:c)) --> true
    @fact (t=darr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5));sort(t,2)==sort(t,2,:a,:b,:c)) --> true

    @fact (t=larr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5),axis=rand(0:1,3));sort(t)==sort(t,1)) --> true
    @fact (t=larr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5),axis=darr(k=rand(0:1,3)));sort(t)==sort(t,1)) --> true
    @fact (t=larr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5),axis=darr(k=rand(0:1,3)));sort(t)==sort(t,1,:k)) --> true
    @fact (t=larr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5),axis2=darr(k=rand(0:1,5)));sort(t,2)==sort(t,2,:k)) --> true
    @fact (t=larr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5),axis2=darr(k1=rand(0:1,5),k2=rand(0:1,5)));sort(t,2)==sort(t,2,:k1,:k2)) --> true

    @fact (t=darr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5));sortperm(t)==sortperm(t,1)) --> true
    @fact (t=darr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5));sortperm(t)==sortperm(t,1)) --> true
    @fact (t=darr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5));sortperm(t)==sortperm(t,1,:a,:b,:c)) --> true
    @fact (t=darr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5));sortperm(t,2)==sortperm(t,2,:a,:b,:c)) --> true
    @fact (t=darr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5));sortperm(t,2)==sortperm(t,2,:a,:b,:c)) --> true

    @fact (t=larr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5),axis=rand(0:1,3));sortperm(t)==sortperm(t,1)) --> true
    @fact (t=larr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5),axis=darr(k=rand(0:1,3)));sortperm(t)==sortperm(t,1)) --> true
    @fact (t=larr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5),axis=darr(k=rand(0:1,3)));sortperm(t)==sortperm(t,1,:k)) --> true
    @fact (t=larr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5),axis2=darr(k=rand(0:1,5)));sortperm(t,2)==sortperm(t,2,:k)) --> true
    @fact (t=larr(a=rand(0:1,3,5),b=rand(0:1,3,5),c=rand(0:1,3,5),axis2=darr(k1=rand(0:1,5),k2=rand(0:1,5)));sortperm(t,2)==sortperm(t,2,:k1,:k2)) --> true
  end
end

end


text:
function orthogonal_sum(G1::LocalGenusHerm, G2::LocalGenusHerm)
  @req prime(G1) === prime(G2) "Local genera must have the same prime ideal"
  if !G1.isdyadic || !G2.isramified
    return _direct_sum_easy(G1, G2)
  else
    L1 = representative(G1)
    L2 = representative(G2)
    L3, = orthogonal_sum(L1, L2)
    return genus(L3, prime(G1))
  end
end

function _direct_sum_easy(G1::LocalGenusHerm, G2::LocalGenusHerm)
  # We do a merge sort
  i1 = 1
  i2 = 1
  z = Tuple{Int, Int, Int}[]
  while i1 <= length(G1) && i2 <= length(G2)
    if scale(G1, i1) < scale(G2, i2)
      push!(z, G1.data[i1])
      i1 += 1
    elseif scale(G2, i2) < scale(G1, i1)
      push!(z, G2.data[i2])
      i2 += 1
    else
      @assert scale(G1, i1) == scale(G2, i2)
      push!(z, (scale(G1, i1), rank(G1, i1) + rank(G2, i2), det(G1, i1) * det(G2, i2)))
      i1 += 1
      i2 += 1
    end
  end

  if i1 <= length(G1)
    append!(z, G1.data[i1:length(G1)])
  end

  if i2 <= length(G2)
    append!(z, G2.data[i2:length(G2)])
  end

  return genus(HermLat, base_field(G1), prime(G1), z)
end

# Genus representatives for quadratic lattices,
#
# With permission ported from Magma package of Markus Kirschmer:
# http://www.math.rwth-aachen.de/~Markus.Kirschmer/magma/lat.html

################################################################################
#
#  Genus representatives for hermitian lattices
#
################################################################################

@doc Markdown.doc"""
    genus_representatives(L::HermLat; max = inf, use_auto = false)
                                                        -> Vector{HermLat}

Computes representatives for the isometry classes in the genus of $L$.

At most `max` representatives are returned. The use of automorphims can be
disabled by `use_auto = false`. The computation of the mass can be enable
by `use_mass = true`.
"""
function genus_representatives(L::HermLat; max = inf, use_auto::Bool = true,
                                                      use_mass::Bool = false)
  @req rank(L) >= 2 "Lattice must have rank >= 2"
  R = base_ring(L)
  definite = isdefinite(L)
  gens, fl, P0 = genus_generators(L)
  a = involution(L)
  LL = typeof(L)[ L ]
  for g in gens
    if definite && g[1] == P0
      continue
    end
    I = g[1]^Int(g[2] - 1)
    J = a(I)
    N = neighbours_with_ppower(L, g[1], g[2] - 1, use_auto = use_auto)::Vector{typeof(L)}
    inter = typeof(L)[]
    for i in 2:length(LL)
      M = pseudo_matrix(LL[i])
      IM = I * M
      JM = J * M
      inter = append!(inter, lattice(ambient_space(L), meet(IM + pseudo_matrix(x), JM) for x in N))
    end
    LL = vcat(LL, N, inter)
  end
  @assert length(LL) == prod(Int[g[2] for g in gens if !definite || g[1] != P0])
  @assert all(X -> genus(X) == genus(L), LL)

  if use_mass
    mass = Hecke.mass(L)
  else
    mass = zero(fmpq)
  end

  missing_mass = Ref(mass)

  local result::Vector{typeof(L)}

  if definite
    result = typeof(L)[]
    for L in LL
      # Should never happen!
      @assert all(X -> !isisometric(X, L), result)
      neig = iterated_neighbours(L, P0, max = max, use_auto = use_auto, missing_mass = missing_mass)
      append!(result, neig)
      max = max - length(result)
    end
    for i in 1:length(result)
      for j in 1:i-1
        @hassert :GenRep 1 !(isisometric(result[i], result[j])[1])
      end
    end
  else
    result = LL
  end
  return result
end

function _neighbours(L, P, result, max, callback = stdcallback, use_auto = true)
  ok, scale = ismodular(L, P);
  @req ok "The lattice must be locally modular"
  R = base_ring(L)
  K = nf(R)
  a = involution(L)
  T = local_basis_matrix(L, P, type = :submodule)
  # This is a bit tricky, since C does not know enough to construct the residue field ...
  # C = a(P)
  p = minimum(P)
  lp = prime_decomposition(R, p)
  if length(lp) == 1
    C = P
  else
    if lp[1][1] == P
      C = lp[2][1]
    elseif lp[2][1] == P
      C = lp[1][1]
    else
      error("This should not happen.")
    end
  end
  k, h = ResidueField(R, C)
  hext = extend(h, K)
  local form::dense_matrix_type(K)
  form = gram_matrix(ambient_space(L))
  special = false
  if scale != 0
    if isramified(R, minimum(P))
      special = isodd(scale)
      scale = div(scale + 1, 2)
    end
    form = K(elem_in_nf(uniformizer(minimum(P))))^(-scale) * form
  end
  n = rank(L)
  W = vector_space(k, n)

  if use_auto
    G = automorphism_group_generators(L)
    @hassert :GenRep 1 all(g -> g * gram_matrix(ambient_space(L)) * transpose(g) == gram_matrix(ambient_space(L)), G)
    Tinv = inv(T)
    adjust_gens = eltype(G)[T * g * Tinv for g in G]
    @hassert :GenRep 1 all(g -> g * form * transpose(g) == form, adjust_gens)
    adjust_gens_mod_p = dense_matrix_type(k)[map_entries(hext, g) for g in adjust_gens]
    adjust_gens_mod_p = dense_matrix_type(k)[x for x in adjust_gens_mod_p if !isdiagonal(x)]
    @hassert :GenRep 1 all(g -> g * pform * transpose(g) == pform, adjust_gens_mod_p)
    if length(adjust_gens_mod_p) > 0
      _LO = line_orbits(adjust_gens_mod_p)
      LO = Vector{eltype(k)}[x[1] for x in _LO]
      @vprint :GenRep 1 "Checking $(length(LO)) representatives (instead of $(div(order(k)^n - 1, order(k) - 1)))\n"
    else
      @vprint :GenRep 1 "Enumerating lines over $k of length $n\n"
      LO = enumerate_lines(k, n)
    end
  else
    @vprint :GenRep 1 "Enumerating lines over $k of length $n\n"
    LO = enumerate_lines(k, n)
  end

  keep = true
  cont = true

  if P != C
    _G = T * form * _map(transpose(T), a)
    G = map_entries(hext, _G)
    pi = p_uniformizer(P)
    pih = h(pi)
    for w in LO
      x = elem_type(K)[ sum(T[i, j] * (hext\w[i]) for i in 1:n) for j in 1:ncols(T)]
      LL = neighbour(L, T, pih * matrix(k, 1, length(w), w) * G, K(pi) .* x, hext, P, C, true)
      keep, cont = callback(result, LL)
      if keep
        push!(result, LL)
      end
      if !cont || (length(result) >= max)
        break
      end
    end
  elseif special
    pi = uniformizer(P)
    _G = elem_in_nf(pi) * T * form * _map(transpose(T), a)
    G = map_entries(hext, _G)
    for w::Vector{fq} in LO
      Gw = G * matrix(k, length(w), 1, w)
      ok = 0
      for d in 1:n
        if !iszero(Gw[d])
          ok = d
          break
        end
      end
      @assert ok != 0
      x = elem_type(K)[ sum(T[i, j] * (hext\w[i]) for i in 1:n) for j in 1:ncols(T)]
      nrm = _inner_product(form, x, x, a)
      b = hext\(-hext(nrm) // (2*Gw[ok, 1]))
      #x = x + b * pi * B[ok]
      for j in 1:ncols(T)
        x[j] = x[j] + b * elem_in_nf(pi) * T[ok, j]
      end
      nrm = _inner_product(form, x, x, a)
      LL = neighbour(L, T, matrix(k, 1, length(w), w) * G, x, hext, P, P, false)
      keep, cont = callback(result, LL)
      if keep
        push!(result, LL)
      end
      if !cont || length(result) >= max
        break
      end
    end
  else
    _G = T * form * _map(transpose(T), a)
    G = map_entries(hext, _G)
    ram = isramified(R, minimum(P))
    if ram
      pi = uniformizer(P)
      S = [ h\x for x in k ]
    else
      p = minimum(P)
      pi = uniformizer(p)
      kp, hp = ResidueField(order(p), p)
      alpha = h\(degree(k) == 1 ? one(k) : gen(k))
      Tram = matrix(kp, 2, 1, [2, hp(tr(alpha))])
    end
    for w::Vector{fq} in LO
      __w = [ (hext\w[i]) for i in 1:n]
      x = [ sum(T[i, j] * (__w[i]) for i in 1:n if !iszero(w[i])) for j in 1:ncols(T)]
      nrm = _inner_product(form, x, x, a) #(x*Form, v) where v:= Vector([ a(e): e in Eltseq(x) ]);
      if !(nrm in P)
        continue
      end
      wG = matrix(k, 1, length(w), w) * G
      ok = 0
      for j in 1:n
        if !iszero(wG[j])
          ok = j
          break
        end
      end
      @assert ok != 0
      if !ram
        el = order(p)(base_field(K)(nrm)//pi)
        b, s, V = can_solve_with_kernel(Tram, matrix(kp, 1, 1, [hp(-el)]), side = :left)
        @assert b
        @assert s * Tram == matrix(kp, 1, 1, [hp(-el)])
        #s, V = solution(Tram, vector(kp, [hp(-el)]))
        #@show _all_row_span(V)
        _kernel = [ matrix(kp, 1, 2, v) for v in _all_row_span(V)]
        l = a(hext\(inv(wG[ok])))
        S = elem_type(K)[ l * (hext\((s + v)[1]) + (hext\(s + v)[2])*alpha) for v in _kernel ]
      end
      for s in S
        LL = neighbour(L, T, wG, elem_type(K)[x[o] + K(elem_in_nf(pi))*s*T[ok, o] for o in 1:ncols(T)], hext, P, P, false)
        keep, cont = callback(result, LL)
        if keep
          push!(result, LL)
        end
        if !cont || (length(result) >= max)
          break
        end
      end
    end
  end
  return result
end

function neighbour(L, B, xG, x, h, P, CC, split)
  R = order(P)
  K = nf(R)
  n = nrows(B)

  local C::ideal_type(R)

  if CC isa Int
    C = split ? involution(L)(P) : P
  else
    C = CC
  end

  I = Int[ i for i in 1:n if !iszero(xG[i])]
  i = I[1]
  a = involution(L)
  M = zero_matrix(K, n - length(I), ncols(B))
  for (k, nk) in enumerate(setdiff(1:n, I))
    for j in 1:ncols(B)
      M[k, j] = B[nk, j]
    end
  end
  CI = fractional_ideal_type(R)[ K(1) * R for j in 1:(n - length(I)) ]
  for j in I
    if j != i
      M = vcat(M, matrix(K, 1, ncols(B), elem_type(K)[B[j, k] - a(h\(divexact(xG[j], xG[i]))) * B[i, k] for k in 1:ncols(B)]))
      push!(CI, K(1) * R)
    end
  end
  M = vcat(M, sub(B, i:i, 1:ncols(B)))
  push!(CI, fractional_ideal(order(P), P))
  M = vcat(M, matrix(K, 1, ncols(B), x))
  push!(CI, inv(C))
  pm = pseudo_matrix(M, CI)
  _M = _sum_modules(L, _module_scale_ideal((split ? P : P * C), pseudo_matrix(L)), pm)
  LL = lattice(ambient_space(L), _M)

  @assert islocally_isometric(L, LL, P)
  return LL
end

function neighbours(L::HermLat, P, max = inf)
#{The immediate P-neighbours of L}
#  require Order(P) eq BaseRing(L): "Arguments are incompatible";
#  require IsPrime(P): "Second argument must be prime";
#  require not IsRamified(P) or Minimum(P) ne 2: "Second argument cannot be a ramified prime over 2";
#  require IsModular(L, P) : "The lattice must be locally modular";
#  require Rank(L) ge 2: "The rank of the lattice must be at least 2";
#  require IsIsotropic(L, P): "The lattice must be locally isotropic";
#
  return _neighbours(L, P, [], max)
end

function stdcallback(list, L)
  keep = all(LL -> !isisometric(LL, L)[1], list)
#//  keep, #List;
  return keep, true;
end

function iterated_neighbours(L::HermLat, P; use_auto = false, max = inf,
                                            callback = false,
                                            missing_mass = Ref{fmpq}(zero(fmpq)))
  # AutoOrbits, CallBack:= false, Max:= Infinity()) -> []
  #require Order(P) eq BaseRing(L): "Arguments are incompatible";
  #require IsPrime(P): "Second argument must be prime";
  #require not IsRamified(P) or Minimum(P) ne 2: "Second argument cannot be a ramified prime over 2";
  #require IsModular(L, P) : "The lattice must be locally modular";
  #require Rank(L) ge 2: "The rank of the lattice must be at least 2";
  #require IsIsotropic(L, P): "The lattice must be locally isotropic";

  if callback == false && isdefinite(L)
    _callback = stdcallback
  else
    _callback = callback
  end

  result = typeof(L)[ L ]

  use_mass = !iszero(missing_mass[])

  if use_mass
    _mass = missing_mass[] - 1//automorphism_group_order(L)
  end

  i = 1
  oldlength = length(result)
  while length(result) < max && i <= length(result)
    # _Neighbours and the callback only add new lattices if not isometric to known ones and stop if Max is reached.
    # So we have nothing to at all.
    result = _neighbours(result[i], P, result, max, _callback)# : CallBack:= CallBack);
    no_lattices = length(result) - oldlength
    oldlength = length(result)
    if use_mass && no_lattices > 0
      _mass = _mass - sum(fmpq[1//automorphism_group_order(result[i]) for i in (length(result) - no_lattices + 1):length(result)])
      if iszero(_mass)
        break
      end
    end
    if use_mass && _mass < 0
      throw(error("This should not happen"))
    end
    i = i + 1
  end
  if use_mass
    missing_mass[] = _mass
  end
  return result
end


text:
using NeuralQuantum, Test
using NeuralQuantum: set_index!
num_types = [Float32, Float64]

re_machines = Dict()
im_machines = Dict()

ma = (T, N) -> NDM(T, N, 2, 3)
re_machines["NDM"] = ma

ma = (T, N) -> RBMSplit(T, N, 2)
im_machines["RBMSplit"] = ma

all_machines = merge(re_machines, im_machines)

N = 4

@testset "Test Properties $name" for name=keys(all_machines)
    for T=num_types
        if name ∈ keys(im_machines)
            T = Complex{T}
        end
        net = all_machines[name](T,N)
        cnet = cached(net)

        @test NeuralQuantum.input_type(net) == real(T)
        @test NeuralQuantum.out_type(net) == Complex{real(T)}
        @test NeuralQuantum.input_shape(net) == (N, N)
        @test NeuralQuantum.is_analytic(net)
    end
end

@testset "Test Cached Value $name" for name=keys(all_machines)
    for T=num_types
        net = all_machines[name](T,N)
        cnet = cached(net)

        v = state(T, SpinBasis(1//2)^N, net)
        # compute exact
        vals = []; cvals = [];
        for i=1:spacedimension(v)
            set_index!(v, i)
            push!(vals, net(v))
            push!(cvals, cnet(v))
        end
        @test vals ≈ cvals
    end
end

@testset "Test Cached Gradient $name" for name=keys(all_machines)
    for T=num_types
        net = all_machines[name](T,N)
        cnet = cached(net)
        cder = grad_cache(net)

        v = state(T, SpinBasis(1//2)^N, net)
        # compute exact
        grads = [];
        for i=1:spacedimension(v)
            set_index!(v, i)
            der_ad = ∇logψ(net,  v)
            ∇logψ!(cder, cnet, v)
            for f=propertynames(der_ad)
                ∇  = getproperty(der_ad,  f)
                c∇ = getproperty(cder, f)
                push!(grads, ∇ ≈ c∇)
            end
        end
        @test all(grads)
    end
end


text:
"""
    @params struct_def

A macro that changes all fields' types to type parameters while respecting the type bounds specificied by the user. For example:
```
@params struct MyType{T}
    f1::T
    f2::AbstractVector{T}
    f3::AbstractVector{<:Real}
    f4
end
```
will define:
```
struct MyType{T, T1 <: T, T2 <: AbstractVector{T}, T3 <: AbstractVector{<:Real}, T4}
    f1::T1
    f2::T2
    f3::T3
    f4::T4
end
```
The default non-parameteric constructor, e.g. `MyType(f1, f2, f3, f4)`, will always work if all the type parameters used in the type header are used in the field types. Using a type parameter in the type header that is not used in the field types such as:
```
struct MyType{T}
    f1
    f2
end
```
is not recommended.
"""
macro params(struct_expr)
    header = struct_expr.args[2]
    fields = @view struct_expr.args[3].args[2:2:end]
    params = []
    for i in 1:length(fields)
        x = fields[i]
        T = gensym()
        if x isa Symbol
            push!(params, T)
            fields[i] = :($x::$T)
        elseif x.head == :(::)
            abstr = x.args[2]
            var = x.args[1]
            push!(params, :($T <: $abstr))
            fields[i] = :($var::$T)
        end
    end
    if header isa Symbol && length(params) > 0
        struct_expr.args[2] = :($header{$(params...)})
    elseif header.head == :curly
        append!(struct_expr.args[2].args, params)
    elseif header.head == :<:
        if struct_expr.args[2].args[1] isa Symbol
            name = struct_expr.args[2].args[1]
            struct_expr.args[2].args[1] = :($name{$(params...)})
        elseif header.head == :<: && struct_expr.args[2].args[1] isa Expr
            append!(struct_expr.args[2].args[1].args, params)
        else
            error("Unidentified type definition.")
        end
    else
        error("Unidentified type definition.")
    end
    esc(struct_expr)
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Helpers for SGP4 algorithm.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export sgp4

"""
    sgp4(Δt, args...)

Function that initialize the SGP4 structure and propagate the orbit until the time Δt.

# Returns

* The position vector [km].
* The velocity vector [km/s].
* The SGP4 structure (see `SGP4_Structure`).

"""
function sgp4(Δt, args...)
    sgp4d = sgp4_init(args...)
    r,v   = sgp4!(sgp4d, Δt)
    return r,v,sgp4d
end


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

mutable struct ServiceDelegationPropertiesFormat <: SwaggerModel
    serviceName::Any # spec type: Union{ Nothing, String } # spec name: serviceName
    actions::Any # spec type: Union{ Nothing, Vector{String} } # spec name: actions
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function ServiceDelegationPropertiesFormat(;serviceName=nothing, actions=nothing, provisioningState=nothing)
        o = new()
        validate_property(ServiceDelegationPropertiesFormat, Symbol("serviceName"), serviceName)
        setfield!(o, Symbol("serviceName"), serviceName)
        validate_property(ServiceDelegationPropertiesFormat, Symbol("actions"), actions)
        setfield!(o, Symbol("actions"), actions)
        validate_property(ServiceDelegationPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ServiceDelegationPropertiesFormat

const _property_map_ServiceDelegationPropertiesFormat = Dict{Symbol,Symbol}(Symbol("serviceName")=>Symbol("serviceName"), Symbol("actions")=>Symbol("actions"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ServiceDelegationPropertiesFormat = Dict{Symbol,String}(Symbol("serviceName")=>"String", Symbol("actions")=>"Vector{String}", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ ServiceDelegationPropertiesFormat }) = collect(keys(_property_map_ServiceDelegationPropertiesFormat))
Swagger.property_type(::Type{ ServiceDelegationPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ServiceDelegationPropertiesFormat[name]))}
Swagger.field_name(::Type{ ServiceDelegationPropertiesFormat }, property_name::Symbol) =  _property_map_ServiceDelegationPropertiesFormat[property_name]

function check_required(o::ServiceDelegationPropertiesFormat)
    true
end

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


text:
function eos(model::vdWModel, V, T, z=@SVector [1.0])
    n = sum(z)
    return R̄*n*T * (a_ideal(idealmodel(model),V,T,z)+ a_res(model,V,T,z))
end

function a_res(model::ABCubicModel, V, T, z=@SVector [1.0])
    n = sum(z)
    n⁻¹ = 1/n
    x = z.*n⁻¹
    v = V*n⁻¹
    return a_resx(model,v,T,x)
end

#fast shortcut to evaluate cubics, pressure is known
#=
function ∂f∂v(model::ABCubicModel,v,t,z)
    #@info "fast shortcut"
    a,b,p = cubic_abp(model,v,t,z)
    return -p
end
=#
function second_virial_coefficient(model::ABCubicModel,T,z = SA[1.0])
    #@info "fast shortcut"
    a,b = cubic_ab(model,T,z)
    return b-a/(R̄*T)
end

function x0_volume_sc(model::ABCubicModel,p,T,z)
    Zc = cubic_zc(model)
    return Zc*R̄*T/p
end

function volume(model::ABCubicModel,p,T,z=SA[1.0];phase="unknown")
    lb_v   =lb_volume(model,z)
    xx = z/sum(z)
    RTp = R̄*T/p
    _poly = cubic_poly(model,p,T,z)

    #sols = Solvers.poly3(_poly)
    sols = Solvers.roots3(_poly)
    #xa = real.(sols) .* RTp
    #@show xa
    #with this, real sol is first
    function imagfilter(x)
        absx = abs(imag(x))
        return absx <8*eps(typeof(absx))
    end

    isreal = map(imagfilter,sols)
    _sols = sort(real.(sols)) .* RTp

    #@show _sols[isreal]
    #@show model.params.b.values
    #@show sum(isreal)
    if sum(isreal) == 3 #3 roots
        sols = sort(real.(sols))
        vg = last(sols)*RTp
        _vl = first(sols)*RTp
        vl = ifelse(_vl>lb_v,_vl,vg)
    elseif  sum(isreal) == 1
        i = findfirst(imagfilter,sols)
        vl = real(sols[i])*RTp
        vg = real(sols[i])*RTp
    elseif  sum(isreal) == 0
        @show sols
        error("weird")
    end

    #fp(_v) = pressure(model,_v,T,z) - p
    #@show vl
    #@show vg
    function gibbs(v)
        _df,f =  ∂f(model,v,T,z)
        dv,dt = _df
        if abs((p+dv)/p) > 0.03
            return Inf
        else
            return f  +p*v
        end
    end
    #this catches the supercritical phase as well
    if vl ≈ vg
        return vl
    end
    if phase == "unknown"
        gg = gibbs(vg)
        gl = gibbs(vl)
        #@show vg,vl
        #@show gg,gl
        return ifelse(gg<gl,vg,vl)
    elseif is_liquid(phase)
        return vl
    elseif is_vapour(phase)
        return vg
    else
        gg = gibbs(vg)
        gl = gibbs(vl)
        #@show vg,vl
        #@show gg,gl
        return ifelse(gg<gl,vg,vl)
    end
end

function ab_consts(model::CubicModel)
    return ab_consts(typeof(model))
end


text:
using LinearAlgebra
@testset "compress" begin
    A = [3 4.0; 0.1 -0.3]

    ## Test dangling node removal
    graph = Compgraph()
    add_lincomb!(graph, :T1, 1.0, :I, 3.0, :A)
    add_mult!(graph, :T2, :T1, :T1)

    # This node is dangling
    add_mult!(graph, :T3, :T2, :A)

    add_output!(graph, :T2)

    graph1 = graph
    graph2 = deepcopy(graph)

    compress_graph_dangling!(graph2)
    @test eval_graph(graph1, A) == eval_graph(graph2, A)

    ## Test zero coeff removal (and dangling)
    graph = Compgraph()
    add_mult!(graph, :A2, :A, :A)
    add_mult!(graph, :A4, :A2, :A2)
    add_mult!(graph, :A8, :A4, :A4)
    add_sum!(graph, :T, [3; 4; 0.1; 0; 0], [:I; :A; :A2; :A4; :A8])
    add_output!(graph, :T)

    graph1 = graph
    graph2 = deepcopy(graph)
    compress_graph!(graph2)
    @test eval_graph(graph1, A) == eval_graph(graph2, A)

    # Check that we actually removed something
    @test size(get_all_cref(graph2), 1) < size(get_all_cref(graph1), 1)

    ## Test function to detect if graph has linear combination of identities
    graph = Compgraph()
    add_lincomb!(graph, :P1, 1.0, :A, 1.0, :I)
    add_lincomb!(graph, :P2, 1.0, :I, 1.0, :I)
    add_mult!(graph, :P0, :P1, :P2)
    @test has_identity_lincomb(graph) == true

    ## Test trivial node detection and removal
    graph = Compgraph()
    add_mult!(graph, :AI, :A, :I)
    add_output!(graph, :AI)
    @test has_trivial_nodes(graph) == true

    graph = Compgraph()
    add_mult!(graph, :IA, :I, :A)
    add_output!(graph, :IA)
    @test has_trivial_nodes(graph) == true

    graph = Compgraph()
    add_ldiv!(graph, :IinvA, :I, :A)
    add_output!(graph, :IinvA)
    @test has_trivial_nodes(graph) == true

    graph = Compgraph()
    add_mult!(graph, :I2, :I, :I) # This node is trivial
    add_mult!(graph, :AI, :A, :I2) # This node is trivial
    add_mult!(graph, :IA, :I2, :AI) # This node is trivial
    add_ldiv!(graph, :P0, :I, :IA) # This node is trivial
    add_output!(graph, :P0)
    graph1 = graph
    graph2 = deepcopy(graph)
    @test has_trivial_nodes(graph1) == true
    cref = get_all_cref(graph) # Empty as graph has no linear combination nodes.
    compress_graph_trivial!(graph2, cref)
    @test isempty(cref)
    @test eval_graph(graph1, A) == eval_graph(graph2, A)
    @test isempty(graph2.operations)
    @test isempty(graph2.parents)
    @test isempty(graph2.coeffs)
    @test graph2.outputs == [:A]

    graph = Compgraph()
    add_lincomb!(graph, :P1, 1.0, :A, 2.0, :I)
    add_mult!(graph, :P2a, :A, :I)
    add_mult!(graph, :P2b, :I, :A)
    add_mult!(graph, :P2, :P2a, :P2b)
    add_mult!(graph, :P0, :P1, :P2)
    add_output!(graph, :P0)
    graph1 = graph
    graph2 = deepcopy(graph)
    @test has_trivial_nodes(graph1) == true
    cref = get_all_cref(graph) # Empty as graph has no linear combination nodes.
    compress_graph_trivial!(graph2, cref)
    @test cref == [(:P1, 1), (:P1, 2)]
    @test graph2.operations == Dict(:P0 => :mult, :P1 => :lincomb, :P2 => :mult)
    @test graph2.parents ==
          Dict(:P0 => (:P1, :P2), :P1 => (:A, :I), :P2 => (:A, :A))
    @test graph2.outputs == [:P0]

    ## Test function to remove redundant nodes.
    graph = Compgraph()
    add_lincomb!(graph, :AI1, 1.0, :A, 2.0, :I)
    add_lincomb!(graph, :AI2, 1.0, :A, 2.0, :I)
    add_lincomb!(graph, :IA1, 2.0, :I, 1.0, :A)
    add_mult!(graph, :P1, :AI1, :AI2)
    add_mult!(graph, :P2, :AI2, :AI1)
    add_ldiv!(graph, :D1, :AI1, :P2)
    add_ldiv!(graph, :D2, :AI1, :P2)
    add_mult!(graph, :out, :D1, :D2)
    add_output!(graph, :out)

    cref = get_all_cref(graph)

    graph1 = deepcopy(graph)
    cref1 = get_all_cref(graph1)
    compress_graph_redundant!(graph1, cref, compress_lincomb = false)
    @test length(get_sorted_keys(graph1)) == 6
    @test cref == cref1
    @test eval_graph(graph1, A) == eval_graph(graph, A)

    graph1 = deepcopy(graph)
    cref1 = get_all_cref(graph1)
    compress_graph_redundant!(graph1, cref1, compress_lincomb = true)
    @test length(get_sorted_keys(graph1)) == 4
    @test cref1 == [(:AI1, 1), (:AI1, 2)]
    @test eval_graph(graph1, A) == eval_graph(graph, A)

    graph1 = deepcopy(graph)
    cref1 = get_all_cref(graph1)
    compress_graph_redundant!(graph1, cref1)
    @test length(get_sorted_keys(graph1)) == 4
    @test cref1 == [(:AI1, 1), (:AI1, 2)]
    @test eval_graph(graph1, A) == eval_graph(graph, A)

    add_output!(graph, :D2) # Now :D2 cannot be merged with :D1.
    graph1 = deepcopy(graph)
    cref1 = get_all_cref(graph1)
    compress_graph_redundant!(graph1, cref1)
    @test length(get_sorted_keys(graph1)) == 5
    @test cref1 == [(:AI1, 1), (:AI1, 2)]
    @test eval_graph(graph1, A) == eval_graph(graph, A)
end


text:
println("===============================================")
println("Testing definition.jl")
println("===============================================")

function _RheoTimeData_explicit_nolog()
    a = Vector{RheoFloat}(1.0:1.0:3.0)
    b = Vector{RheoFloat}(4.0:1.0:6.0)
    c = Vector{RheoFloat}(7.0:1.0:9.0)

    data = RheoTimeData(a, b, c, nothing)

    data.σ==a && data.ϵ==b && data.t==c && isnothing(data.log)
end
@test _RheoTimeData_explicit_nolog()

function _RheoTimeData_const_nolog()
    eps = Vector{RheoFloat}(1.0:1.0:3.0)
    sig = Vector{RheoFloat}(4.0:1.0:6.0)
    t = Vector{RheoFloat}(7.0:1.0:9.0)

    data = RheoTimeData(ϵ=eps, σ=sig, t=t, savelog=false)

    data.σ==sig && data.ϵ==eps && data.t==t && isnothing(data.log)
end
@test _RheoTimeData_const_nolog()

function _RheoTimeData_const_nostrain()
    sig = Vector{RheoFloat}(4.0:1.0:6.0)
    t = Vector{RheoFloat}(7.0:1.0:9.0)

    data = RheoTimeData(σ=sig, t=t)

    data.σ==sig && data.ϵ==RheoFloat[] && data.t==t && data.log[1].info.type == stress_only::TimeDataType
end
@test _RheoTimeData_const_nostrain()

function _operators_logs()
    d1=strainfunction(timeline(0:0.5:10),t->exp(-t))
    d2=strainfunction(timeline(0:0.5:10),t->1-exp(-t))
    d=2*d1 - (-d2) + d2

    d3=rheologrun(d.log)

    (d3.ϵ == d.ϵ) && all([ abs(e-2.)<=eps(RheoFloat) for e in d.ϵ ])
end
@test _operators_logs()

function _union_strain1stress2()
    time_instance = timeline(t_start=0.0, t_end=15.0, step=0.2)
    imposed_strain = strainfunction(time_instance, sawtooth(offset=5.0, amp=2, period=5))
    imposed_stress = stressfunction(time_instance, stairs(offset=3.0, amp=0.5, width=4.0))

    combined = imposed_strain|imposed_stress

    (combined.σ==imposed_stress.σ) && (combined.ϵ==imposed_strain.ϵ) && (combined.t==imposed_stress.t)
end
@test _union_strain1stress2()

function _union_stress1strain2()
    time_instance = timeline(t_start=0.0, t_end=15.0, step=0.2)
    imposed_strain = strainfunction(time_instance, sawtooth(offset=5.0, amp=2, period=5))
    imposed_stress = stressfunction(time_instance, stairs(offset=3.0, amp=0.5, width=4.0))

    combined = imposed_stress|imposed_strain

    (combined.σ==imposed_stress.σ) && (combined.ϵ==imposed_strain.ϵ) && (combined.t==imposed_stress.t)
end
@test _union_stress1strain2()

function _freeze_params()
    SLS2_mod = freeze_params( SLS2, G₀=2, η₂=3.5)

    relaxmod(SLS2, 1, [2,1,2,3,3.5]) == relaxmod(SLS2_mod, 1, [1,2,3])
end
@test _freeze_params()

#
#   The moduli functions on arrays are extensively tested as part of the model tests
#

function _scalar_moduli()
    m=RheoModel(Spring, k=2)
    relaxmod(m, 1) == relaxmod(Spring, k=2, 1) == 2. && creepcomp(m, 1) == creepcomp(Spring, k=2, 1) == 0.5 &&  storagemod(m, 1) == storagemod(Spring, k=2, 1) == 2. && lossmod(m, 1) == lossmod(Spring, k=2, 1) == 0.0
end
@test _scalar_moduli()


text:
@doc raw"""
    CompressibleEulerEquations1D

The compressible Euler equations for an ideal gas in one space dimension.
"""
struct CompressibleEulerEquations1D <: AbstractCompressibleEulerEquations{1, 3}
  gamma::Float64
end

function CompressibleEulerEquations1D()
  gamma = parameter("gamma", 1.4)

  CompressibleEulerEquations1D(gamma)
end

get_name(::CompressibleEulerEquations1D) = "CompressibleEulerEquations1D"
varnames_cons(::CompressibleEulerEquations1D) = @SVector ["rho", "rho_v1", "rho_e"]
varnames_prim(::CompressibleEulerEquations1D) = @SVector ["rho", "v1", "p"]

# Set initial conditions at physical location `x` for time `t`
function initial_conditions_density_pulse(x, t, equation::CompressibleEulerEquations1D)
  rho = 1 + exp(-(x[1]^2 ))/2
  v1 = 1
  rho_v1 = rho * v1
  p = 1
  rho_e = p/(equation.gamma - 1) + 1/2 * rho * v1^2
  return @SVector [rho, rho_v1, rho_e]
end

"""
    initial_conditions_density_wave(x, t, equation::CompressibleEulerEquations1D)

test case for stability of EC fluxes from paper: https://arxiv.org/pdf/2007.09026.pdf
domain [-1, 1]
mesh = 4x4, polydeg = 5
"""

function initial_conditions_density_wave(x, t, equation::CompressibleEulerEquations1D)
  v1 = 0.1
  rho = 1 + 0.98 * sinpi(2 * (x[1] - t * v1))
  rho_v1 = rho * v1
  p = 20
  rho_e = p / (equation.gamma - 1) + 1/2 * rho * v1^2
  return @SVector [rho, rho_v1, rho_e]
end

function initial_conditions_constant(x, t, equation::CompressibleEulerEquations1D)
  rho = 1.0
  rho_v1 = 0.1
  rho_e = 10.0
  return @SVector [rho, rho_v1, rho_e]
end

function initial_conditions_convergence_test(x, t, equation::CompressibleEulerEquations1D)
  c = 2
  A = 0.1
  L = 2
  f = 1/L
  ω = 2 * pi * f
  ini = c + A * sin(ω * (x[1] - t))

  rho = ini
  rho_v1 = ini
  rho_e = ini^2

  return @SVector [rho, rho_v1, rho_e]
end

function initial_conditions_weak_blast_wave(x, t, equation::CompressibleEulerEquations1D)
  # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3)
  # Set up polar coordinates
  inicenter = SVector(0.0)
  x_norm = x[1] - inicenter[1]
  r = abs(x_norm)
  # The following code is equivalent to
  # phi = atan(0.0, x_norm)
  # cos_phi = cos(phi)
  # in 1D but faster
  cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm)

  # Calculate primitive variables
  rho = r > 0.5 ? 1.0 : 1.1691
  v1  = r > 0.5 ? 0.0 : 0.1882 * cos_phi
  p   = r > 0.5 ? 1.0 : 1.245

  return prim2cons(SVector(rho, v1, p), equation)
end

function initial_conditions_blast_wave(x, t, equation::CompressibleEulerEquations1D)
  # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
  # Set up polar coordinates
  inicenter = SVector(0.0)
  x_norm = x[1] - inicenter[1]
  r = abs(x_norm)
  # The following code is equivalent to
  # phi = atan(0.0, x_norm)
  # cos_phi = cos(phi)
  # in 1D but faster
  cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm)

  # Calculate primitive variables
  rho = r > 0.5 ? 1.0 : 1.1691
  v1  = r > 0.5 ? 0.0 : 0.1882 * cos_phi
  p   = r > 0.5 ? 1.0E-3 : 1.245

  return prim2cons(SVector(rho, v1, p), equation)
end

function initial_conditions_sedov_blast_wave(x, t, equation::CompressibleEulerEquations1D)
  # Set up polar coordinates
  inicenter = SVector(0.0)
  x_norm = x[1] - inicenter[1]
  r = abs(x_norm)

  # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000
  r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
  # r0 = 0.5 # = more reasonable setup
  E = 1.0
  p0_inner = 6 * (equation.gamma - 1) * E / (3 * pi * r0)
  p0_outer = 1.0e-5 # = true Sedov setup
  # p0_outer = 1.0e-3 # = more reasonable setup

  # Calculate primitive variables
  rho = 1.0
  v1  = 0.0
  p   = r > r0 ? p0_outer : p0_inner

  return prim2cons(SVector(rho, v1, p), equation)
end

# Apply boundary conditions
function boundary_conditions_convergence_test(u_inner, orientation, direction, x, t,
                                              surface_flux_function,
                                              equation::CompressibleEulerEquations1D)
  u_boundary = initial_conditions_convergence_test(x, t, equation)

  # Calculate boundary flux
  if direction == 2 # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equation)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equation)
  end

  return flux
end

# Apply source terms
function source_terms_convergence_test(ut, u, x, element_id, t, n_nodes, equation::CompressibleEulerEquations1D)
  # Same settings as in `initial_conditions`
  c = 2
  A = 0.1
  L = 2
  f = 1/L
  ω = 2 * pi * f
  γ = equation.gamma

  for i in 1:n_nodes
    x1 = x[1, i, element_id]

    si, co = sincos((t - x1)*ω)
    tmp = (-((4 * si * A - 4c) + 1) * (γ - 1) * co * A * ω) / 2

    ut[2, i, element_id] += tmp
    ut[3, i, element_id] += tmp

    # Original terms (without performanc enhancements)
    # ut[1, i, element_id] += 0
    # ut[2, i, element_id] += (-(((4 * sin((t - x1) * ω) * A - 4c) + 1)) *
    #                          (γ - 1) * cos((t - x1) * ω) * A * ω) / 2
    # ut[3, i, element_id] += (-(((4 * sin((t - x1) * ω) * A - 4c) + 1)) *
    #                          (γ - 1) * cos((t - x1) * ω) * A * ω) / 2
  end

  return nothing
end

# Calculate 1D flux for a single point
@inline function calcflux(u, orientation, equation::CompressibleEulerEquations1D)
  rho, rho_v1, rho_e = u
  v1 = rho_v1/rho
  p = (equation.gamma - 1) * (rho_e - 1/2 * rho * v1^2)
  # Ignore orientation since it is always "1" in 1D
  f1 = rho_v1
  f2 = rho_v1 * v1 + p
  f3 = (rho_e + p) * v1
  return SVector(f1, f2, f3)
end

"""
    function flux_shima_etal(u_ll, u_rr, orientation, equation::CompressibleEulerEquations1D)

This flux is is a modification of the original kinetic energy preserving two-point flux by
Kuya, Totani and Kawai (2018)
  Kinetic energy and entropy preserving schemes for compressible flows
  by split convective forms
  [DOI: 10.1016/j.jcp.2018.08.058](https://doi.org/10.1016/j.jcp.2018.08.058)
The modification is in the energy flux to guarantee pressure equilibrium and was developed by
Nao Shima, Yuichi Kuya, Yoshiharu Tamaki, Soshi Kawai (JCP 2020)
  Preventing spurious pressure oscillations in split convective form discretizations for
  compressible flows
"""
@inline function flux_shima_etal(u_ll, u_rr, orientation, equation::CompressibleEulerEquations1D)
  # Unpack left and right state
  rho_ll, rho_v1_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll / rho_ll
  v1_rr = rho_v1_rr / rho_rr
  p_ll  = (equation.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * (v1_ll^2))
  p_rr  = (equation.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * (v1_rr^2))

  # Average each factor of products in flux
  rho_avg = 1/2 * (rho_ll + rho_rr)
  v1_avg  = 1/2 * ( v1_ll +  v1_rr)
  p_avg   = 1/2 * (  p_ll +   p_rr)
  kin_avg = 1/2 * (v1_ll*v1_rr)

  # Calculate fluxes
  # Ignore orientation since it is always "1" in 1D
  pv1_avg = 1/2 * (p_ll*v1_rr + p_rr*v1_ll)
  f1 = rho_avg * v1_avg
  f2 = rho_avg * v1_avg * v1_avg + p_avg
  f3 = p_avg*v1_avg/(equation.gamma-1) + rho_avg*v1_avg*kin_avg + pv1_avg

  return SVector(f1, f2, f3)
end

"""
    flux_kennedy_gruber(u_ll, u_rr, orientation, equation::CompressibleEulerEquations1D)

Kinetic energy preserving two-point flux by Kennedy and Gruber (2008)
  Reduced aliasing formulations of the convective terms within the
  Navier-Stokes equations for a compressible fluid
[DOI: 10.1016/j.jcp.2007.09.020](https://doi.org/10.1016/j.jcp.2007.09.020)
"""
@inline function flux_kennedy_gruber(u_ll, u_rr, orientation, equation::CompressibleEulerEquations1D)
  # Unpack left and right state
  rho_ll, rho_v1_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll/rho_ll
  v1_rr = rho_v1_rr/rho_rr

  # Average each factor of products in flux
  rho_avg = 1/2 * (rho_ll + rho_rr)
  v1_avg = 1/2 * (v1_ll + v1_rr)
  p_avg = 1/2 * ((equation.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * (v1_ll^2)) +
                 (equation.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * (v1_rr^2)))
  e_avg = 1/2 * (rho_e_ll/rho_ll + rho_e_rr/rho_rr)

  # Ignore orientation since it is always "1" in 1D
  f1 = rho_avg * v1_avg
  f2 = rho_avg * v1_avg * v1_avg + p_avg
  f3 = (rho_avg * e_avg + p_avg) * v1_avg

  return SVector(f1, f2, f3)
end

"""
    flux_chandrashekar(u_ll, u_rr, orientation, equation::CompressibleEulerEquations2D)

Entropy conserving two-point flux by Chandrashekar (2013)
  Kinetic Energy Preserving and Entropy Stable Finite Volume Schemes
  for Compressible Euler and Navier-Stokes Equations
[DOI: 10.4208/cicp.170712.010313a](https://doi.org/10.4208/cicp.170712.010313a)
"""
@inline function flux_chandrashekar(u_ll, u_rr, orientation, equation::CompressibleEulerEquations1D)
  # Unpack left and right state
  rho_ll, rho_v1_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll/rho_ll
  v1_rr = rho_v1_rr/rho_rr
  p_ll =  (equation.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * (v1_ll^2))
  p_rr =  (equation.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * (v1_rr^2))
  beta_ll = 0.5*rho_ll/p_ll
  beta_rr = 0.5*rho_rr/p_rr
  specific_kin_ll = 0.5*(v1_ll^2)
  specific_kin_rr = 0.5*(v1_rr^2)

  # Compute the necessary mean values
  rho_avg  = 0.5*(rho_ll+rho_rr)
  rho_mean = ln_mean(rho_ll,rho_rr)
  beta_mean = ln_mean(beta_ll,beta_rr)
  beta_avg = 0.5*(beta_ll+beta_rr)
  v1_avg = 0.5*(v1_ll+v1_rr)
  p_mean = 0.5*rho_avg/beta_avg
  velocity_square_avg = specific_kin_ll + specific_kin_rr

  # Calculate fluxes
  # Ignore orientation since it is always "1" in 1D
  f1 = rho_mean * v1_avg
  f2 = f1 * v1_avg + p_mean
  f3 = f1 * 0.5*(1/(equation.gamma-1)/beta_mean - velocity_square_avg)+f2*v1_avg

  return SVector(f1, f2, f3)
end

"""
    flux_ranocha(u_ll, u_rr, orientation, equation::CompressibleEulerEquations2D)

Entropy conserving and kinetic energy preserving two-point flux by Ranocha (2018)
  Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods
  for Hyperbolic Balance Laws
[PhD thesis, TU Braunschweig](https://cuvillier.de/en/shop/publications/7743)
See also Ranocha (2020)
  Entropy Conserving and Kinetic Energy Preserving Numerical Methods for
  the Euler Equations Using Summation-by-Parts Operators
[Proceedings of ICOSAHOM 2018](https://doi.org/10.1007/978-3-030-39647-3_42)
"""
@inline function flux_ranocha(u_ll, u_rr, orientation, equation::CompressibleEulerEquations1D)
  # Unpack left and right state
  rho_ll, rho_v1_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll / rho_ll
  v1_rr = rho_v1_rr / rho_rr
  p_ll =  (equation.gamma - 1) * (rho_e_ll - 0.5 * rho_ll * (v1_ll^2))
  p_rr =  (equation.gamma - 1) * (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 ))

  # Compute the necessary mean values
  rho_mean   = ln_mean(rho_ll, rho_rr)
  rho_p_mean = ln_mean(rho_ll / p_ll, rho_rr / p_rr)
  v1_avg = 0.5 * (v1_ll + v1_rr)
  p_avg  = 0.5 * (p_ll + p_rr)
  velocity_square_avg = 0.5 * (v1_ll*v1_rr)

  # Calculate fluxes
  # Ignore orientation since it is always "1" in 1D
  f1 = rho_mean * v1_avg
  f2 = f1 * v1_avg + p_avg
  f3 = f1 * ( velocity_square_avg + 1 / ((equation.gamma-1) * rho_p_mean) ) + 0.5 * (p_ll*v1_rr + p_rr*v1_ll)

  return SVector(f1, f2, f3)
end

function flux_lax_friedrichs(u_ll, u_rr, orientation, equation::CompressibleEulerEquations1D)
  # Calculate primitive variables and speed of sound
  rho_ll, rho_v1_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll / rho_ll
  v_mag_ll = abs(v1_ll)
  p_ll = (equation.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * v_mag_ll^2)
  c_ll = sqrt(equation.gamma * p_ll / rho_ll)
  v1_rr = rho_v1_rr / rho_rr
  v_mag_rr = abs(v1_rr)
  p_rr = (equation.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * v_mag_rr^2)
  c_rr = sqrt(equation.gamma * p_rr / rho_rr)

  # Obtain left and right fluxes
  f_ll = calcflux(u_ll, orientation, equation)
  f_rr = calcflux(u_rr, orientation, equation)

  λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr)
  f1 = 1/2 * (f_ll[1] + f_rr[1]) - 1/2 * λ_max * (rho_rr    - rho_ll)
  f2 = 1/2 * (f_ll[2] + f_rr[2]) - 1/2 * λ_max * (rho_v1_rr - rho_v1_ll)
  f3 = 1/2 * (f_ll[3] + f_rr[3]) - 1/2 * λ_max * (rho_e_rr  - rho_e_ll)

  return SVector(f1, f2, f3)
end

function flux_hll(u_ll, u_rr, orientation, equation::CompressibleEulerEquations1D)
  # Calculate primitive variables and speed of sound
  rho_ll, rho_v1_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll / rho_ll
  p_ll = (equation.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * v1_ll^2)

  v1_rr = rho_v1_rr / rho_rr
  p_rr = (equation.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * v1_rr^2)

  # Obtain left and right fluxes
  f_ll = calcflux(u_ll, orientation, equation)
  f_rr = calcflux(u_rr, orientation, equation)

  Ssl = v1_ll - sqrt(equation.gamma * p_ll / rho_ll)
  Ssr = v1_rr + sqrt(equation.gamma * p_rr / rho_rr)

  if Ssl >= 0.0 && Ssr > 0.0
    f1 = f_ll[1]
    f2 = f_ll[2]
    f3 = f_ll[3]
  elseif Ssr <= 0.0 && Ssl < 0.0
    f1 = f_rr[1]
    f2 = f_rr[2]
    f3 = f_rr[3]
  else
    f1 = (Ssr*f_ll[1] - Ssl*f_rr[1] + Ssl*Ssr*(rho_rr[1]    - rho_ll[1]))/(Ssr - Ssl)
    f2 = (Ssr*f_ll[2] - Ssl*f_rr[2] + Ssl*Ssr*(rho_v1_rr[1] - rho_v1_ll[1]))/(Ssr - Ssl)
    f3 = (Ssr*f_ll[3] - Ssl*f_rr[3] + Ssl*Ssr*(rho_e_rr[1]  - rho_e_ll[1]))/(Ssr - Ssl)
  end

  return SVector(f1, f2, f3)
end

# Determine maximum stable time step based on polynomial degree and CFL number
function calc_max_dt(u, element_id, invjacobian, cfl,
                     equation::CompressibleEulerEquations1D, dg)
  λ_max = 0.0
  for i in 1:nnodes(dg)
    rho, rho_v1, rho_e = get_node_vars(u, dg, i, element_id)
    v1 = rho_v1 / rho
    v_mag = abs(v1)
    p = (equation.gamma - 1) * (rho_e - 1/2 * rho * v_mag^2)
    c = sqrt(equation.gamma * p / rho)
    λ_max = max(λ_max, v_mag + c)
  end

  dt = cfl * 2 / (nnodes(dg) * invjacobian * λ_max)

  return dt
end

# Convert conservative variables to primitive
@inline function cons2prim(u, equation::CompressibleEulerEquations1D)
  rho, rho_v1, rho_e = u

  v1 = rho_v1 / rho
  p = (equation.gamma - 1) * (rho_e - 0.5 * rho * (v1^2))

  return SVector(rho, v1, p)
end

# Convert conservative variables to entropy
@inline function cons2entropy(u, equation::CompressibleEulerEquations1D)
  rho, rho_v1, rho_e = u

  v1 = rho_v1 / rho
  v_square = v1^2
  p = (equation.gamma - 1) * (rho_e - 0.5 * rho * v_square)
  s = log(p) - equation.gamma*log(rho)
  rho_p = rho / p

  w1 = (equation.gamma - s) / (equation.gamma-1) - 0.5 * rho_p * v_square
  w2 = rho_p * v1
  w3 = -rho_p

  return SVector(w1, w2, w3)
end

# Convert primitive to conservative variables
@inline function prim2cons(prim, equation::CompressibleEulerEquations1D)
  rho, v1, p = prim
  rho_v1 = rho * v1
  rho_e  = p/(equation.gamma-1) + 0.5 * (rho_v1 * v1)
  return SVector(rho, rho_v1, rho_e)
end

@inline function density(u, equation::CompressibleEulerEquations1D)
 rho = u[1]
 return rho
end

@inline function pressure(u, equation::CompressibleEulerEquations1D)
 rho, rho_v1, rho_e = u
 p = (equation.gamma - 1) * (rho_e - 0.5 * (rho_v1^2) / rho)
 return p
end

@inline function density_pressure(u, equation::CompressibleEulerEquations1D)
 rho, rho_v1, rho_e = u
 rho_times_p = (equation.gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2))
 return rho_times_p
end

# Calculate thermodynamic entropy for a conservative state `cons`
@inline function entropy_thermodynamic(cons, equation::CompressibleEulerEquations1D)
  # Pressure
  p = (equation.gamma - 1) * (cons[3] - 1/2 * (cons[2]^2) / cons[1])

  # Thermodynamic entropy
  s = log(p) - equation.gamma*log(cons[1])

  return s
end

# Calculate mathematical entropy for a conservative state `cons`
@inline function entropy_math(cons, equation::CompressibleEulerEquations1D)
  # Mathematical entropy
  S = -entropy_thermodynamic(cons, equation) * cons[1] / (equation.gamma - 1)

  return S
end

# Default entropy is the mathematical entropy
@inline entropy(cons, equation::CompressibleEulerEquations1D) = entropy_math(cons, equation)

# Calculate total energy for a conservative state `cons`
@inline energy_total(cons, ::CompressibleEulerEquations1D) = cons[3]

# Calculate kinetic energy for a conservative state `cons`
@inline function energy_kinetic(cons, equation::CompressibleEulerEquations1D)
  return 0.5 * (cons[2]^2)/cons[1]
end

# Calculate internal energy for a conservative state `cons`
@inline function energy_internal(cons, equation::CompressibleEulerEquations1D)
  return energy_total(cons, equation) - energy_kinetic(cons, equation)
end


text:
using Printf

export isip

function isip(str::AbstractString, version::Int=0)::Bool
    ipv4SegmentFormat = "(?:[0-9]|[1-9][0-9]|1[0-9][0-9]|2[0-4][0-9]|25[0-5])"
    ipv4AddressFormat = @sprintf("(%s[.]){3}%s", ipv4SegmentFormat, ipv4SegmentFormat)
    ipv4AddressReg = Regex("^$(ipv4AddressFormat)\$")

    ipv6SegmentFormat = "(?:[0-9a-fA-F]{1,4})"
    ipv6AddressReg = Regex("^(" *
        "(?:$(ipv6SegmentFormat):){7}(?:$(ipv6SegmentFormat)|:)|" *
        "(?:$(ipv6SegmentFormat):){6}(?:$(ipv4AddressFormat)|:$(ipv6SegmentFormat)|:)|" *
        "(?:$(ipv6SegmentFormat):){5}(?::$(ipv4AddressFormat)|(:$(ipv6SegmentFormat)){1,2}|:)|" *
        "(?:$(ipv6SegmentFormat):){4}(?:(:$(ipv6SegmentFormat)){0,1}:$(ipv4AddressFormat)|(:$(ipv6SegmentFormat)){1,3}|:)|" *
        "(?:$(ipv6SegmentFormat):){3}(?:(:$(ipv6SegmentFormat)){0,2}:$(ipv4AddressFormat)|(:$(ipv6SegmentFormat)){1,4}|:)|" *
        "(?:$(ipv6SegmentFormat):){2}(?:(:$(ipv6SegmentFormat)){0,3}:$(ipv4AddressFormat)|(:$(ipv6SegmentFormat)){1,5}|:)|" *
        "(?:$(ipv6SegmentFormat):){1}(?:(:$(ipv6SegmentFormat)){0,4}:$(ipv4AddressFormat)|(:$(ipv6SegmentFormat)){1,6}|:)|" *
        "(?::((?::$(ipv6SegmentFormat)){0,5}:$(ipv4AddressFormat)|(?::$(ipv6SegmentFormat)){1,7}|:))" *
        ")(%[0-9a-zA-Z-.:]{1,})?\$")

    if version === 4
        if !contains(str, ipv4AddressReg)
            return false
        end

        return maximum(parse.(Int, split(str, '.'))) <= 255
    elseif version === 6
        return contains(str, ipv6AddressReg)
    else
        return isip(str, 4) || isip(str, 6)
    end
end