text:
## Sequence space defintions

# A Fun for SequenceSpace can be an iterator
iterate(::Fun{SequenceSpace}) = 1
iterate(f::Fun{SequenceSpace}, st) = f[st], st+1

getindex(f::Fun{SequenceSpace}, k::Integer) =
    k ≤ ncoefficients(f) ? f.coefficients[k] : zero(cfstype(f))
getindex(f::Fun{SequenceSpace},K::CartesianIndex{0}) = f[1]
getindex(f::Fun{SequenceSpace},K) = cfstype(f)[f[k] for k in K]

length(f::Fun{SequenceSpace}) = ℵ₀

dotu(f::Fun{SequenceSpace},g::Fun{SequenceSpace}) =
    mindotu(f.coefficients,g.coefficients)
dotu(f::Fun{SequenceSpace},g::AbstractVector) =
    mindotu(f.coefficients,g)
dotu(f::AbstractVector,g::Fun{SequenceSpace}) =
    mindotu(f,g.coefficients)

norm(f::Fun{SequenceSpace}) = norm(f.coefficients)
norm(f::Fun{SequenceSpace},k::Int) = norm(f.coefficients,k)
norm(f::Fun{SequenceSpace},k::Number) = norm(f.coefficients,k)

Fun(cfs::AbstractVector,S::SequenceSpace) = Fun(S,cfs)
coefficients(cfs::AbstractVector,::SequenceSpace) = cfs  # all vectors are convertible to SequenceSpace

## Constant space defintions

# setup conversions for spaces that contain constants
macro containsconstants(SP)
    esc(quote
        ApproxFunBase.union_rule(A::(ApproxFunBase.ConstantSpace),B::$SP) = B
        Base.promote_rule(A::Type{<:(ApproxFunBase.ConstantSpace)},B::Type{<:($SP)}) = B

        Base.promote_rule(::Type{ApproxFunBase.Fun{S,V,VV}},::Type{T}) where {T<:Number,S<:$SP,V,VV} =
        ApproxFunBase.VFun{S,promote_type(V,T)}
        Base.promote_rule(::Type{ApproxFunBase.Fun{S}},::Type{T}) where {T<:Number,S<:$SP} = ApproxFunBase.VFun{S,T}
        Base.promote_rule(::Type{ApproxFunBase.Fun{S,V,VV}},
        ::Type{Fun{ApproxFunBase.ConstantSpace{ApproxFunBase.AnyDomain},T,VT}}) where {T,S<:$SP,V,VV,VT} =
            ApproxFunBase.VFun{S,promote_type(V,T)}
    end)
end

Fun(c::Number) = Fun(ConstantSpace(typeof(c)),[c])
Fun(c::Number,d::ConstantSpace) = Fun(d,[c])

dimension(::ConstantSpace) = 1

#TODO: Change
setdomain(f::Fun{CS},d::Domain) where {CS<:AnyDomain} = Number(f)*ones(d)

canonicalspace(C::ConstantSpace) = C
spacescompatible(a::ConstantSpace,b::ConstantSpace)=domainscompatible(a,b)

ones(S::ConstantSpace) = Fun(S,fill(1.0,1))
ones(S::Union{AnyDomain,UnsetSpace}) = ones(ConstantSpace())
zeros(S::AnyDomain) = zero(ConstantSpace())
zero(S::UnsetSpace) = zero(ConstantSpace())
evaluate(f::AbstractVector,::ConstantSpace,x...)=f[1]
evaluate(f::AbstractVector,::ZeroSpace,x...)=zero(eltype(f))

convert(::Type{T}, f::Fun{CS}) where {CS<:ConstantSpace,T<:Number} =
    convert(T, f.coefficients[1])

Number(f::Fun) = convert(Number, f)

# promoting numbers to Fun
# override promote_rule if the space type can represent constants
Base.promote_rule(::Type{Fun{CS}},::Type{T}) where {CS<:ConstantSpace,T<:Number} = Fun{CS,T}
Base.promote_rule(::Type{Fun{CS,V}},::Type{T}) where {CS<:ConstantSpace,T<:Number,V} =
    Fun{CS,promote_type(T,V)}
Base.promote_rule(::Type{IF},::Type{T}) where {T<:Number,IF<:Fun} = Fun

# we know multiplication by constants preserves types
Base.promote_op(::typeof(*),::Type{Fun{CS,T,VT}},::Type{F}) where {CS<:ConstantSpace,T,VT,F<:Fun} =
    promote_op(*,T,F)
Base.promote_op(::typeof(*),::Type{F},::Type{Fun{CS,T,VT}}) where {CS<:ConstantSpace,T,VT,F<:Fun} =
    promote_op(*,F,T)

Base.promote_op(::typeof(LinearAlgebra.matprod),::Type{Fun{S1,T1,VT1}},::Type{Fun{S2,T2,VT2}}) where {S1<:ConstantSpace,T1,VT1,S2<:ConstantSpace,T2,VT2} =
            VFun{promote_type(S1,S2),promote_type(T1,T2)}
Base.promote_op(::typeof(LinearAlgebra.matprod),::Type{Fun{S1,T1,VT1}},::Type{Fun{S2,T2,VT2}}) where {S1<:ConstantSpace,T1,VT1,S2,T2,VT2} =
            VFun{S2,promote_type(T1,T2)}
Base.promote_op(::typeof(LinearAlgebra.matprod),::Type{Fun{S1,T1,VT1}},::Type{Fun{S2,T2,VT2}}) where {S1,T1,VT1,S2<:ConstantSpace,T2,VT2} =
            VFun{S1,promote_type(T1,T2)}

# When the union of A and B is a ConstantSpace, then it contains a one
conversion_rule(A::ConstantSpace,B::UnsetSpace)=NoSpace()
conversion_rule(A::ConstantSpace,B::Space)=(union_rule(A,B)==B||union_rule(B,A)==B) ? A : NoSpace()

conversion_rule(A::ZeroSpace,B::Space) = A
maxspace_rule(A::ZeroSpace,B::Space) = B

Conversion(A::ZeroSpace,B::ZeroSpace) = ConversionWrapper(ZeroOperator(A,B))
Conversion(A::ZeroSpace,B::Space) = ConversionWrapper(ZeroOperator(A,B))

# TODO: this seems like it needs more thought
union_rule(A::ConstantSpace,B::Space) = ConstantSpace(domain(B))⊕B

## Special Multiplication and Conversion for constantspace

#  TODO: this is a special work around but really we want it to be blocks
Conversion(a::ConstantSpace,b::Space{D}) where {D<:EuclideanDomain{2}} = ConcreteConversion{typeof(a),typeof(b),
        promote_type(real(prectype(a)),real(prectype(b)))}(a,b)

Conversion(a::ConstantSpace,b::Space) = ConcreteConversion(a,b)
bandwidths(C::ConcreteConversion{CS,S}) where {CS<:ConstantSpace,S<:Space} =
    ncoefficients(ones(rangespace(C)))-1,0
function getindex(C::ConcreteConversion{CS,S,T},k::Integer,j::Integer) where {CS<:ConstantSpace,S<:Space,T}
    if j != 1
        throw(BoundsError())
    end
    on=ones(rangespace(C))
    k ≤ ncoefficients(on) ? convert(T,on.coefficients[k]) : zero(T)
end

coefficients(f::AbstractVector,sp::ConstantSpace{Segment{Vec{2,TT}}},
             ts::TensorSpace{SV,DD}) where {TT,SV,DD<:EuclideanDomain{2}} =
    f[1]*ones(ts).coefficients
coefficients(f::AbstractVector,sp::ConstantSpace{<:Domain{<:Number}},
             ts::TensorSpace{SV,DD}) where {TT,SV,DD<:EuclideanDomain{2}} =
    f[1]*ones(ts).coefficients
coefficients(f::AbstractVector, sp::ConstantSpace{<:Domain{<:Number}}, ts::Space) =
    f[1]*ones(ts).coefficients
coefficients(f::AbstractVector, sp::ConstantSpace, ts::Space) =
    f[1]*ones(ts).coefficients

########
# Evaluation
########

#########
# Multiplication
#########

# this is identity operator, but we don't use MultiplicationWrapper to avoid
# ambiguity errors

defaultMultiplication(f::Fun{CS},b::ConstantSpace) where {CS<:ConstantSpace} =
    ConcreteMultiplication(f,b)
defaultMultiplication(f::Fun{CS},b::Space) where {CS<:ConstantSpace} =
    ConcreteMultiplication(f,b)
defaultMultiplication(f::Fun,b::ConstantSpace) = ConcreteMultiplication(f,b)

bandwidths(D::ConcreteMultiplication{CS1,CS2,T}) where {CS1<:ConstantSpace,CS2<:ConstantSpace,T} =
    0,0
getindex(D::ConcreteMultiplication{CS1,CS2,T},k::Integer,j::Integer) where {CS1<:ConstantSpace,CS2<:ConstantSpace,T} =
    k==j==1 ? convert(T,D.f.coefficients[1]) : one(T)

rangespace(D::ConcreteMultiplication{CS1,CS2,T}) where {CS1<:ConstantSpace,CS2<:ConstantSpace,T} =
    D.space

rangespace(D::ConcreteMultiplication{F,UnsetSpace,T}) where {F<:ConstantSpace,T} =
    UnsetSpace()
bandwidths(D::ConcreteMultiplication{F,UnsetSpace,T}) where {F<:ConstantSpace,T} =
    (ℵ₀,ℵ₀)
getindex(D::ConcreteMultiplication{F,UnsetSpace,T},k::Integer,j::Integer) where {F<:ConstantSpace,T} =
    error("No range space attached to Multiplication")

bandwidths(D::ConcreteMultiplication{CS,F,T}) where {CS<:ConstantSpace,F<:Space,T} = 0,0
blockbandwidths(D::ConcreteMultiplication{CS,F,T}) where {CS<:ConstantSpace,F<:Space,T} = 0,0
subblockbandwidths(D::ConcreteMultiplication{CS,F,T}) where {CS<:ConstantSpace,F<:Space,T} = 0,0
subblockbandwidths(D::ConcreteMultiplication{CS,F,T}, k) where {CS<:ConstantSpace,F<:Space,T} = 0
isbandedblockbanded(D::ConcreteMultiplication{CS,F,T}) where {CS<:ConstantSpace,F<:Space,T} = true
isblockbanded(D::ConcreteMultiplication{CS,F,T}) where {CS<:ConstantSpace,F<:Space,T} = true
getindex(D::ConcreteMultiplication{CS,F,T},k::Integer,j::Integer) where {CS<:ConstantSpace,F<:Space,T} =
    k==j ? convert(T, D.f) : zero(T)
rangespace(D::ConcreteMultiplication{CS,F,T}) where {CS<:ConstantSpace,F<:Space,T} = D.space

bandwidths(D::ConcreteMultiplication{F,CS,T}) where {CS<:ConstantSpace,F<:Space,T} = ncoefficients(D.f)-1,0
function getindex(D::ConcreteMultiplication{F,CS,T},k::Integer,j::Integer) where {CS<:ConstantSpace,F<:Space,T}
    k≤ncoefficients(D.f) && j==1 ? convert(T,D.f.coefficients[k]) : zero(T)
end
rangespace(D::ConcreteMultiplication{F,CS,T}) where {CS<:ConstantSpace,F<:Space,T} = space(D.f)

# functionals always map to Constant space
function promoterangespace(P::Operator,A::ConstantSpace,cur::ConstantSpace)
    @assert isafunctional(P)
    domain(A)==domain(cur) ? P : SpaceOperator(P,domainspace(P),A)
end

for op = (:*,:/)
    @eval $op(f::Fun,c::Fun{CS}) where {CS<:ConstantSpace} = f*convert(Number,c)
end

## Multivariate case
union_rule(a::TensorSpace,b::ConstantSpace{AnyDomain})=TensorSpace(map(sp->union(sp,b),a.spaces))
## Special spaces

function convert(::Type{T},f::Fun{TS}) where {TS<:TensorSpace,T<:Number}
    if all(sp->isa(sp,ConstantSpace),space(f).spaces)
        convert(T,f.coefficients[1])
    else
        error("Cannot convert $f to type $T")
    end
end

convert(::Type{T},
            f::Fun{TensorSpace{Tuple{CS1,CS2},DD,RR}}) where {CS1<:ConstantSpace,CS2<:ConstantSpace,T<:Number,DD,RR} =
    convert(T,f.coefficients[1])

isconstspace(sp::TensorSpace) = all(isconstspace,sp.spaces)

# Supports constants in operators
promoterangespace(M::ConcreteMultiplication{CS,UnsetSpace},
                             ps::UnsetSpace) where {CS<:ConstantSpace} = M
promoterangespace(M::ConcreteMultiplication{CS,UnsetSpace},
                             ps::Space) where {CS<:ConstantSpace} =
                        promoterangespace(Multiplication(M.f,space(M.f)),ps)

# Possible hack: we try uing constant space for [1 Operator()] \ z.
choosedomainspace(M::ConcreteMultiplication{D,UnsetSpace},sp::UnsetSpace) where {D<:ConstantSpace} = space(M.f)
choosedomainspace(M::ConcreteMultiplication{D,UnsetSpace},sp::Space) where {D<:ConstantSpace} = space(M.f)

Base.isfinite(f::Fun{CS}) where {CS<:ConstantSpace} = isfinite(Number(f))


text:
MP.promote_rule_constant(::Type{T}, ::Type{<:DMonomialLike{C}}) where {C, T} = Term{C, promote_type(T, Int)}
MP.promote_rule_constant(::Type{S}, ::Type{Term{C, T}}) where {S, C, T} = Term{C, promote_type(S, T)}
MP.promote_rule_constant(::Type{S}, ::Type{<:TermPoly{C, T}}) where {S, C, T} = Polynomial{C, promote_type(S, T)}

ext: jl


text:
###############################################################################
#
#   Integer.jl : Additional Nemo functionality for Julia Integer
#
###############################################################################

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

JuliaZZ = Integers{BigInt}()

zz = Integers{Int}()

parent(a::T) where T <: Integer = Integers{T}()

elem_type(::Type{Integers{T}}) where T <: Integer = T

parent_type(::Type{T}) where T <: Integer = Integers{T}

base_ring(a::Integer) = Union{}

isexact_type(::Type{T}) where T <: Integer = true

isdomain_type(::Type{T}) where T <: Integer = true

###############################################################################
#
#   Basic manipulation
#
###############################################################################

zero(::Integers{T}) where T <: Integer = T(0)

one(::Integers{T}) where T <: Integer = T(1)

if VERSION < v"0.7.0-DEV.1144"
isone(a::Integer) = a == 1
end

isunit(a::Integer) = a == 1 || a == -1

canonical_unit(a::T) where T <: Integer = a < 0 ? T(-1) : T(1)

###############################################################################
#
#   String I/O
#
###############################################################################

function show(io::IO, R::Integers)
   print(io, "Integers")
end

needs_parentheses(::Integer) = false

isnegative(a::Integer) = a < 0

show_minus_one(::Type{T}) where T <: Integer = false

###############################################################################
#
#   Modular arithmetic
#
###############################################################################

function powmod(a::T, b::Int, c::T) where T <: Integer
   b < 0 && throw(DomainError())
   # special cases
   if a == 0
      return T(0)
   elseif b == 0
      return T(1)
   else
      bit = ~((~UInt(0)) >> 1)
      while (UInt(bit) & b) == 0
         bit >>= 1
      end
      z = a
      bit >>= 1
      while bit != 0
         z = mod(z*z, c)
         if (UInt(bit) & b) != 0
            z = mod(z*a, c)
         end
         bit >>= 1
      end
      return z
   end
end

###############################################################################
#
#   Divides
#
###############################################################################

function divides(a::T, b::T) where T <: Integer
   q, r = divrem(a, b)
   return r == 0, q
end

###############################################################################
#
#   Exact division
#
###############################################################################

divexact(a::T, b::T) where T <: Integer = div(a, b)

divexact(a::BigInt, b::T) where T <: Integer = div(a, b)

###############################################################################
#
#   GCD
#
###############################################################################

function gcdinv(a::T, b::T) where T <: Integer
   g, s, t = gcdx(a, b)
   return g, s
end

###############################################################################
#
#   Exponential
#
###############################################################################

function exp(a::T) where T <: Integer
   a != 0 && throw(DomainError())
   return T(1)
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function zero!(a::Integer)
   return 0
end

function zero!(a::BigInt)
   ccall((:__gmpz_set_si, :libgmp), Void, (Ref{BigInt}, Int), a, 0)
   return a
end

function mul!(a::T, b::T, c::T) where T <: Integer
   return b*c
end

function mul!(a::BigInt, b::BigInt, c::BigInt)
   ccall((:__gmpz_mul, :libgmp), Void, (Ref{BigInt}, Ref{BigInt}, Ref{BigInt}), a, b, c)
   return a
end

function add!(a::T, b::T, c::T) where T <: Integer
   return b + c
end

function add!(a::BigInt, b::BigInt, c::BigInt)
   ccall((:__gmpz_add, :libgmp), Void, (Ref{BigInt}, Ref{BigInt}, Ref{BigInt}), a, b, c)
   return a
end

function addeq!(a::T, b::T) where T <: Integer
   return a + b
end

function addeq!(a::BigInt, b::BigInt)
   ccall((:__gmpz_add, :libgmp), Void, (Ref{BigInt}, Ref{BigInt}, Ref{BigInt}), a, a, b)
   return a
end

function addmul!(a::T, b::T, c::T, d::T) where T <: Integer
   return a + b*c
end

function addmul!(a::BigInt, b::BigInt, c::BigInt, d::BigInt)
   ccall((:__gmpz_addmul, :libgmp), Void, (Ref{BigInt}, Ref{BigInt}, Ref{BigInt}), a, b, c)
   return a
end

function addmul!(a::T, b::T, c::T) where T <: Integer # special case, no temporary required
   return a + b*c
end

function addmul!(a::BigInt, b::BigInt, c::BigInt) # special case, no temporary required
   ccall((:__gmpz_addmul, :libgmp), Void, (Ref{BigInt}, Ref{BigInt}, Ref{BigInt}), a, b, c)
   return a
end

###############################################################################
#
#   Random generation
#
###############################################################################

function rand(R::Integers, n::UnitRange{Int})
   return R(rand(n))
end

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (a::Integers{T})() where T <: Integer
   return T(0)
end

function (a::Integers{T})(b::Integer) where T <: Integer
   return T(b)
end


text:
using BrainFlow
using BrainFlowViz

nchannels = 8
nsamples = 256

# start a live plotting task
data_func = ()->rand(nsamples, nchannels)
t = @task BrainFlowViz.plot_data(data_func, nsamples, nchannels; y_lim = (0, 1))
schedule(t)

sleep(5)

### to stop the task:
schedule(t, InterruptException(), error=true)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Zenith_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Zenith")
JLLWrappers.@generate_main_file("Zenith", UUID("ab215d34-98c8-565c-96d8-6589c1ef91a8"))
end  # module Zenith_jll


text:
"""
    WeightNorm(layer; dim::Union{Tuple,Vector,Int,Nothing}=nothing, ps=Flux.params(layer))

Applies weight normalization to `ps` parameters in the given layer:

``w = g\\frac{v}{||v||}``

!!!note:
    This assumes that `layer` cannot contain non-trainable states.

Ref: [`https://arxiv.org/abs/1602.07868`](@ref)
"""
struct WeightNorm{Re,P,S,D} <: AbstractFluxLayer
    layer_re::Re
    parameters::P
    states::S
    dims::D
end

Flux.@functor WeightNorm (parameters, states)

# function Base.show(io::IO, wn::WeightNorm)
#     ps = Flux.params(wn)
#     p = update_parameters(wn)
#     l = wn.layer_re(p, wn.states)
#     print(io, "WeightNorm(")
#     print(io, l)
#     return print(") ", string(sum(length.(ps))), " Trainable Parameters")
# end

function WeightNorm(layer; dim::Union{Tuple,Vector,Int,Nothing}=nothing, ps=Flux.params(layer))
    if ps isa Array
        ps = Flux.params(ps)
    end
    dim = dim === nothing ? [ndims(p) for p in ps] : (dim isa Int ? [dim for _ in 1:length(ps)] : dim)

    p_, s_, layer_re = destructure_parameters_states(layer, ps.params)

    parameters = []
    for (i, p) in enumerate(ps)
        g_val = _norm(p, dim[i])
        v_val = copy(p)
        push!(parameters, (g_val, v_val))
    end

    return WeightNorm(layer_re, tuple(parameters...), s_, dim)
end

@inline compute_normed_weight(v, g, dim) = v .* (g ./ _norm(v, dim))

@inline function update_parameters(wn::WeightNorm)
    return vcat(ntuple(i -> vec(compute_normed_weight(wn.parameters[i][2], wn.parameters[i][1], wn.dims[i])),
                       length(wn.dims))...)
end

(wn::WeightNorm)(args...; kwargs...) = wn.layer_re(update_parameters(wn),wn.states)(args...; kwargs...)

@inline _norm(x; dims=Colon()) = sqrt.(sum(abs2, x; dims=dims))

# Compute norm over all dimensions except `except_dim`
@inline _norm(x::AbstractArray{T,N}, except_dim) where {T,N} = _norm(x; dims=filter(i -> i != except_dim, 1:N))


text:
using Test
using LinearAlgebra
using Random
using StaticArrays
using Rotations
# using BenchmarkTools

using RecoverPose

Random.seed!(0)

include("triangulate.jl")
include("p3p.jl")

"""
Compute essential matrix with R & t.
"""
function get_target_E(R, t)
    tc = SMatrix{3, 3, Float64}(
        0, t[3], -t[2],
        -t[3], 0, t[1],
        t[2], -t[1], 0,
    )
    E = tc * R
    E /= E[3, 3]
end

@testset "Five Point: Perfect Solution" begin
    res = 1024
    pp = res ÷ 2
    K = SMatrix{3, 3, Float64}(
        910, 0, 0,
        0, 910, 0,
        pp, pp, 1,
    )
    K_inv = inv(K)

    θ = π / 8
    R = RotXYZ(rand() * θ, rand() * θ, rand() * θ)
    t = SVector{3, Float64}(rand(), rand() * 2, rand() * 3)

    P_target = RecoverPose.get_transformation(R, t)
    E_target = get_target_E(R, t)
    P2 = K * P_target

    x1, x2 = SVector{2, Float64}[], SVector{2, Float64}[]

    n_points = 1000
    for i in 1:n_points
        p1 = SVector{3, Float64}(rand() * res, rand() * res, 1.0)
        p1h = SVector{4, Float64}((K_inv * p1)..., 1.0)
        p2 = P2 * p1h
        p2 = p2[1:2] .* (1 / p2[3])

        push!(x1, p1[[1, 2]])
        push!(x2, p2[[1, 2]])
    end

    cache = GEEV4x4Cache()
    pd1, pd2 = RecoverPose.pre_divide(x1, x2, K, K)

    n_inliers, model = five_point_ransac(x1, x2, pd1, pd2, K, K, cache)
    E_res, P_res, inliers, repr_error = model
    E = E_res ./ E_res[3, 3]
    @test repr_error < 1e-5
    # @test all(isapprox.(E, E_target; atol=1e-1))
    # @test all(isapprox.(P_res[1:3, 1:3], P_target[1:3, 1:3]; atol=1e-1))
    # @test norm(P_res[1:3, 4] .- P_target[1:3, 4]) < 1.0
    @test sum(inliers) == n_points
    @test n_inliers == n_points

    # _, model = essential_ransac(pd1, pd2, K[1, 1] + K[2, 2])
    # E_res, _, _ = model
    # E = E_res ./ E_res[3, 3]
    # n_inliers, P_res, inliers, e_error = recover_pose(E_res, x1, x2, K, K, cache)
    # @test e_error < 1e-5
    # @test all(isapprox.(E, E_target; atol=1e-1))
    # @test all(isapprox.(P_res[1:3, 1:3], P_target[1:3, 1:3]; atol=1e-1))
    # @test norm(P_res[1:3, 4] .- P_target[1:3, 4]) < 1.0
    # @test sum(inliers) == n_points
    # @test n_inliers == n_points

    # @btime five_point_ransac($x1, $x2, $pd1, $pd2, $K, $K, $cache)
    # @btime essential_ransac($x1, $x2, $pd1, $pd2, $K, $K)
end


text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

using StringEncodings
using StaticArrays
using Unitful
import Unitful: Length, Temperature, Quantity, Units

"""
Generates .jl files: a `mainfile` and several catalog files.

Each catalog file is a module representing a distinct glass catalog (e.g. NIKON, SCHOTT), generated from corresponding
AGF files in `sourcedir`. These are then included and exported in `mainfile`.
"""
function generate_jls(
    sourcenames::Vector{<:AbstractString}, mainfile::AbstractString, jldir::AbstractString, sourcedir::AbstractString
)
    id = 1
    catalogfiles = []
    glassnames = []

    # generate several catalog files (.jl)
    for catalogname in sourcenames
        # parse the sourcefile (.agf) into a catalog (native Julia dictionary)
        sourcefile = joinpath(sourcedir, "$(catalogname).agf")
        catalog = sourcefile_to_catalog(sourcefile)

        # parse the catalog into a module string and write it to a catalog file (.jl)
        id, modstring = catalog_to_modstring(id, catalogname, catalog)
        push!(catalogfiles, joinpath(jldir, "$(catalogname).jl"))
        open(catalogfiles[end], "w") do io
            write(io, modstring)
        end

        # track glass names for lookup purposes
        append!(glassnames, ["$(catalogname).$(glassname)" for glassname in keys(catalog)])
    end

    # generate the parent main file (.jl)
    agfstrings = [
        "export $(join(sourcenames, ", "))",
        "",
        ["include(raw\"$(catalogfile)\")" for catalogfile in catalogfiles]...,
        "",
        "const AGF_GLASS_NAMES = [$(join(repr.(glassnames), ", "))]",
        "const AGF_GLASSES = [$(join(glassnames, ", "))]",
        ""
    ]
    open(mainfile, "w") do io
        write(io, join(agfstrings, "\n"))
    end
end

"""
Parse a `sourcefile` (.agf) into a native dictionary, where each `kvp = (glassname, glassinfo)` is a glass.
"""
function sourcefile_to_catalog(source_file::AbstractString)
    catalog_dict = Dict{String,Dict{String}}()

    # check whether the file is UTF8 or UTF16 encoded
    is_utf8 = isvalid(readuntil(source_file, " "))
    fo = is_utf8 ? open(source_file, "r") : open(source_file, enc"UTF-16LE", "r")

    # store glass_name and transmission_data as persistent variables in between loops
    glass_name = ""
    transmission_data = nothing

    # read the file
    for line in readlines(fo)
        if strip(line) == "" || length(strip(line)) == 0 || startswith(line, "CC ") || startswith(line, "GC ")
            continue
        end
        if startswith(line, "NM ")
            nm = split(line)

            glass_name = make_valid_name(nm[2])
            transmission_data = Vector{SVector{3,Float64}}(undef, 0)

            catalog_dict[glass_name] = Dict{String,Any}()
            catalog_dict[glass_name]["raw_name"] = nm[2]
            catalog_dict[glass_name]["dispform"] = Int(parse(Float64, nm[3]))
            catalog_dict[glass_name]["Nd"] = parse(Float64, nm[5])
            catalog_dict[glass_name]["Vd"] = parse(Float64, nm[6])
            catalog_dict[glass_name]["exclude_sub"] = length(nm) < 7 ? 0 : Int(parse(Float64, nm[7]))
            catalog_dict[glass_name]["status"] = length(nm) < 8 ? 0 : Int(parse(Float64, nm[8]))
            catalog_dict[glass_name]["meltfreq"] = length(nm) < 9 || "-" ∈ nm ? 0 : Int(parse(Float64, nm[9]))
        elseif startswith(line, "ED ")
            ed = split(line)
            catalog_dict[glass_name]["TCE"] = parse(Float64, ed[2])
            catalog_dict[glass_name]["p"] = parse(Float64, ed[4])
            catalog_dict[glass_name]["ΔPgF"] = parse(Float64, ed[5])
            catalog_dict[glass_name]["ignore_thermal_exp"] = length(ed) < 6 ? 0 : Int(parse(Float64, ed[6]))
        elseif startswith(line, "CD ")
            cd = parse.(Float64, split(line)[2:end])
            catalog_dict[glass_name]["C1"] = get(cd, 1, NaN)
            catalog_dict[glass_name]["C2"] = get(cd, 2, NaN)
            catalog_dict[glass_name]["C3"] = get(cd, 3, NaN)
            catalog_dict[glass_name]["C4"] = get(cd, 4, NaN)
            catalog_dict[glass_name]["C5"] = get(cd, 5, NaN)
            catalog_dict[glass_name]["C6"] = get(cd, 6, NaN)
            catalog_dict[glass_name]["C7"] = get(cd, 7, NaN)
            catalog_dict[glass_name]["C8"] = get(cd, 8, NaN)
            catalog_dict[glass_name]["C9"] = get(cd, 9, NaN)
            catalog_dict[glass_name]["C10"] = get(cd, 10, NaN)
        elseif startswith(line, "TD ")
            td = parse.(Float64, split(line)[2:end])
            catalog_dict[glass_name]["D₀"] = get(td, 1, 0.0)
            catalog_dict[glass_name]["D₁"] = get(td, 2, 0.0)
            catalog_dict[glass_name]["D₂"] = get(td, 3, 0.0)
            catalog_dict[glass_name]["E₀"] = get(td, 4, 0.0)
            catalog_dict[glass_name]["E₁"] = get(td, 5, 0.0)
            catalog_dict[glass_name]["λₜₖ"] = get(td, 6, 0.0)
            catalog_dict[glass_name]["temp"] = get(td, 7, 20.0)
        elseif startswith(line, "OD ")
            od = stringlist_to_floatlist(split(line)[2:end])
            catalog_dict[glass_name]["relcost"] = get(od, 1, -1)
            catalog_dict[glass_name]["CR"] = get(od, 2, -1)
            catalog_dict[glass_name]["FR"] = get(od, 3, -1)
            catalog_dict[glass_name]["SR"] = get(od, 4, -1)
            catalog_dict[glass_name]["AR"] = get(od, 5, -1)
            catalog_dict[glass_name]["PR"] = get(od, 6, -1)
        elseif startswith(line, "LD ")
            ld = parse.(Float64, split(line)[2:end])
            catalog_dict[glass_name]["λmin"] = ld[1]
            catalog_dict[glass_name]["λmax"] = ld[2]
        elseif startswith(line, "IT ")
            it_row = parse.(Float64, split(line)[2:end])
            if length(it_row) == 3 && it_row[1] != 0.0
                entry = SVector(it_row[1], it_row[2], it_row[3])
                push!(transmission_data, entry)
            end
            catalog_dict[glass_name]["transmission"] = transmission_data
        end
    end
    return catalog_dict
end

"""
Make a valid Julia variable name from an arbitrary string.
"""
function make_valid_name(name::AbstractString)
    # remove invalid characters
    name = replace(name, "*" => "_STAR")
    name = replace(name, r"""[ ,.:;?!()&-]""" => "_")
    # cant have module names which are just numbers so add a _ to the start
    if tryparse(Int, "$(name[1])") !== nothing
        name = "_" * name
    end
    return name
end

function stringlist_to_floatlist(x::Vector{<:AbstractString})
    npts = length(x)
    if (npts == 0) || ((npts == 1) && (strip(x[1]) == "-"))
        return (repeat([-1.0], 10))
    end
    res = []
    for a in x
        if (strip(a) == "-")
            push!(res, -1.0)
        else
            try
                push!(res, parse(Float64, a))
            catch
                push!(res, NaN)
            end
        end
    end
    return (res)
end

"""
Convert a `catalog` dict into a `modstring` which can be written to a Julia source file.
"""
function catalog_to_modstring(start_id::Integer, catalogname::AbstractString, catalog::Dict{<:AbstractString})
    id = start_id
    isCI = haskey(ENV, "CI")

    modstrings = [
        "module $catalogname",
        "using ..GlassCat: Glass, GlassID, AGF",
        # "using StaticArrays: SVector",
        "export $(join(keys(catalog), ", "))",
        ""
    ]
    for (glassname, glassinfo) in catalog
        # skip docstrings for CI builds to avoid 'missing docstring' warnings in makedocs
        docstring = isCI ? "" : glassinfo_to_docstring(glassinfo, id, catalogname, glassname)
        argstring = glassinfo_to_argstring(glassinfo, id)
        append!(modstrings, [docstring, "const $glassname = Glass($argstring)", ""])
        id += 1
    end
    append!(modstrings, ["end #module", ""]) # last "" is for \n at EOF

    return id, join(modstrings, "\n")
end

"""
Convert a `glassinfo` dict into a `docstring` to be prepended to a `Glass` const.
"""
function glassinfo_to_docstring(
    glassinfo::Dict{<:AbstractString}, id::Integer, catalogname::AbstractString, glassname::AbstractString
)
    raw_name = glassinfo["raw_name"] == glassname ? "" : " ($(glassinfo["raw_name"]))"
    pad(str, padding=25) = rpad(str, padding)
    getinfo(key, default=0.0) = get(glassinfo, key, default)

    return join([
        "\"\"\"    $catalogname.$glassname$raw_name",
        "```",
        "$(pad("ID:"))AGF:$id",
        "$(pad("RI @ 587nm:"))$(getinfo("Nd"))",
        "$(pad("Abbe Number:"))$(getinfo("Vd"))",
        "$(pad("ΔPgF:"))$(getinfo("ΔPgF"))",
        "$(pad("TCE (÷1e-6):"))$(getinfo("TCE"))",
        "$(pad("Density:"))$(getinfo("p"))g/m³",
        "$(pad("Valid wavelengths:"))$(getinfo("λmin"))μm to $(getinfo("λmax"))μm",
        "$(pad("Reference Temp:"))$(getinfo("temp", 20.0))°C",
        "```",
        "\"\"\""
    ], "\n")
end

"""
Convert a `glassinfo` dict into an `argstring` to be passed into a `Glass` constructor.
"""
function glassinfo_to_argstring(glassinfo::Dict{<:AbstractString}, id::Integer)
    argstrings = []
    for fn in string.(fieldnames(Glass))
        if fn == "ID"
            push!(argstrings, "GlassID(AGF, $id)")
        elseif fn in ["D₀", "D₁", "D₂", "E₀", "E₁", "λₜₖ"]
            push!(argstrings, repr(get(glassinfo, fn, 0.0)))
        elseif fn == "temp"
            push!(argstrings, repr(get(glassinfo, fn, 20.0)))
        elseif fn == "transmission"
            v = get(glassinfo, "transmission", nothing)
            if isnothing(v)
                push!(argstrings, repr(nothing))
            else
                str = join(["($(join(a, ", ")))" for a in v], ", ")
                push!(argstrings, "[$str]")
            end
        elseif fn == "transmissionN"
            continue
        else
            push!(argstrings, repr(get(glassinfo, fn, NaN)))
        end
    end
    return join(argstrings, ", ")
end


text:
# two arg clamp:
# takes a Type arg, and uses typemin and typemax to get the limits
# AND ALSO converts to this type since after clamping, this is safe

# since I'm effectively adding to Base here, get warnings after reloading
# see https://github.com/JuliaLang/julia/issues/7860#issuecomment-51340412

import Base.clamp
clamp{T}( x, ::Type{T} ) = clamp( x, typemin(T), typemax(T) ) % T

# copied and modified from base/math.jl
# this makes two arg clamp work with arrays
clamp{T,T2}(x::AbstractArray{T,1}, ::Type{T2} ) = [clamp(xx, T2 ) for xx in x]
clamp{T,T2}(x::AbstractArray{T,2}, ::Type{T2}) =
    [clamp(x[i,j], T2) for i in 1:size(x,1), j in 1:size(x,2)]  # fixme (iter): change to `eachindex` when #15459 is ready
clamp{T,T2}(x::AbstractArray{T}, ::Type{T2}) =
    reshape([clamp(xx, T2) for xx in x], size(x))

function clamp!{T,T2}(x::AbstractArray{T}, ::Type{T2})
    @inbounds for i in eachindex(x)
        x[i] = clamp(x[i], T2 )
    end
    x
end

import Base.clamp
# make two arg clamp more efficient for BInt RHS
@generated function Base.clamp{Tx<:Integer,R}( x::Tx, ::Type{BInt{R}} )
        l = typemin(BInt{R})
        u = typemax(BInt{R})
        Tw = default_int_type( range_union( R, Tx ) )
        # converting x to Tw makes this work for Tx being either
        # an ordinary Int, or a Bint
        :( UncheckedBInt{R}( Base.clamp( x % $Tw, $l , $u ) ) )
end
@generated function Base.clamp{Tx<:Integer,R,T}( x::Tx, ::Type{BInt{R,T}} )
        l = typemin(BInt{R,T})
        u = typemax(BInt{R,T})
        Tw = default_int_type( range_union( R, Tx ) )
        :( UncheckedBInt{R,T}( Base.clamp( x % $Tw, $l , $u ) ) )
end

# this works for arrays
# clamp( BInt{Int16}[ 10, 2000, 30 ], BInt{Int8} )
# ideally by adding Base.clamp( n, T ) and associated Array code


text:
## Exercise 9-3
# Write a function named avoids that takes a word and a string of forbidden letters, and that returns true if the word doesn’t use any of the forbidden letters.
println("Ans:\n[Part 1]")

function avoids(string, forbidden_letters)
    for letter in forbidden_letters
        if letter ∈ string
            return false
        end
    end
    return true
end

println(avoids("hello world", "hwe"))
println(avoids("printf", "xyz"))

# Modify your program to prompt the user to enter a string of forbidden letters and then print the number of words that don’t contain any of them. Can you find a combination of 5 forbidden letters that excludes the smallest number of words?
println("[Part 2]")

print("Please Enter forbidden letters: ")
forbidden_letters = strip(readline())

count = 0
for line in eachline("words.txt")
    if avoids(line, forbidden_letters)
        global count += 1
    end
end

println("There are total of $count letters which avoids the forbidden letters.")

println("End.")


text:
using Test
using TimeseriesSurrogates
ENV["GKSwstype"] = "100"

N = 1000
ts = cumsum(randn(N))
ts_nan = cumsum(randn(N))
ts_nan[1] = NaN
x = cos.(range(0, 20π, length = N)) .+ randn(N)*0.05

@testset "WLS" begin
    wts = WLS()
    wts = WLS(AAFT())
    s = surrogate(x, wts)

    @test length(s) == length(x)
    @test all([s[i] ∈ x for i = 1:N])
end

@testset "Periodic" begin
	pp = PseudoPeriodic(3, 25, 0.05)
	s = surrogate(x, pp)
	@test length(s) == length(ts)
	@test all(s[i] ∈ x for i in 1:N)
	# Perhaps a more advanced test, e.g. that both components have Fourier peak at
	# the same frequency, should be considered.

	#TODO: Test for noiseradius
end

@testset "BlockShuffle" begin
    bs1 = BlockShuffle()
    bs2 = BlockShuffle(4)
    s1 = surrogate(x, bs1)
    s2 = surrogate(x, bs2)

    @test length(s1) == length(x)
    @test length(s2) == length(x)
    @test all([s1[i] ∈ x for i = 1:N])
    @test all([s2[i] ∈ x for i = 1:N])
end

@testset "RandomShuffle" begin
    rs = RandomShuffle()
    s = surrogate(x, rs)

    @test length(s) == length(x)
    @test all([s[i] ∈ x for i = 1:N])
end

@testset "AutoRegressive" begin
	y = TimeseriesSurrogates.AR1(2000, 0.1, 0.5)
    sg = surrogenerator(y, AutoRegressive(1))
	@test 0.4 ≤ abs(sg.init.φ[1]) ≤ 0.6
	s = sg()
    @test length(s) == length(y)
end

@testset "AAFT" begin
    aaft = AAFT()
    s = surrogate(x, aaft)

    @test length(s) == length(x)
    @test all([s[i] ∈ x for i = 1:N])
end

@testset "IAAFT" begin
    iaaft = IAAFT()
    s = surrogate(x, iaaft)

    @test length(s) == length(x)
    @test all([s[i] ∈ x for i = 1:N])
end

@testset "TFTS" begin
    method_preserve_lofreq = TFTS(0.05)
    method_preserve_hifreq = TFTS(-0.05)

    s = surrogate(x, method_preserve_lofreq)
    @test length(s) == length(x)

    s = surrogate(x, method_preserve_hifreq)
    @test length(s) == length(x)

    @test_throws ArgumentError TFTS(0)
end

@testset "TAAFT" begin
    method_preserve_lofreq = TAAFT(0.05)
    method_preserve_hifreq = TAAFT(-0.05)

    s = surrogate(x, method_preserve_lofreq)
    @test length(s) == length(x)
    @test all([s[i] ∈ x for i = 1:N])

    s = surrogate(x, method_preserve_hifreq)
    @test length(s) == length(x)
    @test all([s[i] ∈ x for i = 1:N])

    @test_throws ArgumentError TAAFT(0)
end

@testset "RandomFourier" begin
    @testset "random phases" begin
        phases = true
        rf = RandomFourier(phases)
        s = surrogate(x, rf)

        @test length(s) == length(x)
    end

    @testset "random amplitudes" begin
        phases = false
        rf = RandomFourier(phases)
        s = surrogate(x, rf)

        @test length(s) == length(x)
    end
end

@testset "Circ/Cycle shuffle" begin
	x = random_cycles()
	s = surrogate(x, CycleShuffle())
	for a in s
		@test a ∈ x
	end
	s = surrogate(x, CircShift(1:length(x)))
	for a in s
		@test a ∈ x
	end
end

using DelayEmbeddings
@testset "ShufleDims" begin
	X = Dataset(rand(100, 3))
	Y = surrogate(X, ShuffleDimensions())
	for i in 1:100
		@test sort(X[i]) == sort(Y[i])
	end
end

#=
@testset "IAAFT" begin
    # With pre-planning
    method = IAAFT(ts)
    surr = surrogate(ts, method)
    @test length(ts) == length(surr)
    @test all(sort(ts) .== sort(surr))

    # Without pre-planning
    method = IAAFT()
    surr = surrogate(ts, method)
    @test length(ts) == length(surr)
    @test all(sort(ts) .== sort(surr))
end
=#


text:
convertStringToVisualMaterial(visualMaterial::Shapes.VisualMaterial) = visualMaterial

convertStringToVisualMaterial(visualMaterial::AbstractString) = Shapes.visualMaterialPalette[visualMaterial]

"""
    Visual(;
        shape = Sphere(diameter = 0.5),
        visualMaterial = VisualMaterial(color = "Green"))

Generate a new visual shape with a visual material. It is only for visualizing, without physical behavior. In this case it is a greenish sphere.

# Arguments
- `shape`: it is possible to use all [Shapes](@ref) from both sets visual and solid, e.g. [`Sphere`](@ref), [`Box`](@ref), [`CoordinateSystem`](@ref). The default value is `nothing`, this means it will not be visualized.
- `visualMaterial`: use a predefined [Visual material](@ref) defined in palettes/visualMaterials.json, or define it with `VisualMaterial(color = "Blue")`. The default value is `VisualMaterial()`.
"""
struct Visual <: Modia3D.AbstractObject3DFeature
    shape::Union{Modia3D.AbstractShape,Nothing}
    visualMaterial::Shapes.VisualMaterial

    function Visual(;
        shape::Union{Modia3D.AbstractShape,Nothing} = nothing,
        visualMaterial::Union{Shapes.VisualMaterial,AbstractString,Nothing} = Shapes.VisualMaterial())
        new(shape, convertStringToVisualMaterial(visualMaterial) )
    end
end


text:
using Parameters
import Base: length

abstract type AbstractSpecies end
export AbstractSpecies

struct EmptySpecies <: AbstractSpecies end

@with_kw struct Species{T<:AbstractThermo,N<:AbstractDiffusivity,N1<:AbstractHenryLawConstant,N2<:AbstractLiquidVolumetricMassTransferCoefficient} <: AbstractSpecies
    name::String
    index::Integer
    inchi::String = ""
    smiles::String = ""
    adjlist::String = ""
    thermo::T
    atomnums::Dict{String,Int64} = Dict()
    bondnum::Int64=-1
    diffusion::N = EmptyDiffusivity()
    radius::Float64 = 0.0
    radicalelectrons::Int64 = -100
    molecularweight::Float64 = 0.0
    henrylawconstant::N1 = EmptyHenryLawConstant()
    liquidvolumetricmasstransfercoefficient::N2 = EmptyLiquidVolumetricMassTransferCoefficient()
end
export Species

length(r::T) where {T<:AbstractSpecies}= 1
export length


text:
using Test
using MPI

MPI.Init()
comm = MPI.COMM_WORLD
nnodes = MPI.Comm_size(comm)
ndims = 3
reorder = 1
periods = [0,1,0]
dims = [0,0,0]
MPI.Dims_create!(nnodes, dims)
comm_cart = MPI.Cart_create(comm, dims, periods, reorder)
coords = MPI.Cart_coords(comm_cart)
disp = 1

for i in 0:2
    neighbors = MPI.Cart_shift(comm_cart, i, disp)
    @test all(rank ->  0 <= rank < nnodes || rank == MPI.MPI_PROC_NULL, neighbors)
end

comm_cart = nothing
GC.gc()

MPI.Finalize()
@test MPI.Finalized()


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct AndersonFouadMachine <: Machine
        R::Float64
        Xd::Float64
        Xq::Float64
        Xd_p::Float64
        Xq_p::Float64
        Xd_pp::Float64
        Xq_pp::Float64
        Td0_p::Float64
        Tq0_p::Float64
        Td0_pp::Float64
        Tq0_pp::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        internal::InfrastructureSystemsInternal
    end

Parameters of 6-states synchronous machine: Anderson-Fouad model

# Arguments
- `R::Float64`: Resistance after EMF in machine per unit, validation range: `(0, nothing)`
- `Xd::Float64`: Reactance after EMF in d-axis per unit, validation range: `(0, nothing)`
- `Xq::Float64`: Reactance after EMF in q-axis per unit, validation range: `(0, nothing)`
- `Xd_p::Float64`: Transient reactance after EMF in d-axis per unit, validation range: `(0, nothing)`
- `Xq_p::Float64`: Transient reactance after EMF in q-axis per unit, validation range: `(0, nothing)`
- `Xd_pp::Float64`: Sub-Transient reactance after EMF in d-axis per unit, validation range: `(0, nothing)`
- `Xq_pp::Float64`: Sub-Transient reactance after EMF in q-axis per unit, validation range: `(0, nothing)`
- `Td0_p::Float64`: Time constant of transient d-axis voltage, validation range: `(0, nothing)`
- `Tq0_p::Float64`: Time constant of transient q-axis voltage, validation range: `(0, nothing)`
- `Td0_pp::Float64`: Time constant of sub-transient d-axis voltage, validation range: `(0, nothing)`
- `Tq0_pp::Float64`: Time constant of sub-transient q-axis voltage, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states are:
	ψq: q-axis stator flux,
	ψd: d-axis stator flux,
	eq_p: q-axis transient voltage,
	ed_p: d-axis transient voltage,
	eq_pp: q-axis subtransient voltage,
	ed_pp: d-axis subtransient voltage
- `n_states::Int`: AndersonFouadMachine has 6 states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct AndersonFouadMachine <: Machine
    "Resistance after EMF in machine per unit"
    R::Float64
    "Reactance after EMF in d-axis per unit"
    Xd::Float64
    "Reactance after EMF in q-axis per unit"
    Xq::Float64
    "Transient reactance after EMF in d-axis per unit"
    Xd_p::Float64
    "Transient reactance after EMF in q-axis per unit"
    Xq_p::Float64
    "Sub-Transient reactance after EMF in d-axis per unit"
    Xd_pp::Float64
    "Sub-Transient reactance after EMF in q-axis per unit"
    Xq_pp::Float64
    "Time constant of transient d-axis voltage"
    Td0_p::Float64
    "Time constant of transient q-axis voltage"
    Tq0_p::Float64
    "Time constant of sub-transient d-axis voltage"
    Td0_pp::Float64
    "Time constant of sub-transient q-axis voltage"
    Tq0_pp::Float64
    ext::Dict{String, Any}
    "The states are:
	ψq: q-axis stator flux,
	ψd: d-axis stator flux,
	eq_p: q-axis transient voltage,
	ed_p: d-axis transient voltage,
	eq_pp: q-axis subtransient voltage,
	ed_pp: d-axis subtransient voltage"
    states::Vector{Symbol}
    "AndersonFouadMachine has 6 states"
    n_states::Int
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function AndersonFouadMachine(R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, ext=Dict{String, Any}(), )
    AndersonFouadMachine(R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, ext, [:ψq, :ψd, :eq_p, :ed_p, :eq_pp, :ed_pp], 6, InfrastructureSystemsInternal(), )
end

function AndersonFouadMachine(; R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, ext=Dict{String, Any}(), states=[:ψq, :ψd, :eq_p, :ed_p, :eq_pp, :ed_pp], n_states=6, internal=InfrastructureSystemsInternal(), )
    AndersonFouadMachine(R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, ext, states, n_states, internal, )
end

# Constructor for demo purposes; non-functional.
function AndersonFouadMachine(::Nothing)
    AndersonFouadMachine(;
        R=0,
        Xd=0,
        Xq=0,
        Xd_p=0,
        Xq_p=0,
        Xd_pp=0,
        Xq_pp=0,
        Td0_p=0,
        Tq0_p=0,
        Td0_pp=0,
        Tq0_pp=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`AndersonFouadMachine`](@ref) `R`."""
get_R(value::AndersonFouadMachine) = value.R
"""Get [`AndersonFouadMachine`](@ref) `Xd`."""
get_Xd(value::AndersonFouadMachine) = value.Xd
"""Get [`AndersonFouadMachine`](@ref) `Xq`."""
get_Xq(value::AndersonFouadMachine) = value.Xq
"""Get [`AndersonFouadMachine`](@ref) `Xd_p`."""
get_Xd_p(value::AndersonFouadMachine) = value.Xd_p
"""Get [`AndersonFouadMachine`](@ref) `Xq_p`."""
get_Xq_p(value::AndersonFouadMachine) = value.Xq_p
"""Get [`AndersonFouadMachine`](@ref) `Xd_pp`."""
get_Xd_pp(value::AndersonFouadMachine) = value.Xd_pp
"""Get [`AndersonFouadMachine`](@ref) `Xq_pp`."""
get_Xq_pp(value::AndersonFouadMachine) = value.Xq_pp
"""Get [`AndersonFouadMachine`](@ref) `Td0_p`."""
get_Td0_p(value::AndersonFouadMachine) = value.Td0_p
"""Get [`AndersonFouadMachine`](@ref) `Tq0_p`."""
get_Tq0_p(value::AndersonFouadMachine) = value.Tq0_p
"""Get [`AndersonFouadMachine`](@ref) `Td0_pp`."""
get_Td0_pp(value::AndersonFouadMachine) = value.Td0_pp
"""Get [`AndersonFouadMachine`](@ref) `Tq0_pp`."""
get_Tq0_pp(value::AndersonFouadMachine) = value.Tq0_pp
"""Get [`AndersonFouadMachine`](@ref) `ext`."""
get_ext(value::AndersonFouadMachine) = value.ext
"""Get [`AndersonFouadMachine`](@ref) `states`."""
get_states(value::AndersonFouadMachine) = value.states
"""Get [`AndersonFouadMachine`](@ref) `n_states`."""
get_n_states(value::AndersonFouadMachine) = value.n_states
"""Get [`AndersonFouadMachine`](@ref) `internal`."""
get_internal(value::AndersonFouadMachine) = value.internal

"""Set [`AndersonFouadMachine`](@ref) `R`."""
set_R!(value::AndersonFouadMachine, val) = value.R = val
"""Set [`AndersonFouadMachine`](@ref) `Xd`."""
set_Xd!(value::AndersonFouadMachine, val) = value.Xd = val
"""Set [`AndersonFouadMachine`](@ref) `Xq`."""
set_Xq!(value::AndersonFouadMachine, val) = value.Xq = val
"""Set [`AndersonFouadMachine`](@ref) `Xd_p`."""
set_Xd_p!(value::AndersonFouadMachine, val) = value.Xd_p = val
"""Set [`AndersonFouadMachine`](@ref) `Xq_p`."""
set_Xq_p!(value::AndersonFouadMachine, val) = value.Xq_p = val
"""Set [`AndersonFouadMachine`](@ref) `Xd_pp`."""
set_Xd_pp!(value::AndersonFouadMachine, val) = value.Xd_pp = val
"""Set [`AndersonFouadMachine`](@ref) `Xq_pp`."""
set_Xq_pp!(value::AndersonFouadMachine, val) = value.Xq_pp = val
"""Set [`AndersonFouadMachine`](@ref) `Td0_p`."""
set_Td0_p!(value::AndersonFouadMachine, val) = value.Td0_p = val
"""Set [`AndersonFouadMachine`](@ref) `Tq0_p`."""
set_Tq0_p!(value::AndersonFouadMachine, val) = value.Tq0_p = val
"""Set [`AndersonFouadMachine`](@ref) `Td0_pp`."""
set_Td0_pp!(value::AndersonFouadMachine, val) = value.Td0_pp = val
"""Set [`AndersonFouadMachine`](@ref) `Tq0_pp`."""
set_Tq0_pp!(value::AndersonFouadMachine, val) = value.Tq0_pp = val
"""Set [`AndersonFouadMachine`](@ref) `ext`."""
set_ext!(value::AndersonFouadMachine, val) = value.ext = val


text:
# # A model with homeowners and renters

using Test #jl
using Aiyagari
using QuantEcon, Parameters
using Plots #md

PKG_HOME = joinpath(dirname(pathof(Aiyagari)), "..")

# ## Setting up the state space

# Exogenous states (incomes)

z_grid = [0.5; 1.0; 1.5]
z_prob = [0.7 0.15 0.15;
          0.2 0.6 0.2;
          0.15 0.15 0.7]
z_MC = MarkovChain(z_prob, z_grid, :z)

exo = ExogenousStateSpace([z_MC])

# ## Define utility function for the case of two goods (consumption good and housing)

u(c; γ) = c^(1-γ) / (1-γ)

function u(c,h; ξ=0.8159, ρ=map(s -> (s-1)/s, 0.13), γ=2.0)
  C = (ξ * h^ρ + (1-ξ) * c^ρ)^(1/ρ)
  u(C, γ=γ)
end

include(joinpath(PKG_HOME, "test", "housing-simple-nlopt.jl"))
include(joinpath(PKG_HOME, "test", "renting-nlopt.jl"))

# ## Aggregate state

mutable struct HousingAS{T1,T2,T3,T4} <: AggregateState
  r::T1
  p::T2 # house price
  ρ::T3 # rent
  dist::T4 # the distribution over idiosynchratic states
end

function HousingAS(r, p, w_grid, exo, param; ρ=p * (param.δ + r))
  dist_proto = zeros((length(w_grid), length(exo)))
  HousingAS(r, p, ρ, dist_proto)
end

# # Case 1: Forever home owners

w_grid_own = LinRange(0.0, 3.0, 40)
endo_own = EndogenousStateSpace((w=w_grid_own,))

r_own = 0.15
p_own = 0.9
param_own  = (β = 0.9, θ = 0.9, δ = 0.1, h_thres = eps())
agg_state_own = HousingAS(r_own, p_own, w_grid_own, exo, param_own)

@unpack policies_full, val = solve_bellman(endo_own, exo, agg_state_own, param_own, Owner(), rtol=√eps())

dist = stationary_distribution(endo_own, exo, policies_full.w_next)

using DelimitedFiles #jl

#

plt_w = plot(w_grid_own, policies_full.w_next, title="wealth next period", xlab="wealth") #md
plt_m_zoomed = plot(w_grid_own[1:15], policies_full.m[1:15,:], title="mortgage (zoomed)", xlab="wealth") #md

plt_m = plot(w_grid_own, policies_full.m, title="mortgage", xlab="wealth") #md

plt_dist = plot(w_grid_own, dist, xlab="wealth", title="stationary distribution") #md

plot(plt_w, plt_dist, plt_m, plt_m_zoomed, legend=false) #md

# using DelimitedFiles #jl
#writedlm(joinpath(PKG_HOME, "test/matrices", "housing-simple-value.txt"), val) #src
#writedlm(joinpath(PKG_HOME, "test/matrices", "housing-simple-dist.txt"), dist) #src

@testset "regression test housing-simple" begin #jl
  value_test = readdlm(joinpath(PKG_HOME, "test/matrices", "housing-simple-value.txt")) #jl
  dist_test = readdlm(joinpath(PKG_HOME, "test/matrices", "housing-simple-dist.txt")) #jl
  @test all(val .== value_test) #jl
  @test maximum(abs, dist .- dist_test) < 1e-12 #jl
end #jl

# ## Equilibrium

# First, let's compute excess demand

function excess_demand(r, p)
  agg_state_own = HousingAS(r, p, w_grid_own, exo, param_own)

  @unpack policies_full, val = solve_bellman(endo_own, exo, agg_state_own, param_own, Owner(), rtol=√eps())

  dist = stationary_distribution(endo_own, exo, policies_full.w_next)

  (m = sum(dist .* policies_full.m), h=sum(dist .* policies_full.h))
end

excess_demand(r_own, p_own)

# TODO: Housing supply

# # Own big, rent small

# Now we want to give households to choice whether to buy or to rent

function combined_policies(policies_full, owner, own, rent=own; f_own =identity, f_rent=identity)
  f_own.(getproperty(policies_full[1], own)) .* owner .+ f_rent.(getproperty(policies_full[2], rent)) .* .!(owner)
end

function excess_demand(r, p, endo, exo, param, hh; maxiter_bellman=200)
  agg_state = HousingAS(r, p, endo.grids.w, exo, param, ρ= p * (param[1].δ + r))

  out = solve_bellman(endo, exo, agg_state, param, hh, maxiter=300, rtol=√eps())

  @unpack val, policy, policies_full, policy_hh = out
  owner = policy_hh[1] .== 1

  w_next_all = combined_policies(policies_full, owner, :w_next)

  dist = stationary_distribution(endo, exo, w_next_all)
  a_all = combined_policies(policies_full, owner, :m, :w_next, f_own = x -> -x)
  h_all = combined_policies(policies_full, owner, :h)

  (a = sum(dist .* a_all), h=sum(dist .* h_all), dist=dist, out...)
end

r = 0.01
p = 1.7
w_grid = LinRange(0.0, 4.0, 40)
endo = EndogenousStateSpace((w=w_grid,))

param_both = (β = 0.93, θ = 0.9, δ = 0.1, h_thres = 1.2)
param = [param_both, param_both]

@unpack val, policy, policies_full, policy_hh, a, h, dist = excess_demand(r, p, endo, exo, param, OwnOrRent())

@show r => a, p => h
owner = policy_hh[1] .== 1
value = val[1]

plt_own = plot(w_grid, owner, title="Who owns?") #md

w_next_all = combined_policies(policies_full, owner, :w_next)
h_all = combined_policies(policies_full, owner, :h)
c_all = combined_policies(policies_full, owner, :c)
a_all = combined_policies(policies_full, owner, :m, :w_next, f_own = x -> -x)

plt_h = plot(w_grid, h_all, legend=:false, title="House size", markerstrokewidth=0, xlab="wealth") #md

plt_a = plot(w_grid, a_all, title = "asset/mortgage") #md

plt_w = plot(w_grid, w_next_all, title = "wealth next period") #md

plt_dist = plot(w_grid, dist, title = "stationary distribution") #md

plot(plt_h, plt_a, plt_w, plt_dist, legend=false) #md
#- #md

#writedlm(joinpath(PKG_HOME, "test/matrices", "own-rent-value.txt"), value) #src
#writedlm(joinpath(PKG_HOME, "test/matrices", "own-rent-dist.txt"), dist) #src

@testset "own-rent" begin #jl
  value_test = readdlm(joinpath(PKG_HOME, "test/matrices", "own-rent-value.txt")) #jl
  dist_test = readdlm(joinpath(PKG_HOME, "test/matrices", "own-rent-dist.txt")) #jl
  @test all(value .≈ value_test) #jl
  @test maximum(abs, dist .- dist_test) < 1e-12 #jl
end #jl


text:
using Test, ImageFeatures, Images

@testset "Types" begin
    img = zeros(10, 10)
    ids = map(CartesianIndex{2}, [(1, 1), (3, 4), (4, 6), (7, 5), (5, 3)])

    img[ids] .= 1

    keypoints = Keypoints(img)
    @test sort(keypoints) == sort(ids)
    @test Keypoint(1, 2) == CartesianIndex{2}(1, 2)

    features = Features(img)
    keypoints = Keypoints(features)
    @test sort(keypoints) == sort(ids)
end

@testset "Keypoint Matching" begin
    @test hamming_distance([false], [false]) == 0
    @test hamming_distance([false, true], [false, false]) == 0.5
    @test hamming_distance([false, false, true], [false, false, true]) == 0
    @test hamming_distance([false, false, true], [false, true, false]) == 2 / 3
    @test hamming_distance([false], [true]) == 1.0

    k1 = map(CartesianIndex{2}, [(1, 1), (2, 2), (3, 3), (4, 4)])
    k2 = map(CartesianIndex{2}, [(5, 5)])
    d1 = [[false, false, false],
          [false, false, true],
          [false, true, false],
          [true, true, true]]
    d2 = [[false, true, false]]
    matches = match_keypoints(k1, k2, d1, d2)
    #expected_matches = [[CartesianIndex(3, 3), CartesianIndex(5, 5)]]
    #@test all(expected_matches .== matches)

    k2 = map(CartesianIndex{2}, [(5, 5), (6, 6)])
    d2 = [[false, true, false],
          [false, false, false]]
    #matches = match_keypoints(k1, k2, d1, d2)
    expected_matches = [[CartesianIndex(3, 3), CartesianIndex(5, 5)],
                        [CartesianIndex(1, 1), CartesianIndex(6, 6)]]
    #@test all(expected_matches .== matches)
end

@testset "Grade Matches" begin
    k1 = map(CartesianIndex{2}, [(1,1),(2,1)])
    k2 = map(CartesianIndex{2}, [(1,1),(2,1)])
    @test grade_matches(k1, k2, 1) == 1
    @test grade_matches(k1, k2, 1, (i,j)->(abs(i[1]-j[1]) + abs(i[2]-j[2]))) == 1
    k2 = map(CartesianIndex{2}, [(1,0),(0,1)])
    @test grade_matches(k1, k2, 1) == 0
    @test grade_matches(k1, k2, 2) == 0.5
    k2 = map(CartesianIndex{2}, [(0,0),(2,1)])
    @test grade_matches(k1, k2, 1.2) == 0.5
end


text:
"""
Staggered Bernoulli scheme is a generalization of Bernoulli scheme were each positive outcome is possible
for different number of steps.

The constructor takes probabilities of succeses p_1, ... p_N-1 and
n_1, ... n_N-1, number of steps outcome i is possible.
"""
struct StaggeredBernoulliScheme{N}<:DiscreteCountingProcess
    n::SVector{N, Int64}
    p_suc::SVector{N, Float64}

    function StaggeredBernoulliScheme(n::AbstractArray{<:Number}, p_suc::AbstractArray{<:Number})
        @assert 0<sum(p_suc)<=1
        # @assert issorted(n)
        new{length(p_suc)}(SA[n...], SA[p_suc...])
    end
end

# # Distributions.jl interface (https://juliastats.org/Distributions.jl/stable/extends/)
# TODO: implment also for staggered case
# Base.length(bs::BernoulliScheme) = 2

# """
# Get step-shift of the next success and the next succes type.
# Support for step-shift is 1,2, ....
# Support for success type is 1, 2, ..., N  (N is the number of success results)
# """
# function Distributions._rand!(rng::AbstractRNG, bs::BernoulliScheme, x::AbstractVector{T}) where T <: Real
#     x[2] = Distributions.rand(rng, Categorical(bs.p_suc./sum(bs.p_suc)))
#     x[1] = Distributions.rand(rng, Geometric(sum(bs.p_suc)))+1
#     return x
# end

"""
Get parameters of the process
"""
Distributions.params(bs::StaggeredBernoulliScheme) = (bs.n, bs.p_suc)

"""
Get outcome distribution for given index.
Support for the outcome distribution is 1...N+1. N+1 means failure.
Index can be in the range 1, ..., n (n is the number of random variables in the scheme).
"""
function Base.getindex(bs::StaggeredBernoulliScheme, i::Integer)
    if 1<=i<=maximum(bs.n)
        active = bs.n .>= i
        active_p = bs.p_suc .* active
        return Categorical([active_p..., 1-sum(active_p)])
    else
        throw(BoundsError(bs, i))
    end
end


text:
@testset "Algebraic Connectivity: Optimal solution tests" begin

    file_path = joinpath(@__DIR__,"..", "examples/solver.jl")
    include(file_path)

    params = Dict{String, Any}(
        "num_nodes" => 5,
        "instance" => 1,
        "data_type" => "old",

        "solution_type" => "exact",
        "presolve" => true,
        "optimizer_log" => true,
        "relax_integrality" => false
    )

    result_lo = LaplacianOpt.run_LOpt_model(params, glpk_optimizer)

    @test result_lo["termination_status"] == MOI.OPTIMAL
    @test result_lo["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_lo["objective"], 11.1592, atol=1E-4)
    @test isapprox(result_lo["solution"]["z_var"][1,4], 1.0)
    @test isapprox(result_lo["solution"]["z_var"][4,1], 1.0)
    @test isapprox(result_lo["solution"]["z_var"][2,4], 1.0)
    @test isapprox(result_lo["solution"]["z_var"][4,2], 1.0)
    @test isapprox(result_lo["solution"]["z_var"][3,4], 1.0)
    @test isapprox(result_lo["solution"]["z_var"][4,3], 1.0)
    @test isapprox(result_lo["solution"]["z_var"][3,5], 1.0)
    @test isapprox(result_lo["solution"]["z_var"][5,3], 1.0)

    @test LA.issymmetric(result_lo["solution"]["W_var"])

    L_val = result_lo["solution"]["W_var"]

    for i=1:params["num_nodes"]
        L_val[i,i] += result_lo["solution"]["γ_var"] * (params["num_nodes"] - 1)/params["num_nodes"]
    end

    for i=1:(params["num_nodes"] - 1)
        for j=(i+1):params["num_nodes"]
            L_val[i,j] -= result_lo["solution"]["γ_var"] / params["num_nodes"]
            L_val[j,i] = L_val[i,j]
        end
    end

    @test isapprox(LA.eigvals(L_val)[2], result_lo["objective"])
    @test isapprox(LA.eigvals(L_val)[2], result_lo["solution"]["γ_var"])
end

@testset "Max Span Tree: Optimal solution tests" begin

    file_path = joinpath(@__DIR__,"..", "examples/solver.jl")
    include(file_path)

    params = Dict{String, Any}(
        "num_nodes" => 5,
        "instance" => 1
        )

    result_mst = LaplacianOpt.run_MaxSpanTree_model(params, glpk_optimizer)

    @test result_mst["termination_status"] == MOI.OPTIMAL
    @test result_mst["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_mst["objective"], 113.7588, atol=1E-4)
    @test isapprox(result_mst["solution"]["z_var"][1,2], 1.0)
    @test isapprox(result_mst["solution"]["z_var"][2,1], 1.0)
    @test isapprox(result_mst["solution"]["z_var"][2,5], 1.0)
    @test isapprox(result_mst["solution"]["z_var"][5,2], 1.0)
    @test isapprox(result_mst["solution"]["z_var"][3,4], 1.0)
    @test isapprox(result_mst["solution"]["z_var"][4,3], 1.0)
    @test isapprox(result_mst["solution"]["z_var"][3,5], 1.0)
    @test isapprox(result_mst["solution"]["z_var"][5,3], 1.0)

end

@testset "Max Span Tree: Lazy callback tests" begin

    file_path = joinpath(@__DIR__,"..", "examples/solver.jl")
    include(file_path)

    params = Dict{String, Any}(
        "num_nodes" => 15,
        "instance" => 5,
        "data_type" => "old",
        "eigen_cuts_full" => false,
        "topology_flow_cuts" => true
        )

    result_mst = LaplacianOpt.run_MaxSpanTree_model(params, glpk_optimizer, lazy_callback=true)

    @test result_mst["termination_status"] == MOI.OPTIMAL
    @test result_mst["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_mst["objective"], 3576.56202331, atol=1E-6)
    @test isapprox(result_mst["solution"]["z_var"][1,4], 1.0)
    @test isapprox(result_mst["solution"]["z_var"][2,4], 1.0)
    @test isapprox(result_mst["solution"]["z_var"][3,4], 1.0)
    @test isapprox(result_mst["solution"]["z_var"][4,15], 1.0)
    @test isapprox(result_mst["solution"]["z_var"][7,15], 1.0)
    @test isapprox(result_mst["solution"]["z_var"][8,15], 1.0)
end

@testset "SOC relaxations - 1: constraint_soc_linearized_cuts_on_2minors" begin

    file_path = joinpath(@__DIR__,"..", "examples/solver.jl")
    include(file_path)

    params_1 = Dict{String, Any}(
        "num_nodes" => 5,
        "instance" => 5,
        "eigen_cuts_full" => false,
        "soc_linearized_cuts" => true,
        "topology_flow_cuts" => true
    )

    result_1 = LaplacianOpt.run_LOpt_model(params_1, glpk_optimizer)

    @test result_1["termination_status"] == MOI.OPTIMAL
    @test result_1["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_1["objective"], 15.992792245, atol=1E-6)
    @test isapprox(result_1["solution"]["z_var"][1,5], 1.0)
    @test isapprox(result_1["solution"]["z_var"][2,5], 1.0)
    @test isapprox(result_1["solution"]["z_var"][3,4], 1.0)
    @test isapprox(result_1["solution"]["z_var"][3,5], 1.0)
    @test isapprox(result_1["solution"]["z_var"][4,3], 1.0)

    params_2 = Dict{String, Any}(
        "num_nodes" => 5,
        "instance" => 5,
        "eigen_cuts_full" => true,
        "soc_linearized_cuts" => true
    )

    result_2 = LaplacianOpt.run_LOpt_model(params_2, glpk_optimizer)

    @test result_2["termination_status"] == MOI.OPTIMAL
    @test result_2["primal_status"] == MOI.FEASIBLE_POINT
    @test result_2["objective"] <= result_1["objective"]
    @test isapprox(result_2["objective"], 13.065372790, atol=1E-6)
    @test isapprox(result_1["solution"]["z_var"], result_2["solution"]["z_var"], atol = 1E-6)
end


text:
module ReferenceStats

using SparseArrays
using Statistics
using Interpolations
using LinearAlgebra
using DocStringExtensions

using ..ReferenceModels
using ..ModelTypes
using ..LESUtils
using ..HelperFuncs
using ..DistributionUtils

export ReferenceStatistics
export pca_length, full_length
export get_obs, get_profile, obs_PCA, pca
export pca_inds, full_inds

"""
    ReferenceStatistics{FT <: Real, IT <: Integer}

A structure containing statistics from the reference model used to
define a well-posed inverse problem.

# Fields

$(TYPEDFIELDS)

# Constructors

    ReferenceStatistics(y::Vector{FT}, args...) where {FT <: Real}
    ReferenceStatistics(
        RM::Vector{ReferenceModel};
        perform_PCA::Bool = true,
        normalize::Bool = true,
        variance_loss::FT = 0.1,
        tikhonov_noise::FT = 0.0,
        tikhonov_mode::String = "absolute",
        dim_scaling::Bool = false,
        y_type::ModelType = LES(),
        Σ_type::ModelType = LES(),
        Δt::FT = 6 * 3600.0,
    ) where {FT <: Real}

Constructs the ReferenceStatistics defining the inverse problem.

Inputs:
 - `RM`               :: Vector of `ReferenceModel`s.
 - `perform_PCA`      :: Boolean specifying whether to perform PCA.
 - `normalize`        :: Boolean specifying whether to normalize the data.
 - `variance_loss`    :: Fraction of variance loss when performing PCA.
 - `tikhonov_noise`   :: Tikhonov regularization factor for covariance matrices.
 - `tikhonov_mode`    :: If "relative", tikhonov_noise is scaled by the maximum
    eigenvalue in the covariance matrix considered, having the interpretation of
    the inverse of the desired condition number. This value is enforced to be
    larger than the sqrt of the machine precision for stability.
 - `dim_scaling`      :: Whether to scale covariance blocks by their size.
 - `y_type`           :: Type of reference mean data. Either LES() or SCM().
 - `Σ_type`           :: Type of reference covariance data. Either LES() or SCM().
 - `Δt`               :: [LES last time - SCM start time (LES timeframe)] for `LES_driven_SCM` cases.
"""
Base.@kwdef struct ReferenceStatistics{FT <: Real, IT <: Integer}
    "Reference data, length: nSim * n_vars * n_zLevels (possibly reduced by PCA)"
    y::Vector{FT}
    "Data covariance matrix, dims: (y,y) (possibly reduced by PCA)"
    Γ::Matrix{FT}
    "Vector (length: `nSim`) of normalizing factors (length: `n_vars`)"
    norm_vec::Vector{Vector{FT}}
    "Vector (length: `nSim`) of PCA projection matrices with leading eigenvectors as columns"
    pca_vec::Vector{Union{Matrix{FT}, UniformScaling}}
    "Full reference data vector, length: `nSim * n_vars * n_zLevels`"
    y_full::Vector{FT}
    "Full covariance matrix, dims: (y,y)"
    Γ_full::SparseMatrixCSC{FT, Int64}
    "Degrees of freedom per case (possibly reduced by PCA)"
    ndof_case::Vector{IT}
    "Full degrees of freedom per case"
    ndof_full_case::Vector{IT}

    # Constructors

    ReferenceStatistics(y::Vector{FT}, args...) where {FT <: Real} = new{FT, Int64}(y, args...)

    function ReferenceStatistics(
        RM::Vector{ReferenceModel};
        perform_PCA::Bool = true,
        normalize::Bool = true,
        variance_loss::FT = 0.1,
        tikhonov_noise::FT = 0.0,
        tikhonov_mode::String = "absolute",
        dim_scaling::Bool = false,
        y_type::ModelType = LES(),
        Σ_type::ModelType = LES(),
        Δt::FT = 6 * 3600.0,
    ) where {FT <: Real}
        IT = Int64
        # Init arrays
        y = FT[]
        Γ_vec = Matrix{FT}[]
        y_full = FT[]
        Γ_full_vec = Matrix{FT}[]
        pca_vec = []
        norm_vec = Vector[]
        ndof_case = IT[]
        ndof_full_case = IT[]

        for m in RM
            model = m.case_name == "LES_driven_SCM" ? time_shift_reference_model(m, Δt) : m
            # Get (interpolated and pool-normalized) observations, get pool variance vector
            y_, y_var_, pool_var = get_obs(model, y_type, Σ_type, normalize, z_scm = get_z_obs(model))
            push!(norm_vec, pool_var)
            if perform_PCA
                y_pca, y_var_pca, P_pca = obs_PCA(y_, y_var_, variance_loss)
                append!(y, y_pca)
                push!(Γ_vec, y_var_pca)
                push!(pca_vec, P_pca)
                push!(ndof_case, length(y_pca))
            else
                append!(y, y_)
                push!(Γ_vec, y_var_)
                push!(pca_vec, 1.0I)
                push!(ndof_case, length(y_))
            end
            # Save full dimensionality (normalized) output for error computation
            append!(y_full, y_)
            push!(Γ_full_vec, y_var_)
            push!(ndof_full_case, length(y_))
        end

        # Construct global observational covariance matrix, original space
        Γ_full = sparse(cat(Γ_full_vec..., dims = (1, 2)))

        # Scale by number of dimensions (averaging loss per dimension)
        Γ_vec = dim_scaling ? length(Γ_vec) .* map(x -> size(x, 1) * x, Γ_vec) : Γ_vec
        # Construct global observational covariance matrix, PCA
        Γ = cat(Γ_vec..., dims = (1, 2))
        # Condition global covariance matrix, PCA
        if tikhonov_mode == "relative"
            @assert perform_PCA "Relative Tikhonov mode only available after PCA change of basis."
            tikhonov_noise = max(tikhonov_noise, 10 * sqrt(eps(FT)))
            Γ = Γ + tikhonov_noise * maximum(diag(Γ)) * I
        else
            Γ = Γ + tikhonov_noise * I
        end

        @assert isposdef(Γ) "Covariance matrix Γ is ill-conditioned, consider regularization."

        return new{FT, IT}(y, Γ, norm_vec, pca_vec, y_full, Γ_full, ndof_case, ndof_full_case)
    end

end

"Returns dimensionality of the ReferenceStatistics in low-dimensional latent space"
pca_length(RS::ReferenceStatistics) = length(RS.y)

"Returns full dimensionality of the ReferenceStatistics, before latent space encoding"
full_length(RS::ReferenceStatistics) = length(RS.y_full)
pca_length(RS::ReferenceStatistics, case_ind) = RS.ndof_case[case_ind]
full_length(RS::ReferenceStatistics, case_ind) = RS.ndof_full_case[case_ind]

"""Fetch the dof indices (possibly reduced by PCA) for the a case, specified by its in config["reference"]["case_name"]."""
pca_inds(RS::ReferenceStatistics, case_ind) = dof_inds(RS.ndof_case, case_ind)
"""Fetch the dof indices for the a case, specified by its in config["reference"]["case_name"]."""
full_inds(RS::ReferenceStatistics, case_ind) = dof_inds(RS.ndof_full_case, case_ind)

"""Given a vector of `ndofs` for some ordered collection, fetch the dof indices for one of the elements by its index in the ordered collection."""
function dof_inds(ndofs::Vector{IT}, ind::IT) where {IT <: Integer}
    !(1 ≤ ind ≤ length(ndofs)) && throw(ArgumentError("Index `$ind` must be between 1 and $(length(ndofs))."))
    start_ind = (1, 1 .+ cumsum(ndofs)...)[ind]
    stop_ind = cumsum(ndofs)[ind]
    start_ind:stop_ind
end

"""
    get_obs(
        m::ReferenceModel,
        y_names::Vector{String},
        Σ_names::Vector{String},
        normalize::Bool;
        z_scm::Union{Vector{FT}, Nothing} = nothing,
    )
    get_obs(
        m::ReferenceModel,
        y_type::Union{LES, SCM},
        Σ_type::Union{LES, SCM},
        normalize::Bool;
        z_scm::Union{Vector{FT}, Nothing},
    )

Get observations for variables `y_names`, interpolated to
`z_scm` (if given), and possibly normalized with respect to the pooled variance.

Inputs:
 - `m`            :: Reference model
 - `y_names`      :: Names of observed fields from the ReferenceModel `m`.
 - `Σ_names`      :: Names of fields used to construct covariances, may be different than `y_names`
    if there are LES/SCM name discrepancies.
 - `normalize`    :: Whether to normalize the observations.
 - `z_scm`        :: If given, interpolate LES observations to given levels.
 - `y_type`       :: ModelType used to construct observations, ::LES or ::SCM.
 - `y_type`       :: ModelType used to construct noise covariances, ::LES or ::SCM.
Outputs:
 - `y`            :: Mean of observations, possibly interpolated to z_scm levels.
 - `Σ`            :: Observational covariance matrix, possibly pool-normalized.
 - `pool_var`     :: Vector of vertically averaged time-variance, one entry for each variable
"""
function get_obs(
    m::ReferenceModel,
    y_names::Vector{String},
    Σ_names::Vector{String},
    normalize::Bool;
    z_scm::Union{Vector{FT}, Nothing} = nothing,
) where {FT <: Real}
    # time covariance
    Σ, pool_var = get_time_covariance(m, Σ_names, z_scm = z_scm)
    # normalization
    norm_vec = normalize ? pool_var : ones(size(pool_var))
    # Get true observables
    y = get_profile(m, y_nc_file(m), y_names, z_scm = z_scm)
    # normalize
    y = normalize_profile(y, num_vars(m), norm_vec)
    return y, Σ, norm_vec
end

function get_obs(
    m::ReferenceModel,
    y_type::Union{LES, SCM},
    Σ_type::Union{LES, SCM},
    normalize::Bool;
    z_scm::Union{Vector{FT}, Nothing},
) where {FT <: Real}
    y_names = isa(y_type, LES) ? get_les_names(m, y_nc_file(m)) : m.y_names
    Σ_names = isa(Σ_type, LES) ? get_les_names(m, Σ_nc_file(m)) : m.y_names
    get_obs(m, y_names, Σ_names, normalize, z_scm = z_scm)
end

"""
    obs_PCA(y_mean, y_var, allowed_var_loss = 1.0e-1)

Perform dimensionality reduction using principal component analysis on
the variance `y_var`. Only eigenvectors with eigenvalues that contribute
to the leading 1-`allowed_var_loss` variance are retained.
Inputs:

 - `y_mean`           :: Mean of the observations.
 - `y_var`            :: Variance of the observations.
 - `allowed_var_loss` :: Maximum variance loss allowed.

Outputs:

 - `y_pca`            :: Projection of `y_mean` onto principal subspace spanned by eigenvectors.
 - `y_var_pca`        :: Projection of `y_var` on principal subspace.
 - `P_pca`            :: Projection matrix onto principal subspace, with leading eigenvectors as columns.
"""
function obs_PCA(y_mean::Vector{FT}, y_var::Matrix{FT}, allowed_var_loss::FT = 1.0e-1) where {FT <: Real}
    λ_pca, P_pca = pca(y_var, allowed_var_loss)
    # Project mean
    y_pca = P_pca' * y_mean
    y_var_pca = Diagonal(λ_pca)
    return y_pca, y_var_pca, P_pca
end

"""
    pca(covmat::AbstractMatrix{FT}, allowed_var_loss::FT) where {FT <: Real}

Perform dimensionality reduction using principal component analysis on
the variance covmat.

Inputs:
 - `covmat`           :: Variance of the observations.
 - `allowed_var_loss` :: Maximum variance loss allowed.
Outputs:
 - `λ_pca`            :: Principal eigenvalues, ordered in increasing value order.
 - `P_pca`            :: Projection matrix onto principal subspace, with leading eigenvectors as columns.
"""
function pca(covmat::AbstractMatrix{FT}, allowed_var_loss::FT) where {FT <: Real}
    eigvals, eigvecs = eigen(covmat)
    # Get index of leading eigenvalues, eigvals are ordered from low to high in julia
    # This expression recovers 1 extra eigenvalue compared to threshold
    leading_eigs = findall(<(1.0 - allowed_var_loss), -cumsum(eigvals) / sum(eigvals) .+ 1)
    P_pca = eigvecs[:, leading_eigs]
    λ_pca = eigvals[leading_eigs]
    return λ_pca, P_pca
end

"""
    get_profile(
        filename::String,
        y_names::Vector{String};
        ti::Real = 0.0,
        tf = nothing,
        z_scm::Union{Vector{FT}, Nothing} = nothing,
    )
    get_profile(m::ReferenceModel, filename::String; z_scm::Union{Vector{T}, T} = nothing) where {T}
    get_profile(
        m::ReferenceModel,
        filename::String,
        y_names::Vector{String};
        z_scm::Union{Vector{T}, T} = nothing,
    ) where {T}

Get time-averaged profiles for variables `y_names`, interpolated to
`z_scm` (if given), and concatenated into a single output vector.

Inputs:

 - `filename`    :: nc filename
 - `y_names`   :: Names of variables to be retrieved.
 - `ti`          :: Initial time of averaging window.
 - `tf`          :: Final time of averaging window.
 - `z_scm`       :: If given, interpolate LES observations to given levels.
 - `m`           :: ReferenceModel from which to fetch profiles, implicitly defines `ti` and `tf`.

Outputs:

 - `y` :: Output vector used in the inverse problem, which concatenates the requested profiles.
"""
function get_profile(
    filename::String,
    y_names::Vector{String};
    ti::Real = 0.0,
    tf::Union{Real, Nothing} = nothing,
    z_scm::Union{Vector{T}, T} = nothing,
) where {T}

    t = nc_fetch(filename, "t")
    dt = length(t) > 1 ? mean(diff(t)) : 0.0
    y = zeros(0)

    # Check that times are contained in simulation output
    Δt_start, ti_index = findmin(broadcast(abs, t .- ti))
    # If simulation does not contain values for ti or tf, return high value (penalization)
    if t[end] < ti
        @warn string(
            "Note: t_end < ti, which means that simulation stopped before reaching the requested t_start.",
            "Requested t_start = $ti s. However, the last time available is $(t[end]) s.",
            "Defaulting to penalized profiles...",
        )
        for i in 1:length(y_names)
            var_ = isnothing(z_scm) ? get_height(filename) : z_scm
            append!(y, 1.0e5 * ones(length(var_[:])))
        end
        return y
    end
    if !isnothing(tf)
        Δt_end, tf_index = findmin(broadcast(abs, t .- tf))
        if t[end] < tf - dt
            @warn string(
                "Note: t_end < tf - dt, which means that simulation stopped before reaching the requested t_end.",
                "Requested t_end = $tf s. However, the last time available is $(t[end]) s.",
                "Defaulting to penalized profiles...",
            )
            for i in 1:length(y_names)
                var_ = isnothing(z_scm) ? get_height(filename) : z_scm
                append!(y, 1.0e5 * ones(length(var_[:])))
            end
            return y
        end
    end

    # Return time average for non-degenerate cases
    for var_name in y_names
        var_ = fetch_interpolate_transform(var_name, filename, z_scm)
        var_mean = !isnothing(tf) ? mean(var_[:, ti_index:tf_index], dims = 2) : var_[:, ti_index]
        append!(y, var_mean)
    end
    return y
end

function get_profile(m::ReferenceModel, filename::String; z_scm::Union{Vector{T}, T} = nothing) where {T}
    get_profile(m, filename, m.y_names, z_scm = z_scm)
end

function get_profile(
    m::ReferenceModel,
    filename::String,
    y_names::Vector{String};
    z_scm::Union{Vector{T}, T} = nothing,
) where {T}
    get_profile(filename, y_names, ti = get_t_start(m), tf = get_t_end(m), z_scm = z_scm)
end

"""
    get_time_covariance(m::ReferenceModel, y_names::Vector{String}; z_scm::Vector{FT}) where {FT <: Real}

Obtain the covariance matrix of a group of profiles, where the covariance
is obtained in time.

Inputs:
 - `m`            :: Reference model.
 - `y_names`    :: List of variable names to be included.
 - `z_scm`        :: If given, interpolates covariance matrix to this locations.
"""
function get_time_covariance(m::ReferenceModel, y_names::Vector{String}; z_scm::Vector{FT}) where {FT <: Real}
    filename = Σ_nc_file(m)
    t = nc_fetch(filename, "t")
    # Find closest interval in data
    ti_index = argmin(broadcast(abs, t .- get_t_start_Σ(m)))
    tf_index = argmin(broadcast(abs, t .- get_t_end_Σ(m)))
    N_samples = length(ti_index:tf_index)
    ts_vec = zeros(0, N_samples)
    num_outputs = length(y_names)
    pool_var = zeros(num_outputs)

    for (i, var_name) in enumerate(y_names)
        var_ = fetch_interpolate_transform(var_name, filename, z_scm)
        # Store pooled variance
        pool_var[i] = mean(var(var_[:, ti_index:tf_index], dims = 2)) + eps(FT) # vertically averaged time-variance of variable
        # Normalize timeseries
        ts_var_i = var_[:, ti_index:tf_index] ./ sqrt(pool_var[i])
        ts_vec = cat(ts_vec, ts_var_i, dims = 1)  # dims: (Nz*num_outputs, Nt)
    end
    cov_mat = cov(ts_vec, dims = 2)  # covariance, w/ samples across time dimension (t_inds).
    return cov_mat, pool_var
end

end # module


text:
using Heuristics, StaticArrays

# Particle State
N = 5
M = 5
T = typeof(10.0)
s = Heuristics.Swarm{T,N,M}(undef)

# Initialization
@test s.particles isa SizedArray{Tuple{M},Heuristics.Particle{T,N},1,1,Vector{Heuristics.Particle{T,N}}}

# Initializing temp variables
tempXs  = rand(N,M)
tempPs  = rand(N,M)
tempVs  = rand(N,M)
tempD   = rand(N)
tempFxs = rand(M)
tempFps = rand(M)
tempN   = 10.0
tempB   = rand()
tempW   = rand()
tempC   = 11.0
tempY₁  = rand()
tempY₂  = rand()

# Setting values
s.b  = tempB
s.n  = tempN
s.w  = tempW
s.c  = tempC
s.y₁ = tempY₁
s.y₂ = tempY₂
s.d .= tempD
@inbounds for i in 1:M
    s[i].fx = tempFxs[i]
    s[i].fp = tempFps[i]
    for j in 1:N
        s[i].x[j] = tempXs[j]
        s[i].p[j] = tempPs[j]
        s[i].v[j] = tempVs[j]
    end
end

# Testing set values
@test s.particles isa SizedArray{Tuple{M},Heuristics.Particle{T,N},1,1,Vector{Heuristics.Particle{T,N}}}
@test s.b isa T
@test s.d isa SizedArray{Tuple{N},T,1,1,Vector{T}}
@test s.n isa UInt16
@test s.w isa T
@test s.c isa UInt16
@test s.y₁ isa T
@test s.y₂ isa T
@test s.b  == T(tempB)
@test s.n  == UInt16(tempN)
@test s.w  == T(tempW)
@test s.c  == UInt16(tempC)
@test s.y₁ == T(tempY₁)
@test s.y₂ == T(tempY₂)
@inbounds for i in 1:N
    @test s.d[i] == T(tempD[i])
end
@inbounds for i in 1:M
    @test s[i].fx == T(tempFxs[i])
    @test s[i].fp == T(tempFps[i])
    for j in 1:N
        @test s[i].x[j] == T(tempXs[j])
        @test s[i].p[j] == T(tempPs[j])
        @test s[i].v[j] == T(tempVs[j])
    end
end

# Testing length(S::Swarm)
@test length(s) == M

# Testing setindex!
pNew = Heuristics.Particle{T,N}(undef)
s[1] = pNew
@test s[1] == pNew


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/LibUnwind_jll.jl

baremodule LibUnwind_jll
using Base, Libdl
Base.Experimental.@compiler_options compile=min optimize=0 infer=false

const PATH_list = String[]
const LIBPATH_list = String[]

export libunwind

# These get calculated in __init__()
const PATH = Ref("")
const LIBPATH = Ref("")
artifact_dir = ""
libunwind_handle = C_NULL
libunwind_path = ""

const libunwind = "libunwind.so.8"

function __init__()
    # We only do something on Linux/FreeBSD
    @static if Sys.islinux() || Sys.isfreebsd()
        global libunwind_handle = dlopen(libunwind)
        global libunwind_path = dlpath(libunwind_handle)
        global artifact_dir = dirname(Sys.BINDIR)
        LIBPATH[] = dirname(libunwind_path)
        push!(LIBPATH_list, LIBPATH[])
    end
end

# JLLWrappers API compatibility shims.  Note that not all of these will really make sense.
# For instance, `find_artifact_dir()` won't actually be the artifact directory, because
# there isn't one.  It instead returns the overall Julia prefix.
is_available() = @static (Sys.islinux() || Sys.isfreebsd()) ? true : false
find_artifact_dir() = artifact_dir
dev_jll() = error("stdlib JLLs cannot be dev'ed")
best_wrapper = nothing
get_libunwind_path() = libunwind_path

end  # module LibUnwind_jll

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
function gpu_track(
    signal,
    state::GNSSBenchmarks.gpuTrackingState{S, C, CALF, COLF, CN, DS},
    prn::Integer,
    sampling_frequency;
    post_corr_filter = get_default_post_corr_filter(gpu_get_correlator(state)),
    intermediate_frequency = 0.0Hz,
    max_integration_time::typeof(1ms) = 1ms,
    min_integration_time::typeof(1.0ms) = 0.75ms,
    early_late_sample_shift = get_early_late_sample_shift(S,
        gpu_get_correlator(state), sampling_frequency, 0.5),
    carrier_loop_filter_bandwidth = 18Hz,
    code_loop_filter_bandwidth = 1Hz,
    velocity_aiding = 0Hz,
    carrier_amplitude_power::Val{N} = Val(5)
) where {
S <: AbstractGNSSSystem,
C <: AbstractCorrelator,
CALF <: AbstractLoopFilter,
COLF <: AbstractLoopFilter,
CN <: AbstractCN0Estimator,
DS <: CuArray,
N
}
correlator = gpu_get_correlator(state)
downconverted_signal = resize!(gpu_get_downconverted_signal(state), size(signal, 1))
carrier_replica = resize!(gpu_get_carrier(state), size(signal, 1))
code_replica = resize!(gpu_get_code(state), size(signal, 1) + 2 * maximum(early_late_sample_shift))
init_carrier_doppler = GNSSBenchmarks.gpu_get_init_carrier_doppler(state)
init_code_doppler = GNSSBenchmarks.gpu_get_init_code_doppler(state)
carrier_doppler = GNSSBenchmarks.gpu_get_carrier_doppler(state)
code_doppler = GNSSBenchmarks.gpu_get_code_doppler(state)
carrier_phase = GNSSBenchmarks.gpu_get_carrier_phase(state)
code_phase = GNSSBenchmarks.gpu_get_code_phase(state)
sc_bit_detector = GNSSBenchmarks.gpu_get_sc_bit_detector(state)
carrier_loop_filter = GNSSBenchmarks.gpu_get_carrier_loop_filter(state)
code_loop_filter = GNSSBenchmarks.gpu_get_code_loop_filter(state)
prompt_accumulator = GNSSBenchmarks.gpu_get_prompt_accumulator(state)
integrated_samples = GNSSBenchmarks.gpu_get_integrated_samples(state)
cn0_estimator = GNSSBenchmarks.gpu_get_cn0_estimator(state)
signal_start_sample = 1
bit_buffer = Tracking.BitBuffer()
valid_correlator = zero(correlator)
got_correlator = false
while true
    num_samples_left_to_integrate = get_num_samples_left_to_integrate(
        S,
        max_integration_time,
        sampling_frequency,
        code_doppler,
        code_phase,
        found(sc_bit_detector)
    )
    signal_samples_left = get_num_samples(signal) - signal_start_sample + 1
    num_samples_left = min(num_samples_left_to_integrate, signal_samples_left)
    carrier_frequency = get_current_carrier_frequency(
        intermediate_frequency,
        carrier_doppler
    )
    code_frequency = get_current_code_frequency(S, code_doppler)
    code_replica = GNSSBenchmarks.gpu_gen_code_replica!(
        code_replica,
        Tracking.GPSL1,
        code_frequency,
        sampling_frequency,
        code_phase,
        signal_start_sample,
        num_samples_left,
        early_late_sample_shift,
        prn
    )
    carrier_replica = GNSSBenchmarks.gpu_gen_carrier_replica!(
        carrier_replica,
        carrier_frequency,
        sampling_frequency,
        carrier_phase,
        carrier_amplitude_power,
        signal_start_sample,
        num_samples_left
    )
    downconverted_signal = GNSSBenchmarks.gpu_downconvert!(
        downconverted_signal,
        signal,
        carrier_replica,
        signal_start_sample,
        num_samples_left
    )
    correlator = GNSSBenchmarks.gpu_correlate(
        correlator,
        downconverted_signal,
        code_replica,
        early_late_sample_shift,
        signal_start_sample,
        num_samples_left,
        1,
        2,
        Val(7)
    )
    integrated_samples += num_samples_left
    carrier_phase = Tracking.update_carrier_phase(
        num_samples_left,
        carrier_frequency,
        sampling_frequency,
        carrier_phase,
        carrier_amplitude_power
    )
    prev_code_phase = code_phase
    code_phase = Tracking.update_code_phase(
        S,
        num_samples_left,
        code_frequency,
        sampling_frequency,
        code_phase,
        found(sc_bit_detector)
    )
    integration_time = integrated_samples / sampling_frequency
    if num_samples_left == num_samples_left_to_integrate &&
            integration_time >= min_integration_time
        got_correlator = true

        correlator = Tracking.normalize(correlator, integrated_samples)
        valid_correlator = correlator
        filtered_correlator = Tracking.filter(post_corr_filter, correlator)
        pll_discriminator = Tracking.pll_disc(S, filtered_correlator)
        dll_discriminator = Tracking.dll_disc(
            S,
            filtered_correlator,
            early_late_sample_shift,
            code_frequency / sampling_frequency
        )
        carrier_freq_update, carrier_loop_filter = filter_loop(
            carrier_loop_filter,
            pll_discriminator,
            integration_time,
            carrier_loop_filter_bandwidth
        )
        code_freq_update, code_loop_filter = filter_loop(
            code_loop_filter,
            dll_discriminator,
            integration_time,
            code_loop_filter_bandwidth
        )
        carrier_doppler, code_doppler = Tracking.aid_dopplers(
            S,
            init_carrier_doppler,
            init_code_doppler,
            carrier_freq_update,
            code_freq_update,
            velocity_aiding
        )
        cn0_estimator = Tracking.update(cn0_estimator, Tracking.get_prompt(filtered_correlator))
        bit_buffer, prompt_accumulator = Tracking.buffer(
            S,
            bit_buffer,
            prompt_accumulator,
            found(sc_bit_detector),
            prev_code_phase,
            code_phase,
            max_integration_time,
            Tracking.get_prompt(filtered_correlator)
        )
        sc_bit_detector = Tracking.find(S, sc_bit_detector, Tracking.get_prompt(filtered_correlator))
        correlator = zero(correlator)
        integrated_samples = 0
    end

    num_samples_left == signal_samples_left && break
    signal_start_sample += num_samples_left
end
next_state = GNSSBenchmarks.gpuTrackingState{S, C, CALF, COLF, CN, DS}(
    init_carrier_doppler,
    init_code_doppler,
    carrier_doppler,
    code_doppler,
    carrier_phase,
    code_phase,
    correlator,
    carrier_loop_filter,
    code_loop_filter,
    sc_bit_detector,
    integrated_samples,
    prompt_accumulator,
    cn0_estimator,
    downconverted_signal,
    carrier_replica,
    code_replica
)
estimated_cn0 = Tracking.estimate_cn0(cn0_estimator, max_integration_time)
gpuTrackingResults(next_state, valid_correlator, got_correlator, bit_buffer, estimated_cn0)
end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

# Compute X <- a X + b Y.
function matfun_axpby!(X,a,b,Y)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:m
        @inbounds for j=1:n
            if (b isa ValueOne)
                X[i,j]+=Y[i,j]
            else
                X[i,j]+=b*Y[i,j]
            end
        end
    end
end

@inline function exp_bbcs_m2(A)
    T=promote_type(eltype(A),ComplexF64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_bbcs_m2!(A_copy)
end

@inline function exp_bbcs_m2!(A)
    T=promote_type(eltype(A),ComplexF64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=4
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: B2 Bb3 B3 T2k5
    # Computing B2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Computing Bb3 = x*A+x*B2
    coeff1=0.0 + 0.16666657785001893im
    coeff2=0.04166664890333649 + 0.0im
    memslots3 .= coeff1.*memslots1 .+ coeff2.*memslots2
    # Computing B3 with operation: mult
    mul!(memslots4,memslots2,memslots3)
    # Deallocating Bb3 in slot 3
    # Computing T2k5 = x*I+x*A+x*B2+x*B3
    coeff1=1.0 + 0.0im
    coeff2=0.0 - 0.9999999999998107im
    coeff3=-0.4999999999999432 + 0.0im
    coeff4=1.0 + 0.0im
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots4
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 2
    # Deallocating B3 in slot 4
    return memslots1 # Returning T2k5
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
function relevant(interpreter::Interpreter, domain::Domain, state::State)
    actions = Term[]
    for act in values(get_actions(domain))
        act_name = get_name(act)
        act_vars, act_types = get_argvars(act), get_argtypes(act)
        # Compute postconditions from the action's effect
        diff = effect_diff(domain, state, get_effect(act))
        addcond = [Compound(:or, diff.add)]
        delcond = [@julog(not(:t)) for t in diff.del]
        typecond = [@julog($ty(:v)) for (v, ty) in zip(act_vars, act_types)]
        # Include type conditions when necessary for correctness
        if any(has_func(c, domain) ||
               has_quantifier(c) for c in [addcond; delcond])
            conds = [typecond; addcond; delcond]
        else
            conds = [addcond; typecond; delcond]
        end
        # Find all substitutions that satisfy the postconditions
        subst = satisfiers(interpreter, domain, state, conds)
        if isempty(subst) continue end
        for s in subst
            args = [get(s, var, var) for var in act_vars]
            if any(!is_ground(a) for a in args) continue end
            term = isempty(args) ? Const(act_name) : Compound(act_name, args)
            push!(actions, term)
        end
    end
    return actions
end

function relevant(interpreter::Interpreter,
                  domain::Domain, state::State, action::Action, args)
   if any(!is_ground(a) for a in args)
       error("Not all arguments are ground.")
   end
   act_vars, act_types = get_argvars(action), get_argtypes(action)
   subst = Subst(var => val for (var, val) in zip(act_vars, args))
   # Compute postconditions from the action's effect
   diff = effect_diff(domain, state, substitute(get_effect(action), subst))
   postcond = Term[Compound(:or, diff.add); [@julog(not(:t)) for t in diff.del]]
   # Construct type conditions of the form "type(val)"
   typecond = (all(ty == :object for ty in act_types) ? Term[] :
               [@julog($ty(:v)) for (v, ty) in zip(args, act_types)])
   # Check whether postconditions hold
   return satisfy(interpreter, domain, state, [postcond; typecond])
end


text:
# convenience and indirect construction

# conversions from CPU arrays rely on constructors
Base.convert(::Type{T}, a::AbstractArray) where {T<:AbstractGPUArray} = a isa T ? a : T(a)
# TODO: can we implement constructors to and from ::AbstractArray here? by calling the undef
#       constructor and doing a `copyto!`. this is tricky, due to ambiguities, and no easy
#       way to go from <:AbstractGPUArray{T,N} to e.g. CuArray{S,N}

## convenience constructors

function Base.fill!(A::AnyGPUArray{T}, x) where T
    length(A) == 0 && return A
    gpu_call(A, convert(T, x)) do ctx, a, val
        idx = @linearidx(a)
        @inbounds a[idx] = val
        return
    end
    A
end

## identity matrices

function identity_kernel(ctx::AbstractKernelContext, res::AbstractArray{T}, stride, val) where T
    i = linear_index(ctx)
    i > stride && return
    ilin = (stride * (i - 1)) + i
    @inbounds res[ilin] = val
    return
end

function (T::Type{<: AnyGPUArray{U}})(s::UniformScaling, dims::Dims{2}) where {U}
    res = similar(T, dims)
    fill!(res, zero(U))
    gpu_call(identity_kernel, res, size(res, 1), s.λ; elements=minimum(dims))
    res
end

(T::Type{<: AnyGPUArray})(s::UniformScaling{U}, dims::Dims{2}) where U = T{U}(s, dims)

(T::Type{<: AnyGPUArray})(s::UniformScaling, m::Integer, n::Integer) = T(s, Dims((m, n)))

function Base.copyto!(A::AbstractGPUMatrix{T}, s::UniformScaling) where T
    fill!(A, zero(T))
    gpu_call(identity_kernel, A, size(A, 1), s.λ; elements=minimum(size(A)))
    A
end

function _one(unit::T, x::AbstractGPUMatrix) where {T}
    m,n = size(x)
    m==n || throw(DimensionMismatch("multiplicative identity defined only for square matrices"))
    I = similar(x, T)
    fill!(I, zero(T))
    gpu_call(identity_kernel, I, m, unit; elements=m)
    I
end

Base.one(x::AbstractGPUMatrix{T}) where {T} = _one(one(T), x)
Base.oneunit(x::AbstractGPUMatrix{T}) where {T} = _one(oneunit(T), x)