text:
module SCoBA

using StaticArrays
using Hungarian
using LinearAlgebra
using Distributions
using Parameters
using DataStructures
using Random
using ArgParse
using Distances
using NearestNeighbors
using GeometryTypes
using TOML

using IterTools
using Reexport

# POMDP stuff
using POMDPs
using MCTS
using POMDPModelTools

# Include submodule files
include("solver/SCoBASolver.jl")
include("domains/SCoBADomains.jl")

@reexport using SCoBA.SCoBASolver
@reexport using SCoBA.SCoBADomains

end # module


text:
scale =2; shearLevel = 2; inputSize = (25,25,4,5,3); useGpu = false; σm =abs; dType=Float32
cpuShears = 3; singleRes = 3; shears = 3
# standard input size
@testset "shearing constructors tiny" begin
    inputSizes = [(25,25,1,2), (25, 25, 4, 5, 3)]
    scalesShears = [(1,1),(2,1), (2,2),(4,1), (4,2), (4,3), (4,4)]#, (8,1), (8,2), (8,3), (8,4), (8,6), (8,8)]
    shearLevels = [1,2,4,5]
    if CUDA.functional()
        useGpus = [true, false]
    else
        useGpus = [false]
    end
    dType = Float32
    σs =[identity, abs, relu]
    for inputSize in inputSizes, (scale, shearLevel) in scalesShears, useGpu in useGpus, σm in σs
        @testset "inputSize=$inputSize, (scale, shearLevel) =$((scale,shearLevel)),useGpu=$useGpu, σm=$(σm)" begin
            shears=3
            with_logger(ConsoleLogger(stderr,Logging.Error)) do
                #global shears
                shears = shearingLayer(inputSize, scale=scale, shearLevel=shearLevel,
                                       dType = dType, σ=σm)
            end
            @test size(shears.fftPlan) ==
                inputSize .+ (2 .* shears.bc.padBy...,
                              fill(0, length(inputSize)-2)...)
            @test shears.σ == σm
            @test shears.bias == nothing
            @test ndims(shears.weight)==3
            @test outType(shears)<:Real
            if useGpu
                init = gpu(randn(dType, inputSize));
                shears = shears |> gpu
            else
                init = randn(dType, inputSize);
            end

            fullResult = shears(init);
            @test size(fullResult) == (inputSize[1:2]...,size(shears.weight,3),
                                       inputSize[3:end]...)
            if σm!=relu
                # are there only a very small number of zeros?
                @test count(abs.(fullResult) .== 0) < .001 *prod(size(fullResult))
            end
            # compare with the result from Shearlab itself
            effectiveShears = ceil.(Int, (1:shearLevel)/2)
            cpuShears=3; res = 2
            with_logger(ConsoleLogger(stderr,Logging.Error)) do
                #global cpuShears, res
                cpuShears = Shearlab.getshearletsystem2D(inputSize[1:2]...,
                                                         scale,
                                                         effectiveShears,
                                                         typeBecomes=Float32,
                                                         padded=true)
                res = Shearlab.sheardec2D(init, cpuShears)
            end
            @test σm.(res) ≈ cpu(fullResult)
            # just the averaging filter
            with_logger(ConsoleLogger(stderr,Logging.Error)) do
                #global shears
                shears = shearingLayer(inputSize, scale=scale, shearLevel=shearLevel,
                                       dType = dType, σ=σm, averagingLayer=true)
            end
            @test size(shears.weight, 3) ==1
            if useGpu
                shears = shears |> gpu
            end
            singleRes = 3
            with_logger(ConsoleLogger(stderr,Logging.Error)) do
                #global singleRes
                singleRes = shears(init)
            end
            ax = axes(fullResult)
            @test fullResult[:,:,end:end, ax[4:end]...] ≈ singleRes
        end
    end
end

ext: jl


text:
function crra_utility(;
    β = 0.9,
    σ = 2.0,
    γ = 2.0,
    Π = 0.5 * ones(2, 2),
    G = [0.1, 0.2],
    Θ = ones(Float64, 2),
    transfers = false
    )
    function U(c, n)
        if σ == 1.0
            U = log(c)
        else
            U = (c.^(1.0 - σ) - 1.0) / (1.0 - σ)
        end
        return U - n.^(1 + γ) / (1 + γ)
    end
    # Derivatives of utility function
    Uc(c,n) =  c.^(-σ)
    Ucc(c,n) = -σ * c.^(-σ - 1.0)
    Un(c,n) = -n.^γ
    Unn(c,n) = -γ * n.^(γ - 1.0)
    n_less_than_one = false
    return Model(β, Π, G, Θ, transfers,
                 U, Uc, Ucc, Un, Unn, n_less_than_one)
end


text:
# Create a hash of the log for an untainted webapp container
# It is important that the log is not stateful, i.e. that the
# log output is always the same after the container is fully initialized.

using SHA

file = "swarm_vdma"    # hashfile
ID = "4f880d0e7d68"    # ID of an untainted fully executed app
hashlog(ID) = (read(`docker logs $(ID)`) |> sha256 |> bytes2hex)

open(f -> write(f, hashlog(ID)), file, "w")


text:
function solve_disc(model::State_Space_Form{T}, obj::State_Space_Objective{T}, tol::T, maxiters::S) where {T <: AbstractFloat, S <: Int}

  nx    = copy(model.nx)
  ny    = copy(model.ny)
  beta  = copy(obj.beta)
  sigma = copy(model.sigma)

  a = copy(model.a)
  b = copy(model.b)
  c = copy(model.c)

  q = copy(obj.q)
  u = copy(obj.u)
  r = copy(obj.r)

  n = nx+ny

  a11 = a[1:nx,1:nx]
  a12 = a[1:nx,(nx+1):n]
  a21 = a[(nx+1):n,1:nx]
  a22 = a[(nx+1):n,(nx+1):n]

  b1 = b[1:nx,:]
  b2 = b[(nx+1):n,:]

  q11 = q[1:nx,1:nx]
  q12 = q[1:nx,(nx+1):n]
  q21 = q[(nx+1):n,1:nx]
  q22 = q[(nx+1):n,(nx+1):n]

  u1 = u[1:nx,:]
  u2 = u[(nx+1):n,:]

  j = (I-a22)\a21   # I = eye(ny)
  k = (I-a22)\b2    # I = eye(ny)
  f = (k'q22*k+k'u2+u2'k+r)\(k'q21+k'q22*j+u1'+u2'j)
  h = j-k*f
  v = q11+2*q12*h+h'q22*h-2*u1*f-2*h'u2*f+f'r*f

  len = Inf
  iters = 0
  while len > tol && iters < maxiters

    j    = (a22-h*a12)\(h*a11-a21)
    k    = (a22-h*a12)\(h*b1-b2)
    abar = a11+a12*j
    bbar = b1+a12*k
    qbar = q11+j'q21+q12*j+j'q22*j
    ubar = q12*k+j'q22*k+u1+j'u2
    rbar = r+k'q22*k+u2'k+k'u2
    fn   = (rbar+beta*bbar'v*bbar)\(ubar'+beta*bbar'v*abar)
    vn   = beta*(abar-bbar*fn)'v*(abar-bbar*fn)+qbar-fn'ubar'-ubar*fn+fn'rbar*fn
    hn   = j-k*fn

    len = maximum(maximum(abs,vn-v), maximum(abs,vn-v), maximum(abs,fn-f))
    h = hn
    v = vn
    f = fn

    iters += 1

  end

  if iters <= maxiters
    retcode = true
  else
    retcode = false
  end

  f = -f
  p = a11+a12*h+b1*f
  h = [eye(nx); h]
  k = c[1:nx,:]

  soln = State_Space_Soln(p,k,h,f,v,sigma,retcode)

  return soln

end

function solve_disc(model::Generalized_State_Space_Form{T}, obj::State_Space_Objective{T}, tol::T, maxiters::S) where {T <: AbstractFloat, S <: Int}

  nx    = copy(model.nx)
  ny    = copy(model.ny)
  beta  = copy(obj.beta)
  sigma = copy(model.sigma)

  a0 = copy(model.a0)
  a  = copy(model.a)
  b  = copy(model.b)
  c  = copy(model.c)

  q = copy(obj.q)
  u = copy(obj.u)
  r = copy(obj.r)

  n = nx+ny

  a11 = a[1:nx,1:nx]
  a12 = a[1:nx,(nx+1):n]
  a21 = a[(nx+1):n,1:nx]
  a22 = a[(nx+1):n,(nx+1):n]

  b1 = b[1:nx,:]
  b2 = b[(nx+1):n,:]

  q11 = q[1:nx,1:nx]
  q12 = q[1:nx,(nx+1):n]
  q21 = q[(nx+1):n,1:nx]
  q22 = q[(nx+1):n,(nx+1):n]

  u1 = u[1:nx,:]
  u2 = u[(nx+1):n,:]

  j = (a0-a22)\a21
  k = (a0-a22)\b2
  f = (k'q22*k+k'u2+u2'k+r)\(k'q21+k'q22*j+u1'+u2'j)
  h = j-k*f
  v = q11+2*q12*h+h'q22*h-2*u1*f-2*h'u2*f+f'r*f

  len = Inf
  iters = 0
  while len > tol && iters < maxiters

    j    = (a22-a0*h*a12)\(a0*h*a11-a21)
    k    = (a22-a0*h*a12)\(a0*h*b1-b2)
    abar = a11+a12*j
    bbar = b1+a12*k
    qbar = q11+j'q21+q12*j+j'q22*j
    ubar = q12*k+j'q22*k+u1+j'u2
    rbar = r+k'q22*k+u2'k+k'u2
    fn   = (rbar+beta*bbar'v*bbar)\(ubar'+beta*bbar'v*abar)
    vn   = beta*(abar-bbar*fn)'v*(abar-bbar*fn)+qbar-fn'ubar'-ubar*fn+fn'rbar*fn
    hn   = j-k*fn

    len = maximum(maximum(abs,vn-v), maximum(abs,vn-v), maximum(abs,fn-f))
    h = hn
    v = vn
    f = fn

    iters += 1

  end

  if iters <= maxiters
    retcode = true
  else
    retcode = false
  end

  f = -f
  p = a11+a12*h+b1*f
  h = [eye(nx); h]
  k = c[1:nx,:]

  soln = State_Space_Soln(p,k,h,f,v,sigma,retcode)

  return soln

end

function solve_disc(model::Structural_Form{T}, obj::Structural_Objective{T}, tol::T, maxiters::S) where {T <: AbstractFloat, S <: Int}

  a0    = copy(model.a0)
  a1    = copy(model.a1)
  a2    = copy(model.a2)
  a3    = copy(model.a3)
  a5    = copy(model.a5)
  sigma = copy(model.sigma)

  q = copy(obj.q)
  r = copy(obj.r)

  beta = copy(obj.beta)

  h1 = a0\a1
  h2 = a0\a5
  f1 = zeros(size(a3'))
  f2 = zeros(size(a3,2),size(a5,2))
  v = q

  a3_bar = Array(T,size(a3))
  a5_bar = Array(T,size(a5))

  len = Inf
  iters = 0
  while len > tol && iters < maxiters

    a1_bar = (a0-a2*h1)\a1
    a3_bar = (a0-a2*h1)\a3
    a5_bar = (a0-a2*h1)\a5

    vn = q+beta*(a1_bar+a3_bar*f1)'v*(a1_bar+a3_bar*f1)

    f1n = -(r+a3_bar'vn*a3_bar)\(a3_bar'vn*a1_bar)
    h1n = a1_bar+a3_bar*f1

    len = maximum(maximum(abs,vn-v), maximum(abs,h1n-h1), maximum(abs,f1n-f1))
    h1 = h1n
    f1 = f1n
    v  = vn

    iters += 1

  end

  if iters <= maxiters
    retcode = true
  else
    retcode = false
  end

  f2 = -(r+a3_bar'v*a3_bar)\(a3_bar'v*a5_bar)
  h2 = a5_bar+a3_bar*f2

  h = [[h1 zeros(size(a3))];
       [f1 zeros(size(r))]]

  g = [h2; f2]

  soln = Structural_Soln(h,g,v,sigma,retcode)

  return soln

end

function solve_disc(model::Generalized_Structural_Form{T}, obj::Structural_Objective{T}, tol::T, maxiters::S) where {T <: AbstractFloat, S <: Int}

  a0    = copy(model.a0)
  a1    = copy(model.a1)
  a2    = copy(model.a2)
  a3    = copy(model.a3)
  a4    = copy(model.a4)
  a5    = copy(model.a5)
  sigma = copy(model.sigma)

  q = copy(obj.q)
  r = copy(obj.r)

  beta = copy(obj.beta)

  h1 = a0\a1
  f1 = zeros(size(a3'))
  di = Array(T,size(a0))
  v  = q

  len = Inf
  iters = 0
  while len > tol && iters < maxiters

    (v,retcode_doubling) = doubling(beta*h1',h1,q+beta*f1'r*f1,tol,maxiters)
    di  = inv(a0-a2*h1-a4*f1)
    f1n = -(r+a3'di'v*di*a3)\(a3'di'v*di*a1)
    h1  = di*(a1+a3*f1n)

    len = maximum(abs,f1n-f1)
    f1 = f1n

    iters += 1

  end

  if iters <= maxiters
    retcode = true
  else
    retcode = false
  end

  f2 = -(r+a3'di'v*di*a3)\(a3'di'v*di*a5)
  h2 = (a0-a2*h1-a4*f1)\(a5+a3*f2)

  h = [h1 zeros(size(a3));
       f1 zeros(size(r))]

  g = [h2; f2]

  soln = Structural_Soln(h,g,v,sigma,retcode)

end


text:
@testset "Get Class Name Function Tests" begin

  @test isdefined(Julz, :get_class_name) == true

  @test_throws MethodError Julz.get_class_name(1)

  @test Julz.get_class_name("first") == "First"

  @test Julz.get_class_name("first class") == "FirstClass"

  @test Julz.get_class_name("first_class") == "FirstClass"

  @test Julz.get_class_name("first-class") == "First-class"

end


text:
import Symbolics: derivative

function linearize(::Star, ctx, F, X)
    top, n, (τ⁻, τ⁺) = getproperty.(Ref(ctx), (:top, :n, :τ))

    T = promote_type(typeof(F), typeof(X))
    diags = Dict{Int,T}()

    for i in eachindex(top)
        p = prod(n[1:i-1])
        ∂ = derivative.(F, τ⁻[i] * X)
        if isperiodic(top[i])
            D = map(CartesianIndices(n), reshape(∂, n...)) do ind, el
                ind[i] == 1 ? zero(el) : el
            end
            diags[-p] = D[begin+p:end]
            D = map(CartesianIndices(n), reshape(∂, n...)) do ind, el
                ind[i] != 1 ? zero(el) : el
            end
            q = (n[i]-1) * p
            diags[q] = D[begin:end-q]
        else
            diags[-p] = ∂[begin+p:end]
        end
    end

    diags[0] = derivative.(F, X)

    for i in eachindex(top)
        p = prod(n[1:i-1])
        ∂ = derivative.(F, τ⁺[i] * X)
        if isperiodic(top[i])
            D = map(CartesianIndices(n), reshape(∂, n...)) do ind, el
                ind[i] == n[i] ? zero(el) : el
            end
            diags[+p] = D[begin:end-p]
            D = map(CartesianIndices(n), reshape(∂, n...)) do ind, el
                ind[i] != n[i] ? zero(el) : el
            end
            q = (n[i]-1) * p
            diags[-q] = D[begin+q:end]
        else
            diags[+p] = ∂[begin:end-p]
        end
    end

#    spdiagm(diags...)
    diags
end


text:
@testset "783.minimum-distance-between-bst-nodes.jl" begin
    t = TreeNode(4,
           TreeNode(2,
               TreeNode(1),
               TreeNode(3)
           ),
           TreeNode(6)
       )
    @test min_diff_in_bst(t) == 1
end


text:
"""
    Object

Overarching type for objects displayable on `Axes`.
"""
abstract type Object end

"""
    Object2D

Overarching type for objects displayable on `Axes2D`.
"""
abstract type Object2D <: Object end

"""
    Text2D <: Object2D

Place text somewhere relative to current axes.
"""
mutable struct Text2D <: Object2D
    anchor::T2F            # where the text is located (x,y) position, relative to the axes
    text  ::String         # what to write
    #
    position ::String      # position with respect to anchor (by default centered on anchor)
    textstyle::TextStyle
end
Text2D(a, t) = default(Text2D, a, t)

"""
    StraightLine2D <: Object2D

Place either a vertical or horizontal straightline at given anchor (from axis
to axis).
"""
mutable struct StraightLine2D <: Object2D
    anchor::F64
    horiz ::Bool
    #
    linestyle::LineStyle
end
StraightLine2D(a, h) = default(StraightLine2D, a, h)

"""
    Box2D <: Object2D

Place a 2D filled box.
"""
mutable struct Box2D <: Object2D
    anchor::T2F  # where the box is (in graph units)
    size  ::T2F  # width, hei (in graph units)
    #
    position ::String             # position of the anchor with respect to the box
    nobox    ::Bool               # show an edge or not
    linestyle::LineStyle          # style of the box edge
    fillstyle::Option{FillStyle}  # if none --> transparent box
    # text::String
    # textStyle::TextStyle = TextStyle()
end
Box2D(a, s) = default(Box2D, a, s)

"""
    Colorbar <: Object2D

Add a colorbar.
"""
mutable struct Colorbar <: Object2D
    zmin::F64
    zmax::F64
    cmap::Vector{Color}
    #
    ticks   ::Ticks        # constructed
    size    ::Option{T2F}  # (width, height)
    pixels  ::Int          # resolution for the color bar
    nobox   ::Bool         #
    position::String       # left, right, bottom, top
    offset  ::T2F          # h, v distance from closest corner of the graph
end
Colorbar(zmin, zmax, cmap) = default(Colorbar,
    zmin, zmax, cmap, Ticks(collect(range(zmin, stop=zmax, length=5))[2:end-1])
)


text:
using MakieRipserer
using Aqua

# Don't test for ambiguities. There are too many in imported packages.
Aqua.test_all(MakieRipserer; ambiguities=false)


text:
module brainflow

include("brainflow_logs.jl")
include("exit_codes.jl")
include("board_shim.jl")
include("data_filter.jl")
include("ml_model.jl")

end # module

max_stars_repo_name: Cognixion-inc/brainflow

max_issues_repo_name: Cognixion-inc/brainflow

max_forks_repo_name: Cognixion-inc/brainflow


text:
@reexport module Blaze

#export BlazePolynomial, BlazePolynomialList, BlazeMap1D, BlazeMap2D

# Default functions for AbstractBlazeTypes
min_log_lambda(blaze::AbstractBlazeType) = blaze.log_lambda_lo
max_log_lambda(blaze::AbstractBlazeType) = blaze.log_lambda_hi
order(blaze::AbstractBlazeType) = blaze.order
pixel_range(blaze::AbstractBlazeType) = blaze.pixel_range
export order, min_log_lambda, max_log_lambda, pixel_range

#= immutable
function update_log_lambda_min_max!(blaze::AbstractBlazeType, data::DataT) where DataT <: AbstractEchelleObservationOrSet
   r = EchelleTools.make_model_wavelength_range(data)
   blaze.log_lambda_lo = r[1]
   blaze.log_lambda_hi = r[end]
end
=#

#@reexport module BlazePolynomial

using Polynomials     # For blaze function

""" Type for storing a polynomial model for one order's blaze function. """
struct BlazePolynomial{T1,T2} <: AbstractBlaze1DType where {T1<:Integer, T2<:Real}
   order::T1
   log_lambda_lo::T2
   log_lambda_hi::T2
   pixel_range::UnitRange{T1}
   model::Poly{T2}
end
export BlazePolynomial

struct BlazePolynomialList{T1,T2} <: AbstractBlaze2DType where {T1<:Integer, T2<:Real}
   blaze_list::Vector{BlazePolynomial{T1,T2}}
end
export BlazePolynomialList

function fit_blaze_1d(x::AbstractRange{T1}, flux::AbstractArray{T2,1}; degree::Int = default_blaze_degree
      ) where {T1<:Real, T2<:Real}
   polyfit(x,flux,degree) # replaced x.-mean(x)
end
export fit_blaze_1d

function make_flat_blaze_polynomial(;num_pixels_per_order::Integer = max_num_pixels_per_order,
      pixel_range::AbstractRange = 1:num_pixels_per_order,
      order::Integer = 0, degree::Integer = 0, log_lambda_lo::Real = 0.0, log_lambda_hi::Real = 0.0)
      @assert 0 <= degree <= max_degree_blaze
      coeff = vcat([1.0],zeros(degree))
      return BlazePolynomial(order, log_lambda_lo, log_lambda_hi, pixel_range, Poly(coeff) )
end

function make_flat_blaze_polynomial(data::EOST, order::Integer;
      pixel_range::AbstractRange = 1:num_pixels_per_order(data), degree::Integer = 0
      ) where { EOST <: Union{AbstractEchelleObservationSet,AbstractEchelleObservation} }
   @assert 1 <= order <= max_num_orders
   @assert 0 <= degree <= max_degree_blaze
   pixel_lo = max(1,pixel_range[1] - 1)
   pixel_hi = min(num_pixels_per_order(data), pixel_range[end]+1)
   log_lambda_lo = log(minimum(data.lambda[pixel_lo,order,:]))
   log_lambda_hi = log(maximum(data.lambda[pixel_hi,order,:]))
   make_flat_blaze_polynomial(order=order, log_lambda_lo=log_lambda_lo, log_lambda_hi=log_lambda_hi,
                  pixel_range=pixel_range, degree=degree)
end
export make_flat_blaze_polynomial

function make_flat_blaze_polynomial_2d(data::EOST;
      pixel_range::AbstractRange = 1:num_pixels_per_order(data), degree::Integer = 0,
      ) where { EOST <: Union{AbstractEchelleObservationSet,AbstractEchelleObservation} }
   map(order-> make_flat_blaze_polynomial(pixel_range=pixel_range, order=order, degree=degree), 1:num_orders(data) )
end
export make_flat_blaze_polynomial_2d

function make_blaze(data::EOST, order::Integer;
      pixel_range::AbstractRange = 1:size(data.lambda,1),
      pixel_buffer_lo::Integer = default_pixel_buffer_lo,
      pixel_buffer_hi::Integer = default_pixel_buffer_hi,
      degree::Integer = default_blaze_degree,
      verbose::Integer = 0
      ) where { EOST <: AbstractEchelleObservationSet }
      @assert 0 <= degree <= max_degree_blaze
      @assert 1 <= order <= num_orders(data)
      max_itterations_fit_blaze = 20
      # TODO: Insert code to truncate pixel_range to enforce positive blaze
      # blaze = fit_blaze_1d(pixel_range, vec(mean(data.flux[pixel_range,order,:],dims=2)), degree = degree)
      #println("# Order = ", order, " pixel_range = ", pixel_range)
      new_pixel_range = pixel_range
      last_pixel_range = 0:0
      local blaze
      i = 0
      while !((2*abs(new_pixel_range[1] - last_pixel_range[1]) <= pixel_buffer_lo) && (2*abs(new_pixel_range[end] - last_pixel_range[end]) <= pixel_buffer_hi )) &&
                  (i<max_itterations_fit_blaze)
            last_pixel_range = new_pixel_range
            blaze = fit_blaze_1d(new_pixel_range, vec(mean(data.flux[new_pixel_range,order,:],dims=2)), degree = degree)
            pixel_lo = max(pixel_range[1],  findfirst(x->x>0,polyval(blaze,pixel_range)) ) +pixel_buffer_lo
            pixel_hi = min(pixel_range[end],findlast( x->x>0,polyval(blaze,pixel_range)) ) -pixel_buffer_hi
            if verbose > 1
                  println("# order: ", order, " it=", i," range: ", pixel_lo, ":", pixel_hi, " blaze = ", polyval(blaze,max(pixel_lo-1,1)), " ", polyval(blaze,pixel_lo), " ... ", polyval(blaze, pixel_hi) )
            end
            new_pixel_range = pixel_lo:pixel_hi
            i += 1
      end
      if verbose > -1 && !((2*abs(new_pixel_range[1] - last_pixel_range[1]) <= pixel_buffer_lo) && (2*abs(new_pixel_range[end] - last_pixel_range[end]) <= pixel_buffer_hi ))
            println("# Warning: Reached maximum number of itterations fitting blaze for order ", order)
            println("# last_pixel_range = ", last_pixel_range)
            println("# new_pixel_range = ", new_pixel_range)
      end
      # blaze = blaze / mean(polyval(blaze,new_pixel_range))
      pixel_lo = max(1,new_pixel_range[1] - 1)
      pixel_hi = min(num_pixels_per_order(data), new_pixel_range[end]+1)
      log_lambda_lo = log(minimum(data.lambda[pixel_lo,order,:]))
      log_lambda_hi = log(maximum(data.lambda[pixel_hi,order,:]))
      return BlazePolynomial(order, log_lambda_lo, log_lambda_hi, new_pixel_range, blaze )
end
export make_blaze

function make_blaze_polynomial(data::EOST;
      pixel_range::AbstractRange = 1:size(data.lambda,1), degree::Integer = default_blaze_degree
      ) where { EOST <: AbstractEchelleObservationSet}
      map(order-> make_blaze_polynomial(data, order, pixel_range=pixel_range, degree=degree), 1:num_orders(data) )
end
export make_blaze_polynomial

# TODO: Change calc_blaze to eval_blaze
function calc_blaze(b::BlazePolynomial, pixel::Real)
      @assert b.pixel_range[1] <= pixel <= b.pixel_range[end]
      polyval(b.model, pixel)
end
export calc_blaze

function calc_blaze(b::BlazePolynomial, pixel_range::AbstractRange{T2}) where {T2<:Real}
      @assert b.pixel_range[1]-0.5 <= pixel_range[1] <= b.pixel_range[end]+0.5
      @assert b.pixel_range[1]-0.5 <= pixel_range[end] <= b.pixel_range[end]+0.5
      blaze = polyval(b.model, pixel_range)
end

function calc_blaze(b::BlazePolynomial) where {T2<:Real}
      blaze = polyval(b.model, b.pixel_range)
end

function calc_blaze(blaze_list::BlazeLT, data::EOT;
      output::AbstractArray{T2,2} = zeros(eltype(data.flux),size(data.lambda))
      ) where { T2<:Real,
      BlazeLT <: AbstractBlaze2DType,
      #BlazeT <: AbstractBlazeType, BlazeLT <: AbstractArray{BlazeT,1},
      EOT <: AbstractEchelleObservation  }
      @assert length(blaze_list) == num_orders(data)
      for blaze in blaze_list
            output[blaze.pixel_range,blaze.order] .= calc_blaze(blaze)
      end
      return output
end

function calc_blaze(blaze_list::BlazeLT, data::EOST;
      output::AbstractArray{T2,3} = zeros(eltype(data.flux),size(data.lambda))
      ) where { T2<:Real,
      BlazeLT <: AbstractBlaze2DType,
      #BlazeT <: AbstractBlazeType, BlazeLT <: AbstractArray{BlazeT,1},
      EOST <: AbstractEchelleObservationSet  }
      @assert length(blaze_list) == num_orders(data)
      for blaze in blaze_list
            output[blaze.pixel_range,blaze.order,:] .= calc_blaze(blaze)
      end
      return output
end

function calc_normalizations(blaze_list::BlazeLT, data::EOST;
      output::OutputT = zeros(eltype(data.flux),num_obs(data))
      ) where { T2<:Real,
      BlazeLT <: AbstractBlaze2DType,
      # BlazeT <: AbstractBlazeType, BlazeLT <: AbstractArray{BlazeT,1},
      EOST <: AbstractEchelleObservationSet, OutputT <: AbstractArray{T2,1}  }
   @assert length(blaze_list) == num_orders(data)
   @assert length(output) == num_obs(data)
   output .= map(obsid->
      mean(map(order->
            mean(data.flux[blaze_list[order].pixel_range,order,obsid] ./
                  calc_blaze(blaze_list[order]) )
            , 1:num_orders(data) ))
      , 1:num_obs(data))
end
export calc_normalizations

function renormalize(blaze::BlazePolynomial, norm::Real)
      return BlazePolynomial(blaze.order, blaze.log_lambda_lo, blaze.log_lambda_hi,
                  blaze.pixel_range, blaze.model*norm )
end
export renormalize

#end # module BlazePolynomial

#@reexport module BlazeMap

struct BlazeMap1D{T1,T2,T3} <: AbstractBlaze1DType where {T1<:Integer, T2<:Real, T3<:Real}
   order::T1
   log_lambda_lo::T2
   log_lambda_hi::T2
   pixel_range::UnitRange{T1}
   model::AbstractArray{T3,1}
end
export BlazeMap1D

function make_flat_blaze_map_1d(;num_pixels_per_order::Integer = max_num_pixels_per_order,
      pixel_range::AbstractRange = 1:num_pixels_per_order,
      order::Integer = 0, log_lambda_lo::Real = 0.0, log_lambda_hi::Real = 0.0)
      return BlazeMap1D(order, log_lambda_lo, log_lambda_hi, pixel_range, ones(length(pixel_range)) )
end
export make_flat_blaze_map_1d

# TODO: Change calc_blaze to eval_blaze
function calc_blaze(b::BT, pixel::Real) where {BT<:AbstractBlaze1DType }
      @assert b.pixel_range[1]-1 <= pixel <= b.pixel_range[end] +1
      pixel_lo = maximum(1, floor(Int64,pixel))
      pixel_hi = minimum(pixel_lo+1, b.pixel_range[end])
      w_hi = (pixel-pixel_lo)
      w_lo = 1-w_hi
      w_lo*b.model[pixel_lo] + w_hi*b.model[pixel_hi]
end

function calc_blaze(b::BT, pixel::Integer) where {BT<:AbstractBlaze1DType }
      @assert b.pixel_range[1] <= pixel <= b.pixel_range[end]
      b.model[pixel]
end

function calc_blaze(b::BT, pixel_range::AbstractRange{T2}) where {BT<:AbstractBlaze1DType, T2<:Real}
      map(p->calc_blaze(b,p), pixel_range)
end

function calc_blaze(b::BT, pixel_range::AbstractRange{T2}) where {BT<:AbstractBlaze1DType, T2<:Integer}
      @assert b.pixel_range[1]-0.5 <= pixel_range[1] <= b.pixel_range[end]+0.5
      @assert b.pixel_range[1]-0.5 <= pixel_range[end] <= b.pixel_range[end]+0.5
      b.model[pixel_range]
end

function calc_blaze(b::BT) where {BT<:AbstractBlaze1DType, T2<:Real}
      calc_blaze(b, b.pixel_range)
end
export calc_blaze

struct BlazeMap2D{T1,T2,T3} <: AbstractBlaze2DType where {T1<:Integer, T2<:Real, T3<:Real}
   blaze_list::Vector{BlazeMap1D{T1,T2,T3}}
end
export BlazeMap2D

function make_flat_blaze_map_2d(;num_pixels_per_order::Integer = max_num_pixels_per_order,
      num_orders::T1 = max_num_orders,
      pixel_range::AbstractArray{AbstractRange} = fill(1:num_pixels_per_order, num_orders),
      orders::AbstractArray{T1,1} = collect(1:num_orders),
      log_lambda_lo::AbstractArray{T2,1} = zeros(num_orders), log_lambda_hi::AbstractArray{T2,1} = zeros(num_orders)
      ) where { T1<:Integer, T2<:Real }
      @assert length(pixel_range) == length(orders) == length(log_lambda_lo) == length(log_lambda_hi) >= 1
      blaze_map = ones(length(pixel_range[1]))
      return BlazeMap2D( [
            BlazeMap1D(orders[i], log_lambda_lo[i], log_lambda_hi[i], pixel_range[i], blaze_map )
                  for i in 1:length(orders) ] )
end
export make_flat_blaze_map_2d

function covert_blaze_map_into_blaze_polynomial(b::BMT; degree = default_blaze_degree
   ) where { BMT <: BlazeMap1D }
   poly_model = fit_blaze_1d(b.pixel_range, b.model[b.pixel_range], degree=degree)
   BlazePolynomial(b.order,b.lambda_lo,b.lambda_hi,b.pixel_range,poly_model)
end
export covert_blaze_map_into_blaze_polynomial

function covert_blaze_map_into_blaze_polynomial(blaze_list::BMT; degree = default_blaze_degree
      ) where { BMT <: BlazeMap2D }
   BlazeMap2D( [ covert_blaze_map_into_blaze_polynomial(blaze_list[i],degree=degree) for i in 1:length(blaze_list) ] )
end
export covert_blaze_map_into_blaze_polynomial

#end # module BlazeMap

end # module Blaze


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.DLR_Visualization (Modia3D/renderer/DLR_Visualization/_module.jl)
#

mutable struct SimVis_Renderer
    velements::Vector{Composition.Object3D}  # Objects to be visualized
    ids::Vector{Ptr{Nothing}}                # ids[i] is the SimVis id of velements[i]
    isInitialized::Bool                      # = true, if SimVis is initialized (SimVis_init was called)

    host::String
    port::Int
    sync::Bool

    directory::String            # Directory Visualization/Extras
    dll_name::String             # Absolute path of SimVis DLL as string
    isProfessionalEdition::Bool  # = true, if DLL is professional SimVis edition otherwise community edition
    isNoRenderer::Bool           # = true, if SimVis DLL is not available and the NoRenderer renderer is used
    dll_handle::Ptr{Nothing}     # Handle for SimVis DLL

    # SimVis functions
    init::Ptr{Nothing}
    shutdown::Ptr{Nothing}
    getObjectID::Ptr{Nothing}
    freeObjectID::Ptr{Nothing}
    setTime::Ptr{Nothing}
    setBaseObject::Ptr{Nothing}
    setFileObject::Ptr{Nothing}
    setTextObject::Union{Ptr{Nothing}, Nothing}

    function SimVis_Renderer(info; host="127.0.0.1",port=11000,sync=false, professionalEdition=true)
        (directory, dll_name, isProfessionalEdition, isNoRenderer) = info

        dll_handle    = Libdl.dlopen(dll_name)
        init          = Libdl.dlsym(dll_handle, :SimVis_init)
        shutdown      = Libdl.dlsym(dll_handle, :SimVis_shutdown)
        getObjectID   = Libdl.dlsym(dll_handle, :SimVis_getObjectID)
        freeObjectID  = Libdl.dlsym(dll_handle, :SimVis_freeObjectID)
        setTime       = Libdl.dlsym(dll_handle, :SimVis_setTime)
        setBaseObject = Libdl.dlsym(dll_handle, :SimVis_setBaseObject)
        setFileObject = Libdl.dlsym(dll_handle, :SimVis_setFileObject)
        if professionalEdition
            setTextObject = Libdl.dlsym(dll_handle, :SimVis_setTextObject)
        else
            setTextObject = nothing
        end

        new(Any[], Ptr{Nothing}[], false, host, port, sync,
            directory, dll_name, isProfessionalEdition, isNoRenderer,
            dll_handle, init, shutdown, getObjectID, freeObjectID, setTime,
            setBaseObject, setFileObject, setTextObject)
    end
end

struct ProfessionalEdition <: Modia3D.AbstractDLR_VisualizationRenderer
    simVis::SimVis_Renderer

    ProfessionalEdition(info; host="127.0.0.1",port=11000,sync=false) =
        new( SimVis_Renderer(info, host=host,port=port,sync=sync) )
end

struct CommunityEdition <: Modia3D.AbstractDLR_VisualizationRenderer
    simVis::SimVis_Renderer

    CommunityEdition(info; host="127.0.0.1",port=11000,sync=false) =
        new( SimVis_Renderer(info, host=host,port=port,sync=sync,professionalEdition=false) )
end

const UndefinedDirectory = "???"

function getSimVisInfo()
    dll_name = ""
    isProfessionalEdition = false

    # Get directory of SimVis2.exe
    if haskey(ENV, "DLR_VISUALIZATION")
        directory = ENV["DLR_VISUALIZATION"]
    else
        directory = UndefinedDirectory
        println("\n... Info: ENV[\"DLR_VISUALIZATION\"] is not defined.\n" *
                "          Therefore, no online renderer is used in Modia3D.\n")
        return (directory,"???",false,true)
    end

    # Check for 64 bit
    if Base.Sys.WORD_SIZE != 64
        @warn "DLR Visualization library only supported for 64-bit system but not for $(Base.Sys.WORD_SIZE) bit.\n" *
              "\nNo Renderer is used in Modia3D (so, animation is switched off)."
        return (directory,"???",false,true)
    end

    # Check whether professional or community edition or on windows or on linux
    if ISWINDOWS()
        dll_name1 = joinpath(directory, "windows", "SimVisInterface_ProfessionalEdition.dll")
        if isfile( dll_name1 )
            dll_name = dll_name1
            isProfessionalEdition = true
        else
            dll_name2 = joinpath(directory, "windows", "SimVisInterface_CommunityEdition.dll")
            if isfile( dll_name2 )
                dll_name = dll_name2
                isProfessionalEdition = false
            else
                @warn "\nModia3D: DLL of DLR-Visualization library not found. Neither of these files\n" *
                      "   $dll_name1, \n" *
                      "   $dll_name2, \n" *
                      "exists. Check whether ENV[\"DLR_VISUALIZATION\"] is correct." *
                      "\nNo Renderer is used in Modia3D (so, animation is switched off)."
                return (directory,"???",false,true)
            end
        end

    elseif ISLINUX()
        dll_name1 = joinpath(directory, "linux", "SimVisInterface_ProfessionalEdition.so")
        if isfile( dll_name1 )
            dll_name = dll_name1
            isProfessionalEdition = true
        else
            dll_name2 = joinpath(directory, "linux", "SimVisInterface_CommunityEdition.so")
            if isfile( dll_name2 )
                dll_name = dll_name2
                isProfessionalEdition = false
            else
                @warn "\nModia3D: *.so of DLR-Visualization library not found. Neither of these files\n" *
                      "   $dll_name1, \n" *
                      "   $dll_name2, \n" *
                      "exists. Check whether ENV[\"DLR_VISUALIZATION\"] is correct." *
                      "\nNo Renderer is used in Modia3D (so, animation is switched off)."
                return (directory, "???", false, true)
            end
        end
    else
        @warn "\nModia3D: DLR Visualization library only supported for Windows or Linux.\n" *
              "\nNo Renderer is used in Modia3D (so, animation is switched off)."
        return (directory, "???", false, true)
    end

    # Try to open the found DLL/SO
    dll = Libdl.dlopen_e(dll_name)
    if dll != C_NULL
        Libdl.dlclose(dll)
    else
        @warn "\nModia3D: DLR Visualization interface library:" *
              "\n   $dll_name" *
              "\nexists, but could not be opened with Libdl.dlopen_e." *
              "\nNo Renderer is used in Modia3D (so, animation is switched off)."
        return (directory, dll_name, false, true)
    end

    # Print info message
    if isProfessionalEdition
        println("   Renderer: Professional edition of the DLR_Visualization library.\n",
                "             (", dll_name, ")")

    else
        println("   Renderer: Community edition of the DLR_Visualization library",
                "\n             (", dll_name, ",",
                "\n              -> the renderer supports only a subset of the Modia3D functionality).\n")
    end

    return (directory, dll_name, isProfessionalEdition, false)
end


text:
# This file was generated, do not modify it.

using MLJ
using UrlDownload
using PyPlot
ioff() # hide
import DataFrames: DataFrame, describe, names, select!
using Statistics

@load LinearRegressor pkg=MLJLinearModels

data_repo = "https://raw.githubusercontent.com/tlienart/DataScienceTutorialsData.jl/master/data"

url_power   = data_repo * "/power_syst/DE_power_hourly.csv"
url_weather = data_repo * "/power_syst/DE_weather_data.csv"

power   = DataFrame(urldownload(url_power))
weather = DataFrame(urldownload(url_weather));

describe(power, :mean, :nmissing)

first(describe(weather, :mean, :nmissing), 20)

select!(power, [
    :utc_timestamp,
    :DE_solar_generation_actual,
    :DE_wind_generation_actual]);

colnames = names(weather)

filter_by_name(name, cols) =
    filter(cn -> occursin(name, String(cn)), cols)

wind   = weather[:, filter_by_name("windspeed", colnames)]
temp   = weather[:, filter_by_name("temperature", colnames)]
raddir = weather[:, filter_by_name("radiation_direct", colnames)]
raddif = weather[:, filter_by_name("radiation_diffuse", colnames)];

dfs = [wind, temp, raddir, raddif]
col_mean = [:windspeed_mean, :temp_mean, :raddir_mean, :raddif_mean];

nrows = size(first(dfs), 1)
for (df, name) in zip(dfs, col_mean)
    df[!, name] = zeros(nrows)
    for (i, row) in enumerate(eachrow(df))
      df[i, name] = mean(row)
    end
end;

data = DataFrame(
    Timestamp     = weather.utc_timestamp,
    Solar_gen     = power.DE_solar_generation_actual,
    Wind_gen      = power.DE_wind_generation_actual,
    Windspeed     = wind.windspeed_mean,
    Temperature   = temp.temp_mean,
    Radiation_dir = raddir.raddir_mean,
    Radiation_dif = raddif.raddif_mean);

describe(data, :mean, :median, :nmissing)

schema(data)

coerce!(data, :Wind_gen => Continuous)
coerce!(data, :Solar_gen => Continuous)

schema(data)

figure(figsize=(8, 6))
hist(data.Solar_gen, color="blue", edgecolor="white", bins=100,
     density=true, alpha=0.5)
xlabel("Solar power generation (MWh)", fontsize=14)
ylabel("Frequency", fontsize=14)
xticks(fontsize=12)
yticks([0, 1e-3, 2e-3], fontsize=12)
savefig(joinpath(@OUTPUT, "hist_solar.svg")) # hide

figure(figsize=(8, 6))
hist(data.Wind_gen, color="blue", edgecolor = "white", bins=50,
     density=true, alpha=0.5)
xlabel("Wind power generation (MWh)", fontsize=14)
ylabel("Frequency", fontsize=14)

savefig(joinpath(@OUTPUT, "hist_wind.svg")) # hide

fig = figure(figsize=(15, 15))

subplot(221)
scatter(data.Solar_gen, data.Radiation_dir)
xlabel("Solar power (kW)", fontsize=14)
ylabel("Solar radiation - directional", fontsize=14)

subplot(222)
scatter(data.Solar_gen, data.Radiation_dif)
xlabel("Solar power (kW)", fontsize=14)
ylabel("Solar radiation - diffuse", fontsize=14)

subplot(223)
scatter(data.Solar_gen, data.Windspeed)
xlabel("Solar power (kW)", fontsize=14)
ylabel("Wind speed (m/s)", fontsize=14)

subplot(224)
scatter(data.Solar_gen, data.Temperature)
xlabel("Solar power (kW)", fontsize=14)
ylabel("Temperature (C)", fontsize=14)

savefig(joinpath(@OUTPUT, "solar_scatter.png"), bbox_inches="tight") # hide

fig = figure(figsize=(15, 15))

subplot(221)
scatter(data.Wind_gen, data.Radiation_dir)
xlabel("Wind power (kW)", fontsize=14)
ylabel("Solar radiation - directional", fontsize=14)

subplot(222)
scatter(data.Wind_gen, data.Radiation_dif)
xlabel("Wind power (kW)", fontsize=14)
ylabel("Solar radiation - diffuse", fontsize=14)

subplot(223)
scatter(data.Wind_gen, data.Windspeed)
xlabel("Wind power (kW)", fontsize=14)
ylabel("Wind speed (m/s)", fontsize=14)

subplot(224)
scatter(data.Wind_gen, data.Temperature)
xlabel("Wind power (kW)", fontsize=14)
ylabel("Temperature (C)", fontsize=14)

savefig(joinpath(@OUTPUT, "wind_scatter.png"), bbox_inches="tight") # hide

y_wind = data.Wind_gen
X = data[:, [:Windspeed, :Temperature, :Radiation_dir, :Radiation_dif]];

train, test = partition(eachindex(y_wind), 0.7, shuffle=true, rng=5);

linReg = LinearRegressor()
m_linReg = machine(linReg, X, y_wind)
fit!(m_linReg, rows=train);

y_hat = predict(m_linReg, rows=test);

figure(figsize=(8, 6))
plot(y_hat, color="blue", label="Predicted")
plot(y_wind[test], color="red", label="Observed")
xlabel("Time", fontsize=14)
ylabel("Power generation", fontsize=14)
xticks([])
yticks(fontsize=12)
xlim(0, 100)
legend(fontsize=14)

savefig(joinpath(@OUTPUT, "obs_v_pred.svg")) # hide

rms(y_wind[train], predict(m_linReg, rows=train))

rms(y_wind[test], y_hat)

res = y_hat .- y_wind[test];

figure(figsize=(12, 6))
stem(res)
xlim(0, length(res))

savefig(joinpath(@OUTPUT, "residuals.png")) # hide

figure(figsize=(8, 6))
hist(res, color="blue", edgecolor="white", bins=50,
     density=true, alpha=0.5)

savefig(joinpath(@OUTPUT, "hist_residuals.svg")) # hide


text:
"""
    VaR(v::AbstractArray,p::Real;rev::Bool=false)

The `p`th quantile of the vector `v` is the Value at Risk. Assumes more positive values are higher risk measures, so a higher p will return a more positive number, but this can be reversed if `rev` is `true`.

Also can be called with `ValueAtRisk(...)`.
"""
function VaR(v,p;rev=false)
    if rev
        return StatsBase.quantile(v,1-p)
    else
        return StatsBase.quantile(v,p)
    end
end
"""
[VaR](@ref)
"""
ValueAtRisk = VaR

"""
    CTE(v::AbstractArray,p::Real;rev::Bool=false)

The average of the values ≥ the `p`th percentile of the vector `v` is the Conditiona Tail Expectation. Assumes more positive values are higher risk measures, so a higher p will return a more positive number, but this can be reversed if `rev` is `true`.

May also be called with `ConditionalTailExpectation(...)`.

Also known as Tail Value at Risk (TVaR), or Tail Conditional Expectation (TCE)
"""
function CTE(v,p;rev=false)
    # filter has the "or approximately equalt to quantile" because
    # of floating point path might make the quantile slightly off from the right indexing
    # e.g. if values should capture <= q, where q should be 10 but is calculated to be
    # 9.99999...
    if rev
        q = StatsBase.quantile(v,1-p)
        filter = (v .<= q) .| (v .≈ q)
    else
        q = StatsBase.quantile(v,p)
        filter = (v .>= q) .| (v .≈ q)
    end

    return sum(v[filter]) / sum(filter)

end

"""
[CTE](@ref)
"""
ConditionalTailExpectation = CTE

max_stars_repo_name: alecloudenback/ActuaryUtilities

max_issues_repo_name: alecloudenback/ActuaryUtilities

max_forks_repo_name: alecloudenback/ActuaryUtilities


text:
# Showcasing the available analysis tools for the Bayes (Junction) tree.
# using Revise

using IncrementalInference
using DistributedFactorGraphs # For `isSolvable` function.
using Combinatorics # For creating the variable ordering `permutations`.
using SuiteSparse.CHOLMOD: SuiteSparse_long # For CCOLAMD constraints.
using Gadfly # For histogram and scatter plots.
Gadfly.set_default_plot_size(35cm,25cm)

latex_fonts = Theme(major_label_font="CMU Serif", major_label_font_size=16pt,
                    minor_label_font="CMU Serif", minor_label_font_size=14pt,
                    key_title_font="CMU Serif", key_title_font_size=12pt,
                    key_label_font="CMU Serif", key_label_font_size=10pt)
Gadfly.push_theme(latex_fonts)

# Get tree for each variable ordering in a factor graph.
fg = generateGraph_Kaess(graphinit=false)
all_trees = getAllTrees(deepcopy(fg))

# scores stores: (tree key ID, nnz, cost fxn 1, cost fxn 2).
unsorted_scores = Vector{Tuple{Int, Float64, Float64, Float64}}()
for key in keys(all_trees)
    e = all_trees[key] # (Bayes tree, var order, nnz
    tree = e[1] # Get the Bayes tree.
    cost1 = getTreeCost_01(tree)
    cost2 = getTreeCost_02(tree)
    push!(unsorted_scores, (key, e[3], cost1, cost2))
end

# Sort them to make sure the keys are in order.
scores = sort(unsorted_scores)

# Separate scores into vectors for plotting.
all_nnzs = (x->(x[2])).(scores)
costs_01 = (x->(x[3])).(scores)
costs_02 = (x->(x[4])).(scores)

min_ids_02 = findall(x->x == minimum(costs_02), costs_02)
max_ids_02 = findall(x->x == maximum(costs_02), costs_02)

min_ids_nnz = findall(x->x == minimum(all_nnzs), all_nnzs)
max_ids_nnz = findall(x->x == maximum(all_nnzs), all_nnzs)

# Find the intersection between best on both rubrics (lower left quadrant).
best_ids = findall(x->x in min_ids_02, min_ids_nnz)
# Find good factorizations but bad trees (upper left quadrant).
bad_trees_good_mats_ids = findall(x->x in max_ids_02, min_ids_nnz)
# Find good trees with bad matrix factorizations (lower right quadrant).
good_trees_bad_mats_ids = min_ids_02[findall(x->x == maximum(all_nnzs[min_ids_02]), all_nnzs[min_ids_02])]

# Get AMDs variable ordering.
amd_ordering = getEliminationOrder(fg)
amd_tree = buildTreeReset!(deepcopy(fg), amd_ordering)
amd_tree_nnz = nnzTree(amd_tree)
amd_tree_cost02 = getTreeCost_02(amd_tree)

# Get CCOLAMD variable ordering. First bring in CCOLAMD.
include(normpath(Base.find_package("IncrementalInference"), "..", "ccolamd.jl"))
A, varsym, fctsym = getBiadjacencyMatrix(fg)
colamd_ordering = varsym[Ccolamd.ccolamd(A)]
colamd_tree = buildTreeReset!(deepcopy(fg), colamd_ordering)
colamd_tree_nnz = nnzTree(colamd_tree)
colamd_tree_cost02 = getTreeCost_02(colamd_tree)

# Now add the iSAM2 constraint.
cons = zeros(SuiteSparse_long, length(A.colptr) - 1)
cons[findall(x->x == :x3, varsym)[1]] = 1 # NOTE(tonioteran) hardcoded for Kaess' example.
ccolamd_ordering = varsym[Ccolamd.ccolamd(A, cons)]
ccolamd_tree = buildTreeReset!(deepcopy(fg), ccolamd_ordering)
ccolamd_tree_nnz = nnzTree(ccolamd_tree)
ccolamd_tree_cost02 = getTreeCost_02(ccolamd_tree)

# Plot data points and underlying histogram.
bincnt = 20
layers = []
push!(layers, Gadfly.layer(x=[amd_tree_nnz],
                           y=[amd_tree_cost02],
                           Theme(default_color=colorant"green")))
push!(layers, Gadfly.layer(x=[colamd_tree_nnz],
                           y=[colamd_tree_cost02],
                           Theme(default_color=colorant"blue")))
push!(layers, Gadfly.layer(x=[ccolamd_tree_nnz],
                           y=[ccolamd_tree_cost02],
                           Theme(default_color=colorant"red")))
push!(layers, Gadfly.layer(x=all_nnzs,
                           y=costs_02,
                           Geom.hexbin(xbincount=bincnt, ybincount=bincnt)))
pl = Gadfly.plot(layers...,
            Guide.xlabel("Number of non zeros [int]"),
            Guide.ylabel("Tree cost [cfxn2]"),
            Guide.manual_color_key("",
             ["AMD", "COLAMD", "iSAM2"],
             ["green", "blue", "red"]))

img = SVG("vo_cost_canon_kaess.svg", 6inch, 6inch)
Gadfly.draw(img, pl)


text:
# Based on src/http/ngx_http_parse.c from NGINX copyright Igor Sysoev
#
# Additional changes are licensed under the same terms as NGINX and
# copyright Joyent, Inc. and other Node contributors. All rights reserved.
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
# IN THE SOFTWARE.
#

module Parsers

export Parser, Header, Headers, ByteView, nobytes,
       reset!,
       parseheaders, parsebody,
       messagestarted, headerscomplete, bodycomplete, messagecomplete,
       messagehastrailing,
       ParsingError

import MbedTLS.SSLContext

import ..@debug, ..@debugshow, ..DEBUG_LEVEL
import ..@require, ..precondition_error
import ..compat_findfirst, ..compat_findprev, ..compat_findnext

include("consts.jl")
include("parseutils.jl")

const strict = false # See macro @errifstrict

const nobytes = view(UInt8[], 1:0)
const ByteView = typeof(nobytes)
const Header = Pair{String,String}
const Headers = Vector{Header}

"""
 - `method::String`: the HTTP method
   [RFC7230 3.1.1](https://tools.ietf.org/html/rfc7230#section-3.1.1)
 - `major` and `minor`: HTTP version
   [RFC7230 2.6](https://tools.ietf.org/html/rfc7230#section-2.6)
 - `target::String`: request target
   [RFC7230 5.3](https://tools.ietf.org/html/rfc7230#section-5.3)
 - `status::Int`: response status
   [RFC7230 3.1.2](https://tools.ietf.org/html/rfc7230#section-3.1.2)
"""

mutable struct Message
    method::String
    major::Int
    minor::Int
    target::String
    status::Int

    Message() = reset!(new())
end

function reset!(m::Message)
    m.method = ""
    m.major = 0
    m.minor = 0
    m.target = ""
    m.status = 0
    return m
end

"""
The parser separates a raw HTTP Message into its component parts.

If the input data is invalid the Parser throws a `ParsingError`.

The parser processes a single HTTP Message. If the input stream contains
multiple Messages the Parser stops at the end of the first Message.
The `parseheaders` and `parsebody` functions return a `SubArray` containing the
unuses portion of the input.

The Parser does not interpret the Message Headers except as needed
to parse the Message Body. It is beyond the scope of the Parser to deal
with repeated header fields, multi-line values, cookies or case normalization.

The Parser has no knowledge of the high-level `Request` and `Response` structs
defined in `Messages.jl`. The Parser has it's own low level
[`Message`](@ref) struct that represents both Request and Response
Messages.
"""

mutable struct Parser

    # state
    state::UInt8
    chunk_length::UInt64
    trailing::Bool
    fieldbuffer::IOBuffer
    valuebuffer::IOBuffer

    # output
    message::Message

    function Parser()
        p = new()
        p.fieldbuffer = IOBuffer()
        p.valuebuffer = IOBuffer()
        p.message = Message()
        return reset!(p)
    end
end

"""
    reset!(::Parser)

Revert `Parser` to unconfigured state.
"""

function reset!(p::Parser)
    p.state = s_start_req_or_res
    p.chunk_length = 0
    p.trailing = false
    truncate(p.fieldbuffer, 0)
    truncate(p.valuebuffer, 0)
    reset!(p.message)
    return p
end

"""
    messagestarted(::Parser)

Has the `Parser` begun processng a Message?
"""

messagestarted(p::Parser) = p.state != s_start_req_or_res

"""
    headerscomplete(::Parser)

Has the `Parser` processed the entire Message Header?
"""

headerscomplete(p::Parser) = p.state > s_headers_done

"""
    bodycomplete(::Parser)

Has the `Parser` processed the Message Body?
"""

bodycomplete(p::Parser) = p.state == s_message_done ||
                          p.state == s_trailer_start

"""
    messagecomplete(::Parser)

Has the `Parser` processed the entire Message?
"""

messagecomplete(p::Parser) = p.state >= s_message_done

"""
    messagehastrailing(::Parser)

Is the `Parser` ready to process trailing headers?
"""
messagehastrailing(p::Parser) = p.trailing

isrequest(p::Parser) = p.message.status == 0

"""
The [`Parser`] input was invalid.

Fields:
 - `code`, internal error code
 - `state`, internal parsing state.
 - `status::Int`, HTTP response status.
 - `msg::String`, error message.
"""

struct ParsingError <: Exception
    code::Symbol
    state::UInt8
    status::Int
    msg::String
    bytes::ByteView
    p::Int
end

function ParsingError(parser::Parser, code::Symbol, bytes, p)
    ParsingError(code, parser.state, parser.message.status, "", bytes, p)
end

function Base.show(io::IO, e::ParsingError)
    println(io, string("HTTP.ParsingError: ",
                       get(ERROR_MESSAGES, e.code, "?"), ", ",
                       ParsingStateCode(e.state), ", ",
                       e.status,
                       e.msg == "" ? "" : "\n",
                       e.msg))
    s = String(e.bytes)
    l = compat_findprev(s, '\n', e.p)
    r = compat_findnext(s, '\n', e.p-1)
    error_line = SubString(s, l == 0 ? 1 : l+1, r == 0 ? length(e.bytes) : r)
    println(io, "\"", chomp(error_line), "\"\n",
                lpad("", e.p - l, " "), "^")
end

macro err(code)
    esc(:(parser.state = p_state; throw(ParsingError(parser, $code, bytes, p))))
end

macro errorif(cond, err)
    esc(:($cond && @err($err)))
end

macro errorifstrict(cond)
    strict ? esc(:(@errorif($cond, :HPE_STRICT))) : :()
end

macro passert(cond)
    DEBUG_LEVEL > 1 ? esc(:(@assert $cond)) : :()
end

macro methodstate(meth, i, char)
    return esc(:(Int($meth) << Int(16) | Int($i) << Int(8) | Int($char)))
end

function parse_token(bytes, len, p, buffer; allowed='a')
    start = p
    while p <= len
        @inbounds ch = Char(bytes[p])
        if !istoken(ch) && ch != allowed
            break
        end
        p += 1
    end
    @passert p <= len + 1

    write(buffer, view(bytes, start:p-1))

    if p > len
        return len, false
    else
        return p, true
    end
end

"""
    parseheaders(::Parser, bytes) do h::Pair{String,String} ... -> excess

Read headers from `bytes`, passing each field/value pair to `f`.
Returns a `SubArray` containing bytes not parsed.

e.g.
```
excess = parseheaders(p, bytes) do (k,v)
    println("\$k: \$v")
end
```
"""

function parseheaders(f, p, bytes)
    v = Vector{UInt8}(bytes)
    parseheaders(f, p, view(v, 1:length(v)))
end

function parseheaders(onheader::Function #=f(::Pair{String,String}) =#,
                      parser::Parser, bytes::ByteView)::ByteView

    @require !isempty(bytes)
    @require messagehastrailing(parser) || !headerscomplete(parser)

    len = length(bytes)
    p_state = parser.state
    @debug 3 "parseheaders(parser.state=$(ParsingStateCode(p_state))), " *
             "$len-bytes:\n" * escapelines(String(collect(bytes))) * ")"

    p = 0
    while p < len && p_state <= s_headers_done

        @debug 4 string("top of while($p < $len) \"",
                        Base.escape_string(string(Char(bytes[p+1]))), "\" ",
                        ParsingStateCode(p_state))
        p += 1
        @inbounds ch = Char(bytes[p])

        if p_state == s_start_req_or_res
            (ch == CR || ch == LF) && continue

            p_state = s_start_req
            p -= 1

        elseif p_state == s_res_first_http_major
            @errorif(!isnum(ch), :HPE_INVALID_VERSION)
            parser.message.major = Int(ch - '0')
            p_state = s_res_http_major

        # major HTTP version or dot
        elseif p_state == s_res_http_major
            if ch == '.'
                p_state = s_res_first_http_minor
                continue
            end
            @errorif(!isnum(ch), :HPE_INVALID_VERSION)
            parser.message.major *= Int(10)
            parser.message.major += Int(ch - '0')
            @errorif(parser.message.major > 999, :HPE_INVALID_VERSION)

        # first digit of minor HTTP version
        elseif p_state == s_res_first_http_minor
            @errorif(!isnum(ch), :HPE_INVALID_VERSION)
            parser.message.minor = Int(ch - '0')
            p_state = s_res_http_minor

        # minor HTTP version or end of request line
        elseif p_state == s_res_http_minor
            if ch == ' '
                p_state = s_res_first_status_code
                continue
            end
            @errorif(!isnum(ch), :HPE_INVALID_VERSION)
            parser.message.minor *= Int(10)
            parser.message.minor += Int(ch - '0')
            @errorif(parser.message.minor > 999, :HPE_INVALID_VERSION)

        elseif p_state == s_res_first_status_code
            if !isnum(ch)
                ch == ' ' && continue
                @err(:HPE_INVALID_STATUS)
            end
            parser.message.status = Int(ch - '0')
            p_state = s_res_status_code

        elseif p_state == s_res_status_code
            if !isnum(ch)
                if ch == ' '
                    p_state = s_res_status_start
                elseif ch == CR
                    p_state = s_res_line_almost_done
                elseif ch == LF
                    p_state = s_header_field_start
                else
                    @err(:HPE_INVALID_STATUS)
                end
            else
                parser.message.status *= Int(10)
                parser.message.status += Int(ch - '0')
                @errorif(parser.message.status > 999, :HPE_INVALID_STATUS)
            end

        elseif p_state == s_res_status_start
            if ch == CR
                p_state = s_res_line_almost_done
            elseif ch == LF
                p_state = s_header_field_start
            else
                p_state = s_res_status
            end

        elseif p_state == s_res_status
            if ch == CR
                p_state = s_res_line_almost_done
            elseif ch == LF
                p_state = s_header_field_start
            end

        elseif p_state == s_res_line_almost_done
            @errorifstrict(ch != LF)
            p_state = s_header_field_start

        elseif p_state == s_start_req
            (ch == CR || ch == LF) && continue

            @errorif(!istoken(ch), :HPE_INVALID_METHOD)

            p_state = s_req_method
            p -= 1

        elseif p_state == s_req_method

            p, complete = parse_token(bytes, len, p, parser.valuebuffer)

            if complete
                parser.message.method = take!(parser.valuebuffer)
                @inbounds ch = Char(bytes[p])
                if parser.message.method == "HTTP" && ch == '/'
                    p_state = s_res_first_http_major
                elseif ch == ' '
                    p_state = s_req_spaces_before_target
                else
                    @err(:HPE_INVALID_METHOD)
                end
            end

        elseif p_state == s_req_spaces_before_target
            ch == ' ' && continue
            p_state = s_req_target
            p -= 1

        elseif (p_state ==  s_req_target)
            start = p
            while p <= len
                @inbounds ch = Char(bytes[p])
                if @anyeq(ch, ' ', CR, LF)
                    if ch == ' '
                        p_state = s_req_http_start
                    else
                        parser.message.major = Int(0)
                        parser.message.minor = Int(9)
                        p_state = ifelse(ch == CR, s_req_line_almost_done,
                                                   s_header_field_start)
                    end
                    break
                end
                p += 1
            end
            @passert p <= len + 1

            write(parser.valuebuffer, collect(view(bytes, start:p-1)))

            if p_state >= s_req_http_start
                target = String(take!(parser.valuebuffer))
                @debugshow 4 target
                parser.message.target = target
                @errorif(isempty(target) ||
                         target[1] == '.' ||
                         startswith(target, "HTTP/"),
                         :HPE_INVALID_TARGET)
                parser.message.target = target
            end

            p = min(p, len)

        elseif p_state == s_req_http_start
            if ch == 'H'
                p_state = s_req_http_H
            elseif ch == ' '
            else
                @err(:HPE_INVALID_CONSTANT)
            end

        elseif p_state == s_req_http_H
            @errorifstrict(ch != 'T')
            p_state = s_req_http_HT

        elseif p_state == s_req_http_HT
            @errorifstrict(ch != 'T')
            p_state = s_req_http_HTT

        elseif p_state == s_req_http_HTT
            @errorifstrict(ch != 'P')
            p_state = s_req_http_HTTP

        elseif p_state == s_req_http_HTTP
            @errorifstrict(ch != '/')
            p_state = s_req_first_http_major

        # first digit of major HTTP version
        elseif p_state == s_req_first_http_major
            @errorif(ch < '1' || ch > '9', :HPE_INVALID_VERSION)
            parser.message.major = Int(ch - '0')
            p_state = s_req_http_major

        # major HTTP version or dot
        elseif p_state == s_req_http_major
            if ch == '.'
                p_state = s_req_first_http_minor
            elseif !isnum(ch)
                @err(:HPE_INVALID_VERSION)
            else
                parser.message.major *= Int(10)
                parser.message.major += Int(ch - '0')
                @errorif(parser.message.major > 999, :HPE_INVALID_VERSION)
            end

        # first digit of minor HTTP version
        elseif p_state == s_req_first_http_minor
            @errorif(!isnum(ch), :HPE_INVALID_VERSION)
            parser.message.minor = Int(ch - '0')
            p_state = s_req_http_minor

        # minor HTTP version or end of request line
        elseif p_state == s_req_http_minor
            if ch == CR
                p_state = s_req_line_almost_done
            elseif ch == LF
                p_state = s_header_field_start
            else
                # FIXME allow spaces after digit?
                @errorif(!isnum(ch), :HPE_INVALID_VERSION)
                parser.message.minor *= Int(10)
                parser.message.minor += Int(ch - '0')
                @errorif(parser.message.minor > 999, :HPE_INVALID_VERSION)
            end

        # end of request line
        elseif p_state == s_req_line_almost_done
            @errorif(ch != LF, :HPE_LF_EXPECTED)
            p_state = s_header_field_start

        elseif p_state == s_header_field_start ||
               p_state == s_trailer_start
            if ch == CR
                p_state = s_headers_almost_done
            elseif ch == LF
                # they might be just sending \n instead of \r\n so this would be
                # the second \n to denote the end of headers
                p_state = s_headers_almost_done
                p -= 1
            else
                c = (!strict && ch == ' ') ? ' ' : tokens[Int(ch)+1]
                @errorif(c == Char(0), :HPE_INVALID_HEADER_TOKEN)
                p_state = s_header_field
                p -= 1
            end

        elseif p_state == s_header_field

            p, complete = parse_token(bytes, len, p, parser.fieldbuffer;
                                      allowed = ' ')
            if complete
                @inbounds ch = Char(bytes[p])
                @errorif(ch != ':', :HPE_INVALID_HEADER_TOKEN)
                p_state = s_header_value_discard_ws
            end

        elseif p_state == s_header_value_discard_ws
            (ch == ' ' || ch == '\t') && continue
            if ch == CR
                p_state = s_header_value_discard_ws_almost_done
                continue
            end
            if ch == LF
                p_state = s_header_value_discard_lws
                continue
            end
            p_state = s_header_value_start
            p -= 1
        elseif p_state == s_header_value_start
            p_state = s_header_value
            c = lower(ch)

            write(parser.valuebuffer, bytes[p])

        elseif p_state == s_header_value
            start = p
            while p <= len
                @inbounds ch = Char(bytes[p])
                @debug 4 Base.escape_string(string('\'', ch, '\''))
                if ch == CR
                    p_state = s_header_almost_done
                    break
                elseif ch == LF
                    p_state = s_header_value_lws
                    break
                elseif strict && !isheaderchar(ch)
                    @err(:HPE_INVALID_HEADER_TOKEN)
                end

                c = lower(ch)

                @debugshow 4 h
                crlf = compat_findfirst(x->(x == bCR || x == bLF), view(bytes, p:len))
                p = crlf == 0 ? len : p + crlf - 2

                p += 1
            end
            @passert p <= len + 1

            write(parser.valuebuffer, view(bytes, start:p-1))

            if p_state != s_header_value
                onheader(String(take!(parser.fieldbuffer)) =>
                         String(take!(parser.valuebuffer)))
            end

            p = min(p, len)

        elseif p_state == s_header_almost_done
            @errorif(ch != LF, :HPE_LF_EXPECTED)
            p_state = s_header_value_lws

        elseif p_state == s_header_value_lws
            p -= 1
            if ch == ' ' || ch == '\t'
                p_state = s_header_value_start
            else
                # finished the header
                p_state = s_header_field_start
            end

        elseif p_state == s_header_value_discard_ws_almost_done
            @errorifstrict(ch != LF)
            p_state = s_header_value_discard_lws

        elseif p_state == s_header_value_discard_lws
            if ch == ' ' || ch == '\t'
                p_state = s_header_value_discard_ws
            else
                # header value was empty
                p_state = s_header_field_start
                onheader(String(take!(parser.fieldbuffer)) => "")
                p -= 1
            end

        elseif p_state == s_headers_almost_done
            @errorifstrict(ch != LF)
            p -= 1
            if parser.trailing
                # End of a chunked request
                p_state = s_message_done
            else
                p_state = s_headers_done
            end

        elseif p_state == s_headers_done
            @errorifstrict(ch != LF)

            p_state = s_body_start
        else
            @err :HPE_INVALID_INTERNAL_STATE
        end
    end

    @assert p <= len
    @assert p == len ||
            p_state == s_message_done ||
            p_state == s_body_start

    # Consume trailing end of line after message.
    if p_state == s_message_done
        while p < len
            ch = Char(bytes[p + 1])
            if ch != CR && ch != LF
                break
            end
            p += 1
        end
    end

    @debug 3 "parseheaders() exiting $(ParsingStateCode(p_state))"

    parser.state = p_state
    return view(bytes, p+1:len)
end

"""
    parsebody(::Parser, bytes) -> data, excess

Parse body data from `bytes`.
Returns decoded `data` and `excess` bytes not parsed.
"""

function parsebody(p, bytes)
    v = Vector{UInt8}(bytes)
    parsebody(p, view(v, 1:length(v)))
end

function parsebody(parser::Parser, bytes::ByteView)::Tuple{ByteView,ByteView}

    @require !isempty(bytes)
    @require headerscomplete(parser)

    if parser.state == s_body_start
        parser.state = s_chunk_size_start
    end

    len = length(bytes)
    p_state = parser.state
    @debug 3 "parsebody(parser.state=$(ParsingStateCode(p_state))), " *
             "$len-bytes:\n" * escapelines(String(collect(bytes))) * ")"

    result = nobytes

    p = 0
    while p < len && result == nobytes && p_state != s_trailer_start

        @debug 4 string("top of while($p < $len) \"",
                        Base.escape_string(string(Char(bytes[p+1]))), "\" ",
                        ParsingStateCode(p_state))
        p += 1
        @inbounds ch = Char(bytes[p])

        if p_state == s_chunk_size_start

            unhex_val = unhex[Int(ch)+1]
            @errorif(unhex_val == -1, :HPE_INVALID_CHUNK_SIZE)

            parser.chunk_length = unhex_val
            p_state = s_chunk_size

        elseif p_state == s_chunk_size
            if ch == CR
                p_state = s_chunk_size_almost_done
            else
                unhex_val = unhex[Int(ch)+1]
                if unhex_val == -1
                    if ch == ';' || ch == ' '
                        p_state = s_chunk_parameters
                        continue
                    end
                    @err(:HPE_INVALID_CHUNK_SIZE)
                end
                t = parser.chunk_length
                t *= UInt64(16)
                t += UInt64(unhex_val)

                # Overflow? Test against a conservative limit for simplicity.
                @debugshow 4 Int(parser.chunk_length)
                if div(typemax(UInt64) - 16, 16) < t
                    @err(:HPE_INVALID_CONTENT_LENGTH)
                end
                parser.chunk_length = t
            end

        elseif p_state == s_chunk_parameters
            # just ignore this?. FIXME check for overflow?
            if ch == CR
                p_state = s_chunk_size_almost_done
            end

        elseif p_state == s_chunk_size_almost_done
            @errorifstrict(ch != LF)

            if parser.chunk_length == 0
                parser.trailing = true
                p_state = s_trailer_start
            else
                p_state = s_chunk_data
            end

        elseif p_state == s_chunk_data
            to_read = Int(min(parser.chunk_length, len - p + 1))

            @passert parser.chunk_length != 0

            @passert result == nobytes
            result = view(bytes, p:p + to_read - 1)
            parser.chunk_length -= to_read
            p += to_read - 1

            if parser.chunk_length == 0
                p_state = s_chunk_data_almost_done
            end

        elseif p_state == s_chunk_data_almost_done
            @passert parser.chunk_length == 0
            @errorifstrict(ch != CR)
            p_state = s_chunk_data_done

        elseif p_state == s_chunk_data_done
            @errorifstrict(ch != LF)
            p_state = s_chunk_size_start

        else
            @err :HPE_INVALID_INTERNAL_STATE
        end
    end

    @assert p <= len
    @assert p == len ||
            result != nobytes ||
            p_state == s_trailer_start

    @debug 3 "parsebody() exiting $(ParsingStateCode(p_state))"

    parser.state = p_state
    return result, view(bytes, p+1:len)
end

const ERROR_MESSAGES = Dict(
    :HPE_INVALID_VERSION => "invalid HTTP version",
    :HPE_INVALID_STATUS => "invalid HTTP status code",
    :HPE_INVALID_METHOD => "invalid HTTP method",
    :HPE_INVALID_TARGET => "invalid HTTP request target",
    :HPE_LF_EXPECTED => "LF character expected",
    :HPE_INVALID_HEADER_TOKEN => "invalid character in header",
    :HPE_INVALID_CONTENT_LENGTH => "invalid character in content-length header",
    :HPE_INVALID_CHUNK_SIZE => "invalid character in chunk size header",
    :HPE_INVALID_CONSTANT => "invalid constant string",
    :HPE_INVALID_INTERNAL_STATE => "encountered unexpected internal state",
    :HPE_STRICT => "strict mode assertion failed",
)

"""
Tokens as defined by rfc 2616. Also lowercases them.
        token       = 1*<any CHAR except CTLs or separators>
     separators     = "(" | ")" | "<" | ">" | "@"
                    | "," | ";" | ":" | "\" | <">
                    | "/" | "[" | "]" | "?" | "="
                    | "{" | "}" | SP | HT
"""

const tokens = Char[
#=   0 nul    1 soh    2 stx    3 etx    4 eot    5 enq    6 ack    7 bel  =#
        0,       0,       0,       0,       0,       0,       0,       0,
#=   8 bs     9 ht    10 nl    11 vt    12 np    13 cr    14 so    15 si   =#
        0,       0,       0,       0,       0,       0,       0,       0,
#=  16 dle   17 dc1   18 dc2   19 dc3   20 dc4   21 nak   22 syn   23 etb =#
        0,       0,       0,       0,       0,       0,       0,       0,
#=  24 can   25 em    26 sub   27 esc   28 fs    29 gs    30 rs    31 us  =#
        0,       0,       0,       0,       0,       0,       0,       0,
#=  32 sp    33  !    34  "    35  #    36  $    37  %    38  &    39  '  =#
        0,      '!',      0,      '#',     '$',     '%',     '&',    '\'',
#=  40  (    41  )    42  *    43  +    44  ,    45  -    46  .    47  /  =#
        0,       0,      '*',     '+',      0,      '-',     '.',      0,
#=  48  0    49  1    50  2    51  3    52  4    53  5    54  6    55  7  =#
       '0',     '1',     '2',     '3',     '4',     '5',     '6',     '7',
#=  56  8    57  9    58  :    59  ;    60  <    61  =    62  >    63  ?  =#
       '8',     '9',      0,       0,       0,       0,       0,       0,
#=  64  @    65  A    66  B    67  C    68  D    69  E    70  F    71  G  =#
        0,      'a',     'b',     'c',     'd',     'e',     'f',     'g',
#=  72  H    73  I    74  J    75  K    76  L    77  M    78  N    79  O  =#
       'h',     'i',     'j',     'k',     'l',     'm',     'n',     'o',
#=  80  P    81  Q    82  R    83  S    84  T    85  U    86  V    87  W  =#
       'p',     'q',     'r',     's',     't',     'u',     'v',     'w',
#=  88  X    89  Y    90  Z    91  [    92  \    93  ]    94  ^    95  _  =#
       'x',     'y',     'z',      0,       0,       0,      '^',     '_',
#=  96  `    97  a    98  b    99  c   100  d   101  e   102  f   103  g  =#
       '`',     'a',     'b',     'c',     'd',     'e',     'f',     'g',
#= 104  h   105  i   106  j   107  k   108  l   109  m   110  n   111  o  =#
       'h',     'i',     'j',     'k',     'l',     'm',     'n',     'o',
#= 112  p   113  q   114  r   115  s   116  t   117  u   118  v   119  w  =#
       'p',     'q',     'r',     's',     't',     'u',     'v',     'w',
#= 120  x   121  y   122  z   123  {   124  |   125  }   126  ~   127 del =#
       'x',     'y',     'z',      0,      '|',      0,      '~',       0 ]

istoken(c) = tokens[UInt8(c)+1] != Char(0)

const unhex = Int8[
     -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1
    ,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1
    ,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1
    , 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,-1,-1,-1,-1,-1,-1
    ,-1,10,11,12,13,14,15,-1,-1,-1,-1,-1,-1,-1,-1,-1
    ,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1
    ,-1,10,11,12,13,14,15,-1,-1,-1,-1,-1,-1,-1,-1,-1
    ,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1
]

Base.show(io::IO, p::Parser) = print(io, "Parser(",
    "state=", ParsingStateCode(p.state), ", ",
    "trailing=", p.trailing, ", ",
    "message=", p.message, ")")

end # module Parsers


text:
input = 240920:789857

part1() = count(validnum, input)

function validnum(n)
    ds = digits(n)
    ns = zip(ds[begin:end-1], ds[begin+1:end])
    all(a >= b for (a, b) in ns) && any(a == b for (a, b) in ns)
end

function validnum2(n)
    ds = digits(n)
    ns = zip(ds[begin:end-1], ds[begin+1:end])
    all(a >= b for (a, b) in ns) && any(count(==(d), ds) == 2 for d in unique(ds))
end

part2() = count(validnum2, input)

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
@testset "Fitting" begin
  TI = training_image("WalkerLake")[1:20,1:20,1]
  xwalker = Float64[i for i=1:20 for j=1:20]
  ywalker = Float64[j for i=1:20 for j=1:20]
  zwalker = Float64[TI[i,j] for i=1:20 for j=1:20]
  γwalker = EmpiricalVariogram(hcat(xwalker,ywalker)', zwalker, maxlag=15.)

  # variogram types to fit
  Vs = (GaussianVariogram, SphericalVariogram,
        ExponentialVariogram, MaternVariogram)

  # all fits lead to same sill
  γs = [Variography.fit(V, γwalker) for V in Vs]
  for γ in γs
    @test isapprox(sill(γ), 0.054, atol=1e-3)
  end

  # best fit is a Gaussian variogram
  γbest = Variography.fit(Variogram, γwalker)
  @test γbest isa GaussianVariogram
  @test isapprox(sill(γbest), 0.054, atol=1e-3)

  if visualtests
    @plottest begin
      plts = []
      for γ in γs
        plt = plot(γwalker, legend=false)
        plot!(γ, maxlag=15.)
        push!(plts, plt)
      end
      plot(plts...)
    end joinpath(datadir,"Fitting.png") !istravis
  end
end


text:
# Python buffer protocol: this is a NumPy-array-like facility to get raw
# pointers to contiguous buffers of data underlying other objects,
# with support for describing multiple dimensions, strides, etc.
#     (thanks to @jakebolewski for his work on this)

#############################################################################
# mirror of Py_buffer struct in Python Include/object.h

immutable Py_buffer
    buf::Ptr{Void}
    obj::PyPtr
    len::Cssize_t
    itemsize::Cssize_t

    readonly::Cint
    ndim::Cint
    format::Ptr{Cchar}
    shape::Ptr{Cssize_t}
    strides::Ptr{Cssize_t}
    suboffsets::Ptr{Cssize_t}

    # some opaque padding fields to account for differences between
    # Python versions (the structure changed in Python 2.7 and 3.3)
    internal0::Ptr{Void}
    internal1::Ptr{Void}
    internal2::Ptr{Void}
end

type PyBuffer
    buf::Py_buffer
    PyBuffer() = begin
        b = new(Py_buffer(C_NULL, C_NULL, 0, 0,
                          0, 0, C_NULL, C_NULL, C_NULL, C_NULL,
                          C_NULL, C_NULL, C_NULL))
        finalizer(b, pydecref)
        return b
    end
end

function pydecref(o::PyBuffer)
    # note that PyBuffer_Release sets o.obj to NULL, and
    # is a no-op if o.obj is already NULL
    ccall(@pysym(:PyBuffer_Release), Void, (Ptr{PyBuffer},), &o)
    o
end

#############################################################################
# Array-like accessors for PyBuffer.

Base.ndims(b::PyBuffer) = UInt(b.buf.ndim)

# from the Python docs: If shape is NULL as a result of a PyBUF_SIMPLE
# or a PyBUF_WRITABLE request, the consumer must disregard itemsize
# and assume itemsize == 1.
Base.length(b::PyBuffer) = b.buf.shape == C_NULL ? Int(b.buf.len) : Int(div(b.buf.len, b.buf.itemsize))

Base.sizeof(b::PyBuffer) = Int(b.buf.len)
Base.pointer(b::PyBuffer) = b.buf.buf

function Base.size(b::PyBuffer)
    b.buf.ndim <= 1 && return (length(b),)
    @assert b.buf.shape != C_NULL
    return tuple(Int[unsafe_load(b.buf.shape, i) for i=1:b.buf.ndim]...)
end
# specialize size(b, d) for efficiency (avoid tuple construction)
function Base.size(b::PyBuffer, d::Integer)
    d > b.buf.ndim && return 1
    d < 0 && throw(BoundsError())
    b.buf.ndim <= 1 && return length(b)
    @assert b.buf.shape != C_NULL
    return Int(unsafe_load(b.buf.shape, d))
end

# stride in bytes for i-th dimension
function Base.stride(b::PyBuffer, d::Integer)
    d > b.buf.ndim && return length(b) # as in base
    d < 0 && throw(BoundsError())
    if b.buf.strides == C_NULL
        if b.buf.ndim == 1
            return b.buf.shape == C_NULL ? 1 : Int(b.buf.itemsize)
        else
            error("unknown buffer strides")
        end
    end
    return Int(unsafe_load(b.buf.strides, d))
end

iscontiguous(b::PyBuffer) =
    1 == ccall((@pysym :PyBuffer_IsContiguous), Cint,
               (Ptr{PyBuffer}, Cchar), &b, 'A')

#############################################################################
# pybuffer constant values from Include/object.h
const PyBUF_MAX_NDIM = convert(Cint, 64)
const PyBUF_SIMPLE    = convert(Cint, 0)
const PyBUF_WRITABLE  = convert(Cint, 0x0001)
const PyBUF_FORMAT    = convert(Cint, 0x0004)
const PyBUF_ND        = convert(Cint, 0x0008)
const PyBUF_STRIDES        = convert(Cint, 0x0010) | PyBUF_ND
const PyBUF_C_CONTIGUOUS   = convert(Cint, 0x0020) | PyBUF_STRIDES
const PyBUF_F_CONTIGUOUS   = convert(Cint, 0x0040) | PyBUF_STRIDES
const PyBUF_ANY_CONTIGUOUS = convert(Cint, 0x0080) | PyBUF_STRIDES
const PyBUF_INDIRECT       = convert(Cint, 0x0100) | PyBUF_STRIDES

# construct a PyBuffer from a PyObject, if possible
function PyBuffer(o::Union{PyObject,PyPtr}, flags=PyBUF_SIMPLE)
    b = PyBuffer()
    @pycheckz ccall((@pysym :PyObject_GetBuffer), Cint,
                     (PyPtr, Ptr{PyBuffer}, Cint), o, &b, flags)
    return b
end

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

# recursive function to write buffer dimension by dimension, starting at
# dimension d with the given pointer offset (in bytes).
function writedims(io::IO, b::PyBuffer, offset, d)
    n = 0
    s = stride(b, d)
    if d < b.buf.ndim
        for i = 1:size(b,d)
            n += writedims(io, b, offset, d+1)
            offset += s
        end
    else
        @assert d == b.buf.ndim
        p = convert(Ptr{UInt8}, pointer(b)) + offset
        for i = 1:size(b,d)
            # would be nicer not to write this one byte at a time,
            # but the alternative seems to be to create an Array
            # object on each loop iteration.
            for j = 1:b.buf.itemsize
                n += write(io, unsafe_load(p))
                p += 1
            end
            p += s
        end
    end
    return n
end

function Base.write(io::IO, b::PyBuffer)
    b.buf.obj != C_NULL || error("attempted to write NULL buffer")

    if iscontiguous(b)
        # (note that 0-dimensional buffers are always contiguous)
        return write(io, pointer_to_array(convert(Ptr{UInt8}, pointer(b)),
                                          sizeof(b)))
    else
        return writedims(io, b, 0, 1)
    end
end

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


text:
using Tracker, Test, Random, Statistics

Random.seed!(0)

@testset "Tracker" begin

include("tracker.jl")

using Tracker: jacobian

@testset "Jacobian" begin
  A = param(randn(2,2))
  x = randn(2)
  m(x) = A*x
  y = m(x)
  J = jacobian(m,x)
  @test J ≈ A.data
end

end


text:
using Test
using QuantumOptics, CollectiveSpins
using LinearAlgebra

@testset "interaction" begin

cs = CollectiveSpins

r = [1.0,1.0,0.0]
dips = [1.0,0.0,0.0]

pos = cs.geometry.chain(0.4, 3)
dips = [[0.0,0.0,1.0],[0.0,0.0,1.0],[0.0,1.0,0.0]]

S = cs.SpinCollection(pos, dips)
Ω = cs.interaction.OmegaMatrix(S)
Γ = cs.interaction.GammaMatrix(S)

for i=1:length(pos), j=1:length(pos)
    @test isapprox(Gamma_ij(pos[i],pos[j],dips[i],dips[j]), Γ[i,j])
    @test isapprox(Omega_ij(pos[i],pos[j],dips[i],dips[j]), Ω[i,j])
end

N = length(pos)
for i = 1:N
    for j = 1:N
        @test Gamma_ij(pos[1],pos[2],dips[1],dips[2]) <= 1.0
        @test Omega_ij(pos[1],pos[2],dips[1],dips[2]) != NaN

    end
end

# Collective effects
function F(ri, rj, µi_, µj_)
    rij = ri - rj
    μi = normalize(µi_)
    μj = normalize(µj_)
    rij_norm = norm(rij)
    rijn = rij./rij_norm
    if rij_norm == 0
        2/3.
    else
        β = 2π*rij_norm
        µi'*(µj*(sin(β)/β + cos(β)/β^2 - sin(β)/β^3) + (rijn'*µj)*rijn*(-sin(β)/β - 3*cos(β)/β^2 + 3*sin(β)/β^3))
    end
end
function G(ri, rj, µi_, µj_)
    rij = ri - rj
    μi = normalize(µi_)
    μj = normalize(µj_)
    rij_norm = norm(rij)
    rijn = rij./rij_norm
    if rij_norm == 0
        0.0
    else
        β = 2π*rij_norm
        µi'*(µj*(-cos(β)/β + sin(β)/β^2 + cos(β)/β^3) + (rijn'*µj)*rijn*(cos(β)/β - 3*sin(β)/β^2 - 3*cos(β)/β^3))
    end
end

# Test circular polarizations (complex μ)
N = 10
pos = geometry.ring(0.1,N;distance=true)

dips = normalize.([[1, im^i, 0] for i=1:N])

Γmat1 = [CollectiveSpins.interaction.Gamma(pos[i], pos[j], dips[i], dips[j]) for i=1:N, j=1:N]
Ωmat1 = [CollectiveSpins.interaction.Omega(pos[i], pos[j], dips[i], dips[j]) for i=1:N, j=1:N]

Γmat2 = [CollectiveSpins.interaction.Gamma_ij(pos[i], pos[j], dips[i], dips[j]) for i=1:N, j=1:N]
Ωmat2 = [CollectiveSpins.interaction.Omega_ij(pos[i], pos[j], dips[i], dips[j]) for i=1:N, j=1:N]

for i=1:N, j=1:N
    @test isapprox(Γmat1[i,j], Γmat2[i,j])
    @test isapprox(Ωmat1[i,j], Ωmat2[i,j])
end

end #testset


text:
# Generate a Agent with PPOPolicy based on inner env
function create_agent(
    ::Val{:PPO};
    update_freq=256,
    n_env=1,
    seed=37,
    kwargs...
)
    rng = StableRNG(seed)
    inner_env = create_env()
    ns = length(state_space(inner_env))
    na = length(action_space(inner_env))

    init = glorot_uniform(rng)
    return Agent(
        policy = PPOPolicy(
            approximator = ActorCritic(
                actor = GaussianNetwork(
                    pre = Chain(
                        Dense(ns, 64, relu; init = init),
                        Dense(64, 64, relu; init = init),
                    ),
                    μ = Chain(
                        Dense(64, na, tanh; init = init),
                    ),
                    logσ = Chain(
                        Dense(64, na; init = init),
                        x -> clamp.(x, eltype(x)(-10), eltype(x)(1)), # Constrain variance, don't want zero and don't want too big
                    )
                ),
                critic = Chain(
                    Dense(ns, 64, relu; init = init),
                    Dense(64, 64, relu; init = init),
                    Dense(64, 1; init = init),
                ),
                optimizer = ADAM(3e-4),
            ) |> cpu,
            update_freq = update_freq,
            n_random_start = 0,
            update_step = 0,
            γ = 0.99f0,
            λ = 0.95f0,
            clip_range = 0.2f0,
            max_grad_norm = 0.5f0,
            n_microbatches = 32,
            n_epochs = 10, # Collect update_freq*n_env data points, loop over them n_epochs times and randomly split into n_microbatches smaller updates each loop
            actor_loss_weight = 1.0f0,
            critic_loss_weight = 0.5f0,
            entropy_loss_weight = 0.00f0,
            dist = Normal,
            rng = rng,
        ),
        trajectory = PPOTrajectory(;
            capacity = update_freq,
            state = Matrix{Float32} => (ns, n_env),
            action = Matrix{Float32} => (na, n_env),
            action_log_prob = Vector{Float32} => (n_env,),
            reward = Vector{Float32} => (n_env,),
            terminal = Vector{Bool} => (n_env,),
        ),
    )
end

max_stars_repo_name: albheim/ServiceMeshControl

max_issues_repo_name: albheim/ServiceMeshControl

max_forks_repo_name: albheim/ServiceMeshControl


text:
using LinearAlgebra

"""
    implementation of the Limiter
"""

# functions to be called when there is a need of Limiting the slopes
# data= dofs of particular cell
# k = varible number . which column in dofs to be modify
# basis =basis of a problem , node informaton
# a = a in linear eqaution, of that particular k
# b = b in linear eqaution, of that particular k
# c = c in linear eqaution, of that particular k
function update_solution(data, k, basis, a, b, c, cell)

    linearindex = vec(collect(Iterators.product(1:basis.order, 1:basis.order))) #linear indexing
    cellcenter = cell.center

    for i in 1:length(basis)
        x, y = globalposition(cell, (basis.quadpoints[linearindex[i][1]], basis.quadpoints[linearindex[i][2]]))
        data[i,k] = a + b * ( x - cellcenter[1] ) + c * ( y - cellcenter[2] )
    end

end

# Caluations required to decide limiting
function minmod_limiter(eq, scenario, dofs, grid, basis, ndofs, ndims)
    # Travelling each cell in a grid
    for i in eachindex(grid.cells)
        @views cell = grid.cells[i]
        facetypes = cell.facetypes
        @views data = dofs[:,:, cell.dataidx]
        weights = kron(basis.quadweights, basis.quadweights)

        # calculation of cell (weighted) averages of own cell and its neighbouring cells
        # calculation of a
        # a = ∑_n w_n u_n
        uown = vec(data' * weights)

        # For all neighbouring cells
        u_mean_neighbour = zeros(4,ndofs)
        for (i, neigh) in enumerate(cell.neighbors)
            @views dofsneigh = dofs[:,:,neigh.dataidx]
            u_mean_neighbour[i,:] =  dofsneigh' * weights
        end

        needs_limiting = BitArray(undef, ndofs) #initialise boolean for whether this cell will use limiting
        s_limit = Array{Float64}(undef, ndofs, ndims)

        for j in 1:ndims
            # points for (x - 1/2) or (y - 1/2)
            if j == 1
                points = repeat(basis.quadpoints, length(basis.quadpoints))
            else
                points = vcat(fill.(basis.quadpoints, length(basis.quadpoints))...)
            end

            #for x-dim:  s1 = 12* ∑_n w_n (x_n - 1/2) u_n
            s1 = 12 * vec(sum(weights .* (points .- 0.5) .* data, dims=1))
            s1 /= cell.size

            s2 = (uown - u_mean_neighbour[j, :]) / cell.size
            if facetypes[j] == 2  #at boundaries set slope to 0
                s2 = 0
            end
            s3 = (u_mean_neighbour[j + 2, :] - uown) / cell.size
            if facetypes[j + 2] == 2 #at boundaries set slope to 0
                s3 = 0
            end

            for k in 1:ndofs
                if s1[k]*s2[k] > 0 && s1[k]*s3[k] > 0  #if all have same sign
                    s_limit[k, j] = sign(s1[k])* min(abs(s1[k]),abs(s2[k]), abs(s3[k]))
                else
                    s_limit[k, j] = 0
                end

                if s_limit[k, j] != s1[k]
                    needs_limiting[k] = true
                end

            end

        end
        for k in 1:ndofs
            if needs_limiting[k]
                update_solution(data, k, basis, uown[k], s_limit[k,1], s_limit[k,2], cell)
            end
        end

    end
end