text:
module Atmosphere

using StaticArrays, StructArrays, ComponentArrays

using Flight.Geodesy

import Flight.Modeling: SystemDescriptor, System, init_x0, init_y0, init_u0, init_d0, f_cont!, f_disc!

export TunableISA
export TunableWind
export AtmosphereDescriptor, AtmosphericData, AtmosphericSystem

const R = 287.05287 #gas constant for dry air
const γ = 1.40 #heat capacity ratio for dry air
const βs = 1.458e-6 #Sutherland's empirical constant for dynamic viscosity
const S = 110.4 #Sutherland's empirical constant for dynamic viscosity

const T_std = 288.15
const p_std = 101325.0
const ρ_std = p_std / (R * T_std)
const g_std = 9.80665

const ISA_layers = StructArray(
    β =      SVector{7,Float64}([-6.5e-3, 0, 1e-3, 2.8e-3, 0, -2.8e-3, -2e-3]),
    h_ceil = SVector{7,Float64}([11000, 20000, 32000, 47000, 51000, 71000, 84852]))

@inline density(p,T) = p/(R*T)
@inline speed_of_sound(T) = √(γ*R*T)
@inline dynamic_viscosity(T) = (βs * T^1.5) / (T + S)

@inline ISA_temperature_law(h::Real, T_b, h_b, β)::Float64 = T_b + β * (h - h_b)

@inline function ISA_pressure_law(h::Real, g0, p_b, T_b, h_b, β)::Float64
    if β != 0.0
        p_b * (1 + β/T_b * (h - h_b)) ^ (-g0/(β*R))
    else
        p_b * exp(-g0/(R*T_b) * (h - h_b))
    end
end

Base.@kwdef struct SLConditions
    p::Float64 = p_std
    T::Float64 = T_std
    g::Float64 = g_std
end

struct ISAData
    p::Float64
    T::Float64
    ρ::Float64
    a::Float64
    μ::Float64
end

ISAData() = ISAData(AltGeop(0))

@inline function ISAData(h_geo::AltGeop, sl::SLConditions = SLConditions())

    h = Float64(h_geo)
    h_base = 0; T_base = sl.T; p_base = sl.p; g_base = sl.g

    for i in 1:length(ISA_layers)
        β, h_ceil = ISA_layers[i]
        if h < h_ceil
            T = ISA_temperature_law(h, T_base, h_base, β)
            p = ISA_pressure_law(h, g_base, p_base, T_base, h_base, β)
            return ISAData(p, T, density(p, T), speed_of_sound(T), dynamic_viscosity(T) )
        end
        T_ceil = ISA_temperature_law(h_ceil, T_base, h_base, β)
        p_ceil = ISA_pressure_law(h_ceil, g_base, p_base, T_base, h_base, β)
        h_base = h_ceil; T_base = T_ceil; p_base = p_ceil
    end

    throw(ArgumentError("Altitude out of bounds"))

end

##################### AbstractISA ###########################

#a System{<:AbstractISA} may have an output type of its own, as any other
#System. however, this output generally will not be an ISAData instance. the
#ISAData returned by a ISA System depends on the location specified in the
#query. instead, the output from an ISA System may hold quantities of interest
#related to its own internal state, if it has one due to it being a dynamic ISA
#implementation.

#AbstractISA subtypes: ConstantUniformISA (TunableISA), ConstantFieldISA,
#DynamicUniformISA, DynamicFieldISA.

abstract type AbstractISA <: SystemDescriptor end

function SLConditions(::T, ::Abstract2DLocation) where {T<:System{<:AbstractISA}}
    error("SLConditions constructor not implemented for $T")
end

@inline function ISAData(sys::System{<:AbstractISA}, p::Geographic)

    h_geop = Altitude{Geopotential}(p.alt, p.l2d)
    sl = SLConditions(sys, p.l2d)
    ISAData(h_geop, sl)

end

##################### TunableISA ###########################

#a TunableISA System does not have a state and therefore cannot evolve on its
#own. but its input vector can still be used to manually tune the SL conditions
#during simulation.

struct TunableISA <: AbstractISA end #Constant, Uniform ISA

Base.@kwdef mutable struct UTunableISA #only allocates upon System instantiation
    T_sl::Float64 = T_std
    p_sl::Float64 = p_std
end

init_u0(::TunableISA) = UTunableISA()
f_cont!(::System{<:TunableISA}, args...) = nothing
f_disc!(::System{<:TunableISA}, args...) = false

function SLConditions(s::System{<:TunableISA}, ::Abstract2DLocation)
    SLConditions(T = s.u.T_sl, p = s.u.p_sl, g = g_std)
    #alternative using actual local SL gravity:
    # return (T = s.u.T_sl, p = s.u.p_sl, g = gravity(Geographic(l2d, AltOrth(0.0))))
end

######################## AbstractWind ############################

abstract type AbstractWind <: SystemDescriptor end

Base.@kwdef struct WindData
    v_ew_n::SVector{3,Float64} = zeros(SVector{3})
end

function WindData(::T, ::Abstract3DLocation) where {T<:System{<:AbstractWind}}
    error("WindData constructor not implemented for $T")
end

######################## TunableWind ############################

struct TunableWind <: AbstractWind end

Base.@kwdef mutable struct USimpleWind
    v_ew_n::MVector{3,Float64} = zeros(MVector{3}) #MVector allows changing single components
end

init_u0(::TunableWind) = USimpleWind()
f_cont!(::System{<:TunableWind}, args...) = nothing
f_disc!(::System{<:TunableWind}, args...) = false

function WindData(wind::System{<:TunableWind}, ::Abstract3DLocation)
    wind.u.v_ew_n |> SVector{3,Float64} |> WindData
end

#################### AtmosphereDescriptor ############################

Base.@kwdef struct AtmosphereDescriptor{I <: AbstractISA, W <: AbstractWind} <: SystemDescriptor
    isa_::I = TunableISA()
    wind::W = TunableWind()
end

#this one's convenient for dispatching on plots() methods, but for U and D we
#can use simply NamedTuples
struct AtmosphereCmpY{I, W}; isa_::I; wind::W; end

init_x0(atm::AtmosphereDescriptor) = ComponentVector(isa_ = init_x0(atm.isa_), wind = init_x0(atm.wind))
init_y0(atm::AtmosphereDescriptor) = AtmosphereCmpY(init_y0(atm.isa_), init_y0(atm.wind))
init_u0(atm::AtmosphereDescriptor) = (isa_ = init_u0(atm.isa_), wind = init_u0(atm.wind))
init_d0(atm::AtmosphereDescriptor) = (isa_ = init_d0(atm.isa_), wind = init_d0(atm.wind))

function System(atm::AtmosphereDescriptor, ẋ = init_x0(atm), x = init_x0(atm),
                    y = init_y0(atm), u = init_u0(atm), d = init_d0(atm), t = Ref(0.0))

    isa_sys = System(atm.isa_, ẋ.isa_, x.isa_, y.isa_, u.isa_, d.isa_, t)
    wind_sys = System(atm.wind, ẋ.wind, x.wind, y.wind, u.wind, d.wind, t)
    params = nothing
    subsystems = (isa_ = isa_sys, wind = wind_sys,)
    System{map(typeof, (atm, x, y, u, d, params, subsystems))...}(
                                ẋ, x, y, u, d, t, params, subsystems)
end

const AtmosphericSystem = System{<:AtmosphereDescriptor}

Base.@kwdef struct AtmosphericData
    isa_::ISAData = ISAData()
    wind::WindData = WindData()
end

function AtmosphericData(a::AtmosphericSystem, pos::Geographic)
    AtmosphericData(
        ISAData(a.subsystems.isa_, pos),
        WindData(a.subsystems.wind, pos))
end

function f_cont!(sys::AtmosphericSystem)
    f_cont!(sys.isa_)
    f_cont!(sys.wind)
end

#we can create a World <: SystemDescriptor to compose Aircraft and
#AtmosphericSystem, and step them together in a single Model. if the
#AtmosphericSystem is dynamic but still we don't want to compose it with the
#Aircraft, we can always put it in a separate model, step it separately, and
#simply pass a reference to its underlying AtmosphericSystem to the
#AircraftSystem, which will be wrapped in its own Model

# #top-down / recursive implementation
# @inline function get_tp(h::Real, T0::Real = T0_std, p0::Real = p0_std, g0::Real = g0_std)

#     h == 0 && return (T0, p0)
#     (h_b, β) = layer_parameters(h)
#     (T_b, p_b) = get_tp(h_b, T0, p0, g0) #get pt at the layer base
#     T = ISA_temperature_law(h, T_b, h_b, β)
#     p = ISA_pressure_law(h, g0, p_b, T_b, h_b, β)
#     return (T, p)
# end

end


text:
""" Fixes the unlabeled leaves of an MIOTree.  """
function fix_labeling!(mt::MIOTree)
    queue = allleaves(mt)
    while !isempty(queue)
        leaf = popfirst!(queue)
        if isnothing(leaf.label) && !isnothing(leaf.parent)
            par = leaf.parent
            labels = [child.label for child in children(leaf.parent) if !isnothing(child.label)]
            if length(labels) == 1
                set_classification_label!(par, labels[1])
                delete_children!(par)
                if !isnothing(par.parent)
                    append!(queue, children(par.parent))
                end
            end
        end
    end
    return
end

ext: jl


text:
using DiffEqPy
using Base.Test

# The test is just that it builds.
# Most of this is tested from the Python side!


text:
const results = ["Balanced","Not Balanced"]

function balancednum(number)
  if number < 100
    # one and two digits are balanced
    return results[1]
  end
  d = digits(number)
  m = (length(d)-1) ÷ 2
  n = (length(d)+4) ÷ 2
  return results[sum(d[1:m]) == sum(d[n:end]) ? 1 : 2]
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using FixedSizeArrays
using FactCheck

immutable Normal{N, T} <: FixedVector{N, T}
    _::NTuple{N, T}
end
immutable RGB{T} <: FixedVectorNoTuple{3, T}
    r::T
    g::T
    b::T
end

typealias Vec1d Vec{1, Float64}
typealias Vec2d Vec{2, Float64}
typealias Vec3d Vec{3, Float64}
typealias Vec4d Vec{4, Float64}
typealias Vec3f Vec{3, Float32}
typealias Mat2d Mat{2,2, Float64}
typealias Mat3d Mat{3,3, Float64}
typealias Mat4d Mat{4,4, Float64}

function test()
a = 1
a1 = @fsa([a,2,3])
a2 = @fsa([a 2 3])
a3 = @fsa([a;2;3])
a4 = @fsa([a 2;3 4])
a5 = @fsa([a 2 3;4 5 6])
a6 = @fsa([a 2;3 4;5 6])

a1 == Vec(a,2,3)
a2 == Mat((a,),(2,),(3,))
a3 == Mat(((a,2,3),))
a4 == Mat(((a,3),(2,4)))
a5 == Mat(((a,4),(2,5),(3,6)))
a6 == Mat(((a,3,5),(2,4,6)))

N = 100
a = Point{3, Float32}[Point{3, Float32}(0.7132) for i=1:N]
b = RGB{Float32}[RGB{Float32}(52.293) for i=1:N]

c = Point{3, Float64}[Point{3, Float64}(typemin(Float64)), a..., Point{3, Float64}(typemax(Float64))]
d = RGB{Float64}[RGB(typemin(Float64)), b..., RGB(typemax(Float64))]

sa = sum(a)
ma = mean(a)
sb = sum(b)
mb = mean(b)

maximum(c) == Point{3, Float32}(typemax(Float64))
minimum(c) == Point{3, Float32}(typemin(Float64))

maximum(d) == RGB(typemax(Float64))
minimum(d) == RGB(typemin(Float64))
af = a + 1f0
bf = b + 1f0
for elem in af
    a[1] + 1f0 == elem
end
for elem in bf
    b[1] + 1f0 == elem
end

for T=[Float32, Float64, Int, UInt, UInt32, UInt8]
    r = rand(T)
    x = RGB{Int}[RGB(1) for i=1:10]
    RGB{Float32}(["0.222", "9.8822", "29.999"]) == RGB{Float32}(0.222, 9.8822, 29.999)
    typeof(map(RGB{Float32}, x))  == Vector{RGB{Float32}}
    RGB{T}(r)                     == RGB(r,r,r)
    RGB{T}([r,r,r])               == RGB(r,r,r)
    length(RGB{T}([r,r,r]))       == 3
    length(RGB{T})                == 3
    eltype(RGB{T}([r,r,r]))       == T
    eltype(RGB{T})                == T
    typeof(RGB(r,r,r))            == RGB{T}
    typeof(RGB{T}(1))             == RGB{T}
    typeof(RGB{T}(1,2,3))         == RGB{T}
    ndims(RGB{T}(1,2,3))          == 1

    typeof(RGB{T}(1f0))           == RGB{T}
    typeof(RGB{T}(1f0,2f0,3f0))   == RGB{T}
    typeof(RGB{T}(1f0, 2, 3.0))   == RGB{T}
    typeof(RGB(1f0, 2, 3.0))      == RGB{Float64}
    typeof(RGB{Int}(1f0, 2, 3.0)) == RGB{Int}
end

# A little brutal, but hey.... Better redudantant tests, than not enough tests
for N=1:3:10
    for VT=[Point, Vec, Normal], VT2=[Normal, Vec, Point], ET=[Float32, Int, UInt64, Float64], ET2=[Float64, UInt64, Int, Float32]
        rand_range  = ET(1):ET(10)
        rand_range2 = ET2(1):ET2(10)
        rn = rand(rand_range, N)
        v0 = VT(rn)
        # parse constructor:
        VT{N, ET}(map(string, rn)) == v0
        # multi constructor
        v1 = VT{N, ET}(rn...)
        v1 == v0
        typeof(v1) == VT{N, ET}
        length(v1) == N
        eltype(v1) == ET
        ndims(v1) == 1

        length(typeof(v1)) == N
        eltype(typeof(v1)) == ET

        for i=1:N
            v1[i] == rn[i]
        end
        # from other FSA without parameters
        v2 = VT2(v1)

        typeof(v2) == VT2{N, ET}
        length(v2) == N
        eltype(v2) == ET
        for i=1:N
            v2[i] == v1[i]
        end
        # from other FSA with parameters
        for i=1:N
            v3 = VT2{i, ET2}(v1)
            typeof(v3) == VT2{i, ET2}
            length(v3) == i
            eltype(v3) == ET2
            for i=1:i
                v3[i] == ET2(v2[i])
            end
        end
        # from single
        r  = rand(rand_range)
        r2 = rand(rand_range2)
        v1 = VT{N, ET}(r)
        v2 = VT{N, ET2}(r)
        v3 = VT{N, ET}(r2)
        v4 = VT{N, ET2}(r2)

        for i=1:N
            v1[i] == r
            v2[i] == ET2(r)
            v3[i] == r2
            v4[i] == ET2(r2)
        end
        x = VT{N, ET}[VT{N, ET}(1) for i=1:10]
        x1 = VT2{N, ET}[VT{N, ET}(1) for i=1:10]
        x2 = map(VT2, x)
        x3 = map(VT, x2)
        typeof(x)  == Vector{VT{N, ET}}
        typeof(x1) == Vector{VT2{N, ET}}
        typeof(x2) == Vector{VT2{N, ET}}
        typeof(x3) == Vector{VT{N, ET}}
        x3         == x
    end
end

typeof(Vec3f(1,1,1))     == Vec{3, Float32}
typeof(Vec3f(1,1f0,1))   == Vec{3, Float32}
typeof(Vec3f(1f0,1,1.0)) == Vec{3, Float32}

typeof(Vec3f(1))  	== Vec{3, Float32}
typeof(Vec3f(0))  	== Vec{3, Float32}
Vec3f(1.0) 			== Vec(1f0,1f0,1f0)
Vec3f(1.0f0) 			== Vec(1f0,1f0,1f0)
Vec3f(1.0f0) 			== Vec3f(1)
Vec(1.0, 1.0, 1.0) 	== Vec3d(1)
Vec2d(Vec3d(1)) 		== Vec(1.0, 1.0)
Vec(Vec3d(1), 1.0) 	== Vec4d(1)
Vec(Vec3d(1), 1) 		== Vec4d(1)
Vec3d(Vec3f(1.0)) 	== Vec3d(1.0)

v2 = Vec(6.0,5.0,4.0)
v1 = Vec(1.0,2.0,3.0)
v2 = Vec(6.0,5.0,4.0)

setindex(v1, 88.9, 1) == Vec(88.9,2.0,3.0)
v1[1] == 1.0
v1[2] == 2.0
v1[3] == 3.0
v1[1:3] == (1.0, 2.0, 3.0)
v1[1:2] == (1.0, 2.0)
v1[1:1] == (1.0,)
v1[(1,2)] == (1.0,2.0)
v1[(2,1)] == (2.0,1.0)
m = Mat{4,4,Int}(
(1,2,3,4),
(5,6,7,8),
(9,10,11,12),
(13,14,15,16)
)
setindex(m, 42.0, 2,2) == Mat{4,4,Int}(
(1,2,3,4),
(5,42.0,7,8),
(9,10,11,12),
(13,14,15,16)
)
m[1] == 1
m[2] == 2
m[10] == 10
m[2,2] == 6
m[3,4] == 15
m[1:4, 1] == (1,5,9,13)
m[1, 1:4] == (1,2,3,4)

-v1 == Vec(-1.0,-2.0,-3.0)
isa(-v1, Vec3d) == true
v1+v2 == Vec3d(7.0,7.0,7.0)
v2-v1 == Vec3d(5.0,3.0,1.0)
v1.*v2 == Vec3d(6.0,10.0,12.0)
v1 ./ v1 == Vec3d(1.0,1.0,1.0)
1.0 + v1 == Vec3d(2.0,3.0,4.0)
1.0 .+ v1 == Vec3d(2.0,3.0,4.0)
v1 + 1.0 == Vec3d(2.0,3.0,4.0)
v1 .+ 1.0 == Vec3d(2.0,3.0,4.0)
1 + v1 == Vec3d(2.0,3.0,4.0)
1 .+ v1 == Vec3d(2.0,3.0,4.0)
v1 + 1 == Vec3d(2.0,3.0,4.0)
v1 .+ 1 == Vec3d(2.0,3.0,4.0)

v1 - 1.0 == Vec3d(0.0,1.0,2.0)
v1 .- 1.0 == Vec3d(0.0,1.0,2.0)
1.0 - v1 == Vec3d(0.0,-1.0,-2.0)
1.0 .- v1 == Vec3d(0.0,-1.0,-2.0)
v1 - 1 == Vec3d(0.0,1.0,2.0)
v1 .- 1 == Vec3d(0.0,1.0,2.0)
1 - v1 == Vec3d(0.0,-1.0,-2.0)
1 .- v1 == Vec3d(0.0,-1.0,-2.0)

2.0 * v1 == Vec3d(2.0,4.0,6.0)
2.0 .* v1 == Vec3d(2.0,4.0,6.0)
v1 * 2.0 == Vec3d(2.0,4.0,6.0)
v1 .* 2.0 == Vec3d(2.0,4.0,6.0)
2 * v1 == Vec3d(2.0,4.0,6.0)
2 .* v1 == Vec3d(2.0,4.0,6.0)
v1 * 2 == Vec3d(2.0,4.0,6.0)
v1 .* 2 == Vec3d(2.0,4.0,6.0)

v1 / 2.0 == Vec3d(0.5,1.0,1.5)
v1 ./ 2.0 == Vec3d(0.5,1.0,1.5)
v1 / 2 == Vec3d(0.5,1.0,1.5)
v1 ./ 2 == Vec3d(0.5,1.0,1.5)

12.0 ./ v1 == Vec3d(12.0,6.0,4.0)
12 ./ v1 == Vec3d(12.0,6.0,4.0)

(v1 .^ 2) == Vec3d(1.0,4.0,9.0)
(v1 .^ 2.0) == Vec3d(1.0,4.0,9.0)
(2.0 .^ v1) == Vec3d(2.0,4.0,8.0)
(2 .^ v1) == Vec3d(2.0,4.0,8.0)
norm(Vec3d(1.0,2.0,2.0)) == 3.0

# cross product
cross(v1,v2) == Vec3d(-7.0,14.0,-7.0)
isa(cross(v1,v2),Vec3d)  == true
a = Vec{2,Int}(1,2)
b = Vec{2,Float64}(1.,2.)
hypot(a) == 2.23606797749979
hypot(b) == 2.23606797749979
hypot(a) == hypot(b) == true

# type conversion
isa(convert(Vec3f,v1), Vec3f)  == true

isa(convert(Vector{Float64}, v1), Vector{Float64})  == true
convert(Vector{Float64}, v1) == [1.0,2.0,3.0]

# matrix operations

#typealias Mat1d Matrix1x1{Float64}

zeromat = Mat2d((0.0,0.0),(0.0,0.0))

length(Mat2d) == 4
length(zeromat) == 4

size(Mat2d) == (2,2)
size(zeromat) == (2,2)

zero(Mat2d) == zeromat

for i=1:4, j=1:4
    x1 = rand(i,j)
    Mat(x1') == Mat(x1)'
end

v = Vec(1.0,2.0,3.0,4.0)
r = row(v)
c = column(v)

#r' == c
#c' == r

a = c*r

b = Mat(
	(1.0,2.0,3.0,4.0),
	(2.0,4.0,6.0,8.0),
	(3.0,6.0,9.0,12.0),
	(4.0,8.0,12.0,16.0)
)

length(b) == 16

a==b
mat30 = Mat(((30.0,),))
r*c == mat30

#row(r, 1) == v
#column(c,1) == v
#row(r+c',1) == 2*v
sum(r) == sum(v)
prod(c) == prod(v)

eye(Mat3d) == Mat((1.0,0.0,0.0),
	(0.0,1.0,0.0),
	(0.0,0.0,1.0))
#v*eye(Mat4d)*v == 30.0
x = -r
y = -1.0*r
jm = rand(4,4)
lowl = Mat(jm)
jm2 = convert(Array{Float64,2}, lowl)

#Single valued constructor
Mat4d(0.0)
a = Vec4d(0)
b = Vec4d(0,0,0,0)

v = rand(4)
m = rand(4,4)
vfs = Vec(v)
mfs = Mat(m)

for i=1:4, j=1:4
	vfs = rand(Vec{j, Float64})
	mfs = rand(Mat{i,j, Float64})
	vm = m * v
	fsvm = mfs * vfs
	if i == j
		fmm = mfs * mfs
		fmm = det(mfs)
		fmm = inv(mfs)
		fmm = expm(mfs)
    end
end
mm = m'
fmm = mfs'

ac = rand(3)
bc = rand(3)

a = rand(4)
b = rand(4)
c = rand(4,4)

d = cross(ac, bc)
d2 = a+b
f = c*a
g = c*b
h = c*f
i = dot(f, a)
j = dot(a, g)
k = abs(f)
l = abs(-f)

acfs = Vec(ac)
bcfs = Vec(bc)

afs = Vec(a)
bfs = Vec(b)
cfs = Mat(c)

dfs = cross(acfs, bcfs)
d2fs = afs+bfs
ffs = cfs*afs
gfs = cfs*bfs
hfs = cfs*ffs
ifs = dot(ffs, afs)
jfs = dot(afs, gfs)
kfs = abs(ffs)
lfs = abs(-ffs)

isapprox(acfs, ac)  == true
isapprox(bcfs, bc)  == true

isapprox(afs, a) == true
isapprox(bfs, b) == true
isapprox(cfs, c) == true

isapprox(dfs, d) == true
isapprox(d2fs, d2) == true
isapprox(ffs, f) == true
isapprox(gfs, g) == true
isapprox(hfs, h) == true
isapprox(ifs, i) == true
isapprox(jfs, j) == true
isapprox(kfs, k) == true
isapprox(lfs, l) == true

Vec{3, Int}(1) == Vec{3, Float64}(1)
Vec(1,2,3) == Vec(1.0,2.0,3.0)
Vec(1,2,3) == [1,2,3]
Mat((1,2),(3,4)) == Mat((1,2),(3,4))

const unaryOps = (
    -, ~, conj, abs,
    sin, cos, tan, sinh, cosh, tanh,
    asin, acos, atan, asinh, acosh, atanh,
    sec, csc, cot, asec, acsc, acot,
    sech, csch, coth, asech, acsch, acoth,
    sinc, cosc, cosd, cotd, cscd, secd,
    sind, tand, acosd, acotd, acscd, asecd,
    asind, atand, rad2deg, deg2rad,
    log, log2, log10, log1p, exponent, exp,
    exp2, expm1, cbrt, sqrt, erf,
    erfc, erfcx, erfi, dawson,

    #trunc, round, ceil, floor, #see JuliaLang/julia#12163
    significand, lgamma,
    gamma, lfact, frexp, modf, airy, airyai,
    airyprime, airyaiprime, airybi, airybiprime,
    besselj0, besselj1, bessely0, bessely1,
    eta, zeta, digamma
)

# vec-vec and vec-scalar
const binaryOps = (
    .+, .-,.*, ./, .\, /,
    .==, .!=, .<, .<=, .>, .>=, +, -,
    min, max,

    atan2, besselj, bessely, hankelh1, hankelh2,
    besseli, besselk, beta, lbeta
)

test1 = (Vec(1,2,typemax(Int)), Mat((typemin(Int),2,5), (2,3,5), (-2,3,6)), Vec{4, Float32}(0.777))
test2 = (Vec(1,0,typemax(Int)), Mat((typemin(Int),77,1), (2,typemax(Int),5), (-2,3,6)), Vec{4, Float32}(-23.2929))
for op in binaryOps
    for i=1:length(test1)
        v1 = test1[i]
        v2 = test2[i]
        try # really bad tests, but better than nothing...
            if applicable(op, v1[1], v2[1]) && typeof(op(v1[1], v2[1])) == eltype(v1)
                r = op(v1, v2)
                for j=1:length(v1)
                    r[j] == op(v1[j], v2[j])
                end

            end
        end
    end
end

test = (Vec(1,2,typemax(Int)), Mat((typemin(Int),2,5), (2,3,5), (-2,3,6)), Vec{4, Float32}(0.777))
for op in unaryOps
    for t in test
        try
            if applicable(op, t[1]) && typeof(op(t[1])) == eltype(t)
                v = op(t)
                for i=1:length(v)
                    v[i] == op(t[i])
                end
            end
        end
    end
end

end

using SnoopCompile

snoop_on()
SnoopCompile.@snoop "/tmp/fsa_compiles.csv" begin
   test()
end
snoop_off()

data = SnoopCompile.read("/tmp/fsa_compiles.csv")
pc, discards = SnoopCompile.parcel(data[end:-1:1,2])
SnoopCompile.write(pc, "/tmp/precompile")


text:
using Complementarity, JuMP
using Base.Test

#########################################################################
m = MCPModel()

@variable(m, x3 >= 0)
@variable(m, x4 >= 0)
@variable(m, x1 >= 0)
@variable(m, x2 >= 0)

@NLexpression(m, F2, x3-2x4 +2)
@NLexpression(m, F3, x1-x2+2x3-2x4 -2)
@NLexpression(m, F4, x1+2x2-2x3+4x4 -6)
@NLexpression(m, F1, -x3-x4 +2)

complements(m, F4, x4)
complements(m, F1, x1)
complements(m, F3, x3)
complements(m, F2, x2)

PATHSolver.path_options(
                "convergence_tolerance 1e-8",
                "output no",
                "time_limit 3600"
                )
status = solveMCP(m)

z = [getvalue(x1), getvalue(x2), getvalue(x3), getvalue(x4)]
@show z

@test isapprox(z, [2.8, 0.0, 0.8, 1.2])

#########################################################################

println("------------------------------------------------------------------")

#########################################################################
m = MCPModel()

@variable(m, x3 >= 0)
@variable(m, x4 >= 0)
@variable(m, x1 >= 0)
@variable(m, x2 >= 0)

@NLexpression(m, F2, x3-2x4 +2)
@NLexpression(m, F3, x1-x2+2x3-2x4 -2)
@NLexpression(m, F4, x1+2x2-2x3+4x4 -6)
@NLexpression(m, F1, -x3-x4 +2)

complements(m, F2, x2)
complements(m, F3, x3)
complements(m, F1, x1)
complements(m, F4, x4)

PATHSolver.path_options(
                "convergence_tolerance 1e-8",
                "output no",
                "time_limit 3600"
                )
status = solveMCP(m)

z = [getvalue(x1), getvalue(x2), getvalue(x3), getvalue(x4)]
@show z
@test isapprox(z, [2.8, 0.0, 0.8, 1.2])
#########################################################################

println("------------------------------------------------------------------")

#########################################################################
m = MCPModel()

M = [0  0 -1 -1 ;
     0  0  1 -2 ;
     1 -1  2 -2 ;
     1  2 -2  4 ]

q = [2; 2; -2; -6]

lb = zeros(4)
ub = Inf*ones(4)

@variable(m, lb[i] <= x[i in 1:4] <= ub[i])
@NLexpression(m, F[i=1:4], sum{M[i,j]*x[j], j=1:4} + q[i])
complements(m, F, x)

PATHSolver.path_options(
                "convergence_tolerance 1e-8",
                "output no",
                "time_limit 3600"
                )

status = solveMCP(m)

z = getvalue(x)
@test isapprox(z[1], 2.8)
@test isapprox(z[2], 0.0)
@test isapprox(z[3], 0.8)
@test isapprox(z[4], 1.2)
#########################################################################

println("------------------------------------------------------------------")

#########################################################################
m = MCPModel()

M = [0  0 -1 -1 ;
     0  0  1 -2 ;
     1 -1  2 -2 ;
     1  2 -2  4 ]

q = [2; 2; -2; -6]

lb = zeros(4)
ub = Inf*ones(4)

items = 1:4

@variable(m, lb[i] <= x[i in items] <= ub[i])
# @variable(m, x[i in items] >= 0)
@NLexpression(m, F[i in items], sum{M[i,j]*x[j], j in items} + q[i])
complements(m, F, x)

PATHSolver.path_options(
                "convergence_tolerance 1e-8",
                "output no",
                "time_limit 3600"
                )

status = solveMCP(m)

z = getvalue(x)
@test isapprox(z[1], 2.8)
@test isapprox(z[2], 0.0)
@test isapprox(z[3], 0.8)
@test isapprox(z[4], 1.2)
#########################################################################

println("------------------------------------------------------------------")

println("-------[Testing NLSolve]------------------------------------------")

#########################################################################
m = nothing
m = MCPModel()

M = [0  0 -1 -1 ;
     0  0  1 -2 ;
     1 -1  2 -2 ;
     1  2 -2  4 ]

q = [2; 2; -2; -6]

lb = zeros(4)
ub = Inf*ones(4)

items = 1:4

# @variable(m, lb[i] <= x[i in items] <= ub[i])
@variable(m, x[i in items] >= 0)
@NLexpression(m, F[i in items], sum{M[i,j]*x[j], j in items} + q[i])
complements(m, F, x)

status = solveMCP(m, solver=:NLsolve)
@show status

z = getvalue(x)
@show z
# @test isapprox(z[1], 2.8)
# @test isapprox(z[2], 0.0)
# @test isapprox(z[3], 0.8)
# @test isapprox(z[4], 1.2)
#########################################################################

println("------------------------------------------------------------------")

#########################################################################
m = nothing
m = MCPModel()

lb = zeros(4)
ub = Inf*ones(4)
items = 1:4
 @variable(m, lb[i] <= x[i in items] <= ub[i])

@NLexpression(m, F1, 3*x[1]^2+2*x[1]*x[2]+2*x[2]^2+x[3]+3*x[4]-6)
@NLexpression(m, F2, 2*x[1]^2+x[1]+x[2]^2+3*x[3]+2*x[4]-2)
@NLexpression(m, F3, 3*x[1]^2+x[1]*x[2]+2*x[2]^2+2*x[3]+3*x[4]-1)
@NLexpression(m, F4, x[1]^2+3*x[2]^2+2*x[3]+3*x[4]-3)

complements(m, F1, x[1])
complements(m, F2, x[2])
complements(m, F3, x[3])
complements(m, F4, x[4])

status = solveMCP(m, solver=:NLsolve, method=:trust_region)
@show status

z = getvalue(x)
@show z

# @test isapprox(z[1], 1.2247, atol=1e-4)
# @test isapprox(z[2], 0.0, atol=1e-4)
# @test isapprox(z[3], 0.0, atol=1e-4)
# @test isapprox(z[4], 0.5, atol=1e-4)

#########################################################################


text:
module JupyterSite

end


text:
using DistributedFactorGraphs

import DistributedFactorGraphs: compare

mutable struct SASDebug
  beams::Vector{Vector{Float64}}
  azi_smpls::Vector{Float64}
  SASDebug() = new()
end

mutable struct SASREUSE
    CBFFull::Array{Complex{Float64}}
    CBFLIE::Array{Complex{Float64}}
    BFOutFull::Array{Complex{Float64}}
    BFOutLIE::Array{Complex{Float64}}
    BFOutLIETemp::Array{Complex{Float64}}
    phaseshiftLOO::Array{Complex{Float64}}
    arrayPos::Array{Float64}
    arrayPosLIE::Array{Float64}
    hackazi::Tuple{Vector{Int}, Vector{Float64}} # [idx], [azi]
    oncebackidx::Vector{Int} # [idx]

    workvarlist::Vector{Symbol}
    looelement::Int
    elementset::Vector{Int}
    waveformsLIE::Array{Complex{Float64},2}
    waveformsLOOc::Array{Complex{Float64},1}

    BFtemp::Array{Complex{Float64}}

    LIEtemp::Array{Complex{Float64}}
    temp1::Array{Complex{Float64},1}
    temp2::Array{Complex{Float64},2}

    temp::Array{Complex{Float64}}
    sourceXY::Vector{Float64}

    dbg::SASDebug
    SASREUSE() = new()
end

mutable struct SASBearing2D <: AbstractRelativeMinimize
  cfgTotal::CBFFilterConfig
  cfgLIE::CBFFilterConfig
  waveformsIn::Array{Complex{Float64}}

  # to allow multithreaded behaviour
  threadreuse::Vector{SASREUSE}

  rangemodel::Union{Distributions.Rayleigh, Vector{IIF.AliasingScalarSampler}} # , Vector{<:StatsBase.AbstractWeights}
  ranget::Vector{Float64}
  cfg::Dict
  waveformsRaw::Array{Float64,2}
  debugging::Bool

  SASBearing2D() = new()
  # SASBearing2D(z1) = new(z1){Float64,2}
end

# untested manifold use
getManifold(::SASBearing2D) = BearingRange_Manifold

function getSample(cfo::CalcFactor{<:SASBearing2D}, N::Int=1)
  sas2d = cfo.factor
  # TODO bring the solvefor index tp getSample function
  if isa(sas2d.rangemodel, Distributions.Rayleigh)
    return (reshape(rand(sas2d.rangemodel, N),1,N), )
  else
    # TODO have to use the first waveform, since we do not know who is being solved for
    return (reshape(rand(sas2d.rangemodel[1], N), 1,N), )
  end
end

function forwardsas(sas2d::SASBearing2D,
                    thread_data::SASREUSE,
                    res::AbstractVector{<:Real},
                    idx::Int,
                    meas,
                    L1,
                    XX...)
  # dx, dy, azi = 0.0, 0.0, 0.0
  # thread_data = sas2d.threadreuse[Threads.threadid()]

  if thread_data.hackazi[1][1] != idx
    for i in 1:length(XX)
      thread_data.arrayPos[i,1:2] = XX[i][1:2]-XX[1][1:2]
    end
    constructCBFFilter2D!(sas2d.cfgTotal, thread_data.arrayPos, thread_data.CBFFull,thread_data.sourceXY, thread_data.BFtemp)

    #fill!(thread_data.temp1,0)
    #fill!(thread_data.temp2,0)

    CBF2D_DelaySum!(sas2d.cfgTotal, sas2d.waveformsIn, thread_data.BFOutFull,thread_data.temp1,thread_data.temp2,thread_data.CBFFull)

    thread_data.hackazi[1][1] = idx

    azi_weights = (repeat(norm.(thread_data.BFOutFull), 1,8)')[:]
    azi_domain = collect(LinRange(0,2pi-1e-10,length(azi_weights)))
    bss = AliasingScalarSampler(azi_domain, azi_weights, SNRfloor=sas2d.cfg["azi_snr_floor"])

    thread_data.hackazi[2][1] = wrapRad(rand(bss)[1])

    if sas2d.debugging
        push!(thread_data.dbg.beams, deepcopy(thread_data.BFOutFull))
        push!(thread_data.dbg.azi_smpls, thread_data.hackazi[2][1])
    end
  end

  dx, dy = L1[1]-XX[1][1] ,  L1[2]-XX[1][2]
  res[1] = (thread_data.hackazi[2][1] - atan(dy, dx))^2 +
            (meas[1] - norm([dx; dy]))^2
  nothing
end

function (cfo::CalcFactor{<:SASBearing2D})( meas,
                                            L1,
                                            XX...  )
  #
  sas2d = cfo.factor
  userdata = cfo.metadata
  idx = cfo._sampleIdx
  #
  dx, dy, azi = 0.0, 0.0, 0.0
  thread_data = sas2d.threadreuse[Threads.threadid()]
  res = [0.0;]

  if string(userdata.solvefor)[1] == 'l'
    # Solving all poses to Landmark
    # Reference the filter positions locally (first element)

    # looks like res is ignored later, but used in this function??
    forwardsas(sas2d, thread_data, res, idx,meas,L1,XX...)

  else

    if thread_data.oncebackidx[1] != idx
      # ccall(:jl_, Void, (Any,), "idx=$(idx)\n")
      thread_data.oncebackidx[1] = idx
      #run for eacn new idx
      thread_data.workvarlist = userdata.variablelist[2:end] # skip landmark element
      thread_data.looelement = findfirst(thread_data.workvarlist .== userdata.solvefor)
      thread_data.elementset = setdiff(Int[1:length(thread_data.workvarlist);],   [thread_data.looelement;])

      # Reference first or second element
      for i in 1:length(XX)
        if thread_data.looelement==1
            thread_data.arrayPos[i,1:2] = XX[i][1:2]-XX[2][1:2]
            dx, dy = L1[1]-XX[2][1] ,  L1[2]-XX[2][2]
        else
            thread_data.arrayPos[i,1:2] = XX[i][1:2]-XX[1][1:2]
            dx, dy = L1[1]-XX[1][1] ,  L1[2]-XX[1][2]
        end
      end

      for i in 1:sas2d.cfgLIE.dataLen
          for j in 1:length(thread_data.elementset)
            thread_data.waveformsLIE[i,j] = sas2d.waveformsIn[i,thread_data.elementset[j]]
          end
          thread_data.waveformsLOOc[i] = sas2d.waveformsIn[i,thread_data.looelement]
      end

    end #end of idx

    # Reference first or second element
    for i in 1:length(XX)
      if thread_data.looelement==1
          thread_data.arrayPos[i,1:2] = XX[i][1:2]-XX[2][1:2]
          dx, dy = L1[1]-XX[2][1] ,  L1[2]-XX[2][2]
      else
          thread_data.arrayPos[i,1:2] = XX[i][1:2]-XX[1][1:2]
          dx, dy = L1[1]-XX[1][1] ,  L1[2]-XX[1][2]
      end
    end

    for i in 1:length(thread_data.elementset), j in 1:2
        elem = thread_data.elementset[i]
        thread_data.arrayPosLIE[i,j] = thread_data.arrayPos[elem,j]
    end
    azi = atan(dy,dx)
    sas2d.cfgLIE.azimuths = [azi;]

    constructCBFFilter2D!(sas2d.cfgLIE, thread_data.arrayPosLIE, thread_data.CBFLIE,thread_data.sourceXY,thread_data.LIEtemp)

    liebf!(thread_data.BFOutLIE, thread_data.waveformsLIE, thread_data.CBFLIE, 1, thread_data.temp, normalize=true)

    copy!(thread_data.phaseshiftLOO,thread_data.waveformsLOOc)
    phaseShiftSingle!(thread_data.sourceXY, sas2d.cfgLIE, azi, thread_data.arrayPos[thread_data.looelement,1:2], thread_data.phaseshiftLOO)

    res[1] = 0.0
    copy!(thread_data.BFOutLIETemp,thread_data.BFOutLIE)
    @fastmath @inbounds @simd for i in 1:length(thread_data.phaseshiftLOO)
      thread_data.BFOutLIETemp[i] += thread_data.phaseshiftLOO[i]
      res[1] -= norm(thread_data.BFOutLIETemp[i])
    end
  end
  res[1]
end

function reset!(dbg::SASDebug)
  dbg.beams = [Float64[];]
  dbg.azi_smpls = Float64[]
  nothing
end

function DFG.compare(a::SASDebug,b::SASDebug)
  TP = true
  TP = TP && a.beams == b.beams
  TP = TP && a.azi_smpls == b.azi_smpls
  return TP
end

mutable struct PackedSASBearing2D <: DFG.AbstractPackedFactor
    rangemodel::String
    totalPhones::Int
    wavedataRawV::Vector{Float64}
    wavedataRawV_dim::Int
    cfgJson::String
    cfgTotal::String
    cfgLIE::String
    waveformsInReal::Vector{Float64}
    waveformsInIm::Vector{Float64}
    PackedSASBearing2D() = new()
    PackedSASBearing2D(rangemodel::String, tp::Int, wavedataRawV::Vector{Float64}, wavedataRawV_dim::Int, cfgJson::String, cfgTotal::String, cfgLIE::String, waveformsInReal::Vector{Float64}, waveformsInIm::Vector{Float64}) = new(rangemodel, tp, wavedataRawV, wavedataRawV_dim, cfgJson, cfgTotal, cfgLIE, waveformsInReal, waveformsInIm)
    PackedSASBearing2D(rangemodel::String, tp::Int, wd::Array{Float64,2}, cf::Dict, cfgTotal::String, cfgLIE::String, waveformsInReal::Vector{Float64}, waveformsInIm::Vector{Float64}) = new(rangemodel, tp, wd[:], size(wd,1), JSON2.write(cf), cfgTotal, cfgLIE, waveformsInReal, waveformsInIm)
end

# Note: Need this otherwise it uses a default converter and causes type conversion issues.
# Safest to just do it here to guarantee it's done.
import Base.convert

function convert(::Type{<:PackedSASBearing2D}, d::SASBearing2D)
  # range model is vector, packing it cleanly
  rangeString = ""
  if typeof(d.rangemodel) == Vector{AliasingScalarSampler}
    rangeString = JSON.json(map(rm -> string(rm), d.rangemodel))
  else
    error("Can't pack range model of type $(typeof(d.rangemodel)) yet.")
  end

  cfgTotalJson = JSON2.write(d.cfgTotal)
  cfgLIEJson = JSON2.write(d.cfgLIE)
  waveformsInReal = map(r -> real(r), d.waveformsIn[:])
  waveformsInIm = map(r -> imag(r), d.waveformsIn[:])
  totalPhones = size(d.waveformsRaw,2)
  pPacked = PackedSASBearing2D(rangeString, totalPhones, d.waveformsRaw, d.cfg, cfgTotalJson, cfgLIEJson, waveformsInReal, waveformsInIm)
  return pPacked
end

function convert(::Type{SASBearing2D}, d::PackedSASBearing2D)::SASBearing2D
  rangePacked = JSON2.read(d.rangemodel)
  rangeModel = map(r -> extractdistribution(r), rangePacked)
  cfgJson = JSON2.read(d.cfgJson, Dict)
  cfgTotal = JSON2.read(d.cfgTotal, CBFFilterConfig)
  cfgLIE = JSON2.read(d.cfgLIE, CBFFilterConfig)
  waveformsIn = map(i -> Complex{Float64}(d.waveformsInReal[i], d.waveformsInIm[i]), 1:(length(d.waveformsInReal)))
  waveformsIn = reshape(waveformsIn, Int(length(waveformsIn)/d.totalPhones), :)

  sas2d = SASBearing2D()
  sas2d.rangemodel = rangeModel
  sas2d.cfg = cfgJson
  sas2d.cfgTotal = cfgTotal
  sas2d.cfgLIE = cfgLIE
  sas2d.waveformsIn = waveformsIn
  sas2d.debugging = false
  sas2d.waveformsRaw = reshape(d.wavedataRawV, d.wavedataRawV_dim, :)
  prepareSASMemory(sas2d)

  return sas2d
end

function DFG.compare(a::SASBearing2D, b::SASBearing2D)
  TP = true
  TP &= compare(a.cfgTotal, b.cfgTotal)               # CBFFilterConfig
  @debug "cfgTotal: $(compare(a.cfgTotal, b.cfgTotal))"
  TP &= compare(a.cfgLIE, b.cfgLIE)                   # CBFFilterConfig
  @debug "cfgLIE: $(compare(a.cfgLIE, b.cfgLIE))"
  TP &= a.waveformsIn == b.waveformsIn         # Array{Complex{Float64}}
  @debug "waveformsIn: $(a.waveformsIn == b.waveformsIn)"

  # Rangemodel
  # Utility dist formula
  calcArrayDist = (a,b, percAcceptable=1) -> begin
    sum(a) == 0 && sum(b) == 0 && return false
    sum(a) == 0 && return false
    sum(abs.(a-b))/sum(abs.(a)) * 100.0 < percAcceptable && return true
    return false
  end
  TP &= typeof(a.rangemodel) == typeof(b.rangemodel)
  @debug "typeof(rangemodel): $((typeof(a.rangemodel) == typeof(b.rangemodel)))"
  if typeof(a.rangemodel) == typeof(b.rangemodel) == Array{AliasingScalarSampler,1}
    TP &= length(a.rangemodel) == length(b.rangemodel)
    @debug "rangemodel (length): $(length(a.rangemodel) == length(b.rangemodel))"

    if length(a.rangemodel) == length(b.rangemodel)
      for i in 1:length(a.rangemodel)
        # Domain
        TP &= calcArrayDist(a.rangemodel[i].domain, b.rangemodel[i].domain)
        @debug "rangemodel (model[$i].domain): $(calcArrayDist(a.rangemodel[i].domain, b.rangemodel[i].domain))"
        # Weights
        TP &= calcArrayDist(a.rangemodel[i].weights, b.rangemodel[i].weights)
        @debug "rangemodel (model[$i].weights): $(calcArrayDist(a.rangemodel[i].weights, b.rangemodel[i].weights))"
      end
    end
  elseif typeof(a.rangemodel) == typeof(b.rangemodel) == Rayleigh
    TP &= a.rangemodel == b.rangemodel
    @debug "rangemodel (Raleigh): $(a.rangemodel == b.rangemodel)"
  end

  TP &= a.cfg == b.cfg                         # Dict
  @debug "cfg: $(a.cfg == b.cfg)"
  TP &= a.waveformsRaw == b.waveformsRaw       # Array{Float64,2}
  @debug "waveformsRaw: $(a.waveformsRaw == b.waveformsRaw)"
  TP &= a.debugging == b.debugging                       # SASDebug
  @debug "debugging: $(a.debugging == b.debugging)"
  return TP
end


text:
module OscillatoryIntegralsODE

using StaticArrays
using OrdinaryDiffEq
import SpecialFunctions: besselj, sphericalbesselj
import LinearAlgebra: dot

# T is the output type
# Tpar is the type of the parameters (not necessarily the same as output)
# M is the dimension of the Levin ODE, changes for each function
abstract type OscillatoryIntegral{T, Tpar, M} end

"""
    levintegrate(oi, f, a, b, alg=Vern9(); ode_kwargs...)

This computes ``\\int_a^b f(x) w(x) dx`` for ``0 < a < b``, where ``f(x)`` is
a smooth, univariate function, and ``w(x)`` is a highly-oscillatory function.

# Arguments:
- `oi::OscillatoryIntegral`: parameters of the highly-oscillatory integral
- `f`: the function to be integrated
- `a`: lower integration bound, must be greater than zero
- `b`: upper integration bound, must be greater than `a`
- `alg`: algorithm used in the ODE solve, you probably want something high order

# Keywords
- Any keyword arguments are passed to the `solve` call of the Levin ODE

# Returns:
- the result of the integral

# Examples
```julia-repl
julia> f(x) = exp(-x^2/16)
       bi = BesselJIntegral{Float64}(100., 100.)  # nu, r
       levintegrate(bi, f, 1.0, 5.0, Vern9(); abstol=1e-6, reltol=1e-6)

0.006311599451652101
```
"""
function levintegrate end

# M=2 case, where the Levin ODE is two-dimensional, like Bessel
function levintegrate(oi::OscillatoryIntegral{T, Tpar, 2}, f::F, a, b, alg=Vern9();
                      ode_kwargs...) where {T, Tpar, F}
    levin_static(u, p_oi, t) =
        SA[zero(T); f(t)] .- transpose(osc_kernel(p_oi, t)) * u  # Levin94 eq 2.6
    va = SA[zero(T); zero(T)]  # arbitrary, Levin prescribes no BC
    prob = ODEProblem{false}(levin_static, va, (a, b), oi)  # set up ODE from a to b
    sol = solve(prob, alg, save_start=false, save_everystep=false; ode_kwargs...)
    vb = last(sol)
    return dot(osc_funcs(oi, b), vb) # - dot(osc_funcs(oi, a), va)  # Levin94 eq 2.9, we chose va = 0
end

# OSCILLATORY FUNCTION DEFINITIONS ==================================

# kernel and functions for the 𝐁𝐞𝐬𝐬𝐞𝐥 function J_ν(rx)
struct BesselJIntegral{T, Tpar} <: OscillatoryIntegral{T, Tpar, 2}
    ν::Tpar
    r::Tpar
end
BesselJIntegral{T}(ν::Tpar, r::Tpar) where {T, Tpar} = BesselJIntegral{T,Tpar}(ν, r)
osc_kernel(oi::BesselJIntegral{T, Tpar}, x) where {T, Tpar} =
    SA[(oi.ν-one(Tpar))/x  (-oi.r);
       (oi.r)            -oi.ν/x]  # Levin94 eq 3.3 (the matrix A(x))
osc_funcs(oi::BesselJIntegral{T, Tpar}, x) where {T, Tpar} =
    SA[besselj(oi.ν-one(Tpar), oi.r * x);  besselj(oi.ν, oi.r * x)]  # Levin94 eq 3.1-3.3, v(x)

"""
    BesselJIntegral{T}(ν, r)

This constructs a Bessel integration problem which solves
```math
I = \\int_a^b f(x) J_{\\nu}(r x)\\, dx
```
and produces output of type `{T}`.
"""
function BesselJIntegral end

# kernel and functions for the 𝐬𝐩𝐡𝐞𝐫𝐢𝐜𝐚𝐥 𝐁𝐞𝐬𝐬𝐞𝐥 function j_ν(rx)
struct SphericalBesselJIntegral{T, Tpar} <: OscillatoryIntegral{T, Tpar, 2}
    ν::Tpar
    r::Tpar
end
SphericalBesselJIntegral{T}(ν::Tpar, r::Tpar) where {T, Tpar} =
    SphericalBesselJIntegral{T,Tpar}(ν, r)
osc_kernel(oi::SphericalBesselJIntegral{T, Tpar}, x) where {T, Tpar} =
    SA[((oi.ν-one(Tpar))/x)  (-oi.r);
        (oi.r)            ((-oi.ν-one(Tpar))/x)]  # the matrix A(x)
osc_funcs(oi::SphericalBesselJIntegral{T, Tpar}, x) where {T, Tpar} =  # the vector v(x)
    SA[sphericalbesselj(oi.ν-one(Tpar), oi.r*x);  sphericalbesselj(oi.ν, oi.r*x)]

"""
    SphericalBesselJIntegral{T}(ν, r)

This constructs a spherical Bessel integration problem which solves
```math
I = \\int_a^b f(x) j_{\\nu}(r x)\\, dx
```
and produces output of type `{T}`.

"""
function SphericalBesselJIntegral end

include("harmonic.jl")
export levintegrate, BesselJIntegral, SphericalBesselJIntegral, HarmonicIntegral

end

max_stars_repo_name: xzackli/OscillatoryIntegralsODE

max_issues_repo_name: xzackli/OscillatoryIntegralsODE

max_forks_repo_name: xzackli/OscillatoryIntegralsODE


text:
"""
    add(x, y)

Adds `x` and `y` together
"""
function add(x, y)
    x + y
end


text:
###############################################################################
#                       Script to generate Maple matrices and equations       #
#                       To be use jointly with Procs.mpl                      #
###############################################################################
# To use it, create a file defining your networks, say networks.jl, and execute:
#
# julia networks.jl
#
# Examples of files defining networks are M1.jl, M2.jl and M1&2.jl

using Nemo
using Polymake
import LinearAlgebra: I, dot

# Nemo.isnonzero(x) = (x) != (zero(x))
# use !iszero instead (it is performance equivalent and more clean)

function findfirstnonzero(A)
    return findfirst((!iszero), A)
end

function findpivots(W)
    # return findfirstnonzero.(eachrow(W)) ## it works just for matrices, code below is more generic.
    # return mapslices(findfirstnonzero, W, dims=collect(2:ndims(W)))
    return n -> findfirstnonzero.(eachslice(W, dims=n))
end

anynonzero(A) = any((!iszero).(A))

function nonzeroslices(M)
    # return n->vec(mapslices(anynonzero, M, dims=deleteat!(collect(1:ndims(M)),n)))
    return n -> anynonzero.(eachslice(M, dims=n))
end

function dropzeroslices(M)
    indeces = nonzeroslices(M).(1:ndims(M))
    return M[indeces...]
end

function stoichiometricsources(net::AbstractMatrix, xs::AbstractVector, getcoefficient)
    return Int.([getcoefficient(p, x) for x in xs, p in net[:,1]])
end

function stoichiometrictargets(net::AbstractMatrix, xs::AbstractVector, getcoefficient)
    return Int.([getcoefficient(p, x) for x in xs, p in net[:,2]])
end

function stoichiometriccoeffs(net::AbstractMatrix, xs::AbstractVector, getcoefficient)
    return hcat(stoichiometricsources(net, xs, getcoefficient), stoichiometrictargets(net, xs, getcoefficient))
end

function stoichiometriccoeffs(net::AbstractMatrix{T}, xs::AbstractVector) where {T <: RingElem}
    return hcat(stoichiometricsources(net, xs, Nemo.coeff), stoichiometrictargets(net, xs, Nemo.coeff))
end

function stoichiometriccoeffs(net::AbstractMatrix, xs::AbstractVector)
    R = parent(xs[1])
    return stoichiometriccoeffs(R.(net), xs)
end

function stoichiometricsources(net::AbstractMatrix)
    return net[:,1:div(end, 2)]
end

function stoichiometrictargets(net::AbstractMatrix)
    return net[:,div(end, 2) + 1:end]
end

# Matrix Y
function kineticorder(net::AbstractMatrix)
    return dropzeroslices(stoichiometricsources(net))
end

# Matrix N
function stoichiometricmatrix(net::AbstractMatrix)
    return dropzeroslices(stoichiometrictargets(net) - stoichiometricsources(net))
end

function conservativelaws(N::AbstractMatrix{T}) where {T <: Integer}
    Nnemo = matrix(FlintIntegerRing(), N)
    nTs, W = left_kernel(Nnemo)
    W = hnf(W)
    return nTs, Int.(Array(W))
end

function cone_positivenullspace(N::AbstractMatrix{T}) where {T <: Integer}
    ## Migrating to Nemo
    Nnemo = Nemo.matrix(Nemo.FlintZZ, N)
    ## Computing nullsapce
    r, U = Nemo.nullspace_right_rational(Nnemo)
    ## Comming back to julia
    nullsp = T.(Array(U[:,1:r]))
    ## The vector space generated by nullsp as a cone
    rays1 = transpose(hcat(nullsp, -nullsp))
    c1 = Polymake.polytope.Cone(INPUT_RAYS=rays1)
    ## The nonnegative orthant as a cone
    rays2 = Matrix{T}(I, size(N, 2), size(N, 2))
    c2 = Polymake.polytope.Cone(INPUT_RAYS=rays2)
    ## Compute intersection of c1 c2 and save RAYS
    d = Polymake.polytope.intersection(c1, c2).RAYS
    ## Convert to integers
    d = Polymake.@convert_to Matrix{Integer} d
    ## Return the intersection as a Base.Array matrix of integers
    return T.(transpose(Array(d)))
end

function matrixtoMaple(io, M, name)
    write(io, "\n\n$(name) := Matrix$(size(M)):\n\n")
    for i in CartesianIndices(M)
        if M[i] != 0
            write(io, "$(name)[$(i[1]),$(i[2])] := $(M[i]):\n")
        end
    end
end

function xstoMaple(io, stoichiometricsources)
    xs = (1:(size(stoichiometricsources, 1)))[nonzeroslices(stoichiometricsources)(1)]
    write(io, "\n\nnxs := $(size(xs, 1)):\n")
    write(io, "xs := [seq(x[i], i = [")
    for x in xs[1:(end - 1)]
        write(io, "$x, ")
    end
    write(io, "$(xs[end])])];\n\n")
    write(io, "depenvars := [seq(x[i], i = [")
    for x in xs[2:(end - 2)]
        write(io, "$x, ")
    end
    write(io, "$(xs[end - 1])])];\n\n")
end

function kstoMaple(io, stoichiometricsources)
    ks = (1:(size(stoichiometricsources, 2)))[nonzeroslices(stoichiometricsources)(2)]
    write(io, "\n\nnks := $(size(ks, 1)):\n")
    write(io, "ks := [seq(k[i], i = [")
    for k in ks[1:(end - 1)]
        write(io, "$k, ")
    end
    write(io, "$(ks[end])])];\n\n")
end

function systemtoMaple(io)
    write(io, "v := Velocities(Y,xs):\n")
    write(io, "digK := Matrix(nks):\n")
    write(io, "for i to nks do digK[i,i] := ks[i] end do:\n")
    write(io, "S := N.digK.v:\n")
    write(io, "Seq:= equfy(S);\n")
end

function WsystemtoMaple(io, nts, W)
    write(io, "\n\nSw := copy(S):\n")
    write(io, "Wx := (W.(Vector[column](xs))) - Vector[column]([seq(T[i], i = 1 .. $(nts))]):\n")
    for (i, p) in enumerate(findpivots(W)(1))
        write(io, "Sw[$p] := Wx[$i]:\n")
    end
    write(io, "Sweq:= equfy(Sw):\n")
    write(io, "J := VectorCalculus[Jacobian](Sw, xs):\n")
    write(io, "DJ := (-1)^(Rank(N))*Determinant(J):\n")
end

## It needs Y and E be defined in Maple
function ConvexparamtoMaple(io)
    write(io, "\n\ndigL := Matrix(convert(E.(Vector[column]([seq(lambda[i], i=1..LinearAlgebra[ColumnDimension](E))])), Vector[row]), shape = diagonal):\n")
    write(io, "digH := DiagonalMatrix([seq(h[i], i = 1..LinearAlgebra[ColumnDimension](LinearAlgebra[Transpose](Y)))]):\n")
    write(io, "Jconv := N.digL.LinearAlgebra[Transpose](Y).digH:\n")
end

macro matrixtoMaple(io, M)
    return quote
        matrixtoMaple($(esc(io)), $(esc(M)), $(esc((string(M)))))
    end
end

function toMaple(net, nxs, file::String)
    S = stoichiometriccoeffs(net, nxs) ## S may have zero rows and cols
    N = stoichiometricmatrix(S) ## Has no zero row or col
    nts, W = conservativelaws(N)
    Y = kineticorder(S)
    E = cone_positivenullspace(N)
    open(file, "w") do io
        write(io, "read(\"ModelsMatricies/Procs.mpl\"):\n")
        @matrixtoMaple io Y
        @matrixtoMaple io N
        @matrixtoMaple io W
        @matrixtoMaple io E
        xstoMaple(io, stoichiometricsources(S))
        kstoMaple(io, stoichiometricsources(S))
        systemtoMaple(io)
        WsystemtoMaple(io, nts, W)
        ConvexparamtoMaple(io)
    end
end

# function toMaple(net, nxs, file::String)
#     S = stoichiometriccoeffs(net, nxs) ## S may have zero rows and cols
#     N = stoichiometricmatrix(S) ## Has no zero row or col
#     nts, W = conservativelaws(N)
#     open(file, "w") do io
#         write(io, "read(\"ModelsMatricies/Procs.mpl\"):\n")
#         matrixtoMaple(io, kineticorder(S), "Y")
#         matrixtoMaple(io, N, "N")
#         matrixtoMaple(io, W, "W")
#         matrixtoMaple(io, cone_positivenullspace(N), "E")
#         xstoMaple(io, stoichiometricsources(S))
#         kstoMaple(io, stoichiometricsources(S))
#         systemtoMaple(io)
#         WsystemtoMaple(io, nts, W)
#         ConvexparamtoMaple(io)
#     end
# end

max_stars_repo_name: LauraBMo/CRNTscripts

max_issues_repo_name: LauraBMo/CRNTscripts

max_forks_repo_name: LauraBMo/CRNTscripts


text:
using LinearAlgebra

LinearAlgebra.diag(a::AbstractJaxArray, k::Integer=0) = numpy.diag(a, -k)

LinearAlgebra.kron(a::AbstractJaxArray, b::AbstractJaxArray) = numpy.kron(a, b)

LinearAlgebra.tr(a::AbstractJaxMatrix, diag=0) = numpy.trace(a, offset=-diag)

LinearAlgebra.dot(a::AbstractJaxVector, b::AbstractJaxVector) = numpy.dot(a, b)


text:
module OutputRelaxationmod

using OutputALASmod

type OutputRelaxation

 solved::Int64 #renvoi 0 si réussi
 solve_message::AbstractString #message de sortie
 xtab::Array #tableau des itérés
 inner_output_alas::Any #en fait une liste de OutputALAS

 rtab::Array{Float64,1}
 stab::Array{Float64,1}
 ttab::Array{Float64,1}
 prectab::Array{Float64,1}

 realisabilite::Array{Float64,1}
 objtab::Array{Float64,1}

end

#Initialisation:
function OutputRelaxation(x0::Vector,realisabilite::Float64,obj::Float64)
 solved=0
 solve_message="Success"
 xtab=Array{Float64,2}
 xtab=collect(x0)
 inner_output=[]
 rtab=[]
 stab=[]
 ttab=[]
 prectab=[]

   if length(x0)<=4
    print_with_color(:green, "j: 0 (r,s,t)=(-,-,-) eps=- xj=$(x0) f(xj)=$(obj) \|c(xj)\|=$(realisabilite)\n\n")
   else
    print_with_color(:green, "j: 0 (r,s,t)=(-,-,-) eps=- f(xj)=$(obj) \|c(xj)\|=$(realisabilite)\n\n")
   end

 return OutputRelaxation(solved,solve_message,xtab,inner_output,rtab,stab,ttab,prectab,[realisabilite],[obj])
end

#Après une itération :
function UpdateOR(or::OutputRelaxation,xk::Vector,stat::Int64,r::Float64,s::Float64,t::Float64,prec::Float64,realisabilite::Float64,outputalas::OutputALASmod.OutputALAS,obj::Float64)
 if stat==0
  or.solved=0
  or.solve_message="Success"
 else #stat !=0
  or.solved=-1
  or.solve_message="Fail"
 end

 or.xtab=[or.xtab xk]
 or.inner_output_alas=[or.inner_output_alas;outputalas]
 or.rtab=[or.rtab;r]
 or.stab=[or.stab;s]
 or.ttab=[or.ttab;t]
 or.prectab=[or.prectab;prec]
 or.realisabilite=[or.realisabilite;realisabilite]
 or.objtab=[or.objtab;obj]

Print(or,length(xk),1,j=length(or.rtab))

 return or
end

function UpdateOR(or::OutputRelaxation,xk::Vector,stat::Int64,r::Float64,s::Float64,t::Float64,prec::Float64,realisabilite::Float64,output,obj::Float64)
 if stat==0
  or.solved=0
  or.solve_message="?"
 else #stat !=0
  or.solved=-1
  or.solve_message="Fail"
 end

 or.xtab=[or.xtab xk]
 or.inner_output_alas=[or.inner_output_alas;output]
 or.rtab=[or.rtab;r]
 or.stab=[or.stab;s]
 or.ttab=[or.ttab;t]
 or.prectab=[or.prectab;prec]
 or.realisabilite=[or.realisabilite;realisabilite]
 or.objtab=[or.objtab;obj]

 return or
end

function UpdateFinalOR(or::OutputRelaxation,mess::String)
 or.solve_message=mess
 return or
end

function Print(or::OutputRelaxation,n::Int64,verbose::Int64;j::Int64=0)
 if verbose>0
  println("")
  if j==0
   if n<=4
    print_with_color(:green, "j: 0 (r,s,t)=(-,-,-) eps=- xj=$(or.xtab[1:n,j]) f(xj)=$(or.objtab[j]) \|c(xj)\|=$(or.realisabilite[j])\n\n")
   else
    print_with_color(:green, "j: 0 (r,s,t)=(-,-,-) eps=- f(xj)=$(or.objtab[j]) \|c(xj)\|=$(or.realisabilite[j])\n\n")
   end
  else
   while j<= length(or.rtab)
    if n<=4
     print_with_color(:green, "j: $j (r,s,t)=($(or.rtab[j]),$(or.stab[j]),$(or.ttab[j])) eps=$(or.prectab[j]) \n xj=$(or.xtab[1:n,j+1]) f(xj)=$(or.objtab[j+1]) \|c(xj)\|=$(or.realisabilite[j+1])\n\n")
    else
     print_with_color(:green, "j: $j (r,s,t)=($(or.rtab[j]),$(or.stab[j]),$(or.ttab[j])) eps=$(or.prectab[j]) \n f(xj)=$(or.objtab[j+1]) \|c(xj)\|=$(or.realisabilite[j+1])\n\n")
    end
    OutputALASmod.Print(or.inner_output_alas[j],n,verbose)
    j+=1
   end
  end
 end

 return true
end

#end of module
end


text:
export cellvolume

"""
    cellvolume(a, b, c, α, β, γ)

Calculates the cell volume from 6 cell parameters.
"""
cellvolume(a, b, c, α, β, γ) =
    a * b * c * sqrt(sind(α)^2 - cosd(β)^2 - cosd(γ)^2 + 2 * cosd(α) * cosd(β) * cosd(γ))
"""
    cellvolume(l::Lattice)
    cellvolume(c::Cell)

Calculates the cell volume from a `Lattice` or a `Cell`.
"""
cellvolume(lattice::AbstractLattice) = abs(det(lattice.data))
cellvolume(cell::Cell) = cellvolume(cell.lattice)
"""
    cellvolume(g::MetricTensor)

Calculates the cell volume from a `MetricTensor`.
"""
cellvolume(g::MetricTensor) = sqrt(det(g.data))  # `sqrt` is always positive!


text:
function BHAdata(pth_gml::String)
    pth_tile = split(pth_gml,"mastermap-topo")[1]

    dir = readdir(pth_tile)
    pth_bha = joinpath(pth_tile,dir[occursin.("building_heights",dir)][1])
    dir = readdir(pth_bha)

    iDocs = findall(dir.=="docs")[1]
    iBha = intersect(findall(dir.!="docs"),findall(occursin.(".pdf",dir).==false))[1]
    pth_bha_docs = joinpath(pth_bha,dir[iDocs])
    pth_bha = joinpath(pth_bha,dir[iBha])

    fname = readdir(pth_bha)[1]

    df_headers = DataFrame(load(joinpath(pth_bha_docs,"BHA_Header.csv"),header_exists=false))
    df = DataFrame(load(joinpath(pth_bha,fname),header_exists=false))

    colnames = [Symbol(str) for str in df_headers[1,:]]

    rename!(df,colnames)

    df.OS_TOPO_TOID = parse.(Int64,replace.(df.OS_TOPO_TOID,"osgb"=>""))

    return df
end

function ExtractFeatures(features,children,tag::String;T::Type=String)
    # Init
    lTag = length(tag)+1
    # Find feature
    tagContents = [features[i][children[i].==tag] for i in 1:length(features)]
    # Remove tags
    tagContents = map(tagContent->map(tc->tc[lTag:(end-lTag)],tagContent),tagContents)
    # Change type if specified
    if T!=String
        tagContents = map(tc->parse.(T,tc),tagContents)
    end
    # Convert to Vector if single valued results throughout
    if length(tagContents) == sum(length.(tagContents))
        tagContents = vcat(tagContents...)
    end
    return tagContents
end


text:
module TransferMatrixMethod

using Statistics
using LinearAlgebra
using ..CommonStructures: _get_beam_parameters, PlaneWave, LayerTMMO

export TMMOptics, tmm_optics

# Wrap the output into different types
abstract type Output end
struct Spectra{T1, T2} <: Output where {T1<:Float64, T2<:ComplexF64}
    Rp::Array{T1}; Rs::Array{T1}
    Tp::Array{T1}; Ts::Array{T1}
    ρp::Array{T2}; ρs::Array{T2}
    τp::Array{T2}; τs::Array{T2}
end
struct Field{T1} <: Output where {T1<:Float64}
    emfp::Array{T1}
    emfs::Array{T1}
end
struct Bloch{T1, T2} <: Output where {T1<:ComplexF64, T2<:Float64}
    κp::Array{T1}
    κs::Array{T1}
    ω::Array{T2}
    Λ::T2
    ωl::T2
    ωh::T2
    qz::Array
end
struct Misc{T1, T2} <: Output where {T1<:Float64, T2<:ComplexF64}
    d::Array{T1}
    ℓ::Array{T1}
    nλ0::Array{T1}
    layers_n::Array{T2}
    λ0::T1
end
struct AdmPhase{T1} <: Output where {T1<:ComplexF64}
    ηp::Array{T1}
    ηs::Array{T1}
    δ::Array{T1}
end
struct TMMOptics{T1, T2, T3, T4, T5, T6, T7} <: Output where {T1<:Spectra, T2<:Field, T3<:Bloch, T4<:Misc, T5<:AdmPhase, T6<:PlaneWave, T7<:LayerTMMO}
    Spectra::T1
    Field::T2
    Bloch::T3
    Misc::T4
    AdmPhase::T5
    Beam::T6
    Layers::Array{T7}
end

## Data conditioning and main calculations.
function tmm_optics(
    beam::T1, layers::Array{T2,N2};
    emfflag::T3=false, h::T4=10, pbgflag::T3=false, λ0::T5=-eps(),
) where {T1<:PlaneWave, T2<:LayerTMMO, N2, T3<:Bool, T4<:Int64, T5<:Float64}
    # Check if λ0 ∈ λ
    λ0 = _checkλ0(λ0, beam)
    # Find beam.λ closest to λ0
    idxλ0 = findmin(abs.(beam.λ .- λ0))[2][1]
    # Build the sequence of n and d depending on the input
    nLen = length(vec(layers))
    d = Array{Float64,1}(undef, nLen)
    nλ0 = Array{Float64,1}(undef, nLen)
    layers_n = Array{ComplexF64,2}(undef, length(beam.λ), nLen)
    _build_arrays!(d, layers_n, nλ0, idxλ0, λ0, vec(layers), beam)
    # Get the parameters for beam
    beam_, _ = _get_beam_parameters(beam)
    # Call transfer matrix method
    if emfflag
        tmmout = _transfer_matrix_emf(layers_n, d, beam_, h)
    else
        tmmout = transfer_matrix(layers_n, d, beam_)
    end
    # Photonic band gap for crystals without defects
    if pbgflag & (nLen > 3)
        ω = Array{Float64,1}(undef, length(beam_.λ))
        κp = Array{ComplexF64,2}(undef, (length(beam_.λ), length(beam_.θ)))
        κs = similar(κp)
        qz::Array{Float64,1} = sin.(beam_.θ).*π./2.0 # parallel wavevector qz
        Λ::Float64 = sum(d[2:3])
        d1, d2 = d[2], d[3]
        n0, n1, n2 = nλ0[1], nλ0[2], nλ0[3]
        ωh::Float64 = Λ/π/(d1*n1 + d2*n2)*acos(-abs(n1 - n2)/(n1 + n2))
        ωl::Float64 = Λ/π/(d2*sqrt(n2^2 - n0^2) + d1*sqrt(n1^2 - n0^2))*acos(abs((n1^2*sqrt(n2^2 - n0^2) - n2^2*sqrt(n1^2 - n0^2))/(n1^2*sqrt(n2^2 - n0^2) + n2^2*sqrt(n1^2 - n0^2))))
        _photonic_dispersion!(κp, κs, ω, Λ, beam_, layers_n[idxλ0, 2:3], d[2:3])
    else
        κp = []; κs = []; ω = []; Λ = 0.0; ωl = 0.0; ωh = 0.0; qz = []
    end
    # Provide the multilayer depth considering the h division
    _ℓ::Array{Float64,2} = (d[2:end-1]/h)*ones.(1, h) # outer product
    ℓ = cumsum([0; _ℓ[:]], dims=1)[1:end-1] # remove last from cumsum
    # Return results
    sol = TMMOptics(
            tmmout[1],
            tmmout[2],
            Bloch(κp, κs, ω, Λ, ωl, ωh, qz),
            Misc(d, ℓ, nλ0, layers_n, λ0),
            tmmout[3],
            beam,
            layers
    )
    return sol
end

"""

    Check input λ0 againts beam.λ.

"""
function _checkλ0(λ0, beam::PlaneWave)
    x = isequal(λ0, -eps()) ? mean(beam.λ) : λ0
    return x
end

"""

    Build the sequence of indices of refraction and thicknesses depending on the input. It follows some logic depending on whether nλ0 was input.

"""
function _build_arrays!(
    d::Array{T0,1},
    layers_n::Array{T1,2},
    nλ0::Array{T0,1},
    idxλ0::T4,
    λ0::T0,
    layers::Array{T2,1},
    beam::T5,
) where {T0<:Float64, T1<:ComplexF64, T2<:LayerTMMO, T3<:Bool, T4<:Int64, T5<:PlaneWave}
    @inbounds for s in eachindex(layers)
        # Refractive index
        @views layers_n[:, s] = layers[s].n
        # nλ0 depending on the input
        if !isequal(layers[s].nλ0[1], -eps()) # nλ0 specified
            nλ0[s] = real(layers[s].nλ0[1])
        else # nλ0 not specified
            nλ0[s] = real(layers_n[idxλ0, s])
        end
        # Build thickness depending on the input
        d[s] = layers[s].d
        if layers[s].type == :OT
            d[s] *= λ0/nλ0[s]
        end
    end
    return nothing
end

"""

    Computes the reflection and transmission coefficients and spectra with the transfer matrix method.

"""
function transfer_matrix(
    layers_n::Array{T1,2}, d::Array{T2,1}, beam::T3;
    λLen::T4=length(beam.λ), θLen::T4=length(beam.θ),
) where {T1<:ComplexF64, T2<:Float64, T3<:PlaneWave, T4<:Int64}
    # λ, λLen, θ, θLen = _get_beam_parameters(beam)
    nLen::Int64 = size(layers_n, 2)
    τs = Array{ComplexF64,2}(undef, (λLen, θLen))
    τp = similar(τs); ρs = similar(τs); ρp = similar(τs)
    δ = Array{ComplexF64,3}(undef, (λLen, θLen, nLen))
    ηs = similar(δ); ηp = similar(ηs)
    # Calculation of complex coefficients of reflection, transmission and emf
    _fresnel_coefficients!(ηs, ηp, δ, ρs, ρp, τs, τp, layers_n, d, beam, nLen)
    x = (
         Spectra(abs2.(ρp),
                 abs2.(ρs),
                 real(@view(ηp[:, :, 1]).*@view(ηp[:, :, end])).*abs2.(τp),
                 real(@view(ηs[:, :, 1]).*@view(ηs[:, :, end])).*abs2.(τs), ρp, ρs, τp, τs),
         Field([], []),
         AdmPhase(ηp, ηs, δ),
    )
    return x
end

"""

    Compute the Fresnel coefficients.

"""
function _fresnel_coefficients!(
    ηs::Array{T1,3}, ηp::Array{T1,3},
    δ::Array{T1,3},
    ρs::Array{T1,2}, ρp::Array{T1,2},
    τs::Array{T1,2}, τp::Array{T1,2},
    layers_n::Array{T1,2},
    d::Array{T2,1},
    beam::T3,
    nLen::T4,
) where {T1<:ComplexF64, T2<:Float64, T3<:PlaneWave, T4<:Int64}
    cosϕ = Array{ComplexF64,1}(undef, nLen)
    I2d = Array{ComplexF64,2}(I,2,2) # Identity 2x2 matrix
    @inbounds for a in eachindex(beam.θ), l in eachindex(beam.λ)
        # Angle of incidence inside each layer according to the cosine Snell law, to avoid cutoff of total internal reflection with complex angles
        cosϕ[1] = cos(beam.θ[a])
        for c = 2:nLen
            cosϕ[c] = cosϑ(layers_n[l, c-1], layers_n[l, c], cosϕ[c-1])
        end
        # Admittance of the first medium for both polarizations
        @. ηp[l, a, :] = ζp(@view(layers_n[l, :]), cosϕ)
        @. ηs[l, a, :] = ζs(@view(layers_n[l, :]), cosϕ)
        # Phase shift of the thickness
        @. δ[l, a, :] = 2.0*π*@view(layers_n[l, :])*d*cosϕ/beam.λ[l]
        # Total transfer matrix
        Ψp::Array{ComplexF64,2} = reduce(*, [[I2d]; Φ.(@view(δ[l, a, 2:end-1]), @view(ηp[l, a, 2:end-1]))])
        Ψs::Array{ComplexF64,2} = reduce(*, [[I2d]; Φ.(@view(δ[l, a, 2:end-1]), @view(ηs[l, a, 2:end-1]))])
        # Compute the Fresnell coefficients
        ρs[l, a] = ρ(ηs[l, a, 1], ηs[l, a, end], Ψs)
        ρp[l, a] = ρ(ηp[l, a, 1], ηp[l, a, end], Ψp)
        τs[l, a] = τ(ηs[l, a, 1], ηs[l, a, end], Ψs)
        τp[l, a] = τ(ηp[l, a, 1], ηp[l, a, end], Ψp)
    end
    return nothing
end

"""

    Computes the reflection and transmission coefficients, spectra and the electromagnetic field with the transfer matrix method.

"""
function _transfer_matrix_emf(
    layers_n::Array{T1,2}, d::Array{T2,1}, beam::T3, h::T4;
    λLen::T4=length(beam.λ), θLen::T4=length(beam.θ),
) where {T1<:ComplexF64, T2<:Float64, T3<:PlaneWave, T4<:Int64}
    # λ, λLen, θ, θLen = _get_beam_parameters(beam)
    nLen::Int64 = size(layers_n, 2)
    τs = Array{ComplexF64,2}(undef, (λLen, θLen))
    τp = similar(τs); ρs = similar(τs); ρp = similar(τs)
    δ = Array{ComplexF64,3}(undef, (λLen, θLen, nLen))
    ηs = similar(δ); ηp = similar(ηs)
    emfs = Array{Float64,3}(undef, (λLen, θLen, (length(d) - 2)*h)); emfp = similar(emfs)
    # Calculation of complex coefficients of reflection, transmission and emf
    _fresnel_coefficients_emf!(ηs, ηp, δ, ρs, ρp, τs, τp, emfs, emfp, layers_n, d, beam, nLen, h)
    x = (
         Spectra(abs2.(ρp),
                 abs2.(ρs),
                 real(@view(ηp[:, :, 1]).*@view(ηp[:, :, end])).*abs2.(τp),
                 real(@view(ηs[:, :, 1]).*@view(ηs[:, :, end])).*abs2.(τs), ρp, ρs, τp, τs),
         Field(emfp, emfs),
         AdmPhase(ηp, ηs, δ),
    )
    return x
end

"""

    Compute the Fresnel coefficients and the electromagnetic field.

"""
function _fresnel_coefficients_emf!(
    ηs::Array{T1,3}, ηp::Array{T1,3},
    δ::Array{T1,3},
    ρs::Array{T1,2}, ρp::Array{T1,2},
    τs::Array{T1,2}, τp::Array{T1,2},
    emfs::Array{T2,3}, emfp::Array{T2,3},
    layers_n::Array{T1,2},
    d::Array{T2,1},
    beam::T3,
    nLen::T4,
    h::T4,
) where {T1<:ComplexF64, T2<:Float64, T3<:PlaneWave, T4<:Int64}
    cosϕ = Array{ComplexF64,1}(undef, nLen)
    I2d = Array{ComplexF64,2}(I, 2, 2) # Identity 2x2 matrix
    @inbounds for a in eachindex(beam.θ), l in eachindex(beam.λ)
        # Angle of incidence inside each layer according to the cosine Snell law, to avoid cutoff of total internal reflection with complex angles
        cosϕ[1] = cos(beam.θ[a])
        for c = 2:nLen
            cosϕ[c] = cosϑ(layers_n[l, c-1], layers_n[l, c], cosϕ[c-1])
        end
        # Admittance of the first medium for both polarizations
        @. ηp[l, a, :] = ζp(@view(layers_n[l, :]), cosϕ)
        @. ηs[l, a, :] = ζs(@view(layers_n[l, :]), cosϕ)
        # Phase shift of the thickness
        @. δ[l, a, :] = 2.0*π*@view(layers_n[l, :])*d*cosϕ/beam.λ[l]
        # Total transfer matrix
        Ψp::Array{ComplexF64,2} = reduce(*, [[I2d]; Φ.(@view(δ[l, a, 2:end-1]), @view(ηp[l, a, 2:end-1]))])
        Ψs::Array{ComplexF64,2} = reduce(*, [[I2d]; Φ.(@view(δ[l, a, 2:end-1]), @view(ηs[l, a, 2:end-1]))])
        # Compute the Fresnell coefficients
        ρs[l, a] = ρ(ηs[l, a, 1], ηs[l, a, end], Ψs)
        ρp[l, a] = ρ(ηp[l, a, 1], ηp[l, a, end], Ψp)
        τs[l, a] = τ(ηs[l, a, 1], ηs[l, a, end], Ψs)
        τp[l, a] = τ(ηp[l, a, 1], ηp[l, a, end], Ψp)
        # Compute the _emfield
        emfs[l, a, :] .= _emfield(@view(layers_n[l, :]), d, @view(δ[l, a, :]), @view(ηs[l, a, :]), Ψs, nLen, h)
        emfp[l, a, :] .= _emfield(@view(layers_n[l, :]), d, @view(δ[l, a, :]), @view(ηp[l, a, :]), Ψp, nLen, h)
    end
    return nothing
end

"""

    Calculates the optical transfer matrix of a layer. (2x2 optical tranfer matrix.)
        δ: phase shift of the layer
        η: admittance of the layer

"""
Φ(δ::T1, η::T1) where {T1<:ComplexF64} = [cos(δ) (-im*sin(δ)/η); (-im*sin(δ)*η) cos(δ)]

"""

    Computes the reflection and transmission coefficients given the admittance and transfer matrix of the whole structure per wavelenth and angle of incidence.

"""
# Reflection coefficient
function ρ(η0::T1, ηm::T1, Ψ::Array{T1,2}) where {T1<:ComplexF64}
    r = (η0*Ψ[1,1] - Ψ[2,1] + η0*ηm*Ψ[1,2] - ηm*Ψ[2,2])/(η0*Ψ[1,1] + Ψ[2,1] + η0*ηm*Ψ[1,2] + ηm*Ψ[2,2])
    return r
end
# Transmission coefficient
function τ(η0::T1, ηm::T1, Ψ::Array{T1,2}) where {T1<:ComplexF64}
    t = 2.0/(η0*Ψ[1,1] + Ψ[2,1] + η0*ηm*Ψ[1,2] + ηm*Ψ[2,2])
    return t
end

"""

    Computes the inverse total transfer matrix for the whole structure at each wavelenth and angle of incidence and return the field.

"""
function _emfield(
    N::SubArray{T1,1},
    d::Array{T2,1},
    δ::SubArray{T1,1},
    η::SubArray{T1,1},
    Ψ::Array{T1,2},
    nLen::T3,
    h::T3,
) where {T1<:ComplexF64, T2<:Float64, T3<:Int64}
    m0 = Array{ComplexF64,2}(undef, 2, 2)
    m1 = Array{ComplexF64,2}(I, 2, 2) # Identity 2x2 matrix
    g11 = Array{ComplexF64,1}(undef, (nLen - 2)*h); g12 = similar(g11)
    # Divide the phase shift by h but keep η as is for each layer
    mδ::Array{ComplexF64} = δ/h
    _m1::Array = Ξ.(mδ, η)
    _g_elements!(g11, g12, m0, _m1, m1, Ψ, nLen, h)
    fi = _field_intensity(g11, g12, η[1], η[end], Ψ)
    return fi
end

"""

	Compute the elements of the G matrix.

"""
function _g_elements!(
    g11::Array{T1,1}, g12::Array{T1,1},
    m0::Array{T1,2},
    _m1::Array{Array{T1,2},1},
    m1::Array{T1,2},
    Ψ::Array{T1,2},
    nLen::T2,
    h::T2,
) where {T1<:ComplexF64, T2<:Int64}
    _h = collect(1:h)
    @inbounds for c =2:nLen-1, j in eachindex(_h)
        k::Int64 = h*(c - 2) + j
        m1 = _m1[c]*m1
        m0 = m1*Ψ
        g11[k] = m0[1, 1]
        g12[k] = m0[1, 2]
    end
    return nothing
end

"""

	Compute the field intensity.

"""
function _field_intensity(
    g11::Array{T1,1}, g12::Array{T1,1}, η1::T1, ηf::T1, Ψ::Array{T1,2},
) where {T1<:ComplexF64}
    fi = abs2.((g11 .+ ηf.*g12)./(0.25*(η1*Ψ[1,1] + Ψ[2,1] + η1*ηf*Ψ[1,2] + ηf*Ψ[2,2])))
    return fi
end

"""

	Calculates the inverse of optical transfer matrix of a layer.

        T = Ξ(φ, η)

            φ:  phase shift of the layer
            η: admittance of the layer

            T: 2x2 optical tranfer matrix.

"""
Ξ(φ::T1, η::T1) where {T1<:ComplexF64} = [cos(φ) (im*sin(φ)/η); (im*sin(φ)*η) cos(φ)]

"""

    Computes the photonic dispersion of binary structures (crystals only) alternating
    two different dielectric layers.

"""
function _photonic_dispersion!(
    κp::Array{T1,2}, κs::Array{T1,2},
    ω::Array{T2,1},
    Λ::T2,
    beam::T3,
    n::Array{T1,1},
    d::Array{T2,1},
) where {T1<:ComplexF64, T2<:Float64, T3<:PlaneWave}
    @. ω = 2*π/beam.λ # Angular frequency
    # Angle of incidence of the second layer with Snell's law of cosine
    cosθ1::Vector{ComplexF64} = cos.(beam.θ)
    cosθ2::Vector{ComplexF64} = cosϑ.(n[1], n[2], cosθ1)
    # Prefactor for Bloch wavevector
    factor_s = _adm_factor.(ζs.(n[1], cosθ1), ζs.(n[2], cosθ2))
    factor_p = _adm_factor.(ζp.(n[1], cosθ1), ζp.(n[2], cosθ2))
    # Bloch wavevectors
    @inbounds for a in eachindex(cosθ1), b in eachindex(ω)
        κp[b, a] = acos(cosκ(d[1]*ω[b]*n[1]*cosθ1[a], d[2]*ω[b]*n[2]*cosθ2[a], factor_p[a]))/Λ
        κs[b, a] = acos(cosκ(d[1]*ω[b]*n[1]*cosθ1[a], d[2]*ω[b]*n[2]*cosθ2[a], factor_s[a]))/Λ
    end
    return nothing
end

"""

	Prefactor for bloch wavevector.

"""
_adm_factor(η1::T1, η2::T1) where {T1<:ComplexF64} = 0.5*(η1^2 + η2^2)/η1/η2

"""

	Bloch wavevector.

"""
cosκ(a1::T1, a2::T1, f::T1) where {T1<:ComplexF64} = cos(a1)*cos(a2) - f*sin(a1)*sin(a2)

"""

	Snell's law in cosine form. Returns the cosine already.

"""
function cosϑ(n1::T1, n2::T1, cosθ::T2) where {T1<:ComplexF64, T2<:Number}
    x = sqrt(1.0 - (n1/n2)^2*(1.0 - cosθ^2))
    return x
end

"""

	Admittance of the medium for p and s polarizations.

"""
ζp(n::T1, cosθ::T2) where {T1<:ComplexF64, T2<:Number} = n/cosθ
ζs(n::T1, cosθ::T2) where {T1<:ComplexF64, T2<:Number} = n*cosθ

end # module


text:
#!/usr/bin/julia

# Trizen
# Date: 15th October 2013
# https://github.com/trizen
# https://trizenx.blogspot.ro/2013/11/smart-word-wrap.html

# Smart word wrap algorithm
# See: https://en.wikipedia.org/wiki/Word_wrap#Minimum_raggedness

# This is the ugliest method! It, recursively,
# prepares the words for the combine() function.
function prepare_words(array, width, callback, depth=0)

    root, len, i = [], 0, 0
    limit = length(array)

    while ((i+=1) <= limit)
        len += (word_len = length(array[i]))

        if len > width
             if word_len > width
                    len -= word_len
                    value = array[i]
                    a = collect((m.match for m = eachmatch(Regex(".{1,$width}"), value)))
                    splice!(array, i, map(x -> convert(String, x), a))
                    limit = length(array)
                    i -= 1
                    continue
                end
            break
        end

        push!(root, Any[
            join(array[1:i], ' '),
            prepare_words(array[i+1 : limit], width, callback, depth+1)
        ])

        if depth == 0
            callback(root[1])
            root = []
        end

        ((len += 1) >= width) && break
    end

    root
end

# This function combines the
# the parents with the childrens.
function combine(root, path, block)
    key = popfirst!(path)
    for value in path
        push!(root, key)
        if isempty(value)
            block(root)
        else
            for item in value
                combine(root, item, block)
            end
        end
        pop!(root)
    end
end

# This is the main function of the algorithm
# which calls all the other functions and
# returns the best possible wrapped string.
function smart_wrap(text, width)

    words = isa(text, AbstractString) ? split(text) : text

    best = Dict{String, Any}(
        "score" => Inf,
        "value" => [],
    )

    prepare_words(words, width, function (path)
        combine([], path, function (combination)
            score = 0
            for line in combination[1:length(combination)-1]
                score += (width - length(line))^2
            end

            if score < best["score"]
                best["score"] = score
                best["value"] = copy(combination)
            end
        end)
    end)

    join(best["value"], "\n")
end

#
## Usage examples
#

text = "aaa bb cc ddddd"
println(smart_wrap(text, 6))

println("-" ^ 80)

text = "As shown in the above phases (or steps), the algorithm does many useless transformations"
println(smart_wrap(text, 20))

println("-" ^ 80)

text = "Lorem ipsum dolor sit amet, consectetur adipiscing elit."
println(smart_wrap(text, 20))

println("-" ^ 80)

text = "Lorem ipsum dolor ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ amet, consectetur adipiscing elit."
println(smart_wrap(text, 20))

println("-" ^ 80)

text = "Will Perl6 also be pre-installed on future Mac/Linux operating systems? ... I can\'t predict the success of the project...";
println(smart_wrap(text, 20))

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts


text:
mutable struct Structure_Ion <: Str
    atom_name
    atom_charge
    para_mass
    num_atom_types
    vec_type_id
    res_name
end


text:
"""
    strel([T], X::AbstractArray)

Convert structuring element (SE) `X` to appropriate presentation format with element type `T`.
This is a useful tool to generate SE that most ImageMorphology functions understand.

ImageMorphology currently supports two commonly used representations:

- `T=CartesianIndex`: offsets to its center point. The output type is
  `Vector{CartesianIndex{N}}`.
- `T=Bool`: connectivity mask where `true` indicates connected to its center point. The
  output type is `BitArray{N}`.

```jldoctest; setup=:(using ImageMorphology)
julia> se_mask = centered(Bool[1 1 0; 1 1 0; 0 0 0]) # connectivity mask
3×3 OffsetArray(::Matrix{Bool}, -1:1, -1:1) with eltype Bool with indices -1:1×-1:1:
 1  1  0
 1  1  0
 0  0  0

julia> se_offsets = strel(CartesianIndex, se_mask) # displacement offsets to its center point
3-element Vector{CartesianIndex{2}}:
 CartesianIndex(-1, -1)
 CartesianIndex(0, -1)
 CartesianIndex(-1, 0)

julia> se = strel(Bool, se_offsets)
3×3 OffsetArray(::BitMatrix, -1:1, -1:1) with eltype Bool with indices -1:1×-1:1:
 1  1  0
 1  1  0
 0  0  0
```

See also [`strel_diamond`](@ref) and [`strel_box`](@ref) for SE constructors for two
special cases.
"""
function strel end

strel(se) = strel(strel_type(se), se)

# convenient user interface without exporting MorphologySE
function strel(::Type{ET}, se::AbstractArray) where {ET<:CartesianIndex}
    return strel(SEOffset{strel_ndims(se)}(), se)
end
strel(::Type{ET}, se::AbstractArray) where {ET<:Bool} = strel(SEMask{strel_ndims(se)}(), se)

# conversion between different SE arrays
function strel(SET::MorphologySE, se::T) where {T}
    strel_type(se) == SET || error("unsupported conversion from type $T to $SET")
    return se
end

function strel(::SEMask{N}, offsets::AbstractArray{CartesianIndex{N}}) where {N}
    isempty(offsets) && return centered(trues(ntuple(_ -> 1, N)))
    mn, mx = extrema(offsets)
    r = ntuple(N) do i
        max(abs(mn.I[i]), abs(mx.I[i]))
    end
    sz = @. 2r + 1
    se = centered(falses(sz))
    se[offsets] .= true
    se[zero(eltype(offsets))] = true # always set center point to true
    return centered(BitArray(OffsetArrays.no_offset_view(se)))
end
strel(::SEMask{N}, mask::AbstractArray{Bool,N}) where {N} = mask

function strel(::SEOffset{N}, connectivity::AbstractArray{Bool,N}) where {N}
    all(isodd, size(connectivity)) || error("`connectivity` must be odd-sized")
    ax = axes(connectivity)
    is_symmetric = all(r -> first(r) == -last(r), ax)
    if !is_symmetric && all(first.(axes(connectivity)) .== 1)
        # To keep consistent with the "kernel" concept in ImageFiltering, we require
        # the connectivity mask to be centered as well.
        # This is a perhaps permanent depwarn to throw friendly message to the user
        # if they're used to use, e.g., `trues(3, 3)` as the input.
        msg = "connectivity mask is expected to be a centered bool array"
        hint = "Do you mean `centered(connectivity)`"
        Base.depwarn("$msg. $hint?", :strel)
    elseif !is_symmetric
        throw(ArgumentError("`connectivity` must be symmetric bool array"))
    end
    connectivity = centered(connectivity)
    # always skip center point
    return [i for i in CartesianIndices(connectivity) if connectivity[i] && !iszero(i)]
end

"""
    strel_type(x)

Infer the structuring element type for `x`.

!!! note "developer note"
    This function is used to dispatch special SE types, e.g., [`SEBoxArray`](@ref), to
    optimized implementation of particular morphology filter. In this sense it is required
    for custom SE array types to define this method.
"""
strel_type(se::MorphologySE) = se
strel_type(::AbstractArray{Bool,N}) where {N} = SEMask{N}()
strel_type(::AbstractVector{CartesianIndex{N}}) where {N} = SEOffset{N}()
strel_type(::CartesianIndices{N}) where {N} = SEOffset{N}()
strel_type(::T) where {T} = error("invalid structuring element data type: $T")

"""
    strel_size(x)

Calculate the minimal block size that contains the structuring element. The result
will be a tuple of odd integers.

```jldoctest; setup=:(using ImageMorphology; using ImageMorphology.StructuringElements)
julia> se = strel_diamond((5, 5); r=1)
5×5 SEDiamondArray{2, 2, UnitRange{$Int}, 0} with indices -2:2×-2:2:
 0  0  0  0  0
 0  0  1  0  0
 0  1  1  1  0
 0  0  1  0  0
 0  0  0  0  0

julia> strel_size(se) # is not (5, 5)
(3, 3)

julia> strel(Bool, strel(CartesianIndex, se)) # because it only checks the minimal enclosing block
3×3 OffsetArray(::BitMatrix, -1:1, -1:1) with eltype Bool with indices -1:1×-1:1:
 0  1  0
 1  1  1
 0  1  0

julia> se = [CartesianIndex(1, 1), CartesianIndex(-2, -2)];

julia> strel_size(se) # is not (4, 4)
(5, 5)

julia> strel(Bool, se) # because the connectivity mask has to be odd size
5×5 OffsetArray(::BitMatrix, -2:2, -2:2) with eltype Bool with indices -2:2×-2:2:
 1  0  0  0  0
 0  0  0  0  0
 0  0  1  0  0
 0  0  0  1  0
 0  0  0  0  0

julia> se = strel_diamond((5, 5), (1, ); r=1)
5×5 SEDiamondArray{2, 1, UnitRange{$Int}, 1} with indices -2:2×-2:2:
 0  0  0  0  0
 0  0  1  0  0
 0  0  1  0  0
 0  0  1  0  0
 0  0  0  0  0

julia> strel_size(se)
(3, 1)
```
"""
strel_size(se) = size(strel(Bool, strel(CartesianIndex, se)))

"""
    strel_ndims(x)::Int

Infer the dimension of the structuring element `x`
"""
strel_ndims(se) = strel_ndims(strel_type(se))
strel_ndims(::MorphologySE{N}) where {N} = N


text:
OLDARGS = copy(ARGS)
atreplinit() do args...
  append!(ARGS,OLDARGS)
  repl = Base.active_repl
  println(STDOUT,"="^Base.Terminals.width(Base.active_repl.t))
  println(STDOUT," "^(div(Base.Terminals.width(Base.active_repl.t),2)-12),"Launching Gallium")
  println(STDOUT,"="^Base.Terminals.width(Base.active_repl.t))
  if (isdefined(repl,:mi) && !isdefined(repl.mi,:interface))
    repl.mi.interface = Base.REPL.setup_interface(repl)
  elseif !isdefined(repl,:interface)
    repl.interface = Base.REPL.setup_interface(repl)
  end
  include("call.jl")
  empty!(Core.ARGS)
end
empty!(Core.ARGS)
Base._start()
exit()


text:
"""
todo
"""
abstract Visualization

Base.show(::Visualization) = error("implement!")


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

mutable struct VirtualMachineCaptureResult <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    schema::Any # spec type: Union{ Nothing, String } # spec name: \$schema
    contentVersion::Any # spec type: Union{ Nothing, String } # spec name: contentVersion
    parameters::Any # spec type: Union{ Nothing, Any } # spec name: parameters
    resources::Any # spec type: Union{ Nothing, Vector{Any} } # spec name: resources

    function VirtualMachineCaptureResult(;id=nothing, schema=nothing, contentVersion=nothing, parameters=nothing, resources=nothing)
        o = new()
        validate_property(VirtualMachineCaptureResult, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(VirtualMachineCaptureResult, Symbol("\$schema"), schema)
        setfield!(o, Symbol("schema"), schema)
        validate_property(VirtualMachineCaptureResult, Symbol("contentVersion"), contentVersion)
        setfield!(o, Symbol("contentVersion"), contentVersion)
        validate_property(VirtualMachineCaptureResult, Symbol("parameters"), parameters)
        setfield!(o, Symbol("parameters"), parameters)
        validate_property(VirtualMachineCaptureResult, Symbol("resources"), resources)
        setfield!(o, Symbol("resources"), resources)
        o
    end
end # type VirtualMachineCaptureResult

const _property_map_VirtualMachineCaptureResult = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"), Symbol("\$schema")=>Symbol("schema"), Symbol("contentVersion")=>Symbol("contentVersion"), Symbol("parameters")=>Symbol("parameters"), Symbol("resources")=>Symbol("resources"))
const _property_types_VirtualMachineCaptureResult = Dict{Symbol,String}(Symbol("id")=>"String", Symbol("\$schema")=>"String", Symbol("contentVersion")=>"String", Symbol("parameters")=>"Any", Symbol("resources")=>"Vector{Any}")
Base.propertynames(::Type{ VirtualMachineCaptureResult }) = collect(keys(_property_map_VirtualMachineCaptureResult))
Swagger.property_type(::Type{ VirtualMachineCaptureResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineCaptureResult[name]))}
Swagger.field_name(::Type{ VirtualMachineCaptureResult }, property_name::Symbol) =  _property_map_VirtualMachineCaptureResult[property_name]

function check_required(o::VirtualMachineCaptureResult)
    true
end

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


text:
using Turing: Tracker

"""
$(TYPEDEF)
Defines a callable that returns the log density for given parameter values when called with
a `NamedTuple` `(parameters = ..., σ = ...)` where `parameters` is a vector of parameters,
and `σ` is the vector of noise scales.
For a common use case, see [`dynamichmc_inference`](@ref).
# Fields
$(FIELDS)
"""
Base.@kwdef struct DynamicHMCPosterior{TA,TP,TD,TT,TR,TS,TK,TI}
    "Algorithm for the ODE solver."
    algorithm::TA
    "An ODE problem definition (`DiffEqBase.DEProblem`)."
    problem::TP
    "Time values at which the simulated path is compared to `data`."
    t::TT
    "Data, as a matrix with each time value in a column."
    data::TD
    "Priors for parameters, an iterable with the same length as the number of parameters."
    parameter_priors::TR
    """
    Priors for the noise scales (currently the standard deviation of a normal distribution),
    one for each variable.
    """
    σ_priors::TS
    "Keyword arguments passed on the the ODE solver `solve`."
    solve_kwargs::TK
    sample_u0::Bool
    save_idxs::TI
end

function (P::DynamicHMCPosterior)(θ)
    @unpack parameters, σ = θ
    @unpack algorithm, problem, data, t, parameter_priors = P
    @unpack σ_priors, solve_kwargs, sample_u0, save_idxs = P
    T = eltype(parameters)
    nu = save_idxs == nothing ? length(problem.u0) : length(save_idxs)
    u0 = convert.(T, sample_u0 ? parameters[1:nu] : problem.u0)
    p = convert.(T, sample_u0 ? parameters[(nu + 1):end] : parameters)
    if length(u0) < length(problem.u0)
      # assumes u is ordered such that the observed variables are in the begining, consistent with ordered theta
      for i in length(u0):length(problem.u0)
          push!(u0, convert(T,problem.u0[i]))
      end
    end
    _saveat = t === nothing ? Float64[] : t
    sol = solve(problem, algorithm; u0=u0, p=p, saveat = _saveat, save_idxs = save_idxs, solve_kwargs...)
    failure = size(sol, 2) < length(_saveat)
    failure && return T(0) * sum(σ) + T(-Inf)
    log_likelihood = sum(sum(map(logpdf, Normal.(0.0, σ), sol[:, i] .- data[:, i])) for (i, t) in enumerate(t))
    log_prior_parameters = sum(map(logpdf, parameter_priors, parameters))
    log_prior_σ = sum(map(logpdf, σ_priors, σ))
    log_likelihood + log_prior_parameters + log_prior_σ
end

# function (P::DynamicHMCPosterior)(θ)
#   @unpack parameters, σ = θ
#   @unpack algorithm, problem, data, t, parameter_priors, σ_priors, solve_kwargs = P
#   prob = remake(problem, u0 = convert.(eltype(parameters), problem.u0), p = parameters)
#   solution = solve(prob, algorithm; solve_kwargs...)
#   any((s.retcode ≠ :Success && s.retcode ≠ :Terminated) for s in solution) && return -Inf
#   log_likelihood = sum(sum(logpdf.(Normal.(0.0, σ), solution(t) .- data[:, i]))
#                        for (i, t) in enumerate(t))
#   log_prior_parameters = sum(map(logpdf, parameter_priors, parameters))
#   log_prior_σ = sum(map(logpdf, σ_priors, σ))
#   log_likelihood + log_prior_parameters + log_prior_σ
# end

"""
$(SIGNATURES)
Run MCMC for an ODE problem. Return a `NamedTuple`, which is similar to the one returned by
`DynamicHMC.mcmc_with_warmup`, with an added field `posterior` which contains a vector of
posterior values (transformed from `ℝⁿ`).
# Arguments
- `problem` is the ODE problem
- `algorithm` is the ODE algorithm
- `t` is the time values at which the solution is compared to `data`
- `data` is a matrix of data, with one column for each element in `t`
- `parameter_priors` is an iterable with the length of the number of paramers, and is used
  as a prior on it, should have comparable structure.
- `parameter_transformations`: a `TransformVariables` transformation to mapping `ℝⁿ` to the
  vector of valid parameters.
# Keyword arguments
- `rng` is the random number generator used for MCMC. Defaults to the global one.
- `num_samples` is the number of MCMC draws (default: 1000)
- `AD_gradient_kind` is passed on to `LogDensityProblems.ADgradient`, make sure to `import`
  the corresponding library.
- `solve_kwargs` is passed on to `solve`
- `mcmc_kwargs` are passed on as keyword arguments to `DynamicHMC.mcmc_with_warmup`
"""
function dynamichmc_inference(problem::DiffEqBase.DEProblem, algorithm, t, data,
                              parameter_priors, parameter_transformations=as(Vector, asℝ₊, length(parameter_priors));
                              σ_priors = fill(Normal(0, 5), size(data, 1)),sample_u0 = false, rng = Random.GLOBAL_RNG,
                              num_samples = 1000, AD_gradient_kind = Val(:ForwardDiff), save_idxs = nothing,solve_kwargs = (),
                              mcmc_kwargs = (initialization = (q = zeros(length(parameter_priors) + (save_idxs === nothing ? length(data[:,1]) : length(save_idxs))),),))
    P = DynamicHMCPosterior(; algorithm = algorithm, problem = problem, t = t, data = data,
                            parameter_priors = parameter_priors, σ_priors = σ_priors,
                            solve_kwargs = solve_kwargs, sample_u0 = sample_u0, save_idxs = save_idxs)
    trans = as((parameters = parameter_transformations,
                σ = as(Vector, asℝ₊, length(σ_priors))))
    ℓ = TransformedLogDensity(trans, P)
    ∇ℓ = LogDensityProblems.ADgradient(AD_gradient_kind, ℓ)
    results = mcmc_with_warmup(rng, ∇ℓ, num_samples; mcmc_kwargs...)
    merge((posterior = transform.(Ref(trans), results.chain), ), results)
end


text:
#! format: off
struct RangeReserve <: AbstractReservesFormulation end
struct StepwiseCostReserve <: AbstractReservesFormulation end
############################### Reserve Variables #########################################

get_variable_sign(_, ::Type{<:PSY.Reserve}, ::AbstractReservesFormulation) = NaN
############################### ActiveServiceVariable, Reserve #########################################

get_variable_binary(::ActiveServiceVariable, ::Type{<:PSY.Reserve}, ::AbstractReservesFormulation) = false
get_variable_upper_bound(::ActiveServiceVariable, ::PSY.Reserve, ::PSY.Component, _) = nothing
get_variable_lower_bound(::ActiveServiceVariable, ::PSY.Reserve, ::PSY.Component, _) = 0.0

############################### ServiceRequirementVariable, ReserveDemandCurve ################################

get_variable_binary(::ServiceRequirementVariable, ::Type{<:PSY.ReserveDemandCurve}, ::AbstractReservesFormulation) = false
get_variable_upper_bound(::ServiceRequirementVariable, ::PSY.ReserveDemandCurve, ::PSY.Component, ::AbstractReservesFormulation) = nothing
get_variable_lower_bound(::ServiceRequirementVariable, ::PSY.ReserveDemandCurve, ::PSY.Component, ::AbstractReservesFormulation) = 0.0

#! format: on
################################## Reserve Requirement Constraint ##########################
function service_requirement_constraint!(
    optimization_container::OptimizationContainer,
    service::SR,
    ::ServiceModel{SR, T},
) where {SR <: PSY.Reserve, T <: AbstractReservesFormulation}
    parameters = model_has_parameters(optimization_container)
    use_forecast_data = model_uses_forecasts(optimization_container)
    initial_time = model_initial_time(optimization_container)
    @debug initial_time
    time_steps = model_time_steps(optimization_container)
    name = PSY.get_name(service)
    constraint =
        get_constraint(optimization_container, make_constraint_name(REQUIREMENT, SR))
    reserve_variable = get_variable(optimization_container, name, SR)
    use_slacks = get_services_slack_variables(optimization_container.settings)

    ts_vector = get_time_series(optimization_container, service, "requirement")

    use_slacks && (slack_vars = reserve_slacks(optimization_container, name))

    requirement = PSY.get_requirement(service)
    if parameters
        container = get_parameter_container(
            optimization_container,
            UpdateRef{SR}(SERVICE_REQUIREMENT, "requirement"),
        )
        param = get_parameter_array(container)
        multiplier = get_multiplier_array(container)
        for t in time_steps
            param[name, t] = add_parameter(optimization_container.JuMPmodel, ts_vector[t])
            multiplier[name, t] = 1.0
            if use_slacks
                resource_expression = sum(reserve_variable[:, t]) + slack_vars[t]
            else
                resource_expression = sum(reserve_variable[:, t])
            end
            mul = (requirement * multiplier[name, t])
            constraint[name, t] = JuMP.@constraint(
                optimization_container.JuMPmodel,
                resource_expression >= param[name, t] * mul
            )
        end
    else
        for t in time_steps
            constraint[name, t] = JuMP.@constraint(
                optimization_container.JuMPmodel,
                sum(reserve_variable[:, t]) >= ts_vector[t] * requirement
            )
        end
    end
    return
end

function service_requirement_constraint!(
    optimization_container::OptimizationContainer,
    service::SR,
    ::ServiceModel{SR, T},
) where {SR <: PSY.StaticReserve, T <: AbstractReservesFormulation}
    parameters = model_has_parameters(optimization_container)
    initial_time = model_initial_time(optimization_container)
    @debug initial_time
    time_steps = model_time_steps(optimization_container)
    name = PSY.get_name(service)
    constraint =
        get_constraint(optimization_container, make_constraint_name(REQUIREMENT, SR))
    reserve_variable = get_variable(optimization_container, name, SR)
    use_slacks = get_services_slack_variables(optimization_container.settings)

    use_slacks && (slack_vars = reserve_slacks(optimization_container, name))

    requirement = PSY.get_requirement(service)
    for t in time_steps
        resource_expression = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}()
        JuMP.add_to_expression!(resource_expression, sum(reserve_variable[:, t]))
        if use_slacks
            resource_expression += slack_vars[t]
        end
        constraint[name, t] = JuMP.@constraint(
            optimization_container.JuMPmodel,
            resource_expression >= requirement
        )
    end

    return
end

function cost_function!(
    optimization_container::OptimizationContainer,
    service::SR,
    ::ServiceModel{SR, T},
) where {SR <: PSY.Reserve, T <: AbstractReservesFormulation}
    reserve = get_variable(optimization_container, PSY.get_name(service), SR)
    for r in reserve
        JuMP.add_to_expression!(
            optimization_container.cost_function,
            r,
            DEFAULT_RESERVE_COST,
        )
    end
    return
end

function service_requirement_constraint!(
    optimization_container::OptimizationContainer,
    service::SR,
    ::ServiceModel{SR, StepwiseCostReserve},
) where {SR <: PSY.ReserveDemandCurve}
    initial_time = model_initial_time(optimization_container)
    @debug initial_time
    time_steps = model_time_steps(optimization_container)
    name = PSY.get_name(service)
    constraint =
        get_constraint(optimization_container, make_constraint_name(REQUIREMENT, SR))
    reserve_variable = get_variable(optimization_container, name, SR)
    requirement_variable = get_variable(optimization_container, SERVICE_REQUIREMENT, SR)

    for t in time_steps
        constraint[name, t] = JuMP.@constraint(
            optimization_container.JuMPmodel,
            sum(reserve_variable[:, t]) >= requirement_variable[name, t]
        )
    end

    return
end

function AddCostSpec(
    ::Type{T},
    ::Type{StepwiseCostReserve},
    optimization_container::OptimizationContainer,
) where {T <: PSY.Reserve}
    return AddCostSpec(;
        variable_type = ServiceRequirementVariable,
        component_type = T,
        has_status_variable = false,
        has_status_parameter = false,
        variable_cost = PSY.get_variable,
        start_up_cost = nothing,
        shut_down_cost = nothing,
        fixed_cost = nothing,
        sos_status = SOSStatusVariable.NO_VARIABLE,
    )
end

function add_to_cost!(
    optimization_container::OptimizationContainer,
    spec::AddCostSpec,
    service::SR,
    component_name::String,
) where {SR <: PSY.Reserve}
    time_steps = model_time_steps(optimization_container)
    use_forecast_data = model_uses_forecasts(optimization_container)
    if !use_forecast_data
        error("StepwiseCostReserve is only supported with forecast")
    end
    variable_cost_forecast =
        get_time_series(optimization_container, service, "variable_cost")
    variable_cost_forecast = map(PSY.VariableCost, variable_cost_forecast)
    for t in time_steps
        variable_cost!(
            optimization_container,
            spec,
            component_name,
            variable_cost_forecast[t],
            t,
        )
    end
    return
end

function cost_function!(
    optimization_container::OptimizationContainer,
    service::SR,
    model::ServiceModel{SR, StepwiseCostReserve},
) where {SR <: PSY.ReserveDemandCurve}
    spec = AddCostSpec(SR, model.formulation, optimization_container)
    @debug SR, spec
    add_to_cost!(optimization_container, spec, service, PSY.get_name(service))
    return
end

function modify_device_model!(
    devices_template::Dict{Symbol, DeviceModel},
    service_model::ServiceModel{<:PSY.Reserve, <:AbstractReservesFormulation},
    contributing_devices::Vector{<:PSY.Device},
)
    device_types = unique(typeof.(contributing_devices))
    for dt in device_types
        for (device_model_name, device_model) in devices_template
            # add message here when it exists
            get_component_type(device_model) != dt && continue
            service_model in device_model.services && continue
            push!(device_model.services, service_model)
        end
    end

    return
end

function include_service!(
    constraint_info::T,
    services,
    ::ServiceModel{SR, <:AbstractReservesFormulation},
) where {
    T <: Union{AbstractRangeConstraintInfo, AbstractRampConstraintInfo},
    SR <: PSY.Reserve{PSY.ReserveUp},
}
    for (ix, service) in enumerate(services)
        push!(
            constraint_info.additional_terms_ub,
            make_constraint_name(PSY.get_name(service), SR),
        )
    end
    return
end

function include_service!(
    constraint_info::T,
    services,
    ::ServiceModel{SR, <:AbstractReservesFormulation},
) where {
    T <: Union{AbstractRangeConstraintInfo, AbstractRampConstraintInfo},
    SR <: PSY.Reserve{PSY.ReserveDown},
}
    for (ix, service) in enumerate(services)
        push!(
            constraint_info.additional_terms_lb,
            make_constraint_name(PSY.get_name(service), SR),
        )
    end
    return
end

function add_device_services!(
    constraint_info::T,
    device::D,
    model::DeviceModel,
) where {
    T <: Union{AbstractRangeConstraintInfo, AbstractRampConstraintInfo},
    D <: PSY.Device,
}
    for service_model in get_services(model)
        if PSY.has_service(device, service_model.component_type)
            services = (
                s for s in PSY.get_services(device) if isa(s, service_model.component_type)
            )
            @assert !isempty(services)
            include_service!(constraint_info, services, service_model)
        end
    end
    return
end

function add_device_services!(
    constraint_data_in::AbstractRangeConstraintInfo,
    constraint_data_out::AbstractRangeConstraintInfo,
    device::D,
    model::DeviceModel{D, <:AbstractStorageFormulation},
) where {D <: PSY.Storage}
    for service_model in get_services(model)
        if PSY.has_service(device, service_model.component_type)
            services = (
                s for s in PSY.get_services(device) if isa(s, service_model.component_type)
            )
            @assert !isempty(services)
            if service_model.component_type <: PSY.Reserve{PSY.ReserveDown}
                for service in services
                    push!(
                        constraint_data_in.additional_terms_ub,
                        make_constraint_name(
                            PSY.get_name(service),
                            service_model.component_type,
                        ),
                    )
                end
            elseif service_model.component_type <: PSY.Reserve{PSY.ReserveUp}
                for service in services
                    push!(
                        constraint_data_out.additional_terms_ub,
                        make_constraint_name(
                            PSY.get_name(service),
                            service_model.component_type,
                        ),
                    )
                end
            end
        end
    end
    return
end


text:
numconnections(network::EmptyNetwork) = 0
isdirected(network::EmptyNetwork) = false
connectivity(network::EmptyNetwork) = 0
shortestpath(network::EmptyNetwork) = 0

adjVet(network::EmptyNetwork) = Int[]

function adjMat(network::EmptyNetwork; sparse::Bool=false)
	return sparse ?
		SparseArrays.spzeros(Bool, network.N, network.N) :
		BitArray(zeros(Bool, network.N, network.N))
end

#---------- Calculators ----------
calcConnectivity(network::EmptyNetwork, idx_node::Integer; dir_behaviour::Symbol) =
	dir_behaviour == :both ? (0,0) : 0