text:
# ------------------------------------------------------------------------------------------
# # Cadenas (Strings)
#
# Temas:
# 1. Cómo conseguir una cadena
# 2. Interpolación de cadenas
# 3. Concatenación de cadenas
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Cómo conseguir una cadena
#
# Enclose your characters in " " or """ """!
# ------------------------------------------------------------------------------------------

s1 = "Yo soy una cadena."

s2 = """Yo también soy una cadena. """

# ------------------------------------------------------------------------------------------
# Existen algunas diferencias funcionales entre expresas cadenas con una sola y triple
# comillas.

# Una diferencia es que en la segunda opción puedes usar comillas dentro de tu cadena.
# ------------------------------------------------------------------------------------------

"Aquí me sale un "error" porque es ambiguo dónde se acaba la cadena. "

"""Mira Mamá, sin "errors"!!! """

# ------------------------------------------------------------------------------------------
# Ojo, '' define a un caracter, ¡ NO una cadena!
# ------------------------------------------------------------------------------------------

typeof('a')

'Esto va a dar un error'

# ------------------------------------------------------------------------------------------
# NOTA: En Julia todas las cadenas por default son del tipo UTF-8. Esto significa que
# podemos usar nuestros querídisimos acentos del habla hispana y nuestros signos de apertura
# de interrogación sin miedo a que nuestro código no corra en una máquina ajena. ¡Inténtalo!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Interpolación de cadenas
#
# Podemos usar el símbolo $ para insertar variables existentes a una cadena y para evaluar
# expresiones en una cadena.

# Abajo hay un ejemplo con información súmamente sensible.
# ------------------------------------------------------------------------------------------

nombre = "Jane"
num_dedos = 10
num_dedos_de_los_pies = 10

println("Hola, me llamo $nombre.")
println("Yo tengo $num_dedos dedos y $num_dedos_de_los_pies. ¡Esos son $(num_dedos + num_dedos_de_los_pies) dígitos en total!!")

# ------------------------------------------------------------------------------------------
# ## Concatenación de cadenas
#
# ¡Abajo hay tres maneras de concatenar cadenas!

# La primera es usando la función de `string()`

# `string()` convierte no-cadenas en cadenas.
# ------------------------------------------------------------------------------------------

string("¿Cuántos gatos ", "son demasiados?")

string("No lo sé, pero ", 10, " son muy pocos.")

# ------------------------------------------------------------------------------------------
# We can also use `*` or string interpolation!
# ------------------------------------------------------------------------------------------

s3 = "¿Cuántos gatos ";
s4 = "son demasiados?";

s3*s4

"$s3$s4"

# ------------------------------------------------------------------------------------------
# ### Ejercicios
#
# 2.1 Crea una cadena que diga "hola" 1000 veces.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 2.2 Agrega dos números dentro de una cadena.
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using CALCEPH
using Test

CALCEPH.disableCustomHandler()
CALCEPH.setCustomHandler(s::String->Nothing )
testpath = joinpath(dirname(pathof(CALCEPH)), "..", "test")

# NAIF ID tests

for (name,id) ∈ naifId.id
    @test name ∈ naifId.names[id]
end

for (id,names) ∈ naifId.names
    for name ∈ names
       @test naifId.id[name] == id
   end
end

@test naifId.id[:ssb] == naifId.id[:solar_system_barycenter] == 0
@test naifId.id[:sun] == 10
@test naifId.id[:mercury_barycenter] == 1
@test naifId.id[:mercury] == 199
@test naifId.id[:venus_barycenter] == 2
@test naifId.id[:venus] == 299
@test naifId.id[:emb] == naifId.id[:earth_barycenter] == 3
@test naifId.id[:moon] == 301
@test naifId.id[:earth] == 399
@test naifId.id[:mars_barycenter] == 4
@test naifId.id[:phobos] == 401
@test naifId.id[:deimos] == 402
@test naifId.id[:mars] == 499
@test naifId.id[:jupiter_barycenter] == 5
@test naifId.id[:io] == 501
@test naifId.id[:europa] == 502
@test naifId.id[:ganymede] == 503
@test naifId.id[:callisto] == 504
@test naifId.id[:jupiter] == 599
@test naifId.id[:saturn_barycenter] == 6
@test naifId.id[:titan] == 606
@test naifId.id[:saturn] == 699
@test naifId.id[:uranus_barycenter] == 7
@test naifId.id[:uranus] == 799
@test naifId.id[:neptune_barycenter] == 8
@test naifId.id[:triton] == 801
@test naifId.id[:neptune] == 899
@test naifId.id[:pluto_barycenter] == 9
@test naifId.id[:charon] == 901
@test naifId.id[:pluto] == 999

# test error case: changing name->id mapping
@test_throws CALCEPHException CALCEPH.add!(naifId,:jupiter,1)
# test error case: parsing wrongly formatted body id input file
bid = CALCEPH.BodyId()
@test_throws CALCEPHException CALCEPH.loadData!(bid,joinpath(testpath,"badIds.txt"))

# check memory management
eph1 = Ephem(joinpath(testpath,"example1.dat"))
eph2 = Ephem([joinpath(testpath,"example1.bsp"),
                     joinpath(testpath,"example1.tpc")])

@test eph1.data != C_NULL
finalize(eph1)
@test eph1.data == C_NULL

@test eph2.data != C_NULL
finalize(eph2)
@test eph2.data == C_NULL
finalize(eph2)
CALCEPH._ephemDestructor(eph2)

# Opening invalid ephemeris
@test_throws CALCEPHException Ephem(String[])

# check constants
eph1 = Ephem(joinpath(testpath,"example1.dat"))
eph2 = Ephem([joinpath(testpath,"example1.bsp"),
                     joinpath(testpath,"example1.tpc")])
eph3 = Ephem([joinpath(testpath,"checktpc_11627.tpc")])
eph4 = Ephem([joinpath(testpath,"checktpc_str.tpc")])
con1 = constants(eph1)
con2 = constants(eph2)
con3 = constants(eph3)
con4 = constants(eph4)

@test isa(con1,Dict{Symbol,Any})
@test length(con1) == 402
@test con1[:EMRAT] ≈ 81.30056
@test isa(con2,Dict{Symbol,Any})
@test length(con2) == 313
@test con2[:AU] ≈ 1.49597870696268e8
@test isa(con3,Dict{Symbol,Any})
@test length(con3) == 3
@test con3[:BODY000_GMLIST4] == [ 199.0 ; 299.0 ; 301.0 ; 399.0 ]
@test con3[:BODY000_GMLIST2] == [ 499 ; 599 ]
@test con3[:BODY000_GMLIST1] == 699
@test isa(con4,Dict{Symbol,Any})
@test length(con4) == 4
@test con4[:MESSAGE] == "You can't always get what you want."
@test con4[:DISTANCE_UNITS] == "KILOMETERS"
@test con4[:MISSION_UNITS] == [ "KILOMETERS" ; "SECONDS" ; "KILOMETERS/SECOND" ]
@test con4[:CONTINUED_STRINGS] == ["This //", "is //", "just //", "one long //", "string.", "Here's a second //", "continued //", "string."]

# Retrieving constants from closed ephemeris
finalize(eph2)
@test_throws CALCEPHException constants(eph2)

# test compute*
# test data and thresholds from CALCEPH C library tests
inpop_files = [joinpath(testpath,"example1.dat")]
spk_files = [joinpath(testpath,"example1.bsp"),
             joinpath(testpath,"example1.tpc"),
             joinpath(testpath,"example1.tf"),
             joinpath(testpath,"example1.bpc"),
             joinpath(testpath,"example1spk_time.bsp")]

testfile = joinpath(testpath,"example1_tests.dat")

test_data = [
    (inpop_files,false),
    (spk_files,false),
    (inpop_files,true),
    (spk_files,true)
]

include("testfunction1.jl")

for (ephFiles,prefetch) in test_data
    testFunction1(testfile,ephFiles,prefetch)
end

testfile2 = joinpath(testpath,"example1_tests_naifid.dat")

include("testfunction2.jl")

for (ephFiles,prefetch) in test_data
    testFunction2(testfile,testfile2,ephFiles,prefetch)
end

# test error case wrong order
eph1 = Ephem(joinpath(testpath,"example1.bsp"))
@test_throws CALCEPHException compute(eph1,0.0,0.0,1,0,0,4)
@test_throws CALCEPHException compute(eph1,0.0,0.0,1,0,0,-1)

# test error case:
@test_throws CALCEPHException compute(eph1,0.0,0.0,-144,0,0)

# Five-Point Stencil

f(x)=x^8
@test_throws ErrorException CALCEPH.fivePointStencil(f,1.5,5,0.001)
@test_throws ErrorException CALCEPH.fivePointStencil(f,1.5,-1,0.001)
@test_throws ErrorException CALCEPH.fivePointStencil(f,1.5,4,0.0)
val = CALCEPH.fivePointStencil(f,1.5,4,0.001)
ref = [25.62890625,136.6875,637.875,2551.5,8505.0]
@test ref[1] ≈ val[1] atol=1e-10
@test ref[2] ≈ val[2] atol=1e-8
@test ref[3] ≈ val[3] atol=1e-5
@test ref[4] ≈ val[4] atol=1e-2
@test ref[5] ≈ val[5] atol=1e-2

# introspection
eph = Ephem(inpop_files)
@test timeScale(eph) == 1

records = positionRecords(eph)
@test length(records) == 12

records = orientationRecords(eph)
@test length(records) == 1

@test timespan(eph) == (2.442457e6, 2.451545e6, 1)

# rotangmom
eph = Ephem(joinpath(testpath,"example2_rotangmom.dat"))
a = rotAngMom(eph,2.4515e6,0.0,399,useNaifId+unitSec)
b = rotAngMom(eph,2.4515e6,0.0,399,useNaifId+unitSec,1)
@test a == b
@test length(a) == 6

ext: jl


text:
export SymmetricGroup,
       infinity_norm,
       infinity_distance,
       fubini_study,
       PrecisionOption,
       PRECISION_ADAPTIVE,
       PRECISION_FIXED_64,
       PRECISION_FIXED_128

"""
    unpack(a::Union{Nothing, T}, b::T)

Returns `a` if it is not `nothing`, otherwise `b`.
"""
unpack(a::Union{Nothing,T}, b::T) where {T} = a === nothing ? b : a

"""
    SymmetricGroup(n)

Group action of the symmetric group S(n).
"""
struct SymmetricGroup{N}
    permutations::Vector{SVector{N,Int}}
end
SymmetricGroup(N::Int) = SymmetricGroup(permutations(Val(N)))
function permutations(::Val{N}) where {N}
    s = StaticArrays.MVector{N}(1:N)
    perms = [SVector(s)]
    while true
        i = N - 1
        while i >= 1 && s[i] >= s[i+1]
            i -= 1
        end
        if i > 0
            j = N
            while j > i && s[i] >= s[j]
                j -= 1
            end
            s[i], s[j] = s[j], s[i]
            reverse!(s, i + 1)
        else
            s[1] = N + 1
        end

        s[1] > N && break
        push!(perms, SVector(s))
    end
    perms
end

Base.eltype(::Type{SymmetricGroup{N}}) where {N} = SVector{N,Int}
Base.length(p::SymmetricGroup{N}) where {N} = length(p.permutations)
Base.iterate(p::SymmetricGroup) = iterate(p.permutations)
Base.iterate(p::SymmetricGroup, s) = iterate(p.permutations, s)

"""
    deprecatekwarg(oldkw, newkw)

Logs a deprecation warning and assigns `oldkw` to `newkw` if `oldkw` is not `nothing`.
"""
macro deprecatekwarg(oldkw, newkw)
    quote
        if $(esc(oldkw)) !== nothing
            old = $(Expr(:quote, oldkw))
            new = $(Expr(:quote, newkw))
            @warn(
                "`$(old)=$($(esc(oldkw)))` is deprecated, use `$(new)=$($(esc(oldkw)))` instead.",
            )
            $(esc(newkw)) = $(esc(oldkw))
        end
    end
end

"""
    invalid_keywords(kwargs, allowed)

Look for invalid keywords in the `kwargs`.
"""
function invalid_kwargs(kwargs, allowed)
    invalids = []
    for kwarg in kwargs
        kw = first(kwarg)
        if !any(isequal(kw), allowed)
            push!(invalids, kwarg)
        end
    end
    invalids
end

"""
    nthroot(x, n)

Compute the `n`-th root of `x`.
"""
function nthroot(x::Real, N::Integer)
    if N == 4
        √(√(x))
    elseif N == 2
        √(x)
    elseif N == 3
        cbrt(x)
    elseif N == 1
        x
    elseif N == 0
        one(x)
    else
        x^(1 / N)
    end
end

"""
    log₁₀(x)

Same as `log10` but faster.
"""
log₁₀(x) = log(x) / 2.302585092994046

"""
    @modulenum(name, block)

This is a modification of `@enum` and creates an intermediate enum.

## Example

The definition

```julia
@moduleenum Car begin
    audi
    volkswagen
    bmw
end
```

expands into

```julia
module Car
    @enum t begin
        audi
        volkswagen
        bmw
    end
end
import .Car
"""
macro moduleenum(name, content...)
    mname = esc(name)
    quote
        @eval module $name
        @enum t $(content...)
        end
        import .$name
    end
end

"""
     print_fieldnames(io::IO, obj)

 A better default printing for structs.
 """
function print_fieldnames(io::IO, obj)
    println(io, typeof(obj), ":")
    for name in fieldnames(typeof(obj))
        if getfield(obj, name) !== nothing
            val = getfield(obj, name)
            print(io, " • ", name, " → ")
            if val isa AbstractFloat
                println(io, round(val; sigdigits = 5))
            else
                println(io, val)
            end
        end
    end
end

"""
    check_kwargs_empty(kwargs, [allowed_kwargs])

Chack that the list of `kwargs` is empty. If not, print all unsupported keywords
with their arguments.
"""
function check_kwargs_empty(kwargs, allowed_kwargs = [])
    if !isempty(kwargs)
        msg = "Unexpected keyword argument(s): "
        first_el = true
        for kwarg in kwargs
            if !first_el
                msg *= ", "
            end
            msg *= "$(first(kwarg))=$(last(kwarg))"
            first_el = false
        end
        if !isempty(allowed_kwargs)
            msg *= "\nAllowed keywords are\n"
            msg *= join(allowed_kwargs, ", ")
        end
        throw(ArgumentError(msg))
    end
end

"""
    is_real_vector(v::AbstractVector, tol=1e-6)

Check whether the 2-norm of the imaginary part of `v` is at most `tol`.
"""
is_real_vector(z::AbstractVector{<:Real}, tol = 1e-6) = true
function is_real_vector(z::AbstractVector{<:Complex}, tol = 1e-6)
    total = zero(real(eltype(z)))
    for zᵢ in z
        total += abs2(imag(zᵢ))
    end
    sqrt(total) < tol
end
is_real_vector(z::NTuple{N,T}, tol = 1e-6) where {N,T} = is_real_vector(SVector{N}(z), tol)
is_real_vector(v::PVector, tol = 1e-6) = isreal(v, tol)

"""
    real_vector(v)

Obtain the real part of a vector.
"""
function real_vector(v::PVector{T}) where {T}
    w = LA.normalize(v)
    λ = map(dimension_indices(w)) do rᵢ
        maxᵢ = zero(real(T))
        max_ind = 0
        for i in rᵢ
            wᵢ = abs(w[i])
            if wᵢ > maxᵢ
                max_ind = i
                maxᵢ = wᵢ
            end
        end
        cis(-angle(w[max_ind]))
    end
    LA.rmul!(w, λ)
    real.(w)
end
real_vector(v::AbstractVector) = real.(v)

"""
    randseed(range=1_000:1_000_000)

Return a random seed in the range `range`.
"""
randseed(range = 1_000:1_000_000) = rand(range)

"""

    infinity_distance(u, v)

Compute the ∞-norm of `u-v`.
"""
function infinity_distance(z₁::AbstractVector, z₂::AbstractVector)
    m = abs(z₁[1] - z₂[1])
    n₁, n₂ = length(z₁), length(z₂)
    if n₁ ≠ n₂
        return convert(typeof(m), Inf)
    end
    @inbounds for k = 2:n₁
        @fastmath m = max(m, abs(z₁[k] - z₂[k]))
    end
    m
end

"""

    infinity_norm(z)

Compute the ∞-norm of `z`. If `z` is a complex vector this is more efficient
than `norm(z, Inf)`.
"""
infinity_norm(z::AbstractVector{<:Complex}) = sqrt(maximum(abs2, z))

function randomish_gamma()
    # Usually values near 1, i, -i, -1 are not good randomization
    # Therefore we artificially constrain the choices
    theta = rand() * 0.30 + 0.075 + (rand(Bool) ? 0.0 : 0.5)
    cis(2π * theta)
end

"""
    splitkwargs(kwargs, supported_keywords)

Split the vector of `kwargs` in two vectors, the first contains all `kwargs`
whose keywords appear in `supported_keywords` and the rest the other one.
"""
function splitkwargs(kwargs, supported_keywords)
    supported = []
    rest = []
    for kwarg in kwargs
        if any(kw -> kw == first(kwarg), supported_keywords)
            push!(supported, kwarg)
        else
            push!(rest, kwarg)
        end
    end
    supported, rest
end

start_solution_sample(xs) = promote_start_solution(first(xs))
start_solution_sample(x::AbstractVector{<:Number}) = promote_start_solution(x)
promote_start_solution(x::AbstractVector{ComplexF64}) = x
promote_start_solution(x::SVector) = map(xᵢ -> first(promote(xᵢ, 0.0im)), x)
function promote_start_solution(x)
    x_new = similar(x, promote_type(eltype(x), ComplexF64), length(x))
    copyto!(x_new, x)
    x_new
end

"""
    ComplexSegment(start, target)

Represents the line segment from `start` to `finish`.
Supports indexing of the values `t ∈ [0, length(target-start)]` in order to get the
corresponding point on the line segment.
"""
struct ComplexSegment
    start::ComplexF64
    target::ComplexF64
    # derived
    Δ_target_start::ComplexF64
    abs_target_start::Float64
end
function ComplexSegment(start::Number, target::Number)
    Δ_target_start = convert(ComplexF64, target) - convert(ComplexF64, start)
    abs_target_start = abs(Δ_target_start)

    ComplexSegment(start, target, Δ_target_start, abs_target_start)
end

function Base.getindex(segment::ComplexSegment, t::Real)
    Δ = t / segment.abs_target_start
    segment.start + Δ * segment.Δ_target_start
end
Base.length(segment::ComplexSegment) = segment.abs_target_start

function Base.show(io::IO, segment::ComplexSegment)
    print(io, "ComplexSegment($(segment.start), $(segment.target))")
end
Base.show(io::IO, ::MIME"application/prs.juno.inline", opts::ComplexSegment) = opts

# Parallelization

"""
    single_thread_blas()

We set the number of BLAS threads to 1 since we multithread by ourself.
But even if we would not, the overhead of the threading is so large
that the single threaded version is around two times faster (at least on Mac Julia v0.6.2).
Returns the previous number of BLAS threads.
"""
function single_thread_blas()
    nblas_threads = get_num_BLAS_threads()
    set_num_BLAS_threads(1)
    nblas_threads
end

set_num_BLAS_threads(n) = LinearAlgebra.BLAS.set_num_threads(n)
get_num_BLAS_threads() = convert(Int, _get_num_BLAS_threads())
# This is into 0.7 but we need it for 0.6 as well
const _get_num_BLAS_threads = function () # anonymous so it will be serialized when called
    blas = LinearAlgebra.BLAS.vendor()
    # Wrap in a try to catch unsupported blas versions
    try
        if blas == :openblas
            return ccall((:openblas_get_num_threads, Base.libblas_name), Cint, ())
        elseif blas == :openblas64
            return ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())
        elseif blas == :mkl
            return ccall((:MKL_Get_Max_Num_Threads, Base.libblas_name), Cint, ())
        end

        # OSX BLAS looks at an environment variable
        if Sys.isapple()
            return ENV["VECLIB_MAXIMUM_THREADS"]
        end
    catch
    end

    return nothing
end


text:
function calculateDLMP(dispatch, buses, generators, lines, SMP, gensetP, gensetU)
   lineset = 1:length(lines)
   busset = 1:length(buses)
   genset = 1:length(generators)

   B_g = []
   for g in genset
      push!(B_g, generators[g].location)
   end
   B_d = setdiff(busset,B_g)

   B_gn = Array{Array{Int64}}(length(buses),1)
   for ii in 1:length(buses)
      B_gn[ii] = Int64[]
   end
   for g in genset
      push!(B_gn[generators[g].location], g)
   end
   #########################
   m = Model(solver = IpoptSolver())

   @variable(m, v[b in busset] >= 0) # variable for voltage square, voltage^2
   @variable(m, a[l in lineset] >= 0) # variable for current square, current^2
   @variable(m, fp[l in lineset])
   @variable(m, fq[l in lineset])
   @variable(m, pg[g in genset])
   @variable(m, qg[g in genset])
   @variable(m, OC >= 0)
   @variable(m, OCgen[g in genset] >= 0)
   #@variable(m, pgextra[g in genset]) ##################

   @objective(m, Max, -OC)

   @constraint(m, OperatingCost, OC == sum(OCgen[g] for g in genset))
   #### Polynomial Generation Cost
   #@constraint(m, GenCost[g in genset], OCgen[g] == sum(generators[g].cost[k]*pg[g]^(length(generators[g].cost)-k) for k in 1:length(generators[g].cost)))
   #### Linear Generation Cost
   #@constraint(m, GenCost[g in setdiff(genset,geninter)], OCgen[g] == generators[g].cost[2]*pg[g] + generators[g].cost[3])
   @constraint(m, GenCost[g in genset], OCgen[g] == generators[g].cost[2]*pg[g] + generators[g].cost[3])

   @NLconstraint(m, LineCapFW[l = 1:length(lineset)], (fp[l])^2 + (fq[l])^2 <= ( lines[l].u)^2 )
   #@constraint(m, LineCapFW[l = lineset], norm([fp[l], fq[l]]) <= lines[l].u )

   @NLconstraint(m, LineCapBW[l = 1:length(lineset)], (fp[l] - a[l]*lines[l].r)^2 + (fq[l] - a[l]*lines[l].x)^2 <= ( lines[l].u)^2 )
   #@constraint(m, LineCapBW[l = lineset], norm([fp[l] - a[l]*lines[l].r, (fq[l] - a[l]*lines[l].x)]) <= lines[l].u )

   @constraint(m, BetweenNodes[l = lineset], v[lines[l].fbus] - 2*(lines[l].r*fp[l] + lines[l].x*fq[l]) + a[l]*(lines[l].r^2 + lines[l].x^2) == v[lines[l].tbus])

   @constraint(m, SOCP[l = lineset], ((fp[l])^2 + (fq[l])^2) <= v[lines[l].fbus] * a[l] )
   #@constraint(m, SOCP[l = lineset], norm([2*fp[l], 2*fq[l], v[lines[l].fbus]-a[l]]) <= v[lines[l].fbus] + a[l])

   if testsystem == "AP15busDN" || testsystem == "IEEE33busDN" || testsystem == "ISONE8busTN" || testsystem == "IEEE33busDN2"
      @constraint(m, v[1] == 1) #Voltage constraint for root node
   end

   @constraint(m, PBalance[b = busset], sum(fp[l] for l in buses[b].outline)  - sum(fp[l] - lines[l].r*a[l] for l in buses[b].inline) - sum(pg[g] for g in B_gn[b]) + buses[b].Pd + v[b]*buses[b].Gs == 0)
   @constraint(m, QBalance[b = busset], sum(fq[l] for l in buses[b].outline)  - sum(fq[l] - lines[l].x*a[l] for l in buses[b].inline) - sum(qg[g] for g in B_gn[b])+ buses[b].Qd - v[b]*buses[b].Bs == 0)

   @constraint(m, Vmax[b = busset], v[b] <= buses[b].Vmax^2) #upper limit constraint for voltage square
   @constraint(m, Vmin[b = busset], v[b] >= buses[b].Vmin^2) #lower limit constraint for voltage square

   @constraint(m, PGmin[g = genset], pg[g] >= generators[g].Pmin)
   @constraint(m, PGmax[g = genset], pg[g] <= generators[g].Pmax)
   @constraint(m, QGmin[g = genset], qg[g] >= generators[g].Qmin)
   @constraint(m, QGmax[g = genset], qg[g] <= generators[g].Qmax)

   @constraint(m, genPower[g in gensetP], dispatch[g] == pg[g])
   #@constraint(m, genPower[g in gensetP], dispatch[g] + pgextra[g] == pg[g])
   #@constraint(m, extraPower[g in gensetP], pgextra[g]^2 <= (1e-3 * dispatch[g])^2)

   #tic()
   status = solve(m)
   if status != Symbol(:Optimal)
      return status, [], [], [], [], []
   end

   obj_value = getobjectivevalue(m)
   v = getvalue(v)[:]
   a = getvalue(a)[:]
   fp = getvalue(fp)[:]
   fq = getvalue(fq)[:]
   pg = getvalue(pg)[:]
   qg = getvalue(qg)[:]

   OC = getvalue(OC)
   OCgen = getvalue(OCgen)

   pgb = zeros(length(buses))
   qgb = zeros(length(buses))
   for b in 1:length(buses)
      if ~isempty(B_gn[b])
         pgb[b] = sum(pg[g] for g in B_gn[b])
         qgb[b] = sum(qg[g] for g in B_gn[b])
      end
   end

   dual_OC = getdual(OperatingCost)

   dual_LineCapFW = getdual(LineCapFW)
   dual_LineCapBW = getdual(LineCapBW)
   dual_BetweenNodes = getdual(BetweenNodes)
   # dual_SOCP = getdual(SOCP)

   dual_PBalance = getdual(PBalance)
   dual_QBalance = getdual(QBalance)

   dual_Vmax = getdual(Vmax)
   dual_Vmin = getdual(Vmin)

   dual_PGmax = getdual(PGmax)
   dual_PGmin = getdual(PGmin)
   dual_QGmax = getdual(QGmax)
   dual_QGmin = getdual(QGmin)
   A1 = Float64[]
   for l in lineset
       temp = ((fp[l]^2+fq[l]^2)*lines[l].x+a[l]*fq[l]*(lines[l].r^2-lines[l].x^2)-2*a[l]*fp[l]*lines[l].r*lines[l].x)/(
       (fp[l]^2+fq[l]^2)*lines[l].x-a[l]*fq[l]*(lines[l].r^2+lines[l].x^2))
       push!(A1, temp)
   end

   A2 = Float64[]
   for l in lineset
       temp = ((fp[l]^2+fq[l]^2)*lines[l].r-a[l]*fp[l]*(lines[l].r^2+lines[l].x^2))/(
       (fp[l]^2+fq[l]^2)*lines[l].x-a[l]*fq[l]*(lines[l].r^2+lines[l].x^2))
       push!(A2, temp)
   end

   A3 = Float64[]
   for l in lineset
       temp = (-(fp[l]^2+fq[l]^2)*lines[l].r+a[l]*fp[l]*(lines[l].r^2-lines[l].x^2)+2*a[l]*fq[l]*lines[l].r*lines[l].x)/(
       (fp[l]^2+fq[l]^2)*lines[l].x-a[l]*fq[l]*(lines[l].r^2+lines[l].x^2))
       push!(A3, temp)
   end

   A4 = Float64[]
   for l in lineset
       temp = (2*(fq[l]^3*lines[l].r-fp[l]^3*lines[l].x)+2*fp[l]*fq[l]*(fp[l]*lines[l].r-fq[l]*lines[l].x))/(
       (fp[l]^2+fq[l]^2)*lines[l].x-a[l]*fq[l]*(lines[l].r^2+lines[l].x^2))
       push!(A4, temp)
   end

   A5 = Float64[]
   for l in lineset
       temp = (2*(fq[l]^3*lines[l].r-2*fp[l]^3*lines[l].x)+2*fp[l]*fq[l]*(fp[l]*lines[l].r-fq[l]*lines[l].x))/(
       (fp[l]^2+fq[l]^2)*lines[l].x-a[l]*fq[l]*(lines[l].r^2+lines[l].x^2)) +
       (2*a[l]^2*(fq[l]*lines[l].r^3-fp[l]*lines[l].x^3)-4*a[l]*fp[l]*fq[l]*(lines[l].r^2-lines[l].x^2))/(
       (fp[l]^2+fq[l]^2)*lines[l].x-a[l]*fq[l]*(lines[l].r^2+lines[l].x^2)) +
       (4*a[l]*lines[l].r*lines[l].x*(fp[l]^2-fq[l]^2))/(
       (fp[l]^2+fq[l]^2)*lines[l].x-a[l]*fq[l]*(lines[l].r^2+lines[l].x^2)) +
       (-2*a[l]^2*lines[l].r*lines[l].x*(fp[l]*lines[l].r-fq[l]*lines[l].x))/(
       (fp[l]^2+fq[l]^2)*lines[l].x-a[l]*fq[l]*(lines[l].r^2+lines[l].x^2))
       push!(A5, temp)
   end
#=
   DLMP = Array{Float64,1}(length(buses))
   DLMP[1] = SMP
   Eq40 = zeros(length(buses))
   Eq41 = zeros(length(buses))
   Eq42 = zeros(length(buses))
   Eq43 = zeros(length(buses))
   Eq44 = zeros(length(buses))
   for b in 2:length(buses)
       DLMP[b] = A1[b-1]*DLMP[lines[b-1].tbus] + A2[b-1]*dual_QBalance[b] + A3[b-1]*dual_QBalance[lines[b-1].tbus] + A4[b-1]*dual_LineCapFW[b-1] + A5[b-1]*dual_LineCapBW[b-1]
       #DLMP[b] = A1[b-1]*dual_PBalance[lines[b-1].tbus] + A2[b-1]*dual_QBalance[b] + A3[b-1]*dual_QBalance[lines[b-1].tbus] + A4[b-1]*dual_LineCapFW[b-1] + A5[b-1]*dual_LineCapBW[b-1]
       Eq40[b] = A1[b-1]*DLMP[lines[b-1].tbus]
       #Eq40[b] = A1[b-1]*dual_PBalance[lines[b-1].tbus]
       Eq41[b] = A2[b-1]*dual_QBalance[b]
       Eq42[b] = A3[b-1]*dual_QBalance[lines[b-1].tbus]
       Eq43[b] = A4[b-1]*dual_LineCapFW[b-1]
       Eq44[b] = A5[b-1]*dual_LineCapBW[b-1]
   end
=#
   DLMP = Array{Float64,1}(length(buses))
   #DLMP[1] = SMP
   DLMP[1] = dual_PBalance[1]
   Eq40 = zeros(length(buses))
   Eq41 = zeros(length(buses))
   Eq42 = zeros(length(buses))
   Eq43 = zeros(length(buses))
   Eq44 = zeros(length(buses))
   for l in lineset
       DLMP[lines[l].fbus] = A1[l]*dual_PBalance[lines[l].tbus] + A2[l]*dual_QBalance[lines[l].fbus] + A3[l]*dual_QBalance[lines[l].tbus] + A4[l]*dual_LineCapFW[l] + A5[l]*dual_LineCapBW[l]
       Eq40[lines[l].fbus] = A1[l]*dual_PBalance[lines[l].tbus]
       Eq41[lines[l].fbus] = A2[l]*dual_QBalance[lines[l].fbus]
       Eq42[lines[l].fbus] = A3[l]*dual_QBalance[lines[l].tbus]
       Eq43[lines[l].fbus] = A4[l]*dual_LineCapFW[l]
       Eq44[lines[l].fbus] = A5[l]*dual_LineCapBW[l]
   end

   println()
   println("::Network Info::")
   NodeInfo = DataFrame(pg=round.(pgb,3),qg=round.(qgb,3),v=round.(v,3))
   LineInfo = DataFrame(flow=round.(sqrt.(fp.^2+fq.^2),3),i=round.(a,3),fp=round.(fp,3),fq=round.(fq,3))
   display(NodeInfo)
   println()
   println()
   display(LineInfo)

   #println()
   #println("::Network Info::")
   #NetworkInfo = DataFrame(pg=round.(pg,3),qg=round.(qg,3),flow=round.(fp.^2+fq.^2,3),v=round.(v,3),i=round.(a,3),fp=round.(fp,3),fq=round.(fq,3))
   #display(NetworkInfo)
   #println()
   #println()
   println("::DLMP Info::")
   DLMPInfo = DataFrame(λ_i=round.(DLMP,2), Eq40=round.(Eq40,2), Eq41=round.(Eq41,2), Eq42=round.(Eq42,2), Eq43=round.(Eq43,2), Eq44=round.(Eq44,2))
   display(DLMPInfo)
   println()

   return status, DLMP, pg, NodeInfo, LineInfo, DLMPInfo
end


text:
function solve()
    x = parse(Int, readline())
    x ≥ 0 ? x : 0
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using Test, HTTP, JSON

const dir = joinpath(dirname(pathof(HTTP)), "..", "test")
include("resources/TestRequest.jl")

@testset "HTTP" begin
    for f in ["ascii.jl",
              "chunking.jl",
              "utils.jl",
              "client.jl",
              "multipart.jl",
              "parsemultipart.jl",
              "sniff.jl",
              "uri.jl",
              "url.jl",
              "cookies.jl",
              "parser.jl",
              "loopback.jl",
              "websockets.jl",
              "messages.jl",
              "handlers.jl",
              "server.jl",
              "async.jl",
              "aws4.jl",
              "insert_layers.jl",
              "mwe.jl"]
        file = joinpath(dir, f)
        println("Running $file tests...")
        if isfile(file)
            include(file)
        else
            @show readdir(dirname(file))
        end
    end
end


text:
import ChainRulesCore.Tangent
import Yota: update!, setfield_nested!

mutable struct StructA
    t::Array{Float64, 1}
    s::Float64
end

mutable struct StructB
    a::StructA
    s::Float64
end

@testset "update" begin

    # setfield_nested
    b = StructB(StructA([1, 2, 3], 4), 5)
    setfield_nested!(b, (:a, :t), [-1.0, -2.0, -3.0])
    setfield_nested!(b, (:a, :s), -4.0)
    setfield_nested!(b, (:s,), -5.0)
    @test b.a.t == [-1.0, -2.0, -3.0]
    @test b.a.s == -4.0
    @test b.s == -5.0

    # struct update
    g = Dict((:a, :t) => [1.0, 2.0, 3.0],
             (:a, :s) => 4.0,
             (:s,) => 5.0)
    g = Tangent{StructB}(s=5.0, a=Tangent{StructA}(t=[1.0, 2.0, 3.0], s=4.0))
    update!(b, g)
    @test b.a.t == [-2.0, -4.0, -6.0]
    @test b.a.s == -8.0
    @test b.s == -10.0

    # struct update with function
    update!(b, g, (x, g) -> x .+ 2g)
    @test b.a.t == [0, 0, 0]
    @test b.a.s == 0
    @test b.s == 0

    # array update
    x = rand(3,4)
    update!(x, x)
    @test x == zero(x)

    # array update with function
    x = rand(3,4)
    xo = copy(x)
    update!(x, x, (x, g) -> x - 2g)
    @test x == -xo

    # with actual grad update
    b = StructB(StructA(ones(4), 1.0), 1.0)
    _, g = grad(b -> sum(b.a.t) + (b.a.s + b.s), b)
    update!(b, g[2], (x, gx) -> x .- 0.5gx)
    @test b.a.t == [0.5, 0.5, 0.5, 0.5]
    @test b.a.s == 0.5
    @test b.s == 0.5
end


text:
"""
    SamplingStats

A struct that tracks sampling information.

The fields available are:

- `start`: A `Float64` Unix timestamp indicating the start time of sampling.
- `stop`: A `Float64` Unix timestamp indicating the stop time of sampling.
- `duration`: The sampling time duration, defined as `stop - start`.
"""
struct SamplingStats
    start::Float64
    stop::Float64
    duration::Float64
end


text:
using LinearAlgebra
using Statistics

using InteractiveUtils

const MOD_METHODS = Dict{Module, Vector{Symbol}}()

methods = methodswith(AbstractArray)

for method in methods
    mod = method.module
    names = get!(MOD_METHODS, mod, Vector{Symbol}())
    push!(names, method.name)
end

for (k, v) in MOD_METHODS
    print(k)
    print(":\n\t")
    show(v)
    print("\n\n")
end


text:
module PedestrianAvoidancePOMDP

using Random
using Distributions
using AutomotiveDrivingModels
using AutoViz
using AutomotiveSensors
using AutomotivePOMDPs
using Parameters
using StaticArrays
using GridInterpolations
using POMDPs
using POMDPModelTools
using POMDPPolicies
using QMDP
using JLD2
using Reel
using EmergencyBrakingSystem

include("pomdp_types.jl")
include("spaces.jl")
include("transition.jl")
include("observation.jl")
include("belief.jl")
include("policy.jl")
include("pedestrian_pomdp.jl")
include("helpers.jl")
include("rendering.jl")
include("pedestrian_pomdp_emergency_system.jl")
include("evaluation_scenario.jl")

export
    PedestrianAvoidancePOMDPFrenet,
    PedestrianAvoidanceSystem,

    SingleOCFPOMDP,
    SingleOCFUpdater,
    SingleOCFBelief,
    SingleOCFAction,
    SingleOCFState,
    CrosswalkEnv,
    ObservationCallback,
    animate_record,
    generate_scenario,
    evaluateScenarioMetric,
    evaluateScenariosMetric

end


text:
######
# This example shows how to do function approximation in 2D using Padua points
# and in 3D by hand
######

using ApproxFun
    import ApproxFun: transform
f = (x,y) -> exp(x*y)

S = Chebyshev(0..1)^2

N = 10^2; p = points(S, N) # 105 padua points
f̃ = (xy) -> f(xy[1], xy[2])
F = f̃.(p)
F̌ = transform(S, F)  # 105 Chebyshev^2 coefficients ordered by polynomial degree
f̌ = Fun(S, F̌)
f̌(0.1,0.2) ≈ f(0.1,0.2)

S = Chebyshev(0..1)
p_x = points(S, 10)
p_y = points(S, 10)'
p_z = reshape(points(S, 10), 1, 1, 10)

f = (x,y,z) -> exp(x*y+z)

F = f.(p_x, p_y, p_z)

for k = 1:size(F,2),j=1:size(F,3)
    F[:,k,j] = transform(S,F[:,k,j])
end
for k = 1:size(F,1),j=1:size(F,3)
    F[k,:,j] = transform(S,F[k,:,j])
end
for k = 1:size(F,1),j=1:size(F,2)
    F[k,j,:] = transform(S,F[k,j,:])
end

f̌ = function(x,y,z)
    ret1 = Array{Float64}(undef, size(F,2), size(F,3))
    for j=1:size(F,2),l=1:size(F,3)
        ret1[j,l] = Fun(S,F[:,j,l])(x)
    end
    ret2 = Array{Float64}(undef, size(F,3))
    for l=1:size(F,3)
        ret2[l] = Fun(S,ret1[:,l])(y)
    end
    Fun(S, ret2)(z)
end

f̌(0.1,0.2,0.3) ≈ f(0.1,0.2,0.3)

max_stars_repo_path: Extras/Manual 3D interpolation.jl

max_issues_repo_path: Extras/Manual 3D interpolation.jl

max_forks_repo_path: Extras/Manual 3D interpolation.jl


text:
type TopographicPoint2D <: Singleton
  Heatmap::Array{Float64,2}
  radiusOfInterest::Float64
end

function evalPotential(topo::TopographicPoint2D, Xi::Array{Graphs.ExVertex,1}; N::Int64=100)
  warn("evalPotential TopographicPoint2D is not finished yet")
  return topo.Heatmap
end

function evalPotential(obs::PriorPose2, Xi::Array{Graphs.ExVertex,1}; N::Int64=200)
    cov = diag(obs.Cov)
    ret = zeros(3,N)
    for j in 1:N
      for i in 1:size(obs.Zi,1)
        ret[i,j] += obs.Zi[i,1] + (cov[i]*randn())
      end
    end
    return ret
end

# MOVED TO TRANSFORM UTILS
# function wrapRad(th::Float64)
#   if th >= pi
#     th -= 2.0*pi
#   end
#   if th < -pi
#     th += 2.0*pi
#   end
#   return th
# end
#
#
# R(th) = [[cos(th);-sin(th)]';[sin(th);cos(th)]'];
#
# function SE2(X::Array{Float64,1})
#     T = eye(3)
#     T[1:2,1:2] = R(X[3])
#     T[1,3] = X[1]
#     T[2,3] = X[2]
#     return T
# end
#
# function se2vee!(retval::Array{Float64,1}, T::Array{Float64,2})
#     retval[1] = T[1,3]
#     retval[2] = T[2,3]
#     retval[3] = wrapRad(atan2(-T[1,2], T[1,1]))
#     nothing
# end
#
# function se2vee(T::Array{Float64,2})
#     retval = zeros(3)
#     #retval[1:2,1] = T[1:2,3]
#     #retval[3,1] = wrapRad(atan2(-T[1,2], T[1,1]))
#     se2vee!(retval, T)
#     return retval
# end

# DX = [transx, transy, theta]
function addPose2Pose2(x::Array{Float64,1}, dx::Array{Float64,1})
    X = SE2(x)
    DX = SE2(dx)
    return se2vee(X*DX)
end

function addPose2Pose2!(retval::Array{Float64,1}, x::Array{Float64,1}, dx::Array{Float64,1})
  X = SE2(x)
  DX = SE2(dx)
  se2vee!(retval, X*DX)
  nothing
end

function evalPotential(odom::Pose2Pose2, Xi::Array{Graphs.ExVertex,1}, Xid::Int64)
    rz,cz = size(odom.Zij)
    Xval = Array{Float64,2}()
    # implicit equation portion -- bi-directional pairwise function
    if Xid == Xi[1].index #odom.
        #Z = (odom.Zij\eye(rz)) # this will be used for group operations
        Z = se2vee(SE2(vec(odom.Zij)) \ eye(3))
        Xval = getVal(Xi[2])
    elseif Xid == Xi[2].index
        Z = odom.Zij
        Xval = getVal(Xi[1])
    else
        error("Bad evalPairwise Pose2Pose2")
    end

    r,c = size(Xval)
    RES = zeros(r,c*cz)
    # cov = diag(odom.Cov)

    # increases the number of particles based on the number of modes in the measurement Z
    # dd = c*cz
    # ENT = randn(dd, 3) # this should be the covariate error from Distributions
    # for d in 1:dd
    #   ENT[:,d] = ENT[:,d].*cov[d]
    # end
    #di = 1.0/(c+1)
    # for i in 1:dd
    #     # ent = 1.0*[randn()*cov[1]; randn()*cov[2]; randn()*cov[3]] # TODO -- should be mvnormal, but more expensive
    #     idx = floor(Int,i*di+1) #floor(Int,i/(c+1)+1)
    #     RES[:,i] = addPose2Pose2(Xval[:,i], (Z[:,idx]+ENT[:,i]) )
    # end

    # TODO -- this should be the covariate error from Distributions, only using diagonals here (approxmition for speed in first implementation)
    # dd = size(Z,1) == r
    ENT = randn(r,c)
    for d in 1:r
       @fastmath @inbounds ENT[d,:] = ENT[d,:].*odom.Cov[d,d]
    end
    # Repeat ENT values for new modes from meas
    for j in 1:cz
      for i in 1:c
          z = Z[1:r,j].+ENT[1:r,i]
          # addPose2Pose2!(RES[1:r,i*j], Xval[1:r,i], z ) # doesn't seem to work between modules
          RES[1:r,i*j] = addPose2Pose2(Xval[1:r,i], z )
          #sum(abs(RES[1:r,i*j])) < 1e-14 ? error("evalPotential(odom::Pose2Pose2..) -- RET col is zero, z=$(z), X=$(Xval[1:r,i])") : nothing
      end
    end
    #sum(abs(RES)) < 1e-14 ? error("evalPotential(odom::Pose2Pose2..) -- an input is zero") : nothing

    return RES
end

# Bearing-Range residual function seems to work just fine
# but not through the NLsolve optimizer. Wrap around issue
# function bearrang!(residual::Array{Float64,1}, Z::Array{Float64,1}, X::Array{Float64,1}, L::Array{Float64,1})
#     #@show X, L
#     dx = L[1] - X[1]
#     dy = L[2] - X[2]
#     b = atan2(dy,dx)
#     r = sqrt(dx^2 + dy^2)
#     residual[1] = Z[1] - b + X[3]
#     residual[2] = Z[2] - r
#     nothing
# end

function bearrang!(residual::Array{Float64,1}, Z::Array{Float64,1}, X::Array{Float64,1}, L::Array{Float64,1})
  wTb = SE2(X)
  bTl = wTb\[L[1:2];1.0]
  b = atan2(bTl[2],bTl[1])
  residual[1] = Z[2]-norm(bTl[1:2])
  residual[2] = Z[1]-b
  nothing
end

function pack3(xL1, xL2, p1, p2, p3, xF3)
    X = zeros(3)
    X[p1] = xL1
    X[p2] = xL2
    X[p3] = xF3
    return X
end

function solveLandm(Zbr::Array{Float64,1}, par::Array{Float64,1}, init::Array{Float64,1})
    return (nlsolve(   (l, res) -> bearrang!(res, Zbr, par, l), init )).zero
end

function solvePose2(Zbr::Array{Float64,1}, par::Array{Float64,1}, init::Array{Float64,1})
    p = collect(1:3);
    shuffle!(p);
    p1 = p.==1; p2 = p.==2; p3 = p.==3
    #@show init, par
    r = nlsolve(    (x, res) -> bearrang!(res, Zbr,
                    pack3(x[1], x[2], p1, p2, p3, init[p3]), par),
                    [init[p1];init[p2]] )
    return pack3(r.zero[1], r.zero[2], p1, p2, p3, init[p3]);
end

function solveSetSeps(fnc::Function, Zbr::Array{Float64,1}, CovZ::Array{Float64,2},
                      pars::Array{Float64,2}, inits::Array{Float64,2})
    out = zeros(size(inits))
    for i in 1:size(pars,2)
        ent = 1.0*[CovZ[1,1]*randn(); CovZ[2,2]*randn()]
        out[:,i] = fnc((Zbr+ent), vec(pars[:,i]), vec(inits[:,i]) )
    end
    return out
end

# Xid is the one you want to get back
function evalPotential(brpho::Pose2DPoint2DBearingRange, Xi::Array{Graphs.ExVertex,1}, Xid::Int64)
    # TODO -- add null hypothesis here, might even be done one layer higher in call stack
    val = Array{Float64,2}()
    ini = Array{Graphs.ExVertex,1}()
    par = Array{Graphs.ExVertex,1}()
    oth = Array{Graphs.ExVertex,1}()
    ff::Function = +
    nullhyp = 0.0
    # mmodes = 1 < length(brpho.W)
    # implicit equation portion -- multi-dependent function
    if Xid == Xi[1].index # brpho. ## find the pose

        ff = solvePose2
        # ini = brpho.Xi[1]
        par = Xi[2:end]
        for j in 1:(length(Xi)-1)
            push!(ini, Xi[1])
        end
        #println("Xid == brpho.Xi[1].index inits=", size(inits), ", par=", size(pars))
    elseif Xid == Xi[2].index # find landmark
        if length(Xi) > 2
            nullhyp = 0.5 # should be variable weight
            oth = Xi[3]
        end
        ff = solveLandm
        ini = Xi[2]
        par = Xi[1]
    elseif length(Xi) > 2
        nullhyp = 0.5 # should be variable weight
        if Xid == Xi[3].index # find second mode landmark
            ff = solveLandm
            ini = Xi[3]
            oth = Xi[2]
            par = Xi[1]
        end
    end
    if ff == +
        error("Bad evalPotential Pose2Point2DBearingRange")
    end

    # Gamma = Categorical(brpho.W)

    inits = getVal(ini)
    pars = getVal(par)
    others =  length(Xi) > 2 && Xid != Xi[1].index ? getVal(oth) : Union{}
    # add null hypothesis case
    len = length(Xi) > 2 && Xid != Xi[1].index ? size(others,2) : 0
    # gamma = mmodes ? rand(Gamma) : 1
    numnh = floor(Int, 2*nullhyp*len) # this doubles the value count for null cases
    nhvals = zeros(size(inits,1),numnh)
    for i in 1:numnh
        idx = floor(Int,len*rand()+1)
        # nhvals[:,i] = inits[:,idx] # WRONG! this should be the other landmark, not current landmark
        nhvals[:,i] = others[:,idx]
    end

    val = solveSetSeps(ff, vec(brpho.Zij[:,1]), brpho.Cov, pars, inits)

    return [val';nhvals']'
end

function evalPotentialNew(brpho::Pose2DPoint2DBearingRange, Xi::Array{Graphs.ExVertex,1}, Xid::Int64)
    # TODO -- add null hypothesis here, might even be done one layer higher in call stack
    val = Array{Float64,2}()
    ini = Array{Graphs.ExVertex,1}()
    par = Array{Graphs.ExVertex,1}()
    ff::Function = +
    mmodes = 1 < length(brpho.W)
    x = getVal(brpho.Xi[1])
    l1 = getVal(brpho.Xi[2])
    l2 = mmodes ? getVal(brpho.Xi[3]) : nothing
    nPts = size(x,2)

    pars = Array{Float64,2}(size(x))

    # discrete option, marginalized out before message is sent
    Gamma = mmodes ? rand(Categorical(brpho.W),nPts) : ones(Int,nPts)
    L1s = Gamma .== 1
    nl1s = sum(map(Int, L1s))
    nl2s = nPts-nl1s
    # implicit equation portion -- multi-dependent function

    if Xid == brpho.Xi[1].index # find the pose
        ff = solvePose2
        # par = brpho.Xi[2:end]
        push!(ini, brpho.Xi[1])
        for j in 1:nPts
            pars[:,j] = L1s[j] ? l1[:,j] : l2[:,j]
        end
    elseif Xid == brpho.Xi[2].index # find landmark
        ff = solveLandm
        pars = x
        push!(ini,brpho.Xi[2])
    elseif mmodes
        if Xid == brpho.Xi[3].index # find second mode landmark
            ff = solveLandm
            push!(ini,brpho.Xi[3])
            pars = x
        end
    end
    if ff == +
        error("Bad evalPotential Pose2Point2DBearingRange")
    end

    inits = getVal(ini)

    L1 = sample(getKDE(brpho.Xi[2]),numsampls) #?? # TODO you where here

    len = size(inits,2)
    numnh = floor(Int, 2*nullhyp*len) # this doubles the value count for null cases
    nhvals = zeros(size(inits,1),numnh)
    for i in 1:numnh
        idx = floor(Int,len*rand()+1)
        nhvals[:,i] = inits[:,idx]
    end

    val = solveSetSeps(ff, vec(brpho.Zij[:,1]), brpho.Cov, pars, inits)

    return [val';nhvals']'
end


text:
# test for variations in DOF solving,  #227 #316 #430

# using Revise

using IncrementalInference
using Test

@testset "Test for variations in N while solving" begin

fg = generateCanonicalFG_CaesarRing1D(graphinit=true)

pts = approxConv(fg, :x0x1f1, :x1, N=101)
if size(pts,2) != 101
  @warn "approxConv not adhering to N=101 != $(size(pts,2)), see issue #105"
end

# 150 is non-standard
getSolverParams(fg).N = 150
# getSolverParams(fg).multiproc = false
# getSolverParams(fg).async = false
tree, smt, hist = solveTree!(fg)

@test getKDE(fg, :x1) |> getPoints |> size == (1,150)

# test with change for second solve
getSolverParams(fg).N = 200
tree, smt, hist = solveTree!(fg)

@test getKDE(fg, :x1) |> getPoints |> size == (1,200)

getSolverParams(fg).N = 99
tree, smt, hist = solveTree!(fg)

@test getKDE(fg, :x1) |> getPoints |> size == (1,99)

end

#


text:
using Fst

a = Wfst()
add_arc!(a, 0, 1, "x", "x", 1.0)
add_arc!(a, 1, 2, "x", "x", 1.0)
add_arc!(a, 1, 3, "x", "x", 1.0)
add_arc!(a, 3, 4, "x", "x", 1.0)
add_arc!(a, 4, 3, "x", "x", 1.0)
add_arc!(a, 1, 5, "x", "x", 1.0)
add_arc!(a, 5, 1, "x", "x", 1.0)
add_arc!(a, 6, 1, "x", "x", 1.0)
add_arc!(a, 6, 7, "x", "x", 1.0)
add_arc!(a, 7, 6, "x", "x", 1.0)
add_initial_state!(a, 0, 1.0)
add_final_state!(a, 2, 1.0)
create_pdf(a, "a.pdf")
println(a)
create_pdf(accessible(a), "A(a).pdf")
create_pdf(coaccessible(accessible(a)), "Co(A(a)).pdf")
create_pdf(coaccessible(accessible(a)), "trim(a).pdf")


text:
# script that outputs some statistics of the model fits

using Printf, Random, DataFrames, CSV

include("common.jl")
include("models.jl")
include("behaviorstats.jl")
include("simmodel.jl")
include("modelfitting.jl")

# settings (when manually evaluating the objective function)
const blockreps = 10
const simtrials = 10^5
const burnin = 1000
const nrats = 1 # do not divide by nrats = 4
const δt = 0.0001
const datapath = "../fitdata/fits"

# models
const modeldict = Dict("opt" => OptimalModel,
                       "optlapse" => OptimalLapseModel,
                       "optcollapse" => OptimalColLapseModel,
                       "optcolbiaslapse" => OptimalColBiasLapseModel,
                       "adflapse" => ADFLapseModel,
                       "adfcollapse" => ADFColLapseModel,
                       "adfbiaslapse" => ADFBiasLapseModel,
                       "adfcolbiaslapse" => ADFColBiasLapseModel,
                       "randomcolbiaslapse" => RandomColBiasLapseModel,
                       "randomcollapse" => RandomColLapseModel,
                       "gammalapse" => GammaLapseModel,
                       "gammacollapse" => GammaColLapseModel,
                       "gammabiaslapse" => GammaBiasLapseModel,
                       "gammacolbiaslapse" => GammaColBiasLapseModel,
                       "deltalapse" => DeltaLapseModel,
                       "deltacollapse" => DeltaColLapseModel,
                       "deltabiaslapse" => DeltaBiasLapseModel,
                       "deltacolbiaslapse" => DeltaColBiasLapseModel,
                       "deltaexpcolbiaslapse" => DeltaExpColBiasLapseModel,
                       "lcalapse" => LCALapseModel,
                       "lcacollapse" => LCAColLapseModel)
# separate list to preserve order that keys(modeldict) wouldn't have
const modelnames = ["opt", "optlapse", "optcollapse", "optcolbiaslapse",
                    "adflapse", "adfcollapse", "adfbiaslapse", "adfcolbiaslapse",
                    "randomcolbiaslapse", "randomcollapse",
                    "gammalapse", "gammacollapse", "gammabiaslapse", "gammacolbiaslapse",
                    "deltalapse", "deltacollapse", "deltabiaslapse", "deltacolbiaslapse",
                    "deltaexpcolbiaslapse", "lcalapse", "lcacollapse"]
@assert issubset(Set(modelnames), Set(keys(modeldict))) # unordered comparison

# objective functions
# second tuple component is trials, 1 = ident, 2 = categ, 3 = itlvd
const objfndict = Dict("psychchron1" => (PsychChron1ObjFn, [1, 2]),
                       "psychchron1ident" => (PsychChron1IdentObjFn, [1]),
                       "psychchron1categ" => (PsychChron1CategObjFn, [2]),
                       "psychchron2" => (PsychChron2ObjFn, [1, 2]),
                       "psychchron3" => (PsychChron3ObjFn, [1, 2]),
                       "psychchronseq1" => (PsychChronSeq1ObjFn, [1, 2]),
                       "psychchronseq2" => (PsychChronSeq2ObjFn, [1, 2]),
                       "psychchronseq3" => (PsychChronSeq3ObjFn, [1, 2]),
                       "psychchron1itlvd" => (PsychChron1ItlvdObjFn, [3]))
#const objfns = ["psychchron1", "psychchron1ident", "psychchron2", "psychchron3",
#                "psychchronseq1", "psychchronseq2", "psychchronseq3"]
const objfns = ["psychchron1", "psychchron1ident", "psychchron1categ", "psychchronseq1"]
const objfnsitlvd = ["psychchron1itlvd"]
const allobjfns = [objfns; objfnsitlvd]
@assert issubset(Set(objfns), Set(keys(objfndict)))

# model fit criteria, converted to log-evidence in log(10)
BIC(avgllh, ϕn, trials) = (trials * avgllh - 0.5ϕn * log(trials)) / log(10)
AIC(avgllh, ϕn, trials) = (trials * avgllh - ϕn) / log(10)
AICc(avgllh, ϕn, trials) = AIC(avgllh, ϕn, trials) - (ϕn * (ϕn+1) / (trials-ϕn-1)) / log(10)

# returns number of trials in categ/ident dataset as vector
datatrials() = [size(d, 1) / nrats for d in (load_data(:categ), load_data(:ident), load_data(:itlvd))]

# returns dataframe with best-fitting parameters and avg. log-likelihood
function loadmodelfits(modelname, objfnname)
    fitfilename = "$(datapath)/fit_$(modelname)_$(objfnname).csv"
    if !isfile(fitfilename)
        error("Could not find fit file $fitfilename")
    end
    d = CSV.read(fitfilename)
    # identify max-llh parameters
    maxllh, maxi = findmax(d[:llh])
    if maxllh != d[end,:llh]
        warning("Last element in $fitfilename is not best-fitting")
    end
    return d[maxi,:]
end

# prints the fit parameters for all objective functions for a given model
function printfitparams(modelname)
    # use first objective function to get parameter names
    d = loadmodelfits(modelname, objfns[1])
    ϕnames = names(d)
    # print parameters
    objfnlen = maximum([length(s) for s in allobjfns])
    println(repeat("-", objfnlen), "  ",
            join([@sprintf("%7s", string(ϕi)) for ϕi in ϕnames], "  "))
    # print fits
    for objfn in allobjfns
        d = loadmodelfits(modelname, objfn)
        println(repeat(" ", objfnlen - length(objfn)), objfn, "  ",
                join([@sprintf("%7.2f", x) for x in Vector(d)], "  "))
    end
end

# prints the fit quality for all models for the given objective function
function printfitquality(objfn, trials)
    modellen = maximum([length(s) for s in modelnames])
    println(repeat("-", modellen), "  ",
            join([@sprintf("%10s", s) for s in ["avg. llh", "BIC", "AIC", "AICc"]], "  "))
    for modelname in modelnames
        d = loadmodelfits(modelname, objfn)
        avgllh = d[:llh]
        ϕn = length(d) - 1
        fitq = [avgllh; [f(avgllh, ϕn, trials) for f in [BIC, AIC, AICc]]]
        println(repeat(" ", modellen - length(modelname)), modelname, "  ",
                join([@sprintf("%10.2f", x) for x in fitq], "  "))
    end
end

# plots the fit quality for all models using the parameters corresponding to
# paramobjfn and uses the given objective function type objfn to compute the
# performance for all models
function printsepfitquality(paramobjfn, evalobjfn)
    # trial sequence to fit data on
    d_ident = load_data(:ident)
    d_categ = load_data(:categ)
    Random.seed!(0)
    evidence, ids = catevidence(ident_cs, categ_cs, 10, 10^5)
    # iterate over all models
    modellen = maximum([length(s) for s in modelnames])
    println(repeat("-", modellen), "  ",
            join([@sprintf("%10s", s) for s in ["avg. llh", "BIC", "AIC", "AICc"]], "  "))
    evalobjfntype, evalobjfntrials = objfndict[evalobjfn]
    evalobjfntrials = sum([size(d_ident, 1), size(d_categ, 1)][evalobjfntrials]) / nrats
    for modelname in modelnames
        # instantiate model and extract parameters
        m = modeldict[modelname]()
        ϕ = Vector(loadmodelfits(modelname, paramobjfn)[ϕnames(m)])
        ϕn = length(ϕ)
        # evaluate objective function
        o = evalobjfntype(m, evidence, ids, burnin, d_ident, d_categ, nrats)
        avgllh = llh(o, ϕ)
        fitq = [avgllh; [f(avgllh, ϕn, evalobjfntrials) for f in [BIC, AIC, AICc]]]
        println(repeat(" ", modellen - length(modelname)), modelname, "  ",
                join([@sprintf("%10.2f", x) for x in fitq], "  "))
    end
end

# print model parameters
for modelname in modelnames
    println("Model $modelname")
    printfitparams(modelname)
    println()
end
println()

# print model fit statistics
const trials = datatrials()
for objfn in allobjfns
    println("Objective function $objfn")
    printfitquality(objfn, sum(trials[objfndict[objfn][2]]))
    println()
end

# print model fit statistics when evaluated with psychchronseq1 objective
for objfn in objfns
    println("Fitted on $objfn, evaluated on psychchronseq1")
    printsepfitquality(objfn, "psychchronseq1")
    println()
end


text:
import Base: getindex, size, summary, show
using Base: dims2string
using AutoGrad: Value

# https://docs.julialang.org/en/v1/manual/types/#man-custom-pretty-printing-1
# show(io::IO, z): single line format used in show, print, inside other objects.
# show(io::IO, ::MIME"text/plain", z): multi-line format used by display.
# show(io::IO, ::MIME"text/html", z): multi-line format for html output.
# get(io, :compact, false), show(IOContext(stdout, :compact=>true),z) for compact (array) printing.
# summary(io::IO, x) = print(io, typeof(x))
# string(z): uses print_to_string.

# Hack for printing without copying the whole KnetArray and without inheriting AbstractArray:
struct KnetDisplay{T,N} <: AbstractArray{T,N}; a::KnetArray{T,N}; end
getindex(a::KnetDisplay, i...) = getindex(a.a, i...)
size(a::KnetDisplay) = size(a.a)
summary(io::IO, a::KnetDisplay) = summary(io, a.a)
summary(io::IO, a::KnetArray) = print(io, dims2string(size(a)), " ", typeof(a))
show(io::IO, m::MIME"text/plain", a::KnetArray) = show(io, m, KnetDisplay(a))
summary(io::IO, x::Value{A}) where {A<:KnetArray} = print(io, dims2string(size(x)), " ", typeof(x))

function show(io::IO, a::KnetArray) # Compact display used by print
    T = eltype(a)
    print(io, T <: AbstractFloat ? "K$(sizeof(T)*8)" : "K{$T}")
    print(io, "($(join(size(a),',')))")
    print(io, isempty(a) ? "[]" : "[$(a[1])⋯]")
end


text:
# HellwarthEffectiveFrequency.jl
#   - use Hellwarth et al. 1999 PRB method to reduce multiple phonon modes to a single effective frequency

@testset "HellwarthEffectiveFrequency" begin

# ((freq THz)) ((IR Activity / e^2 amu^-1))
# These data from MAPbI3-Cubic_PeakTable.csv
# https://github.com/WMD-group/Phonons/tree/master/2015_MAPbI3/SimulatedSpectra
# Data published in Brivio2015 (PRB)
# https://doi.org/10.1103/PhysRevB.92.144308
MAPI= [
96.20813558773261 0.4996300522819191
93.13630357703363 1.7139631746083817
92.87834578121567 0.60108592692181
92.4847918585963 0.0058228799414729
92.26701437594754 0.100590086574602
89.43972834606603 0.006278895133832249
46.89209141511332 0.2460894564364346
46.420949316788 0.14174282581124137
44.0380222871706 0.1987196948553428
42.89702947649343 0.011159939465770681
42.67180170168193 0.02557751102757614
41.46971205834201 0.012555230726601503
37.08982543385215 0.00107488277468418
36.53555265689563 0.02126940080871224
30.20608114002676 0.009019481779712388
27.374810898415028 0.03994453721421388
26.363055017011728 0.05011922682554448
9.522966890022039 0.00075631870522737
4.016471586720514 0.08168931020200264
3.887605410774121 0.006311654262282101
3.5313112232401513 0.05353548710183397
2.755392921480459 0.021303020776321225
2.4380741812443247 0.23162784335484837
2.2490917637719408 0.2622203718355982
2.079632190634424 0.23382298607799906
2.0336707697261187 0.0623239656843172
1.5673011873879714 0.0367465760261409
1.0188379384951798 0.0126328938653956
1.0022960504442775 0.006817361620021601
0.9970130778462072 0.0103757951973341
0.9201781906386209 0.01095811116040592
0.800604081794174 0.0016830270365341532
0.5738689505255512 0.00646428491253749
#0.022939578929507105 8.355742795827834e-05   # Acoustic modes!
#0.04882611767873102 8.309858592685e-06
#0.07575149723846182 2.778248540373041e-05
]

MAPbBr3 = [
#v / THz IR Activity / e^2 amu^-1
96.25494581101505 0.4856105419754306
93.36827101597564 1.5237045178260684
93.19828072025739 0.54951133126244
92.59570097621418 0.0017497942636829
92.36075758752216 0.054483693922451705
89.49998024367731 0.003950754031148851
47.175408954380885 0.21156964115237242
46.57576693654983 0.11627478349663156
44.3361503231313 0.19163912283684897
42.967352451897554 0.011987727827091217
42.71066583604171 0.027797225985699996
41.490751582679316 0.01234463397562453
37.22033177431264 0.001649578352399309
36.57010484641483 0.019161684306162618
30.438999700507342 0.007365167772350685
27.625465570185764 0.035161502048866464
26.458785163031635 0.038174133532805965
9.947074760143552 0.0012355225072173002
4.493935108724231 0.07714630445769342
4.215763207131703 0.01724042408702847
4.021795301113039 0.06167913225780508
3.2393471899404047 0.040335730459227545
2.914057478007975 0.27361583921530497
2.6223962950449407 0.23493646364819373
2.400751729577614 0.2753271857365242
2.370809478597967 0.04428317547154772
2.0388405406936783 0.0517120158657721
1.4017373572275915 0.0046510016948825
1.3634089519825427 0.0018248126743814707
1.336413930635821 0.009190430969212338
1.2144969737450764 0.013885351107636132
1.0763905048534899 0.0014676450708402083
0.7427782546143306 0.008829646844819504
#0.012830037764855392 6.62960832415985e-06
#0.02156203260482768 4.708045983230229e-05
#0.05156888885824344 2.602755859934495e-05
]

MAPbCl3= [
#v / THz IR Activity / e^2 amu^-1
96.65380905798915 0.4542931751264273
93.30434573145013 1.4046841552487999
93.27126600254645 0.5286008698609284
92.719992888923 0.0015896642905404556
92.45885414527524 0.0866990370170258
89.61967951587117 0.0031392395308660263
47.44029120639781 0.18901278025443116
46.764089990499585 0.09717407330258415
44.60198174803584 0.19436643102794884
43.04971569182655 0.012986533289171532
42.7637228646608 0.02998334671150052
41.49924134521671 0.010714577213688365
37.34275112339727 0.0028788161278648656
36.62941263763097 0.018016387094969555
30.651112964494686 0.006550727272302003
27.89469865012452 0.03515737503204058
26.603746656189923 0.03497660426639104
10.493140658721968 9.69032373922839e-05
5.242459815070028 0.12301715153349158
4.8348539006445845 0.11307729502447525
4.633936235943856 0.009371996258413128
4.1582548301810665 0.08368265471427451
3.9021555976797266 0.44406366947684967
3.3808364269416287 0.3379317982808455
3.099250159114845 0.38405867852912545
2.8466202448045217 0.025463767673713028
2.5900539956188333 0.07754012506771049
2.0931563078880493 0.008533555182635471
1.869999049657922 0.011347668353507898
1.8382471614901006 0.0005371757988
1.829843559303894 0.019270325513936128
1.6503257694697964 0.000950380841742176
1.0962868606064022 0.018510433944796113
#0.032159306112378404 2.774125107895896e-05 #Acoustic
#0.04846756590412458 7.796048392205337e-05
#0.05253089286287104 2.8730156133e-06
]

MAPI_low=MAPI[19:33,:] # Just inorganic components, everything below 10THz; modes 3-18

println("\n\nMAPI: BScheme (athermal)")
println("\t MAPI: (all values)")
HellwarthBScheme(MAPI)
println("\t MAPI: (low-frequency, non molecular IR, only)")

Bscheme=HellwarthBScheme(MAPI_low)
println(Bscheme, " ~= 2.25 THz")
@test Bscheme ≈ 2.25 atol=0.01

println("\n\nHellwarth1999 (party like it's)")

Hellwarth=
[
106.23 8.86
160.51 9.50
180.33 20.85
206.69 10.05
252.76 27.00
369.64 61.78
501.71 52.87
553.60 86.18
585.36 75.41
607.29 98.15
834.53 89.36
]

println("\nHellwarth B scheme... Hellwarth data...")
Bscheme=HellwarthBScheme(Hellwarth)
println(Bscheme, " ~= W 196.9 cm^-1  Ω = 500 cm^-1")

println("\nHellwarth A scheme T=0... Hellwarth data...")
Ascheme=HellwarthAScheme(Hellwarth,T=0)

println("\nHellwarth A scheme T=295... Hellwarth data...")
Ascheme=HellwarthAScheme(Hellwarth,T=295)
println(Ascheme, " ~= W 196.9 cm^-1   Ω = 504 cm^-1")

end


text:
# Autogenerated wrapper script for bitwarden_cli_jll for x86_64-w64-mingw32
export bw

JLLWrappers.@generate_wrapper_header("bitwarden_cli")
JLLWrappers.@declare_executable_product(bw)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        bw,
        "bin\\bw.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
struct ChoiceRecord{T}
    retval::T
    score::Float64
end

struct CallRecord{T}
    subtrace::T
    score::Float64
    noise::Float64
end

struct ChoiceOrCallRecord{T}
    subtrace_or_retval::T
    score::Float64
    noise::Float64 # if choice then NaN
    is_choice::Bool
end

function ChoiceRecord(record::ChoiceOrCallRecord)
    if !record.is_choice
        error("Found call but expected choice")
    end
    ChoiceRecord(record.subtrace_or_retval, record.score)
end

function CallRecord(record::ChoiceOrCallRecord)
    if record.is_choice
        error("Found choice but expected call")
    end
    CallRecord(record.subtrace_or_retval, record.score, record.noise)
end

mutable struct DynamicDSLTrace{T} <: Trace
    gen_fn::T
    trie::Trie{Any,ChoiceOrCallRecord}
    isempty::Bool
    score::Float64
    noise::Float64
    args::Tuple
    retval::Any
    function DynamicDSLTrace{T}(gen_fn::T, args) where {T}
        trie = Trie{Any,ChoiceOrCallRecord}()
        # retval is not known yet
        new(gen_fn, trie, true, 0, 0, args)
    end
end

DynamicDSLTrace(gen_fn::T, args) where {T} = DynamicDSLTrace{T}(gen_fn, args)

set_retval!(trace::DynamicDSLTrace, retval) = (trace.retval = retval)

function has_choice(trace::DynamicDSLTrace, addr)
    haskey(trace.trie, addr) && trace.trie[addr].is_choice
end

function has_call(trace::DynamicDSLTrace, addr)
    haskey(trace.trie, addr) && !trace.trie[addr].is_choice
end

function get_choice(trace::DynamicDSLTrace, addr)
    choice = trace.trie[addr]
    if !choice.is_choice
        throw(KeyError(addr))
    end
    ChoiceRecord(choice)
end

function get_call(trace::DynamicDSLTrace, addr)
    call = trace.trie[addr]
    if call.is_choice
        throw(KeyError(addr))
    end
    CallRecord(call)
end

function add_choice!(trace::DynamicDSLTrace, addr, retval, score)
    if haskey(trace.trie, addr)
        error("Value or subtrace already present at address $addr.
            The same address cannot be reused for multiple random choices.")
    end
    trace.trie[addr] = ChoiceOrCallRecord(retval, score, NaN, true)
    trace.score += score
    trace.isempty = false
end

function add_call!(trace::DynamicDSLTrace, addr, subtrace)
    if haskey(trace.trie, addr)
        error("Value or subtrace already present at address $addr.
            The same address cannot be reused for multiple random choices.")
    end
    score = get_score(subtrace)
    noise = project(subtrace, EmptySelection())
    submap = get_choices(subtrace)
    trace.isempty = trace.isempty && isempty(submap)
    trace.trie[addr] = ChoiceOrCallRecord(subtrace, score, noise, false)
    trace.score += score
    trace.noise += noise
end

###############
# GFI methods #
###############

get_args(trace::DynamicDSLTrace) = trace.args
get_retval(trace::DynamicDSLTrace) = trace.retval
get_score(trace::DynamicDSLTrace) = trace.score
get_gen_fn(trace::DynamicDSLTrace) = trace.gen_fn

## get_choices ##

function get_choices(trace::DynamicDSLTrace)
    if !trace.isempty
        DynamicDSLChoiceMap(trace.trie) # see below
    else
        EmptyChoiceMap()
    end
end

struct DynamicDSLChoiceMap <: ChoiceMap
    trie::Trie{Any,ChoiceOrCallRecord}
end

get_address_schema(::Type{DynamicDSLChoiceMap}) = DynamicAddressSchema()
Base.isempty(::DynamicDSLChoiceMap) = false # TODO not necessarily true
has_value(choices::DynamicDSLChoiceMap, addr::Pair) = _has_value(choices, addr)
get_value(choices::DynamicDSLChoiceMap, addr::Pair) = _get_value(choices, addr)
get_submap(choices::DynamicDSLChoiceMap, addr::Pair) = _get_submap(choices, addr)

function get_submap(choices::DynamicDSLChoiceMap, addr)
    trie = choices.trie
    if has_leaf_node(trie, addr)
        # leaf node, must be a call
        call = trie[addr]
        if call.is_choice
            throw(KeyError(addr))
        end
        get_choices(call.subtrace_or_retval)
    elseif has_internal_node(trie, addr)
        # internal node
        subtrie = get_internal_node(trie, addr)
        DynamicDSLChoiceMap(subtrie) # see below
    else
        EmptyChoiceMap()
    end
end

function has_value(choices::DynamicDSLChoiceMap, addr)
    trie = choices.trie
    has_leaf_node(trie, addr) && trie[addr].is_choice
end

function get_value(choices::DynamicDSLChoiceMap, addr)
    trie = choices.trie
    choice = trie[addr]
    if !choice.is_choice
        throw(KeyError(addr))
    end
    choice.subtrace_or_retval
end

function get_values_shallow(choices::DynamicDSLChoiceMap)
    ((key, choice.subtrace_or_retval)
     for (key, choice) in get_leaf_nodes(choices.trie)
     if choice.is_choice)
end

function get_submaps_shallow(choices::DynamicDSLChoiceMap)
    calls_iter = ((key, get_choices(call.subtrace_or_retval))
        for (key, call) in get_leaf_nodes(choices.trie)
        if !call.is_choice)
    internal_nodes_iter = ((key, DynamicDSLChoiceMap(trie))
        for (key, trie) in get_internal_nodes(choices.trie))
    Iterators.flatten((calls_iter, internal_nodes_iter))
end

## Base.getindex ##

function _getindex(trace::DynamicDSLTrace, trie::Trie, addr::Pair)
    (first, rest) = addr
    if haskey(trie.leaf_nodes, first)
        choice_or_call = trie.leaf_nodes[first]
        if choice_or_call.is_choice
            error("Unknown address $addr; random choice at $first")
        else
            subtrace = choice_or_call.subtrace_or_retval
            return subtrace[rest]
        end
    elseif haskey(trie.internal_nodes, first)
        return _getindex(trace, trie.internal_nodes[first], rest)
    end
end

function _getindex(trace::DynamicDSLTrace, trie::Trie, addr)
    if haskey(trie.leaf_nodes, addr)
        choice_or_call = trie.leaf_nodes[addr]
        if choice_or_call.is_choice
            # the value of the random choice
            return choice_or_call.subtrace_or_retval
        else
            # the return value of the generative function call
            return get_retval(choice_or_call.subtrace_or_retval)
        end
    elseif haskey(trie.internal_nodes, addr)
        error("No random choice or generative function call at address $addr")
    end
end

function Base.getindex(trace::DynamicDSLTrace, addr)
    _getindex(trace, trace.trie, addr)
end


text:
##############################################
# This file has been automatically generated #
# Do not edit !                              #
##############################################

# TODO: pending https://github.com/JuliaLang/julia/issues/29420

const PNGCAPI = nothing
const PNG_BYTES_TO_CHECK = 8
const tm = Libc.TmStruct
const time_t = UInt
const jmp_buf = Nothing


text:
# This file was generated, do not modify it. # hide
function WtoV(W)
    Q = Matrix(first(qr(W)))
    L = first(lq(Q))
    V = L * Diagonal(sign.(diag(L)))
    return V
end

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
exp_bits(n) = Dict(32 => 8, 64 => 11)[n]
exp_shift(n) = n - exp_bits(n) - 1
exp_bias(n) = (1 << (exp_bits(n) - 1)) - 1
exp_mask(n) = top_bits(n) - (top_bits(n) >> exp_bits(n))

frac_mask(n)  = (top_bits(n) >> 1) - one(bits_to_UInt[n])
guard_mask(n, fpsize) = top_bits(fpsize) >> n
inner_mask(n, fpsize) = top_bits(fpsize) >> (n - 1)
summ_mask(n, fpsize) = guard_mask(n, fpsize) - one(bits_to_UInt[n])

function floatconvert(n, es)
  maximum_exponent = (n - 2) * (2 ^ es)
  minimum_exponent = -((n - 1) * (2 ^ es))

  ftype = c_ftype(float_for_posit_size(n))

  fcall = (es == 0) ? "$(ftype)_to_p$(n)_zero_es(fval)" : "$(ftype)_to_p$(n)(fval, $(es), $maximum_exponent, $minimum_exponent)"

  """
$(f2pheader(n,es)){
  $(p(n,es)) res;
  res.udata = $(fcall);
  return res;
}
  """
end

#printf specifiers
const ps = Dict(32 => "%08x", 64 => "%016lx")

function floatconvert(n; zero_es = false)

  fpsize = float_for_posit_size(n)  #for now.
  ftype = c_ftype(fpsize)

  es_fn = zero_es ? "_zero_es" : ""
  es_hd = zero_es ? "" : ", int16_t es, int16_t maximum_exponent, int16_t minimum_exponent"

  max_exp = zero_es ? n - 2 : "maximum_exponent"
  min_exp = zero_es ? 1 - n : "minimum_exponent"

  es_expfrc = zero_es ? """
  //use an uint$(fpsize)_t value as an intermediary store for
    //raw fraction bits.  Shift all the way to the right and then

    int16_t regime = exponent;
    uint$(fpsize)_t frac = (*ival) << ($(exp_bits(fpsize) - 1));
    //there are no exponent bits.
    frac &= $(c_literal((top_bits(fpsize) >> 1) - one(bits_to_UInt[fpsize])));
  """ : """
  //divide up this exponent into a proper exponent and regime.
    int16_t regime = exponent >> es;
    exponent = exponent & ((1 << es) - 1);
    //use an uint$(fpsize)_t value as an intermediary store for
    //all off the fraction bits, initially backing off by es.  Mask out the top
    //two bits.
    uint$(fpsize)_t frac = ((*ival) << ($(exp_bits(fpsize) - 1) - es)) & ($(c_literal(frac_mask(fpsize))) >> es);
    //append the exponent bits to frac representation.
    frac |= ((uint$(fpsize)_t) exponent) << ($(fpsize - 2) - es);
  """

  """
static uint$(n)_t $(ftype)_to_p$(n)$(es_fn)($(ftype) fval$(es_hd)){
  //create a result value
  //infinity and NaN checks:
  if (isinf(fval)) {return $(c_literal(top_bits(n)));};
  if (fval == 0)   {return $(c_literal(zero_bits(n)));};
  if (isnan(fval)) {
    if (POSIT_ENV.nanmode)
      { return P$(n)INF; }
    else
      { printf(\"bad: %f \\n\", fval); errno = EDOM; return P$(n)INF; }
  }
  //do a surreptitious conversion from $(ftype) precision to UInt$(n)
  uint$(fpsize)_t *ival = (uint$(fpsize)_t *) &fval;
  bool signbit = (($(c_literal(top_bits(fpsize))) & (*ival)) != 0);
  //capture the exponent value
  int16_t exponent = ((($(c_literal(exp_mask(fpsize))) & (*ival)) >> $(exp_shift(fpsize))) - $(exp_bias(fpsize)));

  //pin the exponent.
  exponent = (exponent > $(max_exp)) ? $(max_exp) : exponent;
  //underflow behavior is defined by the POSIT_ENV.underflows setting.
  exponent = (exponent >= $(min_exp) || POSIT_ENV.underflows) ? exponent : $(min_exp);

  $(es_expfrc)

  //next, append the appropriate shift prefix in front of the value.
  int16_t shift;

  if (regime >= 0) {
    shift = 1 + regime;
    frac |= $(c_literal(top_bits(fpsize)));
  } else {
    shift = -regime;
    frac |= $(c_literal(top_bits(fpsize) >> 1));
  }

  //perform an *arithmetic* shift; convert back to unsigned.
  frac = (uint$(fpsize)_t)(((int$(fpsize)_t) frac) >> shift);

  //mask out the top bit of the fraction, which is going to be the
  //basis for the result.
  frac = frac & $(c_literal(top_bits(fpsize) - one(bits_to_UInt[fpsize])));

  bool guard = (frac & $(c_literal(guard_mask(n,fpsize)))) != 0;
  bool summ  = (frac & $(c_literal(summ_mask(n,fpsize))) ) != 0;
  bool inner = (frac & $(c_literal(inner_mask(n,fpsize)))) != 0;

  //round the frac variable in the event it needs be augmented.

  frac += ((guard && inner) || (guard && summ)) ? $(c_literal(inner_mask(n, fpsize))) : $(c_literal(zero_bits(fpsize)));
  //shift as necessary

  //shift further, as necessary, to match sizes
  frac = frac >> $(fpsize - n);

  //check to recast zeros to the smallest value
  frac = (frac != 0 || POSIT_ENV.underflows) ? frac : $(c_literal(one(bits_to_UInt[fpsize])));

  uint$(n)_t sfrac = (uint$(n)_t) frac;

  return (signbit ? -sfrac : sfrac);
}
"""
end

function positconvert(n, es)

  fpsize = float_for_posit_size(n)
  ftype = c_ftype(fpsize)

  """
$(p2fheader(n,es)){
  return p$(n)_to_$(ftype)(pval.udata, $(es), $(c_literal(-(top_bits(n) >> (es - 1)))));
}
  """
end

function positconvert(n)

  fpsize = float_for_posit_size(n)
  ftype = c_ftype(fpsize)
  clz = (n == 64) ? "__builtin_clzl" : "__builtin_clz"

  """
static $(ftype) p$(n)_to_$(ftype)(uint$(n)_t pval, int16_t es, uint$(n)_t es_mask){
  //check for infs and zeros, which do not necessarily play nice with our algorithm.
  if (pval == $(c_literal(top_bits(n))))  return (POSIT_ENV.nanmode) ? NAN : INFINITY;
  if (pval == $(c_literal(zero_bits(n)))) return ($(ftype)) 0;
  //next, determine the sign of the posit value
  bool negative = ((pval & $(c_literal(top_bits(n)))) != 0);
  uint$(n)_t pos_val = negative ? -pval : pval;
  //ascertain if it's inverted.
  bool inverted = (pos_val & $(c_literal(top_bits(n) >> 1))) == 0;
  //note that the clz/clo intrinsics operate on 32-bit data types.
  uint16_t u_regime;
  int16_t s_regime;
  if (inverted){
    //just count the leading zeros, which will tell you the regime.
    u_regime = $(clz)(pos_val)$(n < 32 ? n - 32 : "") - 1;
    s_regime = - u_regime;
  } else {
    //there's no "clo" intrinsic in standard c (whether or not there is a
    //machine opcode) so we have to do this very awkward transformation first.
    uint$(n)_t z_posit = ~pos_val & $(c_literal(~top_bits(n)));
    //__builtin_clz has "undefined" state for a value of 0.  W.T.F, C??
    u_regime = (z_posit == 0) ? $(n - 1) : $(clz)(z_posit)$(n < 32 ? n - 32 : "") - 1;
    s_regime = u_regime - 1;
  }
  //next create the proper exp/frac value by shifting the pos_val, based on the
  //unsigned regime value.
  uint$(n)_t p_exp_frac = (pos_val << (u_regime + 2));

  //extract the exponent value by grabbing the top bits.
  int16_t exponent = (p_exp_frac & es_mask) >> ($(n) - es);

  //append the (signed) regime value to this.
  exponent |= (s_regime << es);
  //finally, add the bias value
  exponent += $(exp_bias(fpsize));

  //shift the fraction again to obliterate the exponent section.
  uint$(n)_t p_frac = p_exp_frac << es;

  uint$(fpsize)_t result;
  result = (negative ? $(c_literal(top_bits(fpsize))) : $(c_literal(zero_bits(fpsize)))) |
    (((uint$(fpsize)_t) exponent) << $(exp_shift(fpsize))) |
    (((uint$(fpsize)_t) p_frac) << ($(exp_shift(fpsize) - n)));

  //convert the result floatint to the desired float type
  return *(($ftype *) &result);
}
"""
end

function generate_posit_conv_cpp(io, posit_defs)
  #generates "posit.h" based on the posit_definitions
  write(io, """
  #include <errno.h>
  #include <stdio.h>
  #include \"include/posit.h\"
  """)

  for n in sort(collect(keys(posit_defs)))

    write(io, """
    /*************************************************************/
    /*  posit_$(n) section, general forms                           */
    /*************************************************************/

    $(floatconvert(n))

    $(floatconvert(n, zero_es = true))

    $(positconvert(n))

    /*************************************************************/
    /*  posit_$(n) section, variable ES adapters                    */
    /*************************************************************/
    """)

    for es in posit_defs[n]
      write(io, """
      extern "C" $(floatconvert(n, es))
      extern "C" $(positconvert(n, es))

      """)
    end
  end
end


text:
mutable struct CompoundLikelihood <: MyXMLElement
    el::XMLOrNothing
    id::String
    likelihoods::Vector{<:MyXMLElement}
    already_made::Vector{Bool}
    function CompoundLikelihood(id::String, likelihoods::Vector{<:MyXMLElement};
                                already_made::Vector{Bool} = fill(true, length(likelihoods)))

        @assert length(already_made) == length(likelihoods)
        return new(nothing, id, likelihoods, already_made)
    end
end

function name(cl::CompoundLikelihood)
    return bn.LIKELIHOOD
end

function get_loggables(cl::CompoundLikelihood)
    make_xml(cl)
    return cl.el
end

function add_likelihood!(cp::CompoundLikelihood, xml::MyXMLElement; already_made::Bool = true)
    push!(cp.likelihoods, xml)
    push!(cp.already_made, already_made)
end

mutable struct MCMCXMLElement <: MyXMLElement
    el::XMLOrNothing
    likelihoods::CompoundLikelihood
    priors::Vector{<:MyXMLElement}
    operators::OperatorsXMLElement
    chain_length::Int
    file_logEvery::Int
    screen_logEvery::Int
    loggables::LoggablesXMLElement
    filename::String
    log_files::Bool
    attrs::Dict{String, String}
    extras::Vector{<:MyXMLElement}

    function MCMCXMLElement(tl_el::TraitLikelihoodXMLElement,
                    mbd_el::MBDXMLElement,
                    rm_el::RepeatedMeasuresXMLElement,
                    os_el::OperatorsXMLElement; chain_length = 100)
        fle = max(1, div(chain_length, 1000))
        sle = 1000
        filename = "defaultFile"
        lg_el = LoggablesXMLElement([mbd_el, rm_el], [true, true])
        attrs = Dict(bn.AUTO_OPTIMIZE => bn.TRUE)

        return new(nothing, CompoundLikelihood("likelihood", [tl_el]), [mbd_el, rm_el], os_el, chain_length,
                    fle, sle, lg_el, filename, true, attrs,
                    MyXMLElement[])
    end

    function MCMCXMLElement(tl_el::TraitLikelihoodXMLElement,
                    mbd_el::MBDXMLElement,
                    rm_el::RepeatedMeasuresXMLElement,
                    os_el::OperatorsXMLElement,
                    loggables::LoggablesXMLElement;
                    chain_length = 100)
        fle = max(1, div(chain_length, 1000))
        sle = 1000
        filename = "defaultFile"
        attrs = Dict(bn.AUTO_OPTIMIZE => bn.TRUE)

        return new(nothing, [tl_el], CompoundLikelihood("likelihood", [mbd_el, rm_el]), os_el, chain_length,
                    fle, sle, loggables, filename, true, attrs,
                    MyXMLElement[])
    end

    function MCMCXMLElement(tl_el::TraitLikelihoodXMLElement,
                    mbd_el::MBDXMLElement,
                    if_el::Union{IntegratedFactorsXMLElement, LatentFactorModelXMLElement},
                    os_el::OperatorsXMLElement; chain_length = 100)
        fle = max(1, div(chain_length, 1000))
        sle = 1000
        filename = "defaultFile"
        lg_el = LoggablesXMLElement([mbd_el, if_el], [true, true])
        attrs = Dict(bn.AUTO_OPTIMIZE => bn.TRUE)

        return new(nothing, CompoundLikelihood("likelihood", [tl_el, if_el]), [mbd_el, if_el], os_el, chain_length,
                    fle, sle, lg_el, filename, true, attrs,
                    MyXMLElement[])
    end
end

function set_options!(mcmc::MCMCXMLElement, options::MCMCOptions)
    mcmc.chain_length = options.chain_length
    mcmc.file_logEvery = options.file_log_every
    mcmc.screen_logEvery = options.screen_log_every
    set_full_eval!(mcmc, options.likelihood_check_count)
end

function set_chain_length!(mcmc::MCMCXMLElement, cl::Int)
    mcmc.chain_length = cl
end

function set_full_eval!(mcmc::MCMCXMLElement, n_eval::Int)
    mcmc.attrs[bn.FULL_EVALUATION] = string(n_eval)
end

function make_xml(mc_el::MCMCXMLElement)
    el = new_element(bn.MCMC)
    set_attribute(el, bn.ID, bn.MCMC)

    mc_el.attrs[bn.CHAINLENGTH] = string(mc_el.chain_length)
    set_attributes(el, mc_el.attrs)

    posterior_el = new_child(el, bn.POSTERIOR)
    set_attribute(posterior_el, bn.ID, bn.POSTERIOR)
    prior_el = new_child(posterior_el, bn.PRIOR)
    set_attribute(prior_el, bn.ID, bn.PRIOR)

    for prior in mc_el.priors
        add_ref_els(prior_el, get_priors(prior))
    end
    like_el = make_xml(mc_el.likelihoods)
    add_child(posterior_el, like_el)

    add_ref_el(el, mc_el.operators.el)

    screen_log_el = new_child(el, bn.LOG)
    attrs = [(bn.ID, bn.SCREEN_LOG_ID), (bn.LOGEVERY, string(mc_el.screen_logEvery))]
    set_attributes(screen_log_el, attrs)

    default_logs = [posterior_el, prior_el, like_el]
    add_screenLog_cols(screen_log_el, default_logs)

    if mc_el.log_files

        file_log_el = new_child(el, bn.LOG)
        attrs = [(bn.ID, bn.FILE_LOG_ID), (bn.LOGEVERY, string(mc_el.file_logEvery)),
                (bn.FILENAME, "$(mc_el.filename).log"), (bn.OVERWRITE, bn.TRUE)]
        set_attributes(file_log_el, attrs)

        file_logs = Vector{XMLElement}(undef, 0)
        for myxml in mc_el.loggables.els
            file_logs = [file_logs; get_loggables(myxml)]
        end
        # file_logs = [get_loggable(loggable_el) for loggable_el in mc_el.loggables.els]

        final_file_logs = [default_logs; file_logs]
        for log in final_file_logs
            add_ref_el(file_log_el, log)
        end

    end

    for ext_el in mc_el.extras
        make_xml(ext_el)
        add_child(el, ext_el.el)
    end

    mc_el.el = el
    return el
end

function add_screenLog_cols(pel::XMLElement, cols::Array{XMLElement};
            width::Int = bn.WIDTH_DEFAULT, dp::Int = bn.DP_DEFAULT)
    for el in cols
        label = attribute(el, bn.ID, required = true)
        col_el = new_element(bn.COLUMN)
        attrs = [(bn.LABEL, label), (bn.DP, string(dp)), (bn.WIDTH, string(width))]
        set_attributes(col_el, attrs)
        add_ref_el(col_el, el)
        add_child(pel, col_el)
    end
end

function set_screen_logEvery(mcmc::MCMCXMLElement, sle::Int)
    mcmc.screen_logEvery = sle
end

function set_file_logEvery(mcmc::MCMCXMLElement, fle::Int)
    mcmc.file_logEvery = fle
end

function set_filename(mcmc::MCMCXMLElement, filename::String)
    mcmc.filename = filename
end

function merge_mcmc!(mcmc::MCMCXMLElement, pmcmc::ParsedMCMCXMLElement;
                     separate_logs::Bool = false)

    trait_id = "traitLikelihoods"
    trait_likelihood = mcmc.likelihoods
    trait_likelihood.id = trait_id

    mcmc.likelihoods = CompoundLikelihood("likelihood",
                            [mcmc.likelihoods; pmcmc.likelihoods],
                            already_made = [false; fill(true, length(pmcmc.likelihoods))])
    add_loggable(mcmc.loggables, trait_likelihood, already_made = true)

    mcmc.priors = [mcmc.priors; pmcmc.priors]

    if separate_logs
        push!(mcmc.extras, pmcmc.file_log)
    else
        push!(mcmc.loggables, LoggablesXMLElement(pmcmc.file_loggables))
    end

    push!(mcmc.extras, pmcmc.tree_log)
end

function add_prior!(mcmc::MCMCXMLElement, xml::MyXMLElement)
    push!(mcmc.priors, xml)
end

function add_likelihood!(mcmc::MCMCXMLElement, xml::MyXMLElement)
    add_likelihood!(mcmc.likelihoods, xml)
end

function add_loggable!(mcmc::MCMCXMLElement, xml::MyXMLElement)
    add_loggable(mcmc.loggables, xml)
end

function make_xml(cl::CompoundLikelihood)
    el = new_element(bn.LIKELIHOOD)
    set_attribute(el, bn.ID, cl.id)
    for i = 1:length(cl.likelihoods)
        if cl.already_made[i]
            add_ref_el(el, cl.likelihoods[i])
        else
            c_el = make_xml(cl.likelihoods[i])
            add_child(el, c_el)
        end
    end
    cl.el = el
    return el
end


text:
# Auxiliary functions for TimeModeling module
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: September 2016
#

export ricker_wavelet, get_computational_nt, smooth10, damp_boundary, calculate_dt, setup_grid, setup_3D_grid
export convertToCell, limit_model_to_receiver_area, extend_gradient, remove_out_of_bounds_receivers
export time_resample, remove_padding, backtracking_linesearch, subsample
export generate_distribution, select_frequencies, process_physical_parameter
export load_pymodel, load_acoustic_codegen, load_numpy

function limit_model_to_receiver_area(srcGeometry::Geometry, recGeometry::Geometry, model::Model, buffer; pert=[])
    # Restrict full velocity model to area that contains either sources and receivers
    ndim = length(model.n)
    println("N orig: ", model.n)

    # scan for minimum and maximum x and y source/receiver coordinates
    min_x = minimum([vec(recGeometry.xloc[1]); vec(srcGeometry.xloc[1])])
    max_x = maximum([vec(recGeometry.xloc[1]); vec(srcGeometry.xloc[1])])
    if ndim == 3
        min_y = minimum([vec(recGeometry.yloc[1]); vec(srcGeometry.yloc[1])])
        max_y = maximum([vec(recGeometry.yloc[1]); vec(srcGeometry.yloc[1])])
    end

    # add buffer zone if possible
    min_x = max(model.o[1], min_x-buffer)
    max_x = min(model.o[1] + model.d[1]*(model.n[1]-1), max_x+buffer)
    if ndim == 3
        min_y = max(model.o[2], min_y-buffer)
        max_y = min(model.o[2] + model.d[2]*(model.n[2]-1), max_y+buffer)
    end

    # extract part of the model that contains sources/receivers
    nx_min = Int(round(min_x/model.d[1])) + 1
    nx_max = Int(round(max_x/model.d[1])) + 1
    if ndim == 2
        ox = (nx_min - 1)*model.d[1]
        oz = model.o[2]
    else
        ny_min = Int(round(min_y/model.d[2])) + 1
        ny_max = Int(round(max_y/model.d[2])) + 1
        ox = (nx_min - 1)*model.d[1]
        oy = (ny_min - 1)*model.d[2]
        oz = model.o[3]
    end

    # Extract relevant model part from full domain
    n_orig = model.n
    if ndim == 2
        model.m = model.m[nx_min: nx_max, :]
        typeof(model.rho) <: Array && (model.rho = model.rho[nx_min: nx_max, :])
        model.o = (ox, oz)
    else
        model.m = model.m[nx_min:nx_max,ny_min:ny_max,:]
        typeof(model.rho) <: Array && (model.rho = model.rho[nx_min:nx_max,ny_min:ny_max,:])
        model.o = (ox,oy,oz)
    end
    model.n = size(model.m)
    println("N new: ", model.n)
    if isempty(pert)
        return model
    else
        if ndim==2
            pert = reshape(pert,n_orig)[nx_min: nx_max, :]
        else
            pert = reshape(pert,n_orig)[nx_min: nx_max,ny_min: ny_max, :]
        end
        return model,vec(pert)
    end
end

function extend_gradient(model_full::Model,model::Model,gradient::Array)
    # Extend gradient back to full model size
    ndim = length(model.n)
    full_gradient = zeros(Float32,model_full.n)
    nx_start = Int((model.o[1] - model_full.o[1])/model.d[1] + 1)
    nx_end = nx_start + model.n[1] - 1
    if ndim == 2
        full_gradient[nx_start:nx_end,:] = gradient
    else
        ny_start = Int((model.o[2] - model_full.o[2])/model.d[2] + 1)
        ny_end = ny_start + model.n[2] - 1
        full_gradient[nx_start:nx_end,ny_start:ny_end,:] = gradient
    end
    return full_gradient
end

function remove_out_of_bounds_receivers(recGeometry::Geometry, model::Model)

    # Only keep receivers within the model
    xmin = model.o[1]
    if typeof(recGeometry.xloc[1]) <: Array
        idx_xrec = findall(x -> x >= xmin, recGeometry.xloc[1])
        recGeometry.xloc[1] = recGeometry.xloc[1][idx_xrec]
        recGeometry.zloc[1] = recGeometry.zloc[1][idx_xrec]
    end

    # For 3D shot records, scan also y-receivers
    if length(model.n) == 3 && typeof(recGeometry.yloc[1]) <: Array
        ymin = model.o[2]
        idx_yrec = findall(x -> x >= ymin, recGeometry.yloc[1])
        recGeometry.yloc[1] = recGeometry.yloc[1][idx_yrec]
        recGeometry.zloc[1] = recGeometry.zloc[1][idx_yrec]
    end
    return recGeometry
end

function remove_out_of_bounds_receivers(recGeometry::Geometry, recData::Array, model::Model)

    # Only keep receivers within the model
    xmin = model.o[1]
    if typeof(recGeometry.xloc[1]) <: Array
        idx_xrec = findall(x -> x > xmin, recGeometry.xloc[1])
        recGeometry.xloc[1] = recGeometry.xloc[1][idx_xrec]
        recGeometry.zloc[1] = recGeometry.zloc[1][idx_xrec]
        recData[1] = recData[1][:, idx_xrec]
    end

    # For 3D shot records, scan also y-receivers
    if length(model.n) == 3 && typeof(recGeometry.yloc[1]) <: Array
        ymin = model.o[2]
        idx_yrec = findall(x -> x > ymin, recGeometry.yloc[1])
        recGeometry.yloc[1] = recGeometry.yloc[1][idx_yrec]
        recGeometry.zloc[1] = recGeometry.zloc[1][idx_yrec]
        recData[1] = recData[1][:, idx_yrec]
    end
    return recGeometry, recData
end

"""
    convertToCell(x)

Convert an array `x` to a cell array (`Array{Any,1}`) with `length(x)` entries,\\
where the i-th cell contains the i-th entry of `x`.

"""
function convertToCell(x)
    n = length(x)
    y = Array{Any}(undef, n)
    for j=1:n
        y[j] = x[j]
    end
    return y
end

# 1D source time function
"""
    source(tmax, dt, f0)

Create seismic Ricker wavelet of length `tmax` (in milliseconds) with sampling interval `dt` (in milliseonds)\\
and central frequency `f0` (in kHz).

"""
function ricker_wavelet(tmax, dt, f0)
    t0 = 0.
    nt = Int(trunc((tmax - t0)/dt + 1))
    t = range(t0,stop=tmax,length=nt)
    r = (pi * f0 * (t .- 1 / f0))
    q = zeros(Float32,nt,1)
    q[:,1] = (1f0 .- 2f0 .* r.^2f0) .* exp.(-r.^2f0)
    return q
end

function load_pymodel()
    pushfirst!(PyVector(pyimport("sys")["path"]), joinpath(JUDIPATH, "Python"))
    return pyimport("PyModel")
end

function load_numpy()
    pushfirst!(PyVector(pyimport("sys")["path"]), joinpath(JUDIPATH, "Python"))
    return pyimport("numpy")
end

function load_acoustic_codegen()
    pushfirst!(PyVector(pyimport("sys")["path"]), joinpath(JUDIPATH, "Python"))
    return pyimport("JAcoustic_codegen")
end

function calculate_dt(n,d,o,v,rho; epsilon=0)
    pm = load_pymodel()
    length(n) == 2 ? pyDim = [n[2], n[1]] : pyDim = [n[3],n[2],n[1]]
    modelPy = pm["Model"](o, d, pyDim, PyReverseDims(v))
    dtComp = modelPy[:critical_dt]
end

"""
    get_computational_nt(srcGeometry, recGeoemtry, model)

Estimate the number of computational time steps. Required for calculating the dimensions\\
of the matrix-free linear modeling operators. `srcGeometry` and `recGeometry` are source\\
and receiver geometries of type `Geometry` and `model` is the model structure of type \\
`Model`.

"""
function get_computational_nt(srcGeometry, recGeometry, model::Model)
    # Determine number of computational time steps
    if typeof(srcGeometry) == GeometryOOC
        nsrc = length(srcGeometry.container)
    else
        nsrc = length(srcGeometry.xloc)
    end
    nt = Array{Any}(undef, nsrc)
    dtComp = calculate_dt(model.n, model.d, model.o, sqrt.(1f0 ./ model.m), model.rho)
    for j=1:nsrc
        ntRec = recGeometry.dt[j]*(recGeometry.nt[j]-1) / dtComp
        ntSrc = srcGeometry.dt[j]*(srcGeometry.nt[j]-1) / dtComp
        nt[j] = max(Int(ceil(ntRec)), Int(ceil(ntSrc)))
    end
    return nt
end

function setup_grid(geometry,n, origin)
    # 3D grid
    if length(n)==3
        if length(geometry.xloc[1]) > 1
            source_coords = Array{Float32,2}([vec(geometry.xloc[1]) vec(geometry.yloc[1]) vec(geometry.zloc[1])])
        else
            source_coords = Array{Float32,2}([geometry.xloc[1] geometry.yloc[1] geometry.zloc[1]])
        end
        orig = Array{Float32}([origin[1] origin[2] origin[3]])
    else
    # 2D grid
        if length(geometry.xloc[1]) > 1
            source_coords = Array{Float32,2}([vec(geometry.xloc[1]) vec(geometry.zloc[1])])
        else
            source_coords = Array{Float32,2}([geometry.xloc[1] geometry.zloc[1]])
        end
        orig = Array{Float32}([origin[1] origin[2]])
    end
    return source_coords .- orig
end

function setup_3D_grid(xrec::Array{Any,1},yrec::Array{Any,1},zrec::Array{Any,1})
    # Take input 1d x and y coordinate vectors and generate 3d grid. Input are cell arrays
    nsrc = length(xrec)
    xloc = Array{Any}(undef, nsrc)
    yloc = Array{Any}(unfef, nsrc)
    zloc = Array{Any}(undef, nsrc)
    for i=1:nsrc
        nxrec = length(xrec[i])
        nyrec = length(yrec[i])

        xloc[i] = zeros(nxrec*nyrec)
        yloc[i] = zeros(nxrec*nyrec)
        zloc[i] = zeros(nxrec*nyrec)

        idx = 1

        for k=1:nyrec
            for j=1:nxrec
                xloc[i][idx] = xrec[i][j]
                yloc[i][idx] = yrec[i][k]
                zloc[i][idx] = zrec[i]
                idx += 1
            end
        end
    end
    return xloc, yloc, zloc
end

function setup_3D_grid(xrec,yrec,zrec)
# Take input 1d x and y coordinate vectors and generate 3d grid. Input are arrays/ranges
    nxrec = length(xrec)
    nyrec = length(yrec)

    xloc = zeros(nxrec*nyrec)
    yloc = zeros(nxrec*nyrec)
    zloc = zeros(nxrec*nyrec)
    idx = 1
    for k=1:nyrec
        for j=1:nxrec
            xloc[idx] = xrec[j]
            yloc[idx] = yrec[k]
            zloc[idx] = zrec
            idx += 1
        end
    end
    return xloc, yloc, zloc
end

function smooth10(velocity,shape)
    # 10 point smoothing function
    out = ones(Float32,shape)
    nz = shape[end]
    if length(shape)==3
        out[:,:,:] = velocity[:,:,:]
        for a=5:nz-6
            out[:,:,a] = sum(velocity[:,:,a-4:a+5], dims=3) / 10
        end
    else
        out[:,:] = velocity[:,:]
        for a=5:nz-6
            out[:,a] = sum(velocity[:,a-4:a+5], dims=2) / 10
        end
    end
    return out
end

function remove_padding(gradient::Array, nb::Integer; true_adjoint::Bool=false)
    if ndims(gradient) == 2
        if true_adjoint
            gradient[nb+1,:] = sum(gradient[1:nb,:], dims=1)
            gradient[end-nb,:] = sum(gradient[end-nb+1:end,:], dims=1)
            gradient[:,nb+1] = sum(gradient[:,1:nb], dims=2)
            gradient[:,end-nb] = sum(gradient[:,end-nb+1:end], dims=2)
        end
        return gradient[nb+1:end-nb,nb+1:end-nb]
    elseif ndims(gradient)==3
        if true_adjoint
            gradient[nb+1,:,:] = sum(gradient[1:nb,:,:], dims=1)
            gradient[end-nb,:,:] = sum(gradient[end-nb+1:end,:,:], dims=1)
            gradient[:,nb+1,:] = sum(gradient[:,1:nb,:], dims=2)
            gradient[:,end-nb,:] = sum(gradient[:,end-nb+1:end,:], dims=2)
            gradient[:,:,nb+1] = sum(gradient[:,:,1:nb], dims=3)
            gradient[:,:,end-nb] = sum(gradient[:,:,end-nb+1:end], dims=3)
        end
        return gradient[nb+1:end-nb,nb+1:end-nb,nb+1:end-nb]
    else
        error("Gradient must have 2 or 3 dimensions")
    end
end

# Vectorization of single variable (not defined in Julia)
vec(x::Float64) = x;
vec(x::Float32) = x;
vec(x::Int64) = x;
vec(x::Int32) = x;

function time_resample(data::Array,geometry_in::Geometry,dt_new;order=2)

    if dt_new==geometry_in.dt[1]
        return data, geometry_in
    else
        geometry = deepcopy(geometry_in)
        numTraces = size(data,2)
        timeAxis = 0:geometry.dt[1]:geometry.t[1]
        timeInterp = 0:dt_new:geometry.t[1]
        dataInterp = zeros(Float32,length(timeInterp),numTraces)
        for k=1:numTraces
            spl = Spline1D(timeAxis,data[:,k];k=order)
            dataInterp[:,k] = spl(timeInterp)
        end
        geometry.dt[1] = dt_new
        geometry.nt[1] = length(timeInterp)
        geometry.t[1] = (geometry.nt[1] - 1)*geometry.dt[1]
        return dataInterp, geometry
    end
end

function time_resample(data::Array,dt_in, geometry_out::Geometry;order=2)

    if dt_in==geometry_out.dt[1]
        return data
    else
        geometry = deepcopy(geometry_out)
        numTraces = size(data,2)
        timeAxis = 0:dt_in:geometry_out.t[1]
        timeInterp = 0:geometry_out.dt[1]:geometry_out.t[1]
        dataInterp = zeros(Float32,length(timeInterp),numTraces)
        for k=1:numTraces
            spl = Spline1D(timeAxis,data[:,k];k=order)
            dataInterp[:,k] = spl(timeInterp)
        end
        return dataInterp
    end
end

#subsample(x::Nothing) = x

function generate_distribution(x; src_no=1)
	# Generate interpolator to sample from probability distribution given
	# from spectrum of the input data

	# sampling information
	nt = x.geometry.nt[src_no]
	dt = x.geometry.dt[src_no]
	t = x.geometry.t[src_no]

	# frequencies
	fs = 1/dt	# sampling rate
	fnyq = fs/2	# nyquist frequency
	df = fnyq/nt	# frequency interval
	f = 0:2*df:fnyq	# frequencies

	# amplitude spectrum of data (serves as probability density function)
	ns = convert(Integer,ceil(nt/2))
	amp = abs.(fft(x.data[src_no]))[1:ns]	# get first half of spectrum

	# convert to cumulative probability distribution (integrate)
	pd = zeros(ns)
	pd[1] = dt*amp[1]
	for j=2:ns
		pd[j] = pd[j-1] + amp[j]*df
	end
	pd /= pd[end]	# normalize

	return Spline1D(pd, f)
end

function select_frequencies(q_dist; fmin=0f0, fmax=Inf, nf=1)
	freq = zeros(Float32, nf)
	for j=1:nf
		while (freq[j] <= fmin) || (freq[j] > fmax)
			freq[j] = q_dist(rand(1)[1])[1]
		end
	end
	return freq
end

function process_physical_parameter(param, dims)
    if length(param) ==1
        return param
    else
        return PyReverseDims(permutedims(param, dims))
    end
end