text:
# Read saved R datasets in the RDA2 or RDX2 format
# Files written in this format often have the extension .rda or .RData

##############################################################################
##
## A table of R's SEXPREC tags.  No longer used but good to have around.
##
##############################################################################

const SXPtab = @compat Dict(      # Defined in Rinternals.h
    0x00=>"NULL",
    0x01=>"Symbol",
    0x02=>"Pairlist",
    0x03=>"Closure",
    0x04=>"Environment",
    0x05=>"Promise",
    0x06=>"Lang",
    0x07=>"Special",
    0x08=>"Builtin",
    0x09=>"Char",          # "scalar" string type (internal only)
    0x0a=>"Logical",       # almost BitArray but allows NA's
    0x0d=>"Integer",       # Array{Int32, 1}
    0x0e=>"Real",          # Array{Float64, 1}
    0x0f=>"Complex",       # Array{Complex128, 1}
    0x10=>"String",        # Array{ASCIIString, 1}
    0x11=>"Dot",           # dot-dot-dot object
    0x12=>"Any",           # make "any" args work
    0x13=>"List",          # generic vector, {} but with names
    0x14=>"Expr",          # expressions vectors
    0x15=>"ByteCode",
    0x16=>"XPtr",
    0x17=>"WeakRef",
    0x18=>"Raw",
    0x19=>"S4",
    0x1e=>"New",           # fresh node created in new page
    0x1f=>"Free",          # node released by GC
    0x63=>"Function",      # closure or builtin
                           # Defined in serialize.c
    0xf1=>"BaseEnv",
    0xf2=>"EmptyEnv",
    0xf3=>"BCREPREF",
    0xf4=>"BCREPDEF",
    0xf5=>"GenericRef",
    0xf6=>"ClassRef",
    0xf7=>"Persist",
    0xf8=>"Package",
    0xf9=>"Namespace",
    0xfa=>"BaseNamespace",
    0xfb=>"MissingArg",
    0xfc=>"UnboundValue",
    0xfd=>"GlobalEnv",
    0xfe=>"NilValue",      # terminates a pairs list?
    0xff=>"Ref"
)

##############################################################################
##
## Constants used as NA patterns in R.
## (I assume 1954 is the year of Ross's birth or something like that.)
##
##############################################################################

if ENDIAN_BOM == 0x01020304
    const R_NA_FLOAT64 = reinterpret(Float64, [0x7ff00000, @compat(UInt32(1954))])[1]
else
    const R_NA_FLOAT64 = reinterpret(Float64, [@compat(UInt32(1954)), 0x7ff00000])[1]
end
const R_NA_INT32 = typemin(Int32)
const R_NA_STRING = "NA"

##############################################################################
##
## Julia types encapsulating various types of R SEXPRECs
##
##############################################################################

typealias Hash Dict{ASCIIString, Any}
const nullhash = Hash()

abstract RSEXPREC{S}             # Basic R object - symbolic expression

type RSymbol <: RSEXPREC{0x01}   # Not quite the same as a Julia symbol
    displayname::ASCIIString
end

abstract ROBJ{S} <: RSEXPREC{S}  # R object that can have attributes

abstract RVEC{T, S} <: ROBJ{S}   # abstract R vector (actual storage implementation may differ)

type RVector{T, S} <: RVEC{T, S} # R vector object
    data::Array{T, 1}
    attr::Hash                   # collection of R object attributes
end

typealias RLogical RVector{Int32, 0x0a}
typealias RInteger RVector{Int32, 0x0d}
typealias RNumeric RVector{Float64, 0x0e}
typealias RComplex RVector{Complex128, 0x0f}

type RNullableVector{T, S} <: RVEC{T, S} # R vector object with explicit NA values
    data::Array{T, 1}
    na::BitArray{1}             # mask of NA elements
    attr::Hash                  # collection of R object attributes
end

typealias RString RNullableVector{ASCIIString,0x10}
typealias RList RVector{Any, 0x13}  # "list" in R == Julia cell array

##############################################################################
##
## R objects in the file are preceded by a UInt32 giving the type and
## some flags.  These functiona unpack bits in the flags.  The isobj
## bit might be useful for distinguishing an RInteger from a factor or
## an RList from a data.frame.
##
##############################################################################

typealias RDATag UInt32

isobj(fl::RDATag) = (fl & 0x00000100) != 0
hasattr(fl::RDATag) = (fl & 0x00000200) != 0
hastag(fl::RDATag) = (fl & 0x00000400) != 0

if VERSION < v"0.4-"
    sxtype = uint8
else
    sxtype(fl::RDATag) = fl % UInt8
end

##############################################################################
##
## Utilities for reading a single data element.
## The read<type>orNA functions are needed because the ASCII format
## stores the NA as the string 'NA'.  Perhaps it would be easier to
## wrap the conversion in a try/catch block.
##
##############################################################################

LONG_VECTOR_SUPPORT = (WORD_SIZE > 32) # disable long vectors support on 32-bit machines

if LONG_VECTOR_SUPPORT
    typealias RVecLength Int64
else
    typealias RVecLength Int
end

# abstract RDA format IO stream wrapper
abstract RDAIO

type RDAXDRIO{T<:IO} <: RDAIO # RDA XDR(binary) format IO stream wrapper
    sub::T             # underlying IO stream
    buf::Vector{UInt8} # buffer for strings

    RDAXDRIO( io::T ) = new( io, Array(UInt8, 1024) )
end
RDAXDRIO{T <: IO}(io::T) = RDAXDRIO{T}(io)

readint32(io::RDAXDRIO) = ntoh(read(io.sub, Int32))
readuint32(io::RDAXDRIO) = ntoh(read(io.sub, UInt32))
readfloat64(io::RDAXDRIO) = ntoh(read(io.sub, Float64))

readintorNA(io::RDAXDRIO) = readint32(io)
readintorNA(io::RDAXDRIO, n::RVecLength) = map!(ntoh, read(io.sub, Int32, n))

readfloatorNA(io::RDAXDRIO) = readfloat64(io)
readfloatorNA(io::RDAXDRIO, n::RVecLength) = map!(ntoh, read(io.sub, Float64, n))

function readnchars(io::RDAXDRIO, n::Int32)  # a single character string
    readbytes!(io.sub, io.buf, n)
    bytestring(pointer(io.buf), n)::ASCIIString
end

type RDAASCIIIO{T<:IO} <: RDAIO # RDA ASCII format IO stream wrapper
    sub::T              # underlying IO stream

    RDAASCIIIO( io::T ) = new( io )
end
RDAASCIIIO{T <: IO}(io::T) = RDAASCIIIO{T}(io)

readint32(io::RDAASCIIIO) = parse(Int32, readline(io.sub))
readuint32(io::RDAASCIIIO) = parse(UInt32, readline(io.sub))
readfloat64(io::RDAASCIIIO) = parse(Float64, readline(io.sub))

function readintorNA(io::RDAASCIIIO)
    str = chomp(readline(io.sub));
    str == R_NA_STRING ? R_NA_INT32 : parse(Int32, str)
end
readintorNA(io::RDAASCIIIO, n::RVecLength) = Int32[readintorNA(io) for i in 1:n]

function readfloatorNA(io::RDAASCIIIO)
    str = chomp(readline(io.sub));
    str == R_NA_STRING ? R_NA_FLOAT64 : parse(Float64, str)
end
readfloatorNA(io::RDAASCIIIO, n::RVecLength) = Float64[readfloatorNA(io) for i in 1:n]

function readnchars(io::RDAASCIIIO, n::Int32)  # reads N bytes-sized string
    if (n==-1) return "" end
    str = unescape_string(chomp(readline(io.sub)))
    length(str) == n || error("Character string length mismatch")
    str
end

type RDANativeIO{T<:IO} <: RDAIO # RDA native binary format IO stream wrapper (TODO)
    sub::T               # underlying IO stream

    RDANativeIO( io::T ) = new( io )
end
RDANativeIO{T <: IO}(io::T) = RDANativeIO{T}(io)

function rdaio(io::IO, formatcode::AbstractString)
    if formatcode == "X" RDAXDRIO(io)
    elseif formatcode == "A" RDAASCIIIO(io)
    elseif formatcode == "B" RDANativeIO(io)
    else error( "Unrecognized RDA format \"$formatcode\"" )
    end
end

if LONG_VECTOR_SUPPORT
    # reads the length of any data vector from a stream
    # from R's serialize.c
    function readlength(io::RDAIO)
        len = convert(RVecLength, readint32(io))
        if (len < -1) error("negative serialized length for vector")
        elseif (len >= 0)
            return len
        else # big vectors, the next 2 ints encode the length
            len1, len2 = convert(RVecLength, readint32(io)), convert(RVecLength, readint32(io))
            # sanity check for now
            if (len1 > 65536) error("invalid upper part of serialized vector length") end
            return (len1 << 32) + len2
        end
    end
else
    # reads the length of any data vector from a stream
    # fails when long (> 2^31-1) vector encountered
    # from R's serialize.c
    function readlength(io::RDAIO)
        len = convert(RVecLength, readint32(io))
        if (len >= 0)
            return len
        elseif (len < -1)
            error("negative serialized length for vector")
        else
            error("negative serialized vector length:\nperhaps long vector from 64-bit version of R?")
        end
    end
end

immutable CHARSXProps # RDA CHARSXP properties
  levs::UInt32       # level flags (encoding etc) TODO process
  nchar::Int32       # string length, -1 for NA strings
end

function readcharsxprops(io::RDAIO) # read character string encoding and length
    fl = readuint32(io)
    @assert sxtype(fl) == 0x09
### watch out for levs in here.  Generally it has the value 0x40 so that fl = 0x00040009 (262153)
### if levs == 0x00 then the next line should be -1 to indicate the NA_STRING
    CHARSXProps(fl >> 12, readint32(io))
end

function readcharacter(io::RDAIO)  # a single character string
    props = readcharsxprops(io)
    props.nchar==-1 ? "" : readnchars(io, props.nchar)
end

function readcharacter(io::RDAIO, n::RVecLength)  # a single character string
    res = fill("", n)
    na = falses(n)
    for i in 1:n
        props = readcharsxprops(io)
        if (props.nchar==-1) na[i] = true
        else res[i] = readnchars(io, props.nchar)
        end
    end
    return res, na
end

##############################################################################
##
## Utilities for reading compound RDA items: lists, arrays etc
##
##############################################################################

type RDAContext{T <: RDAIO}    # RDA reading context
    io::T                      # R input stream

    # RDA properties
    fmtver::UInt32             # RDA format version
    Rver::VersionNumber        # R version that has written RDA
    Rmin::VersionNumber        # R minimal version to read RDA

    # behaviour
    convertdataframes::Bool    # if R dataframe objects should be automatically converted into DataFrames

    # intermediate data
    symtab::Array{RSymbol,1}   # symbols array

    function RDAContext(io::T, kwoptions::Array{Any})
        fmtver = readint32(io)
        rver = readint32(io)
        rminver = readint32(io)
        kwdict = Dict{Symbol,Any}(kwoptions)
        new(io,
            fmtver,
            VersionNumber( div(rver,65536), div(rver%65536, 256), rver%256 ),
            VersionNumber( div(rminver,65536), div(rminver%65536, 256), rminver%256 ),
            get(kwdict,:convertdataframes,false),
            Array(RSymbol,0))
    end
end

RDAContext{T <: RDAIO}(io::T, kwoptions::Array{Any}) = RDAContext{T}(io, kwoptions)

function readnamedobjects(ctx::RDAContext, fl::RDATag)
    if !hasattr(fl) return nullhash end
    res = Hash()
    fl = readuint32(ctx.io)
    while sxtype(fl) != 0xfe
        ## need to call RSymbol here b/c of symbol reference table
        nm = readsymbol(ctx, readuint32(ctx.io)).displayname
        setindex!(res, readitem(ctx), nm)
        fl = readuint32(ctx.io)
    end
    res
end

readattrs(ctx::RDAContext, fl::RDATag) = readnamedobjects(ctx, fl)

function readnumeric(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == 0x0e
    RNumeric(readfloatorNA(ctx.io, readlength(ctx.io)),
             readattrs(ctx, fl))
end

function readinteger(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == 0x0d
    RInteger(readintorNA(ctx.io, readlength(ctx.io)),
             readattrs(ctx, fl))
end

function readlogical(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == 0x0a # excluding this check, the method is the same as readinteger()
    RLogical(readintorNA(ctx.io, readlength(ctx.io)),
             readattrs(ctx, fl))
end

function readcomplex(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == 0x0f
    n = readlength(ctx.io)
    data = readfloatorNA(ctx.io, 2n)
    RComplex(Complex128[@compat(Complex128(data[i],data[i+1])) for i in 2(1:n)-1],
             readattrs(ctx, fl))
end

function readstring(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == 0x10
    RString(readcharacter(ctx.io, readlength(ctx.io))...,
            readattrs(ctx, fl))
end

function readlist(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == 0x13
    n = readlength(ctx.io)
    res = RList([readitem(ctx) for i in 1:n],
                readattrs(ctx, fl))
    if ctx.convertdataframes && isdataframe(res)
        DataFrame(res)
    else
        res
    end
end

function readsymbol(ctx::RDAContext, fl::RDATag)
    # check if its a reference to an already defined string
    if sxtype(fl) == 0xff return ctx.symtab[fl >> 8] end
    # read the new strings and put it into symbols table
    @assert sxtype(fl) == 0x01
    res = RSymbol(readcharacter(ctx.io))
    push!(ctx.symtab, res)
    res
end

const SXReaders = @compat Dict( # maps RDA type to function that read it
    0x01 => readsymbol,
    0x0a => readlogical,
    0x0d => readinteger,
    0x0e => readnumeric,
    0x0f => readcomplex,
    0x10 => readstring,
    0x13 => readlist
);

function readitem(ctx::RDAContext)
    ff = readuint32(ctx.io)
    fl = sxtype(ff)
    if haskey(SXReaders, fl) return SXReaders[fl](ctx, ff) end
### Should not occur at the top level
###    if fl == 0xf3 return nothing end      # terminates dotted pair lists
###    if fl == 0x09 return readcharacter(ctx.io, ff) end
    error("Encountered flag $ff corresponding to type $(SXPtab[fl])")
end

function read_rda(io::IO, kwoptions::Array{Any})
    header = chomp(readline(io))
    @assert header[1] == 'R' # readable header (or RDX2)
    @assert header[2] == 'D'
    @assert header[4] == '2'
    ctx = RDAContext(rdaio(io, chomp(readline(io))), kwoptions)
    @assert ctx.fmtver == 2    # format version
#    println("Written by R version $(ctx.Rver)")
#    println("Minimal R version: $(ctx.Rmin)")
    return readnamedobjects(ctx, 0x00000200)
end

read_rda(io::IO; kwoptions...) = read_rda(io, kwoptions)

read_rda(fnm::AbstractString; kwoptions...) = gzopen(fnm) do io read_rda(io, kwoptions) end

##############################################################################
##
## Utilities for working with basic properties of R objects:
##    attributes, class inheritance, etc
##
##############################################################################

const emptystrvec = Array(ASCIIString,0)

getattr{T}(ro::ROBJ, attrnm::ASCIIString, default::T) = haskey(ro.attr, attrnm) ? ro.attr[attrnm].data : default;

Base.names(ro::ROBJ) = getattr(ro, "names", emptystrvec)

class(ro::ROBJ) = getattr(ro, "class", emptystrvec)
class(x) = emptystrvec
inherits(x, clnm::ASCIIString) = any(class(x) .== clnm)

isdataframe(rl::RList) = inherits(rl, "data.frame")
isfactor(ri::RInteger) = inherits(ri, "factor")

Base.length(rl::RVEC) = length(rl.data)
Base.size(rv::RVEC) = length(rv.data)
Base.size(rl::RList) = isdataframe(rl) ? (length(rl.data[1]), length(rl.data)) : length(rl.data)

row_names(ro::ROBJ) = getattr(ro, "row.names", emptystrvec)

##############################################################################
##
## Conversion of intermediate R objects into DataArray and DataFrame objects
##
##############################################################################

namask(rl::RLogical) = bitpack(rl.data .== R_NA_INT32)
namask(ri::RInteger) = bitpack(ri.data .== R_NA_INT32)
namask(rn::RNumeric) = bitpack([rn.data[i] === R_NA_FLOAT64 for i in 1:length(rn.data)])
namask(rc::RComplex) = bitpack([rc.data[i].re === R_NA_FLOAT64 ||
                                rc.data[i].im === R_NA_FLOAT64 for i in 1:length(rc.data)])
namask(rv::RNullableVector) = rv.na

DataArrays.data(rv::RVEC) = DataArray(rv.data, namask(rv))

function DataArrays.data(ri::RInteger)
    if !isfactor(ri) return DataArray(ri.data, namask(ri)) end
    # convert factor into PooledDataArray
    pool = getattr(ri, "levels", emptystrvec)
    sz = length(pool)
    REFTYPE = sz <= typemax(UInt8)  ? UInt8 :
              sz <= typemax(UInt16) ? UInt16 :
              sz <= typemax(UInt32) ? UInt32 :
                                      UInt64
    dd = ri.data
    dd[namask(ri)] = 0
    refs = convert(Vector{REFTYPE}, dd)
    return PooledDataArray(DataArrays.RefArray(refs), pool)
end

function DataFrame(rl::RList)
    DataFrame(map(data, rl.data),
              Symbol[identifier(x) for x in names(rl)])
end


text:
push!(LOAD_PATH,"../src/")

using Documenter, SymEngine

makedocs(sitename="Symengine Julia API Docs")

ext: jl


text:
using Documenter

makedocs(
    sitename = "JuliaAtoms",
    pages = [
        "Home" => "index.md",
    ],
    format = Documenter.HTML(
        edit_link = "source",
    ),
)

deploydocs(
    repo = "github.com/JuliaAtoms/juliaatoms.github.io.git",
    branch = "master",
    devbranch = "source",
)


text:
module TestOptimization
# testing helpers for optimization of CovarianceFunctions hyper-parameters
using Test
using CovarianceFunctions
using CovarianceFunctions: EQ, RQ, Lengthscale
using LinearAlgebra
using Flux

# testing taking gradients of kernels and incrementing parameters with them
@testset "optimization" begin
    l = [1/2]
    k = Lengthscale(EQ(), l)
    p = Flux.params(k)
    x, y = randn(), randn()
    g = gradient(()->k(x, y), p) # should we use @forward?
    @test g[l][1] isa Real
end

end # TestOptimization


text:
TrueFX = Val{:TrueFX}

TrueFX_username = "TrueFX_username"
TrueFX_password = "TrueFX_password"
TrueFX_qualifier = "eurates"
TrueFX_format = "csv" # or "html"
TrueFX_currencypairs = "EUR/USD,GBP/USD,CAN/USD,AUD/USD"
ENV["TrueFX_sessionid"] = ""

try
   TrueFX_username = ENV[TrueFX_username]
   TrueFX_password = ENV[TrueFX_password]
catch  KeyError
    throw(ErrorException(
        """You must register a username and password with TrueFX and put them in the environment.\n\t
           Use your .juliarc.jl file:\n\t\t
           ENV["TrueFX_username"] = "<your username>"\n\t\t
           ENV["TrueFX_password"] = "<your password>"\n
        """))
end

const truefx_connect_url = "http://webrates.truefx.com/rates/connect.html"
const truefx_authentication_url = string(truefx_connect_url,"?u=",TrueFX_username,"&p=",TrueFX_password)
truefx_authorization_url = string(truefx_authentication_url,"&q=",TrueFX_qualifier,"&f=",TrueFX_format)

# e.g. "EUR/USD" or "EUR/USD,EUR/GBP"
function truefx_currency_pairs(currencypairs::V) where V<:AbstractVector{S} where S<:String
   return join(V, ",")
end
truefx_currency_pairs(currencypairs::String) = currencypairs

function truefx_currency_pairs!(currencypairs::String)
    global truefx_authorization_url
    truefx_authorization_url = split(truefx_authorization_url,"&c=")[1]
    truefx_authorization_url = string(truefx_authorization_url,"&c=",currencypairs)
    return nothing
end

truefx_authorization_url = truefx_currency_pairs!(truefx_currency_pairs)

function sessionid(::Type{TrueFX}; force::Bool=false)
    global truefx_authorization_url
    if force || length(ENV["TrueFX_sessionid"]) == 0
        response = Requests.get(truefx_authorization_url)
        newsessionid = chomp(readstring(response))
        verify_authorization(TrueFX, newsessionid, response.status)
        ENV["TrueFX_sessionid"] = newsessionid
    end
    return ENV["TrueFX_sessionid"]
end

function verify_authorization(::Type{TrueFX}, obtained::String, status::Int)
    if status != 200
        throw(ErrorException("TrueFX did not process the authorization: status error."))
    end
    if obtained === "not authorized"
        throw(ErrorException("TrueFX did not authorize the username and password given"))
    end
end


text:
#
# Copyright (c) 2021 Frederic Bruder, Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

function fmiLoad(name::AbstractString, tool::AbstractString, version::AbstractString, fmiversion::AbstractString="2.0"; kwargs...)
    FMI.fmiLoad(get_model_filename(name, tool, version, fmiversion); kwargs...)
end

function fmi2Load(name::AbstractString, tool::AbstractString, version::AbstractString; kwargs...)
    FMIImport.fmi2Load(get_model_filename(name, tool, version, "2.0"); kwargs...)
end

function fmi3Load(name::AbstractString, tool::AbstractString, version::AbstractString; kwargs...)
    FMIImport.fmi3Load(get_model_filename(name, tool, version, "3.0"); kwargs...)
end


text:
using Test
using LightGraphs
using DataStructures

n = 100
threshold = 10
g = erdos_renyi(n, 4/n)

@testset "HDA" begin
    attack_nodes = HDA(g, threshold)
    present = fill(true, n)
    present[attack_nodes] .= false
    @test giant_component_size(g, present) <= threshold
end

@testset "CoreHD" begin
    attack_nodes = CoreHD(g, threshold)
    present = fill(true, n)
    present[attack_nodes] .= false
    @test giant_component_size(g, present) <= threshold
end

@testset "WeakNei" begin
    attack_nodes = WeakNei(g, threshold)
    present = fill(true, n)
    present[attack_nodes] .= false
    @test giant_component_size(g, present) <= threshold
end

@testset "tree_break" begin
    decycling_nodes = decore(g)
    present = fill(true, n)
    present[decycling_nodes] .= false
    tree_break!(present, g, threshold)
    @test giant_component_size(g, present) <= threshold
end


text:
"""
    greenseltype(::Type{DQMC}, m::Model)

Returns the type of the elements of the Green's function matrix. Defaults to `ComplexF64`.
"""
greenseltype(::Type{DQMC}, m::Model) = ComplexF64

"""
    energy(mc::DQMC, m::Model, conf)

Calculate bosonic part (non-Green's function determinant part) of energy for configuration `conf` for Model `m`.
"""
energy_boson(mc::DQMC, m::Model, conf) = 0.

"""
    prepare_observables(mc::DQMC, m::Model) -> Dict{String, Observable}

Initializes observables and returns a `Dict{String, Observable}`. In the latter,
keys are abbreviations for the observables names and values are the observables themselves.

See also [`measure_observables!`](@ref) and [`finish_observables!`](@ref).
"""
prepare_observables(mc::DQMC, m::Model) = Dict{String, Observable}()

"""
    measure_observables!(mc::DQMC, m::Model, obs::Dict{String,Observable}, conf)

Measures observables and updates corresponding `MonteCarloObservable.Observable` objects in `obs`.

See also [`prepare_observables`](@ref) and [`finish_observables!`](@ref).
"""
measure_observables!(mc::DQMC, m::Model, obs::Dict{String,Observable}, conf) = nothing

"""
    finish_observables!(mc::DQMC, m::Model, obs::Dict{String,Observable})

Finish measurements of observables.

See also [`prepare_observables`](@ref) and [`measure_observables!`](@ref).
"""
finish_observables!(mc::DQMC, m::Model, obs::Dict{String,Observable}) = nothing

max_stars_repo_name: yosefcap/MC

max_issues_repo_name: yosefcap/MC

max_forks_repo_name: yosefcap/MC


text:
λ(∇ϕxi, θi, Flow::ZigZag) = pos(∇ϕx*θ)

function event(t, x, θ, Z::ZigZag)
    t, copy(x), copy(θ)
end

function pdmp_inner_sub!(Ξ, ∇ϕ!, ∇ϕx, t, x, θ, c, a, b, t′, τref, (acc, num),
     Flow::ZigZag, args...; factor=1.5, adapt=false)
    while true
        if τref < t′
            error("don't go here")
            t, x, θ = move_forward!(τref - t, t, x, θ, Flow)
            #θ = randn!(θ)
            θ = refresh!(θ, Flow)
            τref = t + waiting_time_ref(Flow)
            a, b = ab(x, θ, c, Flow, args...)
            t′ = t + poisson_time(a, b, rand())
            push!(Ξ, event(t, x, θ, Flow))
            return t, x, θ, (acc, num), c, a, b, t′, τref
        else
            τ = t′ - t
            t, x, θ = move_forward!(τ, t, x, θ, Flow)
            (E_tilde, i) = ∇ϕ!(x, args...)
            l, lb = λ(E_tilde, θ[i], Flow), pos(a + b*τ)
            num += 1
            if rand()*lb <= l
                acc += 1
                if l > lb
                    !adapt && error("Tuning parameter `c` too small.")
                    c *= factor
                end
                θ[i] = -θ[i]
                push!(Ξ, event(t, x, θ, Flow))
                a, b = ab(x, θ, c, Flow)
                t′ = t + poisson_time(a, b, rand())
                return t, x, θ, (acc, num), c, a, b, t′, τref
            end
            a, b = ab(x, θ, c, Flow, args...)
            t′ = t + poisson_time(a, b, rand())
        end
    end
end

function pdmp_sub(∇ϕ!, t0, x0, θ0, T, c, Flow::ZigZag, args...; adapt=false, factor=2.0)
    t, x, θ, ∇ϕx = t0, copy(x0), copy(θ0), copy(θ0)
    Ξ = Trace(t0, x0, θ0, Flow)
    τref = waiting_time_ref(Flow)
    num = acc = 0
    a, b = ab(x, θ, c, Flow, args...)
    t′ = t + poisson_time(a, b, rand())
    while t < T
        t, x, θ, (acc, num), c, a, b, t′, τref = pdmp_inner_sub!(Ξ, ∇ϕ!, ∇ϕx, t, x, θ, c, a, b, t′, τref, (acc, num), Flow, args...; factor=factor, adapt=adapt)
    end
    return Ξ, (t, x, θ), (acc, num), c
end


text:
using Mimi

@defcomp ConnectorCompVector begin
    input1 = Parameter(index = [time])
    input2 = Parameter(index = [time])
    output =  Variable(index = [time])

    first = Parameter() # first year to use the shorter data
    last = Parameter()  # last year to use the shorter data

    function run_timestep(p, v, d, t)
        if gettime(t) >= p.first && gettime(t) <= p.last
            input = p.input1
        else
            input = p.input2
        end

        v.output[t] = @allow_missing(input[t])

    end
end

@defcomp ConnectorCompMatrix begin
    ConnectorCompMatrix_Dim2 = Index()

    input1 = Parameter(index = [time, ConnectorCompMatrix_Dim2])
    input2 = Parameter(index = [time, ConnectorCompMatrix_Dim2])
    output =  Variable(index = [time, ConnectorCompMatrix_Dim2])

    first = Parameter() # first year to use the shorter data
    last = Parameter()  # last year to use the shorter data

    function run_timestep(p, v, d, t)

        if gettime(t) >= p.first && gettime(t) <= p.last
            input = p.input1
        else
            input = p.input2
        end

        for r in d.ConnectorCompMatrix_Dim2
            v.output[t, r] = @allow_missing(input[t, r])
        end
    end
end

#
# TBD: define a version with arbitrary dimensions. The problem currently
# is that we have no way to indicate a Parameter with an indeterminate
# dimensions.
#
# @defcomp ConnectorComp begin
#     input1 = Parameter(index = [time, ...])
#     input2 = Parameter(index = [time, ...])
#     output =  Variable(index = [time, ...])

#     # Allow copying of vars/params with arbitrary dimensions
#     function run_timestep(p, v, d, t)
#         colons = repeat([:], inner=ndims(v.output) - 1)
#         if hasvalue(p.input1, t)
#             v.output[t, colons...] = p.input1[t, colons...]
#         elseif hasvalue(p.input2, t)
#             v.output[t, colons...] = p.input2[t, colons...]
#         else
#             error("Neither of the inputs to ConnectorComp have data for the current timestep: $(gettime(t)).")
#         end
#     end
# end


text:
type Linear end
type Circular end
type General end
type Permuted end

function müller_rec(; ϵ::Float64 = 0.0, ϕ::Float64 = 0.0, ΔG::Float64 = 0.0, α::Float64 = 45.0*π/180, ψ::Float64 = 0.0)
    return [1.0         (-2ϵ*sin(ϕ)*sin(2α) + ΔG/2*cos(2α)) 2ϵ*cos(ϕ) (2ϵ*sin(ϕ)*cos(2α) + ΔG/2*sin(2α))
            ΔG/2                     cos(2α)                   0.0                  sin(2α)
            2ϵ*cos(ϕ + ψ)          sin(2α)*sin(ψ)             cos(ψ)            -cos(2α)*sin(ψ)
            2ϵ*sin(ϕ + ψ)         -sin(2α)*cos(ψ)             sin(ψ)             cos(2α)*cos(ψ)]
end
function müller_rec(::Permuted; ϵ::Float64 = 0.0, ϕ::Float64 = 0.0, ΔG::Float64 = 0.0, α::Float64 = 45.0*π/180, ψ::Float64 = 0.0)
    return [1.0         (-2ϵ*sin(ϕ)*sin(2α) + ΔG/2*cos(2α)) 2ϵ*cos(ϕ) (2ϵ*sin(ϕ)*cos(2α) + ΔG/2*sin(2α))
            -2ϵ*sin(ϕ + ψ)         sin(2α)*cos(ψ)             -sin(ψ)             -cos(2α)*cos(ψ)
            2ϵ*cos(ϕ + ψ)          sin(2α)*sin(ψ)             cos(ψ)            -cos(2α)*sin(ψ)
            ΔG/2                     cos(2α)                   0.0                  sin(2α)]

end

function müller_rec(::Linear; ϵ::Float64 = 0.0, ϕ::Float64 = 0.0, ΔG::Float64 = 0.0, δα::Float64 = 0.0, ψ::Float64 = 0.0, s = +)
    pm = s(1.0)
    mp = -pm
    return [1.0           pm*ΔG/2       2ϵ*cos(ϕ) 2ϵ*sin(ϕ)
            ΔG/2          pm            0.0       pm*2δα
            2ϵ*cos(ϕ + ψ) pm*2δα*sin(ψ) cos(ψ)    mp*sin(ψ)
            2ϵ*sin(ϕ + ψ) mp*2δα*cos(ψ) sin(ψ)    pm*cos(ψ)]
end

function müller_rec(::Circular; ϵ::Float64 = 0.0, ϕ::Float64 = 0.0, ΔG::Float64 = 0.0, δα::Float64 = 0.0, ψ::Float64 = 0.0, s = +)
    pm = s(1.0)
    mp = -pm
    return [1.0           mp*2ϵ*sin(ϕ)  2ϵ*cos(ϕ)  pm*ΔG/2
            ΔG/2          mp*2δα          0.0        pm
            2ϵ*cos(ϕ + ψ) pm*sin(ψ)     cos(ψ)   pm*2δα*sin(ψ)
            2ϵ*sin(ϕ + ψ) mp*cos(ψ)     sin(ψ)   mp*2δα*cos(ψ)]
end

function müller_sky(PAaz::Float64)
    cpa = cos(2*PAaz)
    spa = sin(2*PAaz)
    return [1.0   0.0   0.0  0.0
            0.0   cpa   spa  0.0
            0.0  -spa   cpa  0.0
            0.0   0.0   0.0  1.0]
end

function lin_pol_source(;flux::Float64 = 1.0, p::Float64 = 0.1, pa::Float64 = 0.0)
    return Float64[flux
                   p*flux*cos(2*pa)
                   p*flux*sin(2*pa)
                   0.0]
end

max_stars_repo_name: KajWiik/PolarimetryCalibration

max_issues_repo_name: KajWiik/PolarimetryCalibration

max_forks_repo_name: KajWiik/PolarimetryCalibration


text:
# The controller node
# This file is part of the openBuildNet framework developed at EPFL.
# (C) 2016 LA-EPFL
# Author: Truong X. Nghiem (truong.nghiem@gmail.com, xuan.nghiem@epfl.ch)

# MPC node
# See: http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Examples.StandardMPC

using openbuildnet

# Model data
const A = [2. -1.;1. 0.2]
const B = [1. 0.]'
const x0 = [3., 1.];  # Initial state

# MPC data
nx = size(A,2); # Number of states
nu = size(B,2); # Number of inputs

# MPC data
Q = eye(nx);
R = 2.;
N = 7;

# Prepare the optimization model
using GLPKMathProgInterface
using Convex

function create_model()
  local u = Variable(nu, N);
  local x = Variable(nx, N);
  local x0v = Variable(nx);

  local constraints = [-1 <= u, u <= 1, -5 <= x, x <= 5];
  constraints += x[:,1] == A*x0v + B*u[:,1];
  objective = norm(Q*x0v,1) + norm(R*u[:,1],1);
  for k = 2:N
    objective = objective + norm(Q*x[:,k-1],1) + norm(R*u[:,k],1);
    constraints += x[:,k] == A*x[:,k-1] + B*u[:,k];
  end

  controller = minimize(objective, constraints);

  # returns a tuple of the controller, the parameter, and the solution
  (controller, x0v, u)
end

tic()
controller, x0v, u = create_model()
toc()

# Create the node
node = OBNNode("controller", "simplempc")

# Create the ports
input_x = create_input(node, "x", Vector{Float64})  # feedback
output_u = create_output(node, "u", Vector{Float64}) # control values

# Create Callbacks
on_init(node) do
  println("INIT")
end

on_term(node) do
  println("TERMINATE")
end

const MAINBLOCK = 0

US = Vector{Vector{Float64}}()  # To store the control values

on_block_output(node, MAINBLOCK) do
  global controller, x0v, u
  fix!(x0v, get(input_x))
  tic()
  solve!(controller, GLPKSolverLP())
  toc()
  uk = u.value[:,1]
  push!(US, uk)
  set(output_u, uk)
end

println("Ready to start simulation...")
status = runsafe(node, 20)
println("Final status: ", status)

# Plot the results
if status == 2
  using PyPlot
  USmat = hcat(US...)' # convert to matrix
  plot(USmat[:,1])
  show()
end


text:
using ApproxFun, Test, DualNumbers
    import ApproxFun: eigs

@testset "Extras" begin
    @testset "Dual numbers" begin
        @test dual(1.5,1) ∈  Segment(dual(1.0,1),dual(2.0))

        f=Fun(exp,Segment(dual(1.0,1),dual(2.0)),20)
        @test Fun(h->Fun(exp,Segment(1.0+h,2.0)).coefficients[1],0..1)'(0.) ≈ DualNumbers.epsilon(f.coefficients[1])

        ud=let d= Segment(dual(0.0,1.0),1.0)
            B = ldirichlet(d)
            D = Derivative(d)
            a = Fun(exp,d)
            u = [B;D+dual(1.0,4.0)*a] \ [dual(1.0,2.0),0.0]
            u(0.5)
        end

        u0=let d= Segment(0.0,1.0)
            B = ldirichlet(d)
            D = Derivative(d)
            a = Fun(exp,d)
            u = [B;D+a] \ [1.0,0.0]
            u(0.5)
        end
        h=0.00001
        uh=let d= Segment(h,1.0)
            B = ldirichlet(d)
            D = Derivative(d)
            a = Fun(exp,d)
            u = [B;D+(1+4h)*a] \ [1.0+2h,0.0]
            u(0.5)
        end

        @test absdual(ud - dual(u0,(uh-u0)/h)) ≤ h

        let d=Segment(0.0,1.0)
            B = ldirichlet(d)
            D = Derivative(d)
            a = Fun(exp,d)
            u = [B;D+a] \ [dual(1.0,2.0),0.0]
            ur = [B;D+a] \ [1.0,0.0]
            ud = [B;D+a] \ [2.0,0.0]
            @test absdual(u(0.5)  - dual(ur(0.5),ud(0.5))) < 10eps()
        end
    end

    @testset "Eig test #336" begin
        d = 0..π
        A=Derivative(d)^2
        λ=eigvals(Dirichlet(d),A,100)
        @test sort(λ)[end-5:end] ≈ -(-6:-1).^2
        λ=eigvals(Dirichlet(), A, 100)
        @test sort(λ)[end-5:end] ≈ -(-6:-1).^2

        F = x->x^8
        d = Interval(0.0,1.0)
        f = Fun(F,d)
        ginf = Fun(x->exp(-x),d)
        gp = ginf'
        Af = Fun(x->x+f(x),d)
        transport_ = Fun(x-> x - 1,d)
        damping = Fun(x-> 1 - f(x),d)
        A = transport_*Derivative(d) + damping
        P = -DefiniteIntegral(Chebyshev(d))[LowRankFun((x,y)->gp(x)*(y+f(y)),d^2)];
        λ,V = eigs(A,100)
        @test norm(sort(real(filter(x->isreal(x),λ)))[1:5]-(0:4)) ≤ 100000eps()

        λ,V = eigs(A+P,100)
        @test sort(real(filter(x->isreal(x),λ)))[5] ≈ 3.93759261234502 atol=1E-3
    end

    @testset "Sampling" begin
        ff=(x,y)->(x-y)^2*exp(-x^2/2-y^2/2)

        f=Fun(ff,Domain(-4..4)^2)
        r=ApproxFun.sample(f,5000)

        #We can compare the histogram to the 1-point correlation
        g=sum(f,1)/sum(f)
        @test  g(0.1) ≈ 0.2004758624973169

        # check bisection inv
        f = Fun(x -> exp(-x^2/2),-5..5)
        g = cumsum(f)
        @test g(ApproxFun.bisectioninv(g,0.5)) ≈ 0.5
    end

    @testset "findmin/max" begin
        f = Fun(x -> exp(0.25x) + sin(x) + 0.5cos(10x), -4..4)
        @test [findmax(f)...] ≈ [3.0531164509549584, 1.886754631165656]
        @test [findmin(f)...] ≈ [-0.825047261209411, -1.5741041425422948]
    end
end


text:
module SPEEDY

# All these functions need to be updated to accommodate the SPEEDY model

struct Matrices{T,M<:Array{T,3}}
    u0::M
    v0::M
    T0::M
    q0::M
    ps0::M
    rain0::M
end
function setup(nx::Int,
               ny::Int,
               nz::Int,
               T::DataType = Float64)
    # Memory allocation
    u = ones(T, nx, ny, nz)
    # u = reshape(u, nx, ny, nz)
    v = ones(T, nx, ny, nz)
    # v = reshape(v, nx, ny, nz)
    Temp = ones(T, nx, ny, nz)
    # Temp = reshape(Temp, nx, ny, nz)
    q = ones(T, nx, ny, nz)
    # q = reshape(q, nx, ny, nz)
    ps = ones(T, nx, ny, nz)
    rain = ones(T, nx, ny, nz)

    return Matrices(u,v,Temp,q,ps,rain)
end
end# module

max_stars_repo_name: Team-RADDISH/ParticleDA-UseCases

max_issues_repo_name: Team-RADDISH/ParticleDA-UseCases

max_forks_repo_name: Team-RADDISH/ParticleDA-UseCases


text:
################################################################################
#                              BASIC MODEL DEFINITION
################################################################################
"""
    TranscriptionData

A DataType for storing the data mapping an [`InfiniteOpt.InfiniteModel`](@ref)
that has been transcribed to a regular [`JuMP.Model`](@ref) that contains the
transcribed variables. This is stored in the `ext` field of a `JuMP.Model` to
make what is called a `TranscriptionModel` via the [`TranscriptionModel`](@ref)
constructor.

**Fields**
- `infvar_lookup::Dict{InfiniteOpt.GeneralVariableRef, Dict{Vector{Float64}, Int}}`:
   A lookup table of infinite variable transcriptions via support value.
- `infvar_mappings::Dict{InfiniteOpt.GeneralVariableRef, Vector{JuMP.VariableRef}}`:
   Map infinite variables to their transcription variables.
- `infvar_supports::Dict{InfiniteOpt.GeneralVariableRef, Vector{Tuple}}`:
   Map infinite variables to their support values.
- `infvar_support_labels::Dict{InfiniteOpt.GeneralVariableRef, Vector{Set{DataType}}}`:
   Map the infinite variables to their support labels.
- `finvar_mappings::Dict{InfiniteOpt.GeneralVariableRef, JuMP.VariableRef}`:
   Map finite variables to their transcription variables.
- `reduced_vars::Vector{InfiniteOpt.ReducedVariable{InfiniteOpt.GeneralVariableRef}}`:
   Store the core reduced variable objects of reduced variables formed on transcription.
- `last_point_index::Int`: The last internal point variable index added.
- `measure_lookup::Dict{InfiniteOpt.GeneralVariableRef, Dict{Vector{Float64}, Int}}`:
   A lookup table of measure transcriptions via support value.
- `measure_mappings::Dict{InfiniteOpt.GeneralVariableRef, Vector{JuMP.AbstractJuMPScalar}}`:
   Map measures to transcription expressions.
- `measure_supports::Dict{InfiniteOpt.GeneralVariableRef, Vector{Tuple}}`:
   Map measures to their supports values (if the transcribed measure is still infinite).
- `measure_support_labels::Dict{InfiniteOpt.GeneralVariableRef, Vector{Set{DataType}}}`:
   Map measures to their support labels if they have any.
- `constr_mappings::Dict{InfiniteOpt.InfOptConstraintRef, Vector{JuMP.ConstraintRef}}`:
   Map constraints to their transcriptions.
- `constr_supports::Dict{InfiniteOpt.InfOptConstraintRef, Vector{Tuple}}`:
   Map constraints to their support values.
- `constr_support_labels::Dict{InfiniteOpt.InfOptConstraintRef, Vector{Set{DataType}}}`:
   Map constraints to their support labels.
- `supports::Tuple`: Store the collected parameter supports here.
- `support_labels::Tuple`: Store the collected parameter labels here.
- `has_internal_supports::Bool`: Where any internal supports collected?
"""
mutable struct TranscriptionData
    # Variable information
    infvar_lookup::Dict{InfiniteOpt.GeneralVariableRef, Dict{Vector{Float64}, Int}}
    infvar_mappings::Dict{InfiniteOpt.GeneralVariableRef, Vector{JuMP.VariableRef}}
    infvar_supports::Dict{InfiniteOpt.GeneralVariableRef, Vector{Tuple}}
    infvar_support_labels::Dict{InfiniteOpt.GeneralVariableRef, Vector{Set{DataType}}}
    finvar_mappings::Dict{InfiniteOpt.GeneralVariableRef, JuMP.VariableRef}

    # Internal variables (created via internal measure expansions)
    reduced_vars::Vector{InfiniteOpt.ReducedVariable{InfiniteOpt.GeneralVariableRef}}
    last_point_index::Int

    # Measure information
    measure_lookup::Dict{InfiniteOpt.GeneralVariableRef, Dict{Vector{Float64}, Int}}
    measure_mappings::Dict{InfiniteOpt.GeneralVariableRef, Vector{JuMP.AbstractJuMPScalar}}
    measure_supports::Dict{InfiniteOpt.GeneralVariableRef, Vector{Tuple}}
    measure_support_labels::Dict{InfiniteOpt.GeneralVariableRef, Vector{Set{DataType}}}

    # Constraint information
    constr_mappings::Dict{InfiniteOpt.InfOptConstraintRef,
                          Vector{JuMP.ConstraintRef}}
    constr_supports::Dict{InfiniteOpt.InfOptConstraintRef,
                          Vector{Tuple}}
    constr_support_labels::Dict{InfiniteOpt.InfOptConstraintRef,
                               Vector{Set{DataType}}}

    # Collected Supports
    supports::Tuple
    support_labels::Tuple
    has_internal_supports::Bool

    # Default constructor
    function TranscriptionData()
        return new( # variable info
                   Dict{InfiniteOpt.GeneralVariableRef, Dict{Vector{Float64}, Int}}(),
                   Dict{InfiniteOpt.GeneralVariableRef, Vector{JuMP.VariableRef}}(),
                   Dict{InfiniteOpt.GeneralVariableRef, Vector{Tuple}}(),
                   Dict{InfiniteOpt.GeneralVariableRef, Vector{Set{DataType}}}(),
                   Dict{InfiniteOpt.GeneralVariableRef, JuMP.VariableRef}(),
                   # internal variables
                   Vector{InfiniteOpt.ReducedVariable{InfiniteOpt.GeneralVariableRef}}(),
                   0,
                   # measure info
                   Dict{InfiniteOpt.GeneralVariableRef, Dict{Vector{Float64}, Int}}(),
                   Dict{InfiniteOpt.GeneralVariableRef, Vector{JuMP.AbstractJuMPScalar}}(),
                   Dict{InfiniteOpt.GeneralVariableRef, Vector{Tuple}}(),
                   Dict{InfiniteOpt.GeneralVariableRef, Vector{Set{DataType}}}(),
                   # constraint info
                   Dict{InfiniteOpt.InfOptConstraintRef, Vector{JuMP.ConstraintRef}}(),
                   Dict{InfiniteOpt.InfOptConstraintRef, Vector{Vector{Float64}}}(),
                   Dict{InfiniteOpt.InfOptConstraintRef, Vector{Set{DataType}}}(),
                   # support storage
                   (), (), false)
    end
end

"""
    TranscriptionModel([optimizer_constructor;
                       caching_mode::MOIU.CachingOptimizerMode = MOIU.AUTOMATIC,
                       bridge_constraints::Bool = true])::JuMP.Model

Return a [`JuMP.Model`](@ref) with [`TranscriptionData`](@ref) included in the
`ext` data field. Accepts the same arguments as a typical JuMP `Model`.
More detailed variable and constraint naming can be enabled via `verbose_naming`.

**Example**
```julia-repl
julia> TranscriptionModel()
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
```
"""
function TranscriptionModel(; kwargs...)::JuMP.Model
    model = JuMP.Model(; kwargs...)
    model.ext[:TransData] = TranscriptionData()
    return model
end
# Accept optimizer constructors
function TranscriptionModel(optimizer_constructor;
                            kwargs...)::JuMP.Model
    model = JuMP.Model(optimizer_constructor; kwargs...)
    model.ext[:TransData] = TranscriptionData()
    return model
end

################################################################################
#                                BASIC QUERIES
################################################################################
"""
    is_transcription_model(model::JuMP.Model)::Bool

Return true if `model` is a `TranscriptionModel` or false otherwise.

**Example**
```julia-repl
julia> is_transcription_model(model)
true
```
"""
function is_transcription_model(model::JuMP.Model)::Bool
    return haskey(model.ext, :TransData)
end

"""
    transcription_data(model::JuMP.Model)::TranscriptionData

Return the `TranscriptionData` from a `TranscriptionModel`. Errors if it is not
a `TranscriptionModel`.
"""
function transcription_data(model::JuMP.Model)::TranscriptionData
    !is_transcription_model(model) && error("Model is not a transcription model.")
    return model.ext[:TransData]
end

"""
    has_internal_supports(model::JuMP.Model)::Bool

Return a `Bool` whether `model` has any internal supports that were collected.
"""
function has_internal_supports(model::JuMP.Model)::Bool
    return transcription_data(model).has_internal_supports
end

################################################################################
#                              VARIABLE QUERIES
################################################################################
# Define method for checking if the label needs to be accounted for
function _ignore_label(model::JuMP.Model, label::Type{<:InfiniteOpt.AbstractSupportLabel})::Bool
    return label == InfiniteOpt.All || (!has_internal_supports(model) && label == InfiniteOpt.PublicLabel)
end

"""
    transcription_variable(model::JuMP.Model,
        vref::InfiniteOpt.GeneralVariableRef;
        [label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
         ndarray::Bool = false])

Return the transcribed variable reference(s) corresponding to `vref`. Errors
if no transcription variable is found. Also can query via the syntax:
```julia
transcription_variable(vref::InfiniteOpt.GeneralVariableRef;
    [label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
     ndarray::Bool = false])
```
If the infinite model contains a built transcription model. By default, this
method returns only transcribed variables associated with public supports. All the
variables can be returned by setting `label = All`.

If `vref` is infinite and `ndarray = true` then an n-dimensional array will be
returned in accordance with the infinite parameters that have unique object
numbers. In this case, `label` will be used to search the intersection of variable
supports that use the label. This is defers from the default behavior which
considers the union.

**Example**
```julia-repl
julia> transcription_variable(trans_model, infvar)
2-element Array{VariableRef,1}:
 infvar(support: 1)
 infvar(support: 2)

julia> transcription_variable(trans_model, hdvar)
hdvar

julia> transcription_variable(infvar)
2-element Array{VariableRef,1}:
 infvar(support: 1)
 infvar(support: 2)

julia> transcription_variable(hdvar)
hdvar
```
"""
function transcription_variable(model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef;
    label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
    ndarray::Bool = false
    )
    return transcription_variable(model, vref, InfiniteOpt._index_type(vref), label, ndarray)
end

# define convenient aliases
const InfVarIndex = Union{InfiniteOpt.InfiniteVariableIndex,
                          InfiniteOpt.ReducedVariableIndex,
                          InfiniteOpt.DerivativeIndex}
const FinVarIndex = Union{InfiniteOpt.HoldVariableIndex,
                          InfiniteOpt.PointVariableIndex}

## Define the variable mapping functions
# FinVarIndex
function transcription_variable(model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{V},
    label::Type{<:InfiniteOpt.AbstractSupportLabel},
    ndarray::Bool
    )::JuMP.VariableRef where {V <: FinVarIndex}
    var = get(transcription_data(model).finvar_mappings, vref, nothing)
    if var === nothing
        error("Variable reference $vref not used in transcription model.")
    end
    return var
end

# InfVarIndex
function transcription_variable(model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{V},
    label::Type{<:InfiniteOpt.AbstractSupportLabel},
    ndarray::Bool
    )::Array{JuMP.VariableRef} where {V <: InfVarIndex}
    vars = get(transcription_data(model).infvar_mappings, vref, nothing)
    if vars === nothing
        error("Variable reference $vref not used in transcription model.")
    end
    if ndarray
        return make_ndarray(model, vref, vars, label)
    elseif _ignore_label(model, label)
        return vars
    else
        labels = transcription_data(model).infvar_support_labels[vref]
        inds = map(s -> any(l -> l <: label, s), labels)
        return vars[inds]
    end
end

# ParameterFunctionIndex
function transcription_variable(model::JuMP.Model,
    fref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{InfiniteOpt.ParameterFunctionIndex},
    label::Type{<:InfiniteOpt.AbstractSupportLabel},
    ndarray::Bool
    )
    # get the object numbers of the expression and form the support iterator
    obj_nums = InfiniteOpt._object_numbers(fref)
    support_indices = support_index_iterator(model, obj_nums)
    vals = Vector{Float64}(undef, length(support_indices))
    check_labels = length(vals) > 1 && !_ignore_label(model, label)
    label_inds = ones(Bool, length(vals))
    # iterate over the indices and compute the values
    for (i, idx) in enumerate(support_indices)
        supp = index_to_support(model, idx)
        if check_labels && !any(l -> l <: label, index_to_labels(model, idx))
            @inbounds label_inds[i] = false
        end
        @inbounds vals[i] = transcription_expression(model, fref, supp)
    end
    # return the values
    if ndarray
        return make_ndarray(model, fref, vals, label)
    else
        return vals[label_inds]
    end
end

# Fallback
function transcription_variable(model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    index_type,
    label,
    ndarray)
    error("`transcription_variable` not defined for variables with indices of " *
          "type $(index_type) and/or is not defined for labels of type $(label).")
end

# Dispatch for internal models
function transcription_variable(vref::InfiniteOpt.GeneralVariableRef;
    label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
    ndarray::Bool = false
    )
    trans_model = InfiniteOpt.optimizer_model(JuMP.owner_model(vref))
    return transcription_variable(trans_model, vref, label = label,
                                  ndarray = ndarray)
end

"""
    InfiniteOpt.optimizer_model_variable(vref::InfiniteOpt.GeneralVariableRef,
        ::Val{:TransData};
        [label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
        ndarray::Bool = false])

Proper extension of [`InfiniteOpt.optimizer_model_variable`](@ref) for
`TranscriptionModel`s. This simply dispatches to [`transcription_variable`](@ref).
"""
function InfiniteOpt.optimizer_model_variable(vref::InfiniteOpt.GeneralVariableRef,
    ::Val{:TransData};
    label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
    ndarray::Bool = false)
    return transcription_variable(vref, label = label, ndarray = ndarray)
end

"""
    InfiniteOpt.variable_supports(model::JuMP.Model,
        vref::InfiniteOpt.DecisionVariableRef,
        key::Val{:TransData} = Val(:TransData);
        [label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
        ndarray::Bool = false])

Return the support alias mapping associated with `vref` in the transcription model.
Errors if `vref` does not have transcripted variables. See `transcription_variable`
for an explanation of `ndarray`.
"""
function InfiniteOpt.variable_supports(model::JuMP.Model,
    dvref::Union{InfiniteOpt.InfiniteVariableRef, InfiniteOpt.ReducedVariableRef,
                 InfiniteOpt.DerivativeRef},
    key::Val{:TransData} = Val(:TransData);
    label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
    ndarray::Bool = false
    )::Array
    vref = InfiniteOpt._make_variable_ref(JuMP.owner_model(dvref), JuMP.index(dvref))
    if !haskey(transcription_data(model).infvar_mappings, vref)
        error("Variable reference $vref not used in transcription model.")
    elseif !haskey(transcription_data(model).infvar_supports, vref)
        prefs = InfiniteOpt.raw_parameter_refs(dvref)
        lookups = transcription_data(model).infvar_lookup[vref]
        type = typeof(Tuple(first(keys(lookups)), prefs))
        supps = Vector{type}(undef, length(lookups))
        for (s, i) in lookups
            supps[i] = Tuple(s, prefs)
        end
        transcription_data(model).infvar_supports[vref] = supps
    end
    supps = transcription_data(model).infvar_supports[vref]
    if ndarray
        return make_ndarray(model, dvref, supps, label)
    elseif _ignore_label(model, label)
        return supps
    else
        labels = transcription_data(model).infvar_support_labels[vref]
        inds = map(s -> any(l -> l <: label, s), labels)
        return supps[inds]
    end
end

# ParameterFunctionRef
function InfiniteOpt.variable_supports(model::JuMP.Model,
    dvref::InfiniteOpt.ParameterFunctionRef,
    key::Val{:TransData} = Val(:TransData);
    label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
    ndarray::Bool = false
    )::Array
    # get the object numbers of the expression and form the support iterator
    obj_nums = sort(InfiniteOpt._object_numbers(dvref))
    support_indices = support_index_iterator(model, obj_nums)
    supps = Vector{Tuple}(undef, length(support_indices))
    check_labels = length(supps) > 1 && !_ignore_label(model, label)
    param_supps = parameter_supports(model)
    label_inds = ones(Bool, length(supps))
    # iterate over the indices and compute the values
    for (i, idx) in enumerate(support_indices)
        if check_labels && !any(l -> l <: label, index_to_labels(model, idx))
            @inbounds label_inds[i] = false
        end
        @inbounds supps[i] = Tuple(param_supps[j][idx[j]] for j in obj_nums)
    end
    # return the supports
    if ndarray
        return make_ndarray(model, dvref, supps, label)
    else
        return supps[label_inds]
    end
end

"""
    lookup_by_support(model::JuMP.Model,
                      vref::InfiniteOpt.GeneralVariableRef,
                      support::Vector)

Return the transcription expression of `vref` defined at its `support`. This is
intended as a helper method for automated transcription.
"""
function lookup_by_support(model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    support::Vector)
    return lookup_by_support(model, vref, InfiniteOpt._index_type(vref), support)
end

# define error function for not being able to find a variable by its support
_supp_error() = error("Unable to locate transcription variable by support, consider " *
                      "rebuilding the infinite model with less significant digits. " *
                      "Note this might be due to partially evaluating dependent parameters " *
                      "which is not supported by TranscriptionOpt. Such is the case with " *
                      "derivatives/measures that dependent on single dependent parameters.")

# InfiniteIndex
function lookup_by_support(model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{V},
    support::Vector
    )::JuMP.VariableRef where {V <: InfVarIndex}
    if !haskey(transcription_data(model).infvar_lookup, vref)
        error("Variable reference $vref not used in transcription model.")
    end
    idx = get(_supp_error, transcription_data(model).infvar_lookup[vref], support)
    return transcription_data(model).infvar_mappings[vref][idx]
end

# ParameterFunctionIndex
function lookup_by_support(model::JuMP.Model,
    fref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{InfiniteOpt.ParameterFunctionIndex},
    support::Vector
    )
    prefs = InfiniteOpt.raw_parameter_refs(fref)
    func = InfiniteOpt.raw_function(fref)
    return func(Tuple(support, prefs)...)
end

# FiniteVariableIndex
function lookup_by_support(model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{V},
    support::Vector
    )::JuMP.VariableRef where {V <: FinVarIndex}
    if !haskey(transcription_data(model).finvar_mappings, vref)
        error("Variable reference $vref not used in transcription model.")
    end
    return transcription_data(model).finvar_mappings[vref]
end

"""
    InfiniteOpt.internal_reduced_variable(
        vref::InfiniteOpt.ReducedVariableRef,
        ::Val{:TransData}
        )::InfiniteOpt.ReducedVariable{InfiniteOpt.GeneralVariableRef}

Return the internal reduced variable associated with `vref`, assuming it was
added internally during measure expansion at the transcription step. This
extends [`InfiniteOpt.internal_reduced_variable`](@ref) as described in its
docstring. Errors, if no such variable can be found.
"""
function InfiniteOpt.internal_reduced_variable(
    vref::InfiniteOpt.ReducedVariableRef,
    ::Val{:TransData}
    )::InfiniteOpt.ReducedVariable{InfiniteOpt.GeneralVariableRef}
    trans_model = InfiniteOpt.optimizer_model(JuMP.owner_model(vref))
    reduced_vars = transcription_data(trans_model).reduced_vars
    idx = -1 * JuMP.index(vref).value
    if idx in keys(reduced_vars)
        return reduced_vars[idx]
    else
        error("Invalid reduced variable reference, this likely is attributed " *
              "to its being deleted.")
    end
end

################################################################################
#                              MEASURE QUERIES
################################################################################
# MeasureIndex
function transcription_variable(model::JuMP.Model,
    mref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{InfiniteOpt.MeasureIndex},
    label::Type{<:InfiniteOpt.AbstractSupportLabel},
    ndarray::Bool = false
    )
    exprs = get(transcription_data(model).measure_mappings, mref, nothing)
    if exprs === nothing
        error("Measure reference $mref not used in transcription model.")
    end
    if ndarray
        return make_ndarray(model, mref, exprs, label)
    elseif length(exprs) > 1 && _ignore_label(model, label)
        return exprs
    elseif length(exprs) > 1
        labels = transcription_data(model).measure_support_labels[mref]
        inds = map(s -> any(l -> l <: label, s), labels)
        return exprs[inds]
    else
        return first(exprs)
    end
end

# Extend transcription_expression
function lookup_by_support(model::JuMP.Model,
    mref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{InfiniteOpt.MeasureIndex},
    support::Vector
    )::JuMP.AbstractJuMPScalar
    if !haskey(transcription_data(model).measure_lookup, mref)
        error("Measure reference $mref not used in transcription model.")
    end
    idx = get(_supp_error, transcription_data(model).measure_lookup[mref], support)
    return transcription_data(model).measure_mappings[mref][idx]
end

# Extend variable_supports
function InfiniteOpt.variable_supports(model::JuMP.Model,
    dmref::InfiniteOpt.MeasureRef,
    key::Val{:TransData} = Val(:TransData);
    label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
    ndarray::Bool = false
    )
    mref = InfiniteOpt._make_variable_ref(JuMP.owner_model(dmref), JuMP.index(dmref))
    if !haskey(transcription_data(model).measure_mappings, mref)
        error("Measure reference $mref not used in transcription model.")
    elseif !haskey(transcription_data(model).measure_supports, mref)
        lookups = transcription_data(model).measure_lookup[mref]
        prefs = InfiniteOpt.parameter_refs(dmref)
        vt_prefs = InfiniteOpt.Collections.VectorTuple(prefs)
        type = typeof(Tuple(first(keys(lookups)), vt_prefs))
        supps = Vector{type}(undef, length(lookups))
        for (supp, i) in lookups
            supps[i] = Tuple(supp, vt_prefs)
        end
        transcription_data(model).measure_supports[mref] = supps
    end
    supps = transcription_data(model).measure_supports[mref]
    if ndarray
        return make_ndarray(model, dmref, supps, label)
    elseif length(supps) > 1 && _ignore_label(model, label)
        return supps
    elseif length(supps) > 1
        labels = transcription_data(model).measure_support_labels[mref]
        inds = map(s -> any(l -> l <: label, s), labels)
        return supps[inds]
    else
        return first(supps)
    end
end

################################################################################
#                             EXPRESSION QUERIES
################################################################################
"""
    transcription_expression(model::JuMP.Model,
        expr::JuMP.AbstractJuMPScalar;
        [label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
        ndarray::Bool = false])

Return the transcribed expression(s) corresponding to `expr`. Errors
if `expr` cannot be transcribed. Also can query via the syntax:
```julia
transcription_expression(expr::JuMP.AbstractJuMPScalar;
                         [label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
                         ndarray::Bool = false])
```
If the infinite model contains a built transcription model. By default, this
method returns only transcribed expressions associated with public supports. All the
expressions can be returned by setting `label = All`.

If `expr` is infinite and `ndarray = true` then an n-dimensional array will be
returned in accordance with the infinite parameters that have unique object
numbers. In this case, `label` will be used to search the intersection of the
supports that use the label. This is defers from the default behavior which
considers the union.

**Example**
```julia-repl
julia> transcription_expression(trans_model, my_expr)
x(support: 1) - y

julia> transcription_expression(my_expr)
x(support: 1) - y
```
"""
function transcription_expression(model::JuMP.Model,
    expr::Union{JuMP.GenericAffExpr, JuMP.GenericQuadExpr};
    label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
    ndarray::Bool = false
    )
    # get the object numbers of the expression and form the support iterator
    obj_nums = InfiniteOpt._object_numbers(expr)
    support_indices = support_index_iterator(model, obj_nums)
    exprs = Vector{JuMP.AbstractJuMPScalar}(undef, length(support_indices))
    check_labels = length(exprs) > 1 && !_ignore_label(model, label)
    label_inds = ones(Bool, length(exprs))
    # iterate over the indices and compute the values
    for (i, idx) in enumerate(support_indices)
        supp = index_to_support(model, idx)
        if check_labels && !any(l -> l <: label, index_to_labels(model, idx))
            @inbounds label_inds[i] = false
        end
        @inbounds exprs[i] = transcription_expression(model, expr, supp)
    end
    # return the expressions
    if ndarray
        return make_ndarray(model, expr, exprs, label)
    else
        exprs = exprs[label_inds]
        return length(support_indices) > 1 ? exprs : first(exprs)
    end
end

# Define for variables
function transcription_expression(model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef;
    label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
    ndarray::Bool = false)
    return transcription_variable(model, vref, label = label, ndarray = ndarray)
end

# Dispatch for internal models
function transcription_expression(expr::JuMP.AbstractJuMPScalar;
    label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
    ndarray::Bool = false
    )
    model = InfiniteOpt._model_from_expr(expr)
    if model === nothing
        return zero(JuMP.AffExpr) + JuMP.constant(expr)
    else
        trans_model = InfiniteOpt.optimizer_model(model)
    end
    return transcription_expression(trans_model, expr, label = label,
                                    ndarray = ndarray)
end

"""
    InfiniteOpt.optimizer_model_expression(expr::JuMP.AbstractJuMPScalar,
        ::Val{:TransData};
        [label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
        ndarray::Bool = false])

Proper extension of [`InfiniteOpt.optimizer_model_expression`](@ref) for
`TranscriptionModel`s. This simply dispatches to [`transcription_expression`](@ref).
"""
function InfiniteOpt.optimizer_model_expression(
    expr::Union{JuMP.GenericAffExpr, JuMP.GenericQuadExpr},
    ::Val{:TransData};
    label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
    ndarray::Bool = false)
    return transcription_expression(expr, label = label, ndarray = ndarray)
end

"""
    InfiniteOpt.expression_supports(model::JuMP.Model,
        expr::JuMP.AbstractJuMPScalar,
        key::Val{:TransData} = Val(:TransData);
        [label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
        ndarray::Bool = false])

Return the support alias mappings associated with `expr`. Errors if `expr` cannot
be transcribed.
"""
function InfiniteOpt.expression_supports(model::JuMP.Model,
    expr::Union{JuMP.GenericAffExpr, JuMP.GenericQuadExpr},
    key::Val{:TransData} = Val(:TransData);
    label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
    ndarray::Bool = false
    )
    # get the object numbers of the expression and form the support iterator
    obj_nums = sort(InfiniteOpt._object_numbers(expr))
    support_indices = support_index_iterator(model, obj_nums)
    supps = Vector{Tuple}(undef, length(support_indices))
    check_labels = length(supps) > 1 && !_ignore_label(model, label)
    param_supps = parameter_supports(model)
    label_inds = ones(Bool, length(supps))
    # iterate over the indices and compute the values
    for (i, idx) in enumerate(support_indices)
        if check_labels && !any(l -> l <: label, index_to_labels(model, idx))
            @inbounds label_inds[i] = false
        end
        @inbounds supps[i] = Tuple(param_supps[j][idx[j]] for j in obj_nums)
    end
    # return the supports
    if ndarray
        return make_ndarray(model, expr, supps, label)
    else
        supps = supps[label_inds]
        return length(support_indices) > 1 ? supps : first(supps)
    end
end

################################################################################
#                             CONSTRAINT QUERIES
################################################################################
"""
    transcription_constraint(model::JuMP.Model,
        cref::InfiniteOpt.InfOptConstraintRef;
        [label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
        ndarray::Bool = false])

Return the transcribed constraint reference(s) corresponding to `cref`. Errors
if `cref` has not been transcribed. Also can query via the syntax:
```julia
transcription_constraint(cref::InfiniteOpt.InfOptConstraintRef;
    [label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
    ndarray::Bool = false])
```
If the infinite model contains a built transcription model. By default, this
method returns only transcribed constraints associated with public supports. All the
constraints can be returned by setting `label = All`.

If `cref` is infinite and `ndarray = true` then an n-dimensional array will be
returned in accordance with the infinite parameters that have unique object
numbers. In this case, `label` will be used to search the intersection of the
supports that use the label. This is defers from the default behavior which
considers the union.

**Example**
```julia-repl
julia> transcription_constraint(trans_model, fin_con)
fin_con : x(support: 1) - y <= 3.0

julia> transcription_constraint(fin_con)
fin_con : x(support: 1) - y <= 3.0
```
"""
function transcription_constraint(model::JuMP.Model,
    cref::InfiniteOpt.InfOptConstraintRef;
    label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
    ndarray::Bool = false
    )
    constr = get(transcription_data(model).constr_mappings, cref, nothing)
    if constr === nothing
      error("Constraint reference $cref not used in transcription model.")
    end
    if ndarray
        return make_ndarray(model, cref, constr, label)
    elseif length(constr) > 1 && _ignore_label(model, label)
        return constr
    elseif length(constr) > 1
        labels = transcription_data(model).constr_support_labels[cref]
        inds = map(s -> any(l -> l <: label, s), labels)
        return constr[inds]
    else
        return first(constr)
    end
end

# Dispatch for internal models
function transcription_constraint(cref::InfiniteOpt.InfOptConstraintRef;
    label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
    ndarray::Bool = false
    )
    trans_model = InfiniteOpt.optimizer_model(JuMP.owner_model(cref))
    return transcription_constraint(trans_model, cref, label = label,
                                    ndarray = ndarray)
end

"""
    InfiniteOpt.optimizer_model_constraint(cref::InfiniteOpt.InfOptConstraintRef,
        ::Val{:TransData};
        [label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel
        ndarray::Bool = false])

Proper extension of [`InfiniteOpt.optimizer_model_constraint`](@ref) for
`TranscriptionModel`s. This simply dispatches to [`transcription_constraint`](@ref).
"""
function InfiniteOpt.optimizer_model_constraint(
    cref::InfiniteOpt.InfOptConstraintRef,
    ::Val{:TransData};
    label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
    ndarray::Bool = false
    )
    return transcription_constraint(cref, label = label, ndarray = ndarray)
end

"""
    InfiniteOpt.constraint_supports(model::JuMP.Model,
        cref::InfiniteOpt.InfOptConstraintRef,
        key::Val{:TransData} = Val(:TransData);
        [label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
        ndarray::Bool = false])

Return the support alias mappings associated with `cref`. Errors if `cref` is
not transcribed.
"""
function InfiniteOpt.constraint_supports(model::JuMP.Model,
    cref::InfiniteOpt.InfOptConstraintRef,
    key::Val{:TransData} = Val(:TransData);
    label::Type{<:InfiniteOpt.AbstractSupportLabel} = InfiniteOpt.PublicLabel,
    ndarray::Bool = false
    )
    supps = get(transcription_data(model).constr_supports, cref, nothing)
    if supps === nothing
        error("Constraint reference $cref not used in transcription model.")
    end
    if ndarray
        return make_ndarray(model, cref, supps, label)
    elseif length(supps) > 1 && _ignore_label(model, label)
        return supps
    elseif length(supps) > 1
        labels = transcription_data(model).constr_support_labels[cref]
        inds = map(s -> any(l -> l <: label, s), labels)
        return supps[inds]
    else
        return first(supps)
    end
end

################################################################################
#                             OTHER QUERIES
################################################################################
"""
    parameter_supports(model::JuMP.Model)::Tuple

Return the collected parameter support tuple that is stored in
`TranscriptionData.supports`.
"""
function parameter_supports(model::JuMP.Model)::Tuple
    return transcription_data(model).supports
end

"""
    support_index_iterator(model::JuMP.Model, [obj_nums::Vector{Int}])::CartesianIndices

Return the `CartesianIndices` that determine the indices of the unique combinations
of `TranscriptionData.supports` stored in `model`. If `obj_nums` is specified,
then the indices will only include the tuple elements uses indices are included
in the object numbers `obj_nums` and all others will be assigned the last index
which should correspond to an appropriately sized placeholder comprised of `NaN`s.
Note this method assumes that [`set_parameter_supports`](@ref) has already been
called and that the last elements of each support vector contains a placeholder
value.
"""
function support_index_iterator(model::JuMP.Model)::CartesianIndices
    raw_supps = parameter_supports(model)
    return CartesianIndices(ntuple(i -> 1:length(raw_supps[i])-1, length(raw_supps)))
end

# Generate for a subset of object numbers (use last index as placeholder --> support with NaNs)
function support_index_iterator(model::JuMP.Model,
                                obj_nums::Vector{Int})::CartesianIndices
    raw_supps = parameter_supports(model)
    lens = map(i -> length(i), raw_supps)
    # prepare the indices of each support combo
    # note that the actual supports are afrom 1:length-1 and the placeholders are at the ends
    return CartesianIndices(ntuple(i -> i in obj_nums ? (1:lens[i]-1) : (lens[i]:lens[i]),
                                   length(raw_supps)))
end

"""
    index_to_support(model::JuMP.Model, index::CartesianIndex)::Vector{Float64}

Given a particular support `index` generated via [`support_index_iterator`](@ref)
using `model`, return the corresponding support from `TranscriptionData.supports`
using placeholder `NaN`s as appropriate for tuple elements that are unneeded.
"""
function index_to_support(model::JuMP.Model,
                          index::CartesianIndex)::Vector{Float64}
    raw_supps = parameter_supports(model)
    return [j for i in eachindex(index.I) for j in raw_supps[i][index[i]]]
end

"""
    index_to_labels(model::JuMP.Model, index::CartesianIndex)::Set{DataType}

Given a particular support `index` generated via [`support_index_iterator`](@ref)
using `model`, return the corresponding support label set from `TranscriptionData.support_labels`.
"""
function index_to_labels(model::JuMP.Model,
                         index::CartesianIndex)::Set{DataType}
    raw_labels = transcription_data(model).support_labels
    labels = Set{DataType}()
    for (i, j) in enumerate(index.I)
        union!(labels, raw_labels[i][j])
    end
    return labels
end

################################################################################
#                               QUERY FORMATERS
################################################################################
# Helper function for getting the array type T
function _get_array_type(array::Array{T, N}) where {T, N}
    return T
end

## Helper functions to consistently get object numbers
# Fallback
function _get_object_numbers(ref)
    return InfiniteOpt._object_numbers(ref)
end

# Expressions
function _get_object_numbers(expr::Union{JuMP.GenericAffExpr, JuMP.GenericQuadExpr})
    return sort(InfiniteOpt._object_numbers(expr))
end

"""
    make_narray(model::JuMP.Model,
                ref::Union{JuMP.AbstractJuMPScalar, InfiniteOpt.InfOptConstraintRef},
                info::Vector,
                label::Type{<:InfiniteOpt.AbstractSupportLabel})::Array

Take the results`info` associated with `ref` and rearrange them into an
n-dimensional array where the axes correspond to the infinite parameter dependencies
in accordance with their creation. Note that this works by querying the object
numbers. Thus, independent infinite parameters will each get their own dimension
(even if they are defined at the same time in an array) and each dependent infinite
parameter group will have its own dimension.
"""
function make_ndarray(model::JuMP.Model, ref, info::Vector, label::DataType)
    # get the object numbers
    obj_nums = _get_object_numbers(ref)
    # return result if it is from a finite object
    if isempty(obj_nums)
        return info
    end
    # determine the dimensions of the new array
    raw_supps = parameter_supports(model)
    dims = Tuple(length(raw_supps[i]) - 1 for i in eachindex(raw_supps) if i in obj_nums)
    # check that the lengths match (otherwise we'll have some sparse set)
    # TODO add capability to avoid this problem (make reduced array by looking at the supports)
    if length(info) != prod(dims)
        error("Unable to make `ndarray`. This is likely due to the object being " *
              "over a portion of the infinite-domain (e.g., bounded constraints and " *
              "certain reduced infinite variables.")
    end
    # make and populate the array
    narray = Array{_get_array_type(info)}(undef, dims)
    for (i, idx) in enumerate(eachindex(narray))
        narray[idx] = info[i]
    end
    # rearrange the array as needed to match the object number order
    sorted_array = issorted(obj_nums) ? narray : permutedims(narray, sortperm(obj_nums))
    # consider the label specified (this will enforce the intersection of labels)
    if _ignore_label(model, label)
        return sorted_array
    else
        labels = transcription_data(model).support_labels[obj_nums]
        inds = map(sets -> findall(s -> any(l -> l <: label, s), sets), labels)
        return sorted_array[inds...]
    end
end


text:
module Geom

using Colors
using Compose
using DataStructures
using Distributions
using Gadfly
using Measures
using IndirectArrays
using Base.Iterators

import Compose.combine # Prevent DataFrame.combine from taking over.
import Gadfly: render, layers, element_aesthetics, inherit, escape_id,
               default_statistic, default_scales, element_coordinate_type,
               ScaleElement, svg_color_class_from_label, isconcrete,
               concretize, discretize_make_ia, get_stroke_vector, parse_colorant
import IterTools: takestrict

const empty_tag = Symbol("")

function subtags(parent_tag, suffixes...)
    if parent_tag == empty_tag
        return map(s->empty_tag, suffixes)
    end
    return map(s->Symbol(parent_tag, "#", s), suffixes)
end

# Geometry that renders nothing.
struct Nil <: Gadfly.GeometryElement
end

const nil = Nil

render(geom::Nil, theme::Gadfly.Theme, aes::Gadfly.Aesthetics,
        data::Gadfly.Data, scales::Dict{Symbol, ScaleElement},
        subplot_layer_aess::Vector{Gadfly.Aesthetics}) = nothing

# Subplot geometries require some more arguments to render. A simpler render
# function is defined and passed through to here for non-subplot geometries.
render(geom::Gadfly.GeometryElement, theme::Gadfly.Theme, aes::Gadfly.Aesthetics,
                subplot_layer_aess::Union{(Nothing), Vector{Gadfly.Aesthetics}},
                subplot_layer_datas::Union{(Nothing), Vector{Gadfly.Data}},
                scales::Dict{Symbol, ScaleElement}) =
        render(geom, theme, aes)

# Catchall
default_statistic(::Gadfly.GeometryElement) = Gadfly.Stat.identity()

include("geom/bar.jl")
include("geom/blank.jl")
include("geom/boxplot.jl")
include("geom/errorbar.jl")
include("geom/hexbin.jl")
include("geom/hvabline.jl")
include("geom/hvband.jl")
include("geom/label.jl")
include("geom/line.jl")
include("geom/point.jl")
include("geom/rectbin.jl")
include("geom/subplot.jl")
include("geom/ribbon.jl")
include("geom/violin.jl")
include("geom/polygon.jl")
include("geom/beeswarm.jl")
include("geom/segment.jl")

end # module Geom


text:
import ACE.OrthPolys: TransformedPolys

@doc raw"""
`struct Rn1pBasis <: OneParticleBasis`

One-particle basis of the form $R_n(r_{ij})$, i.e.,
no dependence on species or on $l$.

This type basically just translates the `TransformedPolys` into a valid
one-particle basis.
"""
mutable struct Rn1pBasis{T, TT, TJ} <: OneParticleBasis{T}
   R::TransformedPolys{T, TT, TJ}
end

# TODO: allow for the possibility that the symbol where the
#       position is stored is not `rr` but something else!
#
# TODO: Should we drop this type altogether and just
#       rewrite TransformedPolys to become a 1pbasis?

# ---------------------- Implementation of Rn1pBasis

Base.length(basis::Rn1pBasis) = length(basis.R)

# -> TODO : figure out how to do this well!!!
# Base.rand(basis::Ylm1pBasis) =
#       AtomState(rand(basis.zlist.list), ACE.Random.rand_vec(basis.J))

get_spec(basis::Rn1pBasis) =
      [ (n = n) for n = 1:length(basis) ]

==(P1::Rn1pBasis, P2::Rn1pBasis) =  ACE._allfieldsequal(P1, P2)

write_dict(basis::Rn1pBasis{T}) where {T} = Dict(
      "__id__" => "ACE_Rn1pBasis",
          "R" => write_dict(basis.R) )

read_dict(::Val{:ACE_Rn1pBasis}, D::Dict) = Rn1pBasis(read_dict(D["R"]))

fltype(basis::Rn1pBasis{T}) where T = T

gradtype(::Rn1pBasis{T}) where {T} = SVector{3, T}

symbols(::Rn1pBasis) = [:n]

indexrange(basis::Rn1pBasis) = Dict( :n => 1:length(basis) )

isadmissible(b, basis::Rn1pBasis) = (1 <= b.n <= length(basis))

degree(b, basis::Rn1pBasis) = b.n - 1

get_index(basis::Rn1pBasis, b) = b.n

rand_radial(basis::Rn1pBasis) = rand_radial(basis.R)

# ---------------------------  Evaluation code
#

alloc_B(basis::Rn1pBasis, args...) = alloc_B(basis.R)

alloc_dB(basis::Rn1pBasis, args...) =
      zeros( SVector{3, fltype(basis.R)}, length(basis) )

alloc_temp(basis::Rn1pBasis, args...) = alloc_temp(basis.R)

alloc_temp_d(basis::Rn1pBasis, args...) =
      (
      # alloc_temp_d(basis.R)...,
      dRdr = zeros(fltype(basis.R), length(basis.R)),
      )

evaluate!(B, tmp, basis::Rn1pBasis, X::AbstractState) =
      evaluate!(B, tmp, basis.R, norm(X.rr))

function evaluate_d!(dB, tmpd, basis::Rn1pBasis, X::AbstractState)
   r = norm(X.rr)
   r̂ = X.rr / r
   evaluate_d!(tmpd.dRdr, tmpd, basis.R, r)
   dB[:] .= Ref(r̂) .* tmpd.dRdr
   return dB
end

function evaluate_ed!(B, dB, tmpd, basis::Rn1pBasis, X::AbstractState)
   r = norm(X.rr)
   r̂ = X.rr / r
   evaluate!(B, tmpd, basis.R, norm(X.rr))
   evaluate_d!(tmpd.dRdr, tmpd, basis.R, r)
   dB[:] .= Ref(r̂) .* tmpd.dRdr
   return nothing
end


text:
# ========================================================================
#          EXP VARIABLE CELL IMPLEMENTATION
# ========================================================================
#
# F = exp(U) F0
# x =  F * F0^{-1} z  = exp(U) z
#
# ϕ( exp(U+tV) (z+tv) ) ~ ϕ'(x) ⋅ (exp(U) V) + ϕ'(x) ⋅ ( L(U, V) exp(-U) exp(U) z )
#    >>> ∂E(U) : V  =  [S exp(-U)'] : L(U,V)
#                   =  L'(U, S exp(-U)') : V
#                   =  L(U', S exp(-U)') : V
#                   =  L(U, S exp(-U)) : V     (provided U = U')

# define expm for StaticArrays since it is not provided (anymre?)

import StaticArrays

using LinearAlgebra: det

_expm(A::StaticArrays.SMatrix{N,N,T}) where {N,T} =
   StaticArrays.SMatrix{N,N,T}(exp(Array(A)))

mutable struct ExpVariableCell <: AbstractConstraint
   ifree::Vector{Int}
   X0::Vector{JVecF}
   F0::JMatF
   pressure::Float64
   fixvolume::Bool
end

function ExpVariableCell(at::AbstractAtoms;
               free=nothing, clamp=nothing, mask=nothing,
               pressure = 0.0, fixvolume=false)
   if pressure != 0.0 && fixvolume
      warning("the pressure setting will be ignored when `fixvolume==true`")
   end
   return ExpVariableCell( analyze_mask(at, free, clamp, mask),
                        positions(at), defm(at), pressure, fixvolume )
end

function logm_defm(at::AbstractAtoms, cons::ExpVariableCell)
   F = defm(at)
   # remove the reference deformation (F = expm(U) * F0)
   expU = F * inv(cons.F0)
   # expU should be spd, but roundoff (or other things) might mess this up
   # do we need an extra check here?
   U = log(expU |> Array) |> JMat
   U = real(0.5 * (U + U'))
   # check that expm(U) * F0 ≈ F (if not, then something has gone horribly wrong)
   if norm(F - _expm(U) * cons.F0, Inf) > 1e-12
      @show F
      @show _expm(U) * cons.F0
      error("something has gone wrong; U is not symmetric?")
   end
   return U, F
end

expposdofs(x) = x[1:end-6]
expcelldofs(x) = x[end-5:end]
U2dofs(U) = U[(1,2,3,5,6,9)]
dofs2U(x) = JMat(expcelldofs(x)[(1,2,3,2,4,5,3,5,6)])

function position_dofs(at::AbstractAtoms, cons::ExpVariableCell)
   X = positions(at)
   U, _ = logm_defm(at, cons)
   # tranform the positions back to the reference cell (F0)
   broadcast!(x -> _expm(-U) * x, X, X)
   # construct dof vector
   return [mat(X)[cons.ifree]; U2dofs(U)]
end

function set_position_dofs!(at::AbstractAtoms, cons::ExpVariableCell, x::Dofs)
   U = dofs2U(x)
   expU = _expm(U)
   F = expU * cons.F0
   Z = copy(cons.X0)
   mat(Z)[cons.ifree] = expposdofs(x)
   broadcast!(z -> expU * z, Z, Z)
   set_positions!(at, Z)
   set_defm!(at, F)
   return at
end

"""
Directional derivative of matrix exponential. See Theorem 2.1 in
AL-MOHY, HIGHAM SIAM J. Matrix Anal. Appl. 30, 2009.
"""
function dexpm_3x3(X::AbstractMatrix, E::AbstractMatrix)
   @assert size(X) == size(E) == (3, 3)
   return _expm([ X E; zeros(3,3) X ])[1:3, 4:6]
end

function gradient(at::AbstractAtoms, cons::ExpVariableCell)
   U, F = logm_defm(at, cons)
   G = forces(at)
   broadcast!(g -> - (_expm(U) * g), G, G)   # G[n] <- -exp(U) * G[n]
   T = dexpm_3x3(U, - virial(at) * _expm(-U))
   # add the forces acting on the upper diagonal of U to the lower diagonal
   # this way we ensure that T : V = U2dofs(T) : U2dofs(V)
   T[[2,3,6]] += T[[4,7,8]]
   # add the pressure component
   T[[1,5,9]] += cons.pressure * exp(trace(U)) * U[[1,5,9]]
   return [ mat(G)[cons.ifree]; U2dofs(T) ]
end

energy(at::AbstractAtoms, cons::ExpVariableCell) =
               energy(at) + cons.pressure * det(defm(at)) / det(cons.F0)
               #                            >>>>>>  exp(trace(U))  <<<<< (tested)

# TODO: fix this once we implement the volume constraint ??????
project!(at::AbstractAtoms, cons::ExpVariableCell) = at


text:
using Test, InteractiveUtils

using JuMP
using ECOS, SCS
using LightGraphs, MetaGraphs

using Seleroute

struct TestEdgeWise <: EdgeWiseObjectiveFunction end

@testset "Seleroute.jl" begin
    @testset "Data" begin
        include("data_modeltype.jl")
        include("data_parameters.jl")
        include("data_solution.jl")
        include("data_topology.jl")
    end
    @testset "Base models" begin
        include("basemodels_knowncapacities.jl")
        include("basemodels_objectives.jl")
    end
    @testset "Compute" begin
        include("compute_certain.jl")
        include("compute_oblivious_demand.jl")
    end
end


text:
__precompile__(true)
module NNFEM
using PyPlot
using JLD2
using MAT
using PyCall
using Statistics
using LinearAlgebra
using ADCME
using SparseArrays
using ForwardDiff
using Random
using Libdl

STORAGE = Dict{String, Any}()

include("utils/gmsh.jl")
include("utils/options.jl")
include("utils/shapeFunctions.jl")
include("utils/matrix.jl")
include("materials/PlaneStress.jl")
include("materials/PlaneStrain.jl")
include("materials/PlaneStressPlasticity.jl")
include("materials/Elasticity1D.jl")
include("materials/Plasticity1D.jl")
include("materials/Viscoplasticity1D.jl")
include("materials/NeuralNetwork2D.jl")
include("materials/NeuralNetwork1D.jl")
include("materials/PathDependent1D.jl")
include("materials/Scalar1D.jl")
include("materials/ViscoelasticityMaxwell.jl")
include("materials/PlaneStressPlasticityLawBased.jl")
include("materials/ViscoelasticityProny.jl")
include("materials/PlaneStressIncompressibleRivlinSaunders.jl")
include("elements/FiniteStrainContinuum.jl")
include("elements/SmallStrainContinuum.jl")
include("elements/FiniteStrainTruss.jl")
include("fem/fem.jl")
include("fem/fem1D.jl")
include("fem/assembly.jl")
include("utils/Visualize.jl")
include("utils/Visualize2.jl")
include("utils/io.jl")
include("utils/Testsuit.jl")
include("utils/nnconstitutivelaw.jl")
include("solvers/NNSolver.jl")
include("solvers/Solvers.jl")
include("solvers/SolversV2.jl")
include("solvers/AdjointSolver.jl")
include("adutils/helpers.jl")
include("adutils/kernels.jl")
include("adutils/solvers.jl")
include("adutils/solvers_standalone.jl")
include("adutils/representations.jl")
include("adutils/materials.jl")

end


text:
export Resolvent, qresolvent

abstract type MatrixTransformations end

"""
```julia
Resolvent(A::Matrix)
```

- return the resolvent **function** of `A`, ``R(z)= (A-zI)^{-1}``

# Examples

Get the resolvent **function** of a normalized GOE(500).

And evaluate the resolvent at `0`
```julia
R = Resolvent(rand(GOE(500),norm = true)) # this a function

R(0)

500×500 Matrix{Float64}:
  0.460903   -0.384016   0.268931  …  -0.59029    -0.178288
 -0.384016    0.489355  -0.366855      0.452094   -0.190179
  0.268931   -0.366855  -0.697169     -0.221433    0.535448
  0.388846    0.946668  -0.354459     -0.796717    0.207156
  ⋮                                ⋱
 -0.532533    0.950826  -0.366618      0.0926931  -0.288692
 -0.0874979  -0.553262  -0.388048     -0.0156897   0.177695
 -0.59029     0.452094  -0.221433      0.9136      0.343194
 -0.178288   -0.190179   0.535448      0.343194   -0.141701
```
"""
struct Resolvent<:MatrixTransformations
    M::Matrix
end

function (Resolvent::MatrixTransformations)(z)
    (Resolvent.M-UniformScaling(z))^(-1)
end

function Base.show(io::IO, R::Resolvent)
    print("Resolvent{Matrix}")
end

"""
```julia
qresolvent(A)
```
returns the quaternion resolvent **function** of A.
"""
function qresolvent(A::AbstractMatrix)
    return function f(z::Number,η::Number)
    ([UniformScaling(-η) A-UniformScaling(z)
    A'-UniformScaling(z') UniformScaling(-η)])^(-1)
end
end


text:
using Distributed
addprocs(length(Sys.cpu_info())-1) # use all cores
@everywhere using CarrierCapture
@everywhere using NPZ
@everywhere using SharedArrays

# define the necessary parameters for the PES'
@everywhere ħω_i = 8e-3 # 8 meV
@everywhere ħω_f = 8e-3 # 8 meV
@everywhere ΔQs = collect(range(0, stop=25, length=25)) # amu^0.5 Å
@everywhere ΔEs = collect(range(0, stop=2.5, length=10)) # eV

# initialise array
@everywhere ccArray = SharedArray{Float64}(length(ΔQs), length(ΔEs))

# Parallelise over the biggest loop
@sync @distributed for (qi,ΔQ) in collect(enumerate(ΔQs))
    for (ei,ΔE) in enumerate(ΔEs)
        poti, potf = fitHarmonicParams(ħω_i, ħω_f, ΔQ, ΔE)
        capt_coeff, barrier_height = getHarmonicCapture(poti, potf)

        for (i,cc) in enumerate(capt_coeff) # temperature dependence is stored in an array
            ccArray[qi, ei, i] = cc
        end
    end
end

# store output in numpy array for analysis
npzwrite("ccArray.npz", ccArray)