text:
#
# Version 0.1 (IJCNN submission preparation)
#
#     06 Nov 2016: - Added function to add noise to the generated signal
#
# (A lot of functions are copied from Pccf/ repository)
module SignalGenerator

include("recmidpoints.jl")

export genSig, addNoise!, addNoiseRec!

function bisect(v, x, imin::Int, imax::Int)
    # https://docs.python.org/2/library/bisect.html
    # The returned insertion point i partitions the sorted vector v into two halves so that
    # all(val <= x for val in v[1:i]) for the left side and
    # all(val > x  for val in v[i:hi]) for the right side.
    if imin == imax
        x <= v[imin] ? imin : imin + 1
    else
        imid = div(imin+imax, 2)
        if x == v[imid]
            return imid
        elseif x > v[imid]
            bisect(v, x, imid + 1, imax)
        else
            bisect(v, x, imin, max(imid - 1, imin))
        end
    end
end
function bisect(v, x)
    bisect(v,x,1,length(v))
end

"""
Function to generate sequence of elements according to the probabilities proportional to the weights.
"""
function genSamples(values, weights::Vector{Float64}, n::Int64)
    @assert length(values) == length(weights)
    if sum(weights) == 1.0
        probs = weights
    else
        probs = weights/sum(weights)
    end
    cdf = cumsum(probs)
    samples = []
    for i in 1:n
        u   = rand()
        i   = bisect(cdf, u)
        push!(samples, values[i])
    end
    return samples
end

## ;;;;;;;;;;;;;;;;;Pccf/v1.0/src/recurrency_modelling2/sim3.clj;;;;;;;
## (defn gen-vals-in-interval-uniform [n from-v to-v]
##   "Function to generate values within interval inclusive"
##   (let [rng (inc (- to-v from-v))]
##     (map #(+ % from-v) (repeatedly n #(rand-int rng)))))
##
## (defn gen-vals-in-interval-gauss [n from-v to-v]
##   (let [mu (/ (+ to-v from-v) 2.0) sigma (inc (- to-v from-v))]
##     (utils/rnorm n mu sigma)))
## ;; (gen-vals-in-interval 5 7 8)
## ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

# """
# Probably it if naive experiments V.1.1 (After https://github.com/av-maslov/Pccf)
# Function to generate random observations.
# bootom_level, upper_level - boundaries of the interval defining mean and sigma
# Boundaries seemed to be more convenient than setting mu and sigma directly.
# """
# function genValsInIntervalGauss(n, bootom_level, upper_level)
#     # old: mu = (n+st)/2
#     mu = (bootom_level + upper_level)/2
#     sigma = upper_level - bootom_level + 1.0
#     randn(n) * sigma + mu
# end

function calcWidth(mu)
    1 + round(Int, 1.33 * mu)
end

function gaussPdf(x, mu, sigma)
    k = 1.0/(sigma * sqrt(2*pi))
    p = k * exp(-((x-mu)^2/(2*sigma^2)))
end

"Change's locations, all > 0 except the first"
function changesLocations(n, mu, sigma)
    locations = vcat(0, filter(x -> x > 0, cumsum(randn(n)*sigma + mu)))
    locations = map(x -> convert(Int, x), round(locations))
end

"Length of the signal to be generated."
function calcTmax(mu, n)
    round(Int,(mu * (n+1)))
end

function errorsLocations(n,mu,errRate)
    tMax = calcTmax(mu, n)
    avgNumOfErrors = errRate * n
    rand(1:tMax, avgNumOfErrors)
end

"""
Function to generate the signal with changes and noise.
n          : length of the signal
mu         : avg distance between changes
sigma      : std. dev. of the distance
changeStep : absolute value of the signal's level change after each changepoint
sigStd     : signals standard deviation
"""
function genSig(n, mu, sigma, changeStep, sigStd)
    tMax    = calcTmax(mu, n)
    changes = changesLocations(n, mu, sigma)
    n_chps  = length(changes)
    levels  = cumsum(rand([-1, 1], n_chps) * changeStep/2.0)
    @assert length(levels) == n_chps
    sig = []
    for i in 2:n_chps
        sub_len = changes[i] - changes[i-1]
        sig     = vcat(sig, randn(sub_len) * sigStd + levels[i])
    end
    Dict("sig"=>sig, "changes"=>changes)
end

function updVec(vec, ids::Array{Int,1}, newVals)
    v=copy(vec)
    for i in 1:length(ids)
        v[ids[i]]=newVals[i]
    end
    return v
end

"""
Add noise between changepoinst in the artificially generated signal.
Just in the middle between changes.
sig[i] = sig[i] + changeStep
"""
function addNoise!(sig, changes, changeStep; width = 5)
    idx_noise = Array{Int64,1}()
    n_chp = length(changes)
    n     = length(sig)
    @assert n_chp > 3
    for i in 2 : n_chp
        c_j = changes[i-1]
        c_i = changes[i]
        m = convert(Int, div(c_i+c_j, 2))
        idx_noise = vcat(idx_noise, collect(m-div(width,2) : m+div(width,2)))
    end
    for k in idx_noise
        sig[k] = sig[k] + changeStep
    end
    idx_noise
end

"""
Use recursive function to add more noise than just in the middle
"""
function addNoiseRec!(sig, changes, changeStep; width = 5, min_len = 30)
    idx_noise = Array{Int64,1}()
    n_chp = length(changes)
    n     = length(sig)
    @assert n_chp > 3
    for i in 2 : n_chp
        c_j = changes[i-1]
        c_i = changes[i]
        # m = convert(Int, div(c_i+c_j, 2))
        midpoints = recmidpoints(c_j, c_i, min_len)
        for m in midpoints
            idx_noise = vcat(idx_noise, collect(m-div(width,2) : m+div(width,2)))
        end
    end
    sort!(idx_noise)
    for k in idx_noise
        sig[k] = sig[k] + changeStep
    end
    idx_noise
end

end

# """
# Function to generate the signal with changes and noise.
# For changes values are from tvf-tvt; For errors in rvf-evt
# tvf: 'true values from' / tvt: `true values to`
# evf: 'error's values from' / evt: 'error's values to'
# """
# function genSig_old(n, mu, sigma, errRate, tvf, tvt, evf, evt)
#     tMax        = calcTmax(mu,n)
#     changesLoc  = changesLocations(n,mu,sigma)
#     nChps       = length(changesLoc)
#     valsChanges = genValsInIntervalGauss(nChps, tvf, tvt)
#     # ???
#     vals        = zeros(tMax)
#
#     if errRate > 0
#         errorsLoc   = errorsLocations(n,mu,errRate)
#         nErrors     = length(errorsLoc)
#         valsErrors  = genValsInIntervalGauss(nErrors, evf, evt)
#     end
#     changes1    = updVec(vals, errorsLoc, valsErrors)
#     changes2    = updVec(changes1, changesLoc, valsChanges)
#     Dict("changes"=>changesLoc, "sig"=>changes2)
# end

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
#pkg"add CmdStan"

using CmdStan, Statistics
bm = "
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
    y ~ bernoulli(theta);
  }
";
data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])
sm = Stanmodel(name="bernoulli", model=bm, printsummary=false)
sm |> display

rc, samples, cnames = stan(sm, data; CmdStanDir=CMDSTAN_HOME)
if rc == 0
  nt = CmdStan.convert_a3d(samples, cnames, Val(:namedtuple))
  mean(nt.theta, dims=1) |> display
  sdf = read_summary(sm)
  sdf |> display
end

ext: jl


text:
@enum GOPCODES begin
    DISPATCH = 0
    HEARTBEAT
    IDENTIFY
    PRES_UPD
    VOICE_STATE_UPD
    RESUME = 6
    RECONNECT
    REQ_GUILD_MEMB
    INVALID_SES
    HELLO
    HEARTBEAT_ACK
end

const CONN_PROPS = Dict{String,String}(
    "\$os"      => string(Sys.KERNEL),
    "\$browser" => "discord-bot.jl",
    "\$device"  => "discord-bot.jl",
)

# TODO: CHANGE IT COMPLETLY
function dispatch!(c::Client, data::Dict{Symbol,Any})
    return try
        c.seq_num = data[:s]
        handler = get_handler(c.event_h, data[:t])
        handler  === DefaultHandler() && error("no handler found")
        handler(c, data[:d])
        nothing
    catch err
        err_msg = if err isa ErrorException
            err.msg
        else
            sprint(show, err)
        end
        ErrorException("on handling $(data[:t]): $err_msg")
    end
end

function heartbeat(c::Client, args...)
    _, err = write_json(c.conn.io, Dict{Symbol,Any}(
            :op => Int(HEARTBEAT),
            :d => iszero(c.seq_num) ? "null" : c.seq_num
        )
    )
    err
end

function hello(c::Client, args...)
    _, err = write_json(c.conn.io, Dict{Symbol,Any}(
            :op => Int(IDENTIFY),
            :d => Dict{Symbol,Any}(
                :token => c.token,
                :intents => c.intents, # TODO: SET MEANINGFUL VALUE
                :properties => CONN_PROPS
            )
        )
    )
    err
end

heartbeat_ack!(c::Client, args...) = (set_ack!(c.hb); nothing)

function resume_conn(c::Client, args...)
    _, err = write_json(c.conn.io, Dict{Symbol,Any}(
            :op => Int(RESUME),
            :d => Dict{Symbol,Any}(
                :token => c.token,
                :session_id => c.session_id,
                :seq => c.seq_num,
            )
        )
    )
    isnothing(err) || return err
    data, err = read_json(c.conn.io)
    isnothing(err) || return err
    GOPCODES(data[:op]) == INVALID_SES ? ErrorException("on getting invalid session opcode") : nothing
end

function reconnect!(c::Client, args...)
    close!(c) # Do we need to do that?
    return open!(c; resume=true, reconn=true)
end

function invalid_ses!(c::Client, args...)
    c.session_id = ""
    c.seq_num = 0
    ErrorException("on receiving an invalid session signal")
end

# The map of opcodes Client may recieve & corresponding handlers
const HANDLERS = Dict{GOPCODES,Function}(
    DISPATCH => dispatch!,
    HEARTBEAT => heartbeat,
    RECONNECT => reconnect!,
    HELLO => hello,
    HEARTBEAT_ACK => heartbeat_ack!,
    INVALID_SES => invalid_ses!,
)

handle(op::Int) = HANDLERS[GOPCODES(op)]


text:
# Sphere type and related methods

"""
    Sphere <: Shape

Sphere shape model.

"""
struct Sphere <: Shape
    id::Int64
    radius::Float64
    angle::Float64
    obj_to_world::Transformation
    world_to_obj::Transformation
end

Sphere() = Sphere(1, 1.0, π, Transformation(), Transformation())
Sphere(T::Transformation) = Sphere(1, 1.0, π, T, inv(T))
Sphere(r::Real, T::Transformation) = Sphere(1, r, π, T, inv(T))
Sphere(r::Real, angle::Real, T::Transformation) = Sphere(1, r, angle, T, inv(T))
Sphere(id::Int64, r::Real, T::Transformation) = Sphere(id, r, π, T, inv(T))
Sphere(id::Int64, r::Real, angle::Real, T::Transformation) = Sphere(id, r, angle, T, inv(T))

"""
    can_intersect(s::Sphere)

Is the sphere an intersectable object? Always returns true.

"""
can_intersect(s::Sphere) = true

"""
    area(s::Sphere)

Compute surface area of a sphere (doesn't into account scaling...)

TODO: fix for angle (pen and paper required)

"""
area(s::Sphere) = 4π * s.radius^2

"""
    obj_bounds(s::Sphere)

Return bounding box of the sphere in object coordinates.

"""
function obj_bounds(s::Sphere)
    v = Point3(s.radius, s.radius, s.radius)
    return BoundingBox(-v, v)
end

"""
    world_bounds(s::Sphere)

Return bounding box of the sphere in world coordinates.

"""
function world_bounds(s::Sphere)
    s.obj_to_world(obj_bounds(s))
end

"""
    (T::Transformation)(S::Sphere)

Apply a transformation to a sphere.

"""
function (T::Transformation)(S::Sphere)
    T2 = T * S.obj_to_world
    Sphere(S.id, S.radius, S.angle, T2, inv(T2))
end

function sphere_coefs(ray::Ray, sph::Sphere)
    A = dot(ray.direction, ray.direction)
    B = 2.0 * dot(ray.direction, ray.origin)
    C = dot(ray.origin, ray.origin) - sph.radius^2
    return A, B, C
end

"""
    shape_intersect(r::Ray, sph::Sphere)

Compute intersection geometry between given ray and sphere.

"""
function shape_intersect(r::Ray, sph::Sphere)
    ray = sph.world_to_obj(r)

    hit, tnear, tfar = quadratic(sphere_coefs(ray, sph)...)

    # Not hit or both hits out of bounds
    if !hit || tnear > r.tmax || tfar < r.tmin || (tnear < r.tmin && tfar > r.tmax)
        return nothing, NaN, NaN
    end
    P1 = ray(tnear)
    P2 = ray(tfar)
    if P1.z / sph.radius < cos(sph.angle)
        if P2.z / sph.radius < cos(sph.angle)
            return nothing, NaN, NaN
        else
            P = P2
            t = tfar
        end
    else
        P = P1
        t = tnear
    end

    n = normalize(Normal3(P.x, P.y, P.z))
    s = normalize(Vector3(-2π*P.y, 2π*P.x, 0.0))

    n = flip_normal(n, ray.direction)

    DG = DifferentialGeometry(sph.obj_to_world(P), sph.obj_to_world(n), sph.obj_to_world(s))
    return DG, t, 5e-4 * t
end

"""
    intersectP(r::Ray, sph::Sphere)

Quick true/false intersection test between ray and sphere.

"""
function intersectP(r::Ray, sph::Sphere)
    ray = sph.world_to_obj(r)

    hit, tnear, tfar = quadratic(sphere_coefs(ray, sph)...)

    # Not hit or both hits out of bounds
    if !hit || tnear > r.tmax || tfar < r.tmin || (tnear < r.tmin && tfar > r.tmax)
        return false
    end

    P1 = ray(tnear)
    P2 = ray(tfar)
    if P1.z / sph.radius < cos(sph.angle)
        if P2.z / sph.radius < cos(sph.angle)
            return false
        else
            return true
        end
    end

    return true
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
# # Sliding window analysis

# In this example, we will get precipitation data from Québec, and use a sliding
# window analysis to smooth them out. The beginning of the code should now be
# familiar:

using SimpleSDMLayers
using Plots
using Statistics

#-

precipitation = SimpleSDMPredictor(WorldClim, BioClim, 12; left=-80.0, right=-56.0, bottom=44.0, top=62.0)

# The sliding window works by taking all pixels *within a given radius*
# (expressed in kilometres) around the pixel of interest, and then applying the
# function given as the second argument to their values. Empty pixels are
# removed. In this case, we will do a summary across a 100 km radius around each
# pixel:

averaged = slidingwindow(precipitation, Statistics.mean, 100.0; threaded=false)

# We can finally overlap the two layers -- the result of sliding window is a
# little bit smoother than the raw data.

plot(precipitation, c=:alpine)
contour!(averaged, c=:white, lw=2.0)


text:
function run_server()
    # Start health check - necessary for AWS
    @async begin
        http = HttpHandler() do req::Request, res::Response
            Response(string(JSON.json({
                "status" =>1,
                "env" => env,
                "hello" => "world",
            })))
        end
        http.events["listen"] = (port) -> Logging.info("Health check listening on port $port")

        server = Server(http)
        run(server, config["port"])
    end

    Logging.info("Starting server")
    while true
        task = get_task()
        if task != nothing
            run_task(task)
            mark_task_done(task)
        else
            Logging.debug("No task. Sleeping.")
            sleep(config["sleep"])
        end
    end
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
using QuasiArrays, Test
import QuasiArrays: ApplyQuasiMatrix

@testset "Hcat" begin
    A = QuasiArray(randn(2,2), ([0,0.5], Base.OneTo(2)))
    B = [A A]
    @test B isa ApplyQuasiMatrix{Float64,typeof(hcat)}
    @test axes(B) == (axes(A,1), Base.OneTo(4))
    @test B[0.0,1] == B[0.0,3] == A[0.0,1]
    @test_throws BoundsError B[0.1,1]
    @test_throws BoundsError B[0.0,5]
    @test QuasiArray(B) == B
end


text:
"Calculate the largest needed order for the hankel series."
function maximum_hankel_order(ω::Union{T,Complex{T}}, medium::Medium{T}, species::Vector{Specie{T}};
        tol::T=1e-7, verbose::Bool = false) where T <: Number

    # estimation is based on far-field scattering pattern contribution to the primary effective wavenumber^2 and reflection coefficient.
    f0 = ω^(2.0)/(medium.c^2)
    hankel_order = 0
    next_order = -4im*sum(sp.num_density*Zn(ω,sp,medium,hankel_order) for sp in species)

    # increase hankel order until the relative error^2 < tol. Multiplying the tol by 200 has been chosen based on empircal comparisons with other tolerance used for reflection coefficients.
    while abs(next_order/f0) > tol*200
        f0 = f0 + next_order
        hankel_order += 1
        next_order = -4im*sum(sp.num_density*Zn(ω,sp,medium,hankel_order) for sp in species)
    end

    return hankel_order
end

"A t_matrix in the form of a vector, because for now we only deal with diagonal T matrices."
function t_vectors(ω::T, medium::Medium{T}, species::Vector{Specie{T}}; hankel_order = 3, dim = 2) where T <: AbstractFloat
    t_vecs = [ zeros(Complex{T},1+2hankel_order) for s in species]
    for i = 1:length(species), n = 0:hankel_order
        t_vecs[i][n+hankel_order+1] = - Zn(ω,species[i],medium,n; dim = dim)
        t_vecs[i][-n+hankel_order+1] = t_vecs[i][n+hankel_order+1]
    end
    return t_vecs
end

"Pre-calculate a matrix of Zn's"
function Zn_matrix(ω::T, medium::Medium{T}, species::Vector{Specie{T}}; hankel_order = 3) where T <: Number
    Zs = OffsetArray{Complex{T}}(undef, 1:length(species), -hankel_order:hankel_order)
    for i = 1:length(species), n = 0:hankel_order
        Zs[i,n] = Zn(ω,species[i],medium,n)
        Zs[i,-n] = Zs[i,n]
    end
    return Zs
end

Zn(ω::T, p::Specie{T}, med::Medium{T}, m::Int; kws...) where T<: Number = Zn(Complex{T}(ω), p, med, m; kws...)

"Returns a ratio which gives the scattering strength of a particle"
function Zn(ω::Complex{T}, p::Specie{T}, med::Medium{T},  m::Int; dim::Int = 2) where T<: Number
    m = abs(m)
    ak = p.r*ω/med.c
    # check for material properties that don't make sense or haven't been implemented
    if abs(p.c*p.ρ) == T(NaN)
        error("scattering from a particle with density =$(p.ρ) and phase speed =$(p.c) is not defined")
    elseif abs(med.c*med.ρ) == T(NaN)
        error("wave propagation in a medium with density =$(med.ρ) and phase speed =$(med.c) is not defined")
    elseif abs(med.c) == zero(T)
        error("wave propagation in a medium with phase speed =$(med.c) is not defined")
    elseif abs(med.ρ) == zero(T) && abs(p.c*p.ρ) == zero(T)
        error("scattering in a medium with density $(med.ρ) and a particle with density =$(p.ρ) and phase speed =$(p.c) is not defined")
    end
    if dim == 3
        j = sbesselj
        dj(m,x) = diffsbessel(sbesselj,m,x)
        h1 = shankelh1
        dh1(m,x) = diffsbessel(shankelh1,m,x)
    else
        j = besselj
        dj = diffbesselj
        h1 = hankelh1
        dh1 = diffhankelh1
    end
    # set the scattering strength and type
    if abs(p.c) == T(Inf) || abs(p.ρ) == T(Inf) || abs(med.ρ) == zero(T)
        numer = dj(m, ak)
        denom = dh1(m, ak)
    else
        q = (p.c*p.ρ)/(med.c*med.ρ) #the impedance
        if q == zero(T)
          numer =  j(m, ak)
          denom =  h1(m, ak)
        else
          γ = med.c/p.c #speed ratio
          numer = q * dj(m, ak)*j(m, γ*ak) - j(m, ak)*dj(m, γ*ak)
          denom = q * dh1(m, ak)*j(m, γ*ak) - h1(m, ak)*dj(m, γ*ak)
        end
    end

    return numer / denom
end


text:
@doc raw"""
        ub_ff(canonical_model::CanonicalModel,
              cons_name::Symbol,
              param_reference::UpdateRef,
              var_name::Symbol)

Constructs a parametrized upper bound constraint to implement feedforward from other models.
The Parameters are initialized using the uppper boundary values of the provided variables.

# Constraints
``` variable[var_name, t] <= param_reference[var_name] ```

# LaTeX

`` x \leq param^{max}``

# Arguments
* canonical_model::CanonicalModel : the canonical model built in PowerSimulations
* cons_name::Symbol : name of the constraint
* param_reference : Reference to the Parameter used to determine the upperbound
* var_name::Symbol : the name of the continuous variable
"""
function ub_ff(canonical_model::CanonicalModel,
               cons_name::Symbol,
               param_reference::UpdateRef,
               var_name::Symbol)

    time_steps = model_time_steps(canonical_model)
    ub_name = _middle_rename(cons_name, "_", "ub")
    variable = var(canonical_model, var_name)

    axes = JuMP.axes(variable)
    set_name = axes[1]

    @assert axes[2] == time_steps
    _add_param_container!(canonical_model, param_reference, set_name)
    param_ub = par(canonical_model, param_reference)

    _add_cons_container!(canonical_model, ub_name, set_name, time_steps)
    con_ub = con(canonical_model, ub_name)

    for name in axes[1]
        value = JuMP.upper_bound(variable[name, 1])
        param_ub[name] = PJ.add_parameter(canonical_model.JuMPmodel, value)
        for t in axes[2]
            con_ub[name, t] = JuMP.@constraint(canonical_model.JuMPmodel,
                                                variable[name, t] <= param_ub[name])
        end
    end

    return

end

@doc raw"""
        range_ff(canonical_model::CanonicalModel,
                        cons_name::Symbol,
                        param_reference::NTuple{2, UpdateRef},
                        var_name::Symbol)

Constructs min/max range parametrized constraint from device variable to include feedforward.

# Constraints

``` param_reference[1][var_name] <= variable[var_name, t] ```
``` variable[var_name, t] <= param_reference[2][var_name] ```

where r in range_data.

# LaTeX

`` param^{min} \leq x ``
`` x \leq param^{max}``

# Arguments
* canonical_model::CanonicalModel : the canonical model built in PowerSimulations
* param_reference::NTuple{2, UpdateRef} : Tuple with the lower bound and upper bound parameter reference
* cons_name::Symbol : name of the constraint
* var_name::Symbol : the name of the continuous variable
"""
function range_ff(canonical_model::CanonicalModel,
                  cons_name::Symbol,
                  param_reference::NTuple{2, UpdateRef},
                  var_name::Symbol)

    time_steps = model_time_steps(canonical_model)
    ub_name = _middle_rename(cons_name, "_", "ub")
    lb_name = _middle_rename(cons_name, "_", "lb")

    variable = var(canonical_model, var_name)
    axes = JuMP.axes(variable)
    set_name = axes[1]
    @assert axes[2] == time_steps

    #Create containers for the constraints
    _add_param_container!(canonical_model, param_reference[1], set_name)
    param_lb = par(canonical_model, param_reference[1])
    _add_param_container!(canonical_model, param_reference[2], set_name)
    param_ub = par(canonical_model, param_reference[2])

    #Create containers for the parameters
    _add_cons_container!(canonical_model, lb_name, set_name, time_steps)
    con_lb = con(canonical_model, lb_name)
    _add_cons_container!(canonical_model, ub_name, set_name, time_steps)
    con_ub = con(canonical_model, ub_name)

    for name in axes[1]
        param_lb[name] = PJ.add_parameter(canonical_model.JuMPmodel,
                                          JuMP.lower_bound(variable[name, 1]))
        param_ub[name] = PJ.add_parameter(canonical_model.JuMPmodel,
                                          JuMP.upper_bound(variable[name, 1]))
        for t in axes[2]
            con_ub[name, t] = JuMP.@constraint(canonical_model.JuMPmodel,
                                            variable[name, t] <= param_ub[name])
            con_lb[name, t] = JuMP.@constraint(canonical_model.JuMPmodel,
                                            variable[name, t] >= param_lb[name])
        end
    end

    return

end

@doc raw"""
            semicontinuousrange_ff(canonical_model::CanonicalModel,
                                    cons_name::Symbol,
                                    var_name::Symbol,
                                    param_reference::UpdateRef)

Constructs min/max range constraint from device variable with parameter setting.

# Constraints
If device min = 0:

``` variable[var_name, t] <= r[2].max*param_reference[var_name] ```

Otherwise:

``` variable[var_name, t] <= r[2].max*param_reference[var_name] ```

``` variable[var_name, t] >= r[2].min*param_reference[var_name] ```

where r in range_data.

# LaTeX

`` 0.0 \leq x^{var} \leq r^{max} x^{param}, \text{ for } r^{min} = 0 ``

`` r^{min} x^{param} \leq x^{var} \leq r^{min} x^{param}, \text{ otherwise } ``

# Arguments
* canonical_model::CanonicalModel : the canonical model built in PowerSimulations
* cons_name::Symbol : name of the constraint
* var_name::Symbol : the name of the continuous variable
* param_reference::UpdateRef : UpdateRef of the parameter
"""
function semicontinuousrange_ff(canonical_model::CanonicalModel,
                                cons_name::Symbol,
                                param_reference::UpdateRef,
                                var_name::Symbol)

    time_steps = model_time_steps(canonical_model)
    ub_name = _middle_rename(cons_name, "_", "ub")
    lb_name = _middle_rename(cons_name, "_", "lb")

    variable = var(canonical_model, var_name)

    axes = JuMP.axes(variable)
    set_name = axes[1]
    @assert axes[2] == time_steps

    _add_param_container!(canonical_model, param_reference, set_name)
    param = par(canonical_model, param_reference)

    _add_cons_container!(canonical_model, ub_name, set_name, time_steps)
    _add_cons_container!(canonical_model, lb_name, set_name, time_steps)
    con_ub = con(canonical_model, ub_name)
    con_lb = con(canonical_model, lb_name)

    for name in axes[1]
        ub_value = JuMP.upper_bound(variable[name, 1])
        lb_value = JuMP.lower_bound(variable[name, 1])
        param[name] = PJ.add_parameter(canonical_model.JuMPmodel, 1.0)
        for t in axes[2]
            con_ub[name, t] = JuMP.@constraint(canonical_model.JuMPmodel,
                                            variable[name, t] <= ub_value*param[name])
            con_lb[name, t] = JuMP.@constraint(canonical_model.JuMPmodel,
                                        variable[name, t] >= lb_value*param[name])
        end
    end

    # If the variable was a lower bound != 0, not removing the LB can cause infeasibilities
    for v in variable
        if JuMP.has_lower_bound(v)
            JuMP.set_lower_bound(v, 0.0)
        end
    end

    return

end

########################## FeedForward Constraints #########################################

function feedforward!(canonical_model::CanonicalModel,
                     device_type::Type{T},
                     ff_model::Nothing) where {T<:PSY.Component}
    return
end

function feedforward!(canonical_model::CanonicalModel,
                     device_type::Type{I},
                     ff_model::UpperBoundFF) where {I<:PSY.Injection}

    for prefix in get_vars_prefix(ff_model)
        var_name = Symbol(prefix, "_$(I)")
        parameter_ref = UpdateRef{JuMP.VariableRef}(var_name)
        ub_ff(canonical_model,
              Symbol("FF_$(I)"),
                     parameter_ref,
                     var_name)
    end

    return

end

function feedforward!(canonical_model::CanonicalModel,
                     device_type::Type{I},
                     ff_model::SemiContinuousFF) where {I<:PSY.Injection}

    bin_var = Symbol(get_bin_prefix(ff_model), "_$(I)")
    parameter_ref = UpdateRef{JuMP.VariableRef}(bin_var)
    for prefix in get_vars_prefix(ff_model)
        var_name = Symbol(prefix, "_$(I)")
        semicontinuousrange_ff(canonical_model,
                               Symbol("FFbin_$(I)"),
                               parameter_ref,
                               var_name)
    end

    return

end


text:
using Test
import NNlib
import KissConv

to_whcn(arr) = reshape(arr, (size(arr)..., 1, 1))

@testset "Against NNlib" begin
    for item in [
            (arr1=randn(10), arr2=randn(2)),
            (arr1=randn(Float32, 10), arr2=randn(Float32, 2)),
            (arr1=randn(Float32, 2,3), arr2=randn(Float32, 2,1)),
            (arr1=rand(1:10, 2,3), arr2=randn(Float32, 2,1)),
            (arr1=rand(1:100, 2,3,4), arr2=rand(1:100, 2,1, 3)),
        ]

        res_nn = NNlib.conv(
            to_whcn(item.arr1),
            to_whcn(item.arr2),
        )

        res_kiss = @inferred KissConv.conv(item.arr1, item.arr2)
        @test eltype(res_kiss) == eltype(res_nn)
        @test to_whcn(res_kiss) ≈ res_nn
    end
end


text:
module GraphBuilding

export build_graph
export weights_voronoi, weights_cutoff
export inverse_square, exp_decay

using PyCall
using ChemistryFeaturization
using Serialization

# options for decay of bond weights with distance...
# user can of course write their own as well
inverse_square(x) = x^-2.0
exp_decay(x) = exp(-x)

"""
Function to build graph from a file storing a crystal structure (currently supports anything ase.io.read can read in). Returns an AtomGraph object.

# Arguments
## Required Arguments
- `file_path::String`: Path to ASE-readable file containing a molecule/crystal structure

## Keyword Arguments
- `normalize_weights::Bool=true`: Whether to rescale graph weights such that the maximum value is 1.0 (recommended)
- `use_voronoi::bool`: if true, use Voronoi method for neighbor lists, if false use cutoff method

    (The rest of these parameters are only used if `use_voronoi==false`)

- `cutoff_radius::Real=8.0`: cutoff radius for atoms to be considered neighbors (in angstroms)
- `max_num_nbr::Integer=12`: maximum number of neighbors to include (even if more fall within cutoff radius)
- `dist_decay_func::Function=inverse_square`: function to determine falloff of graph edge weights with neighbor distance
"""
function build_graph(
    file_path::String;
    use_voronoi::Bool = false,
    cutoff_radius::Real = 8.0,
    max_num_nbr::Integer = 12,
    dist_decay_func::Function = inverse_square,
    normalize_weights::Bool = true,
)
    aseio = pyimport_conda("ase.io", "ase", "conda-forge")
    atoms_object = aseio.read(file_path)

    # list of atom symbols
    atom_ids = [get(atoms_object, i - 1).symbol for i = 1:length(atoms_object)]

    # check if any nonperiodic BC's
    nonpbc = any(.!atoms_object.pbc)
    local cant_voronoi = false
    if nonpbc & use_voronoi
        @warn "Voronoi edge weights are not supported if any direction in the structure is nonperiodic. Using cutoff weights method..."
        cant_voronoi = true
    end

    if use_voronoi && !cant_voronoi
        s = pyimport_conda("pymatgen.core.structure", "pymatgen", "conda-forge")
        pmgase = pyimport_conda("pymatgen.io.ase", "pymatgen", "conda-forge")
        aa = pmgase.AseAtomsAdaptor()
        struc = aa.get_structure(atoms_object)
        weight_mat = weights_voronoi(struc)
    else
        nl = pyimport_conda("ase.neighborlist", "ase", "conda-forge")
        is, js, dists = nl.neighbor_list("ijd", atoms_object, cutoff_radius)
        weight_mat = weights_cutoff(
            is .+ 1,
            js .+ 1,
            dists;
            max_num_nbr = max_num_nbr,
            dist_decay_func = dist_decay_func,
        )
    end

    if normalize_weights
        weight_mat = weight_mat ./ maximum(weight_mat)
    end

    return weight_mat, atom_ids
end

"""
Build graph using neighbor number cutoff method adapted from original CGCNN. Note that `max_num_nbr` is a "soft" max, in that if there are more of the same distance as the last, all of those will be added.
"""
function weights_cutoff(is, js, dists; max_num_nbr = 12, dist_decay_func = inverse_square)
    # sort by distance
    ijd = sort([t for t in zip(is, js, dists)], by = t -> t[3])

    # initialize neighbor counts
    num_atoms = maximum(is)
    local nb_counts = Dict(i => 0 for i = 1:num_atoms)
    local longest_dists = Dict(i => 0.0 for i = 1:num_atoms)

    # iterate over list of tuples to build edge weights...
    # note that neighbor list double counts so we only have to increment one counter per pair
    weight_mat = zeros(Float32, num_atoms, num_atoms)
    for (i, j, d) in ijd
        # if we're under the max OR if it's at the same distance as the previous one
        if nb_counts[i] < max_num_nbr || isapprox(longest_dists[i], d)
            weight_mat[i, j] += dist_decay_func(d)
            longest_dists[i] = d
            nb_counts[i] += 1
        end
    end

    # average across diagonal, just in case
    weight_mat = 0.5 .* (weight_mat .+ weight_mat')
end

"""
Build graph using neighbors from faces of Voronoi polyedra and weights from areas. Based on the approach from https://github.com/ulissigroup/uncertainty_benchmarking/blob/aabb407807e35b5fd6ad06b14b440609ae09e6ef/BNN/data_pyro.py#L268
"""
function weights_voronoi(struc)
    num_atoms = size(struc)[1]
    sa = pyimport_conda("pymatgen.analysis.structure_analyzer", "pymatgen", "conda-forge")
    vc = sa.VoronoiConnectivity(struc)
    conn = vc.connectivity_array
    weight_mat = zeros(Float32, num_atoms, num_atoms)
    # loop over central atoms
    for atom_ind = 1:size(conn)[1]
        # loop over neighbor atoms
        for nb_ind = 1:size(conn)[2]
            # loop over each possible PBC image for chosen image
            for image_ind = 1:size(conn)[3]
                # only add as neighbor if atom is not current center one AND there is connectivity to image
                if (atom_ind != image_ind) && (conn[atom_ind, nb_ind, image_ind] != 0)
                    weight_mat[atom_ind, nb_ind] +=
                        conn[atom_ind, nb_ind, image_ind] / maximum(conn[atom_ind, :, :])
                end
            end
        end
    end

    # average across diagonal (because neighborness isn't strictly symmetric in the way we're defining it here)
    weight_mat = 0.5 .* (weight_mat .+ weight_mat')
end

# TODO: graphs from SMILES via OpenSMILES.jl

end


text:
# function makeRadauIntegrator(N::Int64, tol::Float64, de_object::T_object, NC::Int64) where {T_object}
#     return RadauIntegrator{T_object, N, 3, NC}(tol, de_object)
# end
#
# function makeRadauIntegrator(x0::Vector{Float64}, tol::Float64, de_object::T_object, NC::Int64) where {T_object}
#     return makeRadauIntegrator(length(x0), tol, de_object, NC)
# end

function makeRadauIntegrator(de_object::T_object, x_or_N::Union{Int64,Vector{Float64}}, tol::Float64=1.0e-16,
    NR::Int64=2, NC::Int64=6) where {T_object}

    NX = ifelse(x_or_N isa Int64, x_or_N, length(x_or_N))
    NSM = radau_rule_to_stage(NR)
    return RadauIntegrator{NX,NC,NR,NSM}(tol, de_object)
end

function get_X_final(rr::RadauIntegrator{T_object, N, n_stage_max}, table::RadauTable{n_stage}) where {n_stage, n_stage_max, N, T_object}
    return rr.ct.X_stage[n_stage] * 1.0
end

@inline get_exponent(table::RadauTable{n_stage}) where {n_stage} = 1 / (1 + n_stage)
@inline get_h(rr::RadauIntegrator) = rr.step.h

@inline function get_table_from_current_s(rr::RadauIntegrator{T_object, N, n_stage_max}) where {n_stage_max, N, T_object}
    return rr.table[rr.rule.s]
end

function print_exit_flag(rr::RadauIntegrator)
    exit_flag = rr.step.exit_flag
    (exit_flag == -9999) && println("never run")
    (exit_flag == 1) && println("iteration limit exceeded")
    (exit_flag == 2) && println("tolerance failure")
    (exit_flag == 3) && println("newton divergence")
    is_exit_flag_success(rr) && println("step successful")
end

@inline set_exit_flag_success(rr::RadauIntegrator) = (rr.step.exit_flag = 0)
@inline set_exit_flag_fail_iter(rr::RadauIntegrator) = (rr.step.exit_flag = 1)
@inline set_exit_flag_fail_tol(rr::RadauIntegrator) = (rr.step.exit_flag = 2)
@inline set_exit_flag_fail_diverge(rr::RadauIntegrator) = (rr.step.exit_flag = 3)

@inline is_exit_flag_success(rr::RadauIntegrator) = (rr.step.exit_flag == 0)


text:
# This file was generated, do not modify it. # hide
datac = select!(datac, Not(:hospital_number));


text:
"""
    AbstractGeoArray <: DimensionalData.AbstractDimArray

Abstract supertype for objects that wrap an array (or location of an array)
and metadata about its contents. It may be memory or hold a `FileArray`, which
holds the filename, and is only opened when required.

`AbstractGeoArray`s inherit from [`AbstractDimArray`]($DDarraydocs)
from DimensionalData.jl. They can be indexed as regular Julia arrays or with
DimensionalData.jl [`Dimension`]($DDdimdocs)s. They will plot as a heatmap in
Plots.jl with correct coordinates and labels, even after slicing with
`getindex` or `view`. `getindex` on a `AbstractGeoArray` will always return
a memory-backed `GeoArray`.
"""
abstract type AbstractGeoArray{T,N,D,A} <: AbstractDimensionalArray{T,N,D,A} end

# Interface methods ###########################################################
"""
    missingval(x)

Returns the value representing missing data in the dataset
"""
function missingval end
missingval(x) = missing
missingval(A::AbstractGeoArray) = A.missingval

filename(A::AbstractGeoArray) = filename(data(A))

cleanreturn(A::AbstractGeoArray) = modify(cleanreturn, A)
cleanreturn(x) = x

isdisk(A::AbstractGeoArray) = parent(A) isa DiskArrays.AbstractDiskArray
isdisk(x) = false
ismem(x) = !isdisk(A)

function Base.:(==)(A::AbstractGeoArray{T,N}, B::AbstractGeoArray{T,N}) where {T,N}
    size(A) == size(B) && all(A .== B)
end

"""
    crs(x)

Get the projected coordinate reference system of a `Y` or `X` `Dimension`,
or of the `Y`/`X` dims of an `AbstractGeoArray`.

For [`Mapped`](@ref) mode this may be `nothing` as there may be no projected
coordinate reference system at all.
"""
function crs end
function crs(obj)
    if hasdim(obj, Y)
        crs(dims(obj, Y))
    elseif hasdim(obj, X)
        crs(dims(obj, X))
    else
        error("No Y or X dimension, crs not available")
    end
end
crs(dim::Dimension) = crs(mode(dim))

"""
    mappedcrs(x)

Get the mapped coordinate reference system for the `Y`/`X` dims of an array.

In [`Projected`](@ref) mode this is used to convert [`Selector`]($DDselectordocs)
values form the mappedcrs defined projection to the underlying projection, and to
show plot axes in the mapped projection.

In `Mapped` mode this is the coordinate reference system of the index values.
"""
function mappedcrs end
function mappedcrs(obj)
    if hasdim(obj, Y)
        mappedcrs(dims(obj, Y))
    elseif hasdim(obj, X)
        mappedcrs(dims(obj, X))
    else
        error("No Y or X dimension, mappedcrs not available")
    end
end
mappedcrs(dim::Dimension) = mappedcrs(mode(dim))

for f in (:mappedbounds, :projectedbounds, :mappedindex, :projectedindex)
    @eval ($f)(A::AbstractGeoArray, dims_) = ($f)(dims(A, dims_))
    @eval ($f)(A::AbstractGeoArray) = ($f)(dims(A))
end

# DimensionalData methods

DD.units(A::AbstractGeoArray) = get(metadata(A), :units, nothing)
# Rebuild all types of AbstractGeoArray as GeoArray
function DD.rebuild(
    A::AbstractGeoArray, data, dims::Tuple, refdims, name,
    metadata, missingval=missingval(A)
)
    GeoArray(data, dims, refdims, name, metadata, missingval)
end
function DD.rebuild(A::AbstractGeoArray;
    data=data(A), dims=dims(A), refdims=refdims(A), name=name(A),
    metadata=metadata(A), missingval=missingval(A)
)
    rebuild(A, data, dims, refdims, name, metadata, missingval)
end

function DD.DimTable(As::Tuple{<:AbstractGeoArray,Vararg{<:AbstractGeoArray}}...)
    DimTable(DimStack(map(read, As...)))
end

# DiskArrays methods

DiskArrays.eachchunk(A::AbstractGeoArray) = DiskArrays.eachchunk(parent(A))
DiskArrays.haschunks(A::AbstractGeoArray) = DiskArrays.haschunks(parent(A))
function DA.readblock!(A::AbstractGeoArray, dst, r::AbstractUnitRange...)
    DA.readblock!(parent(A), dst, r...)
end
function DA.writeblock!(A::AbstractGeoArray, src, r::AbstractUnitRange...)
    DA.writeblock!(parent(A), src, r...)
end

# Base methods

Base.parent(A::AbstractGeoArray) = data(A)

"""
    open(f, A::AbstractGeoArray; write=false)

`open` is used to open any `AbstractGeoArray` and do multiple operations
on it in a safe way. The `write` keyword opens the file in write mode so that it
can be altered on disk using e.g. a broadcast.

`f` is a method that accepts a single argument - an `GeoArray` object
which is just an `AbstractGeoArray` that holds an open disk-based object.
Often it will be a `do` block:

```julia
# A is an `GeoArray` wrapping the opened disk-based object.
open(GeoArray(filepath); write=true) do A
    mask!(A; to=maskfile)
    A[I...] .*= 2
    # ...  other things you need to do with the open file
end
```

By using a do block to open files we ensure they are always closed again
after we finish working with them.
"""
function Base.open(f::Function, A::AbstractGeoArray; kw...)
    # Open FileArray to expose the actual dataset object, even inside nested wrappers
    select = FileArray
    ignore = Union{Dict,Set,Base.MultiplicativeInverses.SignedMultiplicativeInverse}
    fa = Flatten.flatten(parent(A), select, ignore)
    if fa == ()
        f(GeoArray(parent(A), dims(A), refdims(A), name(A), metadata(A), missingval(A)))
    else
        open(fa[1]; kw...) do x
            # Rewrap the opened object where the FileArray was
            d = Flatten.reconstruct(parent(A), (x,), select, ignore)
            f(GeoArray(d, dims(A), refdims(A), name(A), metadata(A), missingval(A)))
        end
    end
end
Base.write(A::T) where T <: AbstractGeoArray = write(filename(A), A)

# Concrete implementation ######################################################

"""
    GeoArray <: AbsractGeoArray

    GeoArray(filepath::AbstractString, dims; kw...)
    GeoArray(A::AbstractArray{T,N}, dims; kw...)
    GeoArray(A::AbstractGeoArray; kw...)

A generic [`AbstractGeoArray`](@ref) for spatial/raster array data. It may hold
memory-backed arrays or [`FileArray`](@ref), that simply holds the `String` path
to an unopened file. This will only be opened lazily when it is indexed with `getindex`
or when `read(A)` is called. Broadcasting, taking a view, reversing and most other
methods _do not_ load data from disk: they are applied later, lazily.

# Keywords

- `data`: can replace the data in an `AbstractGeoArray`
- `dims`: `Tuple` of `Dimension`s for the array.
- `refdims`: `Tuple of` position `Dimension`s the array was sliced from, defaulting to `()`.
- `key`: `Symbol` key to desired layer in a multi-layer dataset, when a `filpath` is used.
- `name`: `Symbol` name for the array. `key` is used by default when a filepath `String` is pased in.
- `missingval`: value reprsenting missing data, normally detected form the file. Set manually
    when you know the value is not specified or is incorrect. This will *not* change any
    values in the raster, it simply assigns which value is treated as missing. To replace all of
    the missing values in the raster, use [`replace_missing`](@ref).
- `metadata`: `ArrayMetadata` object for the array, or `NoMetadata()`.
- `crs`: the coordinate reference system of  the objects `XDim`/`YDim` dimensions.
    Only set this if you know the detected crs is incrorrect, or it is not present in
    the file.
- `mappedcrs`: the mapped coordinate reference system of the objects `XDim`/`YDim` dimensions.
    for `Mapped` lookups these are the actual values of the index. For `Projected` lookups
    this can be used to index in eg. `EPSG(4326)` lat/lon values, having it converted automatically.
    Only set this if the detected `mappedcrs` in incorrect, or the file does not have a `mappedcrs`,
    e.g. a tiff.
"""
struct GeoArray{T,N,D<:Tuple,R<:Tuple,A<:AbstractArray{T,N},Na,Me,Mi} <: AbstractGeoArray{T,N,D,A}
    data::A
    dims::D
    refdims::R
    name::Na
    metadata::Me
    missingval::Mi
end
function GeoArray(A::AbstractArray, dims::Tuple;
    refdims=(), name=Symbol(""), metadata=NoMetadata(), missingval=missing,
)
    GeoArray(A, DD.formatdims(A, dims), refdims, name, metadata, missingval)
end
function GeoArray(A::AbstractArray{<:Any,1}, dims::Tuple{<:Dimension,<:Dimension,Vararg}; kw...)
    GeoArray(reshape(A, map(length, dims)), dims; kw...)
end
function GeoArray(table, dims::Tuple; name, kw...)
    cols = Tables.getcolumns(table)
    reshape(cols[name], map(length, dims))
    GeoStack(layers; name, kw...)
end
GeoArray(A::AbstractArray; dims, kw...) = GeoArray(A, dims; kw...)
function GeoArray(A::AbstractDimArray;
    data=data(A), dims=dims(A), refdims=refdims(A),
    name=name(A), metadata=metadata(A), missingval=missingval(A)
)
    GeoArray(data, dims, refdims, name, metadata, missingval)
end
function GeoArray(filename::AbstractString; key=nothing, kw...)
    _open(filename) do ds
        key = filekey(ds, key)
        GeoArray(ds, filename, key; kw...)
    end
end
function GeoArray(ds, filename::AbstractString, key=nothing;
    crs=nothing, mappedcrs=nothing, dims=nothing, refdims=(),
    name=Symbol(key isa Nothing ? "" : string(key)),
    metadata=metadata(ds), missingval=missingval(ds), write=false,
    source=_sourcetype(filename)
)
    crs = defaultcrs(source, crs)
    mappedcrs = defaultmappedcrs(source, mappedcrs)
    dims = dims isa Nothing ? DD.dims(ds, crs, mappedcrs) : dims
    data = FileArray(ds, filename; key, write)
    GeoArray(data, dims, refdims, name, metadata, missingval)
end

@propagate_inbounds function Base.setindex!(A::GeoArray, x, I::DD.StandardIndices...)
    setindex!(data(A), x, I...)
end

filekey(ds, key) = key
filekey(filename::String) = Symbol(splitext(basename(filename))[1])

# Precompile
precompile(GeoArray, (String,))

@deprecate geoarray(args...; kw...) GeoArray(args...; kw...)


text:
using Test
@isdefined(TestHelpers) || include("../../TestHelpers.jl")

TestHelpers.@timed_testset "generated_weights/" begin
    include("generated_weights/mfc+mfc+softmax.jl")
    include("generated_weights/conv+fc+softmax.jl")
end


text:
immutable Float120 <: AbstractFloat
    hi::Float64
    lo::Float64

    Float120(hi::Float64, lo::Float64) = new(hi,lo)
end

hi(x::Float120) = x.hi
lo(x::Float120) = x.lo

Float120(hi::Float64) = Float120(hi,0.0)

# fully clean a,b then set hi,lo
function F120{T<:Float64}(a::T, b::T) # hi,lo = eftAdd(a,b)
  hi = a + b
  t = hi - a
  lo = (a - (hi - t)) + (b - t)
  Float120(hi,lo)
end


text:
module HumanoidLCMSimTests

module SideTest
using Test, Random
using HumanoidLCMSim

if isdefined(Random, :seed!)
    using Random: seed!
else
    const seed! = srand
end

seed!(1)

@testset "side" begin
    @test -Sides.left == Sides.right
    @test -Sides.right == Sides.left

    sides = [rand(Side) for i = 1 : 10000];
    @test isapprox(count(side -> side == Sides.left, sides) / length(sides), 0.5; atol = 0.05)

    @test Sides.flipsign_if_right(2., Sides.left) == 2.
    @test Sides.flipsign_if_right(2., Sides.right) == -2.
end
end # module

module ControlTest
using Test
using Random
using HumanoidLCMSim
using ValkyrieRobot
using RigidBodyDynamics
using BotCoreLCMTypes
import HumanoidLCMSim: set!

@testset "robot_state_t -> MechanismState" begin
    mechanism = Valkyrie().mechanism
    remove_fixed_tree_joints!(mechanism)
    feet = Dict(Sides.left => findbody(mechanism, "leftFoot"), Sides.right => findbody(mechanism, "rightFoot"))
    hands = Dict(Sides.left => findbody(mechanism, "leftPalm"), Sides.right => findbody(mechanism, "rightPalm"))
    robot_info = HumanoidRobotInfo(mechanism, feet, hands, parse_actuators(mechanism, ValkyrieRobot.urdfpath()))
    msg = robot_state_t()
    for joint in tree_joints(mechanism)
        num_velocities(joint) == 1 || continue
        push!(msg.joint_name, string(joint))
    end
    msg.num_joints = length(msg.joint_name)
    resize!(msg)
    state = MechanismState(mechanism)
    rand!(state)
    result = DynamicsResult(mechanism)
    dynamics!(result, state)
    τprev = similar(velocity(state))
    τprev .= 0
    t = 1.0
    set!(msg, result, robot_info, τprev, t, state)
    state_back = MechanismState(mechanism)
    set!(state_back, msg, robot_info)

    # terrible precision due to float32:
    principal_value!(state)
    principal_value!(state_back)
    @test configuration(state) ≈ configuration(state_back) atol=1e-6
    @test velocity(state) ≈ velocity(state_back) atol=1e-6
end
end # module

module AtlasTest
using Test
using LinearAlgebra
using HumanoidLCMSim
using RigidBodyDynamics
using AtlasRobot

@testset "no controller" begin
    # just make sure we get all the way to simulating:
    @test_throws HumanoidLCMSim.NoCommandError AtlasSim.run(headless = true)
end

@testset "standing controller" begin
    # run controller
    AtlasControl.run(async=true)

    # run simulation
    sol = AtlasSim.run(controlΔt = 1 / 500, final_time = 5.0);

    # checks
    mechanism = AtlasRobot.mechanism(add_flat_ground=true)
    finalstate = MechanismState(mechanism)
    copyto!(finalstate, sol.u[end])
    @test sol.retcode == :Success
    @test center_of_mass(finalstate).v[3] ≈ 1.06 atol=5e-2
    @test norm(velocity(finalstate)) ≈ 0 atol=5e-3
end

end # module

end # module


text:
using QMTK
using QMTK.Consts.Pauli

h = @ham sum(Region{1}, σ₁⊗σ₂) + sum(Region{2}, σ₁⊗σ₂)

chain = Chain(Fixed, 4)

t = kronsum(bonds(chain, 1)) do i, j
    :(σ₁[$i] ⊗ σ₂[$j])
end

t += kronsum(bonds(chain, 2)) do i, j
    :(σ₁[$i] ⊗ σ₂[$j])
end

rhs = rand(Bit, 4)
for (val, lhs) in h(chain, rhs)
    @test t[Int(lhs)+1, Int(rhs)+1] == val
end

# merge factor into local matrix
# @ham 2 * sum(Region{1}, σ₁⊗σ₂) # => sum(Region{1}, 2 * σ₁⊗σ₂)

# mat = σ₁⊗σ₂
# h1 = LocalHamiltonian(Region{1}, mat)
# h2 = LocalHamiltonian(Region{2}, mat)

# FusedHamiltonian(h1, h2)


text:
struct Views
    viewstore::Vector{SequencingView}
end

Views() = Views(Vector{SequencingView}())

function Views(n::Integer)
    return Views(Vector{SequencingView}(undef, n))
end

@inline Base.getindex(vs::Views, i::Integer) = getindex(vs.viewstore, i)

function Base.getindex(vs::Views, is::Vector{<:Integer})
    newviews = Views(length(is))
    for (i, j) in enumerate(is)
        newviews[i] = vs[j]
    end
    return newviews
end

@inline Base.setindex!(vs::Views, sv::SequencingView, i::Integer) = setindex!(vs.viewstore, sv, i)

@inline function Base.iterate(vs::Views, state = 1)
    if state > length(vs)
        nothing
    else
        newstate = state + 1
        return vs[state], newstate
    end
end

@inline Base.length(vs::Views) = length(vs.viewstore)
@inline Base.eachindex(vs::Views) = Base.OneTo(length(vs))

Base.push!(x::Views, y::SequencingView) = push!(x.viewstore, y)

Base.resize!(x::Views, s::Integer) = resize!(x.viewstore, s)

Base.show(io::IO, x::Views) = show(io, "text/plain", x.viewstore)


text:
using FromFile
@from "Core.jl" import Node, Options

"""
    evalTreeArray(tree::Node, cX::AbstractMatrix{T}, options::Options)

Evaluate a binary tree (equation) over a given input data matrix. The
options contain all of the operators used. This function fuses doublets
and triplets of operations for lower memory usage.

# Returns

- `(output, complete)::Tuple{AbstractVector{T}, Bool}`: the result,
    which is a 1D array, as well as if the evaluation completed
    successfully (true/false). A `false` complete means an infinity
    or nan was encountered, and a large loss should be assigned
    to the equation.
"""
function evalTreeArray(tree::Node, cX::AbstractMatrix{T}, options::Options)::Tuple{AbstractVector{T}, Bool} where {T<:Real}
    if tree.degree == 0
        deg0_eval(tree, cX, options)
    elseif tree.degree == 1
        if tree.l.degree == 2 && tree.l.l.degree == 0 && tree.l.r.degree == 0
            deg1_l2_ll0_lr0_eval(tree, cX, Val(tree.op), Val(tree.l.op), options)
        else
            deg1_eval(tree, cX, Val(tree.op), options)
        end
    else
        if tree.l.degree == 0 && tree.r.degree == 0
            deg2_l0_r0_eval(tree, cX, Val(tree.op), options)
        elseif tree.l.degree == 0
            deg2_l0_eval(tree, cX, Val(tree.op), options)
        elseif tree.r.degree == 0
            deg2_r0_eval(tree, cX, Val(tree.op), options)
        else
            deg2_eval(tree, cX, Val(tree.op), options)
        end
    end
end

# Break whenever there is a nan or inf. Since so many equations give
# such numbers, it saves a lot of computation to just skip computation,
# and return a large loss!
macro break_on_check(val, flag)
    :(if isnan($(esc(val))) || !isfinite($(esc(val)))
          $(esc(flag)) = false
          break
    end)
end

macro return_on_false(flag, retval)
    :(if !$(esc(flag))
          return ($(esc(retval)), false)
    end)
end

function deg2_eval(tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, options::Options)::Tuple{AbstractVector{T}, Bool} where {T<:Real,op_idx}
    n = size(cX, 2)
    (cumulator, complete) = evalTreeArray(tree.l, cX, options)
    @return_on_false complete cumulator
    (array2, complete2) = evalTreeArray(tree.r, cX, options)
    @return_on_false complete2 cumulator
    op = options.binops[op_idx]
    finished_loop = true
    @inbounds @simd for j=1:n
        x = op(cumulator[j], array2[j])::T
        @break_on_check x finished_loop
        cumulator[j] = x
    end
    return (cumulator, finished_loop)
end

function deg1_eval(tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, options::Options)::Tuple{AbstractVector{T}, Bool} where {T<:Real,op_idx}
    n = size(cX, 2)
    (cumulator, complete) = evalTreeArray(tree.l, cX, options)
    @return_on_false complete cumulator
    op = options.unaops[op_idx]
    finished_loop = true
    @inbounds @simd for j=1:n
        x = op(cumulator[j])::T
        @break_on_check x finished_loop
        cumulator[j] = x
    end
    return (cumulator, finished_loop)
end

function deg0_eval(tree::Node, cX::AbstractMatrix{T}, options::Options)::Tuple{AbstractVector{T}, Bool} where {T<:Real}
    n = size(cX, 2)
    if tree.constant
        return (fill(convert(T, tree.val), n), true)
    else
        return (cX[tree.feature, :], true)
    end
end

function deg1_l2_ll0_lr0_eval(tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, ::Val{op_l_idx}, options::Options)::Tuple{AbstractVector{T}, Bool} where {T<:Real,op_idx,op_l_idx}
    n = size(cX, 2)
    op = options.unaops[op_idx]
    op_l = options.binops[op_l_idx]
    finished_loop = true
    if tree.l.l.constant && tree.l.r.constant
        val_ll = convert(T, tree.l.l.val)
        val_lr = convert(T, tree.l.r.val)
        x_l = op_l(val_ll, val_lr)::T
        if isnan(x_l) || !isfinite(x_l)
            return (Array{T, 1}(undef, n), false)
        end
        x = op(x_l)::T
        if isnan(x) || !isfinite(x)
            return (Array{T, 1}(undef, n), false)
        end
        return (fill(x, n), true)
    elseif tree.l.l.constant
        val_ll = convert(T, tree.l.l.val)
        feature_lr = tree.l.r.feature
        cumulator = Array{T, 1}(undef, n)
        finished_loop = true
        @inbounds @simd for j=1:n
            x_l = op_l(val_ll, cX[feature_lr, j])::T
            @break_on_check x_l finished_loop
            x = op(x_l)::T
            @break_on_check x finished_loop
            cumulator[j] = x
        end
        return (cumulator, finished_loop)
    elseif tree.l.r.constant
        feature_ll = tree.l.l.feature
        val_lr = convert(T, tree.l.r.val)
        cumulator = Array{T, 1}(undef, n)
        finished_loop = true
        @inbounds @simd for j=1:n
            x_l = op_l(cX[feature_ll, j], val_lr)::T
            @break_on_check x_l finished_loop
            x = op(x_l)::T
            @break_on_check x finished_loop
            cumulator[j] = x
        end
        return (cumulator, finished_loop)
    else
        feature_ll = tree.l.l.feature
        feature_lr = tree.l.r.feature
        cumulator = Array{T, 1}(undef, n)
        finished_loop = true
        @inbounds @simd for j=1:n
            x_l = op_l(cX[feature_ll, j], cX[feature_lr, j])::T
            @break_on_check x_l finished_loop
            x = op(x_l)::T
            @break_on_check x finished_loop
            cumulator[j] = x
        end
        return (cumulator, finished_loop)
    end
    return (cumulator, finished_loop)
end

function deg2_l0_r0_eval(tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, options::Options)::Tuple{AbstractVector{T}, Bool} where {T<:Real,op_idx}
    n = size(cX, 2)
    op = options.binops[op_idx]
    finished_loop = true
    if tree.l.constant && tree.r.constant
        val_l = convert(T, tree.l.val)
        val_r = convert(T, tree.r.val)
        x = op(val_l, val_r)::T
        if isnan(x) || !isfinite(x)
            return (Array{T, 1}(undef, n), false)
        end
        return (fill(x, n), true)
    elseif tree.l.constant
        cumulator = Array{T, 1}(undef, n)
        val_l = convert(T, tree.l.val)
        feature_r = tree.r.feature
        @inbounds @simd for j=1:n
            x = op(val_l, cX[feature_r, j])::T
            @break_on_check x finished_loop
            cumulator[j] = x
        end
    elseif tree.r.constant
        cumulator = Array{T, 1}(undef, n)
        feature_l = tree.l.feature
        val_r = convert(T, tree.r.val)
        @inbounds @simd for j=1:n
            x = op(cX[feature_l, j], val_r)::T
            @break_on_check x finished_loop
            cumulator[j] = x
        end
    else
        cumulator = Array{T, 1}(undef, n)
        feature_l = tree.l.feature
        feature_r = tree.r.feature
        @inbounds @simd for j=1:n
            x = op(cX[feature_l, j], cX[feature_r, j])::T
            @break_on_check x finished_loop
            cumulator[j] = x
        end
    end
    return (cumulator, finished_loop)
end

function deg2_l0_eval(tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, options::Options)::Tuple{AbstractVector{T}, Bool} where {T<:Real,op_idx}
    n = size(cX, 2)
    (cumulator, complete) = evalTreeArray(tree.r, cX, options)
    @return_on_false complete cumulator
    op = options.binops[op_idx]
    finished_loop = true
    if tree.l.constant
        val = convert(T, tree.l.val)
        @inbounds @simd for j=1:n
            x = op(val, cumulator[j])::T
            @break_on_check x finished_loop
            cumulator[j] = x
        end
    else
        feature = tree.l.feature
        @inbounds @simd for j=1:n
            x = op(cX[feature, j], cumulator[j])::T
            @break_on_check x finished_loop
            cumulator[j] = x
        end
    end
    return (cumulator, finished_loop)
end

function deg2_r0_eval(tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, options::Options)::Tuple{AbstractVector{T}, Bool} where {T<:Real,op_idx}
    n = size(cX, 2)
    (cumulator, complete) = evalTreeArray(tree.l, cX, options)
    @return_on_false complete cumulator
    op = options.binops[op_idx]
    finished_loop = true
    if tree.r.constant
        val = convert(T, tree.r.val)
        @inbounds @simd for j=1:n
            x = op(cumulator[j], val)::T
            @break_on_check x finished_loop
            cumulator[j] = x
        end
    else
        feature = tree.r.feature
        @inbounds @simd for j=1:n
            x = op(cumulator[j], cX[feature, j])::T
            @break_on_check x finished_loop
            cumulator[j] = x
        end
    end
    return (cumulator, finished_loop)
end

# Evaluate an equation over an array of datapoints
# This one is just for reference. The fused one should be faster.
function differentiableEvalTreeArray(tree::Node, cX::AbstractMatrix{T}, options::Options)::Tuple{AbstractVector{T}, Bool} where {T<:Real}
    n = size(cX, 2)
    if tree.degree == 0
        if tree.constant
            return (ones(T, n) .* tree.val, true)
        else
            return (cX[tree.feature, :], true)
        end
    elseif tree.degree == 1
        return deg1_diff_eval(tree, cX, Val(tree.op), options)
    else
        return deg2_diff_eval(tree, cX, Val(tree.op), options)
    end
end

function deg1_diff_eval(tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, options::Options)::Tuple{AbstractVector{T}, Bool} where {T<:Real,op_idx}
    n = size(cX, 2)
    (left, complete) = differentiableEvalTreeArray(tree.l, cX, options)
    @return_on_false complete left
    op = options.unaops[op_idx]
    out = op.(left)
    no_nans = !any(x -> (isnan(x) || !isfinite(x)), out)
    return (out, no_nans)
end

function deg2_diff_eval(tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, options::Options)::Tuple{AbstractVector{T}, Bool} where {T<:Real,op_idx}
    n = size(cX, 2)
    (left, complete) = differentiableEvalTreeArray(tree.l, cX, options)
    @return_on_false complete left
    (right, complete2) = differentiableEvalTreeArray(tree.r, cX, options)
    @return_on_false complete2 left
    op = options.binops[op_idx]
    out = op.(left, right)
    no_nans = !any(x -> (isnan(x) || !isfinite(x)), out)
    return (out, no_nans)
end


text:
module CUFFT

import CUDAapi

using ..CuArrays
using ..CuArrays: libcufft

import AbstractFFTs: plan_fft, plan_fft!, plan_bfft, plan_bfft!,
    plan_rfft, plan_brfft, plan_inv, normalization, fft, bfft, ifft, rfft,
    Plan, ScaledPlan
import Base: show, *, convert, unsafe_convert, size, strides, ndims
import Base.Sys: WORD_SIZE

using LinearAlgebra
import LinearAlgebra: mul!

include("libcufft_types.jl")
include("error.jl")

include("libcufft.jl")
include("genericfft.jl")
include("fft.jl")
include("wrappers.jl")
include("highlevel.jl")

version() = VersionNumber(cufftGetProperty(CUDAapi.MAJOR_VERSION),
                          cufftGetProperty(CUDAapi.MINOR_VERSION),
                          cufftGetProperty(CUDAapi.PATCH_LEVEL))

end


text:
export rotation2D

"""
getQ(xc)

builds basis for fast computation of rigid2D
"""
function getQ(xc::Array{Float64})
	n  = round(Int,length(xc)/2)
	xc = reshape(xc,n,2)
	return  kron(speye(2), [xc ones(n)])
end

"""
rotation2D(w,x,doDerivative,Q)

returns rotated points x and derivatives.
"""
function rotation2D(w::Float64,x::Array{Float64};doDerivative::Bool=false,Q::AbstractArray{Float64,2}=getQ(x),c::Vector{Float64}=zeros(2))

R  = [ cos(w) -sin(w);sin(w) cos(w)]
g  = (eye(2)-R)*vec(c)
f  = [R[1,1];R[1,2];g[1];R[2,1];R[2,2];g[2]]
y  = Q*f

if doDerivative
	dR = [-sin(w[1]) -cos(w[1]);cos(w[1]) -sin(w[1])]
	dg = -dR*vec(c)
	df = [dR[1,1];dR[1,2];dg[1];dR[2,1];dR[2,2];dg[2]]
	dy = Q*df
	return y,dy
else
	return y
end
end


text:
mutable struct Lattice{D,F,T,C, B}
    lattice_dims::NTuple{D,Int64}
    basis_vectors::Matrix{T}
    cell_vectors::Vector{Vector{T}}
    increments::Vector{Vector{T}}
    Js::NTuple{3,Float64}
    hs::NTuple{3,Float64}
    Jfunc::Interaction{F}
    cell_size::Int64
    tot_cell_num::Int64
    tot_spin_num::Int64
    basis_name::Symbol
    cell_name::Symbol
    ccell_cartesian::NTuple{D,Int64}
    ccell_linear::Int64
    complex::Val{C}
    function Lattice(lattice_dims::NTuple{D,Int64}, Js::NTuple{3,Float64}, Jfunc::Interaction{F}, hs::NTuple{3,Float64}, basis_vectors::Matrix{T}, cell_vectors::Vector{Vector{T}}, basis_name::Symbol, cell_name::Symbol, B::Bool = true) where {D,F,T}
        Dim = length(lattice_dims)
        assert(det(basis_vectors) != 0)
        if length(cell_vectors)==0
            push!(cell_vectors,zeros(T,Dim))
        end
        cell_size = length(cell_vectors)
        increments = Vector{Vector{Float64}}()
        for incr in product([[-1,0,1] for i in 1:Dim]...)
            increment = zeros(Float64, Dim)
            for i in 1:Dim
                @. increment += basis_vectors[:,i]*incr[i]*lattice_dims[i]
            end
            push!(increments, increment)
        end
        tot_cell_num = prod(lattice_dims)
        tot_spin_num = tot_cell_num*cell_size
        cc = hs[2]!=0.
        ccell_cartesian = div.(lattice_dims,2).+1
        ccell_linear = sub2ind(lattice_dims, ccell_cartesian...)
        new{D,F,T,cc,B}(lattice_dims, basis_vectors, cell_vectors, increments, Js, hs, Jfunc, cell_size, tot_cell_num, tot_spin_num, basis_name, cell_name, ccell_cartesian, ccell_linear, Val{cc}())
    end
end

function Lattice(dims, Js, Jfunc, hs, T::Type = Int64, B::Bool = true)
    basis_vectors = eye(T, length(dims), length(dims))
    cell_vectors = Vector{Vector{T}}()
    Lattice(dims, Js, Jfunc, hs, basis_vectors, cell_vectors, :cubic, :simple, B)
end
Lattice(dims, Js, Jfunc, hs, basis_vector::Matrix{T}, basis_name::Symbol, B::Bool = true) where {T} = Lattice(dims, Js, Jfunc, hs, basis_vector, Vector{Vector{T}}(), basis_name, :simple, B)

function get_vector(L::Lattice, spin::I) where {I}
    L.basis_vectors*collect(ind2sub(L.lattice_dims, spin))
end
function get_vector(L::Lattice, spin::Tuple{I, I}) where {I}
    L.basis_vectors*collect(ind2sub(L.lattice_dims, spin[1])) + L.cell_vectors[spin[2]]
end

function check_spin(L::Lattice, spin::I) where {I}
    assert(0<spin<=L.tot_spin_num)
end
function check_spin(L::Lattice, spin::Tuple{I, I}) where {I}
    assert((0<spin[1]<=prod(L.lattice_dims)) && (0<spin[2]<=L.cell_size))
end

function get_patch_vector(L::Lattice, spin1::S, spin2::S) where {S<:Tuple{Integer,Integer}}
    check_spin(L, spin1)
    check_spin(L, spin2)
    new_s2_cart = ind2sub(L.lattice_dims, spin2[1]) .- ind2sub(L.lattice_dims, spin1[1]) .+ L.ccell_cartesian
    new_s2_cart = mod.(new_s2_cart.-1, L.lattice_dims).+1
    return get_vector(L, (sub2ind(L.lattice_dims, new_s2_cart...), spin2[2])) .- get_vector(L, (L.ccell_linear, spin1[2]))
end

function get_patch_vector(L1::Lattice, L2::Lattice, spin1::S, spin2::S) where {S<:Tuple{Integer,Integer}}
    check_spin(L1, spin1)
    check_spin(L2, spin2)
    new_s2_cart = ind2sub(L1.lattice_dims, spin2[1]) .- ind2sub(L1.lattice_dims, spin1[1]) .+ L1.ccell_cartesian
    new_s2_cart = mod.(new_s2_cart.-1, L1.lattice_dims).+1
    return get_vector(L2, (sub2ind(L2.lattice_dims, new_s2_cart...), spin2[2])) .- get_vector(L1, (L1.ccell_linear, spin1[2]))
end

@inline get_patch_vector(L::Lattice, spin1::S, spin2::S) where {S<:Integer} = get_patch_vector(L, (spin1, 1), (spin2, 1))

@inline get_patch_vector(L1::Lattice, L2::Lattice, spin1::S, spin2::S) where {S<:Integer} = get_patch_vector(L1, L2, (spin1, 1), (spin2, 1))

function get_simple_vector(L::Lattice, spin1::S, spin2::S) where {S<:Union{Integer,Tuple{Integer,Integer}}}
    check_spin(L, spin1)
    check_spin(L, spin2)
    return get_vector(L, spin1) .- get_vector(L, spin2)
end

function get_smallest_vector(L::Lattice, spin1::S, spin2::S) where {S<:Union{Integer,Tuple{Integer,Integer}}}
    check_spin(L, spin1)
    check_spin(L, spin2)
    vector = get_vector(L, spin2) .- get_vector(L, spin1)
    smallest = deepcopy(vector)
    temp_vector = deepcopy(vector)
    r = norm(vector)
    for incr in L.increments
        temp_vector .= vector .+ incr
        temp_r = norm(temp_vector)
        if temp_r<r
            smallest = deepcopy(temp_vector)
            r = temp_r
        end
    end
    return smallest
end

function get_value(L::Lattice{D,F,T,C,true}, spin1::S, spin2::S) where {D, F, T, C, S<:Union{Integer,Tuple{Integer,Integer}}}
    return L.Jfunc(get_patch_vector(L, spin1, spin2))
end

function get_value(L::Lattice{D,F,T,C,false}, spin1::S, spin2::S) where {D, F, T, C, S<:Union{Integer,Tuple{Integer,Integer}}}
    return L.Jfunc(get_simple_vector(L, spin1, spin2))
end

function get_value(L1::Lattice, L2::Lattice, spin1::S, spin2::S) where {S<:Union{Integer,Tuple{Integer,Integer}}}
    return L1.Jfunc(get_patch_vector(L1, L2, spin1, spin2))
end

function get_string(L::Ltype) where {Ltype<:Lattice}
    @sprintf("D%d d%d S:%s C:%s Js(%.3f,%.3f,%.3f) Jf%s H(%.3f,%.3f,%.3f)", length(L.lattice_dims), L.cell_size, string(L.basis_name), string(L.cell_name), L.Js[1], L.Js[2], L.Js[3], func_type(L.Jfunc), L.hs[1], L.hs[2], L.hs[3])
end

function get_string(vec::Vector)
    assert(length(vec)>0)
    str = "("
    i = 1
    str *= "$(vec[i])"
    i += 1
    while i<= length(vec)
        str *= ", $(vec[i])"
        i += 1
    end
    str *= ")"
    return str
end

function translate_indices(L::Lattice, indices::Vector{Tuple{Int64, Int64}})
    sort!(indices)
    return indices
end
function translate_indices(L::Lattice, indices::AbstractVector{Int64})
    nIndices = Vector{Tuple{Int64,Int64}}()
    for index in indices for c in 1:L.cell_size
        push!(nIndices, (index, c))
    end end
    sort!(nIndices)
    return nIndices
end
function translate_indices(L::Lattice{D,F,T,C}, indices::CartesianRange{CartesianIndex{D2}}) where {D,F,T,C,D2}
    assert(length(first(indices).I)==length(L.lattice_dims)+1)
    nIndices = Vector{NTuple{2, Int64}}()
    for index in indices
        push!(nIndices, (sub2ind(L.lattice_dims, index.I[1:D]...), index.I[end]))
    end
    sort!(nIndices)
    return nIndices
end

function read_indices(L::Lattice{D,F,T,C}, indices::Vector{NTuple{D, Int64}}) where {D,F,T,C}
    fff(index::NTuple{D}) = sub2ind(L.lattice_dims, index...)
    return translate_indices(L, map(fff, indices))
end

get_all_spins(L::Lattice) = NTuple{2,Int64}[(l,c) for l in Base.OneTo(L.tot_cell_num) for c in Base.OneTo(L.cell_size)]

function get_central_spins(L::Lattice, dims::NTuple{N,Int64}) where {N}
    v = Vector{Vector{Int64}}()
    for i in 1:N
        s = div(dims[i],2)
        if dims[i]%2==0
            push!(v, [s,s+1])
        else
            push!(v, [s+1])
        end
    end
    c_cells = product(v...)
    return translate_indices(L, [sub2ind(L.lattice_dims, idxs...) for idxs in product(v...)])
end

change_dims(dims1, dims2, index) = ind2sub(dims2, sub2ind(dims1, index...))

##########################3

struct SpinArray
    name::Symbol
    spins::Vector{NTuple{2,Int64}}
end
@inline SpinArray(spins::Vector{NTuple{2,Int64}}, name = :undef) = SpinArray(name, spins)
@inline SpinArray(L::Lattice{D}, spins::Vector{NTuple{D,Int64}}, name = :undef) where {D} = SpinArray(name, read_indices(L, spins))
function SpinArray(L::Lattice{D}, spins::Vector{T}, name = :undef) where T<:Tuple{NTuple{D, Int64}, Int64} where D
    fff(index::T) = (sub2ind(L.lattice_dims, index[1]...), index[2])
    return SpinArray(name, map(fff, spins))
end

@inline get_string(s::Symbol) = string(s)
@inline get_string(s::SpinArray) = s.name==:undef ? "$(s.spins)" : string(s.name)
get_string(link::Tuple{T1,T2}) where {T1,T2} = @sprintf("%s_%s", get_string(link[1]), get_string(link[2]))

Base.isless(s1::SpinArray, s2::Symbol) = Base.isless(s1.name, s2)
Base.isless(s1::Symbol, s2::SpinArray) = Base.isless(s1, s2.name)
Base.isless(s1::SpinArray, s2::SpinArray) = Base.isless(s1.name, s2.name)


text:
using Test, AdmittanceModels, LinearAlgebra
import .PSOModelExample, .LRCExample, .CircuitCascadeExample
const ae = PSOModelExample
const lrce = LRCExample
const cce = CircuitCascadeExample

@testset "PSOModel to Blackbox" begin
    # test 1
    ω = collect(range(4, stop=6, length=1000)) * 1e9 * 2π
    bbox = Blackbox(ω, ae.pso)
    @test get_ports(bbox) == get_ports(ae.pso)
    @test length(get_Y(bbox)) == length(ω)
    @test size(get_Y(bbox)[1]) == size(ae.pso.K)
    bbox = canonical_gauge(bbox)
    @test size(get_Y(bbox)[1]) == (2,2)
    # test 2
    l, r, c, ω₀ = lrce.l, lrce.r, lrce.c, lrce.ω
    ω = collect(range(.8, stop=1.2, length=1000)) * ω₀
    bbox = Blackbox(ω, lrce.pso)
    z = 1 ./(1 ./(l * 1im * ω) .+ 1/r .+ c * 1im * ω)
    @test length(get_Y(bbox)) == length(ω)
    impedance = impedance_matrices(bbox)
    @test all([size(x) == (1,1) for x in impedance])
    @test [x[1,1] for x in impedance] ≈ z
end

@testset "cascade, unite_ports, open_ports, short_ports" begin
    # cascade
    pso0 = PSOModel(cce.circuit0, [(0,1), (0,2)], [1, 2])
    pso1 = PSOModel(cce.circuit1, [(4,5), (4,6)], [5, 6])
    pso2 = PSOModel(cce.circuit2, [(8,9), (8,10)], [9, 10])
    pso = cascade(pso0, pso1, pso2)
    ω = collect(range(.1, stop=100, length=1000)) * 2π
    bbox_from_pso = Blackbox(ω, pso)
    @test get_ports(bbox_from_pso) == get_ports(pso)
    bboxes0 = [Blackbox(ω, a) for a in [pso0, pso1, pso2]]
    bbox0 = cascade(bboxes0)
    @test bbox0 ≈ bbox_from_pso
    bboxes1 = [canonical_gauge(b) for b in bboxes0]
    bbox1 = cascade(bboxes1)
    @test bbox1 ≈ canonical_gauge(bbox_from_pso)
    # unite_ports
    united_pso = unite_ports(pso, 1, 5, 9)
    bbox_from_pso = Blackbox(ω, united_pso)
    @test get_ports(bbox_from_pso) == get_ports(united_pso)
    united_bbox0 = unite_ports(bbox0, 1, 5, 9)
    @test united_bbox0 ≈ bbox_from_pso
    united_bbox1 = unite_ports(bbox1, 1, 5, 9)
    @test canonical_gauge(united_bbox1) ≈ canonical_gauge(bbox_from_pso)
    # open_ports
    open_pso = open_ports(united_pso, 2, 6)
    bbox_from_pso = Blackbox(ω, open_pso)
    open_bbox = open_ports(united_bbox1, 2, 6)
    @test canonical_gauge(open_bbox) ≈ canonical_gauge(bbox_from_pso)
    # short_ports
    short_pso = short_ports(united_pso, 2, 6)
    bbox_from_pso = Blackbox(ω, short_pso)
    short_bbox = short_ports(united_bbox1, 2, 6)
    @test canonical_gauge(short_bbox) ≈ canonical_gauge(bbox_from_pso)
end

# Pozar p 192
@testset "transfer function conversions" begin
    a, b, c, d = 1.0, 2.0, 3.0, 4.0
    Z0 = 10.0
    Y = [[d, -1] [b * c - a * d, a]] ./ b
    Z = [[a, 1] [a * d - b * c, d]] ./ c
    denom = a + b/Z0 + c * Z0 + d
    S = [[a + b/Z0 - c * Z0 - d, 2] [2 * (a * d - b * c), -a + b/Z0 - c * Z0 + d]] ./ denom
    @test inv(Y) ≈ Z
    Z0s = [Z0, Z0]
    @test AdmittanceModels.admittance_to_scattering(Y, Z0s) ≈ S
    @test AdmittanceModels.scattering_to_admittance(S, Z0s) ≈ Y
    @test AdmittanceModels.impedance_to_scattering(Z, Z0s) ≈ S
    @test AdmittanceModels.scattering_to_impedance(S, Z0s) ≈ Z
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module Simulate_Pendulum

Demonstrate how a pendulum can be defined, simulated and plotted
with ModiaMath.
"""
module Simulate_Pendulum

using ModiaMath
using ModiaMath.StaticArrays   # included via ModiaMath, to avoid requirement to add it in the standard environment

@component Pendulum(;L=1.0, m=1.0, d=0.1, g=9.81, phi0_deg=90.0) begin
    @assert(L > 0.0)
    @assert(m > 0.0)
    @assert(d >= 0.0)

    phi = RealScalar(start=phi0_deg*pi/180, unit="rad"    , fixed=true, nominal=1.0, info="Relative rotation angle",                     numericType=ModiaMath.XD_EXP)
    w   = RealScalar(start=0.0            , unit="rad/s"  , fixed=true, nominal=1.0, info="Relative angular velocity",     integral=phi, numericType=ModiaMath.XD_EXP)
    a   = RealScalar(                       unit="rad/s^2",                          info="Relative angular acceleration", integral=w  , numericType=ModiaMath.DER_XD_EXP)
    r   = RealSVector{2}(                   unit="m"      ,                          info="Tip position of pendulum"     ,               numericType=ModiaMath.WC)
end

function ModiaMath.computeVariables!(p::Pendulum, sim::ModiaMath.SimulationState)
    L = p.L
    m = p.m
    d = p.d
    g = p.g
    phi = p.phi.value
    w   = p.w.value

    p.a.value = (-m*g*L*sin(phi) - d*w) / (m*L^2)

    if ModiaMath.isStoreResult(sim)
        p.r.value = @SVector [L*sin(phi), -L*cos(phi)]
    end
end

simulationModel = ModiaMath.SimulationModel(Pendulum(L=0.8, m=0.5, d=0.2); stopTime=5.0, structureOfDAE=ModiaMath.DAE_ExplicitDerivatives)
result = ModiaMath.simulate!(simulationModel; log=true)

plot(result, [(:phi, :w) :a])
#plot(result, "r[2]", xAxis="r[1]", figure=2)

# ModiaMath.print_ModelVariables(simulationModel)
# println("result variables = ", result)

end


text:
using Test
using GaussQuadrature: legendre, both
using ClimateMachine.VTK: writemesh_highorder, writemesh_raw

@testset "VTK" begin
    for writemesh in (writemesh_highorder, writemesh_raw)
        Ns = writemesh == writemesh_raw ? ((4, 4, 4), (3, 5, 7)) : ((4, 4, 4),)
        for N in Ns
            for dim in 1:3
                nelem = 3
                T = Float64

                Nq = N .+ 1
                r = legendre(T, Nq[1], both)[1]
                s = dim < 2 ? [0] : legendre(T, Nq[2], both)[1]
                t = dim < 3 ? [0] : legendre(T, Nq[3], both)[1]
                Nq1 = length(r)
                Nq2 = length(s)
                Nq3 = length(t)
                Np = Nq1 * Nq2 * Nq3

                x1 = Array{T, 4}(undef, Nq1, Nq2, Nq3, nelem)
                x2 = Array{T, 4}(undef, Nq1, Nq2, Nq3, nelem)
                x3 = Array{T, 4}(undef, Nq1, Nq2, Nq3, nelem)

                for e in 1:nelem, k in 1:Nq3, j in 1:Nq2, i in 1:Nq1
                    xoffset = nelem + 1 - 2e
                    x1[i, j, k, e], x2[i, j, k, e], x3[i, j, k, e] =
                        r[i] - xoffset, s[j], t[k]
                end

                if dim == 1
                    x1 = x1 .^ 3
                elseif dim == 2
                    x1, x2 = x1 + sin.(π * x2) / 5, x2 + exp.(-x1 .^ 2)
                else
                    x1, x2, x3 = x1 + sin.(π * x2) / 5,
                    x2 + exp.(-hypot.(x1, x3) .^ 2),
                    x3 + sin.(π * x1) / 5
                end
                d = exp.(sin.(hypot.(x1, x2, x3)))
                s = copy(d)

                if dim == 1
                    @test "test$(dim)d.vtu" == writemesh(
                        "test$(dim)d",
                        x1;
                        fields = (("d", d), ("s", s)),
                    )[1]
                    @test "test$(dim)d.vtu" == writemesh(
                        "test$(dim)d",
                        x1;
                        x2 = x2,
                        fields = (("d", d), ("s", s)),
                    )[1]
                    @test "test$(dim)d.vtu" == writemesh(
                        "test$(dim)d",
                        x1;
                        x2 = x2,
                        x3 = x3,
                        fields = (("d", d), ("s", s)),
                    )[1]
                elseif dim == 2
                    @test "test$(dim)d.vtu" == writemesh(
                        "test$(dim)d",
                        x1,
                        x2;
                        fields = (("d", d), ("s", s)),
                    )[1]
                    @test "test$(dim)d.vtu" == writemesh(
                        "test$(dim)d",
                        x1,
                        x2;
                        x3 = x3,
                        fields = (("d", d), ("s", s)),
                    )[1]
                elseif dim == 3
                    @test "test$(dim)d.vtu" == writemesh(
                        "test$(dim)d",
                        x1,
                        x2,
                        x3;
                        fields = (("d", d), ("s", s)),
                    )[1]
                end
            end
        end
    end
end

using MPI
MPI.Initialized() || MPI.Init()
using ClimateMachine.Mesh.Topologies: BrickTopology
using ClimateMachine.Mesh.Grids: DiscontinuousSpectralElementGrid
using ClimateMachine.VTK: writevtk_helper

let
    mpicomm = MPI.COMM_SELF
    for FT in (Float64,) #Float32)
        for dim in 2:3
            for _N in ((2, 3, 4), (0, 2, 5), (3, 0, 0), (0, 0, 0))
                N = _N[1:dim]
                if dim == 2
                    Ne = (4, 5)
                    brickrange = (
                        range(FT(0); length = Ne[1] + 1, stop = 1),
                        range(FT(0); length = Ne[2] + 1, stop = 1),
                    )
                    topl = BrickTopology(
                        mpicomm,
                        brickrange,
                        periodicity = (false, false),
                    )
                    warpfun =
                        (x1, x2, _) -> begin
                            (x1 + sin(x1 * x2), x2 + sin(2 * x1 * x2), 0)
                        end
                elseif dim == 3
                    Ne = (3, 4, 5)
                    brickrange = (
                        range(FT(0); length = Ne[1] + 1, stop = 1),
                        range(FT(0); length = Ne[2] + 1, stop = 1),
                        range(FT(0); length = Ne[3] + 1, stop = 1),
                    )
                    topl = BrickTopology(
                        mpicomm,
                        brickrange,
                        periodicity = (false, false, false),
                    )
                    warpfun =
                        (x1, x2, x3) -> begin
                            (
                                x1 + (x1 - 1 / 2) * cos(2 * π * x2 * x3) / 4,
                                x2 + exp(sin(2π * (x1 * x2 + x3))) / 20,
                                x3 + x1 / 4 + x2^2 / 2 + sin(x1 * x2 * x3),
                            )
                        end
                end
                grid = DiscontinuousSpectralElementGrid(
                    topl,
                    FloatType = FT,
                    DeviceArray = Array,
                    polynomialorder = N,
                    meshwarp = warpfun,
                )
                Q = rand(FT, prod(N .+ 1), 3, prod(Ne))
                prefix = "test$(dim)d_raw$(prod(ntuple(i->"_$(N[i])", dim)))"
                @test "$(prefix).vtu" == writevtk_helper(
                    prefix,
                    grid.vgeo,
                    Q,
                    grid,
                    ("a", "b", "c");
                    number_sample_points = 0,
                )[1]
                prefix = "test$(dim)d_high_order$(prod(ntuple(i->"_$(N[i])", dim)))"
                @test "$(prefix).vtu" == writevtk_helper(
                    prefix,
                    grid.vgeo,
                    Q,
                    grid,
                    ("a", "b", "c");
                    number_sample_points = 10,
                )[1]
            end
        end
    end
end

max_forks_repo_name: jm-c/CLIMA


text:
using TetrahedronIntegration
using Test

@testset "TetrahedronIntegration.jl" begin
    include("test_tetrahedron.jl")
    include("test_gaussian.jl")
end

max_stars_repo_name: jaemolihm/TetrahedronIntegration

max_issues_repo_name: jaemolihm/TetrahedronIntegration

max_forks_repo_name: jaemolihm/TetrahedronIntegration


text:
using Snappy
##### set a path
path = "c:/scratch/test.parquet"

##### load metadata
metadata(path)

#### continues to try

io = open(path)

# the dicvtionary page starts here at 4
seek(io, 4); dict_page_header = read_thrift(io, PAR2.PageHeader)

compressed_length = dict_page_header.compressed_page_size

data_compressed = read(io, compressed_length)
data_uncompressed = Snappy.uncompress(b)

data_actual = reinterpret(Int64, data_uncompressed)

next_page = read_thrift(io, PAR2.PageHeader)

Parquet.coltype(page.colchunk)

io = IOBuffer(data_unc)
using Parquet
@which Parquet.read_plain_values(io, num_values, ctype)

typ = ctype
@which Parquet.read_plain(io, typ)

Parquet.PLAIN_JTYPES[typ+1]

@which Parquet.read_fixed(io, Int64)

ret = convert(UInt64, 0)
N = 8
T = UInt64
for n in 0:(N-1)
    global ret
    byte = convert(T, read(io, UInt8))
    ret |= (byte << *(8,n))
end
reinterpret(Int64, ret)

res = [read(io, Int64) for i in 1:389490]

Int(length(data_unc)//8)

pos= position(io)
read_thrift(io, PAR2.PageHeader)

data = read(io, 34966)
data_unc = Snappy.uncompress(data)

reinterpret(Int32, data_unc)

seek(io, 4)
read_thrift(io, PAR2.PageHeader)

num_values
read_plain_values(io, num_values, 6)

using StatsBase
countmap([pg.hdr._type for pg in pgs])

Debugger.@enter read_column(path, 1)


text:
@eval mutable struct Counter
    $((:($(Symbol(op[1], typ[2])) ::Int) for op in ops for typ in typs)...)
    Counter() = new($((0 for _ in 1:length(ops)*length(typs))...))
end

# Relatively inefficient, but there should be no need for performance here...

function flop(c::Counter)
    total = 0
    for (typ, suffix) in typs
        for (name, op, cnt) in ops
            fieldname = Symbol(name, suffix)
            total += cnt * getfield(c, fieldname)
        end
    end
    total
end

import Base: ==, *, show

function Base.show(io::IO, c::Counter)
    type_names  = [typ    for (typ, _)    in typs]
    type_suffix = [suffix for (_, suffix) in typs]
    op_names    = [name   for (name, _)   in ops]

    mat = [getfield(c, Symbol(name, suffix)) for
           name   in op_names,
           suffix in type_suffix]

    fl = flop(c)
    print(io, "Flop Counter: $fl flop")
    fl == 0 && return

    print(io, "\n")
    fc(data, i) = any(data[:,i] .> 0)
    fr(data, i) = any(data[i,:] .> 0)
    pretty_table(io, mat,
                 header = type_names,
                 row_names = op_names,
                 filters_col = (fc,),
                 filters_row = (fr,),
                 newline_at_end = false)
end

function ==(c1::Counter, c2::Counter)
    all(getfield(c1, field)==getfield(c2, field) for field in fieldnames(Counter))
end

function *(n::Int, c::Counter)
    ret = Counter()
    for field in fieldnames(Counter)
        setfield!(ret, field, n*getfield(c, field))
    end
    ret
end


text:
#returns squareroot of sum of squares, i.e. positive number
function  z(x::Matrix{Int64}, q::Vector{Int64}, p::Vector{Int64})
   # S is the number of students, N is the number of courses
   # x is the SxN allocation matrix
   # q is the Nx1 course capacity vector
   # p is the Nx1 price vector
   N = length(q)
   zer = zeros(N)        #Create empty Nx1 vector of clearing errors
   numstud = transpose(sum(x, dims = 1)) #Create Nx1 vector with the number of students enrolled in each course.
   for i in 1:N        #For every course
       if p[i] > 0     #As long as the price in this course is positive.
           zer[i] = numstud[i] - q[i]
       else
           zer[i] = max(numstud[i] - q[i], 0)
       end
   end
   α = sqrt(sum(zer.^2))
   return α
end

export z


text:
using FactCheck

facts("sNES") do

context("mix_with_indices") do
  cs = randn(2,3)
  r = BlackBoxOptim.mix_with_indices(cs)
  @fact length(r) => 2
  @fact r[1][2] => 1
  @fact r[2][2] => 2
  @fact r[1][1] => cs[1,:]
  @fact r[2][1] => cs[2,:]
end

assign_weights(rankedCandidates) = begin
  u = BlackBoxOptim.fitness_shaping_utilities_linear(length(rankedCandidates))
  BlackBoxOptim.assign_weights(rankedCandidates, u)
end

context("assign_weights") do

  context("when indices are already ordered") do
    u = assign_weights([(:dummy, 1.0), (:dummy, 2.0)])
    @fact length(u) => 2
    @fact u[1] => 1.0
    @fact u[2] => 0.0

    u = assign_weights([(:d, 1.0), (:d, 2.0), (:d, 3.0)])
    @fact length(u) => 3
    @fact u[1] => 1.0
    @fact u[2] => 0.0
    @fact u[3] => 0.0

    u = assign_weights([(:d, 1), (:d, 2), (:d, 3), (:d, 4)])
    @fact length(u) => 4
    @fact isapprox(u[1], 2/3) => true
    @fact isapprox(u[2], 1/3) => true
    @fact u[3] => 0.0
    @fact u[4] => 0.0

    u = assign_weights([(:d, 1), (:d, 2), (:d, 3), (:d, 4), (:d, 5)])
    @fact length(u) => 5
    @fact isapprox(u[1], 2/3) => true
    @fact isapprox(u[2], 1/3) => true
    @fact u[3] => 0.0
    @fact u[4] => 0.0
    @fact u[5] => 0.0

    u = assign_weights([(:d, 1), (:d, 2), (:d, 3), (:d, 4), (:d, 5), (:d, 6)])
    @fact length(u) => 6
    @fact isapprox(u[1], (((3/3)/(3/3+2/3+1/3)))) => true
    @fact isapprox(u[2], (((2/3)/(3/3+2/3+1/3)))) => true
    @fact isapprox(u[3], (((1/3)/(3/3+2/3+1/3)))) => true
    @fact u[4] => 0.0
    @fact u[5] => 0.0
    @fact u[6] => 0.0
  end

  context("when indices are not ordered") do
    u = assign_weights([(:d, 4), (:d, 1), (:d, 2), (:d, 3)])
    @fact isapprox(u[4], 2/3) => true
    @fact isapprox(u[1], 1/3) => true
    @fact u[2] => 0.0
    @fact u[3] => 0.0
  end
end

end


text:
with_testdir() do testdir
    registry_dir = joinpath(testdir, "TestRegistry")
    upstream_dir = joinpath(dirname(registry_dir), "upstream")
    upstream_url = "file:///$(upstream_dir)"
    create_registry(registry_dir, upstream_url,
                    description = "For testing purposes only.",
                    uuid = "ed6ca2f6-392d-11ea-3224-d3daf7fee369",
                    gitconfig = TEST_GITCONFIG, push = false)

    reg_path, is_temp = check_git_registry(registry_dir, TEST_GITCONFIG)
    @test reg_path == registry_dir
    @test !is_temp

    mkpath(upstream_dir)
    upstream_git = gitcmd(upstream_dir, TEST_GITCONFIG)
    run(`$(upstream_git) clone -q --bare $(registry_dir) .`)

    reg_path, is_temp = check_git_registry(upstream_url, TEST_GITCONFIG)
    @test isdir(joinpath(reg_path, ".git"))
    @test is_temp

    # Emulate a registry downloaded from a package server.
    downstream_git = gitcmd(registry_dir, TEST_GITCONFIG)
    tree_hash = readchomp(`$(downstream_git) rev-parse HEAD:`)
    rm(joinpath(registry_dir, ".git"), recursive = true)
    write(joinpath(registry_dir, ".tree_info.toml"),
          "git-tree-sha1 = \"$(tree_hash)\"")
    reg_path, is_temp = check_git_registry(registry_dir, TEST_GITCONFIG)
    @test isdir(joinpath(reg_path, ".git"))
    @test is_temp

    # Emulate a registry downloaded from a package server without unpacking.
    tar = read(`$(upstream_git) -c core.autocrlf=false archive $(tree_hash)`)
    gzip = transcode(GzipCompressor, tar)
    write(joinpath(testdir, "TestRegistry.tar.gz"), gzip)
    registry_toml = joinpath(testdir, "TestRegistry.toml")
    write(registry_toml,
          """git-tree-sha1 = "$(tree_hash)"
             uuid = "ed6ca2f6-392d-11ea-3224-d3daf7fee369"
             path = "TestRegistry.tar.gz"
          """)
    if VERSION < v"1.7-"
        @test_throws ErrorException check_git_registry(registry_toml,
                                                       TEST_GITCONFIG)
    else
        reg_path, is_temp = check_git_registry(registry_toml, TEST_GITCONFIG)
        @test isdir(joinpath(reg_path, ".git"))
        @test is_temp
    end
end


text:
struct SimConstructor{T}
    # String or parsed Json
    source :: T
    options :: NamedTuple
end

SimConstructor(source) = SimConstructor(source, NamedTuple())

function addoptions(s::SimConstructor, d)
    options = parse_options(s.options)
    @show options
    d["mapper_options"] = options
    return d
end

function parse(s::SimConstructor{T}) where T
    json_dict = Dict{String,Any}(
        "profile" => deepcopy(s.source),
    )
    return addoptions(s, json_dict)
end

function parse(s::SimConstructor{String})
    profile = open(s.source) do f
        JSON.parse(f)
    end
    json_dict = Dict{String,Any}(
        "profile" => profile,
    )
    return addoptions(s, json_dict)
end

function build_map(sc::SimConstructor)
    json_dict = parse(sc)
    toplevel = build_architecture(sc, json_dict)
    taskgraph = build_taskgraph(sc, json_dict)

    # Get options, do some things
    options = json_dict["mapper_options"]
    rule = getrule(options)

    map = Map(rule, toplevel, taskgraph)
    map.options = options

    return map
end

# Assume a callback is given.
build_architecture(s::SimConstructor, json_dict) = json_dict["mapper_options"][:architecture]()

function build_taskgraph(sc::SimConstructor, json_dict::AbstractDict)
    taskgraph = Taskgraph()
    options = json_dict["mapper_options"]

    parse_input!(sc, taskgraph, json_dict)

    # Run some unpacking transforms on the taskgraph.
    ops = (transform_task_types, compute_edge_metadata, apply_link_weights)
    for op in ops
        taskgraph = op(taskgraph, options)
    end
    return taskgraph
end

function parse_input!(::SimConstructor, taskgraph, json_dict::AbstractDict)
    # Iterate through each node in the dictionary. The keys are the names of
    # the nodes.
    #
    # A second pass will unpack edges.
    for (name, measurements) in json_dict["profile"]
        # Make a metadata dictionary that mimics the structure of what the
        # Project_Manager generates - allowing us to use many of the same passes.
        typestring = titlecase(measurements["Get_Type_String()"], wordsep = isequal('_'))

        metadata = Dict{String,Any}(
            "measurements_dict" => measurements,
            "type" => typestring,
        )
        node = TaskgraphNode(name, metadata = metadata)
        add_node(taskgraph, node)
    end

    # Edges - first pass grabs the explicit edges defined in the profile.
    # Assume all are circuit links.
    #
    # Second pass will unpack any attached memories.
    for node in getnodes(taskgraph)
        # Skip memory nodes, we'll assign those later.
        node.metadata["type"] == "Memory" && continue

        measurements = node.metadata["measurements_dict"]
        # Iterate over output buffers. Must be careful with "nulls"
        for link in get(measurements, "output_buffers", ())
            link == nothing && continue

            reader_core = link["reader_core"]

            # Make a metadata.
            metadata = Dict{String,Any}(
                "pm_class" => "Circuit_Link",
                "measurements_dict" => link,
                "route_link" => true,
                "preserve_dest" => false,
            )
            new_edge = TaskgraphEdge(
                node.name,
                reader_core,
                metadata = metadata
            )

            add_edge(taskgraph, new_edge)
        end
    end

    for node in getnodes(taskgraph)
        attached_memory = get(
            node.metadata["measurements_dict"],
            "attached_memory",
            nothing
        )
        attached_memory == nothing && continue

        # Create an edge to the memory.

        # Proc -> Memory
        metadata = Dict{String,Any}(
            "pm_class" => "Memory_Request_Link",
            "measurements_dict" => Dict{String,Any}(),
            "route_link" => true,
            "preserve_dest" => false,
        )
        new_edge = TaskgraphEdge(node.name, attached_memory, metadata = metadata)
        add_edge(taskgraph, new_edge)

        # Memory -> Proc
        metadata = Dict{String,Any}(
            "pm_class" => "Memory_Response_Link",
            "measurements_dict" => Dict{String,Any}(),
            "route_link" => true,
            "preserve_dest" => false,
        )
        new_edge = TaskgraphEdge(attached_memory, node.name, metadata = metadata)
        add_edge(taskgraph, new_edge)
    end

    return taskgraph
end


text:
using PkgBenchmark
using AMDB

import AMDB: tryparse_base10_tosep, tryparse_base10_fixed

@benchgroup "parsing primitives" begin
    str = b"1980;"
    @bench "parse tosep" tryparse_base10_tosep($str, UInt8(';'), 1)
    @bench "parse fixed" tryparse_base10_fixed($str, 1, 4)
end


text:
# addprocs(6)

using Caesar, RoME
using TransformUtils
using Rotations, CoordinateTransformations
using PyLCM
using PyCall

using KernelDensityEstimate

using JLD, HDF5

@show lcmtpath = joinpath(dirname(@__FILE__),"lcmtypes")

# println("Compile LCM types")
# run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/rome_pose_node_t.lcm`)
# run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/rome_pose_pose_nh_t.lcm`)

println("Adding lcmtypes dir to Python path: $(lcmtpath)")
unshift!(PyVector(pyimport("sys")["path"]),lcmtpath)

@pyimport rome

function whatsthere(channel, msgdata)
  odomsg = rome.pose_pose_nh_t[:decode](msgdata)
  @show odomsg[:utime]
  nothing
end

# MSGDATA = Dict{Symbol, Vector{}}()
# MSGDATA[:pose_node_t] = []
# lastpose = SE3(0)

function lcmposenodes!(vc, slaml::SLAMWrapper, msgdata)
  # push!(MSGDATA[:pose_node_t], msgdata)
  posenode = rome.pose_node_t[:decode](msgdata)
  # @show posenode[:id]
  tA = [posenode[:mean][1];posenode[:mean][2];posenode[:mean][3]]
  EA = Euler(posenode[:mean][6], posenode[:mean][5], posenode[:mean][4])
  q = TransformUtils.convert(Quaternion, EA)
  newpose = SE3(tA, q)
  am = Translation(tA...) ∘ LinearMap(Quat(q.s,q.v...))
  # updaterealtime!(vc,:body, am)
  # visualizerealtime(vc)

  # Dx = slaml.lastpose\newpose # get delta transform between poses
  # copy!(slaml.lastpose.t, newpose.t)
  # copy!(slaml.lastpose.R.R, newpose.R.R)

  @show "here"
  if posenode[:id] == 0
    # first pose prior
    dcovar = [posenode[:covar][1];posenode[:covar][2];posenode[:covar][3];posenode[:covar][4];posenode[:covar][5];posenode[:covar][6]]
    # TODO -- this wont do, must add to existing slaml.fg
    slaml.fg = identitypose6fg(fg=slaml.fg, initpose=newpose, initCov=diagm(dcovar))

    @show length(slaml.fg.g.vertices)
    @show ls(slaml.fg)
    @show lbls, = ls(slaml.fg)
    slaml.usrid2lbl[posenode[:id]] = lbls[1]
    slaml.lbl2usrid[lbls[1]] = posenode[:id]
  end

  rovt( vc, am )
  # change to using pose prior and interpose translation only
  # println("Added node $(v.label)")
  nothing
end

function lcmodomsg!(vc, slaml::SLAMWrapper, msgdata)
  odomsg = rome.pose_pose_nh_t[:decode](msgdata)
  # @show odomsg[:mean_dim]
  # @show odomsg[:node_1_id], odomsg[:node_2_id]
  @show odomsg[:confidence]
  tA = [odomsg[:mean][1];odomsg[:mean][2];odomsg[:mean][3]]
  EA = Euler(odomsg[:mean][6], odomsg[:mean][5], odomsg[:mean][4])
  Dx = SE3(tA, EA)
  dcovar = [odomsg[:covar][1];odomsg[:covar][2];odomsg[:covar][3];odomsg[:covar][4];odomsg[:covar][5];odomsg[:covar][6]]

  if abs(odomsg[:confidence] - 1.0) < 1e-10
    # tA = [odomsg[:mean][1];odomsg[:mean][2];odomsg[:mean][3]]
    # EA = Euler(odomsg[:mean][6], odomsg[:mean][5], odomsg[:mean][4])
    # Dx = SE3(tA, EA)
    # dcovar = [odomsg[:covar][1];odomsg[:covar][2];odomsg[:covar][3];odomsg[:covar][4];odomsg[:covar][5];odomsg[:covar][6]]
    v,f = addOdoFG!(slaml, Pose3Pose3( MvNormal([tA; EA], diagm(dcovar)^2) ) , saveusrid=odomsg[:node_2_id] )
    @show ls(slaml.fg)
  else
    # loop closure case
    println("LOOP CLOSURE")
    @show odomsg[:node_1_id], odomsg[:node_2_id]
    odoc3 = Pose3Pose3NH( MvNormal(veeEuler(Dx), diagm(dcovar)), [0.5;0.5]) # define 50/50% hypothesis
    addFactor!(slam.fg,[Symbol("x$(odomsg[:node_1_id]+1)");Symbol("x$(odomsg[:node_2_id]+1)")],odoc3)
  end
  visualizeallposes!(vc, slaml.fg)

  nothing
end

function runlistener(
      fl::Vector{Bool},
      slam::SLAMWrapper,
      vc;
      until::Symbol=:x999 ) #slam::SLAMWrapper)
  #
  gg = (channel, msg_data) -> lcmposenodes!(vc, slam, msg_data)

  odohdl = (channel, msg_data) -> lcmodomsg!(vc, slam, msg_data)

  print("Preparing LCM...")
  lc = LCM()
  subscribe(lc, "ROME_NODES", gg)
  subscribe(lc, "ROME_FACTORS", odohdl)
  while fl[1]
    handle(lc)
    if slam.lastposesym == until
      break
    end
  end
  println("done")
  nothing
end

function setupSLAMinDB(cloudGraph, addrdict)
  SLAMWrapper(initfg(sessionname=addrdict["session"], cloudgraph=cloudGraph), nothing, 0)
end

function lcmsendpose(lc, vert::Graphs.ExVertex, usrid::Int)
  msg = rome.pose_node_t()
  # p = getBelief(vert)
  p = getKDE(vert)
  mv = getKDEfit(p)
  # pme = getKDEMean(vert)
  msg[:utime] = 0
  msg[:id] = usrid
  msg[:mean_dim] = 6
  msg[:mean] = mv.μ
  msg[:covar_dim] = 6
  msg[:covar] = diag(mv.Σ)

  publish(lc, "MMISAM_POSE_EST_MEAN", msg)
  nothing
end
function lcmsendallposes(lc, slaml::SLAMWrapper)
  for lbluid in slaml.lbl2usrid
    @show lbluid
    lcmsendpose(lc, getVert(slaml.fg, lbluid[1]), lbluid[2])
  end
  nothing
end

vc = startdefaultvisualization()
sleep(3.0)
rovt = loadmodel(:rov)
rovt(vc)

Nparticles = 100
include(joinpath(dirname(@__FILE__),"..","database","blandauthremote.jl"))
addrdict["session"] = "SESSHAUV"
cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict)
# cloudGraph, addrdict = standardcloudgraphsetup()

slam = setupSLAMinDB(cloudGraph, addrdict)

flags = Bool[1]
flags[1] = true

println("LCM listener running...")
runlistener(flags, slam, vc, until=:x5)

# @show 1
#
#
#
# println("AD-HOC DEBUGGING")
#
# lcmposenodes!(vc, slam, MSGDATA[:pose_node_t][1])
# splice!(MSGDATA[:pose_node_t], 2)
#
# @show -1
#
# #
# # # run this to stop the lcm listener loop
# # flags[1] = false
# #
# #
# solve(slam.fg) #, densitymeshes=[:x1;:x33;:x60])
# #
# visualizeallposes!(vc, slam.fg, drawtype=:fit)
#
#
#
#
# lc = LCM()
# lcmsendallposes(lc, slam)
#
#
#
# visualizeDensityMesh!(vc, slam.fg, :x1,meshid=2)

# savefg = deepcopy(slam.fg)
# savefg.registeredModuleFunctions = nothing
# @save joinpath(dirname(@__FILE__),"hauvnh67.jld") savefg

#
# using Gadfly
#
#
# draw(PDF("/home/dehann/Desktop/test.pdf",30cm,20cm),
#         plotKDE( marginal(p60,collect(4:6)),
# 				dimLbls=["x";"y";"z";"phi";"the";"psi"]) )

#