text:
@inline undef_memory_buffer(::Type{T}, ::StaticInt{L}) where {T,L} = MemoryBuffer{L,T}(undef)
@inline undef_memory_buffer(::Type{T}, L) where {T} = Vector{T}(undef, L)
@inline undef_memory_buffer(::Type{Bit}, ::StaticInt{L}) where {L} = MemoryBuffer{(L+7)>>>3,UInt8}(undef)
@inline undef_memory_buffer(::Type{Bit}, L) = Vector{UInt8}(undef, (L+7)>>>3)
# @inline undef_memory_buffer(::Type{Bit}, L) = BitVector(undef, L)

struct StrideArray{S,D,T,N,C,B,R,X,O,A<:Union{MemoryBuffer,AbstractArray}} <: AbstractStrideArray{S,D,T,N,C,B,R,X,O}
  ptr::PtrArray{S,D,T,N,C,B,R,X,O}
  data::A
end
struct StrideBitArray{S,D,N,C,B,R,X,O,A<:Union{MemoryBuffer,AbstractArray}} <: AbstractStrideArray{S,D,Bit,N,C,B,R,X,O}
  ptr::BitPtrArray{S,D,N,C,B,R,X,O}
  data::A
end
StrideArray(ptr::BitPtrArray{S,D,N,C,B,R,X,O}, data::A) where {S,D,N,C,B,R,X,O,A<:Union{MemoryBuffer,AbstractArray}} = StrideBitArray{S,D,N,C,B,R,X,O,A}(ptr, data)
@inline LayoutPointers.stridedpointer(A::StrideArray) = getfield(getfield(A,:ptr),:ptr)
@inline LayoutPointers.stridedpointer(A::StrideBitArray) = getfield(getfield(A,:ptr),:ptr)

const StrideVector{S,D,T,C,B,R,X,O,A} = StrideArray{S,D,T,1,C,B,R,X,O,A}
const StrideMatrix{S,D,T,C,B,R,X,O,A} = StrideArray{S,D,T,2,C,B,R,X,O,A}

@inline StrideArray(A::AbstractArray) = StrideArray(PtrArray(A), A)

@inline function StrideArray{T}(::UndefInitializer, s::Tuple{Vararg{Integer,N}}) where {N,T}
  x, L = calc_strides_len(T,s)
  b = undef_memory_buffer(T, L ÷ static_sizeof(T))
  # For now, just trust Julia's alignment heuristics are doing the right thing
  # to save us from having to over-allocate
  # ptr = LayoutPointers.align(pointer(b))
  StrideArray(reinterpret(Ptr{T}, pointer(b)), s, x, b, all_dense(Val{N}()))
end
@inline function StrideArray(ptr::Ptr{T}, s::S, x::X, b, ::Val{D}) where {S,X,T,D}
  StrideArray(PtrArray(ptr, s, x, Val{D}()), b)
end
@inline StrideArray(f, s::Vararg{Integer,N}) where {N} = StrideArray{Float64}(f, s)
@inline StrideArray(f, s::Tuple{Vararg{Integer,N}}) where {N} = StrideArray{Float64}(f, s)
@inline StrideArray(f, ::Type{T}, s::Vararg{Integer,N}) where {T,N} = StrideArray{T}(f, s)
@inline StrideArray{T}(f, s::Vararg{Integer,N}) where {T,N} = StrideArray{T}(f, s)
@inline function StrideArray(A::PtrArray{S,D,T,N}, s::Tuple{Vararg{Integer,N}}) where {S,D,T,N}
  PtrArray(stridedpointer(A), s, val_dense_dims(A))
end
@inline function StrideArray(A::AbstractArray{T,N}, s::Tuple{Vararg{Integer,N}}) where {T,N}
  StrideArray(PtrArray(stridedpointer(A), s, val_dense_dims(A)), preserve_buffer(A))
end
@inline function StrideArray{T}(f, s::Tuple{Vararg{Integer,N}}) where {T,N}
  A = StrideArray{T}(undef, s)
  @inbounds for i ∈ eachindex(A)
    A[i] = f(T)
  end
  A
end
@inline function StrideArray{T}(::typeof(zero), s::Tuple{Vararg{Integer,N}}) where {T,N}
  ptr = Ptr{T}(Libc.calloc(prod(s), sizeof(T)))
  A = unsafe_wrap(Array{T}, ptr, s; own=true)
  StrideArray(A)
end
mutable struct StaticStrideArray{S,D,T,N,C,B,R,X,O,L} <: AbstractStrideArray{S,D,T,N,C,B,R,X,O}
  data::NTuple{L,T}
  @inline StaticStrideArray{S,D,T,N,C,B,R,X,O}(::UndefInitializer) where {S,D,T,N,C,B,R,X,O} = new{S,D,T,N,C,B,R,X,O,Int(ArrayInterface.reduce_tup(*, to_static_tuple(Val(S))))}()
  @inline StaticStrideArray{S,D,T,N,C,B,R,X,O,L}(::UndefInitializer) where {S,D,T,N,C,B,R,X,O,L} = new{S,D,T,N,C,B,R,X,O,L}()
end

@generated function to_static_tuple(::Val{S}) where {S}
  t = Expr(:tuple)
  for s ∈ S.parameters
    push!(t.args, Expr(:new, s))
  end
  t
end
@inline ArrayInterface.size(::StaticStrideArray{S}) where {S} = to_static_tuple(Val(S))
@inline function ArrayInterface.strides(::StaticStrideArray{S,D,T,N,C,B,R,X}) where {S,D,T,N,C,B,R,X}
  map(Base.Fix2(>>>,intlog2(static_sizeof(T))), to_static_tuple(Val(X)))
end
@inline LayoutPointers.bytestrides(::StaticStrideArray{S,D,T,N,C,B,R,X}) where {S,D,T,N,C,B,R,X} = to_static_tuple(Val(X))
@inline ArrayInterface.offsets(::StaticStrideArray{S,D,T,N,C,B,R,X,O}) where {S,D,T,N,C,B,R,X,O} = to_static_tuple(Val(O))
@inline Base.unsafe_convert(::Type{Ptr{T}}, A::StaticStrideArray{S,D,T}) where {S,D,T} = Base.unsafe_convert(Ptr{T}, pointer_from_objref(A))
@inline Base.pointer(A::StaticStrideArray{S,D,T}) where {S,D,T} = Base.unsafe_convert(Ptr{T}, pointer_from_objref(A))
@inline StrideArray{T}(::UndefInitializer, s::Tuple{Vararg{StaticInt,N}}) where {N,T} = StaticStrideArray{T}(undef, s)
@inline StrideArray{T}(f, s::Tuple{Vararg{StaticInt,N}}) where {N,T} = StaticStrideArray{T}(f, s)
@inline StrideArray{T}(::typeof(zero), s::Tuple{Vararg{StaticInt,N}}) where {N,T} = StaticStrideArray{T}(zero, s) # Eager when static; assumed small
@inline function StaticStrideArray{T}(::UndefInitializer, s::Tuple{Vararg{StaticInt,N}}) where {N,T}
  StaticStrideArray{T}(undef, s, all_dense(Val(N)))
end
@inline function LayoutPointers.bytestrideindex(A::StaticStrideArray{S,D,T,N,C,B,R}) where {S,D,T,N,C,B,R}
  StrideIndex{N,R,C}(LayoutPointers.bytestrides(A), offsets(A))
end

@inline function StaticStrideArray{T}(::UndefInitializer, s::Tuple{Vararg{StaticInt,N}}, ::Val{D}) where {N,T,D}
  x, L = calc_strides_len(T,s)
  R = ntuple(Int, Val(N))
  O = ntuple(_ -> One(), Val(N))
  Tshifter = intlog2(static_sizeof(T))
  StaticStrideArray{typeof(s),D,T,N,1,0,R,typeof(x),typeof(O),Int(L >>> Tshifter)}(undef)
end

function dense_quote(N::Int, b::Bool)
  d = Expr(:tuple)
  for n in 1:N
    push!(d.args, b)
  end
  Expr(:call, Expr(:curly, :Val, d))
end
@generated all_dense(::Val{N}) where {N} = dense_quote(N, true)

@generated function calc_strides_len(::Type{T}, s::Tuple{Vararg{StaticInt,N}}) where {T, N}
  L = Base.allocatedinline(T) ? sizeof(T) : sizeof(Int)
  t = Expr(:tuple)
  for n ∈ 1:N
    push!(t.args, static_expr(L))
    l::Int = (s.parameters[n].parameters[1])::Int
    if (T ≢ Bit) || (n ≠ 1)
      L *= l
    else
      L *= ((l + 7) & -8)
    end
  end
  Expr(:tuple, t, static_expr(L))
end
@generated function calc_strides_len(::Type{T}, s::Tuple{Vararg{Integer,N}}) where {T, N}
  last_sx = :s_0
  st = Base.allocatedinline(T) ? sizeof(T) : sizeof(Int)
  q = Expr(:block, Expr(:meta,:inline), Expr(:(=), last_sx, static_expr(st)))
  t = Expr(:tuple)
  for n ∈ 1:N
    push!(t.args, last_sx)
    new_sx = Symbol(:s_, n)
    sz_expr = Expr(:call, getfield, :s, n, false)
    ((T === Bit) && (n == 1)) && (sz_expr = :(($sz_expr + StaticInt{7}()) & StaticInt{-8}()))
    push!(q.args, Expr(:(=), new_sx, Expr(:call, *, last_sx, sz_expr)))
    last_sx = new_sx
  end
  push!(q.args, Expr(:tuple, t, last_sx))
  q
end
@inline function StaticStrideArray{T}(f, s::Tuple{Vararg{StaticInt,N}}) where {T,N}
  A = StaticStrideArray{T}(undef, s)
  @inbounds for i ∈ eachindex(A)
    A[i] = f(T)
  end
  A
end

@inline LayoutPointers.preserve_buffer(A::MemoryBuffer) = A
@inline LayoutPointers.preserve_buffer(A::StrideArray) = preserve_buffer(getfield(A, :data))

@inline PtrArray(A::Union{StrideArray,StrideBitArray}) = getfield(A, :ptr)

@inline maybe_ptr_array(A) = A
@inline maybe_ptr_array(A::AbstractArray) = maybe_ptr_array(ArrayInterface.device(A), A)
@inline maybe_ptr_array(::ArrayInterface.CPUPointer, A::AbstractArray) = PtrArray(A)
@inline maybe_ptr_array(_, A::AbstractArray) = A

@inline ArrayInterface.size(A::Union{StrideArray,StrideBitArray}) = getfield(getfield(A, :ptr), :size)

@inline LayoutPointers.bytestrides(A::Union{StrideArray,StrideBitArray}) = bytestrides(getfield(getfield(A, :ptr), :ptr))
@inline ArrayInterface.strides(A::Union{StrideArray,StrideBitArray}) = strides(getfield(A, :ptr))
@inline ArrayInterface.offsets(A::Union{StrideArray,StrideBitArray}) = offsets(getfield(getfield(A, :ptr), :ptr))

@inline zeroindex(r::ArrayInterface.OptionallyStaticUnitRange{One}) = CloseOpen(Zero(), last(r))
@inline zeroindex(r::Base.OneTo) = CloseOpen(Zero(), last(r))
@inline zeroindex(r::AbstractUnitRange) = Zero():(last(r)-first(r))

@inline zeroindex(r::CloseOpen{Zero}) = r
@inline zeroindex(r::ArrayInterface.OptionallyStaticUnitRange{Zero}) = r
@inline zeroindex(A::PtrArray{S,D}) where {S,D} = PtrArray(zstridedpointer(A), size(A), Val{D}())
@inline zeroindex(A::Union{StrideArray,StrideBitArray}) = StrideArray(zeroindex(PtrArray(A)), preserve_buffer(A))
@inline zeroindex(A::StaticStrideArray) = StrideArray(zeroindex(PtrArray(A)), A)

@generated rank_to_sortperm_val(::Val{R}) where {R} = :(Val{$(rank_to_sortperm(R))}())
@inline function similar_layout(A::AbstractStrideArray{S,D,T,N,C,B,R}) where {S,D,T,N,C,B,R}
  permutedims(similar(permutedims(A, rank_to_sortperm_val(Val{R}()))), Val{R}())
end
@inline function similar_layout(A::AbstractArray)
  b = preserve_buffer(A)
  GC.@preserve b begin
    similar_layout(PtrArray(A))
  end
end
@inline function Base.similar(A::AbstractStrideArray{S,D,T}) where {S,D,T}
  StrideArray{T}(undef, size(A))
end
@inline function Base.similar(A::AbstractStrideArray, ::Type{T}) where {T}
  StrideArray{T}(undef, size(A))
end

@inline function Base.view(A::AbstractStrideArray, i::Vararg{Union{Integer,AbstractRange,Colon},K}) where {K}
  StrideArray(view(PtrArray(A), i...), preserve_buffer(A))
end
@inline function zview(A::AbstractStrideArray, i::Vararg{Union{Integer,AbstractRange,Colon},K}) where {K}
  StrideArray(zview(PtrArray(A), i...), preserve_buffer(A))
end
@inline function Base.permutedims(A::AbstractStrideArray, ::Val{P}) where {P}
  StrideArray(permutedims(PtrArray(A), Val{P}()), preserve_buffer(A))
end
@inline Base.adjoint(a::AbstractStrideVector) = StrideArray(adjoint(PtrArray(a)), preserve_buffer(a))

function gc_preserve_call(ex, skip=0)
  q = Expr(:block)
  call = Expr(:call, esc(ex.args[1]))
  gcp = Expr(:gc_preserve, call)
  for i ∈ 2:length(ex.args)
    arg = ex.args[i]
    if i+1 ≤ skip
      push!(call.args, arg)
      continue
    end
    A = gensym(:A); buffer = gensym(:buffer);
    if arg isa Expr && arg.head === :kw
      push!(call.args, Expr(:kw, arg.args[1], Expr(:call, :maybe_ptr_array, A)))
      arg = arg.args[2]
    else
      push!(call.args, Expr(:call, :maybe_ptr_array, A))
    end
    push!(q.args, :($A = $(esc(arg))))
    push!(q.args, Expr(:(=), buffer, Expr(:call, :preserve_buffer, A)))
    push!(gcp.args, buffer)
  end
  push!(q.args, gcp)
  q
end
"""
    @gc_preserve foo(A, B, C)

  Apply to a single, non-nested, function call. It will `GC.@preserve` all the arguments, and substitute suitable arrays with `PtrArray`s.
  This has the benefit of potentially allowing statically sized mutable arrays to be both stack allocated, and passed through a non-inlined function boundary.
  """
macro gc_preserve(ex)
  @assert ex.head === :call
  gc_preserve_call(ex)
end

@inline LayoutPointers.zero_offsets(A::AbstractArray) = StrideArray(LayoutPointers.zero_offsets(PtrArray(A)), A)


text:
mutable struct Grid
	obj::Ptr{Cvoid}
end

# destructor to be used in finalizer. Don't call explicitly
function destroy(G::Grid)
    ElError(ccall(("ElGridDestroy", libEl), Cuint,
        (Ptr{Cvoid},), G.obj))
    return nothing
end

# Returns the default Grid. The default grid is finalized when Elemental is finalized
# so we shouldn't register a `destroy` as finalizer.
function Grid()
    obj = Ref{Ptr{Cvoid}}(C_NULL)
        ElError(ccall(("ElDefaultGrid", libEl), Cuint,
            (Ref{Ptr{Cvoid}},), obj))
    return Grid(obj[])
end

const DefaultGrid = Ref{Grid}()

ext: jl


text:
import GLFW

# Create a windowed mode window and its OpenGL context
window = GLFW.CreateWindow(640, 480, "GLFW.jl")

# Make the window's context current
GLFW.MakeContextCurrent(window)

# Loop until the user closes the window
while !GLFW.WindowShouldClose(window)

	# Render here

	# Swap front and back buffers
	GLFW.SwapBuffers(window)

	# Poll for and process events
	GLFW.PollEvents()
end


text:
using ManifoldsBase: LinearAlgebra
include("../utils.jl")

@testset "Tucker" begin
    n⃗ = (4, 5, 6)
    r⃗ = (2, 3, 4)
    r⃗_small = (2, 3, 3)
    M = Tucker(n⃗, r⃗)
    M_small = Tucker(r⃗, (2, 3, 3))

    types = [Float64]
    TEST_FLOAT32 && push!(types, Float32)

    for T in types
        @testset "Type $T" begin
            U = ntuple(d -> Matrix{T}(LinearAlgebra.I, n⃗[d], r⃗[d]), 3)
            C₁ = reshape(
                T[
                    0.0 0.0 1.0 0.0
                    1.0 1/√2 0.0 0.0
                    1.0 -1/√2 0.0 0.0
                    0.0 0.0 0.0 0.0
                    0.0 0.0 0.0 1.0
                    0.0 0.0 0.0 0.0
                ],
                (2, 3, 4),
            )
            C₂ = reshape(
                T[
                    0.0 0.0 -1/√2 1/√2
                    -1/√2 -0.5 0.0 0.0
                    1.0 -1/√2 0.0 0.0
                    0.0 0.0 0.0 0.0
                    0.0 0.0 1/√2 1/√2
                    1/√2 0.5 0.0 0.0
                ],
                (2, 3, 4),
            )
            C₃ = reshape(
                T[
                    1/√2 0.5 1/√2 0.0
                    1/√2 0.5 -1/√2 0.0
                    1/√2 -0.5 0.0 0.0
                    -1/√2 0.5 0.0 0.0
                    0.0 0.0 0.0 1/√2
                    0.0 0.0 0.0 -1/√2
                ],
                (2, 3, 4),
            )

            p = TuckerPoint(C₁, U...)
            p_small = TuckerPoint(C₂, r⃗_small)
            p_ambient = embed(M, p)

            U⊥ = ntuple(d -> nullspace(U[d]') * Matrix{Bool}(I, n⃗[d] - r⃗[d], r⃗[d]), 3)
            v = TuckerTVector(reshape(collect(1.0:prod(r⃗)), r⃗), U⊥)
            U⊥ = ntuple(d -> nullspace(U[d]') * Matrix{Bool}(I, n⃗[d] - r⃗[d], r⃗[d]), 3)
            w = TuckerTVector(reshape((1.0:prod(r⃗)) .^ 2, r⃗), U⊥)

            norm_v = norm(M, p, v)
            @test norm(M, p, 2 * v - v * 2) ≤ √eps(T) * norm_v
            @test norm(M, p, (v / 2) - (2 \ v)) ≤ eps(T) * norm_v
            @test norm(M, p, (v / 2) - 0.5v) ≤ eps(T) * norm_v
            @test +v == v

            @test embed(M, p) ≈ reshape(kron(reverse(U)...) * vec(C₁), n⃗)
            @test embed(M, TuckerPoint(p_ambient, r⃗)) ≈ p_ambient
            @test inner(M, p, v, w) ≈ dot(embed(M, p, v), embed(M, p, w))
            @test inner(M, p, v, v) ≈ norm(M, p, v)^2
            w_ambient = embed(M, p, w)
            @test inner(M, p, v, w_ambient) ≈ inner(M, p, v, w) ≈ inner(M, p, w_ambient, v)

            @test p ≈ TuckerPoint(p_ambient, r⃗)
            @test v == TuckerTVector(v.Ċ, v.U̇)

            @test manifold_dimension(M) ==
                  prod(r⃗) + sum(ntuple(d -> r⃗[d] * (n⃗[d] - r⃗[d]), length(r⃗)))

            @test is_point(M_small, p_small)
            @test is_point(M, p)
            @test is_point(M, embed(M, p))
            @test !is_point(M, zeros(n⃗))
            @test !is_point(M, TuckerPoint(zeros(T, r⃗), U...))
            @test !is_point(M, TuckerPoint(C₁, map(rk -> Matrix{T}(I, rk, rk), r⃗)...))
            @test !is_point(M, p_small)
            @test !is_point(M, embed(M_small, p_small))
            @test !is_point(M_small, p)
            q = allocate(p)
            copyto!(q, p)
            q.hosvd.core .= 2 * p.hosvd.core
            @test !is_point(M, q)

            u = allocate(v)
            u .= v
            @test u == v
            u = allocate(v)
            copyto!(u, v)
            @test u == v

            # broadcasting
            @test axes(v) === ()
            u = copy(v)
            # test that the copy is equal to the original, but represented by
            # a new array
            @test u.U̇ !== v.U̇
            @test u.Ċ == v.Ċ
            x = u .+ v .* 2
            @test x isa TuckerTVector
            @test x == v + u * 2
            x .= 2 .* u .+ v
            @test x == 2 * u + v

            @test is_vector(M, p, v)
            @test !is_vector(M, p_small, v)
            @test !is_vector(
                M,
                p,
                TuckerTVector(ones(T, r⃗), map(u⊥ -> u⊥[:, 1:(end - 1)], U⊥)),
            )

            ℬ = get_basis(M, p, DefaultOrthonormalBasis())
            J = convert(Matrix, ℬ)
            @test J'J ≈ LinearAlgebra.I
            vectors = get_vectors(M, p, DefaultOrthonormalBasis())
            vectors_mtx = hcat(map(ξ -> vec(embed(M, p, ξ)), vectors)...)
            @test J ≈ vectors_mtx

            shiftprint(x) = replace(sprint(show, "text/plain", x), "\n" => "\n ")

            @test sprint(show, "text/plain", M) == "Tucker((4, 5, 6), (2, 3, 4), ℝ)"
            @test sprint(show, "text/plain", p) == """
            $(summary(p))
            U factor 1:
             $(shiftprint(p.hosvd.U[1]))\n
            U factor 2:
             $(shiftprint(p.hosvd.U[2]))\n
            U factor 3:
             $(shiftprint(p.hosvd.U[3]))\n
            Core:
             $(shiftprint(p.hosvd.core))"""

            @test sprint(show, "text/plain", v) == """
            $(summary(v))
            U̇ factor 1:
             $(shiftprint(v.U̇[1]))\n
            U̇ factor 2:
             $(shiftprint(v.U̇[2]))\n
            U̇ factor 3:
             $(shiftprint(v.U̇[3]))\n
            Ċ factor:
             $(shiftprint(v.Ċ))"""

            @test sprint(show, "text/plain", ℬ) == """
            $(summary(ℬ)) ≅ $(shiftprint(J))
            """

            pts = [
                p,
                TuckerPoint(C₂, U[1][[4, 2, 1, 3], :], U[2:end]...),
                TuckerPoint(C₃, U[1][[4, 3, 1, 2], :], U[2:end]...),
            ]
            test_manifold(
                M,
                pts;
                is_mutating=false, # avoid allocations of the wrong type
                basis_types_to_from=(DefaultOrthonormalBasis(),),
                basis_types_vecs=(DefaultOrthonormalBasis(),),
                test_exp_log=false,
                default_inverse_retraction_method=ProjectionInverseRetraction(),
                test_injectivity_radius=false,
                default_retraction_method=PolarRetraction(),
                test_is_tangent=false,
                test_project_tangent=false,
                test_default_vector_transport=false,
                test_vector_spaces=false,
                test_vee_hat=false,
                test_tangent_vector_broadcasting=true,
                test_representation_size=false,
                projection_atol_multiplier=15,
                retraction_methods=[PolarRetraction()],
                inverse_retraction_methods=[ProjectionInverseRetraction()],
                mid_point12=nothing,
            )
        end
    end
end


text:
# This file is part of TaylorSeries.jl, MIT licensed

# HomogeneousPolynomial and TaylorN types, and functions defined on them

## Minimum order of a HomogeneousPolynomial compatible with the vector's length
function orderH{T}(coeffs::Array{T,1})
    ord = 0
    ll = length(coeffs)
    for i = 1:get_order()+1
        @inbounds num_coeffs = size_table[i]
        ll <= num_coeffs && break
        ord += 1
    end
    return ord
end

## HomogeneousPolynomial (homogeneous polynomial) constructors ##
doc"""
    HomogeneousPolynomial{T<:Number} <: Number

DataType for homogenous polynomials in many (>1) independent variables.

**Fields:**

- `coeffs  :: Array{T,1}` Expansion coefficients of the homogeneous
polynomial; the $i$-th component is related to a monomial, where the degrees
of the independent variables are specified by `coeff_table[order+1][i]`.
- `order   :: Int` order (degree) of the homogenous polynomial.
"""
immutable HomogeneousPolynomial{T<:Number} <: Number
    coeffs  :: Array{T,1}
    order   :: Int

    function HomogeneousPolynomial( coeffs::Array{T,1}, order::Int )
        maxOrder = get_order()
        @assert order <= maxOrder
        lencoef = length( coeffs )
        @inbounds num_coeffs = size_table[order+1]
        @assert lencoef <= num_coeffs
        num_coeffs == lencoef && return new(coeffs, order)
        resize!(coeffs, num_coeffs)
        @simd for i = lencoef+1:num_coeffs
            @inbounds coeffs[i] = zero(T)
        end
        new(coeffs, order)
    end
end

HomogeneousPolynomial{T<:Number}(x::HomogeneousPolynomial{T}, order::Int) =
    HomogeneousPolynomial{T}(x.coeffs, order)
HomogeneousPolynomial{T<:Number}(x::HomogeneousPolynomial{T}) = x
HomogeneousPolynomial{T<:Number}(coeffs::Array{T,1}, order::Int) =
    HomogeneousPolynomial{T}(coeffs, order)
HomogeneousPolynomial{T<:Number}(coeffs::Array{T,1}) =
    HomogeneousPolynomial{T}(coeffs, orderH(coeffs))
HomogeneousPolynomial{T<:Number}(x::T, order::Int) =
    HomogeneousPolynomial{T}([x], order)
HomogeneousPolynomial{T<:Number}(x::T) = HomogeneousPolynomial{T}([x], 0)

eltype{T<:Number}(::HomogeneousPolynomial{T}) = T
length(a::HomogeneousPolynomial) = length( a.coeffs )
get_order(a::HomogeneousPolynomial) = a.order

## zero and one ##
function zero{T<:Number}(a::HomogeneousPolynomial{T})
    a.order == 0 && return HomogeneousPolynomial(zero(T), 0)
    @inbounds num_coeffs = size_table[a.order+1]
    return HomogeneousPolynomial{T}( zeros(T,num_coeffs), a.order)
end

function zeros{T<:Number}(::HomogeneousPolynomial{T}, order::Int)
    order == 0 && return [HomogeneousPolynomial(zero(T),0)]
    v = Array(HomogeneousPolynomial{T}, order+1)
    @simd for ord in eachindex(v)
        @inbounds num_coeffs = size_table[ord]
        @inbounds v[ord] = HomogeneousPolynomial(zeros(T,num_coeffs),ord-1)
    end
    return v
end

zeros{T<:Number}(::Type{HomogeneousPolynomial{T}}, order::Int) =
    zeros( HomogeneousPolynomial(zero(T), 0), order)

function one{T<:Number}(a::HomogeneousPolynomial{T})
    a.order == 0 && return HomogeneousPolynomial(one(T), 0)
    @inbounds num_coeffs = size_table[a.order+1]
    return HomogeneousPolynomial{T}( ones(T,num_coeffs), a.order)
end

function ones{T<:Number}(::HomogeneousPolynomial{T}, order::Int)
    order == 0 && return [HomogeneousPolynomial(one(T),0)]
    v = Array(HomogeneousPolynomial{T}, order+1)
    @simd for ord in eachindex(v)
        @inbounds num_coeffs = size_table[ord]
        @inbounds v[ord] = HomogeneousPolynomial(ones(T,num_coeffs),ord-1)
    end
    return v
end

ones{T<:Number}(::Type{HomogeneousPolynomial{T}}, order::Int) =
    ones( HomogeneousPolynomial(one(T), 0), order)

## Conversion and promotion rules ##
convert{T<:Number}(::Type{HomogeneousPolynomial{T}}, a::HomogeneousPolynomial) =
    HomogeneousPolynomial{T}(convert(Array{T,1}, a.coeffs), a.order)
function convert{T<:Integer, S<:AbstractFloat}(
    ::Type{HomogeneousPolynomial{Rational{T}}}, a::HomogeneousPolynomial{S})
    v = Array(Rational{T}, length(a.coeffs))
    for i in eachindex(v)
        # v[i] = convert(Rational{T}, a.coeffs[i])
        v[i] = rationalize(a.coeffs[i], tol=eps(one(S)))
    end
    HomogeneousPolynomial(v, a.order)
end
convert{T<:Number, S<:Number}(::Type{HomogeneousPolynomial{T}}, b::Array{S,1}) =
    HomogeneousPolynomial{T}(convert(Array{T,1}, b), orderH(b))
convert{T<:Number, S<:Number}(::Type{HomogeneousPolynomial{T}}, b::S) =
    HomogeneousPolynomial{T}([convert(T,b)], 0)
convert{T<:Number}(::Type{HomogeneousPolynomial{T}}, b::Array{T,1}) =
    HomogeneousPolynomial{T}(b, orderH(b))
convert{T<:Number}(::Type{HomogeneousPolynomial{T}}, b::T) =
    HomogeneousPolynomial{T}([b], 0)

promote_rule{T<:Number, S<:Number}(::Type{HomogeneousPolynomial{T}},
    ::Type{HomogeneousPolynomial{S}}) = HomogeneousPolynomial{promote_type(T,S)}
promote_rule{T<:Number, S<:Number}(::Type{HomogeneousPolynomial{T}},
    ::Type{Array{S,1}}) = HomogeneousPolynomial{promote_type(T, S)}
promote_rule{T<:Number,S<:Union{Real,Complex}}(::Type{HomogeneousPolynomial{T}},
    ::Type{S}) = HomogeneousPolynomial{promote_type(T,S)}

## Maximum order of a HomogeneousPolynomial vector; used by TaylorN constructor
function maxorderH{T<:Number}(v::Array{HomogeneousPolynomial{T},1})
    m = 0
    @inbounds for i in eachindex(v)
        m = max(m, v[i].order)
    end
    return m
end

doc"""
    TaylorN{T<:Number} <: Number

DataType for polynomial expansions in many (>1) independent variables.

**Fields:**

- `coeffs  :: Array{HomogeneousPolynomial{T},1}` Vector containing the
`HomogeneousPolynomial` entries. The $i$-th component corresponds to the
homogeneous polynomial of degree $i-1$.
- `order   :: Int`  maximum order of the polynomial expansion.
"""
immutable TaylorN{T<:Number} <: Number
    coeffs  :: Array{HomogeneousPolynomial{T},1}
    order   :: Int

    function TaylorN( v::Array{HomogeneousPolynomial{T},1}, order::Int )
        m = maxorderH(v)
        order = max( m, order )
        coeffs = zeros(HomogeneousPolynomial{T}, order)
        @simd for i in eachindex(v)
            @inbounds ord = v[i].order
            @inbounds coeffs[ord+1] += v[i]
        end
        new(coeffs, order)
    end
end

TaylorN{T<:Number}(x::TaylorN{T}, order::Int) = TaylorN{T}(x.coeffs, order)
TaylorN{T<:Number}(x::TaylorN{T}) = x
TaylorN{T<:Number}(x::Array{HomogeneousPolynomial{T},1}, order::Int) =
    TaylorN{T}(x, order)
TaylorN{T<:Number}(x::Array{HomogeneousPolynomial{T},1}) = TaylorN{T}(x,0)
TaylorN{T<:Number}(x::HomogeneousPolynomial{T},order::Int) =
    TaylorN{T}([x], order)
TaylorN{T<:Number}(x::HomogeneousPolynomial{T}) = TaylorN{T}([x], x.order)
TaylorN{T<:Number}(x::T, order::Int) =
    TaylorN{T}([HomogeneousPolynomial(x)],order)
# TaylorN{T<:Number}(x::T) = TaylorN{T}([HomogeneousPolynomial(x)], 0)

## Shortcut to define TaylorN independent variables
"""
    TaylorN(T, nv; [order=get_order()])
    TaylorN(nv; [order=get_order()])

Shortcut to define the `nv`-th independent `TaylorN{T}` variable as a
polynomial. The order is defined through the keyword parameter `order`,
whose default corresponds to `get_order()`. If `T::Type` is ommitted,
`Float64` is assumend.
"""
function TaylorN{T<:Number}(::Type{T}, nv::Int; order::Int=get_order())
    @assert 0 < nv <= get_numvars()
    v = zeros(T, get_numvars())
    @inbounds v[nv] = one(T)
    return TaylorN( HomogeneousPolynomial(v,1), order )
end
TaylorN(nv::Int; order::Int=get_order()) = TaylorN(Float64, nv; order=order)

## get_coeff
"""
    get_coeff(a, v)

Return the coefficient of `a::HomogeneousPolynomial`, specified by
`v::Array{Int,1}` which has the indices of the specific monomial.
"""
function get_coeff(a::HomogeneousPolynomial, v::Array{Int,1})
    @assert length(v) == get_numvars()
    kdic = in_base(get_order(),v)
    @inbounds n = pos_table[a.order+1][kdic]
    a.coeffs[n]
end

"""
    get_coeff(a, v)

Return the coefficient of `a::TaylorN`, specified by
`v::Array{Int,1}` which has the indices of the specific monomial.
"""
function get_coeff(a::TaylorN, v::Array{Int,1})
    order = sum(v)
    get_coeff(a.coeffs[order+1], v)
end

## Type, length ##
eltype{T<:Number}(::TaylorN{T}) = T
length(a::TaylorN) = length( a.coeffs )
get_order(x::TaylorN) = x.order

## zero and one ##
zero{T<:Number}(a::TaylorN{T}) = TaylorN(zero(T), a.order)
one{T<:Number}(a::TaylorN{T}) = TaylorN(one(T), a.order)

## Conversion and promotion rules ##
convert{T<:Number}(::Type{TaylorN{T}}, a::TaylorN) =
    TaylorN{T}( convert(Array{HomogeneousPolynomial{T},1}, a.coeffs), a.order)
convert{T<:Number, S<:Number}(::Type{TaylorN{T}}, b::HomogeneousPolynomial{S}) =
    TaylorN{T}( [convert(HomogeneousPolynomial{T}, b)], b.order)
convert{T<:Number, S<:Number}(::Type{TaylorN{T}},
    b::Array{HomogeneousPolynomial{S},1}) =
    TaylorN{T}( convert(Array{HomogeneousPolynomial{T},1}, b), length(b)-1)
convert{T<:Number, S<:Number}(::Type{TaylorN{T}}, b::S) =
    TaylorN( [HomogeneousPolynomial(convert(T, b))], 0)
convert{T<:Number}(::Type{TaylorN{T}}, b::HomogeneousPolynomial{T}) =
    TaylorN{T}( [b], b.order)
convert{T<:Number}(::Type{TaylorN{T}}, b::Array{HomogeneousPolynomial{T},1}) =
    TaylorN{T}( b, length(b)-1)
convert{T<:Number}(::Type{TaylorN{T}}, b::T) =
    TaylorN( [HomogeneousPolynomial(b)], 0)

promote_rule{T<:Number, S<:Number}(::Type{TaylorN{T}}, ::Type{TaylorN{S}}) =
    TaylorN{promote_type(T, S)}
promote_rule{T<:Number, S<:Number}(::Type{TaylorN{T}},
    ::Type{HomogeneousPolynomial{S}}) = TaylorN{promote_type(T, S)}
promote_rule{T<:Number, S<:Number}(::Type{TaylorN{T}},
    ::Type{Array{HomogeneousPolynomial{S},1}}) = TaylorN{promote_type(T, S)}
promote_rule{T<:Number, S<:Number}(::Type{TaylorN{T}}, ::Type{S}) =
    TaylorN{promote_type(T, S)}

## Auxiliary function ##
fixorder{T<:Number}(a::TaylorN{T}, order::Int64) = TaylorN{T}(a.coeffs, order)

function fixorder(a::TaylorN, b::TaylorN)
    a.order == b.order && return a, b
    a.order < b.order && return TaylorN(a.coeffs, b.order), b
    return a, TaylorN(b.coeffs, a.order)
end

function fixshape(a::HomogeneousPolynomial, b::HomogeneousPolynomial)
    @assert a.order == b.order
    return promote(a, b)
end

function fixshape(a::TaylorN, b::TaylorN)
    if eltype(a) != eltype(b)
        a, b = promote(a, b)
    end
    return fixorder(a, b)
end

## real, imag, conj and ctranspose ##
for TT in (:HomogeneousPolynomial, :TaylorN), f in (:real, :imag, :conj)
    @eval ($f){T<:Number}(a::($TT){T}) = ($TT)(($f)(a.coeffs), a.order)
end

ctranspose{T<:Number}(a::HomogeneousPolynomial{T}) = conj(a)
ctranspose{T<:Number}(a::TaylorN{T}) = conj(a)

## Equality ##
==(a::HomogeneousPolynomial, b::HomogeneousPolynomial) = a.coeffs == b.coeffs
function ==(a::TaylorN, b::TaylorN)
    a, b = fixshape(a, b)
    test = true
    for i in eachindex(a.coeffs)
        @inbounds test = a.coeffs[i] == b.coeffs[i]
        test || return test
    end
    return test
end

function iszero{T}(a::HomogeneousPolynomial{T})
    test = true
    for i in eachindex(a.coeffs)
        @inbounds test = a.coeffs[i] == zero(T)
        test || return test
    end
    return test
end

## Addition and substraction ##
for T in (:HomogeneousPolynomial, :TaylorN), f in (:+, :-)
    @eval begin
        function ($f)(a::($T), b::($T))
            a, b = fixshape(a, b)
            v = Array(eltype(a.coeffs), length(a.coeffs))
            for i in eachindex(v)
                @inbounds v[i] = ($f)(a.coeffs[i], b.coeffs[i])
            end
            return ($T)(v, a.order)
        end
        function ($f)(a::($T))
            v = Array(eltype(a.coeffs), length(a.coeffs))
            for i in eachindex(v)
                @inbounds v[i] = ($f)(a.coeffs[i])
            end
            return ($T)(v, a.order)
        end
    end
end
for f in (:+, :-)
    @eval begin
        function ($f)(a::TaylorN, b::Union{Real,Complex})
            @inbounds aux = ($f)(a.coeffs[1], b)
            S = eltype(aux)
            coeffs = Array(HomogeneousPolynomial{S},length(a.coeffs))
            @simd for i in eachindex(coeffs)
                @inbounds coeffs[i] = a.coeffs[i]
            end
            @inbounds coeffs[1] = aux
            return TaylorN{S}(coeffs, a.order)
        end
        function ($f)(b::Union{Real,Complex}, a::TaylorN)
            @inbounds aux = ($f)(b, a.coeffs[1])
            S = eltype(aux)
            coeffs = Array(HomogeneousPolynomial{S},length(a.coeffs))
            @simd for i in eachindex(coeffs)
                @inbounds coeffs[i] = ($f)(a.coeffs[i])
            end
            @inbounds coeffs[1] = aux
            return TaylorN{S}(coeffs, a.order)
        end
    end
end

## Multiplication ##
function *(a::HomogeneousPolynomial, b::HomogeneousPolynomial)
    T = promote_type( eltype(a), eltype(b) )

    order = a.order + b.order
    order > get_order() && return HomogeneousPolynomial(zero(T), get_order())

    if eltype(a) != eltype(b)
        a, b = promote(a, b)
    end

    res = HomogeneousPolynomial(zero(T), order)
    mul!(res, a, b)
    return res
end

function *(a::TaylorN, b::TaylorN)
    a, b = fixshape(a, b)
    T = eltype(a)
    coeffs = zeros(HomogeneousPolynomial{T}, a.order)

    for ord in eachindex(coeffs)
        for i = 0:ord-1
            @inbounds mul!(coeffs[ord], a.coeffs[i+1], b.coeffs[ord-i])
        end
    end

    return TaylorN{T}(coeffs, a.order)
end

*(a::Bool, b::HomogeneousPolynomial) = *(promote(a,b)...)
*(a::HomogeneousPolynomial, b::Bool) = b * a
function *{T<:Union{Real,Complex}}(a::HomogeneousPolynomial, b::T)
    @inbounds aux = a.coeffs[1] * b
    S = typeof(aux)
    coeffs = Array(S,length(a.coeffs))
    @simd for i in eachindex(coeffs)
        @inbounds coeffs[i] = a.coeffs[i] * b
    end
    return HomogeneousPolynomial{S}(coeffs, a.order)
end
*{T<:Union{Real,Complex}}(b::T, a::HomogeneousPolynomial) = a * b

*(a::Bool, b::TaylorN) = *(promote(a,b)...)
*(a::TaylorN, b::Bool) = b * a
function *{T<:Union{Real,Complex}}(a::TaylorN, b::T)
    @inbounds aux = a.coeffs[1] * b
    S = eltype(aux)
    coeffs = Array(HomogeneousPolynomial{S},length(a.coeffs))
    @simd for i in eachindex(coeffs)
        @inbounds coeffs[i] = a.coeffs[i] * b
    end
    return TaylorN{S}(coeffs, a.order)
end
*{T<:Union{Real,Complex}}(b::T, a::TaylorN) = a * b

"""
```
mul!(c, a, b)
```

Return `c = a*b` with no allocation; all parameters are `HomogeneousPolynomial`.
"""
function mul!(c::HomogeneousPolynomial, a::HomogeneousPolynomial, b::HomogeneousPolynomial)
    (iszero(b) || iszero(a)) && return nothing

    T = eltype(c)
    @inbounds num_coeffs_a = size_table[a.order+1]
    @inbounds num_coeffs_b = size_table[b.order+1]
    @inbounds num_coeffs  = size_table[c.order+1]

    coeffs = c.coeffs
    @inbounds posTb = pos_table[c.order+1]

    @inbounds for na = 1:num_coeffs_a
        ca = a.coeffs[na]
        ca == zero(T) && continue
        inda = index_table[a.order+1][na]

        @inbounds for nb = 1:num_coeffs_b
            cb = b.coeffs[nb]
            cb == zero(T) && continue
            indb = index_table[b.order+1][nb]

            pos = posTb[inda + indb]
            coeffs[pos] += ca * cb
        end
    end

    return nothing
end

"""
    mul!(c, a)

Return `c = a*a` with no allocation; all parameters are `HomogeneousPolynomial`.
"""
function mul!(c::HomogeneousPolynomial, a::HomogeneousPolynomial)
    iszero(a) && return nothing

    T = eltype(c)
    @inbounds num_coeffs_a = size_table[a.order+1]
    @inbounds num_coeffs  = size_table[c.order+1]

    two = convert(T,2)
    coeffs = c.coeffs
    @inbounds posTb = pos_table[c.order+1]

    @inbounds for na = 1:num_coeffs_a
        ca = a.coeffs[na]
        ca == zero(T) && continue
        inda = index_table[a.order+1][na]
        @inbounds pos = posTb[2*inda]
        @inbounds coeffs[pos] += ca * ca
        @inbounds for nb = na+1:num_coeffs_a
            cb = a.coeffs[nb]
            cb == zero(T) && continue
            indb = index_table[a.order+1][nb]
            pos = posTb[inda+indb]
            coeffs[pos] += two * ca * cb
        end
    end

    return nothing
end

## Division ##
/{T<:Union{Real, Complex}}(a::HomogeneousPolynomial, x::T) = a*inv(x)
/{T<:Union{Real, Complex}}(a::TaylorN, x::T) = a*inv(x)
function /(a::TaylorN, b::TaylorN)
    @inbounds b0 = b.coeffs[1].coeffs[1]
    @assert b0 != zero(b0)
    a, b = fixshape(a, b)
    # order and coefficient of first factorized term
    # orddivfact, cdivfact = divfactorization(a, b)
    b0 = inv(b0)
    @inbounds cdivfact = a.coeffs[1] * b0
    T = eltype(cdivfact)
    coeffs = zeros(HomogeneousPolynomial{T}, a.order)
    @inbounds coeffs[1] = cdivfact

    for ord in eachindex(coeffs)
        ord == a.order+1 && continue
        @inbounds for i = 0:ord-1
            mul!(coeffs[ord+1], coeffs[i+1], b.coeffs[ord-i+1])
        end
        @inbounds coeffs[ord+1] = (a.coeffs[ord+1] - coeffs[ord+1]) * b0
    end

    return TaylorN{T}(coeffs, a.order)
end

## TODO: Implement factorization (divfactorization) for TaylorN polynomials

## Division functions: rem and mod
for op in (:mod, :rem)
    @eval begin
        @inbounds function ($op){T<:Real,S<:Real}(a::TaylorN{T}, x::S)
            coeffs = copy(a.coeffs)
            y = ($op)(a.coeffs[1].coeffs[1], x)
            coeffs[1] = HomogeneousPolynomial(y)
            return TaylorN{T}( coeffs, a.order )
        end
    end
end

function mod2pi{T<:Real}(a::TaylorN{T})
    coeffs = copy(a.coeffs)
    @inbounds y = mod2pi(a.coeffs[1].coeffs[1])
    @inbounds coeffs[1] = HomogeneousPolynomial(y)
    return TaylorN{T}( coeffs, a.order )
end

## abs function ##
"""
    abs(a::TaylorN)

Absolute value of a `TaylorN` polynomial, using the 0-th order coefficient.

Return `a` or `-a`, depending on the 0-th order coefficient of `a`.
If it is zero, it throws an `ArgumentError`.
"""
function abs{T<:Real}(a::TaylorN{T})
    if a.coeffs[1].coeffs[1] > zero(T)
        return a
    elseif a.coeffs[1].coeffs[1] < zero(T)
        return -a
    else
        throw(ArgumentError(
        """The 0th order TaylorN coefficient must be non-zero
        (`abs(x)` is not differentiable at zero)."""))
    end
end

## Int power ##
function ^(a::HomogeneousPolynomial, n::Integer)
    n == 0 && return one(a)
    n == 1 && return a
    n == 2 && return square(a)
    n < 0 && throw(DomainError())
    return power_by_squaring(a, n)
end

function ^{T<:Number}(a::TaylorN{T}, n::Integer)
    n == 0 && return one(a)
    n == 1 && return a
    n == 2 && return square(a)
    n < 0 && return inv( a^(-n) )
    return power_by_squaring(a, n)
end

function ^{T<:Integer}(a::TaylorN{T}, n::Integer)
    n == 0 && return one(a)
    n == 1 && return a
    n == 2 && return square(a)
    n < 0 && throw(DomainError())
    return power_by_squaring(a, n)
end

## power_by_squaring; slightly modified from base/intfuncs.jl
## Licensed under MIT "Expat"
for T in (:HomogeneousPolynomial, :TaylorN)
    @eval begin
        function power_by_squaring(x::($T), p::Integer)
            p == 1 && return x# copy(x)
            p == 0 && return one(x)
            p == 2 && return square(x)
            t = trailing_zeros(p) + 1
            p >>= t

            while (t -= 1) > 0
                x = square(x)
            end

            y = x
            while p > 0
                t = trailing_zeros(p) + 1
                p >>= t
                while (t -= 1) >= 0
                    x = square(x)
                end
                y *= x
            end

            return y
        end
    end
end

## Rational power ##
^(a::TaylorN, x::Rational) = a^(x.num/x.den)

^(a::TaylorN, b::TaylorN) = exp( b*log(a) )

## Real power ##
function ^{S<:Real}(a::TaylorN, x::S)
    x == zero(x) && return TaylorN( one(eltype(a)), 0 )
    x == one(x)/2 && return sqrt(a)
    isinteger(x) && return a^round(Int,x)
    @inbounds a0 = a.coeffs[1].coeffs[1]
    @assert a0 != zero(a0)
    aux = ( a0 )^x
    T = typeof(aux)
    coeffs = zeros(HomogeneousPolynomial{T}, a.order)
    @inbounds coeffs[1] = HomogeneousPolynomial( aux )

    for ord in eachindex(coeffs)
        ord == a.order+1 && continue
        @inbounds for i = 0:ord-1
            tt = x*(ord-i)-i
            cpol = coeffs[i+1]
            apol = a.coeffs[ord-i+1]
            (iszero(cpol) || iszero(apol)) && continue
            coeffs[ord+1] += tt * cpol * apol
        end
        coeffs[ord+1] = coeffs[ord+1] / (ord*a0)
    end

    return TaylorN{T}(coeffs, a.order)
end
^{T<:Complex}(a::TaylorN, x::T) = exp( x*log(a) )

## Square ##
function square(a::HomogeneousPolynomial)
    T = eltype(a)

    order = 2*a.order
    order > get_order() && return HomogeneousPolynomial(zero(T), get_order())

    res = HomogeneousPolynomial(zero(T), order)
    mul!(res, a)
    return res
end

function square(a::TaylorN)
    T = eltype(a)
    coeffs = zeros(HomogeneousPolynomial{T}, a.order)
    @inbounds mul!(coeffs[1], a.coeffs[1])

    two = convert(T,2)
    for ord in eachindex(coeffs)
        ord == a.order+1 && continue
        kodd = ord%2
        kord = div(ord-2+kodd, 2)
        for i = 0 : kord
            @inbounds mul!(coeffs[ord+1], a.coeffs[i+1], a.coeffs[ord-i+1])
        end
        @inbounds coeffs[ord+1] = two * coeffs[ord+1]
        kodd == 1 && continue
        kodd = div(ord,2)
        @inbounds mul!(coeffs[ord+1], a.coeffs[kodd+1] )
    end

    return TaylorN{T}(coeffs, a.order)
end

## sqrt ##
function sqrt(a::TaylorN)
    @inbounds p0 = sqrt( a.coeffs[1].coeffs[1] )
    if p0 == zero(p0)
        throw(ArgumentError(
        """The 0-th order TaylorN coefficient must be non-zero
        in order to expand `sqrt` around 0."""))
    end

    T = typeof(p0)
    coeffs = zeros(HomogeneousPolynomial{T}, a.order)
    @inbounds coeffs[1] = HomogeneousPolynomial( p0 )
    two = convert(T,2)

    for ord in eachindex(coeffs)
        ord == a.order+1 && continue
        kodd = ord%2
        kord = div(ord-2+kodd, 2)
        @inbounds for i = 1:kord
            coeffs[ord+1] += coeffs[i+1] * coeffs[ord-i+1]
        end
        @inbounds coeffs[ord+1] = a.coeffs[ord+1] - two * coeffs[ord+1]
        if iseven(ord)
            @inbounds coeffs[ord+1]=coeffs[ord+1]-square( coeffs[div(ord,2)+1])
        end
        @inbounds coeffs[ord+1] = coeffs[ord+1] / (two * p0)
    end

    return TaylorN{T}(coeffs, a.order)
end

## exp ##
function exp(a::TaylorN)
    order = a.order
    @inbounds aux = exp( a.coeffs[1].coeffs[1] )
    T = typeof(aux)
    coeffs = zeros(HomogeneousPolynomial{T}, order)
    @inbounds coeffs[1] = HomogeneousPolynomial(aux, 0)

    @inbounds for ord in eachindex(coeffs)
        ord == order+1 && continue
        @inbounds for j = 0:ord-1
            coeffs[ord+1] += (ord-j) * a.coeffs[ord-j+1] * coeffs[j+1]
        end
        @inbounds coeffs[ord+1] = coeffs[ord+1] / ord
    end

    return TaylorN{T}(coeffs, order)
end

## log ##
function log(a::TaylorN)
    order = a.order
    @inbounds a0 = a.coeffs[1].coeffs[1]
    if a0 == zero(a0)
        throw(ArgumentError(
        """The 0-th order `TaylorN` coefficient must be non-zero
        in order to expand `log` around 0."""))
    end
    l0 = log( a0 )
    T = typeof(l0)
    coeffs = zeros(HomogeneousPolynomial{T}, order)
    @inbounds coeffs[1] = HomogeneousPolynomial(l0)

    @inbounds for ord in eachindex(coeffs)
        ord == order+1 && continue
        @inbounds for j = 1:ord-1
            coeffs[ord+1] += j * a.coeffs[ord-j+1] * coeffs[j+1]
        end
        coeffs[ord+1] = (a.coeffs[ord+1] - coeffs[ord+1] / ord ) / a0
    end

    return TaylorN{T}(coeffs, order)
end

## sin and cos ##
sin(a::TaylorN) = imag( exp(im*a) )
cos(a::TaylorN) = real( exp(im*a) )
# sin(a::TaylorN) = sincos(a)[1]
# cos(a::TaylorN) = sincos(a)[2]
# function sincos(a::TaylorN)
#     order = a.order
#     @inbounds a0 = a.coeffs[1].coeffs[1]
#     s0 = sin( a0 )
#     c0 = cos( a0 )
#     T = typeof(s0)
#     coeffsSin = zeros(HomogeneousPolynomial{T}, order)
#     coeffsCos = zeros(HomogeneousPolynomial{T}, order)
#     @inbounds coeffsSin[1] = HomogeneousPolynomial(s0)
#     @inbounds coeffsCos[1] = HomogeneousPolynomial(c0)
#
#     @inbounds for ord in eachindex(coeffsSin)
#         ord == order+1 && continue
#         @inbounds for j = 0:ord-1
#             coeffsSin[ord+1] += (ord-j) * a.coeffs[ord-j+1] * coeffsCos[j+1]
#             coeffsCos[ord+1] += (ord-j) * a.coeffs[ord-j+1] * coeffsSin[j+1]
#         end
#         @inbounds coeffsSin[ord+1] =  coeffsSin[ord+1] / ord
#         @inbounds coeffsCos[ord+1] = -coeffsCos[ord+1] / ord
#     end
#     return TaylorN{T}(coeffsSin, order), TaylorN{T}(coeffsCos, order)
# end

## tan ##
function tan(a::TaylorN)
    order = a.order
    @inbounds a0 = a.coeffs[1].coeffs[1]
    t0 = tan(a0)
    T = typeof(t0)
    coeffsTan = zeros(HomogeneousPolynomial{T}, order)
    @inbounds coeffsTan[1] = HomogeneousPolynomial(t0)

    @inbounds for ord in eachindex(coeffsTan)
        ord == order+1 && continue
        v = coeffsTan[1:ord]
        tAux = (TaylorN(v, ord))^2
        @inbounds for j = 0:ord-1
            coeffsTan[ord+1] += (ord-j) * a.coeffs[ord-j+1] * tAux.coeffs[j+1]
        end
        coeffsTan[ord+1] = a.coeffs[ord+1] + coeffsTan[ord+1] / ord
    end

    return TaylorN{T}(coeffsTan, order)
end

## Differentiation ##
"""
    derivative(a, r)

Partial differentiation of `a::HomogeneousPolynomial` series with respect
to the `r`-th variable.
"""
function derivative(a::HomogeneousPolynomial, r::Int)
    @assert 1 <= r <= get_numvars()
    T = eltype(a)
    a.order == 0 && return HomogeneousPolynomial(zero(T))
    @inbounds num_coeffs = size_table[a.order]
    coeffs = zeros(T,num_coeffs)
    @inbounds posTb = pos_table[a.order]
    @inbounds num_coeffs = size_table[a.order+1]

    @inbounds for i = 1:num_coeffs
        iind = coeff_table[a.order+1][i]
        n = iind[r]
        n == 0 && continue
        iind[r] -= 1
        kdic = in_base(get_order(),iind)
        pos = posTb[kdic]
        coeffs[pos] = n * a.coeffs[i]
        iind[r] += 1
    end

    return HomogeneousPolynomial{T}(coeffs, a.order-1)
end

"""
    derivative(a, [r=1])

Partial differentiation of `a::TaylorN` series with respect
to the `r`-th variable.
"""
function derivative(a::TaylorN, r=1::Int)
    T = eltype(a)
    coeffs = Array(HomogeneousPolynomial{T},a.order)

    @inbounds for ord in eachindex(coeffs)
        ord == a.order+1 && continue
        coeffs[ord] = derivative( a.coeffs[ord+1], r)
    end
    return TaylorN{T}( coeffs, a.order )
end

## Gradient, jacobian and hessian
"""
```
    gradient(f)
    ∇(f)
```

Compute the gradient of the polynomial `f::TaylorN`.
"""
function gradient(f::TaylorN)
    T = eltype(f)
    numVars = get_numvars()
    grad = Array(TaylorN{T}, numVars)
    @inbounds for nv = 1:numVars
        grad[nv] = derivative(f, nv)
    end
    return grad
end
const ∇ = gradient

"""
```
    jacobian(vf)
    jacobian(vf, [vals])
```

Compute the jacobian matrix of `vf`, a vector of `TaylorN` polynomials,
evaluated at the vector `vals`. If `vals` is ommited, it is evaluated at zero.
"""
function jacobian{T<:Number}(vf::Array{TaylorN{T},1})
    numVars = get_numvars()
    @assert length(vf) == numVars
    jac = Array(T, (numVars,numVars))

    @inbounds for comp = 1:numVars
        jac[:,comp] = vf[comp].coeffs[2].coeffs[1:end]
    end

    return transpose(jac)
end
function jacobian{T<:Number,S<:Number}(vf::Array{TaylorN{T},1},vals::Array{S,1})
    R = promote_type(T,S)
    numVars = get_numvars()
    @assert length(vf) == numVars == length(vals)
    jac = Array(R, (numVars,numVars))

    for comp = 1:numVars
        @inbounds grad = gradient( vf[comp] )
        @inbounds for nv = 1:numVars
            jac[nv,comp] = evaluate(grad[nv], vals)
        end
    end

    return transpose(jac)
end

"""
```
    hessian(f)
    hessian(f, [vals])
```

Return the hessian matrix (jacobian of the gradient) of `f::TaylorN`,
evaluated at the vector `vals`. If `vals` is ommited, it is evaluated at
zero.
"""
hessian{T<:Number,S<:Number}(f::TaylorN{T}, vals::Array{S,1}) =
    (R = promote_type(T,S); jacobian( gradient(f), vals::Array{R,1}) )
hessian{T<:Number}(f::TaylorN{T}) = hessian( f, zeros(T, get_numvars()) )

## TODO: Integration...

## Evaluation
"""
    evaluate(a, vals)

Evaluate a `HomogeneousPolynomial` polynomial using Horner's rule (hand coded)
at `vals`.
"""
function evaluate{T<:Number,S<:Union{Real,Complex}}(a::HomogeneousPolynomial{T},
    vals::Array{S,1} )

    numVars = get_numvars()
    @assert length(vals) == numVars
    R = promote_type(T,S)
    suma = convert(TaylorN{R}, a)

    for nv = 1:numVars
        suma = horner(suma, (nv, vals[nv]))
    end

    return suma.coeffs[1].coeffs[1]
end

"""
    evaluate(a, [vals])

Evaluate a `TaylorN` polynomial using Horner's rule (hand coded) at `vals`.
If `vals` is ommitted, it is evaluated at zero.
"""
function evaluate{T<:Number,S<:Union{Real,Complex}}(a::TaylorN{T},
    vals::Array{S,1} )

    numVars = get_numvars()
    @assert length(vals) == numVars
    R = promote_type(T,S)
    suma = convert(TaylorN{R}, a)

    for nv = 1:numVars
        suma = horner(suma, (nv, vals[nv]))
    end

    return suma.coeffs[1].coeffs[1]
end

evaluate{T<:Number}(a::TaylorN{T}) = evaluate(a, zeros(T, get_numvars()))

## Evaluates HomogeneousPolynomials and TaylorN on a val of the nv variable
## using Horner's rule on the nv variable
function horner{T<:Number,S<:Union{Real,Complex}}(a::HomogeneousPolynomial{T},
    b::Tuple{Int,S} )

    nv, val = b
    numVars = get_numvars()
    @assert 1 <= nv <= numVars
    R = promote_type(T,S)
    @inbounds indTb = coeff_table[a.order+1]
    suma = TaylorN(zero(R), a.order)

    # Horner's rule on the nv variable
    for ord = a.order : -1 : 0
        suma_ord = TaylorN(zero(R), a.order)
        posOrd = order_posTb(a.order, nv, ord)
        neworder = a.order-ord
        for pos in posOrd
            c = a.coeffs[pos]
            iIndices = copy(indTb[pos])
            iIndices[nv] = 0
            kdic = in_base(get_order(), iIndices)
            newpos = pos_table[neworder+1][kdic]
            zhp = HomogeneousPolynomial(zero(R), neworder)
            zhp.coeffs[newpos] = a.coeffs[pos]
            suma_ord += TaylorN(zhp, a.order)
        end
        if ord == a.order
            suma += suma_ord
        else
            suma = suma*val + suma_ord
        end
    end

    return suma
end

function horner{T<:Number,S<:Union{Real,Complex}}(a::TaylorN{T},
    b::Tuple{Int,S} )

    nv, val = b
    @assert 1 <= nv <= get_numvars()
    R = promote_type(T,S)

    suma = TaylorN(zero(R), a.order)
    for ord = a.order:-1:0
        @inbounds polH = a.coeffs[ord+1]
        suma += horner( polH, b)
    end
    suma
end

"""
    order_posTb(order, nv, ord)

Return a vector with the positions, in a `HomogeneousPolynomial` of
order `order`, where the variable `nv` has order `ord`.
"""
function order_posTb(order::Int, nv::Int, ord::Int)
    @assert order <= get_order()
    @inbounds indTb = coeff_table[order+1]
    @inbounds num_coeffs = size_table[order+1]
    posV = Int[]
    for pos = 1:num_coeffs
        @inbounds indTb[pos][nv] != ord && continue
        push!(posV, pos)
    end
    posV
end


text:
include("../julia/ProgramTunnel_fs.jl")

using .ProgramTunnel_fs
using Formatting

PTI = ProgramTunnelInfo(recv_first_sleep=2.0)

println(PTI)

for i=1:100
    println(format("[{:3d}]", i))
    msg = recvText(PTI)
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# This file contains the definition for the Lines struct.

mutable struct Lines
    id
    from
    to
	reactance
	pmax

    function Lines(lines_df::DataFrame)

		# Zbase = (Vbase*1e3)^2/(Sbase*1e6)

		L = lines_df[!, 1]                                    # U = unit
        lines_dict = df_to_dict_with_id(lines_df)

		# lines_dict[:reactance] = lines_dict[:reactance] ./
		# 	lines_dict[:circuits] ./ Zbase

        return new(L,
                   lines_dict[:from],
                   lines_dict[:to],
				   lines_dict[:reactance],
                   lines_dict[:pmax],
				   )
    end
end


text:
module Mathepia

include("OtherTools/utils.jl")

export media_impact_fun

end


text:
module TSX

# https://gcc.gnu.org/onlinedocs/gcc/x86-transactional-memory-intrinsics.html

export @xabort, xbegin, xend, xtest, xabort, xabort_code, XBEGIN_STARTED, XABORT_EXPLICIT

const XBEGIN_STARTED = Int32(-1)
const XABORT_EXPLICIT = Int32(1 << 0)
const XABORT_RETRY = Int32(1 << 1)
const XABORT_CONFLICT = Int32(1 << 2)
const XABORT_CAPACITY = Int32(1 << 3)
const XABORT_DEBUG = Int32(1 << 4)
const XABORT_NESTED = Int32(1 << 5)

is_explicit(code) = (code & XABORT_EXPLICIT) != 0
is_retry(code) = (code & XABORT_RETRY) != 0
is_conflict(code) = (code & XABORT_CONFLICT) != 0

@inline function xbegin()
    Base.llvmcall(
        (
            """
            define i32 @entry() #0 {
            top:
                %0 = call i32 @llvm.x86.xbegin() #1
                ret i32 %0
            }

            declare i32 @llvm.x86.xbegin() #1

            attributes #0 = { alwaysinline }
            attributes #1 = { nounwind }
            """,
            "entry",
        ),
        Int32,
        Tuple{},
    )
end

@inline function xend()
    Base.llvmcall(
        (
            """
            define void @entry() #0 {
            top:
                call void @llvm.x86.xend() #1
                ret void
            }

            declare void @llvm.x86.xend() #1

            attributes #0 = { alwaysinline }
            attributes #1 = { nounwind }
            """,
            "entry",
        ),
        Cvoid,
        Tuple{},
    )
end

@inline function xtest()
    Base.llvmcall(
        (
            """
            define i32 @entry() #0 {
            top:
                %0 = call i32 @llvm.x86.xtest() #1
                ret i32 %0
            }

            declare i32 @llvm.x86.xtest() #1

            attributes #0 = { alwaysinline }
            attributes #1 = { nounwind }
            """,
            "entry",
        ),
        Int32,
        Tuple{},
    )
end

macro xabort(status::Signed)
    typemin(Int8) <= status <= typemax(Int8) || error("require an 8-bit signed integer")
    :(xabort(Val{$status}()))
end

@generated function xabort(::Val{status}) where {status}
    (status isa Signed && typemin(Int8) <= status <= typemax(Int8)) ||
        return :(error("require an 8-bit signed integer"))
    IR = (
        """
        define void @entry() #0 {
        top:
            call void @llvm.x86.xabort(i8 $status) #1
            ret void
        }

        declare void @llvm.x86.xabort(i8) #1

        attributes #0 = { alwaysinline }
        attributes #1 = { nounwind }
        """,
        "entry",
    )
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($IR, Cvoid, Tuple{})
    end
end

@inline xabort_code(status::Int32) = ((status >> 24) & 0xff) % Int8

end  # module TSX


text:
"""
SamplerState for preconditioned stochastic gradient Langevin dynamics with Adagrad following

Li, C., Chen, C., Carlson, D., & Carin., L. (2016). Preconditioned Stochastic Gradient Langevin Dynamics for Deep Neural Networks. AAAI. Retrieved from http://arxiv.org/abs/1512.07666

Fields:

    - x: the current state
    - preconditioner: current preconditioner

    - t: iteration counter
    - stepsize: s.stepsize(s.t) gives the stepsize to be used in iteration s.t
    - minmass

Constructors:

    - construct using `SGLDAdagradState(x,stepsize::Function)` or `SGLDAdagradState(x,stepsize::Float64)` for a fixed stepsize
"""
type SGLDAdagradState <: SamplerState
  x::Array{Float64,1}
  preconditioner::Array{Float64,1}

  t::Int64
  stepsize::Function
  lambda::Float64

  function SGLDAdagradState(x::Array{Float64,1},
    stepsize::Float64;
    preconditioner=zeros(size(x)),
    lambda=1e-8)
    @assert size(x) == size(preconditioner)
    s(niters) = stepsize
    new(x,preconditioner,0,s,lambda)
  end
  function SGLDAdagradState(x::Array{Float64,1},
    stepsize::Function;
    preconditioner=zeros(size(x)),
    lambda=1e-8)
    @assert size(x) == size(preconditioner)
    new(x,preconditioner,0,stepsize,lambda)
  end
end

function sample!(s::SGLDAdagradState, grad::Function)
    # increment iteration
    s.t += 1

    # improve readability
    t = s.t
    x = s.x
    lambda = s.lambda
    precond = s.preconditioner
    stepsize = s.stepsize(t)

    # update preconditioner
    g = grad(x)
    if t == 1
        precond[:] = g.*g
    else
        add_inplace!(1,g.*g,precond)
    end
    m = lambda + sqrt(precond)
    # gradient step (inplace addition)
    add_inplace!(stepsize,g./m,x)

    # injected noise (inplace addition)
    add_inplace!(sqrt(2.*stepsize),randn(length(x))./sqrt(m),x)
    s
end


text:
# alex's funcion
c=0.5
c1mlogc=c*(1-log(c))

function alex(x::Array)
    ind1=x.>-c
    ind2=x.<=-c
    tmp=zeros(size(x))
    tmp[ind1]=x[ind1]
    tmp[ind2]=-log(-x[ind2]).*c-c1mlogc
    ind1=ind2=0
    return tmp
end

function gradalex(x::Array)
    ind1=x.>-c
    ind2=x.<=-c
    tmp=zeros(size(x))
    tmp[ind1]=1.0
    tmp[ind2]=-c./x[ind2]
    ind1=ind2=0
    return tmp
end

Falex(x::Array{Float64,2})=( alex(x) ,[]);

function Falex_inplace(value,auxvalue,x::Array{Float64,2})
    ind1=x.>-c
    ind2=x.<=-c
    tmp=zeros(size(x))
    tmp[ind1]=x[ind1]+c1mlogc
    tmp[ind2]=-log(-x[ind2]).*c
    ind1=ind2=0
    copy!(value,tmp)
end

Dalex(derivativeIDX,f_c,faux_c,grad_c,grad_n,x)=axpy!(1.0,grad_c.*gradalex(x),grad_n)

if PROC=="GPU" # TODO
    #Falex(x::CudaArray)=(alex(x),[])
    Falex_inplace(value,aux,x::CudaArray)=alex!(x,value)

    Dalex(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray)=gradalex!(grad_c,x,grad_n)
end

Derivative[Falex]=Dalex
Inplace[Falex]=Falex_inplace

alex(A::ADnode)=ADnode(Falex,A)
export alex

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
"""
    Module LibDeflate

`LibDeflate` provides Julia bindings for the C library `libdeflate`. The C library,
and the corresponding Julia package, contain highly optimized code for compressing and
uncompressing data using the DEFLATE algorithm, including gzip, or zlib formats.
"""
module LibDeflate

using libdeflate_jll

"""
    Module LibDeflateErrors

Dummy module to contain the variants of the `LibDeflateError` enum.
"""
module LibDeflateErrors

@enum LibDeflateError::UInt8 begin
    deflate_bad_payload
    deflate_output_too_short
    deflate_insufficient_space
    gzip_header_too_short
    gzip_bad_magic_bytes
    gzip_not_deflate
    gzip_string_not_null_terminated
    gzip_null_in_string
    gzip_bad_header_crc16
    gzip_bad_crc32
    gzip_extra_too_long
    gzip_bad_extra
    zlib_input_too_short
    zlib_not_deflate
    zlib_wrong_window_size
    zlib_needs_compression_dict
    zlib_bad_header_check
    zlib_bad_adler32
    zlib_insufficient_space
end

@doc """
    LibDeflateError

A `UInt8` enum representing that LibDeflate encountered an error. The numerical value
of the errors are not stable across non-breaking releases, but their names are.
Code checking for specific errors should check by e.g. ` == LibDeflateErrors.gzip_not_deflate`.
Successful operations will never return a `LibDeflateError`.
"""
LibDeflateError

export LibDeflateError
end # module

using .LibDeflateErrors

const DEFAULT_COMPRESSION_LEVEL = 6

# Must be mutable for the GC to be able to interact with it
"""
    Decompressor()

Create an object which can decompress using the DEFLATE algorithm.
The same decompressor cannot be used by multiple threads at the same time.
To parallelize decompression, create multiple instances of `Decompressor`
and use one for each thread.
Creating this object allocates, so when decompressing multiple blocks, keep
the same decompressor in memory rather than making one for each block.

See also: [`decompress!`](@ref), [`unsafe_decompress!`](@ref)
"""
mutable struct Decompressor
    actual_nbytes_ret::UInt
    ptr::Ptr{Nothing}
end

Base.unsafe_convert(::Type{Ptr{Nothing}}, x::Decompressor) = x.ptr

function Decompressor()
    decompressor = Decompressor(0,
        ccall(
            (:libdeflate_alloc_decompressor, libdeflate),
            Ptr{Nothing},
            ()
        )
    )
    finalizer(free_decompressor, decompressor)
    return decompressor
end

function free_decompressor(decompressor::Decompressor)
    ccall(
        (:libdeflate_free_decompressor, libdeflate),
        Nothing,
        (Ptr{Nothing},),
        decompressor
    )
    return nothing
end

"""
    Compressor(compresslevel::Int=$(DEFAULT_COMPRESSION_LEVEL))

Create an object which can compress using the DEFLATE algorithm. `compresslevel`
can be from 1 (fast) to 12 (slow), and defaults to $(DEFAULT_COMPRESSION_LEVEL).
The same compressor cannot be used by multiple threads at the same time.
To parallelize compression, create multiple instances of `Compressor` and use one for each thread.
Creating this object allocates, so when compressing multiple blocks, keep
the same compressor in memory rather than making one for each block.

See also: [`compress!`](@ref), [`unsafe_compress!`](@ref)
"""
mutable struct Compressor
    level::Int
    ptr::Ptr{Nothing}
end

Base.unsafe_convert(::Type{Ptr{Nothing}}, x::Compressor) = x.ptr

function Compressor(compresslevel::Integer=DEFAULT_COMPRESSION_LEVEL)
    compresslevel in 1:12 || throw(ArgumentError("Compresslevel must be in 1:12"))
    ptr = ccall(
        (:libdeflate_alloc_compressor, libdeflate),
        Ptr{Nothing},
        (Csize_t,),
        compresslevel
    )
    compressor = Compressor(compresslevel, ptr)
    finalizer(free_compressor, compressor)
    return compressor
end

# Called by the garbage collecter, do not use manually
function free_compressor(compressor::Compressor)
    ccall((:libdeflate_free_compressor, libdeflate), Nothing, (Ptr{Nothing},), compressor)
    return nothing
end

# Compression and decompression functions

# Raw C call - do not export this
function _unsafe_decompress!(
    decompressor::Decompressor,
    out_ptr::Ptr,
    out_len::Integer,
    in_ptr::Ptr,
    inlen::Integer,
    nptr::Ptr
)::Union{LibDeflateError, Nothing}
    status = ccall(
        (:libdeflate_deflate_decompress, libdeflate),
        Csize_t,
        (Ptr{Nothing}, Ptr{UInt8}, Csize_t, Ptr{UInt8}, Csize_t, Ptr{UInt}),
        decompressor, in_ptr, inlen, out_ptr, out_len, nptr
    )
    if status == Cint(1)
        return LibDeflateErrors.deflate_bad_payload
    elseif status == Cint(2)
        return LibDeflateErrors.deflate_output_too_short
    elseif status == Cint(3)
        return LibDeflateErrors.deflate_insufficient_space
    else
        return nothing
    end
end

"""
    unsafe_decompress!(
        s::IteratorSize, ::Decompressor,
        out_ptr::Ptr, n_out::Integer,
        in_ptr::Ptr, n_in::Integer
    )::Union{Int, LibDeflateError}

Decompress `n_in` bytes from `in_ptr` to `out_ptr` using the DEFLATE algorithm,
returning the number of decompressed bytes or the error encounted.
`s` should be whether you know the decompressed size or not.

If `s` isa `Base.HasLength`, the number of decompressed bytes is given as `n_out`.
This is more efficient, but will fail if the number is not correct.

If `s` isa `Base.SizeUnknown`, pass the size in bytes of the available space at the output
to `n_out`.

See also: [`decompress!`](@ref)
"""
function unsafe_decompress! end

function unsafe_decompress!(
    ::Base.HasLength,
    decompressor::Decompressor,
    out_ptr::Ptr,
    n_out::Integer,
    in_ptr::Ptr,
    n_in::Integer
)::Union{LibDeflateError, Int}
    y = _unsafe_decompress!(decompressor, out_ptr, n_out, in_ptr, n_in, C_NULL)
    y isa LibDeflateError ? y : Int(n_out)
end

function unsafe_decompress!(
    ::Base.SizeUnknown,
    decompressor::Decompressor,
    out_ptr::Ptr,
    n_out::Integer,
    in_ptr::Ptr,
    n_in::Integer
)::Union{LibDeflateError, Int}
    y = GC.@preserve decompressor begin
        retptr = pointer_from_objref(decompressor)
        _unsafe_decompress!(decompressor, out_ptr, n_out, in_ptr, n_in, retptr)
    end
    y isa LibDeflateError ? y : (decompressor.actual_nbytes_ret % Int)
end

"""
    decompress!(
        ::Decompressor, out_data::Array, in_data, [n_out::Integer]
    )::Union{LibDeflateError, Int}

Use the passed `Decompressor` to decompress the data at `in_data` into the
first bytes of `out_data` using the DEFLATE algorithm,
returning the number of written bytes to the output, or a `LibDeflateError`.
Data must fit in `out_data`. `in_data` must implement `sizeof` and `pointer`.

If the decompressed size is known beforehand, pass it as `n_out`. This will increase
performance, but will fail if it is wrong.
"""
function decompress! end

# Decompress method with length known (preferred)
function decompress!(
    decompressor::Decompressor,
    out_data::Array,
    in_data,
    n_out::Integer
)::Union{LibDeflateError, Int}
    sizeof(out_data) < n_out && return LibDeflateErrors.deflate_insufficient_space
    GC.@preserve out_data in_data unsafe_decompress!(
        Base.HasLength(),
        decompressor,
        pointer(out_data),
        n_out,
        pointer(in_data),
        sizeof(in_data)
    )
end

# Decompress method with length unknown (not preferred)
function decompress!(
    decompressor::Decompressor,
    out_data::Array,
    in_data
)::Union{LibDeflateError, Int}
    GC.@preserve out_data in_data unsafe_decompress!(
        Base.SizeUnknown(),
        decompressor,
        pointer(out_data),
        sizeof(out_data),
        pointer(in_data),
        sizeof(in_data)
    )
end

"""
    unsafe_compress!(
        ::Compressor, out_ptr::Ptr, n_out::Integer, in_ptr::Ptr, n_in::Integer
    )::Union{Int. LibDeflateError}

Use the passed `Compressor` to compress `n_in` bytes from the pointer `in_ptr`
to the pointer `n_out` where there are `n_out` bytes of space to write to.

Return the number of written bytes to the output, or a `LibDeflateError`.

See also: [`compress!`](@ref)
"""
function unsafe_compress!(
    compressor::Compressor,
    out_ptr::Ptr,
    n_out::Integer,
    in_ptr::Ptr,
    n_in::Integer
)::Union{LibDeflateError, Int}
    bytes = ccall(
        (:libdeflate_deflate_compress, libdeflate),
        Csize_t,
        (Ptr{Nothing}, Ptr{UInt8}, Csize_t, Ptr{UInt8}, Csize_t),
        compressor, in_ptr, n_in, out_ptr, n_out
    )
    iszero(bytes) && return LibDeflateErrors.deflate_insufficient_space
    return bytes % Int
end

"""
    compress!(::Compressor, out_data::Array, in_data)::Union{LibDeflateError, Int}

Use the passed `Compressor` to compress `in_data` into the first
bytes of `out_data` using the DEFLATE algorithm.
Data must fit in `out_data`. `in_data` must implement `sizeof` and `pointer`.

Return the number of written bytes to the output, or a `LibDeflateError`.
"""
function compress!(
    compressor::Compressor,
    out_data::Array,
    in_data
)::Union{LibDeflateError, Int}
    GC.@preserve out_data in_data unsafe_compress!(
        compressor,
        pointer(out_data),
        sizeof(out_data),
        pointer(in_data),
        sizeof(in_data)
    )
end

"""
    unsafe_crc32(in_ptr::Ptr, n_in::Integer, start::UInt32)::UInt32

Calculate the crc32 checksum of the first `n_in` bytes of the pointer `in_ptr`,
with seed `start` (default is 0).
Note that crc32 is a different and slower algorithm than the `crc32c` provided
in the Julia standard library.

See also: [`crc32`](@ref)
"""
function unsafe_crc32(in_ptr::Ptr, n_in::Integer, start::UInt32=UInt32(0))
    return ccall(
        (:libdeflate_crc32, libdeflate),
        UInt32, (UInt32, Ptr{UInt8}, Csize_t),
        start, in_ptr, n_in
    )
end

"""
    crc32(data, start=UInt32(0))::UInt32

Calculate the crc32 checksum of `data` and seed `start` (0 by default).
`data` must implement `pointer` and `sizeof`.
Note that crc32 is a different and slower algorithm than the `crc32c` provided
in the Julia standard library.

See also: [`unsafe_crc32`](@ref)
"""
function crc32(data, start::UInt32=UInt32(0))
    GC.@preserve data unsafe_crc32(pointer(data), sizeof(data), start)
end

"""
    unsafe_adler32(data, start=UInt32(1))::UInt32

Calculate the adler32 checksum of the first `n_in` of the pointer `in_ptr`,
with seed `start` (default is 1).

See also: [`adler32`](@ref)
"""
function unsafe_adler32(in_ptr::Ptr, n_in::Integer, start::UInt32=UInt32(1))
    return ccall(
        (:libdeflate_adler32, libdeflate),
        UInt32, (UInt32, Ptr{UInt8}, Csize_t),
        start, in_ptr, n_in
    )
end

"""
    adler32(data, start=UInt32(1))::UInt32

Calculate the adler32 checksum of the byte vector `data` and seed `start` (1 by default).
`data` must implement `pointer` and `sizeof`.

See also: [`unsafe_adler32`](@ref)
"""
function adler32(data, start::UInt32=UInt32(1))
    GC.@preserve data unsafe_adler32(pointer(data), sizeof(data), start)
end

include("gzip.jl")
include("zlib.jl")

export Decompressor,
       Compressor,

       LibDeflateErrors,
       LibDeflateError,

       unsafe_decompress!,
       decompress!,
       unsafe_compress!,
       compress!,

       unsafe_gzip_decompress!,
       gzip_decompress!,
       unsafe_gzip_compress!,
       gzip_compress!,

       unsafe_zlib_decompress!,
       zlib_decompress!,
       unsafe_zlib_compress!,
       zlib_compress!,

       unsafe_crc32,
       crc32,
       unsafe_adler32,
       adler32,

       unsafe_parse_gzip_header,
       parse_gzip_header,
       is_valid_extra_data

end # module


text:
module MonteCarlo

export simMonteCarlo!, simMonteCarloMulti!, runSim
using TypeDef, Tools, Helper, Solve, UserFileIO
using Printf, JLD, Distributed

function simMonteCarlo!(
  b1  ::Int64,
  b2  ::Int64,
  dat ::datType,
  dim ::dimType,
  env ::envType,
  mat ::matType,
  par0::parType,
  set ::setType,
  indx::indxType,
  mask::maskType,
  vInt::Bool=false)

  tic = time()
  eqa = allocEqArray(b1,b2,dim)
  for b=b1:b2
    runSim!(eqa,b,b1,dat,dim,env,mat,par0,set,indx,mask,vInt)
  end
  @printf("Total %.0f sec elapsed.\n",time()-tic);
  saveJLD!(eqa,env,env.fileOutp)

  return eqa
end

function simMonteCarloMulti!(
  b1  ::Int64,
  b2  ::Int64,
  dat ::datType,
  dim ::dimType,
  env ::envType,
  mat ::matType,
  par0::parType,
  set ::setType,
  indx::indxType,
  mask::maskType,
  vInt::Bool=false)

  tic = time()
  eqa = allocEqArray(b1,b2,dim)
  f =
  let dat=deepcopy(dat),mat=deepcopy(mat),indx=deepcopy(indx)
    n->runSim(n,b1,dat,dim,env,mat,par0,set,indx,mask,vInt)
  end
  pout = pmap(f, b1:b2)
  for b in eachindex(pout)
    eqa[b] = pout[b]
  end
  @printf("Total %.0f sec elapsed.\n",time()-tic);
  saveJLD!(eqa,env,env.fileOutp)

  return eqa
end

function runSim(
  b   ::Int64,
  b1  ::Int64,
  dat ::datType,
  dim ::dimType,
  env ::envType,
  mat ::matType,
  par0::parType,
  set ::setType,
  indx::indxType,
  mask::maskType,
  vInt::Bool=false)

  toc = time()
  par = sampleParam!(b,dim,env,mat,par0)
  # par = par0
  sim = allocateSim(dim)
  eq  = allocateEquil(dim)
  initialize!(mat,dat,dim,par,set)
  marginalCost!(eq,dat,dim,env,mat,par,set,sim,indx,mask,vInt)
  bertrandSolve!(eq,dat,dim,env,mat,par,set,sim,indx,mask,vInt)
  printLog!(b,toc,eq,env.fileLog)
  @printf("σₐ[%3d]: %.4f complete in %.0f sec.\n",b,par.θ₂[1],time()-toc)

  return eq

end

function runSim!(
  eqa ::Array{Array{equilType,1},1},
  b   ::Int64,
  b1  ::Int64,
  dat ::datType,
  dim ::dimType,
  env ::envType,
  mat ::matType,
  par0::parType,
  set ::setType,
  indx::indxType,
  mask::maskType,
  vInt::Bool=false)

  toc = time()
  par = sampleParam!(b,dim,env,mat,par0)
  # par = par0
  sim = allocateSim(dim)
  bid = b - b1 + 1
  initialize!(mat,dat,dim,par,set)
  marginalCost!(eqa[bid],dat,dim,env,mat,par,set,sim,indx,mask,vInt)
  bertrandSolve!(eqa[bid],dat,dim,env,mat,par,set,sim,indx,mask,vInt)
  # saveCSV!(b,eqa[bid],dim,env)
  # saveJLD!(eqa,env,env.fileOutp)
  printLog!(b,toc,eqa[bid],env.fileLog)
  @printf("σₐ[%3d]: %.4f complete in %.0f sec.\n",b,par.θ₂[1],time()-toc)

end

end

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand


text:
function isone(a::AlgAssAbsOrdFracIdl)
  a = simplify!(a)
  return isone(denominator(a, copy = false)) && isone(numerator(a, copy = false))
end

function one(S::AlgAssAbsOrdFracIdlSet)
  return frac_ideal(order(S), ideal(order(S), one(order(S))), fmpz(1))
end

function one(I::AlgAssAbsOrdFracIdl)
  return frac_ideal(order(I), ideal(order(I), one(order(I))), fmpz(1))
end

function Base.copy(I::AlgAssAbsOrdFracIdl)
  return I
end

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

function show(io::IO, a::AlgAssAbsOrdFracIdl)
  print(io, "Fractional Ideal of ")
  show(IOContext(io, :compact => true), order(a))
  println(io, " with basis matrix ")
  print(io, basis_matrix(a, copy = false))
end

###############################################################################
#
#  Deepcopy
#
###############################################################################

function Base.deepcopy_internal(a::AlgAssAbsOrdFracIdl, dict::IdDict)
  b = typeof(a)(order(a), numerator(a), denominator(a))
  for i in fieldnames(typeof(a))
    if isdefined(a, i)
      if i != :order && i != :num && i != :den
        setfield!(b, i, Base.deepcopy_internal(getfield(a, i), dict))
      end
    end
  end
  return b
end

###############################################################################
#
#  Field access
#
###############################################################################

order(a::AlgAssAbsOrdFracIdl) = a.order

function numerator(a::AlgAssAbsOrdFracIdl; copy::Bool = true)
  if !isdefined(a, :num)
    if !isdefined(a, :basis_matrix)
      error("No basis matrix defined")
    end
    a.num = ideal(order(a), numerator(basis_matrix(a)))
  end
  if copy
    return deepcopy(a.num)
  else
    return a.num
  end
end

function denominator(a::AlgAssAbsOrdFracIdl; copy::Bool = true)
  if !isdefined(a, :den)
    if !isdefined(a, :basis_matrix)
      error("No basis matrix defined")
    end
    a.den = denominator(basis_matrix(a))
  end
  if copy
    return deepcopy(a.den)
  else
    return a.den
  end
end

function basis_matrix(a::AlgAssAbsOrdFracIdl; copy::Bool = true)
  if !isdefined(a, :basis_matrix)
    if !isdefined(a, :num) || !isdefined(a, :den)
      error("Numerator and/or denominator not defined")
    end
    a.basis_matrix = FakeFmpqMat(basis_matrix(a.num, copy = false), a.den)
  end
  if copy
    return deepcopy(a.basis_matrix)
  else
    return a.basis_matrix
  end
end

function basis_mat_inv(a::AlgAssAbsOrdFracIdl; copy::Bool = true)
  if !isdefined(a, :basis_mat_inv)
    a.basis_mat_inv = inv(basis_matrix(a, copy = false))
  end
  if copy
    return deepcopy(a.basis_mat_inv)
  else
    return a.basis_mat_inv
  end
end

################################################################################
#
#  Basis
#
################################################################################

function basis(a::AlgAssAbsOrdFracIdl)
  B = basis_matrix(a, copy = false)
  O = order(a)
  A = algebra(O)
  d = dim(A)
  Oba = basis(O, copy = false)
  res = Array{elem_type(A)}(undef, d)

  for i in 1:d
    z = A()
    for j in 1:d
      z = z + B.num[i, j]*elem_in_algebra(Oba[j], copy = false)
    end
    z = divexact(z, A(B.den))
    res[i] = z
  end

  return res
end

################################################################################
#
#  Construction
#
################################################################################

function frac_ideal(O::AlgAssAbsOrd{S, T}, M::FakeFmpqMat, M_in_hnf::Bool = false) where {S, T}
  !M_in_hnf ? M = hnf(M, :lowerleft) : nothing

  return AlgAssAbsOrdFracIdl{S, T}(O, M)
end

function frac_ideal(O::AlgAssAbsOrd{S, T}, a::AlgAssAbsOrdIdl{S, T}, d::fmpz) where {S, T}
  return AlgAssAbsOrdFracIdl{S, T}(O, a, d)
end

function frac_ideal(O::AlgAssAbsOrd{S, T}, a::T) where {S, T}
  d = denominator(a, O)
  return frac_ideal(O, ideal(O, O(d*a)), d)
end

*(O::AlgAssAbsOrd{S, T}, a::T) where {S, T} = frac_ideal(O, a)
*(a::T, O::AlgAssAbsOrd{S, T}) where {S, T} = frac_ideal(O, a)

function frac_ideal_from_z_gens(O::AlgAssAbsOrd{S, T}, b::Vector{T}) where {S, T}
  d = degree(O)
  den = lcm([ denominator(bb, O) for bb in b ])
  num = ideal_from_z_gens(O, [ O(den*bb) for bb in b ])
  return frac_ideal(O, num, den)
end

################################################################################
#
#  Arithmetic
#
################################################################################

function +(a::AlgAssAbsOrdFracIdl{S, T}, b::AlgAssAbsOrdFracIdl{S, T}) where {S, T}
  d = lcm(denominator(a, copy = false), denominator(b, copy = false))
  ma = div(d, denominator(a, copy = false))
  mb = div(d, denominator(b, copy = false))
  return frac_ideal(order(a), numerator(a, copy = false)*ma + numerator(b, copy = false)*mb, d)
end

function *(a::AlgAssAbsOrdFracIdl{S, T}, b::AlgAssAbsOrdFracIdl{S, T}) where {S, T}
  d = denominator(a, copy = false)*denominator(b, copy = false)
  return frac_ideal(order(a), numerator(a, copy = false)*numerator(b, copy = false), d)
end

^(A::AlgAssAbsOrdFracIdl, e::Int) = Base.power_by_squaring(A, e)
^(A::AlgAssAbsOrdFracIdl, e::fmpz) = Base.power_by_squaring(A, BigInt(e))

################################################################################
#
#  Ad hoc binary operations
#
################################################################################

function +(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdFracIdl{S, T}) where {S, T}
  c = a*denominator(b, copy = false) + numerator(b, copy = false)
  return frac_ideal(order(a), c, denominator(b))
end

+(a::AlgAssAbsOrdFracIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where {S, T} = b + a

function *(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdFracIdl{S, T}) where {S, T}
  return typeof(b)(order(a), a*numerator(b, copy = false), denominator(b))
end

function *(a::AlgAssAbsOrdFracIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where {S, T}
  return typeof(a)(order(a), numerator(a, copy = false)*b, denominator(a))
end

function *(x::T, a::AlgAssAbsOrdFracIdl{S, T}) where {S, T}
  O = order(a)
  return frac_ideal(O, x)*a
end

function *(a::AlgAssAbsOrdFracIdl{S, T}, x::T) where {S, T}
  O = order(a)
  return a*frac_ideal(O, x)
end

function *(x::T, a::AlgAssAbsOrdIdl{S, T}) where {S, T}
  O = order(a)
  return frac_ideal(O, x)*a
end

function *(a::AlgAssAbsOrdIdl{S, T}, x::T) where {S, T}
  O = order(a)
  return a*frac_ideal(O, x)
end

################################################################################
#
#  Equality
#
################################################################################

function ==(A::AlgAssAbsOrdFracIdl, B::AlgAssAbsOrdFracIdl)
  order(A) !== order(B) && return false
  return basis_matrix(A, copy = false) == basis_matrix(B, copy = false)
end

################################################################################
#
#  Simplification
#
################################################################################

function simplify!(a::AlgAssAbsOrdFracIdl)
  b = basis_matrix(numerator(a, copy = false), copy = false)
  g = gcd(denominator(a, copy = false), content(b))

  if g != 1
    a.num = divexact(a.num, g)
    a.den = divexact(a.den, g)
  end

  return a
end

################################################################################
#
#  Random elements
#
################################################################################

function rand(a::AlgAssAbsOrdFracIdl, B::Int)
  z = rand(numerator(a, copy = false), B)
  return fmpq(1, denominator(a, copy = false))*elem_in_algebra(z, copy = false)
end

################################################################################
#
#  Inclusion of elements
#
################################################################################

function in(x::AbsAlgAssElem, a::AlgAssAbsOrdFracIdl)
  O = order(a)
  M = zero_matrix(FlintZZ, 1, degree(O))
  t = FakeFmpqMat(M)
  elem_to_mat_row!(t.num, 1, t.den, x)
  v = t*basis_mat_inv(O, copy = false)
  v = v*basis_mat_inv(a, copy = false)

  return v.den == 1
end

################################################################################
#
#  Colon
#
################################################################################

function colon(I::AlgAssAbsOrdFracIdl, J::AlgAssAbsOrdFracIdl)
  # Let I = a/c and J = b/d, a and b integral ideals, c, d \in Z, then
  # \{ x \in K | xJ \subseteq I \} = \{ x \in K | xcb \subseteq da \}

  II = numerator(I, copy = false)*denominator(J, copy = false)
  JJ = numerator(J, copy = false)*denominator(I, copy = false)
  return colon(II, JJ)
end

################################################################################
#
#  Move ideals to another order
#
################################################################################

function extend(I::AlgAssAbsOrdFracIdl, O::AlgAssAbsOrd)
  J = extend(numerator(I, copy = false), O)
  return frac_ideal(O, J, denominator(I, copy = false))
end

*(I::AlgAssAbsOrdFracIdl, O::AlgAssAbsOrd) = extend(I, O)
*(O::AlgAssAbsOrd, I::AlgAssAbsOrdFracIdl) = extend(I, O)

function _as_frac_ideal_of_smaller_order(O::AlgAssAbsOrd, I::AlgAssAbsOrdIdl)
  M = basis_matrix(I, copy = false)
  M = M*basis_matrix(order(I), copy = false)*basis_mat_inv(O, copy = false)
  return frac_ideal(O, M)
end

function _as_frac_ideal_of_smaller_order(O::AlgAssAbsOrd, I::AlgAssAbsOrdFracIdl)
  J = _as_frac_ideal_of_smaller_order(O, numerator(I, copy = false))
  return algebra(O)(fmpq(1, denominator(I, copy = false)))*J
end

################################################################################
#
#  Ideal Set
#
################################################################################

function show(io::IO, a::AlgAssAbsOrdFracIdlSet)
  print(io, "Set of fractional ideals of $(order(a))")
end

order(a::AlgAssAbsOrdFracIdlSet) = a.order

parent(I::AlgAssAbsOrdFracIdl) = AlgAssAbsOrdFracIdlSet(order(I))

function FracIdealSet(O::AlgAssAbsOrd)
   return AlgAssAbsOrdFracIdlSet(O)
end

elem_type(::Type{AlgAssAbsOrdFracIdlSet{S, T}}) where {S, T} = AlgAssAbsOrdFracIdl{S, T}

elem_type(::AlgAssAbsOrdFracIdlSet{S, T}) where {S, T} = AlgAssAbsOrdFracIdl{S, T}

parent_type(::Type{AlgAssAbsOrdFracIdl{S, T}}) where {S, T} = AlgAssAbsOrdFracIdlSet{S, T}


text:
export ResidueField

function ResidueField(O::NfMaxOrd, P::NfMaxOrdIdl)
  @assert !is_index_divisor(O, minimum(P))
  @assert has_2_elem(P) && is_prime_known(P)

  gtwo = P.gen_two

  f = nf(O).pol
  g = parent(f)(elem_in_nf(gtwo))

  if nbits(P.gen_one) < 64
    R = ResidueRing(FlintZZ, P.gen_one)

    Zy, y = PolynomialRing(FlintZZ, "y")
    Rx, x = PolynomialRing(R, "x")::Tuple{NmodPolyRing, nmod_poly}

    gmodp = Rx(Zy(g))::nmod_poly
    fmodp = Rx(Zy(f))::nmod_poly

    h = gcd(gmodp,fmodp)

    F = FqNmodFiniteField(h, :$)

    #return F, Mor(O, F, gen(F))
    return F, Mor(O, F, h)
  else
    error("Not yet implemented")
  end
end

function disc_log(x::fq_nmod, g::fq_nmod)
  F = parent(x)
  q = order(F)
  h = g
  for i in 1:Int(q)-2
    if h == x
      return i
    end
    h = h*g
  end
  return i
end

function disc_log(x::fq_nmod, g::fq_nmod, p::Int)
  F = parent(x)
  q = order(F)
  @assert mod(q - 1, p) == 0
  l = div(q - 1, p)
  x = x^l
  g = g^l
  h = deepcopy(g)
  i = 1
  while i < p + 2
    if h == x
      return i
    end
    h = h*g
    i = i + 1
  end
  error("something odd")
end


text:
@block AnshulSinghvi ["geo", "geomakie"] begin

    using GeoMakie

    @cell "Geographical axes" [geoaxis, projection] begin

        projections = [
            Projection("+proj=longlat") => "Standard lon-lat",
            Projection("+proj=robin")   => "Robinson",
            Projection("+proj=moll")    => "Molleweide",
            Projection("+proj=vandg")   => "van der Grinten I",
            Projection("+proj=wintri")  => "Winkel Tripel",
        ]

        geoaxis_scenes = Scene[]

        for (projection, titlestr) in projections
            push!(geoaxis_scenes, title(geoaxis(-180, 180, -90, 90; show_axis = false, crs = (dest = projection,), scale_plot = false), titlestr))
        end

        vbox(geoaxis_scenes...; parent = Scene(resolution = (1440, 200)))

    end

    @cell "Earth and Coastlines" [earth, coastlines, geoaxis] begin
        using GeoMakie

        coastlines(; crs = (dest = WinkelTripel(),), show_axis = false, scale_plot = false)
        earth!(; crs = (dest = WinkelTripel(),))
        geoaxis!(-180, 180, -90, 90; crs = (dest = WinkelTripel(),))
    end

    @cell "Air Particulates" [record, animation] begin
        using GeoMakie
        using GeoMakie: ImageMagick, Glob

        source = LonLat()
        dest = WinkelTripel()

        imgdir = observations("rgb/MYDAL2_M_AER_RA")

        imgpaths = sort(Glob.glob("*.PNG", imgdir)) # change this to your requirements

        img = ImageMagick.load(imgpaths[1])

        re = GeoMakie.date_regex("MYDAL2_M_AER_RA", "PNG")
        titletext = Node(join(match(re, basename(imgpaths[1])).captures, '-'))

        lons = LinRange(-179.5, 179.5, size(img)[2])
        lats = LinRange(89.5, -89.5, size(img)[1])

        xs = [lon for lat in lats, lon in lons]
        ys = [lat for lat in lats, lon in lons]

        Proj4.transform!(source, dest, vec(xs), vec(ys))

        scene = surface(xs, ys, zeros(size(xs)); color = img, shading = false, show_axis = false, scale_plot = false)
        surfplot = scene.plots[end]

        geoaxis!(scene, -180, 180, -90, 90; crs = (src = source, dest = dest,));

        titletext = Node("07/2016")

        fullsc = title(scene, titletext; fontsize = 40);

        record(fullsc, @replace_with_a_path(mp4), imgpaths; framerate = 10) do img
            year, month = match(re, img).captures
            surfplot.color = ImageMagick.load(img)
            titletext[] = "$month/$year"
        end

    end

    @cell "MakieLayout with GeoMakie" [layout] begin

        using MakieLayout

        lons = LinRange(-179.5, 179.5, 360)
        lats = LinRange(-89.5, 89.5, 180)

        field = [exp(cosd(l)) + 3(y/90) for l in lons, y in lats]

        source = LonLat()
        dest = WinkelTripel()

        xs, ys = xygrid(lons, lats)
        Proj4.transform!(source, dest, vec(xs), vec(ys))

        xmin, xmax = extrema(xs)
        ymin, ymax = extrema(ys)

        aspect_ratio = (ymax - ymin) / (xmax - xmin)

        scene, layout = layoutscene(
            resolution = (1300, 900)
        )

        layout[1, 1] = lsc = LScene(scene; scenekw = (show_axis = false, scale_plot = false))

        splot = surface!(lsc, xs, ys; color = field, shading = false, show_axis = false)

        geoaxis!(lsc, -180, 180, -90, 90; crs = (src = source, dest = dest,))

        coastlines!(lsc; crs = (src = source, dest = dest,))

        layout[1, 2] = LColorbar(scene, splot; label = "Arbitrary data", width = 30)

        colsize!(layout, 1, Relative(1))
        rowsize!(layout, 1, Aspect(1, aspect_ratio))

        layout[0, :] = LText(scene, "MakieLayout is cool!"; textsize = 40)

        return scene

    end

end


text:
using DynamicLinearModels
using LinearAlgebra
using RecipesBase
using Test

@testset "DynamicLinearModels" begin
    include("utility.jl")
    include("univariate.jl")
    include("weighted.jl")
    include("dynamic.jl")
end


text:
using GLWindow, GLAbstraction, GLFW, ModernGL, ImmutableArrays, WavefrontObj, FixedPointNumbers, Images, Colorant, Reactive
using GLPlot #toopengl
# window creation
window = createwindow("OBJ-Viewer", 1000, 1000, debugging = false)
cam = PerspectiveCamera(window.inputs, Vec3(2,2,0.5), Vec3(0.0))

# render objects creation
shader = TemplateProgram(Pkg.dir("GLPlot", "src", "shader", "standard.vert"), Pkg.dir("GLPlot", "src", "shader", "phongblinn.frag"))
RGBAU8 = AlphaColorValue{RGB{UFixed8}, UFixed8}
Colorant.rgba(r::Real, g::Real, b::Real, a::Real)    = AlphaColorValue(RGB{Float32}(r,g,b), Float32(a))
rgbaU8(r::Real, g::Real, b::Real, a::Real)  = AlphaColorValue(RGB{UFixed8}(r,g,b), ufixed8(a))
#GLPlot.toopengl{T <: AbstractRGB}(colorinput::Signal{T}) = toopengl(const_lift(x->AlphaColorValue(x, one(T)), RGBA{T}, colorinput))
tohsva(rgba)     = AlphaColorValue(convert(HSV, rgba.c), rgba.alpha)
torgba(hsva)     = AlphaColorValue(convert(RGB, hsva.c), hsva.alpha)
tohsva(h,s,v,a)  = AlphaColorValue(HSV(Float32(h), Float32(s), Float32(v)), Float32(a))
Base.convert{T <: AbstractAlphaColorValue}(typ::Type{T}, x::AbstractAlphaColorValue) = AlphaColorValue(convert(RGB{eltype(typ)}, x.c), convert(eltype(typ), x.alpha))

function unitGeometry{T}(geometry::Vector{Vector3{T}})
    assert(!isempty(geometry))

    xmin = typemax(T)
    ymin = typemax(T)
    zmin = typemax(T)

    xmax = typemin(T)
    ymax = typemin(T)
    zmax = typemin(T)

    for vertex in geometry
        xmin = min(xmin, vertex[1])
        ymin = min(ymin, vertex[2])
        zmin = min(zmin, vertex[3])

        xmax = max(xmax, vertex[1])
        ymax = max(ymax, vertex[2])
        zmax = max(zmax, vertex[3])
    end

    xmiddle = xmin + (xmax - xmin) / 2;
    ymiddle = ymin + (ymax - ymin) / 2;
    zmiddle = zmin + (zmax - zmin) / 2;
    scale = 2 / max(xmax - xmin, ymax - ymin, zmax - zmin);

    result = similar(geometry)

    for i = 1:length(result)
        result[i] = Vector3{T}((geometry[i][1] - xmiddle) * scale,
                               (geometry[i][2] - ymiddle) * scale,
                               (geometry[i][3] - zmiddle) * scale
                    );
    end

    return result
end

objpath = "Butterfly/Butterfly.obj"

assets_path = "Butterfly/"

obj = readObjFile(objpath, faceindextype=GLuint, vertextype=Float32, compute_normals = false, triangulate = false)
computeNormals!(obj, smooth_normals = true, override = false)
triangulate!(obj)

# center geometry
obj.vertices = unitGeometry(obj.vertices)

# load mtl files if present
materials = WavefrontMtlMaterial{Float32}[]

for mtllib in obj.mtllibs
    materials = [materials, readMtlFile( assets_path*mtllib, colortype=Float32 )]
end

render_objects = RenderObject[]

compiled_materials = Dict()

const light         = Vec3[Vec3(1.0,1.0,1.0), Vec3(0.1,0.1,0.1), Vec3(0.9,0.9,0.9), Vec3(20.0,20.0,20.0)]
const material      = Vec3[Vec3(1.0,1.0,1.0), Vec3(1.0,1.0,1.0), Vec3(1.0,1.0,1.0), Vec3(1.0,1.0,1.0)]
const tmaterial      = Vec4[Vec4(1.0,1.0,1.0,1.0), Vec4(1.0,1.0,1.0,1.0), Vec4(1.0,1.0,1.0,1.0), Vec4(1.0,1.0,1.0,1.0)]
const tmaterialused  = GLint[-1,-1,-1,-1]

for material_name in collect(keys(obj.materials))

    (vs, nvs, uvs, fcs) = compileMaterial(obj, material_name)

    # hack: invert normals for glabstraction
    nvs = -nvs

    # holding global references seems necessary here
    # compiled_materials[material_name] = (vs, nvs, uvs, fcs, lines)
    data = [
        :vertex          => GLBuffer(vs),
        :normal          => GLBuffer(nvs),
        :uv              => GLBuffer(uvs),
        :indexes         => indexbuffer(fcs),
        :view            => cam.view,
        :projection      => cam.projection,
        :normalmatrix    => cam.normalmatrix,
        :eyeposition     => cam.eyeposition,
        :model           => eye(Mat4),
        :material        => material,
        :light           => light,
        :textures_used   => tmaterialused,
    ]

    # search for a material with the given name
    texture_array = Any[]
    for mtl in materials
        if mtl.name == material_name
            data[:material]             = Vec3[mtl.diffuse, mtl.ambient, mtl.specular, Vec3(mtl.specular_exponent)]
            println(data[:material])
            if mtl.diffuse_texture != ""
                data[:textures_used][1] = length(texture_array) # insert texture array index
                push!(texture_array, imread(assets_path*"Texture/"*mtl.diffuse_texture).data)
            end
            if mtl.ambient_texture != ""
                data[:textures_used][2] = length(texture_array) # insert texture array index
                push!(texture_array, imread(assets_path*"Texture/"*mtl.ambient_texture).data)
            end
            if mtl.specular_texture != ""
                data[:textures_used][3] = length(texture_array) # insert texture array index
                push!(texture_array, imread(assets_path*"Texture/"*mtl.specular_texture).data)
            end
            #=
            if mtl.bump_texture != ""
                data[:textures_used] = true
                data[:tmaterialused][1] = 0 # insert texture array index
                push!(texture_array, imread(assets_path*mtl.diffuse_texture).data)
            end
            =#
            break
        end
    end

    if !isempty(texture_array)
        data[:textures_used] = true
        data[:texture_maps] = Texture(convert(Vector{Matrix{eltype(first(texture_array))}}, texture_array))
    else
        data[:texture_maps] = Texture(Matrix{RGBA{UFixed8}}[fill(rgbaU8(0,0,0,0), 1,1)])
    end

    ro = RenderObject(data, shader)

    postrender!(ro, render, ro.vertexarray)
    #postrender!(ro, render, ro.vertexarray, GL_LINES)

    push!(render_objects, ro)
end

# OpenGL setup
glClearColor(0.2,0.2,0.2,1)
glDisable(GL_CULL_FACE)
glEnable(GL_DEPTH_TEST)

const_lift(window.inputs[:framebuffer_size]) do wh
    glViewport(0,0,wh...)
end |> preserve
# Loop until the user closes the window
while !GLFW.WindowShouldClose(window.nativewindow)

  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

  # render materials separately
  for ro in render_objects
    render(ro)
  end

  yield() # this is needed for react to work

  GLFW.SwapBuffers(window.nativewindow)
  GLFW.PollEvents()
end

GLFW.Terminate()


text:
function get_energy( p::Param, every::Int , tlimit::Float64)

    #Get variables
    xmin        = p.xmin
    xmax        = p.xmax
    Nx          = p.xnodes

    ymin        = p.ymin
    ymax        = p.ymax
    Ny          = p.ynodes

    #Discretize space
    dx      = (xmax - xmin)/(Nx-1)
    x       = [xmin + (i-1)*dx for i in 1:Nx]

    dy      = (ymax - ymin)/(Ny-1)
    y       = [ymin + (i-1)*dy for i in 1:Ny]

    if(p.Boundaries == :periodic )
        dx      = (xmax - xmin)/(Nx)
        x       = [xmin + (i-1)*dx for i in 1:Nx]

        dy      = (ymax - ymin)/(Ny)
        y       = [ymin + (i-1)*dy for i in 1:Ny]
    end

    #Create Differential Operators
    Dx  =   Diff_Operator_2D(1 ,1, [dx,dy] , [Nx,Ny])
    Dxx =   Diff_Operator_2D(2 ,1, [dx,dy] , [Nx,Ny])

    Dy  =   Diff_Operator_2D(1 ,2, [dx,dy] , [Nx,Ny])
    Dyy =   Diff_Operator_2D(2 ,2, [dx,dy] , [Nx,Ny])

    if(p.Boundaries == :periodic )
        Dx  =   Diff_Operator_2D_Periodic(1 ,1, [dx,dy] , [Nx,Ny])
        Dxx =   Diff_Operator_2D_Periodic(2 ,1, [dx,dy] , [Nx,Ny])

        Dy  =  Diff_Operator_2D_Periodic(1 ,2, [dx,dy] , [Nx,Ny])
        Dyy =  Diff_Operator_2D_Periodic(2 ,2, [dx,dy] , [Nx,Ny])
    end

    #Get max time iteration
    max_sim_iter , tmax = get_t_max(p)
    max_plot_iter = max_sim_iter

    if(tlimit <=  tmax)
        max_plot_iter = floor(Int , max_sim_iter * tlimit/tmax )
    else
        println("ATTENTION: max time input is larger than simulation time, using tmax")
    end
    #Mas user iteration

    #Potential Matrix
    Vpot = potential_matrix(p)

    # Total volume of space
    Volume = ((ymax-ymin)*(xmax-xmin))

    #Storage vectors
    E_vector = Float64[]
    t_vector = Float64[]

    t = 0
    E = 0

    #Loop to get energy
    for i in 1:every:max_plot_iter

        #Time and fields
        t, _ , _ , ψ , dψ = get_fields(p,i)

        #Volume element
        dV = dx*dy
        E = 0

        #Gradient term (possibly wrong)
         vx = Dx * reshape(ψ, Nx*Ny)
         vy = Dy * reshape(ψ, Nx*Ny)

        #Potential term
        v2 = Vpot * reshape(ψ, Nx*Ny)

        #Add all terms
        for i in 1:(Nx*Ny) #length(ψ)
            E+= vx[i]^2 + vy[i]^2 + dψ[i]^2 + (p.μ)^2*ψ[i]^2 + v2[i]^2
        end

        #Add to storage vectors  and normalize
        push!(t_vector, t)
        push!(E_vector, dV * E / Volume )
    end

    return t_vector , E_vector

end

function get_energy( fname::String, every::Int , tlimit::Float64)

    #Get Variables
    fid         = h5open( fname , "r")

    xmin = read_attribute(fid, "xmin")
    xmax = read_attribute(fid, "xmax")
    Nx   = read_attribute(fid, "xnodes")

    ymin = read_attribute(fid, "ymin")
    ymax = read_attribute(fid, "ymax")
    Ny   = read_attribute(fid, "ynodes")

    #Discretize space
    dx      = (xmax - xmin)/(Nx-1)
    x       = [xmin + (i-1)*dx for i in 1:Nx]

    dy      = (ymax - ymin)/(Ny-1)
    y       = [ymin + (i-1)*dy for i in 1:Ny]

    #Create Differential Operators
    Dx  =   Diff_Operator_2D(1 ,1, [dx,dy] , [Nx,Ny])
    Dxx =   Diff_Operator_2D(2 ,1, [dx,dy] , [Nx,Ny])

    Dy  =   Diff_Operator_2D(1 ,2, [dx,dy] , [Nx,Ny])
    Dyy =   Diff_Operator_2D(2 ,2, [dx,dy] , [Nx,Ny])

    #Get max time iteration
    max_sim_iter , tmax = get_t_max(fname)
    max_plot_iter = max_sim_iter

    if(tlimit <=  tmax)
        max_plot_iter = floor(Int , max_sim_iter * tlimit/tmax )
    else
        println("ATTENTION: max time input is larger than simulation time, using tmax")
    end
    #Mas user iteration

    #Potential Matrix
   # Vpot = potential_matrix(p)

    # Total volume of space
    Volume = ((ymax-ymin)*(xmax-xmin))

    #Storage vectors
    E_vector = Float64[]
    t_vector = Float64[]

    t = 0
    E = 0

    #Loop to get energy
    for i in 1:every:max_plot_iter

        #Time and fields
        t, _ , _ , ψ , dψ = get_fields(fname,i)

        #Volume element
        dV = dx*dy
        E = 0

        #Gradient term (possibly wrong)
         vx = Dx * reshape(ψ, Nx*Ny)
         vy = Dy * reshape(ψ, Nx*Ny)

        #Potential term
       # v2 = Vpot * reshape(ψ, Nx*Ny)

        #Add all terms
        for i in 1:(Nx*Ny) #length(ψ)
            E+= vx[i]^2 + vy[i]^2 + dψ[i]^2 # + (p.μ)^2*ψ[i]^2  # + v2[i]^2
        end

        #Add to storage vectors  and normalize
        push!(t_vector, t)
        push!(E_vector, dV * E / Volume )
    end

    return t_vector , E_vector

end


text:
module BrainFlowViz

    using BrainFlow
    using Makie, GLMakie, AbstractPlotting

    function init_scene(xs, ys;
        x_lim = (0, length(xs)),
        y_lim = (0, 1),
        theme = :light,
        color = :green,
    )
        nsamples, nchannels = size(ys)
        ys_n = []
        for n = 1:nchannels
            push!(ys_n, Node(view(ys, :, n)))
        end

        outer_padding = 30

        scene, layout = layoutscene(
            outer_padding,
            resolution = (1200, 700),
            backgroundcolor = RGBf0(0.99, 0.99, 0.99),
        )

        ax = Array{Axis, 1}(undef, nchannels)
        for n = 1:nchannels
            ax[n] = layout[n, 1] = Axis(scene, ylabel = "ch$n")
            lines!(ax[n], xs, ys_n[n], color = color, linewidth = 2)
            ax[n].yticklabelsvisible = false
            limits!(ax[n], x_lim[1], x_lim[2], y_lim[1], y_lim[2])
        end

        # TODO: use AbstractPlotting.Theme
        if theme == :dark
            scene.backgroundcolor[] = RGBf0(0.01, 0.01, 0.01)
            for n = 1:nchannels
                ax[n].backgroundcolor[] = :black
                ax[n].ylabelcolor[] = :white
                ax[n].xlabelcolor[] = :white
                ax[n].yticklabelcolor[] = :white
                ax[n].xticklabelcolor[] = :white
                ax[n].xtickcolor[] = :white
                ax[n].ytickcolor[] = :white
                ax[n].leftspinecolor[] = :white
                ax[n].rightspinecolor[] = :white
                ax[n].topspinecolor[] = :white
                ax[n].bottomspinecolor[] = :white
                ax[n].xgridcolor = :black
                ax[n].ygridcolor = :black
            end
        end
        #linkaxes!(ax...)

        # hide the x stuff
        for n=1:nchannels-1
            hidexdecorations!(ax[n], grid = false)
        end

        return scene, ys_n
    end

    # live plotting function
    function plot_data(data_func::Function, nsamples, nchannels; delay = 0.05, y_lim = :auto, kwargs...)

        xs = collect(1:nsamples)
        #ys = zeros(nsamples, nchannels)
        ys = data_func()

        if y_lim === :auto
            y_lim = (minimum(ys), maximum(ys))
            println("y_lim = $(y_lim)")
        end

        scene, ys_n = init_scene(xs, ys; y_lim = y_lim, kwargs...)

        display(scene)

        while(true)
            sleep(delay)

            # overwrite the original data
            ys = data_func()
            available_nsamples = size(ys, 1)

            # update x-axis without triggering a refresh
            empty!(xs)
            append!(xs, collect(1:available_nsamples))

            # update the scene
            # is this the most efficient? Or could we refresh 4 axes in a scene at the same time
            for n = 1:nchannels
                ys_n[n][] = view(ys, :, n)
            end

            # optional if you want to zoom in + update limits:
            # AbstractPlotting.update_limits!(scene)
            # AbstractPlotting.update!(scene)

            # yield() # allow another task to run, e.g. brainflow data gathering? or use yieldto()
        end
    end

end # module