text:
# initialisation functions for GeMM

"""
    createpop(settings)

Create a new, random individual and replicates it a certain number of times
(depending on metabolic variables) to create a new population of organisms.
Returns an array of individuals.
"""
function createpop(settings::Dict{String, Any})
    # start with an archetypical individual
    archetype = createind(settings, true)
    population = [archetype]
    # calculate the population size
    popsize = 0
    if occursin("metabolic", settings["popsize"]) || occursin("single", settings["popsize"])
        # population size determined by adult size
        popsize = round(settings["fertility"] * archetype.traits["repsize"] ^ (-1 / 4) *
                        exp(-act / (boltz * 298.0)))
    elseif occursin("bodysize", settings["popsize"])
        # population size up to 25% of the maximum possible in this cell
        quarterpopsize = Integer(ceil((settings["cellsize"] / archetype.traits["repsize"]) / 4))
        popsize = rand(1:quarterpopsize) + 1
    elseif occursin("minimal", settings["popsize"]) || popsize == 0
        popsize = 2 #Takes two to tangle ;-)
    else
        simlog("Invalid value for `popsize`: $(settings["popsize"])", 'e')
    end
    # clone & mutate the archetype N times
    locivar = rand()
    for i in 2:popsize
        #TODO avoid `deepcopy` here
        ind = deepcopy(archetype) # XXX this is potentially very slow!
        ind.id = rand(UInt32)
        varyalleles!(ind.genome, settings, locivar)
        ind.traits = gettraitdict(ind.genome, settings["traitnames"])
        if !(settings["indsize"] == "adult" || settings["indsize"] == "seed")
            #XXX sizes uniformally distributed?
            ind.size = ind.traits["seedsize"] + rand() * (ind.traits["repsize"] - ind.traits["seedsize"])
        end
        push!(population, ind)
    end
    population
end

"""
    createind(settings, marked=false)

Create an individual organism of a new species with a random genome.
"""
function createind(settings::Dict{String, Any}, marked::Bool = false)
    id = rand(Int32)
    lineage = randstring(6)
    ngenes = rand(1:settings["maxloci"]) * length(settings["traitnames"])
    settings["degpleiotropy"] == 0 && (ngenes = length(settings["traitnames"]))
    ngenes < 1 && (ngenes = 1)
    traits = createtraits(settings)
    genes = creategenes(ngenes, traits, settings)
    if settings["linkage"] == "none"
        nchrms = length(genes)
    elseif settings["linkage"] == "full"
        nchrms = 1
    else
        linkage = rand(1:length(genes))
        nchrms = Integer(round(ngenes/linkage))
    end
    chromosomes = createchrms(nchrms, genes)
    locivar = rand()
    varyalleles!(chromosomes, settings, locivar)
    traitdict = gettraitdict(chromosomes, settings["traitnames"])
    if settings["indsize"] == "adult"
        indsize = traitdict["repsize"]
    elseif settings["indsize"] == "seed"
        indsize = traitdict["seedsize"]
    else
        #XXX sizes uniformally distributed?
        indsize = traitdict["seedsize"] + rand() * (traitdict["repsize"] - traitdict["seedsize"])
    end
    Individual(lineage, chromosomes, traitdict, marked, 1.0, 1.0, indsize, hermaphrodite, 0, id)
end

"""
    genesis(settings)

Create a new community, composed of random new species populations, for a patch.
Returns an array of individuals.
"""
function genesis(settings::Dict{String, Any})
    community = Individual[]
    totalmass = 0.0
    while true
        population = createpop(settings)
        popsize = length(population)
        # Check the cell capacity
        popmass = sum(map(x -> x.size, population))
        if totalmass + popmass > settings["cellsize"] * settings["overfill"]
            # stop loop if cell is full
            if totalmass >= settings["cellsize"] * 0.9 || occursin("single", settings["popsize"])
                #make sure the cell is full enough
                simlog("Cell is now $(round((totalmass/settings["cellsize"])*100))% full.", 'd') #DEBUG
                break
            else
                continue
            end
        end
        totalmass += popmass
        append!(community, population)
        occursin("single", settings["popsize"]) && break
    end
    simlog("Patch initialized with $(length(community)) individuals.", 'd') #DEBUG
    community
end

"""
    createpatch(patchentry, settings)

Create a new patch and set its parameters by parsing one line of the map file.
"""
function createpatch(patchentry::Array{String,1}, settings::Dict{String, Any})
    size(patchentry,1) < 3 && simlog("please check your map file for incomplete or faulty entries. \n
Each line must contain patch information with at least \n
\t - a unique integer ID, \n
\t - an integer x coordinate, \n
\t - an integer y coordinate, \n
separated by a whitespace character (<ID> <x> <y>).", settings, 'e')
    # create the basic patch
    id = parse(Int, patchentry[1])
    xcord = parse(Int, patchentry[2])
    ycord = parse(Int, patchentry[3])
    capacity = settings["cellsize"]
    # XXX the 'global' here is a hack so that I can use eval() later on
    # (eval() always works on the global scope)
    global newpatch = Patch(id, (xcord, ycord), capacity)
    simlog("Creating patch $id at $xcord/$ycord, size $capacity", 'd') #DEBUG
    # parse other parameter options
    for p in patchentry[4:end]
        varval = split(p, '=')
        var = varval[1]
        if !(var in map(string, fieldnames(Patch)))
            simlog("Unrecognized patch parameter $var.", 'w')
            continue
        elseif length(varval) < 2
            # if no value is specified, assume 'true'
            # (unless the variable is prepended with an exclamation mark, then take 'false')
            if var[1] == '!'
                var = var[2:end]
                val = false
            else
                val = true
            end
        else
            val = Meta.parse(varval[2])
        end
        # check for correct type and modify the new patch
        vartype = typeof(eval(Meta.parse("newpatch."*var)))
        if !isa(val, vartype)
            try
                val = convert(vartype, val)
            catch
                simlog("Invalid patch parameter type $var: $val", 'w')
                continue
            end
        end
        eval(Meta.parse("newpatch."*string(var)*" = $val"))
    end
    newpatch
end

"""
    createworld(maptable, settings)

Use a parsed map file (as returned by `readmapfile`) to create the world.
Initialises each patch with its parameters and a new community, then returns
an array of patches.
"""
function createworld(maptable::Array{Array{String,1},1}, settings::Dict{String, Any})
    simlog("Creating world...")
    world = Patch[]
    for entry in maptable
        newpatch = createpatch(entry, settings)
        if newpatch.initpop
            if settings["mode"] == "zosterops"
                append!(newpatch.community, zgenesis(newpatch, settings))
            elseif settings["indsize"] != "seed"
                append!(newpatch.community, genesis(settings))
            else
                append!(newpatch.seedbank, genesis(settings))
            end
        end
        push!(world, newpatch)
        global newpatch = nothing #clear memory used in `createpatch()` #XXX does this make sense?
    end
    world
end

"""
    updateworld!(world, maptable, cellsize)

Reinitialise the world from another parsed map file. Works analogously to
`createworld`. Intended for use in scenarios where the model world changes
during a run (e.g. through global warming or island ontogeny).
"""
#FIXME this function needs to be folded into createworld(), it duplicates way too much code
function updateworld!(world::Array{Patch,1},maptable::Array{Array{String,1},1}, settings::Dict{String, Any})
    cellsize = settings["cellsize"]
    simlog("Updating world...")
    allids = Int[]
    for entry in maptable
        size(entry,1) < 3 && error("please check your map file for incomplete or faulty entries. \n
                            Each line must contain patch information with at least \n
                            \t - a unique integer ID, \n
                            \t - an integer x coordinate, \n
                            \t - an integer y coordinate, \n
                            separated by a whitespace character (<ID> <x> <y>).")
        id = parse(Int, entry[1])
        push!(allids, id)
        xcord = parse(Int, entry[2])
        ycord = parse(Int, entry[3])
        # XXX the 'global' here is a hack so that I can use eval() later on
        # (this always works on the global scope)
        idx = findall(x -> x.id == id, world)
        if length(idx) == 0
            marked = true
            global newpatch = Patch(id, (xcord, ycord), cellsize)
        else
            marked = false
            world[idx[1]].location = (xcord, ycord)
            global newpatch = world[idx[1]]
        end
        # parse other parameter options
        for p in entry[4:end]
            varval = split(p, '=')
            var = varval[1]
            if !(var in map(string, fieldnames(Patch)))
                simlog("Unrecognized patch parameter $var.", 'w')
                continue
            elseif length(varval) < 2
                val = true # if no value is specified, assume 'true'
            else
                val = Meta.parse(varval[2])
            end
            # check for correct type and modify the new patch
            vartype = typeof(eval(Meta.parse("newpatch."*var)))
            if !isa(val, vartype)
                try
                    val = convert(vartype, val)
                catch
                    simlog("Invalid patch parameter type $var: $val", 'w')
                    continue
                end
            end
            eval(Meta.parse("newpatch."*string(var)*" = $val"))
        end
        if marked
            push!(world, newpatch)
            global newpatch = nothing #clear memory
        end
    end
    filter!(x -> x.id in allids, world)
    world
end

max_stars_repo_name: lleiding/gemm

max_issues_repo_name: lleiding/gemm

max_forks_repo_name: lleiding/gemm


text:
#script for creating Octree grids for POLARIS
#right now a clumpy AGN-Torus is modelled. The underlying continous model is a wedge rotated around the z-axis
#with density(R) = rho_0 * R**-0.5

using PyCall		#using python structs for easy IO compatible with POLARIS
s = pyimport("struct")
using Distributed	#the children of a cell don't depend on each other -> parallelize
addprocs(4)    		#right amount of processes on my machine, should probably be abstracted
using ArgParse
using Random

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "-N"
            help = "Number of Clumps in Torus"
            arg_type = Int
            default = 10^3
        "-S"
            help = "Seed for random generator"
            arg_type = Int
            default = 0
        "-R"
            help = "Clump radius in pc"
            arg_type = Float64
            default = 1.0
        end
    return parse_args(s)
end

@everywhere mutable struct cell
    x::Float64
    y::Float64
    z::Float64

    length::Float64
    level::UInt8
    isleaf::UInt8
    data::Array{Float64}
    branches::Array{cell}
end

@everywhere mutable struct octree
    root::cell
    mass::Float64
    maxlevel::UInt8
	minlevel::UInt8
    cellcounter
    totalmass
end

@everywhere function getcenter(c::cell)::Array{Float64}
    cx = c.x + 0.5 * c.length
    cy = c.y + 0.5 * c.length
    cz = c.z + 0.5 * c.length
    return [cx,cy,cz]
end

@everywhere function density(pos::Vector)::Float64
    dens = 0.0
	#problem with the model: clumps cover the the whole center for certain configurations
	#this ensures that a cylinder in z-direction is left free
	if sqrt(pos[1]^2 + pos[2]^2) < Rᵢ
        return 0.0
    end
    for i =1:length(view(positions,:,1))
        if distance(pos,view(positions,i,:)) <= R_cl
            dens += 1 #whole model gets normalized while writing the tree, so relativ density is fine here
        end
    end
    return dens
end

#current way of deciding if a cell is filled or split up:
#calculate mass in cell, calculate mass that would result if cell is split up again
#if difference is small enough fill cell
#cells near the center always get split up to the maximum level
@everywhere function refinement(c::cell)
	pos = getcenter(c)
	if sqrt(pos[1]^2 + pos[2]^2 + pos[3]^2) <= 3.0 * con_pc
		return false
	else
		mass = density(pos) * (sl/(2^(c.level)))^3
		massNextLevel = 0.
		l = c.length * 0.5
		level = c.level + 1
		v = (sl/(2^(level)))^3
		x = c.x
		y = c.y
		z = c.z
		massNextLevel += density([x + 0.5*l,y + 0.5*l,z + 0.5*l]) * v
		massNextLevel += density([x + l + 0.5*l,y + 0.5*l,z + 0.5*l]) * v
		massNextLevel += density([x + 0.5*l,y + l + 0.5*l,z + 0.5*l]) * v
		massNextLevel += density([x + l + 0.5*l,y + l + 0.5*l,z + 0.5*l]) * v
		massNextLevel += density([x + 0.5*l,y + 0.5*l,z + l + 0.5*l]) * v
		massNextLevel += density([x + l + 0.5*l,y + 0.5*l,z + l + 0.5*l]) * v
		massNextLevel += density([x + 0.5*l,y + l + 0.5*l,z + l + 0.5*l]) * v
		massNextLevel += density([x + l + 0.5*l,y + l + 0.5*l,z + l + 0.5*l]) * v
		if massNextLevel == 0.
			return true
		else
			return abs(massNextLevel - mass) / massNextLevel < 10^-4
		end
	end
end

@everywhere function distance(pos1::Vector,pos2)::Float64
    return sqrt((pos1[1] - pos2[1])^2 + (pos1[2] - pos2[2])^2 + (pos1[3] - pos2[3])^2)
end

@everywhere function fillcell(c::cell)::Float64
	mass = density(getcenter(c)) * c.length^3
    push!(c.data,mass)
    return mass
end

@everywhere function addNextLevel(c::cell,t::octree)
	if c.level == t.maxlevel || (c.level >= t.minlevel && refinement(c))
        c.isleaf = 1
        Threads.atomic_add!(t.cellcounter, 1)
        Threads.atomic_add!(t.totalmass, fillcell(c))
        if t.cellcounter[] % 1000 == 0
            println(t.cellcounter)
        end

    else
        x = c.x
        y = c.y
        z = c.z
        level = c.level + 1
        l = c.length * 0.5
        c1 = cell(x,y,z,l,level,0,Vector{Float64}(),Vector())
        push!(c.branches, c1)
        c2 = cell(x + l,y,z,l,level,0,Vector{Float64}(),Vector())
        push!(c.branches, c2)
        c3 = cell(x,y + l,z,l,level,0,Vector{Float64}(),Vector())
        push!(c.branches, c3)
        c4 = cell(x + l,y + l,z,l,level,0,Vector{Float64}(),Vector())
        push!(c.branches, c4)
        c5 = cell(x,y,z + l,l,level,0,Vector{Float64}(),Vector())
        push!(c.branches, c5)
        c6 = cell(x + l,y,z + l,l,level,0,Vector{Float64}(),Vector())
        push!(c.branches, c6)
        c7 = cell(x,y + l,z + l,l,level,0,Vector{Float64}(),Vector())
        push!(c.branches, c7)
        c8 = cell(x + l,y + l,z + l,l,level,0,Vector{Float64}(),Vector())
        push!(c.branches, c8)

        Threads.@threads for i=1:8
            addNextLevel(c.branches[i],t)
        end
    end
end

function writeTree(t::octree,c::cell,io)
    write(io,s.pack("H",c.isleaf))
    write(io,s.pack("H",c.level))
    if c.isleaf == 1
        gas_m = c.data[1]*t.mass /t.totalmass[]
        gas_d = gas_m / c.length^3
        write(io,s.pack("f",gas_d))
        for i in 1:15
            write(io,s.pack("f",gas_d * 0.01 * dust_fractions[i]))
        end
    else
        for i=1:8
            writeTree(t,c.branches[i],io)
        end
    end
end

###############---MAIN---###############

#Constants
const con_pc = 3.0856775814671916e+16
const con_msun = 1.9884754153381438e+30

#grid constants
const grid_id = 20
const nr_dust_densities = 15
const data_id = ones(Int,nr_dust_densities) .* 29
const data_len = length(data_id) + 1
const dust_fractions = [0.0492868,0.01971472,0.00985736,0.07288324,0.02915329,0.01457665,0.10777664,0.04311066,0.02155533,0.15937553,0.06375021,0.03187511,0.23567779,0.09427111,0.04713556]

#commandline arguments
parsed_args = parse_commandline()
const N = parsed_args["N"]
const S = parsed_args["S"]
const R_cl = parsed_args["R"] * con_pc
const outR = round(Int,parsed_args["R"])
Random.seed!(S)

#define modell parameters here
const sl = 100 * con_pc
const maxlevel = 8
const minlevel = 6
const modellmass = 2236036 * con_msun

#choose positions
const Rₒ = sl * 0.5
const Rᵢ = 0.4 * con_pc
const positions = zeros(Float64,N,3)

for i in 1:N
    θ = acos(2*rand() - 1)
    while true
        if θ >= π/4 && θ <= 3π/4
            break
        else
            θ = acos(2*rand() - 1)
        end
    end
    ϕ = 2π * rand()
    r = (rand() * ( Rₒ^(5/2) - Rᵢ^(5/2)) + Rᵢ^(5/2))^(2/5)

    positions[i,1] = r*sin(θ)cos(ϕ)
    positions[i,2] = r*sin(θ)sin(ϕ)
    positions[i,3] = r*cos(θ)
end

#create initial values
tmpmass = 0
cellcounter = Threads.Atomic{Int}(0)
totalmass = Threads.Atomic{Float64}(0.0)
root = cell(-0.5*sl,-0.5*sl,-0.5*sl,sl,0,0,Vector{Float64}(),Vector())
tree = octree(root,modellmass,maxlevel,minlevel,cellcounter,totalmass)

#build and write tree
addNextLevel(root,tree)
open("bigclumps_$(N).dat", "w") do io
    write(io,s.pack("H",grid_id))
    write(io,s.pack("H",data_len))
    write(io,s.pack("H",28))
    for i=1:data_len - 1
        write(io,s.pack("H",data_id[i]))
    end
    write(io,s.pack("d",sl))
    writeTree(tree,tree.root,io)
end

ext: jl

max_stars_repo_name: tstoertenbecker/octreeJl

max_issues_repo_name: tstoertenbecker/octreeJl

max_forks_repo_name: tstoertenbecker/octreeJl


text:
abstract type AbstractBlockHankel{T<:Number} <: AbstractMatrix{T} end

function size(A::AbstractBlockHankel)
    (size(A, 1), size(A, 2))
end

function getindex(A::AbstractBlockHankel, i::Integer)
    return A[mod(i - 1, size(A, 1)) + 1, div(i - 1, size(A, 1)) + 1]
end
struct BlockHankel{T<:Number} <: AbstractBlockHankel{T}
    w::Vector{Vector{T}}
    N = length(w)                                                                       # Length of the discrete-time vector signal
    m = length(w[1])                                                                    # Dimension of the (discrete-time) signal
    function BlockHankel{T}(w::Vector{Vector{T}},l::S) where {T<:Number,S<:Unsigned}    # w ... vector of vectors, l ... # of block rows
        if first(vc) != first(vr)
            error("First element of the vectors must be the same")
        end
        return new{T}(w)
    end
end

function Matrix(A::AbstractBlockHankel{T}) where T
    m, n = size(A)
    Af = Matrix{T}(undef, m, n)
    for j = 1:n
        for i = 1:m
            Af[i,j] = A[i,j]
        end
    end
    return Af
end

convert(::Type{Matrix}, A::AbstractBlockHankel) = Matrix(A)


text:
using Oceananigans.Operators: index_and_interp_dependencies
using Oceananigans.Utils: tupleit, user_function_arguments

import Oceananigans: location

"""
    ContinuousBoundaryFunction{X, Y, Z, I, F, P, D, N, ℑ} <: Function

A wrapper for the user-defined boundary condition function `func` at location
`X, Y, Z`. `I` denotes the boundary-normal index (`I=1` at western boundaries,
`I=grid.Nx` at eastern boundaries, etc). `F, P, D, N, ℑ` are, respectively, the
user-defined function, parameters, field dependencies, indices of the field dependencies
in `model_fields`, and interpolation operators for interpolating `model_fields` to the
location at which the boundary condition is applied.
"""
struct ContinuousBoundaryFunction{X, Y, Z, I, F, P, D, N, ℑ} <: Function
                          func :: F
                    parameters :: P
            field_dependencies :: D
    field_dependencies_indices :: N
     field_dependencies_interp :: ℑ

    """ Returns a location-less wrapper for `func`, `parameters`, and `field_dependencies`."""
    function ContinuousBoundaryFunction(func, parameters, field_dependencies)

        field_dependencies = tupleit(field_dependencies)

        return new{Nothing, Nothing, Nothing, Nothing,
                   typeof(func), typeof(parameters),
                   typeof(field_dependencies), Nothing, Nothing}(func, parameters, field_dependencies, nothing, nothing)
    end

    function ContinuousBoundaryFunction{X, Y, Z, I}(func, parameters, field_dependencies,
                                                    field_dependencies_indices, field_dependencies_interp) where {X, Y, Z, I}
        return new{X, Y, Z, I,
                   typeof(func),
                   typeof(parameters),
                   typeof(field_dependencies),
                   typeof(field_dependencies_indices),
                   typeof(field_dependencies_interp)}(func,
                                                      parameters,
                                                      field_dependencies,
                                                      field_dependencies_indices,
                                                      field_dependencies_interp)
    end
end

location(bc::ContinuousBoundaryFunction{X, Y, Z}) where {X, Y, Z} = X, Y, Z

#####
##### "Regularization" for IncompressibleModel setup
#####

regularize_boundary_condition(bc, X, Y, Z, I, model_field_names) = bc # fallback

"""
    regularize_boundary_condition(bc::BoundaryCondition{C, <:ContinuousBoundaryFunction},
                                  X, Y, Z, I, model_field_names) where C

Regularizes `bc.condition` for location `X, Y, Z`, boundary index `I`, and `model_field_names`,
returning `BoundaryCondition(C, regularized_condition)`.

The regularization of `bc.condition::ContinuousBoundaryFunction` requries

1. Setting the boundary location to `X, Y, Z`.
   The boundary-normal direction is tagged with `Nothing` location.

2. Setting the boundary-normal index `I` for indexing into `field_dependencies`.
   `I` is either `1` (for left boundaries) or
   `size(grid, n)` for a boundary in the `n`th direction where `n ∈ (1, 2, 3)` corresponds
   to `x, y, z`.

3. Determining the `indices` that map `model_fields` to `field_dependencies`.

4. Determining the `interps` functions that interpolate field_dependencies to the location
   of the boundary.
"""
function regularize_boundary_condition(bc::BoundaryCondition{C, <:ContinuousBoundaryFunction},
                                       LX, LY, LZ, I, model_field_names) where C
    boundary_func = bc.condition

    indices, interps = index_and_interp_dependencies(LX, LY, LZ,
                                                     boundary_func.field_dependencies,
                                                     model_field_names)

    regularized_boundary_func = ContinuousBoundaryFunction{LX, LY, LZ, I}(boundary_func.func,
                                                                       boundary_func.parameters,
                                                                       boundary_func.field_dependencies,
                                                                       indices, interps)

    return BoundaryCondition(C, regularized_boundary_func)
end

#####
##### Kernel functions
#####

@inline function (bc::ContinuousBoundaryFunction{Nothing, LY, LZ, i})(j, k, grid, clock, model_fields) where {LY, LZ, i}
    args = user_function_arguments(i, j, k, grid, model_fields, bc.parameters, bc)
    return bc.func(ynode(LY(), j, grid), znode(LZ(), k, grid), clock.time, args...)
end

@inline function (bc::ContinuousBoundaryFunction{LX, Nothing, LZ, j})(i, k, grid, clock, model_fields) where {LX, LZ, j}
    args = user_function_arguments(i, j, k, grid, model_fields, bc.parameters, bc)
    return bc.func(xnode(LX(), i, grid), znode(LZ(), k, grid), clock.time, args...)
end

@inline function (bc::ContinuousBoundaryFunction{LX, LY, Nothing, k})(i, j, grid, clock, model_fields) where {LX, LY, k}
    args = user_function_arguments(i, j, k, grid, model_fields, bc.parameters, bc)
    return bc.func(xnode(LX(), i, grid), ynode(LY(), j, grid), clock.time, args...)
end

# Don't re-convert ContinuousBoundaryFunctions passed to BoundaryCondition constructor
BoundaryCondition(TBC, condition::ContinuousBoundaryFunction) =
    BoundaryCondition{TBC, typeof(condition)}(condition)

Adapt.adapt_structure(to, bc::BoundaryCondition{C, <:ContinuousBoundaryFunction}) where C =
    BoundaryCondition(C, Adapt.adapt(to, bc.condition))

Adapt.adapt_structure(to, bf::ContinuousBoundaryFunction{LX, LY, LZ, I}) where {LX, LY, LZ, I} =
    ContinuousBoundaryFunction{LX, LY, LZ, I}(Adapt.adapt(to, bf.func),
                                           Adapt.adapt(to, bf.parameters),
                                           nothing,
                                           Adapt.adapt(to, bf.field_dependencies_indices),
                                           Adapt.adapt(to, bf.field_dependencies_interp))


text:
export BinaryTraceHeader, th_byte2sample

mutable struct BinaryTraceHeader
    TraceNumWithinLine                 ::Int32
    TraceNumWithinFile                 ::Int32
    FieldRecord                        ::Int32
    TraceNumber                        ::Int32
    EnergySourcePoint                  ::Int32
    CDP                                ::Int32
    CDPTrace                           ::Int32
    TraceIDCode                        ::Int16
    NSummedTraces                      ::Int16
    NStackedTraces                     ::Int16
    DataUse                            ::Int16
    Offset                             ::Int32
    RecGroupElevation                  ::Int32
    SourceSurfaceElevation             ::Int32
    SourceDepth                        ::Int32
    RecDatumElevation                  ::Int32
    SourceDatumElevation               ::Int32
    SourceWaterDepth                   ::Int32
    GroupWaterDepth                    ::Int32
    ElevationScalar                    ::Int16
    RecSourceScalar                    ::Int16
    SourceX                            ::Int32
    SourceY                            ::Int32
    GroupX                             ::Int32
    GroupY                             ::Int32
    CoordUnits                         ::Int16
    WeatheringVelocity                 ::Int16
    SubWeatheringVelocity              ::Int16
    UpholeTimeSource                   ::Int16
    UpholeTimeGroup                    ::Int16
    StaticCorrectionSource             ::Int16
    StaticCorrectionGroup              ::Int16
    TotalStaticApplied                 ::Int16
    LagTimeA                           ::Int16
    LagTimeB                           ::Int16
    DelayRecordingTime                 ::Int16
    MuteTimeStart                      ::Int16
    MuteTimeEnd                        ::Int16
    ns                                 ::Int16
    dt                                 ::Int16
    GainType                           ::Int16
    InstrumentGainConstant             ::Int16
    InstrumntInitialGain               ::Int16
    Correlated                         ::Int16
    SweepFrequencyStart                ::Int16
    SweepFrequencyEnd                  ::Int16
    SweepLength                        ::Int16
    SweepType                          ::Int16
    SweepTraceTaperLengthStart         ::Int16
    SweepTraceTaperLengthEnd           ::Int16
    TaperType                          ::Int16
    AliasFilterFrequency               ::Int16
    AliasFilterSlope                   ::Int16
    NotchFilterFrequency               ::Int16
    NotchFilterSlope                   ::Int16
    LowCutFrequency                    ::Int16
    HighCutFrequency                   ::Int16
    LowCutSlope                        ::Int16
    HighCutSlope                       ::Int16
    Year                               ::Int16
    DayOfYear                          ::Int16
    HourOfDay                          ::Int16
    MinuteOfHour                       ::Int16
    SecondOfMinute                     ::Int16
    TimeCode                           ::Int16
    TraceWeightingFactor               ::Int16
    GeophoneGroupNumberRoll            ::Int16
    GeophoneGroupNumberTraceStart      ::Int16
    GeophoneGroupNumberTraceEnd        ::Int16
    GapSize                            ::Int16
    OverTravel                         ::Int16
    CDPX                               ::Int32
    CDPY                               ::Int32
    Inline3D                           ::Int32
    Crossline3D                        ::Int32
    ShotPoint                          ::Int32
    ShotPointScalar                    ::Int16
    TraceValueMeasurmentUnit           ::Int16
    TransductionConstnatMantissa       ::Int32
    TransductionConstantPower          ::Int16
    TransductionUnit                   ::Int16
    TraceIdentifier                    ::Int16
    ScalarTraceHeader                  ::Int16
    SourceType                         ::Int16
    SourceEnergyDirectionMantissa      ::Int32
    SourceEnergyDirectionExponent      ::Int16
    SourceMeasurmentMantissa           ::Int32
    SourceMeasurementExponent          ::Int16
    SourceMeasurmentUnit               ::Int16
    Unassigned1                        ::Int32
    Unassigned2                        ::Int32
end

function th_byte2sample()
    Dict{String, Int32}(
        "TraceNumWithinLine"              =>  1  -1,
        "TraceNumWithinFile"              =>  5  -1,
        "FieldRecord"                     =>  9  -1,
        "TraceNumber"                     =>  13 -1,
        "EnergySourcePoint"               =>  17 -1,
        "CDP"                             =>  21 -1,
        "CDPTrace"                        =>  25 -1,
        "TraceIDCode"                     =>  29 -1,
        "NSummedTraces"                   =>  31 -1,
        "NStackedTraces"                  =>  33 -1,
        "DataUse"                         =>  35 -1,
        "Offset"                          =>  37 -1,
        "RecGroupElevation"               =>  41 -1,
        "SourceSurfaceElevation"          =>  45 -1,
        "SourceDepth"                     =>  49 -1,
        "RecDatumElevation"               =>  53 -1,
        "SourceDatumElevation"            =>  57 -1,
        "SourceWaterDepth"                =>  61 -1,
        "GroupWaterDepth"                 =>  65 -1,
        "ElevationScalar"                 =>  69 -1,
        "RecSourceScalar"                 =>  71 -1,
        "SourceX"                         =>  73 -1,
        "SourceY"                         =>  77 -1,
        "GroupX"                          =>  81 -1,
        "GroupY"                          =>  85 -1,
        "CoordUnits"                      =>  89 -1,
        "WeatheringVelocity"              =>  91 -1,
        "SubWeatheringVelocity"           =>  93 -1,
        "UpholeTimeSource"                =>  95 -1,
        "UpholeTimeGroup"                 =>  97 -1,
        "StaticCorrectionSource"          =>  99 -1,
        "StaticCorrectionGroup"           =>  101-1,
        "TotalStaticApplied"              =>  103-1,
        "LagTimeA"                        =>  105-1,
        "LagTimeB"                        =>  107-1,
        "DelayRecordingTime"              =>  109-1,
        "MuteTimeStart"                   =>  111-1,
        "MuteTimeEnd"                     =>  113-1,
        "ns"                              =>  115-1,
        "dt"                              =>  117-1,
        "GainType"                        =>  119-1,
        "InstrumentGainConstant"          =>  121-1,
        "InstrumntInitialGain"            =>  123-1,
        "Correlated"                      =>  125-1,
        "SweepFrequencyStart"             =>  127-1,
        "SweepFrequencyEnd"               =>  129-1,
        "SweepLength"                     =>  131-1,
        "SweepType"                       =>  133-1,
        "SweepTraceTaperLengthStart"      =>  135-1,
        "SweepTraceTaperLengthEnd"        =>  137-1,
        "TaperType"                       =>  139-1,
        "AliasFilterFrequency"            =>  141-1,
        "AliasFilterSlope"                =>  143-1,
        "NotchFilterFrequency"            =>  145-1,
        "NotchFilterSlope"                =>  147-1,
        "LowCutFrequency"                 =>  149-1,
        "HighCutFrequency"                =>  151-1,
        "LowCutSlope"                     =>  153-1,
        "HighCutSlope"                    =>  155-1,
        "Year"                            =>  157-1,
        "DayOfYear"                       =>  159-1,
        "HourOfDay"                       =>  161-1,
        "MinuteOfHour"                    =>  163-1,
        "SecondOfMinute"                  =>  165-1,
        "TimeCode"                        =>  167-1,
        "TraceWeightingFactor"            =>  169-1,
        "GeophoneGroupNumberRoll"         =>  171-1,
        "GeophoneGroupNumberTraceStart"   =>  173-1,
        "GeophoneGroupNumberTraceEnd"     =>  175-1,
        "GapSize"                         =>  177-1,
        "OverTravel"                      =>  179-1,
        "CDPX"                            =>  181-1,
        "CDPY"                            =>  185-1,
        "Inline3D"                        =>  189-1,
        "Crossline3D"                     =>  193-1,
        "ShotPoint"                       =>  197-1,
        "ShotPointScalar"                 =>  201-1,
        "TraceValueMeasurmentUnit"        =>  203-1,
        "TransductionConstnatMantissa"    =>  205-1,
        "TransductionConstantPower"       =>  209-1,
        "TransductionUnit"                =>  211-1,
        "TraceIdentifier"                 =>  213-1,
        "ScalarTraceHeader"               =>  215-1,
        "SourceType"                      =>  217-1,
        "SourceEnergyDirectionMantissa"   =>  219-1,
        "SourceEnergyDirectionExponent"   =>  223-1,
        "SourceMeasurmentMantissa"        =>  225-1,
        "SourceMeasurementExponent"       =>  229-1,
        "SourceMeasurmentUnit"            =>  231-1,
        "Unassigned1"                     =>  233-1,
        "Unassigned2"                     =>  237-1)

end

function BinaryTraceHeader()
    BinaryTraceHeader(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,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,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 show(io::IO, bth::BinaryTraceHeader)

    println("BinaryTraceHeader:")
    for field in fieldnames(BinaryTraceHeader)
        s = @sprintf "     %30s: %9d" String(field) getfield(bth, field)
        println(s)
    end
    println("\n")

end

function show(io::IO, Abth::Array{BinaryTraceHeader,1})

    if length(Abth) == 0
        println("Empty Traceheaders")
    else
        # Show first
        println("BinaryTraceHeader 1:")
        for field in fieldnames(BinaryTraceHeader)
            s = @sprintf "     %30s: %9d" String(field) getfield(Abth[1], field)
            println(s)
        end

        println("\n ... \n")

        # Show last
        println("BinaryTraceHeader $(length(Abth)):")
        for field in fieldnames(BinaryTraceHeader)
            s = @sprintf "     %30s: %9d" String(field) getfield(Abth[end], field)
            println(s)
        end
    end
end


text:
struct InvalidCallbackReturnValue <: Exception
    msg::String
end

function Base.showerror(io::IO, ex::InvalidCallbackReturnValue)
    Base.write(io, "Invalid callback return value: ", ex.msg)
end


text:
# Runs a simple test of kepler_init.jl

const GNEWT = 39.4845/365.242^2
const third = 1.0/3.0

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

@testset "kepler_init" begin

ntrial = 10
for itrial=1:ntrial
t0 = 2.4
mass = 1.0
period = 1.5
elements = ones(Float64,6)
while elements[3]^2+elements[4]^2 >= 0.2^2
  elements = [period,rand()*period,randn(),randn(),rand()*pi,rand()*pi]
end
elements_diff = zeros(Float64,6)
ecc = sqrt(elements[3]^2+elements[4]^2)
ntime = 1000
time = linearspace(t0,t0+period,ntime)
xvec = zeros(Float64,3,ntime)
vvec = zeros(Float64,3,ntime)
vfvec = zeros(Float64,3,ntime)
timebig = big.(time)
#dt = 1e-8
dt = period/ntime
dtbig = big(dt)
jac_init = zeros(Float64,7,7)
# Compute Jacobian in BigFloat precision with analytic formulae:
jac_init_big = zeros(BigFloat,7,7)
# Compute Jacobian in BigFloat precision with finite differences:
jac_init_num = zeros(BigFloat,7,7)
# Variations in the elements for finite differences:
delements = big.([1e-15,1e-15,1e-15,1e-15,1e-15,1e-15,1e-15])
elements_name = ["P","t0","ecosom","esinom","inc","Omega","Mass"]
cartesian_name = ["x","y","z","vx","vy","vz","mass"]
println("Elements: ",elements)
massbig=big(mass)
for i=1:ntime
  x,v = kepler_init(timebig[i],massbig,big.(elements))
  if i==1
    x_ekep,v_ekep = kepler_init(time[i],mass,elements,jac_init)
    # Redo in BigFloat
    elements_big = big.(elements)
    x_ekep_big,v_ekep_big = kepler_init(timebig[i],massbig,elements_big,jac_init_big)
  # Check that these agree:
#    println("x-x_ekep: ",maximum(abs.(x-x_ekep))," v-v_ekep: ",maximum(abs.(v-v_ekep)))
  # Now take some derivatives:
    for j=1:7
      elements_diff = big.(elements)
      if j <= 6
        elements_diff[j] -= delements[j]
      else
        massbig -= delements[7]
      end
      x_minus,v_minus= kepler_init(timebig[i],massbig,elements_diff)
      elements_diff .= elements
      if j <= 6
        elements_diff[j] += delements[j]
      else
        massbig += 2.0*delements[7]
      end
      x_plus,v_plus= kepler_init(timebig[i],massbig,elements_diff)
      if j == 7
        massbig -= delements[7]
      end
      for k=1:3
        jac_init_num[  k,j] = .5*(x_plus[k]-x_minus[k])/delements[j]
        jac_init_num[3+k,j] = .5*(v_plus[k]-v_minus[k])/delements[j]
      end
    end
    jac_init_num[7,7]=1.0
    for j=1:7, k=1:7
      if abs(jac_init[k,j]-jac_init_num[k,j]) > 1e-8
        println(elements_name[j]," ",cartesian_name[k]," ",jac_init[k,j]," ",jac_init_num[k,j]," ",jac_init[k,j]-jac_init_num[k,j])
      end
      if abs(jac_init[k,j]/jac_init_num[k,j]-1) > 1e-12 && jac_init_num[k,j] != 0.0
        println(elements_name[j]," ",cartesian_name[k]," ",jac_init[k,j]," ",jac_init_num[k,j]," ",jac_init[k,j]/jac_init_num[k,j]-1)
      end
    end
    println("Jacobians: difference ",maximum(abs.(jac_init-jac_init_num)))
    println("Jacobians: log diff   ",maximum(abs.(jac_init[jac_init_num .!= 0.0]./jac_init_num[jac_init_num .!= 0.0].-1)))
  end
  xvec[:,i]=x
  vvec[:,i]=v
  # Compute finite difference velocity:
  xf,vf = kepler_init(timebig[i]+dtbig,massbig,big.(elements))
  vfvec[:,i]=(xf-x)/dt
end
period = elements[1]
#omega=atan2(elements[4],elements[3])
omega=atan(elements[4],elements[3])
semi = (GNEWT*mass*period^2/4/pi^2)^third
# Compute the focus:
focus = [semi*ecc,0.,0.]
xfocus = semi*ecc*(-cos(omega)-sin(omega))
zfocus = semi*ecc*sin(omega)
# Check that we have an ellipse (we're assuming that the motion is in the x-z plane; no longer true):
Atot = 0.0
dAdt = zeros(Float64,ntime)
dx = xvec[1,1]-xvec[1,ntime-1]
dy = xvec[2,1]-xvec[2,ntime-1]
dz = xvec[3,1]-xvec[3,ntime-1]
dAvec = 0.5*[xvec[2,1]*dz-xvec[3,1]*dy;xvec[3,1]*dx-xvec[1,1]*dz;xvec[1,1]*dy-xvec[2,1]*dx]
#println("dA: ",norm(dAvec)," ",dAvec/norm(dAvec))
#dA = 0.5*(xvec[3,1]*dx-xvec[1,1]*dz)
dAdt[1] = norm(dAvec)/(time[2]-time[1])
Atot += norm(dAvec)
for i=2:ntime
  dx = xvec[1,i]-xvec[1,i-1]
  dy = xvec[2,i]-xvec[2,i-1]
  dz = xvec[3,i]-xvec[3,i-1]
  dAvec = 0.5*[xvec[2,i]*dz-xvec[3,i]*dy;xvec[3,i]*dx-xvec[1,i]*dz;xvec[1,i]*dy-xvec[2,i]*dx]
#  println("dA: ",norm(dAvec)," ",dAvec/norm(dAvec))
  #dA = 0.5*(xvec[3,i]*dx-xvec[1,i]*dz)
  dAdt[i] = norm(dAvec)/(time[i]-time[i-1])
  Atot += norm(dAvec)
end
println("Total area: ",Atot," ",pi*sqrt(1-ecc^2)*semi^2," ratio: ",Atot/pi/semi^2/sqrt(1-ecc^2))
#using PyPlot
#plot(time,dAdt)
#axis([minimum(time),maximum(time),0.,1.5*maximum(dAdt)])
println("Scatter in dAdt: ",std(dAdt))

#plot(time,vvec[1,:])
#plot(time,vfvec[1,:],".")
#plot(time,vvec[1,:]-vfvec[1,:],".")
#plot(time,vvec[3,:])
#plot(time,vfvec[3,:],".")
#plot(time,vvec[3,:]-vfvec[3,:],".")
# Check that velocities match finite difference values

@test isapprox(jac_init,jac_init_num;norm=maxabs)
@test isapprox(jac_init,convert(Array{Float64,2},jac_init_big);norm=maxabs)
end
end


text:
using DifferentialEquations
using Test
using MRIgeneralizedBloch
using MRIgeneralizedBloch: hamiltonian_linear, propagator_linear_inversion_pulse

##
max_error = 1e-5

## set parameters
α  = π
ω1 = α / 500e-6
ω0 = 200.0
B1 = 1.0
m0s = 0.15
R1f = 0.3
R2f = 1 / 65e-3
R1s = 2
T2s = 10e-6
Rx = 30.0
TRF = 500e-6
N = Inf

## precomputations
(R2sl, dR2sldT2s, dR2sldB1) = precompute_R2sl()

_R2sl = R2sl(TRF, α, B1, T2s)
_dR2sldT2s = dR2sldT2s(TRF, α, B1, T2s)
_dR2sldB1 = dR2sldB1(TRF, α, B1, T2s)

## baseline IDE solution
m0  = [0.5 * (1 - m0s), 0, 0.5 * (1 - m0s), 0, m0s, 1]
m0g = [m0[1:5]; zeros(5); 1]
m = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, _R2sl)) * m0

## derivative wrt. m0s
dm0s = 1e-6
mg = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, _R2sl, _dR2sldT2s, _dR2sldB1, grad_m0s())) * m0g

md = exp(hamiltonian_linear(ω1, B1, ω0, TRF, (m0s + dm0s), R1f, R2f, Rx, R1s, _R2sl)) * m0

mfd = (md - m) / dm0s
@test mg[6:10] ≈ mfd[1:5] rtol = max_error

## test derivative wrt. R1f
dR1f = 1e-6
mg = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, _R2sl, _dR2sldT2s, _dR2sldB1, grad_R1f())) * m0g

md = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, (R1f + dR1f), R2f, Rx, R1s, _R2sl)) * m0

mfd = (md - m) / dR1f
@test mg[6:10] ≈ mfd[1:5] rtol = max_error

## test derivative wrt. R2f
dR2f = 1e-6
mg = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, _R2sl, _dR2sldT2s, _dR2sldB1, grad_R2f())) * m0g

md = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, R1f, (R2f + dR2f), Rx, R1s, _R2sl)) * m0

mfd = (md - m) / dR2f
@test mg[6:10] ≈ mfd[1:5] rtol = max_error

## test derivative wrt. Rx
dRx = 1e-5
mg = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, _R2sl, _dR2sldT2s, _dR2sldB1, grad_Rx())) * m0g

md = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, R1f, R2f, (Rx + dRx), R1s, _R2sl)) * m0

mfd = (md - m) / dRx
@test mg[6:10] ≈ mfd[1:5] rtol = max_error

## test derivative wrt. R1s
dR1s = 1e-6
mg = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, _R2sl, _dR2sldT2s, _dR2sldB1, grad_R1s())) * m0g

md = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, (R1s + dR1s), _R2sl)) * m0

mfd = (md - m) / dR1s
@test mg[6:10] ≈ mfd[1:5] rtol = max_error

## test derivative wrt. T2s
dT2s = 1e-14
mg = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, _R2sl, _dR2sldT2s, _dR2sldB1, grad_T2s())) * m0g

_dR2sl = R2sl(TRF, α, B1, (T2s+dT2s))
md = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, _dR2sl)) * m0

mfd = (md - m) / dT2s
@test mg[6:10] ≈ mfd[1:5] rtol = max_error

## test derivative wrt. ω0
dω0 = 1e-6
mg = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, _R2sl, _dR2sldT2s, _dR2sldB1, grad_ω0())) * m0g

md = exp(hamiltonian_linear(ω1, B1, (ω0 + dω0), TRF, m0s, R1f, R2f, Rx, R1s, _R2sl)) * m0

mfd = (md - m) / dω0
@test mg[6:10] ≈ mfd[1:5] rtol = max_error

## test derivative wrt. B1
dB1 = 1e-6
mg = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, _R2sl, _dR2sldT2s, _dR2sldB1, grad_B1())) * m0g

_dR2sl = R2sl(TRF, α, (B1 + dB1), T2s)
md = exp(hamiltonian_linear(ω1, (B1 + dB1), ω0, TRF, m0s, R1f, R2f, Rx, R1s, _dR2sl)) * m0

mfd = (md - m) / dB1
@test mg[6:10] ≈ mfd[1:5] rtol = max_error

## test derivative wrt. R1a
R1a = 1
dR1a = 1e-6
m = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, R1a, R2f, Rx, R1a, _R2sl)) * m0
mg = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, R1a, R2f, Rx, R1a, _R2sl, _dR2sldT2s, _dR2sldB1, grad_R1a())) * m0g

md = exp(hamiltonian_linear(ω1, B1, ω0, TRF, m0s, (R1a + dR1a), R2f, Rx, (R1a + dR1a), _R2sl)) * m0

mfd = (md - m) / dR1a
@test mg[6:10] ≈ mfd[1:5] rtol = max_error

## ########################################################################
# inversion pulse
###########################################################################
m = propagator_linear_inversion_pulse(ω1, TRF, B1, _R2sl, undef, undef, undef) * m0

## derivative wrt. m0s
dm0s = 1e-6
mg = propagator_linear_inversion_pulse(ω1, TRF, B1, _R2sl, _dR2sldT2s, _dR2sldB1, grad_m0s()) * m0g

md = propagator_linear_inversion_pulse(ω1, TRF, B1, _R2sl, undef, undef, undef) * m0

mfd = (md - m) / dm0s
@test mg[6:10] ≈ mfd[1:5] rtol = max_error

## test derivative wrt. T2s
dT2s = 1e-14
mg = propagator_linear_inversion_pulse(ω1, TRF, B1, _R2sl, _dR2sldT2s, _dR2sldB1, grad_T2s()) * m0g

_dR2sl = R2sl(TRF, α, B1, (T2s+dT2s))
md = propagator_linear_inversion_pulse(ω1, TRF, B1, _dR2sl, undef, undef, undef) * m0

mfd = (md - m) / dT2s
@test mg[6:10] ≈ mfd[1:5] rtol = max_error

## test derivative wrt. B1
dB1 = 1e-9
mg = propagator_linear_inversion_pulse(ω1, TRF, B1, _R2sl, _dR2sldT2s, _dR2sldB1, grad_B1()) * m0g

_dR2sl = R2sl(TRF, α, (B1 + dB1), T2s)
md = propagator_linear_inversion_pulse(ω1, TRF, (B1 + dB1), _dR2sl, undef, undef, undef) * m0

mfd = (md - m) / dB1
@test mg[6:10] ≈ mfd[1:5] rtol = max_error


text:
########################################################################
# Atoms of the block nuclear-norm ball.
########################################################################
"""
    BlkNucBallAtom(m,n,cidx,u,v)

Atom in the BlkNucBall with index cidx.
"""
mutable struct BlkNucBallAtom{T1<:Int64, T2<:CartesianIndex, T3<:NucBallAtom{Int64, Vector{Float64}}} <: AbstractAtom
    m::T1
    n::T1
    cidx::T2
    ba::T3
    function BlkNucBallAtom(m::Int64, n::Int64, cidx::CartesianIndex, ba::NucBallAtom{Int64, Vector{Float64}})
        new{Int64, CartesianIndex, NucBallAtom{Int64, Vector{Float64}}}(m, n, cidx, ba)
    end
end

"""
    *(M::AbstractLinearOp, a::NucBallAtom) -> Y

Multiply a nuclear-norm atom by a linear map, and return the matrix `Y`.
"""
function Base.:(*)(M::AbstractLinearOp, a::BlkNucBallAtom)
    m, n = a.m, a.n
    c1, c2 = a.cidx[1], a.cidx[2]
    u, v = a.ba.u, a.ba.v
    bm = length(u); bn = length(v)
    A = spzeros(m*n)
    @inbounds for bj in 1:bn
        @inbounds for bi in 1:bm
            i = (c1-1)*bm + bi
            j = (c2-1)*bn + bj
            idx = (j-1)*m + i
            val = u[bi] * v[bj]
            push!(A.nzind, idx)
            push!(A.nzval, val)
        end
    end
    return M*A
end
function LinearAlgebra.mul!(Ma::Vector{Float64}, M::LinearOp, a::BlkNucBallAtom)
    m, n = a.m, a.n
    c1, c2 = a.cidx[1], a.cidx[2]
    u, v = a.ba.u, a.ba.v
    bm = length(u); bn = length(v)
    A = spzeros(m*n)
    @inbounds for bj in 1:bn
        @inbounds for bi in 1:bm
            i = (c1-1)*bm + bi
            j = (c2-1)*bn + bj
            idx = (j-1)*m + i
            val = u[bi] * v[bj]
            push!(A.nzind, idx)
            push!(A.nzval, val)
        end
    end
    mul!(Ma, M, A)
    return nothing
end
Base.vec(a::BlkNucBallAtom) = I(a.m*a.n)*a
Base.length(a::BlkNucBallAtom) = a.m * a.n

########################################################################
# Block Nuclear-norm ball.
########################################################################
"""
    BlkNucBall(m, n, bm, bn, maxrank)

"""
struct BlkNucBall{T1<:Int64, T2<:Matrix{Int64}} <: AbstractAtomicSet
    m::T1           # matrix size 1
    n::T1           # matrix size 2
    bm::T1          # block size 1
    bn::T1          # block size 2
    l::T1           # num of blocks 1
    h::T1           # num of blocks 2
    maxrank::T2     # max rank for each block
    function BlkNucBall(n::Int64, m::Int64, bn::Int64, bm::Int64, maxrank::Matrix{Int64})
        l = floor(Int, m/bm)
        h = floor(Int, n/bn)
        new{Int64, Matrix{Int64}}(m, n, bm, bn, l, h, maxrank)
    end
end

"""
    gauge(A::BlkNucBall, x::Matrix)

Gives the block nuclear-norm of `x`.
"""
function gauge(A::BlkNucBall, x::Matrix{Float64})
    bm, bn, l, h = A.bm, A.bn, A.l, A.h
    Ab = NucBall(bm, bn)
    f(i::CartesianIndex) = gauge(Ab, x[(i[1]-1)*bm+1:i[1]*bm, (i[2]-1)*bn+1:i[2]*bn])
    return mapreduce(f, +, CartesianIndices((l, h)))
end
gauge(A::BlkNucBall, x::Vector{Float64}) = gauge(A, reshape(x, A.m, A.n))

"""
    support(A::BlkNucBall, z::Matirx)

Gives the largest singular value of `z`.
"""
function support(A::BlkNucBall, z::Matrix{Float64})
    bm, bn, l, h = A.bm, A.bn, A.l, A.h
    Ab = NucBall(bm, bn)
    f(i::CartesianIndex) = support(Ab, z[(i[1]-1)*bm+1:i[1]*bm, (i[2]-1)*bn+1:i[2]*bn])
    return maximum(pmap(f, CartesianIndices((l, h))))
end
support(A::BlkNucBall, z::Vector{Float64}) = support(A, reshape(z, A.m, A.n))

"""
    expose(A::BlkNucBall, z::Matrix)

"""
function expose!(A::BlkNucBall, z::Matrix{Float64}, a::BlkNucBallAtom)
    bm, bn, l, h = A.bm, A.bn, A.l, A.h
    Ab = NucBall(bm, bn)
    f(i::CartesianIndex) = support(Ab, z[(i[1]-1)*bm+1:i[1]*bm, (i[2]-1)*bn+1:i[2]*bn])
    ci = argmax(pmap(f, CartesianIndices((l, h))))
    a.cidx = ci
    expose!(Ab, z[(ci[1]-1)*bm+1:ci[1]*bm, (ci[2]-1)*bn+1:ci[2]*bn], a.ba)
    return nothing
end
expose!(A::BlkNucBall, z::Vector{Float64}, a::BlkNucBallAtom) = expose!(A, reshape(z, A.m, A.n), a)

function expose(A::BlkNucBall, z::Matrix{Float64})
    bm, bn = A.bm, A.bn
    Ab = NucBall(bm, bn)
    ba = expose(Ab, zeros(bm, bn))
    a = BlkNucBallAtom(A.m, A.n, CartesianIndex(1, 1), ba)
    expose!(A, z, a)
    return a
end
expose(A::BlkNucBall, z::Vector{Float64}) = expose(A, reshape(z, A.m, A.n))

Base.length(A::BlkNucBall) = A.m*A.n
rank(A::BlkNucBall) = A.maxrank
atom_name(A::BlkNucBall) = "block nuclear-norm ball"
atom_description(A::BlkNucBall) = "{ x ∈ ℝⁿ | ||x||₁ ≤ 1 }"
atom_parameters(A::BlkNucBall) = "n×m = $(length(A)); maxrank = $(A.maxrank)"

########################################################################
# Face of the nuclear-norm ball.
########################################################################
mutable struct BlkNucBallFace{T1<:Int64, T2<:Matrix{NucBallFace{Int64, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}}, T3<:SparseMatrixCSC{Float64, Int64}} <: AbstractFace
    m::T1
    n::T1
    nucBallFaces::T2
    S::T3
    function BlkNucBallFace(m::Int64, n::Int64, nucBallFaces::Matrix{NucBallFace{Int64, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}})
        l, h = size(nucBallFaces)
        bm, bn = size(nucBallFaces[1,1])
        I = Vector{Int64}(undef, m*n)
        for k in 1:m*n
            I[k] = indexmatch(k, m, bm, l, bn)
        end
        S = sparse(I, collect(1:m*n), ones(m*n), m*n, m*n)
        new{Int64, Matrix{NucBallFace{Int64, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}},SparseMatrixCSC{Float64, Int64}}(m, n, nucBallFaces, S)
    end
end

function indexmatch(k::Int64, m::Int64, bm::Int64, l::Int64, bn::Int64)
    u, v = divrem(k, bm*bn)
    if v == 0
        u = u - 1; v = bm*bn
    end
    c, d = divrem(u+1, l)
    if d == 0
        c = c - 1; d = l
    end
    e, f = divrem(v, bm)
    if f == 0
        e = e - 1; f = bm
    end
    return (c*bn + e)*m + (d-1)*bm + f
end

"""
    face(A, z) -> F

"""
function face!(A::BlkNucBall, z::Matrix{Float64}, F::BlkNucBallFace)
    bm, bn, l, h, maxrank = A.bm, A.bn, A.l, A.h, A.maxrank
    f(i::CartesianIndex) = face!( NucBall(bn, bm, maxrank=maxrank[i]), z[(i[1]-1)*bm+1:i[1]*bm, (i[2]-1)*bn+1:i[2]*bn], F.nucBallFaces[i])
    pmap(f, CartesianIndices((l, h)))
    return nothing
end
face!(A::BlkNucBall, z::Vector{Float64}, F::BlkNucBallFace) = face!(A, reshape(z, A.m, A.n), F)

function face(A::BlkNucBall, z::Matrix{Float64})
    bm, bn, l, h, maxrank = A.bm, A.bn, A.l, A.h, A.maxrank
    f(i::CartesianIndex) = face( NucBall(bn, bm, maxrank=maxrank[i]), z[(i[1]-1)*bm+1:i[1]*bm, (i[2]-1)*bn+1:i[2]*bn])
    Fs = pmap(f, CartesianIndices((l, h)))
    return BlkNucBallFace(A.m, A.n, Fs)
end
face(A::BlkNucBall, z::Vector{Float64}) = face(A, reshape(z, A.m, A.n))

"""
    *(M::AbstractLinearOp, F::NucBallFace)
"""
function Base.:(*)(M::AbstractLinearOp, F::BlkNucBallFace)
    return M*F.S*Base.cat(map(vec, F.nucBallFaces)...; dims=(1,2))
end
Base.:(*)(λ::Real, F::BlkNucBallFace) = λ*vec(F)
Base.length(F::BlkNucBallFace) = F.n*F.m
rank(F::BlkNucBallFace) = sum(map(rank, F.nucBallFaces))
vec(F::BlkNucBallFace) = I(F.m*F.n)*F
face_name(F::BlkNucBallFace) = "block nuclear-norm ball"
face_parameters(F::BlkNucBallFace) = "rank = $(rank(F)); n = $(length(F))"

"""
Given a face and a set of weights, reveal the corresponding
point on the face.
"""
Base.:(*)(F::BlkNucBallFace, c::Vector) = vec(F)*c


text:
using JSON
using Sockets

function parse_config(filename::String)
    json_data = Dict()
    open(filename, "r") do f
        json_data = JSON.parse(f)
    end
    return parse_config(json_data)
end
function parse_config(json_data::Dict{String, N}) where N<:Any
    if !("strips" in keys(json_data) && "controllers" in keys(json_data))
        error("Config file was invalid")
    end
    known_channels = Array{Int,1}(undef, 0)
    known_controllers = Array{String,1}(undef, 0)
    for strip in json_data["strips"]
        if !(strip["channel"] in known_channels)
            push!(known_channels, strip["channel"])
            sort!(known_channels)
        end
        if !(strip["controller"] in known_controllers)
            push!(known_controllers, strip["controller"])
            sort!(known_controllers)
        end
    end
    num_channels = length(known_channels)
    num_controllers = length(known_controllers)
    channel_to_idx = Dict(num => findfirst(known_channels,num) for num in known_channels)
    controller_to_idx = Dict(num => findfirst(known_controllers,num) for num in known_controllers)
    channels = Array{LEDChannel, 1}(undef, num_channels)
    for i in eachindex(channels)
        channels[i] = LEDChannel()
    end
    controllers = Array{LEDController, 1}(undef, num_controllers)
    for i in eachindex(controllers)
        controller_json = json_data["controllers"][known_controllers[i]]
        controllers[i] = LEDController(controller_json["size"], (IPv4(controller_json["ip"]), controller_json["port"]))
    end
    strips = Array{LEDStrip, 1}(undef, 0)
    for strip in json_data["strips"]
        tmp_channel = channels[channel_to_idx[strip["channel"]]]
        tmp_controller = controllers[controller_to_idx[strip["controller"]]]
        push!(strips, LEDStrip(strip["name"], tmp_channel, tmp_controller, strip["start"], strip["end"]))
    end
    led_array = LEDArray(controllers, Array{LEDChannel,1}(undef, 0), channels, strips)
    return led_array
end

function remote_config(outgoing_socket, main_port, discovery_port, subnet)
    broadcast_packet = getBroadcastPacket()
    discovery_socket = UDPSocket()
    bind(discovery_socket, ip"0.0.0.0", discovery_port)
    send(outgoing_socket, subnet, main_port, broadcast_packet)
    discovered_clients = Array{Tuple{IPAddr,Dict},1}(undef, 0)
    searching = true
    @async begin
        while searching
            try
                packet = recvfrom(discovery_socket)
                json_data = Dict{String,Any}()
                try
                    json_data = JSON.parse(convert(String, packet[2]))
                catch
                    println("Failed to parse json")
                end
                if validate_json(json_data)
                    push!(discovered_clients, (packet[1],json_data))
                end
            catch e
                if !isa(e, EOFError)
                    rethrow(e)
                end
            end
        end
    end

    sleep(3)
    searching = false
    close(discovery_socket)

    output = get_overall_config_template()
    known_channels = Array{Int}(undef, 0)
    for data_tup in discovered_clients
        address, data = data_tup
        try
            output["controllers"][data["name"]] = Dict(
                "ip"=>address,
                "port"=>data["port"],
                "size"=>data["numAddrs"]
            )
            for strip in data["strips"]
                push!(output["strips"], Dict(
                    "name"=>strip["name"],
                    "start"=>strip["startAddr"],
                    "end"=>strip["endAddr"],
                    "channel"=>strip["channel"],
                    "controller"=>data["name"]
                ))
                if !(strip["channel"] in known_channels)
                    push!(known_channels, strip["channel"])
                end
            end
        catch
            println("Got some trash json")
        end
        output["general"]["numChannels"] = length(known_channels)
    end
    if length(output["controllers"]) > 0
        return parse_config(output::Dict)
    else
        error("Failed to find any clients")
    end
end

function getBroadcastPacket()
    #macString = readstring(`cat /sys/class/net/eth0/address`)
    ipString = split(read(`hostname -I`, String))[1]
    output_dict = Dict(
        "ip"=>ipString,
        "mac"=>"mac address here",
        "msg_type"=>"startup"
    )
    return json(output_dict)
end

function get_overall_config_template()
    return Dict(
        "general"=>Dict(
            "numChannels"=>0
        ),
        "controllers"=>Dict{String, Dict{String, Union{Number, String, IPAddr}}}(),
        "strips"=>Array{Dict{String,Union{Number, String}},1}(undef, 0)
    )
end

function validate_json(json_data)
    ref_dict = Dict(
        "name"=>"",
        "ip"=>"",
        "port"=>0,
        "mac"=>"",
        "numStrips"=>0,
        "numAddrs"=>0,
        "strips"=>[
            Dict(
                "name"=>"",
                "startAddr"=>0,
                "endAddr"=>1,
                "channel"=>0
            )
        ]
    )

    return check_json(json_data, ref_dict)
end

check_json(x::T, y::T) where {T<:Dict{String, N} where N<:Any} = check_symmetry(x,y)
check_json(x::T, y::N) where {T,N} = false

function check_symmetry(x::T, y::T) where {T<:Dict{S, N} where {S<:Any, N<:Any}}
    if collect(keys(x)) == collect(keys(y))
        return reduce(check_symmetry.(collect(values(x)), collect(values(y)))) do x, y
            return x && y
        end
    else
        return false
    end
end

function check_symmetry(x::AbstractArray, y::AbstractArray)
    try
        return reduce(check_symmetry.(x,y)) do x, y
            return x && y
        end
    catch
        return false
    end
end

check_symmetry(x::T, y::T) where T<:Union{Number, String, Bool} = true
check_symmetry(x::T, y::N) where {T, N} = false


text:
# Compression
# #+HTML: Source
# #+HTML:

# [[file:~/.julia/dev/MatrixProductStates/README.org::*Compression][Compression:1]]
function compress(ψ::MPS{L, T}, to_the::Right; Dcut::Int=typemax(Int)) where {L, T}
    tensors = Array{T, 3}[]

    B = ψ[1]
    d = length(B[1, 1, :])

    @cast Bm[(σ¹, a⁰), a¹] |= B[a⁰, a¹, σ¹]
    U, S, V = psvd(Bm, rank=Dcut)
    #S = S/√sum(S .^ 2)

    @cast A[a⁰, a¹, σ¹] |= U[(σ¹, a⁰), a¹] (σ¹:d)
    push!(tensors, A)

    for i ∈ 2:L
        B = ψ[i]
        d = length(B[1, 1, :])

        @tensor M[aⁱ⁻¹, aⁱ, σⁱ] := (Diagonal(S)*V')[aⁱ⁻¹, aⁱ⁻¹′] * B[aⁱ⁻¹′, aⁱ, σⁱ]
        @cast   Mm[(σⁱ, aⁱ⁻¹), aⁱ] |= M[aⁱ⁻¹, aⁱ, σⁱ]

        U, S, V = psvd(Mm, rank=Dcut)
        #S = S/√sum(S .^ 2)

        @cast A[aⁱ⁻¹, aⁱ, σⁱ] |= U[(σⁱ, aⁱ⁻¹), aⁱ] (σⁱ:d)
        push!(tensors, A)
    end
    MPS{L, T}(tensors), Left()
end

leftcanonical(ψ) = compress(ψ, right)[1]

function compress(ψ::MPS{L, T}, to_the::Left; Dcut::Int=typemax(Int)) where {L, T}
    tensors = Array{T, 3}[]

    A = ψ[L]
    d = length(A[1, 1, :])
    @cast Am[aᴸ⁻¹, (σᴸ, aᴸ)] |= A[aᴸ⁻¹, aᴸ, σᴸ]

    U, S, V = psvd(Am, rank=Dcut)
    #S = S/√sum(S .^ 2)

    @cast B[aᴸ⁻¹, aᴸ, σᴸ] |= V'[aᴸ⁻¹, (σᴸ, aᴸ)] (σᴸ:d)
    push!(tensors, B)

    for i ∈ (L-1):-1:1
        A = ψ[i]
        d = length(A[1, 1, :])
        @tensor M[aⁱ⁻¹, aⁱ, σⁱ]    := A[aⁱ⁻¹, aⁱ′, σⁱ] * (U * Diagonal(S))[aⁱ′, aⁱ]
        @cast   Mm[aⁱ⁻¹, (σⁱ, aⁱ)] |= M[aⁱ⁻¹, aⁱ, σⁱ]

        U, S, V = psvd(Mm, rank=Dcut)
        #S = S/√sum(S .^ 2)

        @cast B[aⁱ⁻¹, aⁱ, σⁱ] |= V'[aⁱ⁻¹, (σⁱ, aⁱ)] (σⁱ:d)
        push!(tensors, B)
    end
    MPS{L, T}(reverse(tensors)), Right()
end

rightcanonical(ψ) = compress(ψ, left)[1]

#compress(ψ; Dcut) = compress(ψ, left, Dcut=Dcut)[1]

function variational_compress(ψ::MPS{L, T}, to_the::Left, ϕ::MPS{L, T}) where {L, T}
    btensors = Array{T, 3}[]

    A=ψ[L]
    d = length(A[1, 1, :])
    @cast Am[aᴸ⁻¹, (σᴸ, aᴸ)] |= A[aᴸ⁻¹, aᴸ, σᴸ]
    U, S, V = svd(Am)
    @cast B[aᴸ⁻¹, aᴸ, σᴸ] |= V'[aᴸ⁻¹, (σᴸ, aᴸ)] (σᴸ:d)
    push!(btensors, B)

    for i in L-1:-1:2
        prevM=1
        for j in 1:i-1
            M = ϕ[j]
            M̃ = ψ[j]
            prevM=sum([M̃[:,:,k]'*prevM*M[:,:,k] for k in 1:d])
        end
        Lm=prevM

        prevM=1
        temp=0
        for j in L:-1:i+1
            temp += 1
            M = ϕ[j]
            M̃ = btensors[temp]
            prevM = sum([M[:,:,k]*prevM*M̃[:,:,k]' for k in 1:d])
        end
        Rm=prevM

        M = ϕ[i]
        @tensor M̃[aⁱ⁻¹, aⁱ, σⁱ] := Lm[aⁱ⁻¹,aⁱ⁻¹′] * Rm[aⁱ′,aⁱ] * M[aⁱ⁻¹′, aⁱ′, σⁱ]
        #println("error rate : ",round(1.0-real(tr(sum(M̃[:,:,k]*M̃[:,:,k]' for k in 1:2))),digits=10))

        @cast   Mm[aⁱ⁻¹, (σⁱ, aⁱ)] |= M̃[aⁱ⁻¹, aⁱ, σⁱ]
        U, S, V = svd(Mm)
        @cast B[aⁱ⁻¹, aⁱ, σⁱ] |= V'[aⁱ⁻¹, (σⁱ, aⁱ)] (σⁱ:d)

        push!(btensors, B)
    end

    prevM=1
    temp=0
    for j in L:-1:2
        temp += 1
        M = ϕ[j]
        M̃ = btensors[temp]
        prevM=sum([M[:,:,k]*prevM*M̃[:,:,k]' for k in 1:d])
    end
    Rm=prevM
    M = ϕ[1]
    @tensor M̃[aⁱ⁻¹, aⁱ, σⁱ] := Rm[aⁱ′,aⁱ] * M[aⁱ⁻¹, aⁱ′, σⁱ]
    #println("error rate : ",round(1.0-real(tr(sum(M̃[:,:,k]*M̃[:,:,k]' for k in 1:2))),digits=10))
    @cast   Mm[aⁱ⁻¹, (σⁱ, aⁱ)] |= M̃[aⁱ⁻¹, aⁱ, σⁱ]
    U, S, V = svd(Mm)
    @cast B[aⁱ⁻¹, aⁱ, σⁱ] |= V'[aⁱ⁻¹, (σⁱ, aⁱ)] (σⁱ:d)

    push!(btensors, B)

    MPS{L, T}(reverse(btensors)), Right()
end

function variational_compress(ψ::MPS{L, T}, to_the::Right, ϕ::MPS{L, T}) where {L, T}
    atensors = Array{T, 3}[]

    B=ψ[1]
    d = length(B[1, 1, :])
    @cast   Bm[(σⁱ, aⁱ⁻¹), aⁱ] |= B[aⁱ⁻¹, aⁱ, σⁱ]
    U, S, V = svd(Bm)
    @cast A[aⁱ⁻¹, aⁱ, σⁱ] |= U[(σⁱ, aⁱ⁻¹), aⁱ] (σⁱ:d)
    push!(atensors, A)

    for i in 2:L-1
        prevM=1
        temp=0
        for j in 1:i-1
            temp+=1
            M = ϕ[j]
            M̃ = atensors[temp]
            prevM=sum([M̃[:,:,k]'*prevM*M[:,:,k] for k in 1:d])
        end
        Lm=prevM

        prevM=1
        for j in L:-1:i+1
            M = ϕ[j]
            M̃ = ψ[j]
            prevM=sum([M[:,:,k]*prevM*M̃[:,:,k]' for k in 1:d])
        end
        Rm=prevM

        M = ϕ[i]
        @tensor M̃[aⁱ⁻¹, aⁱ, σⁱ] := Lm[aⁱ⁻¹,aⁱ⁻¹′] * Rm[aⁱ′,aⁱ] * M[aⁱ⁻¹′, aⁱ′, σⁱ]
        #println("error rate : ",round(1.0-real(tr(sum(M̃[:,:,k]*M̃[:,:,k]' for k in 1:2))),digits=10))

        @cast   Mm[(σⁱ, aⁱ⁻¹), aⁱ] |= M̃[aⁱ⁻¹, aⁱ, σⁱ]
        U, S, V = svd(Mm)
        @cast A[aⁱ⁻¹, aⁱ, σⁱ] |= U[(σⁱ, aⁱ⁻¹), aⁱ] (σⁱ:d)

        push!(atensors, A)
    end

    prevM=1
    temp=0
    for j in 1:L-1
        temp+=1
        M = ϕ[j]
        M̃ = atensors[temp]
        prevM = sum([M̃[:,:,k]'*prevM*M[:,:,k] for k in 1:d])
    end
    Lm=prevM
    M = ϕ[L]
    @tensor M̃[aⁱ⁻¹, aⁱ, σⁱ] := Lm[aⁱ⁻¹,aⁱ⁻¹'] * M[aⁱ⁻¹', aⁱ, σⁱ]
    #println("error rate : ",round(1.0-real(tr(sum(M̃[:,:,k]*M̃[:,:,k]' for k in 1:2))),digits=10))
    @cast   Mm[(σⁱ, aⁱ⁻¹), aⁱ] |= M̃[aⁱ⁻¹, aⁱ, σⁱ]
    U, S, V = svd(Mm)
    @cast A[aⁱ⁻¹, aⁱ, σⁱ] |= U[(σⁱ, aⁱ⁻¹), aⁱ] (σⁱ:d)

    push!(atensors, A)

    MPS{L, T}(atensors), Left()
end

# Compression:1 ends here


text:
using Test
using MPI
using CLIMA.Mesh.BrickMesh

function main()
  MPI.Init()
  comm = MPI.COMM_WORLD
  rank = MPI.Comm_rank(comm)

  d = 4
  A = rand(1:10, d, 3rank)
  B = BrickMesh.parallelsortcolumns(comm, A, rev=true)

  root = 0
  Acounts = MPI.Allgather(Cint(length(A)), comm)
  A_all = MPI.Gatherv(A, Acounts, root, comm)

  Bcounts = MPI.Allgather(Cint(length(B)), comm)
  B_all = MPI.Gatherv(B, Bcounts, root, comm)

  if MPI.Comm_rank(comm) == root
    A_all = reshape(A_all, d, div(length(A_all),d))
    B_all = reshape(B_all, d, div(length(B_all),d))

    A_all = sortslices(A_all, dims=2, rev=true)

    @test A_all == B_all
  end
end

main()

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
# This file was generated, do not modify it. # hide
plt.figure(figsize=(8,6))
plt.hist(df.price, color = "blue", edgecolor = "white", bins=50,
         density=true, alpha=0.5)
plt.xlabel("Price", fontsize=14)
plt.ylabel("Frequency", fontsize=14)
plt.savefig(joinpath(@OUTPUT, "hist_price.svg")) # hide


text:
_uniqueid(a::AbstractVector)::Vector{Int} = indexin(a,unique(a)) # TODO: replace with indunique() or similar when it is introduced in julia (c.f. sort() and sortperm())

function _partitionpositions(referenceCodons::AbstractMatrix{Symbol})
    # for each position, assign an ID number to each reference s.t. they have the same ID if they coincide at that position
    IDPerPos = Vector{Int}[ _uniqueid(referenceCodons[i,:][:]) for i=1:size(referenceCodons,1) ]
    _uniqueid(IDPerPos) # give each position a number depending on which equivalence class it is in
end

function _partitionmatrices(partitionInd::AbstractVector{Int}, referenceCodons::AbstractMatrix{Symbol}, freqs::AbstractArray{Float64,3}, codons::Vector{Symbol})
    @assert size(referenceCodons,1)==size(freqs,2)

    nbrSamples = size(freqs,3)
    nbrReferences = size(referenceCodons,2)
    nbrPositions = length(partitionInd)

    nbrUniqueCodons = length(unique(referenceCodons[partitionInd[1],:])) # the number of unique codons is the same at all positions within a partition
    nbrUniqueCodons==1 && return zeros(0,nbrReferences), zeros(0,nbrSamples) # skip partition if it only has one unique codon per site (i.e. the partition doesn't help to determine the sample mix)

    f = zeros(nbrUniqueCodons,nbrPositions,nbrSamples) # UniqueCodons x Positions x Samples
    for (i,p) in enumerate(partitionInd) # for each position in the partition
        c = unique(referenceCodons[p,:]) # the codons belonging to the partition at the given site
        cInd = indexin(c, codons)
        f[:,i,:] = freqs[cInd,p,:] # extract slice for UniqueCodons x Samples
    end

    f = f ./ sum(f;dims=1) # normalize freqs to sum to one for each position and sample (i.e. conditional probability given that we have to chose one of the reference codons)

    # output trimmed frequencies (Y matrix)
    Y = zeros(nbrUniqueCodons,nbrSamples)
    for s=1:nbrSamples # for each sample
        # create masks for trimming
        tMask = trues(nbrPositions)
        for k=1:nbrUniqueCodons
            fs = f[k,:,s] # freqs associated to each unqiue group of reference genomes and sample
            q = quantile(fs, [0.25,0.75])
            tMask .&= tMask .& (q[1].<=fs) .& (fs.<=q[2])
        end
        # compute trimmed means
        Y[:,s] = mean(f[:,tMask,s]; dims=2)
    end

    # output reference->codon identification matrix (A matrix)
    # A_ij == 1 iff reference j has codon i and 0 otherwise
    A = zeros(nbrUniqueCodons,nbrReferences)
    setindex!.( (A,), 1, _uniqueid(referenceCodons[partitionInd[1],:][:]), 1:nbrReferences)

    # weight by number of positions in this group
    A *= nbrPositions
    Y *= nbrPositions

    A, Y
end

# TODO: support referenceGenomes::Array{Reference} too
function estimatemix(swarms::AnnotatedArray, referenceGenomes::Vector{T}) where {T<:AbstractString}
    @assert all( map(length,referenceGenomes) .== length(referenceGenomes[1]) )

    nbrReferences = length(referenceGenomes)
    nbrSamples = size(swarms,5)

    pos = swarms[:position][:]
    referenceCodons = hcat(seq2codons(referenceGenomes)...) # assumes all references have equal length
    referenceCodons = referenceCodons[pos,:]

    codons = reshape(swarms[:codon],64)
    freqs = reshape(convert(Array, swarms), 64, size(swarms,4), nbrSamples) # Codons x Positions x Samples

    partitions = _partitionpositions(referenceCodons)
    nbrPartitions = maximum(partitions)

    A = zeros(0,nbrReferences)
    Y = zeros(0,nbrSamples)

    for partition=1:nbrPartitions
        Ap, Yp = _partitionmatrices(findall(partition.==partitions), referenceCodons, freqs, codons)
        # concat
        A = vcat(A,Ap)
        Y = vcat(Y,Yp)
    end

    # solve separately for each sample for increased precision
    mixes = zeros(nbrSamples,nbrReferences)
    X = Variable(nbrReferences)
    for s=1:nbrSamples
        problem = minimize( sumsquares(A*X-Y[:,s]), [X>=0, ones(nbrReferences)'X==1] )
        solve!(problem, SCSSolver(eps=1e-6,max_iters=20000,verbose=false), verbose=false)
        problem.status==:Optimal || error("Failed to compute mix for sample ", swarms[:id][s])
        mixes[s,:] = (z=max.(X.value,0); z/sum(z)) # the solver might return slightly negative values, get rid of those
    end

    # # Fixing Y should improve performance, but the code below doesn't run faster even when we have hundreds of samples...
    # # solve separately for each sample for increased precision
    # mixes = zeros(nbrSamples,nbrReferences)
    # X = Variable(nbrReferences)
    # y = Variable(size(A,1))
    # problem = minimize( sumsquares(A*X-y), [X>=0, ones(nbrReferences)'X==1] )
    # for s=1:nbrSamples
    #     fix!(y,Y[:,s])
    #     solve!(problem, SCSSolver(eps=1e-6,verbose=false), verbose=false)#, warmstart=s>1)
    #     problem.status==:Optimal || error("Failed to compute mix for sample ", swarms[:id][s])
    #     mixes[s,:] = max(X.value,0) # the solver might return slightly negative values, get rid of those
    # end

    mixes
end


text:
struct D2Q37{
    T <: Real,
    Abscissaes <: AbstractMatrix{<:Integer},
    Weights <: AbstractVector{<:Real},
} <: Quadrature
    abscissae::Abscissaes
    weights::Weights
    speed_of_sound_squared::T
end

D2Q37() = D2Q37(
    [
        #  # group 2   # group 3   # group 4   # group 5               # group 6   # group 7   # group 8
        0 1 -1 0 0 1 -1 1 -1 2 -2 0 0 2 -2 -2 2 1 -1 1 -1 2 -2 2 -2 3 -3 0 0 3 -3 3 -3 1 -1 -1 1
        0 0 0 1 -1 1 -1 -1 1 0 0 2 -2 1 -1 1 -1 2 -2 -2 2 2 -2 -2 2 0 0 3 -3 1 -1 -1 1 3 -3 3 -3
    ],
    [
        0.23315066913235250228650

        # group 2
        0.10730609154221900241246
        0.10730609154221900241246
        0.10730609154221900241246
        0.10730609154221900241246

        # group 3
        0.05766785988879488203006
        0.05766785988879488203006
        0.05766785988879488203006
        0.05766785988879488203006

        # group 4
        0.01420821615845075026469
        0.01420821615845075026469
        0.01420821615845075026469
        0.01420821615845075026469

        # group 5
        0.00535304900051377523273
        0.00535304900051377523273
        0.00535304900051377523273
        0.00535304900051377523273
        0.00535304900051377523273
        0.00535304900051377523273
        0.00535304900051377523273
        0.00535304900051377523273

        # group 6
        0.00101193759267357547541
        0.00101193759267357547541
        0.00101193759267357547541
        0.00101193759267357547541

        # group 7
        0.00024530102775771734547
        0.00024530102775771734547
        0.00024530102775771734547
        0.00024530102775771734547

        # group 8
        0.00028341425299419821740
        0.00028341425299419821740
        0.00028341425299419821740
        0.00028341425299419821740
        0.00028341425299419821740
        0.00028341425299419821740
        0.00028341425299419821740
        0.00028341425299419821740
    ],
    # r =
    1.19697977039307435897239^2,
)

order(q::D2Q37) = 9
Base.show(io::IO, q::D2Q37) = show(io, "D2Q37")
Base.string(q::D2Q37) = "D2Q37"


text:
using LCMCore
using Test
using Random
using StaticArrays
using FastIOBuffers: FastReadBuffer, FastWriteBuffer, setdata!

include(joinpath(@__DIR__, "lcmtypes", "lcmtesttypes.jl"))

function test_encode_decode(::Type{T}) where T<:LCMType
    for i = 1 : 100
        in = rand(T)
        bytes = encode(in)::Vector{UInt8}
        out = decode(bytes, T)
        @test in == out
    end
end

function test_in_place_decode_noalloc(::Type{T}) where T<:LCMType
    in = rand(T)
    bytes = encode(in)
    out = deepcopy(in)
    decodestream = FastReadBuffer()
    setdata!(decodestream, bytes)
    decode!(out, decodestream)
    allocs = @allocated begin
        setdata!(decodestream, bytes)
        decode!(out, decodestream)
    end
    @test allocs == 0
end

function test_in_place_encode_noalloc(::Type{T}) where T<:LCMType
    in = rand(T)
    encodestream = FastWriteBuffer()
    encode(encodestream, in)
    take!(encodestream)
    allocs = @allocated encode(encodestream, in)
    @test allocs == 0
end

@testset "LCMType" begin
    using .lcmtesttypes
    Random.seed!(1)

    # Check base hashes (obtained from lcm-gen --debug; note that hash reported by lcm-gen is a **signed** integer (despite the 0x prefix))
    @test LCMCore.basehash(lcm_test_type_1) == 0x6ac981fb91a46d8b
    @test LCMCore.basehash(lcm_test_type_3) == 0x394ce5f0e401c325
    @test LCMCore.basehash(lcm_test_type_2) == 0x5a53eae01a55d4cb

    # Check invertibility
    test_encode_decode(lcm_test_type_1)
    test_encode_decode(lcm_test_type_2)
    test_encode_decode(lcm_test_type_3)
    test_encode_decode(polynomial_t)
    test_encode_decode(polynomial_matrix_t)

    # Check that decoding types without `String`s doesn't allocate
    test_in_place_decode_noalloc(lcm_test_type_1)
    test_in_place_decode_noalloc(lcm_test_type_2)
    test_in_place_decode_noalloc(polynomial_t)
    test_in_place_decode_noalloc(polynomial_matrix_t)

    # Check that encoding types doesn't allocate
    test_in_place_encode_noalloc(lcm_test_type_1)
    test_in_place_encode_noalloc(lcm_test_type_2)
    test_in_place_encode_noalloc(lcm_test_type_3)
    test_in_place_encode_noalloc(polynomial_t)
    test_in_place_encode_noalloc(polynomial_matrix_t)

    # Mismatch between length field and length of corresponding vector
    bad = rand(lcm_test_type_1)
    bad.blength += 1
    @test_throws DimensionMismatch encode(bad)

    # Bad fingerprint
    in = rand(lcm_test_type_1)
    badbytes = encode(in)
    badbytes[1] += 0x01
    @test_throws LCMCore.FingerprintException decode(badbytes, lcm_test_type_1)

    # Test against byte blobs that were encoded using pylcm
    for lcmt in [lcm_test_type_1, lcm_test_type_2, lcm_test_type_3, polynomial_t, polynomial_matrix_t]
        bytes = read(joinpath(@__DIR__, "lcmtypes", string(lcmt.name.name) * "_example_bytes"))
        @test hard_coded_example(lcmt) == decode(bytes, lcmt)
    end

    # resize!
    lcmt2 = lcm_test_type_2()
    lcmt2.dlength = 5
    lcmt2.f_inner_length = 2
    resize!(lcmt2)
    @test length(lcmt2.d) == 5
    @test size(lcmt2.f) == (3, 2)
    d = lcmt2.d
    f = lcmt2.f
    resize!(lcmt2)
    @test d === lcmt2.d
    @test f === lcmt2.f
end

@testset "LCMType: handle" begin
    channel = "CHANNEL_1"
    msg = rand(lcm_test_type_1)

    # start listening
    sublcm = LCM()
    check_msg = let expected = msg
        (channel, msg) -> @test(msg == expected)
    end
    sub = subscribe(sublcm, channel, check_msg, lcm_test_type_1)
    set_queue_capacity(sub, 2)

    # publish two messages
    publcm = LCM()
    for _ = 1 : 2
        publish(publcm, channel, msg)
    end

    # handle
    LCMCore.lcm_handle(sublcm)
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module PyLCM

using PyCall
export LCM, publish, subscribe, handle, @pyimport

const pylcm = PyNULL()

immutable LCM
	lcm_obj::PyObject

	LCM() = new(pylcm[:LCM]())
end

function publish(lc::LCM, channel::AbstractString, msg)
	pycall(lc.lcm_obj[:publish], PyAny, channel, pycall(msg[:encode], PyObject))
end

function subscribe(lc::LCM, channel::AbstractString, handler::Function)
	lc.lcm_obj[:subscribe](channel, pyeval("lambda chan, data, handler=h: handler(chan, bytearray(data))", h=handler))
end

function subscribe(lc::LCM, channel::AbstractString, handler::Function, msg_type::PyObject)
	lc.lcm_obj[:subscribe](channel, pyeval("lambda chan, data, handler=h, msg_type=t: handler(chan, msg_type.decode(data))", h=handler, t=msg_type))
end

function handle(lc::LCM)
	pycall(lc.lcm_obj[:handle], PyObject)
end

function __init__()
	depsjl = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
	isfile(depsjl) ? include(depsjl) : error("PyLCM not properly ",
	    "installed. Please run\nPkg.build(\"PyLCM\")")
    copy!(pylcm, pyimport("lcm"))
end

end # module


text:
VERSION >= v"0.4.0" && __precompile__(true)

module Indicators

export
    runmean, runsum, runvar, runsd, runcov, runcor, runmax, runmin, runmad,
    wilder_sum, mode, diffn, #lag,
    sma, trima, wma, ema, mma, kama, mama, hma, swma, dema, tema, alma,
    mlr_beta, mlr_slope, mlr_intercept, mlr, mlr_se, mlr_ub, mlr_lb, mlr_bands, mlr_rsq,
    momentum, roc, macd, rsi, adx, psar, kst, wpr, cci, stoch, smi,
    bbands, tr, atr, keltner,
    maxima, minima#, support, resistance

include("run.jl")
include("ma.jl")
include("reg.jl")
include("mom.jl")
include("vol.jl")
include("trendy.jl")

using Temporal
include("temporal.jl")

end


text:
module Sesam

using ProgressMeter
using CSV
using DataFrames
using Dates
using JuMP
using Gurobi

println("Using Sesam module.")

# utility functions
include("util.jl")

# power plants
include("powerplants.jl")

# renewables
include("renewables.jl")

# storages
include("storages.jl")

# nodes
include("nodes.jl")

# lines
include("lines_ptdf.jl")

### models
# economic dispatch
include("economic_dispatch.jl")

# congestion management
include("congestion_management_ptdf.jl")

# congestion management with PtG
include("congestion_management_ptg_ptdf.jl")

# helper functions
include("helpers.jl")

export
    CSV,
    DataFrames,
    JuMP,
    Gurobi,
    ProgressMeter,
    ElectronDisplay,
    Nodes,
    Lines,
    PowerPlants,
    Renewables,
    Storages,
    EconomicDispatch,
    CongestionManagementPHS80PtG,
    hoursInWeek,
    hoursInDay,
    hoursInHalfDay,
    has,
    getKeyVector,
    dayslicer,
    calc_h_b,
    CongestionManagementPHS,
    CongestionManagementPHS80,
    EconomicDispatch80,
    CongestionManagementPHS80_LoadShift

end


text:
@enum(Status::Int,
    INITIAL = 0,
    HAS_CONVERGED = 1,
    MAXIMUM_ITERATIONS = 2,
)

"""
    Parameters

This contains the parameters used in ADMM algorithm.
"""
mutable struct Parameters
    mu_max::Float64 # Augmented Lagrangian
    max_auglag::Int # Augmented Lagrangian
    ABSTOL::Float64
    RELTOL::Float64

    rho_max::Float64    # TODO: not used
    rho_min_pq::Float64 # TODO: not used
    rho_min_w::Float64  # TODO: not used
    eps_rp::Float64     # TODO: not used
    eps_rp_min::Float64 # TODO: not used
    rt_inc::Float64     # TODO: not used
    rt_dec::Float64     # TODO: not used
    eta::Float64        # TODO: not used
    verbose::Int

    # Two-Level ADMM
    outer_eps::Float64
    Kf::Int             # TODO: not used
    Kf_mean::Int        # TODO: not used
    MAX_MULTIPLIER::Float64
    DUAL_TOL::Float64
    rho_sigma::Float64  # Penalty for bus

    function Parameters()
        par = new()
        par.mu_max = 1e8
        par.rho_max = 1e6
        par.rho_min_pq = 5.0
        par.rho_min_w = 5.0
        par.eps_rp = 1e-4
        par.eps_rp_min = 1e-5
        par.rt_inc = 2.0
        par.rt_dec = 2.0
        par.eta = 0.99
        par.max_auglag = 50
        par.ABSTOL = 1e-6
        par.RELTOL = 1e-5
        par.verbose = 1
        par.outer_eps = 2*1e-4
        par.Kf = 100
        par.Kf_mean = 10
        par.MAX_MULTIPLIER = 1e12
        par.DUAL_TOL = 1e-8
        par.rho_sigma = 1e8
        return par
    end
end

abstract type AbstractGeneratorModel end

struct GeneratorModel{TD} <: AbstractGeneratorModel
    ngen::Int
    gen_start::Int
    pgmin::TD
    pgmax::TD
    qgmin::TD
    qgmax::TD
    c2::TD
    c1::TD
    c0::TD
end

"""
    Model{T,TD,TI}

This contains the parameters specific to ACOPF model instance.
"""
mutable struct Model{T,TD,TI}
    n::Int
    nline::Int
    nbus::Int
    nvar::Int

    line_start::Int

    gen_mod::AbstractGeneratorModel
    YshR::TD
    YshI::TD
    YffR::TD
    YffI::TD
    YftR::TD
    YftI::TD
    YttR::TD
    YttI::TD
    YtfR::TD
    YtfI::TD
    FrBound::TD
    ToBound::TD
    FrStart::TI
    FrIdx::TI
    ToStart::TI
    ToIdx::TI
    GenStart::TI
    GenIdx::TI
    Pd::TD
    Qd::TD

    # Two-Level ADMM
    nvar_u::Int
    nvar_v::Int
    bus_start::Int # this is for varibles of type v.
    brBusIdx::TI

    function Model{T,TD,TI}(data::OPFData, use_gpu::Bool, use_polar::Bool, use_twolevel::Bool) where {T, TD<:AbstractArray{T}, TI<:AbstractArray{Int}}
        model = new{T,TD,TI}()

        ngen = length(data.generators)
        gen_start = 1
        pgmin, pgmax, qgmin, qgmax, c2, c1, c0 = get_generator_data(data; use_gpu=use_gpu)
        model.gen_mod = GeneratorModel{TD}(ngen, gen_start, pgmin, pgmax, qgmin, qgmax, c2, c1, c0)

        model.n = (use_polar == true) ? 4 : 10
        model.nline = length(data.lines)
        model.nbus = length(data.buses)
        model.nvar = 2*ngen + 8*model.nline
        model.line_start = 2*ngen + 1
        model.YshR, model.YshI, model.YffR, model.YffI, model.YftR, model.YftI, model.YttR, model.YttI, model.YtfR, model.YtfI, model.FrBound, model.ToBound = get_branch_data(data; use_gpu=use_gpu)
        model.FrStart, model.FrIdx, model.ToStart, model.ToIdx, model.GenStart, model.GenIdx, model.Pd, model.Qd = get_bus_data(data; use_gpu=use_gpu)

        # These are only for two-level ADMM.
        model.nvar_u = 2*ngen + 8*model.nline
        model.nvar_v = 2*ngen + 4*model.nline + 2*model.nbus
        model.bus_start = 2*ngen + 4*model.nline + 1
        if use_twolevel
            model.nvar = model.nvar_u + model.nvar_v
            model.brBusIdx = get_branch_bus_index(data; use_gpu=use_gpu)
        end

        return model
    end
end

abstract type AbstractSolution{T,TD} end

"""
    SolutionOneLevel{T,TD}

This contains the solutions of ACOPF model instance, including the ADMM parameter rho.
"""
mutable struct SolutionOneLevel{T,TD} <: AbstractSolution{T,TD}
    status::Status
    u_curr::TD
    v_curr::TD
    l_curr::TD
    u_prev::TD
    v_prev::TD
    l_prev::TD
    rho::TD
    rd::TD
    rp::TD
    objval::T

    function SolutionOneLevel{T,TD}(model::Model) where {T, TD<:AbstractArray{T}}
        sol = new{T,TD}(
            INITIAL,
            TD(undef, model.nvar),
            TD(undef, model.nvar),
            TD(undef, model.nvar),
            TD(undef, model.nvar),
            TD(undef, model.nvar),
            TD(undef, model.nvar),
            TD(undef, model.nvar),
            TD(undef, model.nvar),
            TD(undef, model.nvar),
            Inf,
        )

        sol.u_curr .= 0.0
        sol.v_curr .= 0.0
        sol.l_curr .= 0.0
        sol.u_prev .= 0.0
        sol.v_prev .= 0.0
        sol.l_prev .= 0.0
        sol.rho .= 0.0
        sol.rd .= 0.0
        sol.rp .= 0.0

        return sol
    end
end

function active_power_generation(model::Model, sol::SolutionOneLevel)
    ngen = model.gen_mod.ngen
    return sol.u_curr[1:2:2*ngen]
end
function reactive_power_generation(model::Model, sol::SolutionOneLevel)
    ngen = model.gen_mod.ngen
    return sol.u_curr[2:2:2*ngen]
end
function voltage_magnitude(model::Model, sol::SolutionOneLevel)
    nbus = model.nbus
    return zeros(nbus)
end
function voltage_angle(model::Model, sol::SolutionOneLevel)
    nbus = model.nbus
    return zeros(nbus)
end

function set_active_power_generation!(model::Model, sol::SolutionOneLevel, val)
    ngen = model.gen_mod.ngen
    for g in 1:ngen
        pg_idx = model.gen_mod.gen_start + 2 * (g-1)
        @inbounds sol.v_curr[pg_idx] = val[g]
    end
end
function set_reactive_power_generation!(model::Model, sol::SolutionOneLevel, val)
    ngen = model.gen_mod.ngen
    for g in 1:ngen
        pg_idx = model.gen_mod.gen_start + 2 * (g-1)
        @inbounds sol.v_curr[pg_idx+1] = val[g]
    end
end
function set_voltage_magnitude!(model::Model, sol::SolutionOneLevel, val)
    # TODO
end
function set_voltage_angle!(model::Model, sol::SolutionOneLevel, val)
    # TODO
end

"""
    SolutionTwoLevel{T,TD}

This contains the solutions of ACOPF model instance for two-level ADMM algorithm,
    including the ADMM parameter rho.
"""
mutable struct SolutionTwoLevel{T,TD} <: AbstractSolution{T,TD}
    status::Status
    x_curr::TD
    xbar_curr::TD
    z_outer::TD
    z_curr::TD
    z_prev::TD
    l_curr::TD
    lz::TD
    rho::TD
    rp::TD
    rd::TD
    rp_old::TD
    Ax_plus_By::TD
    wRIij::TD
    s_curr::TD
    objval::T

    function SolutionTwoLevel{T,TD}(model::Model) where {T, TD<:AbstractArray{T}}
        sol = new{T,TD}(
            INITIAL,
            TD(undef, model.nvar),      # x_curr
            TD(undef, model.nvar_v),    # xbar_curr
            TD(undef, model.nvar),      # z_outer
            TD(undef, model.nvar),      # z_curr
            TD(undef, model.nvar),      # z_prev
            TD(undef, model.nvar),      # l_curr
            TD(undef, model.nvar),      # lz
            TD(undef, model.nvar),      # rho
            TD(undef, model.nvar),      # rp
            TD(undef, model.nvar),      # rd
            TD(undef, model.nvar),      # rp_old
            TD(undef, model.nvar),      # Ax_plus_By
            TD(undef, 2*model.nline),   # wRIij
            TD(undef, 2*model.nbus),    # s_curr
            Inf,                        # objval
        )
        sol.x_curr .= 0.0
        sol.xbar_curr .= 0.0
        sol.z_outer .= 0.0
        sol.z_curr .= 0.0
        sol.z_prev .= 0.0
        sol.l_curr .= 0.0
        sol.lz .= 0.0
        sol.rho .= 0.0
        sol.rp .= 0.0
        sol.rd .= 0.0
        sol.rp_old .= 0.0
        sol.Ax_plus_By .= 0.0
        sol.wRIij .= 0.0
        sol.s_curr .= 0.0

        return sol
    end
end

function active_power_generation(model::Model, sol::SolutionTwoLevel)
    ngen = model.gen_mod.ngen
    return sol.xbar_curr[1:2:2*ngen]
end
function reactive_power_generation(model::Model, sol::SolutionTwoLevel)
    ngen = model.gen_mod.ngen
    return sol.xbar_curr[2:2:2*ngen]
end
function voltage_magnitude(model::Model, sol::SolutionTwoLevel)
    bus_start = model.bus_start
    return sol.xbar_curr[bus_start:2:bus_start-1+2*model.nbus]
end
function voltage_angle(model::Model, sol::SolutionTwoLevel)
    bus_start = model.bus_start
    nbus = model.nbus
    return sol.xbar_curr[bus_start+1:2:bus_start-1+2*model.nbus]
end

function set_active_power_generation!(model::Model, sol::SolutionTwoLevel, val)
    ngen = model.gen_mod.ngen
    for g in 1:ngen
        pg_idx = model.gen_mod.gen_start + 2 * (g-1)
        @inbounds sol.xbar_curr[pg_idx] = val[g]
    end
end
function set_reactive_power_generation!(model::Model, sol::SolutionTwoLevel, val)
    ngen = model.gen_mod.ngen
    for g in 1:ngen
        pg_idx = model.gen_mod.gen_start + 2 * (g-1)
        @inbounds sol.xbar_curr[pg_idx+1] = val[g]
    end
end
function set_voltage_magnitude!(model::Model, sol::SolutionTwoLevel, val)
    bus_start = model.bus_start
    for b in 1:model.nbus
        bus_idx = bus_start + 2 * (b-1)
        @inbounds sol.xbar_curr[bus_idx] = val[b]
    end
end
function set_voltage_angle!(model::Model, sol::SolutionTwoLevel, val)
    bus_start = model.bus_start
    for b in 1:model.nbus
        bus_idx = bus_start + 2 * (b-1)
        @inbounds sol.xbar_curr[bus_idx+1] = val[b]
    end
end

"""
    AdmmEnv{T,TD,TI}

This structure carries everything required to run ADMM from a given solution.
"""
mutable struct AdmmEnv{T,TD,TI,TM}
    case::String
    data::OPFData
    use_gpu::Bool
    use_polar::Bool
    use_twolevel::Bool
    allow_infeas::Bool
    gpu_no::Int

    params::Parameters
    model::Model{T,TD,TI}
    solution::AbstractSolution{T,TD}

    membuf::TM # was param

    function AdmmEnv{T,TD,TI,TM}(
        opfdata, rho_pq, rho_va; use_gpu=false, use_polar=true, use_twolevel=false,
        allow_infeas=false, rho_sigma=1e8,
        gpu_no=-1, verbose=1, outer_eps=2e-4,
    ) where {T, TD<:AbstractArray{T}, TI<:AbstractArray{Int}, TM<:AbstractArray{T,2}}
        if use_gpu && gpu_no > 0
            CUDA.device!(gpu_no)
        end
        env = new{T,TD,TI,TM}()

        env.data = opfdata
        env.use_gpu = use_gpu
        env.use_polar = use_polar
        env.allow_infeas = allow_infeas
        env.gpu_no = gpu_no

        env.params = Parameters()
        env.params.verbose = verbose
        env.params.outer_eps = outer_eps
        env.params.rho_sigma = rho_sigma

        env.model = Model{T,TD,TI}(env.data, use_gpu, use_polar, use_twolevel)
        ybus = Ybus{Array{T}}(computeAdmitances(
            env.data.lines, env.data.buses, env.data.baseMVA; VI=Array{Int}, VD=Array{T})...)

        env.solution = ifelse(use_twolevel,
            SolutionTwoLevel{T,TD}(env.model),
            SolutionOneLevel{T,TD}(env.model))
        init_solution!(env, env.solution, ybus, rho_pq, rho_va)

        env.membuf = TM(undef, (31, env.model.nline))
        fill!(env.membuf, 0.0)

        return env
    end
end

function AdmmEnv(case::String, use_gpu::Bool, rho_pq, rho_va; options...)
    opfdata = opf_loaddata(case)
    env = AdmmEnv(opfdata, use_gpu, rho_pq, rho_va; options...)
    env.case = case
    return env
end

function AdmmEnv(opfdata::OPFData, use_gpu::Bool, rho_pq, rho_va; options...)
    VT = use_gpu ? CuArray : Array
    return AdmmEnv{Float64, VT{Float64, 1}, VT{Int, 1}, VT{Float64, 2}}(
        opfdata, rho_pq, rho_va; use_gpu=use_gpu, options...
    )
end

function AdmmEnv(opfdata::OPFData, ::Type{VT}, rho_pq, rho_va; options...) where VT
    use_gpu = VT <: CuArray
    return AdmmEnv{Float64, VT{Float64, 1}, VT{Int, 1}, VT{Float64, 2}}(
        opfdata, rho_pq, rho_va; use_gpu=use_gpu, options...
    )
end

# Getters / setters
active_power_generation(env::AdmmEnv) = active_power_generation(env.model, env.solution)
reactive_power_generation(env::AdmmEnv) = reactive_power_generation(env.model, env.solution)

function set_active_load!(env::AdmmEnv, pd::AbstractArray)
    @assert length(pd) == env.model.nbus
    copyto!(env.model.Pd, pd)
    return
end

function set_reactive_load!(env::AdmmEnv, qd::AbstractArray)
    @assert length(qd) == env.model.nbus
    copyto!(env.model.Qd, qd)
    return
end

set_active_power_generation!(env::AdmmEnv, val) = set_active_power_generation!(env.model, env.solution, val)
set_reactive_power_generation!(env::AdmmEnv, val) = set_reactive_power_generation!(env.model, env.solution, val)
set_voltage_magnitude!(env::AdmmEnv, val) = set_voltage_magnitude!(env.model, env.solution, val)
set_voltage_angle!(env::AdmmEnv, val) = set_voltage_angle!(env.model, env.solution, val)


text:
""" Computing in the category of finite sets and Petri cospans
"""
module AlgebraicPetri
export PetriCospanOb, PetriFunctor, PetriCospan

using Catlab.GAT
using Catlab.Theories: BiproductCategory
using Catlab.CategoricalAlgebra.ShapeDiagrams
using Catlab.CategoricalAlgebra.FinSets
using Petri
using AutoHashEquals

import Catlab.Theories: dom, codom, id, compose, ⋅, ∘, otimes, ⊗, munit,
                        braid, σ, mcopy, Δ, mmerge, ∇, create, □, delete, ◊,
                        pair, copair, proj1, proj2, coproj1, coproj2

""" Finite ordinal (natural number)

An object in the category of Open Petri Nets.
"""
@auto_hash_equals struct PetriCospanOb
    n::Int
end
Base.eachindex(X::PetriCospanOb) = 1:X.n

struct PetriDecorator <: AbstractFunctor end
struct PetriLaxator <: AbstractLaxator end

""" Petri Functor

A functor from FinOrd to Petri defined as a PetriDecorator and a PetriLaxator
"""
const PetriFunctor = LaxMonoidalFunctor{PetriDecorator, PetriLaxator}

id(::Type{PetriFunctor}) = PetriFunctor(PetriDecorator(), PetriLaxator())

""" AlgebraicPetri.PetriDecorator(n::FinOrd)

A functor from FinOrd to Set has an objects part, which given an object n in
FinOrd (a natural number) should return a representation of F(n)::Set, sets can
be represented as a predicate that takes an element and returns true if the
element is in the set. Here we take any julia value and test whether it is a
Petri net on n states.
"""
function (pd::PetriDecorator)(n::FinOrd)
    return p -> typeof(p) <: Petri.Model && length(p.S) == n.n
end

map(f::Function, d::Dict{K,V}) where {K,V} = begin
    out = Dict{K,V}()
    for p in Base.map(x->Pair(f(x[1]), x[2]), collect(d))
        if p[1] in keys(out)
            out[p[1]] += p[2]
        else
            out[p[1]] = p[2]
        end
    end
    out
end

map(f::Function, ts::Vector{Tuple{S,T}}) where {S<:Dict, T<:Dict} = [(map(f, t[1]), map(f, t[2])) for t in ts]

""" AlgebraicPetri.PetriDecorator(f::FinOrdFunction)

A functor from FinOrd to Set has a hom part, which given a hom f in FinOrd
(a function n::Int->m::Int) should return a representation of F(f)::F(n)->F(m),
here we implement this as a function that takes a Petri net of size n to a Petri
net of size m, such that the transitions are mapped appropriately.
"""
function (pd::PetriDecorator)(f::FinOrdFunction)
    return (p::Petri.Model) -> Petri.Model(collect(1:codom(f).n), map(x->f(x), p.Δ))
end

""" AlgebraicPetri.PetriLaxator(p::Petri.Model, q::Petri.Model)

The laxitor takes a pair of decorations and returns the coproduct decoration
For Petri nets, this encodes the idea that you shift the states of q up by the
number of states in p.
"""
function (l::PetriLaxator)(p::Petri.Model, q::Petri.Model)
    return Petri.Model(collect(1:(length(p.S)+length(q.S))),
                       vcat(p.Δ, map(x->x+length(p.S), q.Δ)))
end

""" Petri Cospan

A morphism in the category of Open Petri Nets defined as a decorated cospan with
a [`PetriFunctor`](@ref) as the decorator which maps the category of finite
ordinals to the category Petri and a Petri.Model as the decoration
"""
const PetriCospan = DecoratedCospan{PetriFunctor, Petri.Model}

""" AlgebraicPetri.PetriCospan(l::Vector{Int}, m::Petri.Model, r::Vector{Int})

A constructor for Petri Cospans where `l` is a vector of the input states from
Petri.Model `m`, and `r` is a vector of the output states from Petri.Model `m`

Constructs the cospan: l → m ← r
"""
function (::Type{PetriCospan})(l::AbstractVector, m::Petri.Model, r::AbstractVector)
    return PetriCospan(Cospan(FinOrdFunction(l, length(m.S)),
                              FinOrdFunction(r, length(m.S))),
                       id(PetriFunctor), m)
end

@instance BiproductCategory(PetriCospanOb, PetriCospan) begin
    dom(f::PetriCospan) = PetriCospanOb(dom(left(f)).n)
    codom(f::PetriCospan) = PetriCospanOb(dom(right(f)).n)

    compose(p::PetriCospan, q::PetriCospan) = begin
        # reimplementation of pushout of Span{FinOrdFunc, FinOrdFun}
        # to save the value of coeq
        f, g = right(p), left(q)
        coprod = coproduct(codom(f), codom(g))
        ι1, ι2 = left(coprod), right(coprod)
        coeq = coequalizer(f⋅ι1, g⋅ι2)
        f′, g′ = ι1⋅coeq, ι2⋅coeq
        composite = Cospan(left(p)⋅f′, right(q)⋅g′)
        dpuq = decorator(p).L(decoration(p), decoration(q))
        return PetriCospan(composite, decorator(p), decorator(p).F(coeq)(dpuq))
    end

    id(X::PetriCospanOb) = PetriCospan(
        Cospan(id(FinOrd(X.n)), id(FinOrd(X.n))),
        id(PetriFunctor),
        EmptyPetri(X.n))

    otimes(X::PetriCospanOb, Y::PetriCospanOb) = PetriCospanOb(X.n + Y.n)

    otimes(f::PetriCospan, g::PetriCospan) = begin
        fl, fr = left(f), right(f)
        gl, gr = left(g), right(g)
        # TODO: Replace with universal properties once implemented
        PetriCospan(
            Cospan(
                FinOrdFunction(x->x > dom(fl).n ? gl(x-dom(fl).n)+codom(fl).n : fl(x),
                               FinOrd(dom(fl).n + dom(gl).n),
                               FinOrd(codom(fl).n + codom(gl).n)),
                FinOrdFunction(x->x > dom(fr).n ? gr(x-dom(fr).n)+codom(fr).n : fr(x),
                               FinOrd(dom(fr).n + dom(gr).n),
                               FinOrd(codom(fr).n + codom(gr).n))
            ), decorator(f), decorator(f).L(decoration(f), decoration(g)))
    end

    munit(::Type{PetriCospanOb}) = PetriCospanOb(0)

    braid(X::PetriCospanOb, Y::PetriCospanOb) = begin
        Z = otimes(X, Y)
        PetriCospan(
            Cospan(
                id(FinOrd(Z.n)),
                FinOrdFunction(vcat(X.n+1:Z.n, 1:X.n), Z.n, Z.n)
            ), id(PetriFunctor), EmptyPetri(Z.n))
    end

    mcopy(X::PetriCospanOb) = PetriCospan(
        Cospan(
            id(FinOrd(X.n)),
            FinOrdFunction(vcat(1:X.n,1:X.n), 2*X.n, X.n)
        ), id(PetriFunctor), EmptyPetri(X.n))

    mmerge(X::PetriCospanOb) = PetriCospan(
        Cospan(
            FinOrdFunction(vcat(1:X.n,1:X.n), 2*X.n, X.n),
            id(FinOrd(X.n))
        ), id(PetriFunctor), EmptyPetri(X.n))

    create(X::PetriCospanOb) = PetriCospan(
        Cospan(FinOrdFunction(Int[], 0, X.n), id(FinOrd(X.n))),
        id(PetriFunctor), EmptyPetri(X.n))

    delete(X::PetriCospanOb) = PetriCospan(
        Cospan(id(FinOrd(X.n)), FinOrdFunction(Int[], 0, X.n)),
        id(PetriFunctor), EmptyPetri(X.n))

    pair(f::PetriCospan, g::PetriCospan) = compose(mcopy(dom(f)), otimes(f, g))
    copair(f::PetriCospan, g::PetriCospan) = compose(otimes(f, g), mmerge(codom(f)))

    proj1(A::PetriCospanOb, B::PetriCospanOb) = otimes(id(A), delete(B))
    proj2(A::PetriCospanOb, B::PetriCospanOb) = otimes(delete(A), id(B))

    coproj1(A::PetriCospanOb, B::PetriCospanOb) = otimes(id(A), create(B))
    coproj2(A::PetriCospanOb, B::PetriCospanOb) = otimes(create(A), id(B))
end

include("Epidemiology.jl")

end


text:
module NLoptExamples
import NLopt
import ForwardDiff

include("quadratic.jl")
include("data_fit.jl")

export quadratic_ex
export optim_data
end # module


text:
# Helper function
function convertTo(m::Money, target::Currency)
	if m.currency!=target
		rate=lookup(m.currency, target)
		return rounded(exchange(rate, m))
	end
	return m
end

convertToBase(m::Money)=convertTo(m, baseCurrency(m))

function functionOperator(m1::Money, m2::Money, op::Function)
	if m1.currency==m2.currency
		return op(m1.value, m2.value)
	elseif ConversionType==:BaseCurrencyConversion
		tmp1=convertToBase(m1)
		tmp2=convertToBase(m2)
		return functionOperator(tmp1, tmp2, op)
	elseif ConversionType==:AutomatedConversion
		tmp=convertTo(m2, m1.currency)
		return functionOperator(m1, tmp, op)
	else
		error("currency mismatch and no conversion specified")
	end
end

# Operator functions that return a value
/(m1::Money, m2::Money)=functionOperator(m1, m2, /)
==(m1::Money, m2::Money)=functionOperator(m1, m2, ==)
<(m1::Money, m2::Money)=functionOperator(m1, m2, <)
close(m1::Money, m2::Money)=functionOperator(m1, m2, (x,y)=>abs(x-y)<10.0^(-10))

function functionOperatorMoney(m1::Money, m2::Money, op::Function)
	if m1.currency==m2.currency
		return Money(m1.currency, op(m1.value, m2.value))
	elseif ConversionType==:BaseCurrencyConversion
		tmp1=convertToBase(m1)
		tmp2=convertToBase(m2)
		return functionOperatorMoney(tmp1, tmp2, op)
	elseif ConversionType==:AutomatedConversion
		tmp=convertTo(m2, m1.currency)
		return functionOperatorMoney(m1, tmp, op)
	else
		error("currency mismatch and no conversion specified")
	end
end

# Operator functions that return a Money
+(m1::Money, m2::Money)=functionOperatorMoney(m1, m2, +)
-(m1::Money, m2::Money)=functionOperatorMoney(m1, m2, -)