text:
Base.promote_array_type{FSA <: FixedArray, T}(F, ::Type{T}, ::Type{FSA}) = FSA

# operations
const unaryOps = (:-, :~, :conj, :abs,
                  :sin, :cos, :tan, :sinh, :cosh, :tanh,
                  :asin, :acos, :atan, :asinh, :acosh, :atanh,
                  :sec, :csc, :cot, :asec, :acsc, :acot,
                  :sech, :csch, :coth, :asech, :acsch, :acoth,
                  :sinc, :cosc, :cosd, :cotd, :cscd, :secd,
                  :sind, :tand, :acosd, :acotd, :acscd, :asecd,
                  :asind, :atand, :rad2deg, :deg2rad,
                  :log, :log2, :log10, :log1p, :exponent, :exp,
                  :exp2, :expm1, :cbrt, :sqrt, :erf,
                  :erfc, :erfcx, :erfi, :dawson, :ceil, :floor,
                  :trunc, :round, :significand, :lgamma,
                  :gamma, :lfact, :frexp, :modf, :airy, :airyai,
                  :airyprime, :airyaiprime, :airybi, :airybiprime,
                  :besselj0, :besselj1, :bessely0, :bessely1,
                  :eta, :zeta, :digamma)

# vec-vec and vec-scalar
const binaryOps = (:.+, :.-,:.*, :./, :.\, :.^,
                   :.==, :.!=, :.<, :.<=, :.>, :.>=, :+, :-,
                   :min, :max,
                   :div, :fld, :rem, :mod, :mod1, :cmp,
                   :atan2, :besselj, :bessely, :hankelh1, :hankelh2,
                   :besseli, :besselk, :beta, :lbeta)

const matrixOps = (:*, :/)

const reductions = ((:sum,:+), (:prod,:*), (:minimum,:min), (:maximum,:max))

function gen_functor(func::Symbol, unary::Int)
    functor_name  = gensym()
    arguments     = ntuple(i->symbol("arg$i"), unary)
    functor_expr  = quote
        immutable $functor_name <: Func{$unary} end
        @inline call(::$functor_name, $(arguments...)) = $func($(arguments...))
    end
    return (functor_name, functor_expr)
end

for (callfun, reducefun) in reductions
    functor_name, functor_expr = gen_functor(reducefun, 2)
    eval(quote
        $functor_expr
        @inline $(callfun){T <: FixedArray}(x::T) = reduce($functor_name(), x)
    end)
end

for op in unaryOps
    functor_name, functor_expr = gen_functor(op, 1)
    eval(quote
        $functor_expr
        @inline $(op){T <: FixedArray}(x::T) = map($functor_name(), x)
    end)
end

for op in binaryOps
    functor_name, functor_expr = gen_functor(op, 2)
    eval(quote
        $functor_expr
        @inline $op{T <: FixedArray}(x::T, y::T) = map($functor_name(), x, y)
        @inline $op{T, T2, NDIM, SIZE}(x::FixedArray{T, NDIM, SIZE}, y::FixedArray{T2, NDIM, SIZE}) = $op(promote(x, y)...)
        @inline $op{T <: Number}(x::T, y::FixedArray{T}) = map($functor_name(), x, y)
        @inline $op{T1 <: Number, T2}(x::T1, y::FixedArray{T2}) = $op(promote(x, y)...)
        @inline $op{T <: Number}(x::FixedArray{T}, y::T) = map($functor_name(), x, y)
        @inline $op{T1, T2 <: Number}(x::FixedArray{T1}, y::T2) = $op(promote(x, y)...)
    end)
end

for op in matrixOps
    functor_name, functor_expr = gen_functor(op, 2)
    eval(quote
        $functor_expr
        @inline $op{T <: Number}(x::T, y::FixedArray{T}) = map($functor_name(), x, y)
        @inline $op{T1 <: Number, T2}(x::T1, y::FixedArray{T2}) = $op(promote(x, y)...)
        @inline $op{T <: Number}(x::FixedArray{T}, y::T) = map($functor_name(), x, y)
        @inline $op{T1, T2 <: Number}(x::FixedArray{T1}, y::T2) = $op(promote(x, y)...)
    end)
end

@inline function promote{T1 <: FixedArray, T2 <: FixedArray}(a::T1, b::T2)
    T = promote_type(eltype(T1), eltype(T2))
    map(T, a), map(T, b)
end
function promote{T1, T2 <: Number}(a::FixedArray{T1}, b::T2)
    T = promote_type(T1, T2)
    map(T, a), T(b)
end
function promote{T1 <: Number, T2}(a::T1, b::FixedArray{T2})
    T = promote_type(T1, T2)
    T(a), map(T, b)
end

function promote_rule{N, T, X<:Number}(::Type{Vec{N,T}}, ::Type{X})
    Vec{N, promote_type(T, X)}
end

@inline ctranspose{R, C, T}(a::Mat{R, C, T}) = Mat(ntuple(CRowFunctor(a), Val{R}))
@generated function ctranspose{N,T}(b::Vec{N,T})
    expr = [:(b._[$i]',) for i=1:N]
    return quote
        Mat{1,N,T}($(expr...))
    end
end
@inline transpose{R, C, T}(a::Mat{R, C, T}) = Mat(ntuple(RowFunctor(a), Val{R}))
@generated function transpose{N,T}(b::Vec{N,T})
    expr = [:(transpose(b._[$i]),) for i=1:N]
    return quote
        Mat{1,N,T}($(expr...))
    end
end
@inline Base.hypot{T}(v::FixedVector{2,T}) = hypot(v[1],v[2])

immutable DotFunctor <: Func{2} end
call(::DotFunctor, a, b) = a'*b
@inline dot{T <:  FixedArray}(a::T, b::T) = sum(map(DotFunctor(), a, b))

immutable BilinearDotFunctor <: Func{2} end
call(::BilinearDotFunctor, a, b) = a*b
@inline bilindot{T <: Union{FixedArray, Tuple}}(a::T, b::T) = sum(map(DotFunctor(), a, b))

@inline bilindot{T}(a::NTuple{1,T}, b::NTuple{1,T}) = @inbounds return a[1]*b[1]
@inline bilindot{T}(a::NTuple{2,T}, b::NTuple{2,T}) = @inbounds return (a[1]*b[1] + a[2]*b[2])
@inline bilindot{T}(a::NTuple{3,T}, b::NTuple{3,T}) = @inbounds return (a[1]*b[1] + a[2]*b[2] + a[3]*b[3])
@inline bilindot{T}(a::NTuple{4,T}, b::NTuple{4,T}) = @inbounds return (a[1]*b[1] + a[2]*b[2] + a[3]*b[3]+a[4]*b[4])

#cross{T}(a::FixedVector{2, T}, b::FixedVector{2, T}) = a[1]*b[2]-a[2]*b[1] # not really used!?
@inline cross{T<:Real}(a::FixedVector{3, T}, b::FixedVector{3, T}) = @inbounds return typeof(a)(
    a[2]*b[3]-a[3]*b[2],
    a[3]*b[1]-a[1]*b[3],
    a[1]*b[2]-a[2]*b[1]
)

@inline norm{T, N}(a::FixedVector{T, N})     = sqrt(dot(a,a))
@inline normalize{FSA <: FixedArray}(a::FSA) = a / norm(a)

#Matrix
@inline det{T}(A::FixedMatrix{1, 1, T}) = @inbounds return ( A[1] )
@inline det{T}(A::FixedMatrix{2, 2, T}) = @inbounds return ( A[1,1]*A[2,2] - A[1,2]*A[2,1])
@inline det{T}(A::FixedMatrix{3, 3, T}) = @inbounds return (
    A[1,1]*(A[2,2]*A[3,3]-A[2,3]*A[3,2]) -
    A[1,2]*(A[2,1]*A[3,3]-A[2,3]*A[3,1]) +
    A[1,3]*(A[2,1]*A[3,2]-A[2,2]*A[3,1])
)
@inline det{T}(A::FixedMatrix{4, 4, T}) = @inbounds return (
    A[13] * A[10]  * A[7]  * A[4]  - A[9] * A[14] * A[7]  * A[4]   -
    A[13] * A[6]   * A[11] * A[4]  + A[5] * A[14] * A[11] * A[4]   +
    A[9]  * A[6]   * A[15] * A[4]  - A[5] * A[10] * A[15] * A[4]   -
    A[13] * A[10]  * A[3]  * A[8]  + A[9] * A[14] * A[3]  * A[8]   +
    A[13] * A[2]   * A[11] * A[8]  - A[1] * A[14] * A[11] * A[8]   -
    A[9]  * A[2]   * A[15] * A[8]  + A[1] * A[10] * A[15] * A[8]   +
    A[13] * A[6]   * A[3]  * A[12] - A[5] * A[14] * A[3]  * A[12]  -
    A[13] * A[2]   * A[7]  * A[12] + A[1] * A[14] * A[7]  * A[12]  +
    A[5]  * A[2]   * A[15] * A[12] - A[1] * A[6]  * A[15] * A[12]  -
    A[9]  * A[6]   * A[3]  * A[16] + A[5] * A[10] * A[3]  * A[16]  +
    A[9]  * A[2]   * A[7]  * A[16] - A[1] * A[10] * A[7]  * A[16]  -
    A[5]  * A[2]   * A[11] * A[16] + A[1] * A[6]  * A[11] * A[16]
)

trace(A::FixedMatrix{1,1}) = A[1,1]
trace(A::FixedMatrix{2,2}) = A[1,1] + A[2,2]
trace(A::FixedMatrix{3,3}) = A[1,1] + A[2,2] + A[3,3]
trace(A::FixedMatrix{4,4}) = A[1,1] + A[2,2] + A[3,3] + A[4,4]

\{m,n,T}(mat::Mat{m,n,T}, v::Vec{n, T}) = inv(mat)*v

@inline inv{T}(A::Mat{1, 1, T}) = @inbounds return Mat{1, 1, T}(inv(A[1]))
@inline function inv{T}(A::Mat{2, 2, T})
    determinant = det(A)
    @inbounds return Mat{2, 2, T}(
        (A[2,2] /determinant, -A[2,1]/determinant),
        (-A[1,2]/determinant, A[1,1] /determinant)
    )
end
@inline function inv{T}(A::Mat{3, 3, T})
    determinant = det(A)
    @inbounds return  Mat{3, 3, T}(
        ((A[2,2]*A[3,3]-A[2,3]*A[3,2]) /determinant,
        -(A[2,1]*A[3,3]-A[2,3]*A[3,1])/determinant,
        (A[2,1]*A[3,2]-A[2,2]*A[3,1]) /determinant),

        (-(A[1,2]*A[3,3]-A[1,3]*A[3,2])/determinant,
        (A[1,1]*A[3,3]-A[1,3]*A[3,1]) /determinant,
        -(A[1,1]*A[3,2]-A[1,2]*A[3,1])/determinant),

        ((A[1,2]*A[2,3]-A[1,3]*A[2,2]) /determinant,
        -(A[1,1]*A[2,3]-A[1,3]*A[2,1])/determinant,
        (A[1,1]*A[2,2]-A[1,2]*A[2,1]) /determinant)
    )
end

@inline function inv{T}(A::Mat{4, 4, T})
    determinant = det(A)
    @inbounds return Mat{4, 4, T}(
        ((A[2,3]*A[3,4]*A[4,2] - A[2,4]*A[3,3]*A[4,2] + A[2,4]*A[3,2]*A[4,3] - A[2,2]*A[3,4]*A[4,3] - A[2,3]*A[3,2]*A[4,4] + A[2,2]*A[3,3]*A[4,4]) / determinant,
        (A[2,4]*A[3,3]*A[4,1] - A[2,3]*A[3,4]*A[4,1] - A[2,4]*A[3,1]*A[4,3] + A[2,1]*A[3,4]*A[4,3] + A[2,3]*A[3,1]*A[4,4] - A[2,1]*A[3,3]*A[4,4]) / determinant,
        (A[2,2]*A[3,4]*A[4,1] - A[2,4]*A[3,2]*A[4,1] + A[2,4]*A[3,1]*A[4,2] - A[2,1]*A[3,4]*A[4,2] - A[2,2]*A[3,1]*A[4,4] + A[2,1]*A[3,2]*A[4,4]) / determinant,
        (A[2,3]*A[3,2]*A[4,1] - A[2,2]*A[3,3]*A[4,1] - A[2,3]*A[3,1]*A[4,2] + A[2,1]*A[3,3]*A[4,2] + A[2,2]*A[3,1]*A[4,3] - A[2,1]*A[3,2]*A[4,3]) / determinant),

        ((A[1,4]*A[3,3]*A[4,2] - A[1,3]*A[3,4]*A[4,2] - A[1,4]*A[3,2]*A[4,3] + A[1,2]*A[3,4]*A[4,3] + A[1,3]*A[3,2]*A[4,4] - A[1,2]*A[3,3]*A[4,4]) / determinant,
        (A[1,3]*A[3,4]*A[4,1] - A[1,4]*A[3,3]*A[4,1] + A[1,4]*A[3,1]*A[4,3] - A[1,1]*A[3,4]*A[4,3] - A[1,3]*A[3,1]*A[4,4] + A[1,1]*A[3,3]*A[4,4]) / determinant,
        (A[1,4]*A[3,2]*A[4,1] - A[1,2]*A[3,4]*A[4,1] - A[1,4]*A[3,1]*A[4,2] + A[1,1]*A[3,4]*A[4,2] + A[1,2]*A[3,1]*A[4,4] - A[1,1]*A[3,2]*A[4,4]) / determinant,
        (A[1,2]*A[3,3]*A[4,1] - A[1,3]*A[3,2]*A[4,1] + A[1,3]*A[3,1]*A[4,2] - A[1,1]*A[3,3]*A[4,2] - A[1,2]*A[3,1]*A[4,3] + A[1,1]*A[3,2]*A[4,3]) / determinant),

        ((A[1,3]*A[2,4]*A[4,2] - A[1,4]*A[2,3]*A[4,2] + A[1,4]*A[2,2]*A[4,3] - A[1,2]*A[2,4]*A[4,3] - A[1,3]*A[2,2]*A[4,4] + A[1,2]*A[2,3]*A[4,4]) / determinant,
        (A[1,4]*A[2,3]*A[4,1] - A[1,3]*A[2,4]*A[4,1] - A[1,4]*A[2,1]*A[4,3] + A[1,1]*A[2,4]*A[4,3] + A[1,3]*A[2,1]*A[4,4] - A[1,1]*A[2,3]*A[4,4]) / determinant,
        (A[1,2]*A[2,4]*A[4,1] - A[1,4]*A[2,2]*A[4,1] + A[1,4]*A[2,1]*A[4,2] - A[1,1]*A[2,4]*A[4,2] - A[1,2]*A[2,1]*A[4,4] + A[1,1]*A[2,2]*A[4,4]) / determinant,
        (A[1,3]*A[2,2]*A[4,1] - A[1,2]*A[2,3]*A[4,1] - A[1,3]*A[2,1]*A[4,2] + A[1,1]*A[2,3]*A[4,2] + A[1,2]*A[2,1]*A[4,3] - A[1,1]*A[2,2]*A[4,3]) / determinant),

        ((A[1,4]*A[2,3]*A[3,2] - A[1,3]*A[2,4]*A[3,2] - A[1,4]*A[2,2]*A[3,3] + A[1,2]*A[2,4]*A[3,3] + A[1,3]*A[2,2]*A[3,4] - A[1,2]*A[2,3]*A[3,4]) / determinant,
        (A[1,3]*A[2,4]*A[3,1] - A[1,4]*A[2,3]*A[3,1] + A[1,4]*A[2,1]*A[3,3] - A[1,1]*A[2,4]*A[3,3] - A[1,3]*A[2,1]*A[3,4] + A[1,1]*A[2,3]*A[3,4]) / determinant,
        (A[1,4]*A[2,2]*A[3,1] - A[1,2]*A[2,4]*A[3,1] - A[1,4]*A[2,1]*A[3,2] + A[1,1]*A[2,4]*A[3,2] + A[1,2]*A[2,1]*A[3,4] - A[1,1]*A[2,2]*A[3,4]) / determinant,
        (A[1,2]*A[2,3]*A[3,1] - A[1,3]*A[2,2]*A[3,1] + A[1,3]*A[2,1]*A[3,2] - A[1,1]*A[2,3]*A[3,2] - A[1,2]*A[2,1]*A[3,3] + A[1,1]*A[2,2]*A[3,3]) / determinant)
    )
end

# Matrix
(*){T, M, N, O, K}(a::FixedMatrix{M, N, T}, b::FixedMatrix{O, K, T}) = throw(DimensionMismatch("$N != $O in $(typeof(a)) and $(typeof(b))"))
(*){T, M, N, O}(a::FixedMatrix{M, N, T}, b::FixedVector{O, T}) = throw(DimensionMismatch("$N != $O in $(typeof(a)) and $(typeof(b))"))

@generated function *{T, N}(a::FixedVector{N, T}, b::FixedMatrix{1, N, T})
    expr = Expr(:tuple, [Expr(:tuple, [:(a[$i] * b[$j]) for i in 1:N]...) for j in 1:N]...)
    :( Mat($(expr)) )
end

@generated function *{T, M, N}(a::Mat{M, N, T}, b::Vec{N,T})
    expr = [:(bilindot(row(a, $i), b.(1))) for i=1:M]
    :( Vec($(expr...)) )
end
@generated function *{T, M, N, R}(a::Mat{M, N, T}, b::Mat{N, R, T})
    expr = Expr(:tuple, [Expr(:tuple, [:(bilindot(row(a, $i), column(b,$j))) for i in 1:M]...) for j in 1:R]...)
    :( Mat($(expr)) )
end

function (==)(a::FixedVectorNoTuple, b::FixedVectorNoTuple)
    s_a = size(a)
    s_b = size(b)
    s_a == s_b || return false
    @inbounds for i = 1:length(a)
        a[i] == b[i] || return false
    end
    true
end
(==)(a::FixedArray, b::FixedArray) = a.(1) == b.(1)

(==){R, T, FSA <: FixedVector}(a::FSA, b::Mat{R, 1, T}) = a.(1) == column(b,1)
(==){R, T, FSA <: FixedVector}(a::Mat{R, 1, T}, b::FSA) = column(a,1) == b.(1)
function (==)(a::FixedArray, b::AbstractArray)
    s_a = size(a)
    s_b = size(b)
    s_a == s_b || return false
    @inbounds for i = 1:length(a)
        a[i] == b[i] || return false
    end
    true
end

(==)(a::AbstractArray, b::FixedArray) = b == a

# To support @test_approx_eq
Base.Test.approx_full(a::FixedArray) = a

# UniformScaling

*(J::Base.LinAlg.UniformScaling, A::FixedArray) = J.λ*A
*(A::FixedArray, J::Base.LinAlg.UniformScaling) = A*J.λ
/(A::FixedArray, J::Base.LinAlg.UniformScaling) = A/J.λ

+{m, n, T}(A::Mat{m,n, T}, J::Base.LinAlg.UniformScaling) = A + J.λ*eye(Mat{m,n,T})
+{m, n, T}(J::Base.LinAlg.UniformScaling, A::Mat{m,n, T}) = A + J
-{m, n, T}(A::Mat{m,n, T}, J::Base.LinAlg.UniformScaling) = A + (-J)
-{m, n, T}(J::Base.LinAlg.UniformScaling, A::Mat{m,n, T}) = J.λ*eye(Mat{m,n,T}) - A


text:
#mutable struct MapNarrowPicardGrp{S, T} <: Map{S, T, HeckeMap, MapNarrowPicardGrp}
#  header::MapHeader{S, T}
#
#  picard_group # picard group map of order
#  right_transform::fmpz_mat
#  betas # Vector of factorized algebra elements
#  gammas # the same type as betas
#
#  function MapNarrowPicardGrp{S, T}() where {S, T}
#    return new{S, T}()
#  end
#end
#
#function show(io::IO, mP::MapNarrowPicardGrp)
#  @show_name(io, mP)
#  println(io, "Narrow Picard group map of ")
#  show(IOContext(io, :compact => true), codomain(mP))
#end

################################################################################
#
#  High level functions
#
################################################################################

@doc Markdown.doc"""
      picard_group(O::NfOrd) -> GrpAbFinGen, MapClassGrp

Returns the Picard group of O and a map from the group in the set of
(invertible) ideals of O.
"""
function narrow_picard_group(O::NfOrd)
  U, mU = unit_group(O)

  # determine the totally positive units

  Q, Q_to_elem, elem_to_Q = _principal_ideals_modulo_totally_positive_principal_ideals(O, mU)

  C, mC = picard_group(O)

  @assert issnf(C)

  new_ngens = ngens(Q) + ngens(C)

  R = zero_matrix(FlintZZ, new_ngens, new_ngens)
  for i in 1:ngens(Q)
    R[i, i] = 2
  end

  idealgens = ideal_type(O)[]

  Crels = rels(C)

  for i in 1:ngens(C)
    I0 = mC(C[i])
    I = I0^(order(C[i]))
    push!(idealgens, I0)
    # I is principal
    fl, a = isprincipal(I)
    @assert fl
    q = elem_to_Q(elem_in_nf(a))
    R[ngens(Q) + i,:] = hcat(-q.coeff, Crels[i, :])
  end

  B = abelian_group(R)
  BS, BStoB = snf(B)
  BtoBS = inv(BStoB)

  # If given an element of B, I have to write it in the
  # standard generators [1,...,0], [0,...,1]

  im_Q_gens = elem_type(B)[]
  for i in 1:ngens(Q)
    v = fmpz[0 for j in 1:ngens(B)]
    v[i] = 1
    push!(im_Q_gens, B(v))
  end

  Q_to_B = hom(Q, B, im_Q_gens)

  disclog = function(J)
    d = denominator(J)
    JJ = numerator(J)
    c = mC\(JJ)
    for i in 1:ngens(C)
      JJ = JJ * idealgens[i]^(-c.coeff[i] + order(C[i]))
    end
    fl, b = isprincipal(JJ)
    @assert fl
    @assert b * O == JJ
    q = elem_to_Q(elem_in_nf(b)//d)
    return BtoBS(B(hcat(-q.coeff, c.coeff)))
  end

  _exp = function(el)
    @assert parent(el) == BS
    _el = BStoB(el)
    _elC = C([_el.coeff[1, ngens(Q) + i] for i in 1:ngens(C)])
    JJ = mC(_elC)
    _ell = _el - BStoB(disclog(fractional_ideal(JJ)))
    fl, b = haspreimage(Q_to_B, _ell)
    @assert fl
    J = Q_to_elem(b) * mC(_elC)
    return J
  end

   # A test

   _Q, = units_modulo_totally_positive_units(O)
   r, s = signature(nf(O))
   @assert order(BS) == divexact(order(C) * 2^r, order(_Q))

  return BS, disclog, _exp
end

function _principal_ideals_modulo_totally_positive_principal_ideals(O, mU)
  S, mS, h, rlp = units_modulo_totally_positive_units(O, mU)
  OK = maximal_order(O)
  R, f, g = infinite_primes_map(OK, rlp, 1 * OK)
  # First take the quotient of R by the totally positive units
  RR, mRR = quo(R, elem_type(R)[R(Int[1 for i in 1:ngens(R)])], false)
  SinRR_gen = elem_type(RR)[]
  for i in 1:ngens(S)
    push!(SinRR_gen, mRR(g(elem_in_nf(mU(mS\(S[i]))))))
  end

  Q, mQ = quo(RR, SinRR_gen, false)
  QQ, QQtoQ = snf(Q)
  gg = x -> QQtoQ\(mQ(mRR(g(x))))
  ff = y -> elem_in_nf(f(mRR\(mQ\(QQtoQ(y)))))
  return QQ, ff, gg
end

function units_modulo_totally_positive_units(O::NfOrd)
  U, mU = unit_group(O)
  return units_modulo_totally_positive_units(O, mU)
end

function units_modulo_totally_positive_units(O::NfOrd, mU)
  OK = maximal_order(O) # Will be computed anyway
  U = domain(mU)
  K = nf(O)
  r, = signature(K)
  rlp = real_places(K)

  A = abelian_group([2 for i in 1:r])
  imag = elem_type(A)[]
  for i in 1:ngens(U)
    v = A(Int[ sign(elem_in_nf(mU(U[i])), s) == -1 ? 1 : 0 for s in rlp])
    push!(imag, v)
  end

  h = hom(U, A, imag)

  K, mK = kernel(h, false)

  Q, mQ = quo(U, mK, false)
  S, StoQ = snf(Q)
  return S, mQ * inv(StoQ), h, rlp
end

ext: jl


text:
import Base.del_client
export del_client

# PATCHES FOR STUFF WHAT EFFECTS CLUSTER COMPUTING

# Representation of SharedArray that works for remotely hosted data.
function getrepresentation(S)
    buf = IOBuffer()
    td = TextDisplay(buf)
    Base.Multimedia.display(td, S)
    str = takebuf_string(buf)
    return str
end

function Base.display(S::SharedArray)
    validpid = minimum(S.pids)
    repr = @fetchfrom validpid getrepresentation(S)
    print_with_color(:bold, repr)
end

#workarounds for #14445 - basically doing remotecall can break garbage collection

if VERSION < v"0.5-"

function del_client(pg, id, client)
# As a workaround to issue https://github.com/JuliaLang/julia/issues/14445
# the dict/set updates are executed asynchronously so that they do
# not occur in the midst of a gc. The `@async` prefix must be removed once
# 14445 is fixed.
    @async begin
        rv = get(pg.refs, id, false)
        if rv != false
            delete!(rv.clientset, client)
            if isempty(rv.clientset)
                delete!(pg.refs, id)
            end
        end
    end
    nothing
end

end


text:
module SQL

import MySQL
import DBInterface
import DataFrames: DataFrame
import URIs: URI
import SQLite

export sql, sql_set_default, @sql, @sql_cmd, sql_connect

CC = nothing

function sql_connect(uri::String)
    u = URI(uri)
    if startswith(u.scheme, "mysql")
        return mysql_connect(u)
    end
    if startswith(u.scheme, "sqlite")
        return sqlite_connect(u)
    end
    throw(ArgumentError("not a known schema $(uri)"))
end

function mysql_connect(u::URI)
    user, password = split(u.userinfo, ":")
    # need to be strings not substrings :(
    user = String(user)
    password = String(password)
    host = String(u.host)
    if u.port == ""
        port = 3306
    else
        port = parse(Int, u.port)
    end
    db = String(u.path[2:end])
    query = filter(t -> t[1] == "unix_socket", [split(q, "=") for q in split(u.query, "&")])
    if length(query) >= 1
        unix_socket = String(query[1][2])
    else
        unix_socket = MySQL.API.MYSQL_DEFAULT_SOCKET
    end
    @assert startswith(u.scheme, "mysql")
    global CC
    CC = DBInterface.connect(MySQL.Connection, host, user, password; port = port, db = db, unix_socket = unix_socket)
    CC
end

function sqlite_connect(u::URI)
    @assert startswith(u.scheme, "sqlite") "not an sqlite database"
    db = String(u.path[2:end])
    global CC
    CC = DBInterface.connect(SQLite.DB, db)
    CC
end

"Set default database connection to conn or rest to nothing"
function sql_set_default(conn::Union{DBInterface.Connection,Nothing})::Union{DBInterface.Connection,Nothing}
    global CC
    old = CC
    CC = conn
    old
end
"Execute an SQL query on database conn"
function sql_df(conn, query::String)::DataFrame
    @assert conn !== nothing "no database connection!"
    DataFrame(DBInterface.execute(DBInterface.prepare(conn, query)))
end

function sql(conn, query::String)::DataFrame
    sql_df(conn, query)
end

function sql(query::String)::DataFrame
    sql_df(query)
end

function sql_df(conn::String, query::String)::DataFrame
    conn = mysql_connect(conn)
    sql_df(conn, query)
end

function sql_df(query::String)
    sql_df(CC, query)
end

macro sql_cmd(query::String)

    return quote
        local q0 = @eval $query

        sql_df(CC, q0)
    end
end

macro sql(conn, query)
    return quote
        local c = @eval $conn
        local q2 = @eval $query
        sql_df(c, q2)
    end
end
macro sql(query)
    return quote
        local q = @eval $query
        sql_df(CC, q)
    end
end

end # module


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

using MXNet

#--------------------------------------------------------------------------------
# define lenet with stn layer

# input
data = mx.Variable(:data)

# the localisation network in lenet-stn
# it will increase acc about more than 1%, when num-epoch >=15
# The localization net just takes the data as input and must output a vector in R^n
loc_net = @mx.chain mx.Convolution(data, num_filter=10, kernel=(5, 5), stride=(2,2)) =>
                    mx.Activation(act_type=:relu) =>
                    mx.Pooling( kernel=(2, 2), stride=(2, 2), pool_type=:max) =>
                    mx.Convolution( num_filter=10, kernel=(3, 3), stride=(2,2), pad=(1, 1)) =>
                    mx.Activation(act_type=:relu) =>
                    mx.Pooling( global_pool=true, kernel=(2, 2), pool_type=:avg) =>
                    mx.Flatten() =>
                    mx.FullyConnected(num_hidden=6, name=:stn_loc)

data=mx.SpatialTransformer(data,loc_net, target_shape = (28,28), transform_type="affine", sampler_type="bilinear")

# first conv
conv1 = @mx.chain mx.Convolution(data, kernel=(5,5), num_filter=20)  =>
                  mx.Activation(act_type=:tanh) =>
                  mx.Pooling(pool_type=:max, kernel=(2,2), stride=(2,2))

# second conv
conv2 = @mx.chain mx.Convolution(conv1, kernel=(5,5), num_filter=50) =>
                  mx.Activation(act_type=:tanh) =>
                  mx.Pooling(pool_type=:max, kernel=(2,2), stride=(2,2))

# first fully-connected
fc1   = @mx.chain mx.Flatten(conv2) =>
                  mx.FullyConnected(num_hidden=500) =>
                  mx.Activation(act_type=:tanh)

# second fully-connected
fc2   = mx.FullyConnected(fc1, num_hidden=10)

# softmax loss
lenet = mx.SoftmaxOutput(fc2, name=:softmax)

#--------------------------------------------------------------------------------

# load data
batch_size = 100
include("mnist-data.jl")
train_provider, eval_provider = get_mnist_providers(batch_size; flat=false)

#--------------------------------------------------------------------------------
# fit model
model = mx.FeedForward(lenet, context=mx.cpu())

# optimizer
optimizer = mx.ADAM(η=0.01, λ=0.00001)

# fit parameters
initializer=mx.XavierInitializer(distribution = mx.xv_uniform, regularization = mx.xv_avg, magnitude = 1)
mx.fit(model, optimizer, train_provider, n_epoch=20, eval_data=eval_provider,initializer=initializer)

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet


text:
##############################################################################
#
# LDA
#
##############################################################################

module Lda

mutable struct TopicBasedDocument
    topic::Vector{Int}
    text::Vector{Int}
    topicidcount::Vector{Int}
end
TopicBasedDocument(ntopics) = TopicBasedDocument(Vector{Int}(), Vector{Int}(), zeros(Int, ntopics))

mutable struct Topic
    count::Int
    wordcount::Dict{Int, Int}
end
Topic() = Topic(0, Dict{Int, Int}())

end

"""
    ϕ, θ = lda(dtm::DocumentTermMatrix, ntopics::Int, iterations::Int, α::Float64, β::Float64; kwargs...)

Perform [Latent Dirichlet allocation](https://en.wikipedia.org/wiki/Latent_Dirichlet_allocation).

# Required Positional Arguments
- `α` Dirichlet dist. hyperparameter for topic distribution per document. `α<1` yields a sparse topic mixture for each document. `α>1` yields a more uniform topic mixture for each document.
- `β` Dirichlet dist. hyperparameter for word distribution per topic. `β<1` yields a sparse word mixture for each topic. `β>1` yields a more uniform word mixture for each topic.

# Optional Keyword Arguments
- `showprogress::Bool`. Show a progress bar during the Gibbs sampling. Default value: `true`.

# Return Values
- `ϕ`: `ntopics × nwords` Sparse matrix of probabilities s.t. `sum(ϕ, 1) == 1`
- `θ`: `ntopics × ndocs` Dense matrix of probabilities s.t. `sum(θ, 1) == 1`
"""
function lda(dtm::DocumentTermMatrix, ntopics::Int, iteration::Int,
             alpha::Float64, beta::Float64; showprogress::Bool = true)

    number_of_documents, number_of_words = size(dtm.dtm)
    docs = [Lda.TopicBasedDocument(ntopics) for _ in 1:number_of_documents]
    topics = Vector{Lda.Topic}(undef, ntopics)
    for i in 1:ntopics
        topics[i] = Lda.Topic()
    end

    for wordid in 1:number_of_words
        nzdocs_idxs = nzrange(dtm.dtm, wordid)
        for docid in dtm.dtm.rowval[nzdocs_idxs]
            for _ in 1:dtm.dtm[docid, wordid]
                topicid = rand(1:ntopics)
                update_target_topic = topics[topicid]
                update_target_topic.count += 1
                update_target_topic.wordcount[wordid] = get(update_target_topic.wordcount, wordid, 0) + 1
                topics[topicid] = update_target_topic
                topic_base_document = docs[docid]
                push!(topic_base_document.topic, topicid)
                push!(topic_base_document.text, wordid)
                topic_base_document.topicidcount[topicid] += 1
            end
        end
    end

    probs = Vector{Float64}(undef, ntopics)

    wait_time = showprogress ? 1.0 : Inf

    # Gibbs sampling
    @showprogress wait_time for _ in 1:iteration
        for doc in docs
            for (i, word) in enumerate(doc.text)
                topicid_current = doc.topic[i]
                doc.topicidcount[topicid_current] -= 1
                topics[topicid_current].count -= 1
                topics[topicid_current].wordcount[word] -= 1
                document_lenth = length(doc.text) - 1

                for target_topicid in 1:ntopics
                    topicprob = (doc.topicidcount[target_topicid] + beta) / (document_lenth + beta * ntopics)
                    wordprob = (get(topics[target_topicid].wordcount, word, 0)+ alpha) / (topics[target_topicid].count + alpha * number_of_words)
                    probs[target_topicid] = topicprob * wordprob
                end
                normalize_probs = sum(probs)

                # select new topic
                select = rand()
                sum_of_prob = 0.0
                new_topicid = 0
                for (selected_topicid, prob) in enumerate(probs)
                    sum_of_prob += prob / normalize_probs
                    if select < sum_of_prob
                        new_topicid = selected_topicid
                        break
                    end
                end
                doc.topic[i] = new_topicid
                doc.topicidcount[new_topicid] = get(doc.topicidcount, new_topicid, 0) + 1
                topics[new_topicid].count += 1
                topics[new_topicid].wordcount[word] = get(topics[new_topicid].wordcount, word, 0) + 1
            end
        end
    end

    # ϕ
    # topic x word sparse matrix.
    ϕ = spzeros(ntopics, number_of_words)
    θ = getfield.(docs, :topicidcount)
    θ = Float64.(hcat(θ...))
    θ ./= sum(θ, dims=1)
    for topic in 1:ntopics
        t = topics[topic]
        for (word, count) in t.wordcount
            if 0 < t.count
                ϕ[topic, word] = count / t.count
            end
        end
    end
    return ϕ, θ
end


text:
using CompScienceMeshes
using SauterSchwabQuadrature

pI = point(1,5,3)
pII = point(2,5,3)
pIII = point(7,1,0)

Sourcechart = simplex(pI, pII, pIII)
Testchart = simplex(pII, pIII, pI)

Accuracy = 12

function integrand(x,y)
			return(((x-pI)'*(y-pII))*exp(-im*1*norm(x-y))/(4pi*norm(x-y)))
end

result = sauterschwabintegral(Sourcechart, Testchart, integrand, Accuracy, Accuracy)
println(result)

#=For those who want to test the sauterschwab_nonparameterized() function,
may uncomment the following five lines=#

#sourcechart = simplex(pI,pII,pIII)
#testchart = simplex(pI,pII,pIII)
#cf = CommonFace(Accuracy)
#result2 = sauterschwab_nonparameterized(sourcechart, testchart, integrand, cf)
#println(result2)

#=sourcechart = testchart, hence the required condition is fulfilled.
The user may also compare the two results and see that both are equal=#


text:
using GeoData, Test, Statistics, Dates, Plots
import NCDatasets, ArchGDAL
using GeoData: name, mode, window, DiskStack, bounds

testpath = joinpath(dirname(pathof(GeoData)), "../test/")
include(joinpath(testpath, "test_utils.jl"))
const DD = DimensionalData

maybedownload("https://raw.githubusercontent.com/rspatial/raster/master/inst/external/rlogo.grd", "rlogo.grd")
maybedownload("https://github.com/rspatial/raster/raw/master/inst/external/rlogo.gri", "rlogo.gri")
stem = joinpath(testpath, "data/rlogo")
@test isfile(stem * ".grd")
@test isfile(stem * ".gri")
path = stem * ".gri"

@testset "Grd array" begin
    grdarray = geoarray(path)

    @testset "open" begin
        @test all(open(A -> A[Y=1], grdarray) .=== grdarray[:, 1, :])
    end

    @testset "read" begin
        A = read(grdarray)
        @test A isa GeoArray
        @test parent(A) isa Array
    end

    @testset "array properties" begin
        @test grdarray isa GRDarray{Float32,3}
    end

    @testset "dimensions" begin
        @test length(val(dims(dims(grdarray), X))) == 101
        @test ndims(grdarray) == 3
        @test dims(grdarray) isa Tuple{<:X,<:Y,<:Band}
        @test refdims(grdarray) == ()
        @test bounds(grdarray) == ((0.0, 101.0), (0.0, 77.0), (1, 3))
    end

    @testset "other fields" begin
        @test missingval(grdarray) == -3.4f38
        @test metadata(grdarray) isa Metadata{:GRD}
        @test name(grdarray) == Symbol("red:green:blue")
        @test label(grdarray) == "red:green:blue"
        @test units(grdarray) == nothing
        customgrdarray = GRDarray(path; name=:test, mappedcrs=EPSG(4326));
        @test name(customgrdarray) == :test
        @test label(customgrdarray) == "test"
        @test mappedcrs(dims(customgrdarray, Y)) == EPSG(4326)
        @test mappedcrs(dims(customgrdarray, X)) == EPSG(4326)
        @test mappedcrs(customgrdarray) == EPSG(4326)
        proj = ProjString("+proj=merc +datum=WGS84")
        @test crs(dims(customgrdarray, Y)) == proj
        @test crs(dims(customgrdarray, X)) == proj
        @test crs(customgrdarray) == proj
    end

    @testset "getindex" begin
        @test grdarray[Band(1)] isa GeoArray{Float32,2}
        @test grdarray[Y(1), Band(1)] isa GeoArray{Float32,1}
        @test grdarray[X(1), Band(1)] isa GeoArray{Float32,1}
        @test grdarray[X(50), Y(30), Band(1)] == 115.0f0
        @test grdarray[1, 1, 1] == 255.0f0
        @test grdarray[Y(At(20)), X(At(20)), Band(3)] == 255.0f0
        @test grdarray[Y(Contains(60)), X(Contains(20)), Band(1)] == 255.0f0
    end

    # @testset "setindex" begin
    #     A = grdarray[:, :, :]
    #     temp = grdarray[1, 1, 1]
    #     println(temp)
    #     @test temp != 100.0f0
    #     grdarray[1, 1, 1] = 100.0f0
    #     grdarray[:, :, :] = 100.0f0
    #     @test grdarray[1, 1, 1] == 100.0f0
    #     grdarray[1, 1, 1] = temp
    #     @test grdarray[1, 1, 1] == temp
    #     println("sum: ", sum(A .- grdarray[:, :, :]))
    #     temp = grdarray[X(20), Y(10), Band(3)]
    #     println(temp)
    #     @test temp != 200.0f0
    #     grdarray[X(20), Y(10), Band(3)] = 200.0f0
    #     @test grdarray[20, 10, 3] == 200.0f0
    #     grdarray[X(20), Y(10), Band(3)] = temp
    # end

    @testset "selectors" begin
        geoA = grdarray[Y(Contains(3)), X(:), Band(1)]
        @test geoA isa GeoArray{Float32,1}
        @test grdarray[X(Contains(20)), Y(Contains(10)), Band(1)] isa Float32
    end

    @testset "conversion to GeoArray" begin
        geoA = grdarray[X(1:50), Y(1:1), Band(1)]
        @test size(geoA) == (50, 1)
        @test eltype(geoA) <: Float32
        @time geoA isa GeoArray{Float32,1}
        @test dims(geoA) isa Tuple{<:X,Y}
        @test refdims(geoA) isa Tuple{<:Band}
        @test metadata(geoA) == metadata(grdarray)
        @test missingval(geoA) == -3.4f38
        @test name(geoA) == Symbol("red:green:blue")
    end

    @testset "save" begin
        @testset "2d" begin
            filename2 = tempname() * ".gri"
            write(filename2, grdarray[Band(1)])
            saved = GeoArray(GRDarray(filename2))
            # 1 band is added again on save
            @test size(saved) == size(grdarray[Band(1:1)])
            @test data(saved) == data(grdarray[Band(1:1)])
        end

        @testset "3d with subset" begin
            geoA = GeoArray(grdarray)[1:100, 1:50, 1:2]
            filename = tempname()
            write(filename, GRDarray, geoA)
            saved = GeoArray(GRDarray(filename))
            @test size(saved) == size(geoA)
            @test refdims(saved) == ()
            @test bounds(saved) == bounds(geoA)
            @test size(saved) == size(geoA)
            @test missingval(saved) === missingval(geoA)
            @test metadata(saved) != metadata(geoA)
            @test metadata(saved)["creator"] == "GeoData.jl"
            @test all(metadata.(dims(saved)) .== metadata.(dims(geoA)))
            @test name(saved) == name(geoA)
            @test all(mode.(dims(saved)) .== mode.(dims(geoA)))
            @test dims(saved) isa typeof(dims(geoA))
            @test all(val.(dims(saved)) .== val.(dims(geoA)))
            @test all(mode.(dims(saved)) .== mode.(dims(geoA)))
            @test all(metadata.(dims(saved)) .== metadata.(dims(geoA)))
            @test dims(saved) == dims(geoA)
            @test all(data(saved) .=== data(geoA))
            @test saved isa typeof(geoA)
            @test data(saved) == data(geoA)
        end

        @testset "to netcdf" begin
            filename2 = tempname()
            span(grdarray[Band(1)])
            write(filename2, NCDarray, grdarray[Band(1)])
            saved = GeoArray(NCDarray(filename2; crs=crs(grdarray)))
            @test size(saved) == size(grdarray[Band(1)])
            @test replace_missing(saved, missingval(grdarray)) ≈ reverse(grdarray[Band(1)]; dims=Y)
            @test replace_missing(saved, missingval(grdarray)) ≈ reverse(grdarray[Band(1)]; dims=Y)
            @test index(saved, X) ≈ index(grdarray, X) .+ 0.5
            @test index(saved, Y) ≈ index(grdarray, Y) .+ 0.5
            @test bounds(saved, Y) == bounds(grdarray, Y)
            @test bounds(saved, X) == bounds(grdarray, X)
        end

        @testset "to gdal" begin
            # No Band
            gdalfilename = tempname() * ".tif"
            write(gdalfilename, GDALarray, grdarray[Band(1)])
            size(grdarray)
            gdalarray = GDALarray(gdalfilename)
            # @test convert(ProjString, crs(gdalarray)) == convert(ProjString, EPSG(4326))
            @test val(dims(gdalarray, X)) ≈ val(dims(grdarray, X))
            @test reverse(val(dims(gdalarray, Y))) ≈ val(dims(grdarray, Y))
            @test GeoArray(gdalarray) ≈ permutedims(grdarray[Band(1)], [X(), Y()])
            # 3 Bands
            gdalfilename2 = tempname() * ".tif"
            write(gdalfilename2, GDALarray, grdarray)
            gdalarray2 = GDALarray(gdalfilename2)
            @test all(GeoArray(gdalarray2) .== GeoArray(grdarray))
            @test val(dims(gdalarray2, Band)) == 1:3
        end

    end

    @testset "show" begin
        sh = sprint(show, MIME("text/plain"), grdarray)
        # Test but don't lock this down too much
        @test occursin("GRDarray", sh)
        @test occursin("Y", sh)
        @test occursin("X", sh)
        @test occursin("Band", sh)
    end

    @testset "plot" begin
        grdarray |> plot
        grdarray[Band(1)] |> plot
    end

end

@testset "Grd stack" begin
    grdstack = stack((a=path, b=path))

    @testset "read" begin
        st = read(grdstack)
        @test st isa GeoStack
        @test st.data isa NamedTuple
        @test first(st.data) isa GeoArray
        @test parent(first(st.data)) isa Array
    end

    @testset "indexing" begin
        @test grdstack[:a][Y(20), X(20), Band(3)] == 70.0f0
        @test grdstack[:a][Y([2,3]), X(40), Band(2)] == [240.0f0, 246.0f0]
    end

    @testset "child array properties" begin
        @test size(grdstack[:a]) == size(GeoArray(grdstack[:a])) == (101, 77, 3)
        @test grdstack[:a] isa GeoArray{Float32,3}
    end

    @testset "window" begin
        windowedstack = GRDstack((a=path, b=path); window=(Y(1:5), X(1:5), Band(1)))
        @test window(windowedstack) == (Y(1:5), X(1:5), Band(1))
        windowedarray = windowedstack[:a]
        @test windowedarray isa GeoArray{Float32,2}
        @test length.(dims(windowedarray)) == (5, 5)
        @test size(windowedarray) == (5, 5)
        @test windowedarray[1:3, 2:2] == reshape([255.0f0, 255.0f0, 255.0f0], 3, 1)
        @test windowedarray[1:3, 2] == [255.0f0, 255.0f0, 255.0f0]
        @test windowedarray[1, 2] == 255.0f0
        windowedstack = GRDstack((a=path, b=path); window=(Y(1:5), X(1:5), Band(1:1)))
        windowedarray = windowedstack[:b]
        @test windowedarray[1:3, 2:2, 1:1] == reshape([255.0f0, 255.0f0, 255.0f0], 3, 1, 1)
        @test windowedarray[1:3, 2:2, 1] == reshape([255.0f0, 255.0f0, 255.0f0], 3, 1)
        @test windowedarray[1:3, 2, 1] == [255.0f0, 255.0f0, 255.0f0]
        @test windowedarray[1, 2, 1] == 255.0f0
        windowedstack = GRDstack((a=path, b=path); window=(Band(1),));
        windowedarray = windowedstack[:b]
        @test windowedarray[1:3, 2:2] == reshape([255.0f0, 255.0f0, 255.0f0], 3, 1)
        @test windowedarray[1:3, 2] == [255.0f0, 255.0f0, 255.0f0]
        @test windowedarray[30, 30] == 185.0f0
        windowedarray |> plot
    end

    # Stack Constructors
    @testset "conversion to GeoStack" begin
        geostack = GeoStack(grdstack)
        @test Symbol.(Tuple(keys(grdstack))) == keys(geostack)
        smallstack = GeoStack(grdstack; keys=(:a,))
        @test keys(smallstack) == (:a,)
    end

    if VERSION > v"1.1-"
        @testset "copy" begin
            geoA = zero(GeoArray(grdstack[:a]))
            geoA
            copy!(geoA, grdstack, :a)
            # First wrap with GeoArray() here or == loads from disk for each cell.
            # we need a general way of avoiding this in all disk-based sources
            @test geoA == GeoArray(grdstack[:a])
        end
    end

    @testset "save" begin
        geoA = GeoArray(grdstack[:b])
        filename = tempname() * ".grd"
        write(filename, grdstack)
        base, ext = splitext(filename)
        filename_b = string(base, "_b", ext)
        saved = GeoArray(GRDarray(filename_b))
        @test typeof(saved) == typeof(geoA)
        @test data(saved) == data(geoA)
    end

    @testset "show" begin
        sh = sprint(show, MIME("text/plain"), grdstack)
        # Test but don't lock this down too much
        @test occursin("DiskStack", sh)
        @test occursin("GRDarray", sh)
        @test occursin("Y", sh)
        @test occursin("X", sh)
        @test occursin("Band", sh)
        @test occursin(":a", sh)
        @test occursin(":b", sh)
        @test occursin("logo.gri", sh)
    end

end

@testset "Grd series" begin
    grdseries = series([path, path], (Ti,); childtype=GRDarray, childkwargs=(mappedcrs=EPSG(4326), name=:test))
    @test GeoArray(grdseries[Ti(1)]) ==
        GeoArray(GRDarray(path; mappedcrs=EPSG(4326), name=:test))
    stacks = [DiskStack((a=path, b=path); childtype=GRDarray, childkwargs=(mappedcrs=EPSG(4326), name=:test))]

    grdseries = GeoSeries(stacks, (Ti,))
    @test grdseries[Ti(1)][:a] ==
        GeoArray(GRDarray(path; mappedcrs=EPSG(4326), name=:test))
    modified_ser = modify(Array, grdseries)
    @test typeof(modified_ser) <: GeoSeries{<:GeoStack{<:NamedTuple{(:a,:b),<:Tuple{<:GeoArray{Float32,3,<:Tuple,<:Tuple,<:Array{Float32,3}},Vararg}}}}

    @testset "read" begin
        geoseries = read(grdseries)
        @test geoseries isa GeoSeries{<:GeoStack}
        @test geoseries.data isa Vector{<:GeoStack}
        @test geoseries.data isa Vector{<:GeoStack}
        @test first(geoseries.data[1].data) isa GeoArray
    end
end

nothing


text:
#src # Differentiable programming
using Test: @test #src

# [Julia and its ecosystem provide some tools for differentiable programming](https://sinews.siam.org/Details-Page/scientific-machine-learning-how-julia-employs-differentiable-programming-to-do-it-best).
# Trixi.jl is designed to be flexible, extendable, and composable with Julia's growing ecosystem for
# scientific computing and machine learning. Thus, the ultimate goal is to have fast implementations
# that allow automatic differentiation (AD) without too much hassle for users. If some parts do not
# meet these requirements, please feel free to open an issue or propose a fix in a PR.

# In the following, we will walk through some examples demonstrating how to differentiate through
# Trixi.jl.

# ## Forward mode automatic differentiation

# Trixi integrates well with [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl)
# for forward mode AD.

# ### Computing the Jacobian
# The high-level interface to compute the Jacobian this way is [`jacobian_ad_forward`](@ref).
# First, we load the required packages and compute the Jacobian of a semidiscretization
# of the compressible Euler equations, a system of nonlinear conservation laws.

using Trixi, LinearAlgebra, Plots

equations = CompressibleEulerEquations2D(1.4)

solver = DGSEM(3, flux_central)
mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), initial_refinement_level=2, n_cells_max=10^5)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_density_wave, solver)

J = jacobian_ad_forward(semi);
size(J)
@test size(J) == (1024, 1024) #src

# Next, we compute the eigenvalues of the Jacobian.

λ = eigvals(J)
scatter(real.(λ), imag.(λ), label="central flux")

# As you can see here, the maximal real part is close to zero.

relative_maximum = maximum(real, λ) / maximum(abs, λ)
@test 3.0e-10 < relative_maximum < 8.0e-10 #src

# Interestingly, if we add dissipation by switching to the `flux_lax_friedrichs`
# at the interfaces, the maximal real part of the eigenvalues increases.

solver = DGSEM(3, flux_lax_friedrichs)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_density_wave, solver)

J = jacobian_ad_forward(semi)
λ = eigvals(J)

scatter!(real.(λ), imag.(λ), label="Lax-Friedrichs flux")

# Although the maximal real part is still somewhat small, it's larger than for
# the purely central discretization.

relative_maximum = maximum(real, λ) / maximum(abs, λ)
@test 1.0e-5 < relative_maximum < 5.0e-5 #src

# However, we should be careful when using this analysis, since the eigenvectors
# are not necessarily well-conditioned.

λ, V = eigen(J)
condition_number = cond(V)
@test 1.0e6 < condition_number < 5.0e6 #src

# In one space dimension, the situation is a bit different.

equations = CompressibleEulerEquations1D(1.4)

solver = DGSEM(3, flux_central)
mesh = TreeMesh((-1.0,), (1.0,), initial_refinement_level=6, n_cells_max=10^5)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_density_wave, solver)

J = jacobian_ad_forward(semi)

λ = eigvals(J)

scatter(real.(λ), imag.(λ), label="central flux")

# Here, the maximal real part is basically zero to machine accuracy.

relative_maximum = maximum(real, λ) / maximum(abs, λ)
@test 1.0e-17 < relative_maximum < 1.0e-15 #src

# Moreover, the eigenvectors are not as ill-conditioned as in 2D.

λ, V = eigen(J)
condition_number = cond(V)
@test 200 < condition_number < 4000 #src

# If we add dissipation, the maximal real part is still approximately zero.

solver = DGSEM(3, flux_lax_friedrichs)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_density_wave, solver)

J = jacobian_ad_forward(semi)
λ = eigvals(J)

scatter!(real.(λ), imag.(λ), label="Lax-Friedrichs flux")

# As you can see from the plot generated above, the maximal real part is still
# basically zero to machine precision.

relative_maximum = maximum(real, λ) / maximum(abs, λ)
@test -1.0e-15 < relative_maximum < 1.0e-15 #src

# Let's check the condition number of the eigenvectors.

λ, V = eigen(J)

condition_number = cond(V)
@test 70_000 < condition_number < 200_000 #src

# Note that the condition number of the eigenvector matrix increases but is
# still smaller than for the example in 2D.

# ### Computing other derivatives

# It is also possible to compute derivatives of other dependencies using AD in Trixi. For example,
# you can compute the gradient of an entropy-dissipative semidiscretization with respect to the
# ideal gas constant of the compressible Euler equations as described in the following. This example
# is also available as the elixir
# [examples/special\_elixirs/elixir\_euler\_ad.jl](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/special_elixirs/elixir_euler_ad.jl)

# First, we create a semidiscretization of the compressible Euler equations.

using Trixi, LinearAlgebra, ForwardDiff

equations = CompressibleEulerEquations2D(1.4)

"""
    initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D)

The classical isentropic vortex test case of
- Chi-Wang Shu (1997)
  Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory
  Schemes for Hyperbolic Conservation Laws
  [NASA/CR-97-206253](https://ntrs.nasa.gov/citations/19980007543)
"""
function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D)
  inicenter = SVector(0.0, 0.0) # size and strength of the vortex
  iniamplitude = 0.2

  rho = 1.0 # base flow
  v1 = 1.0
  v2 = 1.0
  vel = SVector(v1, v2)
  p = 10.0

  rt = p / rho                  # ideal gas equation
  cent = inicenter + vel*t      # advection of center
  cent = x - cent               # distance to centerpoint
  cent = SVector(-cent[2], cent[1])
  r2 = cent[1]^2 + cent[2]^2
  du = iniamplitude/(2*π)*exp(0.5*(1-r2)) # vel. perturbation
  dtemp = -(equations.gamma-1)/(2*equations.gamma*rt)*du^2 # isentrop

  rho = rho * (1+dtemp)^(1\(equations.gamma-1))
  vel = vel + du*cent
  v1, v2 = vel
  p = p * (1+dtemp)^(equations.gamma/(equations.gamma-1))

  prim = SVector(rho, v1, v2, p)
  return prim2cons(prim, equations)
end

mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), initial_refinement_level=2, n_cells_max=10^5)

solver = DGSEM(3, flux_lax_friedrichs, VolumeIntegralFluxDifferencing(flux_ranocha))

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_isentropic_vortex, solver)

u0_ode = compute_coefficients(0.0, semi)
size(u0_ode)
@test size(u0_ode) == (1024,) #src

# Next, we compute the Jacobian using `ForwardDiff.jacobian`.

J = ForwardDiff.jacobian((du_ode, γ) -> begin
    equations_inner = CompressibleEulerEquations2D(first(γ))
    semi_inner = Trixi.remake(semi, equations=equations_inner, uEltype=eltype(γ))
    Trixi.rhs!(du_ode, u0_ode, semi_inner, 0.0)
end, similar(u0_ode), [1.4]); # γ needs to be an `AbstractArray`

round.(extrema(J), sigdigits=2)
@test round.(extrema(J), sigdigits=2) == (-5.6, 5.6) #src

# Note that we create a semidiscretization `semi` at first to determine the state `u0_ode` around
# which we want to perform the linearization. Next, we wrap the RHS evaluation inside a closure
# and pass that to `ForwardDiff.jacobian`. There, we need to make sure that the internal caches
# are able to store dual numbers from ForwardDiff.jl by setting `uEltype` appropriately. A similar
# approach is used by [`jacobian_ad_forward`](@ref).

# Note that the ideal gas constant does not influence the semidiscrete rate of change of the
# density, as demonstrated by

norm(J[1:4:end])
@test norm(J[1:4:end]) == 0.0 #src

# Here, we used some knowledge about the internal memory layout of Trixi, an array of structs
# with the conserved variables as fastest-varying index in memory.

# ## Differentiating through a complete simulation

# It is also possible to differentiate through a complete simulation. As an example, let's differentiate
# the total energy of a simulation using the linear scalar advection equation with respect to the
# wave number (frequency) of the initial data.

using Trixi, OrdinaryDiffEq, ForwardDiff, Plots

function energy_at_final_time(k) # k is the wave number of the initial condition
    equations = LinearScalarAdvectionEquation2D(1.0, -0.3)
    mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), initial_refinement_level=3, n_cells_max=10^4)
    solver = DGSEM(3, flux_lax_friedrichs)
    initial_condition = (x, t, equation) -> begin
            x_trans = Trixi.x_trans_periodic_2d(x - equation.advection_velocity * t)
            return SVector(sinpi(k * sum(x_trans)))
    end
    semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                               uEltype=typeof(k))
    ode = semidiscretize(semi, (0.0, 1.0))
    sol = solve(ode, BS3(), save_everystep=false)
    Trixi.integrate(energy_total, sol.u[end], semi)
end

k_values = range(0.9, 1.1, length=101)

plot(k_values, energy_at_final_time.(k_values), label="Energy")

# You see a plot of a curve that resembles a parabola with local maximum around `k = 1.0`.
# Why's that? Well, the domain is fixed but the wave number changes. Thus, if the wave number is
# not chosen as an integer, the initial condition will not be a smooth periodic function in the
# given domain. Hence, the dissipative surface flux (`flux_lax_friedrichs` in this example)
# will introduce more dissipation. In particular, it will introduce more dissipation for "less smooth"
# initial data, corresponding to wave numbers `k` further away from integers.

# We can compute the discrete derivative of the energy at the final time with respect to the wave
# number `k` as follows.

round(ForwardDiff.derivative(energy_at_final_time, 1.0), sigdigits=2)
@test round(ForwardDiff.derivative(energy_at_final_time, 1.0), sigdigits=2) == 1.4e-5 #src

# This is rather small and we can treat it as zero in comparison to the value of this derivative at
# other wave numbers `k`.

dk_values = ForwardDiff.derivative.((energy_at_final_time,), k_values);

plot(k_values, dk_values, label="Derivative")

# If you remember basic calculus, a sufficient condition for a local maximum is that the first derivative
# vanishes and the second derivative is negative. We can also check this discretely.

second_derivative = round(ForwardDiff.derivative(
        k -> Trixi.ForwardDiff.derivative(energy_at_final_time, k), 1.0),
      sigdigits=2)
@test second_derivative ≈ -0.9 #src

# Having seen this application, let's break down what happens step by step.

function energy_at_final_time(k) # k is the wave number of the initial condition
    equations = LinearScalarAdvectionEquation2D(1.0, -0.3)
    mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), initial_refinement_level=3, n_cells_max=10^4)
    solver = DGSEM(3, flux_lax_friedrichs)
    initial_condition = (x, t, equation) -> begin
        x_trans = Trixi.x_trans_periodic_2d(x - equation.advection_velocity * t)
        return SVector(sinpi(k * sum(x_trans)))
    end
    semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                               uEltype=typeof(k))
    ode = semidiscretize(semi, (0.0, 1.0))
    sol = solve(ode, BS3(), save_everystep=false)
    Trixi.integrate(energy_total, sol.u[end], semi)
end

k = 1.0
round(ForwardDiff.derivative(energy_at_final_time, k), sigdigits=2)
@test round(ForwardDiff.derivative(energy_at_final_time, 1.0), sigdigits=2) == 1.4e-5 #src

# When calling `ForwardDiff.derivative(energy_at_final_time, k)` with `k=1.0`, ForwardDiff.jl
# will basically use the chain rule and known derivatives of existing basic functions
# to calculate the derivative of the energy at the final time with respect to the
# wave number `k` at `k0 = 1.0`. To do this, ForwardDiff.jl uses dual numbers, which
# basically store the result and its derivative w.r.t. a specified parameter at the
# same time. Thus, we need to make sure that we can treat these `ForwardDiff.Dual`
# numbers everywhere during the computation. Fortunately, generic Julia code usually
# supports these operations. The most basic problem for a developer is to ensure
# that all types are generic enough, in particular the ones of internal caches.

# The first step in this example creates some basic ingredients of our simulation.
equations = LinearScalarAdvectionEquation2D(1.0, -0.3)
mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), initial_refinement_level=3, n_cells_max=10^4)
solver = DGSEM(3, flux_lax_friedrichs);

# These do not have internal caches storing intermediate values of the numerical
# solution, so we do not need to adapt them. In fact, we could also define them
# outside of `energy_at_final_time` (but would need to take care of globals or
# wrap everything in another function).

# Next, we define the initial condition
initial_condition = (x, t, equation) -> begin
    x_trans = Trixi.x_trans_periodic_2d(x - equation.advection_velocity * t)
    return SVector(sinpi(k * sum(x_trans)))
end;
# as a closure capturing the wave number `k` passed to `energy_at_final_time`.
# If you call `energy_at_final_time(1.0)`, `k` will be a `Float64`. Thus, the
# return values of `initial_condition` will be `SVector`s of `Float64`s. When
# calculating the `ForwardDiff.derivative`, `k` will be a `ForwardDiff.Dual` number.
# Hence, the `initial_condition` will return `SVector`s of `ForwardDiff.Dual`
# numbers.

# The semidiscretization `semi` uses some internal caches to avoid repeated allocations
# and speed up the computations, e.g. for numerical fluxes at interfaces. Thus, we
# need to tell Trixi to allow `ForwardDiff.Dual` numbers in these caches. That's what
# the keyword argument `uEltype=typeof(k)` in
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    uEltype=typeof(k));

# does. This is basically the only part where you need to modify your standard Trixi
# code to enable automatic differentiation. From there on, the remaining steps
ode = semidiscretize(semi, (0.0, 1.0))
sol = solve(ode, BS3(), save_everystep=false)
round(Trixi.integrate(energy_total, sol.u[end], semi), sigdigits=5)
@test round(Trixi.integrate(energy_total, sol.u[end], semi), sigdigits=5) == 0.24986 #src

# do not need any modifications since they are sufficiently generic (and enough effort
# has been spend to allow general types inside these calls).

# ## Propagating errors using Measurements.jl

# [![Error bars by Randall Munroe](https://imgs.xkcd.com/comics/error_bars.png)](https://xkcd.com/2110/)

# Similar to AD, Trixi also allows propagating uncertainties using linear error propagation
# theory via [Measurements.jl](https://github.com/JuliaPhysics/Measurements.jl).
# As an example, let's create a system representing the linear advection equation
# in 1D with an uncertain velocity. Then, we create a semidiscretization using a
# sine wave as initial condition, solve the ODE, and plot the resulting uncertainties
# in the primitive variables.

using Trixi, OrdinaryDiffEq, Measurements, Plots, LaTeXStrings

equations = LinearScalarAdvectionEquation1D(1.0 ± 0.1)

mesh = TreeMesh((-1.0,), (1.0,), n_cells_max=10^5, initial_refinement_level=5)

solver = DGSEM(3)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
                                    solver, uEltype=Measurement{Float64})

ode = semidiscretize(semi, (0.0, 1.5))

sol = solve(ode, BS3(), save_everystep=false);

plot(sol)

# You should see a plot where small error bars are shown around the extrema
# and larger error bars are shown in the remaining parts.
# This result is in accordance with expectations. Indeed, the uncertain propagation
# speed will affect the extrema less since the local variation of the solution is
# relatively small there. In contrast, the local variation of the solution is large
# around the turning points of the sine wave, so the uncertainties will be relatively
# large there.

# All this is possible due to allowing generic types and having good abstractions
# in Julia that allow packages to work together seamlessly.

# ## Finite difference approximations

# Trixi provides the convenience function [`jacobian_fd`](@ref) to approximate the Jacobian
# via central finite differences.

using Trixi, LinearAlgebra

equations = CompressibleEulerEquations2D(1.4)

solver = DGSEM(3, flux_central)

mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), initial_refinement_level=2, n_cells_max=10^5)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_density_wave, solver)

J_fd = jacobian_fd(semi)

J_ad = jacobian_ad_forward(semi)

relative_difference = norm(J_fd - J_ad) / size(J_fd, 1)
@test relative_difference < 7.0e-7 #src

# This discrepancy is of the expected order of magnitude for central finite difference approximations.

# ## Linear systems

# When a linear PDE is discretized using a linear scheme such as a standard DG method,
# the resulting semidiscretization yields an affine ODE of the form

# ```math
# \partial_t u(t) = A u(t) + b,
# ```

# where `A` is a linear operator ("matrix") and `b` is a vector. Trixi allows you
# to obtain this linear structure in a matrix-free way by using [`linear_structure`](@ref).
# The resulting operator `A` can be used in multiplication, e.g. `mul!` from
# LinearAlgebra, converted to a sparse matrix using `sparse` from SparseArrays,
# or converted to a dense matrix using `Matrix` for detailed eigenvalue analyses.
# For example,

using Trixi, LinearAlgebra, Plots

equations = LinearScalarAdvectionEquation2D(1.0, -0.3)

solver = DGSEM(3, flux_lax_friedrichs)

mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), initial_refinement_level=2, n_cells_max=10^5)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver)

A, b = linear_structure(semi)

size(A), size(b)
@test (size(A), size(b)) == ((256, 256), (256,)) #src

# Next, we compute the eigenvalues of the linear operator.

λ = eigvals(Matrix(A))

scatter(real.(λ), imag.(λ))

# As you can see here, the maximal real part is close to machine precision.

λ = eigvals(Matrix(A))
relative_maximum = maximum(real, λ) / maximum(abs, λ)
@test relative_maximum < 1.0e-15 #src


text:
using Flux.Optimise
using Flux.Optimise: runall
using Flux.Tracker
using Test
@testset "Optimise" begin
  w = randn(10, 10)
  @testset for Opt in [ADAMW, ADAGrad, AdaMax, ADADelta, AMSGrad, NADAM, Descent, ADAM, Nesterov, RMSProp, Momentum]
    w′ = param(randn(10, 10))
    loss(x) = Flux.mse(w*x, w′*x)
    opt = Opt(0.001)
    if opt isa Descent || opt isa ADAGrad
      opt = Opt(0.1)
    end
    if opt isa ADADelta
      opt = Opt(0.9)
    end
    for t = 1: 10^5
      l = loss(rand(10))
      back!(l)
      delta = Optimise.update!(opt, w′.data, w′.grad)
      w′.data .-= delta
    end
    @test Flux.mse(w, w′) < 0.01
  end
end

@testset "Optimiser" begin
  w = randn(10, 10)
  @testset for Opt in [InvDecay, WeightDecay, ExpDecay]
    w′ = param(randn(10, 10))
    loss(x) = Flux.mse(w*x, w′*x)
    opt = Optimiser(Opt(), ADAM(0.001))
    for t = 1:10^5
      l = loss(rand(10))
      back!(l)
      delta = Optimise.update!(opt, w′.data, w′.grad)
      w′.data .-= delta
    end
    @test Flux.mse(w, w′) < 0.01
  end
end

@testset "Training Loop" begin
  i = 0
  l = param(1)

  Flux.train!(() -> (sleep(0.1); i += 1; l),
              (),
              Iterators.repeated((), 100),
              Descent(),
              cb = Flux.throttle(() -> (i > 3 && Flux.stop()), 1))

  @test 3 < i < 50

  # Test multiple callbacks
  x = 0
  fs = [() -> (), () -> x = 1]
  cbs = runall(fs)
  cbs()
  @test x == 1
end


text:
const contexts = Ptr{Void}[]

function initctx()
    cuInit(0)
    for dev = 0:devcount()-1
        p = Ptr{Void}[0]
        cuCtxCreate(p, 0, dev)
        push!(contexts, p[1])
    end
    setdevice(0)
end

function infodevices()
    count = devcount()
    if count == 0
        println("No CUDA-capable device found.")
    else
        for dev = 0:count-1
            name = devname(dev)
            cap = capability(dev)
            mem = round(Int, totalmem(dev) / (1024^2))
            info("device[$dev]: $(name), capability $(cap[1]).$(cap[2]), totalmem = $(mem) MB")
        end
    end
end

function capability(dev::Int)
    major, minor = Cint[0], Cint[0]
    cuDeviceComputeCapability(major, minor, dev)
    Int(major[1]), Int(minor[1])
end

function device()
    p = Ptr{Void}[0]
    cuCtxGetDevice(p)
    Int(p[1])
end

function devcount()
    p = Cint[0]
    cuDeviceGetCount(p)
    Int(p[1])
end

function devname(dev::Int)
    p = Array{UInt8}(256)
    cuDeviceGetName(p, length(p), dev)
    unsafe_string(pointer(p))
end

function totalmem(dev::Int)
    p = Csize_t[0]
    cuDeviceTotalMem(p, dev)
    Int(p[1])
end

setdevice(dev::Int) = cuCtxSetCurrent(contexts[dev+1])

synchronize() = cuCtxSynchronize()


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    AbstractShape

A geometrical shape in a `N`-dimensional space with coordinates of type `T`.
"""
abstract type AbstractShape{T,N} end

#------------------
# IMPLEMENTATIONS
#------------------
include("shapes/rectangle.jl")


text:
using Revise
using Plots, Plots.PlotMeasures
using Pkg
# Pkg.activate(".")
using BenchmarkTools
using RadiativeTransfer
using RadiativeTransfer.Absorption
using RadiativeTransfer.Scattering
using RadiativeTransfer.vSmartMOM
using RadiativeTransfer.SolarModel

I_modeled_all = zeros(7, 16);
I_deltas_all  = zeros(7, 16);
Q_modeled_all = zeros(7, 16);
Q_deltas_all  = zeros(7, 16);
U_modeled_all = zeros(7, 16);
U_deltas_all  = zeros(7, 16);

μ = [0.02, 0.06, 0.10, 0.16, 0.20, 0.28, 0.32, 0.40, 0.52, 0.64, 0.72, 0.84, 0.92, 0.96, 0.98, 1.00]
ϕs = collect(0.0:30.0:180.0)
τ = 0.5

include("/home/rjeyaram/RadiativeTransfer/test/benchmarks/natraj_trues.jl")
parameters = parameters_from_yaml("/home/rjeyaram/RadiativeTransfer/test/benchmarks/natraj.yaml");

parameters.spec_bands = [1e7/360.0 (1e7/360.0 + 1)]
parameters.vza = acosd.(μ)
parameters.sza = acosd.(0.2)

for ϕ_i in 1:length(ϕs)
    parameters.vaz = repeat([ϕs[ϕ_i]], 16)
    model = model_from_parameters(parameters);
    model.τ_rayl[1] .= τ

    R = vSmartMOM.rt_run(model, i_band=1)
    @show size(R)

    I_modeled_all[ϕ_i,:] = R[:,1,1]
    Q_modeled_all[ϕ_i,:] = R[:,2,1]
    U_modeled_all[ϕ_i,:] = R[:,3,1]

    I_deltas_all[ϕ_i,:] = abs.(I_trues[:,ϕ_i] - I_modeled_all[ϕ_i,:]) ./ I_trues[:,ϕ_i]
    Q_deltas_all[ϕ_i,:] = abs.(Q_trues[:,ϕ_i] - Q_modeled_all[ϕ_i,:]) ./ Q_trues[:,ϕ_i]
    U_deltas_all[ϕ_i,:] = abs.(U_trues[:,ϕ_i] - U_modeled_all[ϕ_i,:]) ./ U_trues[:,ϕ_i]
end

plot(parameters.vza, I_modeled_all[1,:])
plot!(parameters.vza, I_trues[:,1])

##

plot_phis = [1, 4, 7]
plot_list = []
gr(display_type=:inline)
for i in plot_phis
    p = plot(parameters.vza, Q_modeled_all[i,:], xlabel = "VZA", title = "Q, VAZ = $(ϕs[i])°", legend=:bottomleft, label="model", left_margin = 10mm)
    p = plot!(parameters.vza, Q_trues[:,i], label="natraj")
    push!(plot_list, p)

    p = plot(parameters.vza, Q_modeled_all[i,:] - Q_trues[:,i], xlabel = "VZA", title = "Model – Natraj", legend=:bottomright, label="diff", ylims=(-0.0005, 0.0005))

    push!(plot_list, p)

    p = plot(parameters.vza, Q_deltas_all[i,:] * 100, xlabel = "VZA", title = "% Difference", legend=false, label="diff")

    push!(plot_list, p)
end

plot(plot_list..., layout=(3,3), size=(1200,1200))


text:
using BinaryBuilder

platforms = [
  BinaryProvider.Linux(:i686, :glibc),
  BinaryProvider.Linux(:x86_64, :glibc),
  BinaryProvider.Linux(:aarch64, :glibc),
  BinaryProvider.Linux(:armv7l, :glibc),
  BinaryProvider.Linux(:powerpc64le, :glibc),
  BinaryProvider.MacOS()
]

sources = [
  "https://zlib.net/zlib-1.2.11.tar.gz" =>
  "c3e5e9fdd5004dcb542feda5ee4f0ff0744628baf8ed2dd5d66f8ca1197cb1a1",
]

script = raw"""
cd $WORKSPACE/srcdir
cd zlib-1.2.11/
./configure --prefix=/
make install -j3
"""

products = prefix -> [
  LibraryProduct(prefix,"libz"),
]

autobuild(pwd(), "zlib", platforms, sources, script, products)

# Special-case windows platforms
platforms = [
  BinaryProvider.Windows(:i686),
  BinaryProvider.Windows(:x86_64),
]
script = raw"""
cd $WORKSPACE/srcdir
cd zlib-1.2.11/
./configure --prefix=/
make install SHAREDLIB=libz.dll SHAREDLIBM=libz-1.dll SHAREDLIBV=libz-1.2.11.dll LDSHAREDLIBC='' -j3
"""
autobuild(pwd(), "zlib", platforms, sources, script, products)

max_stars_repo_name: UnofficialJuliaMirror/Keno-ZlibBuilder

max_issues_repo_name: UnofficialJuliaMirror/Keno-ZlibBuilder

max_forks_repo_name: UnofficialJuliaMirror/Keno-ZlibBuilder


text:
@testset "368.largest-divisible-subset.jl" begin
    @test largest_divisible_subset([1,2,3]) == [2, 1]
    @test largest_divisible_subset([1,2,4,8]) == [8, 4, 2, 1]
end


text:
using ProjAssistant
@quickactivate "Chemostat_Rath2017"

@time begin

    import Chemostat
    import Chemostat.MetNets
    import Chemostat.MetLP
    const Ch = Chemostat

    import Chemostat_Rath2017
    const ChR = Chemostat_Rath2017
    const Rd = ChR.RathData
    const H1 = ChR.Human1
    const HG = H1.HumanGEM
end

## ---------------------------------------------------------------------
## Description

# This script prepare the original HumanGEM model for modeling. It include some modifications to it (see comment in code) extracted from several data sources.

## ---------------------------------------------------------------------
# Prepare base base_model
# this will be the base model for all the processing

# Mat Model
@time begin
    mat_model = HG.load_humangem_mat_model()
    base_model = HG.load_humangem_raw_model(mat_model)
    base_model = MetNets.uncompressed_model(base_model)
    M, N = size(base_model)
    println("\nLoaded Mat model")
    println("Base Model: ", (M, N))
end

## ---------------------------------------------------------------------
base_model = ChR.prepare_metnet(HG, base_model; inf_medium = true);

## ---------------------------------------------------------------------
# Saving base_model
HG.save_humangem_base_model(base_model)


text:
mutable struct rigid_transform_t <: LCMType
    utime::Int64
    trans::SVector{3, Float64}
    quat::SVector{4, Float64}
end

@lcmtypesetup(rigid_transform_t)


text:
module DummyModule

export MyType, f

mutable struct MyType
    a::Int
end

f(x) = x^2+1

println("loaded")

end

# @everywhere include("devel/DummyModule.jl") this works


text:
using BenchmarkTools

function timefun(m, p, n, q, k=1, rep=4)
    #q = -expm1(log1p(-qs)/n) / p
    A = sprand(m, n, p)
    B = sprand(n, k, q)
    Corig = SparseWrappers.spmatmul_orig(A, B)
    beorig = @benchmark SparseWrappers.spmatmul_orig($A, $B) samples=rep
    bealt = @benchmark SparseWrappers.spmatmul_alt($A, $B) samples=rep
    torig = minimum(beorig)
    talt = minimum(bealt)
    (m=m, nnzA=nnz(A), n=n, nnzB=nnz(B), k=k, nnzC=nnz(Corig), torig=torig, talt=talt)
end

function sample(cnt::Integer, seed=0, oplimit=10^9, slimit=10^9)
    Random.seed!(seed)
    res = []
    i = 0
    while i < cnt
        m = 10^rand(3:6)
        n = 10^rand(3:6)
        k = 10^rand(0:1)
        p = exp10(rand(-6:-2))
        q = exp10(rand(-6:-1))
        nnzA = Int(ceil(m * n * p))
        nnzB = Int(ceil(n * k * q))
        nnzC = estimate_mulsize(m, nnzA, n, nnzB, k)
        space = (nnzA + nnzB + nnzC)*(sizeof(Int)+sizeof(Float64)) + (n + 2k)*(sizeof(Int))
        op1 = Int(ceil(m * n * k * p * q))
        op2 = Int(ceil(m * k))
        if n * p > 1 && n * q > 1 && op1 <= oplimit && op2 <= oplimit && space <= slimit
            i += 1
            println("$i: m=$m n=$n k=$k p=$p q=$q nnzA=$nnzA nnzB=$nnzB nnzC=$nnzC op1=$op1 op2=$op2 space=$space")
            push!(res, timefun(m, p, n, q, k))
        end
    end
    res
end

op1(r) = r.nnzA * r.nnzB / r.n
op2(r) = r.m * r.k
space(r) = r.nnzC * (sizeof(Int) + sizeof(Float64)) + r.k * sizeof(Int)

nlogn(n::Integer) = ilog2(n) * n

# convert vector of 7-tuples to array
function model_matrix(res::Vector)
    A = vcat((collect(hcat(1.0, op1(r), op2(r), r.nnzC/r.k, r.k*nlogn(r.nnzC÷r.k)) for r = res))...)
    B = vcat((collect(hcat(r.torig, r.talt) for r = res))...)
    W = [r.torig for r = res]
    A, getproperty.(B, :time), getproperty.(W, :time)
end

q(qs, n, p ) = -expm1(log1p(-qs)/n) / p
rhs(qs, k, n, p) = sprand(n, k, q(qs, n, p))
mat(n, p ) = sprandn(n, n, p)

function timefun2(matmul::Function, A::SparseMatrixCSC, qs::Float64, k::Integer)
    n = size(A, 2)
    p = nnz(A) / length(A)
    B = rhs(qs, k, n, p)
    be = @benchmark $matmul($A, $B) evals=10
    minimum(be).time
end

function sample2(matmul::Function, A, qs, kvec::AbstractVector{<:Integer})
    timefun2.(matmul, Ref(A), qs, kvec)
end

function sample1(matmul::Function, n::Integer, p::Float64, kvec::AbstractVector, qsvec::AbstractVector)
    Random.seed!(0)
    A = mat(n, p)
    res = Float64[]
    for qs in qsvec
        append!(res, sample2(matmul, A, qs, kvec))
    end
    reshape(res, length(kvec), length(qsvec))
end

#test environment to find speed of sorting vs. scanning

function sort_or_scan(rowval::Vector, nzval::Vector, xb::Vector, len::Integer, m::Integer, nnz::Integer)

    len >= m >= nnz || throw(ArgumentError("len($len) >= m($m) >= nnz($nnz) requested"))
    println("sort_or_scan($len, $m, $nnz)")
    rowval, nzval, xb, ip0, ip, m = nzsetup(rowval, nzval, xb, len, m, nnz)

    bescan = @benchmark runscan(rowval, nzval, xb, $ip0, $ip, $m) setup=begin
                                rowval, nzval, xb = copy($rowval), copy($nzval), copy($xb)
                            end samples=10 evals=1

    besort = @benchmark runsort(rowval, nzval, xb, $ip0, $ip, $m) setup=begin
                                rowval, nzval, xb = copy($rowval), copy($nzval), copy($xb)
                            end samples=10 evals=1

    getproperty.(minimum.([bescan, besort]), :time)
end

function runsort(rowvalC, nzvalC, xb, ip0, ip, m)
    let ip0=ip0, ip=ip, k0 = ip0 - 1, nnzi = ip - ip0
        alg = nnzi <= 21 ? InsertionSort : QuickSort
        sort!(rowvalC, ip0, ip-1, alg, Base.Order.Forward)
        #sort!(rowvalC, ip0, ip-1, QuickSort, Base.Order.Forward)
        for vp = ip0:ip-1
            k = rowvalC[vp]
            xb[k] = false
            nzvalC[vp] = nzvalC[k+k0]
        end
    end
end

function runscan(rowvalC, nzvalC, xb, ip0, ip, m)
    let ip0=ip0, m=m, k0 = ip0 - 1
        for k = 1:m
            if xb[k]
                xb[k] = false
                rowvalC[ip0] = k
                nzvalC[ip0] = nzvalC[k+k0]
                ip0 += 1
            end
        end
    end
end

function nzsetup(rowval, nzval, xb, len, m, nnz, Tv=Float64, Ti=Int)
    @assert 0 <= nnz <= m <= len
    Random.seed!(65+len)
    length(nzval) < len && resize!(nzval, len)
    length(rowval) < len && resize!(rowval, len)
    length(xb) < m && resize!(xb, m)
    fill!(xb, false)
    ip0 = rand(1:len-m)
    k0 = ip0 - 1
    for j = 1:nnz
        k = rand(1:m)
        xb[k] = true
        nzval[k+k0] = randn(Tv)
        rowval[j+k0] = k
    end
    ip = k0 + nnz
    rowval, nzval, xb, ip0, ip, m
end

samplex(rowval, nzval, xb, m::Integer, nz::Integer) = sort_or_scan(rowval, nzval, xb, max(10^7,2^m*9÷8), 2^m, 2^nz)
samplex(rowval, nzval, xb, m::Integer) = samplex.(Ref(rowval), Ref(nzval), Ref(xb), m, max(0,m-8):m)
samplex(r::UnitRange) = samplex.(Ref(Int[]), Ref(Float64[]), Ref(Bool[]), r)

"""
    samplex_plot(r::UnitRange, res::Vector{Vector{Vector{Float64}}})

usage:
    res = samplex(0:10)
    xscanm, scannnz, sortm, sortnnz, xm, ym = samplex_sort(0:10, res)

    scatter(scanm, scannnz)

    scatter!(sortm, sortnnz)

    plot!(xm, ym)

    See also sort_or_scan.png in data directory.
"""
function samplex_plot(r::UnitRange, res::Vector{<:Vector})
    sortm = Int[]
    sortnnz = Int[]
    scanm = Int[]
    scannnz = Int[]
    j = 0
    for m = r
        cm = res[j+=1]
        k = length(cm)
        for nz = m:-1:m-k+1
            data = cm[k]
            k -= 1
            if data[1] < data[2]
                push!(scanm, m)
                push!(scannnz, nz)
            else
                push!(sortm, m)
                push!(sortnnz, nz)
            end
        end
    end

    ym = collect(1:length(res)-7)
    xm = log2.(3 * nlogn.((^).(2, ym)))

    scanm, scannnz, sortm, sortnnz, xm, ym
end

nlogn(x) = x <= 0 ? zero(x) : ilog2(x) * x

probqs(pq::AbstractFloat, n::Integer) = -expm1(log1p(-pq) * n)
probpq(qs::AbstractFloat, n::Integer) = -expm1(log1p(-qs) / n)

function benchmark01(m::Integer, n::Integer, k::Integer, nz::Integer)
    0 <= nz <= m || throw(ArgumentError("0 <= nnz($nz) <= m($m) required"))

    qs = prevfloat(nz / m) # targeted fill factor of result kill case nz == m
    pq = probpq(qs, n) # product of fill factors of A nd B
    #1 / (k * m) <= pq || throw(ArgumentError("nnz too small"))

    x = 0.5
    p = pq * (1-x) + x
    q = pq + 1 - p
    spdq = sqrt(p/q)
    spq = sqrt(pq)
    p, q = spq * spdq, spq / spdq
    A = sprandn(m, n, p)
    B = sprandn(n, k, q)
    println("nnz(A)=$(nnz(A)) nnz(b)=$(nnz(B))")
    beorig = @benchmark SparseWrappers.spmatmul_orig($A, $B) samples=10
    benew = @benchmark SparseWrappers.spmatmul($A, $B) samples=10
    minimum(beorig).time, minimum(benew).time
end


text:
###
### show as MIME type
###
using FileIO: @format_str
using Base64: Base64EncodePipe
import Base.showable

# This is used by IJulia (for example) to display images

# showable to PNG if 2D colorant array
# issue #12 -- 0 length image isn't showable
showable(::MIME"image/png", img::AbstractMatrix{C}) where {C<:Colorant} = _length1(img) > 0

# Colors.jl turns on SVG display of colors, which leads to poor
# performance and weird spacing if you're displaying images. We need
# to disable that here.
# See https://github.com/JuliaLang/IJulia.jl/issues/229 and Images #548
showable(::MIME"image/svg+xml", img::AbstractMatrix{C}) where {C<:Color} = false

# Really large images can make display very slow, so we shrink big
# images.  Conversely, tiny images don't show up well, so in such
# cases we repeat pixels.
function Base.show(io::IO, mime::MIME"image/png", img::AbstractMatrix{C};
#                   minpixels=10^3, maxpixels=10^5,
                   minpixels=10^4, maxpixels=10^6,
                   # Jupyter seemingly can't handle 16-bit colors:
                   mapi=x->mapc(N0f8, clamp01nan(csnormalize(x)))) where C<:Colorant
    if Base.has_offset_axes(img)
        img = collect(img)
    end
    if !get(io, :full_fidelity, false)
        while _length1(img) > maxpixels
            img = _restrict1(img)  # big images
        end
        npix = _length1(img)
        if npix < minpixels
            # Tiny images
            fac = ceil(Int, sqrt(minpixels/npix))
            r = ones(Int, ndims(img))
            r[[coords_spatial(img)...]] .= fac
            img = repeat(img, inner=r)
        end
        FileIO.save(FileIO.Stream(format"PNG", io), img, mapi=mapi)
    else
        FileIO.save(FileIO.Stream(format"PNG", io), img)
    end
end

Base.show(io::IO, mime::MIME"image/png", img::OffsetArray{C}; kwargs...) where C<:Colorant =
    show(io, mime, parent(img); kwargs...)

# Not all colorspaces are supported by all backends, so reduce types to a minimum
csnormalize(c::AbstractGray) = Gray(c)
csnormalize(c::Color) = RGB(c)
csnormalize(c::Colorant) = RGBA(c)

const ColorantMatrix{T<:Colorant} = AbstractMatrix{T}

function _show_odd(io::IO, m::MIME"text/html", imgs::AbstractArray{T, 1}) where T<:ColorantMatrix
    # display a vector of images in a row
    for j in eachindex(imgs)
        write(io, "")
        show_element(IOContext(io, :thumbnail => true), imgs[j])
        write(io, "	")
    end
end

function _show_odd(io::IO, m::MIME"text/html", imgs::AbstractArray{T, N}) where {T<:ColorantMatrix, N}
    colons = ([Colon() for i=1:(N-1)]...,)
    for i in axes(imgs, N)
        write(io, "")
        _show_even(io, m, view(imgs, colons..., i)) # show even
        write(io, "	")
    end
end

function _show_even(io::IO, m::MIME"text/html", imgs::AbstractArray{T, N}, center=true) where {T<:ColorantMatrix, N}
    colons = ([Colon() for i=1:(N-1)]...,)
    centering = center ? " style='margin: auto'" : ""
    write(io, "
")
    write(io, "
")
    for i in axes(imgs, N)
        write(io, "
")
        _show_odd(io, m, view(imgs, colons..., i)) # show odd
        write(io, "
")
    end
    write(io, "
")
    write(io, "
")
end

function Base.show(io::IO, m::MIME"text/html", imgs::AbstractArray{T, N}) where {T<:ColorantMatrix, N}
    imgs = permutedims(imgs, N:-1:1)
    if N % 2 == 1
        write(io, "
")
        write(io, "
")
        write(io, "
")
        _show_odd(io, m, imgs) # Stack horizontally
        write(io, "
")
        write(io, "
")
        write(io, "
")
        if N == 1
            write(io, "
(a vector displayed as a row to save space)
")
        end
    else
        _show_even(io, m, imgs, false) # Stack vertically
    end
end

function downsize_for_thumbnail(img, w, h)
    a,b=size(img)
    a > 2w && b > 2h ?
        downsize_for_thumbnail(_restrict1(img), w, h) : img
end

function show_element(io::IOContext, img)
    io2=IOBuffer()
    w,h=get(io, :thumbnailsize, (100,100))
    im_resized = downsize_for_thumbnail(img, w, h)
    thumbnail_style = get(io, :thumbnail, false) ? "max-width: $(w)px; max-height:$(h)px;" : ""
    b64pipe=Base64EncodePipe(io2)
    write(io,"")
end

function _restrict1(img)
    if _have_restrict[]
        res = ImageTransformations.restrict(img)
    else
        @info "For better quality inline display of large images or thumbnails, load the Images package." maxlog=1
        res = img[1:2:end,1:2:end]
    end
    res
end

_length1(A::AbstractArray) = length(eachindex(A))
_length1(A) = length(A)


text:
# This file needs the source file to have run with appropriate parameters!!!
# It loads and plots the results!
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using JLD, PyPlot

ds = Float64[4, 3, 2]

a = 9.0
b = 9.8
L = 50.0
N = 250
ic = :random
prefix = "brusselator2D"
suffix = "jld2"
fnames = [savename(prefix, @strdict(a, b, d, L, N, ic), suffix) for d in ds]
fig = figure(figsize=(figx, 1.6*figy))

iplt = 0
nfile = length(fnames)
println(nfile)

for ifile = 1:nfile
    fname = fnames[ifile]
    println(fname)
    data = load(datadir("Brusselator", fname))
    uout = data["u"]
    tvec = data["t"]
    nts = length(tvec)
    for n = 1:nts
        ures = uout[n]
        iplt += 1
        ax = subplot(2,4,iplt)
        im = ax.pcolormesh(ures[1,:,:]; cmap = "inferno")
        ax.set_xticks([])
        ax.set_yticks([])
        plt.title(L"t = "*string(floor(Int,tvec[n])), fontsize = 24)
        ax.set_aspect("equal")
    end
end # of iplt loop

# plt.colorbar(im)
fig.tight_layout(pad=0.4)
add_identifiers!(fig)
wsave(plotsdir("11", "Brusselator_2D_Turing_patterns"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
@testset "fit_tensors.jl" begin

end


text:
using Interpolations # For making initial guess of template
using Polynomials

function make_blaze_and_template_1d(data::EOST; order::Integer, oversample_factor::Real = 1.0, num_knots::Integer=round(Int64,oversample_factor*num_pixels_per_order(data)),
      z_bc::AbstractVector{T3} = zeros(eltype(data.lambda),num_obs(data))
      ) where { EOST <: AbstractEchelleObservationSet, T3<:Real}
      #=
      pixel_range = 1:size(data.lambda,1)
      blaze = fit_blaze_1d(pixel_range, vec(mean(data.flux[:,order,:],dims=2)), degree = 3)
      blaze = blaze / mean(polyval(blaze,pixel_range))
      =#
      log_lambda_range = make_model_wavelength_range(data, order=order)
      blaze = make_blaze_1d(data, order=order)
      # log_lambda_range = range(log(minimum(data.lambda[1,order,:])), stop=log(maximum(data.lambda[end,order,:])), )
      template = make_template_1d(data, log_lambda_range, order, z_bc=z_bc, blaze=blaze) #pixel_range=blaze.pixel_range, blaze=blaze.model)
      #template = make_template_1d(data, log_lambda_range, order, z_bc=z_bc, pixel_range=blaze.pixel_range, blaze=blaze.model).template
      #return (pixel_range=blaze.pixel_range, blaze=blaze.model, log_lambda_range=log_lambda_range, template=template )
      return blaze, template
end

function make_blaze_and_template_2d(data::EOST;
      orders::Union{AbstractVector{T1},AbstractRange{T1}} = 1:size(data.lambda,2),
      oversample_factor::Real = 1.0, num_knots::Integer=round(Int64,oversample_factor*num_pixels_per_order(data)),
      z_bc::AbstractVector{T3} = zeros(eltype(data.lambda),num_obs(data))
      ) where { EOST <: AbstractEchelleObservationSet, T1<:Integer, T3<:Real}
   map(order->make_blaze_and_template_1d(data,order=order,z_bc=z_bc, num_knots=num_knots), orders) # TODO: TEST
end

""" make_stellar_basis(mean_spectrum; num_basis_vectors = 0)
   Returns array with one row containing the mean_spectrum followed by num_basis_vectors rows of zeros. """
function make_stellar_basis(mean_spectrum::AbstractArray{T,1}; num_basis_vectors::Integer = 0) where {T<:Real}
   @assert 1 <= length(mean_spectrum) <= max_template_length
   @assert 0 <= num_basis_vectors <= max_num_basis_vectors
   output = zeros(T,1+num_basis_vectors,length(mean_spectrum))
   output[1,:] .= mean_spectrum
   return output
end

function make_template_1d(data::EOST, log_lambda_out::AbstractRange{T1}, order::Integer;
    template::AbstractArray{T2,1} = zeros(T1,length(log_lambda_out)),
    z_bc::AbstractVector{T3} = zeros(T1,num_obs(data)),
    pixel_range::AbstractRange{T4} = 1:size(data.lambda,1),
    blaze::AbstractBlazeType = make_flat_blaze(data, pixel_range=pixel_range, order=order),
    num_star_basis_vectors::Integer = 0, num_telluric_basis_vectors::Integer = 0
    ) where { EOST <: AbstractEchelleObservationSet, T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Real}
    @assert num_obs(data) == length(z_bc)
    @assert length(log_lambda_out) == length(template)
    @assert 1<=order<=num_orders(data)
    interpolators = [ interpolate( (log.(lambda_ssb_from_obs(data.lambda[:,order,j],z_bc[j])), ),
         data.flux[:,order,j], Gridded(Linear()) ) for j in 1:num_obs(data) ]

    k_start = findfirst(k->(log_lambda_out[k]) >= log(minimum(data.lambda[1,order,:])),   1:length(log_lambda_out)  )
    if isnothing(k_start)
      k_start = 1
    end
    k_stop  = findlast( k->(log_lambda_out[k]) <= log(maximum(data.lambda[end,order,:])), 1:length(log_lambda_out)  )
    if isnothing(k_stop)
      k_stop = length(log_lambda_out)
    end
    weight = vec(sum(data.flux[:,order,:],dims=1))/size(data.flux,1) # average counts per pixel in each observation
    for k in k_start:k_stop
       weight_sum = 0.0
       for j in 1:num_obs(data)
         if interpolators[j].knots[1][1] <= log_lambda_out[k] <=  interpolators[j].knots[end][end]
            template[k] += interpolators[j](log_lambda_out[k]) # normalized flux times weight
            weight_sum  += weight[j]
         end
       end # for j over observations
       if weight_sum > 0
          # TODO: WARNING: Blaze calculation assumes length(log_lambda_out) == num_pixels_per_order
          pixel = blaze.pixel_range.start+length(blaze.pixel_range)*(k-1)/length(log_lambda_out)
          # println("k=",k," pixel=",pixel)
          b = calc_blaze(blaze,pixel) # polyval(blaze,k) #-mean(pixel_range))
          template[k] /= b * weight_sum
       end
   end # for k over log_lambda_out
   return template
end

function make_spectral_order_model(data::EOST, log_lambda_out::AbstractRange{T1}, order::Integer;
   pixel_range::AbstractRange{T4} = 1:size(data.lambda,1),
   blaze::AbstractBlazeType = make_blaze_1d(data, order=order),
   z_bc::AbstractVector{T3} = zeros(T1,num_obs(data)),
    template::AbstractArray{T2,1} = make_template_1d(data,log_lambda_out,order, z_bc=z_bc, pixel_range=pixel_range, blaze=blaze ), # zeros(T1,length(log_lambda_out)),
    num_star_basis_vectors::Integer = 0, num_telluric_basis_vectors::Integer = 0
    ) where { EOST <: AbstractEchelleObservationSet, T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Real}
    @assert num_obs(data) == length(z_bc)
    @assert length(log_lambda_out) == length(template)
    @assert 1<=order<=num_orders(data)
    # template = make_template_1d(data,log_lambda_out,order, z_bc=z_bc, pixel_range=pixel_range, blaze=blaze )
    model = SpectralOrderModel(log_lambda_out,
      star_basis = make_stellar_basis(template,num_basis_vectors=num_star_basis_vectors),
      telluric_basis = default_telluric_model(length(log_lambda_out),num_basis_vectors=num_telluric_basis_vectors),
      pixel_range=blaze.pixel_range, blaze=blaze )
   return model
#   return (pixel_range=pixel_range, blaze=blaze, log_lambda_range=log_lambda_out,template=template)

end

function make_spectral_model(data::EOST; #, log_lambda_out::AbstractRange{T1};
    pixel_range = 1:size(data.lambda,1),
    z_bc::AbstractVector{T3} = zeros(T1,num_obs(data)),
    num_star_basis_vectors::Integer = 0, num_telluric_basis_vectors::Integer = 0
    ) where { EOST <: AbstractEchelleObservationSet, T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Real}
    @assert num_obs(data) == length(z_bc)

    model = map( order->make_spectral_order_model(data,make_model_wavelength_range(data,order),order,
               pixel_range=pixel_range, z_bc=z_bc,
               num_star_basis_vectors=num_star_basis_vectors,
               num_telluric_basis_vectors=num_telluric_basis_vectors), 1:size(data.lambda,2) )
end

#=
import EchelleTools.AbstractEchelleObservationSet, EchelleTools.SpectralOrderModelAbstract
import EchelleTools.min_num_star_basis_vectors
import EchelleTools.max_num_star_basis_vectors
import EchelleTools.min_num_telluric_basis_vectors
import EchelleTools.max_num_telluric_basis_vectors
=#

function make_observations_model(data::DataT, som::AbstractArray{SOMT,1};
   orders::AbstractRange = 1:size(data.lambda,2),
   num_star_basis_vectors::Integer = 1,
   num_telluric_basis_vector::Integer = 0,
   z_bc::AbstractArray{T1,1} = zeros(num_obs),
   )  where {T1<:Real,
      DataT <: AbstractEchelleObservationSet, SOMT <: SpectralOrderModelAbstract
      }
   @assert min_num_star_basis_vectors <= num_star_basis_vectors <= max_num_star_basis_vectors
   @assert min_num_telluric_basis_vectors <= num_telluric_basis_vector <= max_num_telluric_basis_vectors
   @assert length(z_bc) == size(data.lambda,3)
   num_pixels = size(data.lambda,1)
   num_orders = size(data.lambda,2)
   num_obs = size(data.lambda,3)
   star_score = default_star(num_pixels,num_basis_vectors=num_star_basis_vectors)
   telluric_score = default_telluric_model(num_pixels,num_basis_vectors=num_telluric_basis_vector)
   normalization = zeros(num_obs)
   for i in 1:num_obs
      # TODO: specify pixel_range

      normalization[i] = mapreduce(
         order->mean( view(data.flux,:,order,i) ./
                        eval_model_baseline(view(data.lambda,:,order,i), som[order]) ),
                        +, orders) / length(orders)
   end
   ObservationsModel{eltype(star_score),eltype(telluric_score),eltype(z_bc),eltype(normalization)}(star_score,telluric_score,z_bc,normalization)
end

# Rest is no longer used?
function model_spectrum_1d(log_lambdas_template::AbstractRange{T1}, template::AbstractArray{T2,1},
      log_lambdas_out::Union{AbstractRange{T3},AbstractVector{T3}}, z_bc::Real;
      pixel_range::AbstractRange{T4} = 1:length(log_lambdas_out),
      blaze::Poly{T5} = Poly([one(T2)]),
      output::AbstractVector{T6} = zeros(T2,length(log_lambdas_out)) ) where {T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Real, T6<:Real}
   @assert length(log_lambdas_template) == length(template)
   @assert length(log_lambdas_out) == length(pixel_range) == length(output)

   fb = polyval(blaze, pixel_range )
   interpolator = scale(interpolate( template, BSpline(Linear())), log_lambdas_template)
   ll = log_lambda_ssb_from_obs.(log_lambdas_out, -z_bc)
   fm = interpolator.(ll)
   output .= fb.*fm
end

function model_spectrum_2d(log_lambdas_template::AbstractVector{ART1},
      template::AbstractArray{AVT2,1}, log_lambdas_out::AbstractArray{T3,2}, z_bc::Real;
      orders::Union{AbstractVector{T4},AbstractRange{T4}} = 1:size(log_lambdas_out,2),
      pixel_range::AbstractRange{T5} = 1:size(log_lambdas_out,1), # TODO: WARNING Assumes all orders have same pixel_range
      blaze::AbstractVector{Poly{T6}} = fill(Poly([one(eltype(eltype(AVT2)))]),size(log_lambdas_out,1)),
      output::AbstractArray{T7,2} = zeros(eltype(eltype(AVT2)),size(log_lambdas_out))
      ) where {T1<:Real, ART1<:AbstractRange{T1}, T2<:Real, AVT2<:AbstractVector{T2}, T3<:Real, T4<:Integer, T5<:Real, T6<:Real, T7<:Real}

   map(order->model_spectrum_1d(log_lambdas_template[order], template[order],
         view(log_lambdas_out,:,order), z_bc, pixel_range=pixel_range, blaze=blaze[order],
         output=view(output,:,order)), orders) # TODO: TEST
   return output
end

#=
ENV["PYTHON"] = "/usr/bin/python3"
using Pkg
Pkg.activate(".")
using Revise
using EchelleTools
using DataFrames, FileIO # JLD2 #
#JLD2.@load("../../Data/GJ_436.jld2")
data = load("../../Data/GJ_436.jld2")["data"]
res2d = make_blaze_and_template_2d(data)
orders = 1:12
llt = map(i->res2d[i].log_lambda_range,orders)
t = map(i->res2d[i].template,orders)
llo = log.(data.lambda[:,orders,1])
b = map(i->res2d[i].blaze,orders)
out = model_spectrum_2d(llt,t,llo,0.0,orders=orders,blaze=b )

using Polynomials, Interpolations
import EchelleTools.num_obs, EchelleTools.num_orders, EchelleTools.num_pixels_per_order
import EchelleTools.model_spectrum_2d
include("src/analyze.jl")

=#


text:
import HTTP
import JSON3

# url = "http://127.0.0.1:8000"
# q = HTTP.request("GET", url, [], [])

# HTTP.request("POST", "http://localhost:8000/jsonpayload",
#     [("Content-Type", "application/json")], """{"name":"Adrian"}""")

x = randn(3)
y = randn(3)
xq = randn(30)

input_array = Vector{Vector{Float64}}(undef, 2)
input_array[1] = x
input_array[2] = y

q =  """{"input_array":$(input_array)}"""
#q =  """{"x":$(x)}"""
#q = """{"name":"Adrian"}"""

b = JSON3.write(input_array)
c = JSON3.read(b)
d = convert(Vector{Vector{Float64}}, c)
x_rec = d[1]
y_rec = d[2]

z = HTTP.request("POST", "http://localhost:8000/jsonpayload",
    [("Content-Type", "application/json")], q)

h = String(HTTP.payload(z.request))

dic = JSON3.read(h)
Z = convert(Vector{Vector{Float64}}, dic["input_array"])

dic = JSON3.read(q)
Q = convert(Vector{Vector{Float64}}, dic["input_array"])
x_rec2 = Q[1]
y_rec2 = Q[2]

max_stars_repo_name: RoyCCWang/WebAppRKHS

max_issues_repo_name: RoyCCWang/WebAppRKHS

max_forks_repo_name: RoyCCWang/WebAppRKHS


text:
using Test
using Citations

# include("test_bib.jl")

@testset "add citation" begin
    @add_reference("kawczynski2018" => "@article{kawczynski2018characterization,....}")
    @test cite("kawczynski2018") == "@article{kawczynski2018characterization,....}"
end

# @testset "add bib" begin
#     @bibliography("bib_test.bib")
#     @test cite("angenendt") == ""
#     @test cite("bertram") == ""
# end