text:
struct SteadyStateSolution{T,N,uType,R,P,A} <: AbstractSteadyStateSolution{T,N}
  u::uType
  resid::R
  prob::P
  alg::A
  retcode::Symbol
end

function build_solution(prob::AbstractSteadyStateProblem,
                        alg,u,resid;calculate_error = true,
                        retcode = :Default, kwargs...)

  T = eltype(eltype(u))
  N = length((size(prob.u0)...,))

  SteadyStateSolution{T,N,typeof(u),typeof(resid),typeof(prob),typeof(alg)}(u,resid,prob,alg,retcode)
end


text:
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: September 2020
# Copyright: Georgia Institute of Technology, 2020

module FastApproximateInference

using JOLI
using JLD
using Random
using DataFrames
using LinearAlgebra
using Distributions
using Statistics
using InvertibleNetworks
using CUDA

# Utilities
include("./utils/load_experiment.jl")
include("./utils/data_loader.jl")
include("./utils/preprocessing.jl")
include("./utils/savefig.jl")
include("./utils/logpdf.jl")

# Sampling
include("./sampling/pSGLD.jl")

# Data generation
include("./modeling/forward_model.jl")
include("./modeling/rosenbrock.jl")

# Objective functions
include("./objectives/objectives.jl")

# Likelihoods with change of variable formula
include("./objectives/exact_likelihood.jl")

end

ext: jl


text:
const ARCH = @static if string(Sys.ARCH)[1] == 'i'
    "x86"
else
    string(Sys.ARCH)
end

const PLATFORM = @static if Sys.islinux()
    "linux-" * ARCH
elseif Sys.iswindows()
    "windows-" * ARCH
elseif Sys.isapple()
    "darwin-" * ARCH
end

OUTPUT = "$(PLATFORM)-v$(VERSION.major)_$(VERSION.minor)_$(VERSION.patch).jl"
OUTPUT_DEST = joinpath("src/memlayout", OUTPUT)

isfile(OUTPUT_DEST) && exit(0)

const PROJECT_DIR = (@__DIR__) |> dirname |> dirname |> dirname
const INCLUDE = joinpath(dirname(Sys.BINDIR), "include/julia")
const OPTIONS = if string(Sys.ARCH)[1] == 'i'
    Sys.islinux() && run(`sudo apt-get install g++-multilib -y`)
    ["-std=c++11","-march=pentium4", "-m32", "-static-libgcc", "-static-libstdc++"]
else
    ["-std=c++11"]
end

run(`g++ $(OPTIONS) -I$(INCLUDE) $(PROJECT_DIR)/.github/workflows/tasklayout/memlayout.cpp -o memlayout-gen.exe`)
run(`./memlayout-gen.exe`)

isfile(OUTPUT) && Base.Filesystem.mv(OUTPUT, OUTPUT_DEST)


text:
# Laplacian

struct Laplacian{T,D} <: LazyQuasiMatrix{T}
    axis::Inclusion{T,D}
end

Laplacian{T}(axis::Inclusion{<:Any,D}) where {T,D} = Laplacian{T,D}(axis)
# Laplacian{T}(domain) where T = Laplacian{T}(Inclusion(domain))
# Laplacian(axis) = Laplacian{eltype(axis)}(axis)

axes(D::Laplacian) = (D.axis, D.axis)
==(a::Laplacian, b::Laplacian) = a.axis == b.axis
copy(D::Laplacian) = Laplacian(copy(D.axis))

@simplify function *(Δ::Laplacian, P::AbstractSphericalHarmonic)
     # Spherical harmonics are the eigenfunctions of the Laplace operator on the unit sphere
     P * Diagonal(mortar(Fill.((-(0:∞)-(0:∞).^2), 1:2:∞)))
end

# Negative fractional Laplacian (-Δ)^α or equiv. abs(Δ)^α

struct AbsLaplacianPower{T,D,A} <: LazyQuasiMatrix{T}
    axis::Inclusion{T,D}
    α::A
end

AbsLaplacianPower{T}(axis::Inclusion{<:Any,D},α) where {T,D} = AbsLaplacianPower{T,D,typeof(α)}(axis,α)

axes(D:: AbsLaplacianPower) = (D.axis, D.axis)
==(a:: AbsLaplacianPower, b:: AbsLaplacianPower) = a.axis == b.axis && a.α == b.α
copy(D:: AbsLaplacianPower) = AbsLaplacianPower(copy(D.axis), D.α)

abs(Δ::Laplacian) = AbsLaplacianPower(axes(Δ,1),1)
-(Δ::Laplacian) = abs(Δ)

^(D::AbsLaplacianPower, k) = AbsLaplacianPower(D.axis, D.α*k)


text:
using LinearAlgebra, Plots, JSON, Jacobi
include("../src/SFFM.jl")
include("METools.jl")

include("stochasticApproximations.jl")
p = plot!()

T = [-2.0 2.0; 1.0 -1.0]#[-2.0 2.0 0; 0.5 -2.0 1.5; 0 2 -2]#
C = [1.0; -2.0]#; -1]
fn(x) = [ones(size(x)) ones(size(x))]# ones(size(x))]
model = SFFM.Model(;T=T,C=C,r=(r=fn,R=fn),Bounds=[0 10;-Inf Inf])
N₋ = sum(C.<=0)
N₊ = sum(C.>=0)
NPhases = length(C)

tim = 3
# τ = SFFM.FixedTime(T=tim)
# NSim = 500_000
# sims = SFFM.SimSFM(model=model,StoppingTime=τ,InitCondition=(φ=2*ones(Int,NSim),X=zeros(NSim)))

function e(i,order)
    e = zeros(order)
    e[i] = 1
    return e
end

let
    globalerrME = []
    globalerrDG = []
    orders = 1:2:23
    for order in orders
        μ = vecΔ
        ME = MakeME(CMEParams[order], mean = μ)#MakeErlang(order, mean = μ)#
        D = CMEParams[order]["D"]

        # if order==1
        #     tvec = [0]
        #     midpoints = 0
        # elseif order==2
        #     tvec=[0;1]
        # else
        #     h = (1)./(order-2)
        #     # tvec = [0;range(0+h/2,1-h/2,length=order-2);1]
        #     # tvec = range(0,2*(π./ME.Q[3,2]),length=order+1)[1:end-1]
        #     # tvec = ((-cos.(π.*range(0,1,length=order+1)).+1.0)./2*μ)[1:end-1]
        #     tvec = range(0,μ,length=order)[1:end]
        #     # tvec = ((Jacobi.zglj(order+1, 0, 0).+1.0)./2*μ)[1:end-1]#2*(π./ME.Q[3,2]))[1:end-1]
        #     # tvec4 = [0;range(0+h/2,1-h/2,length=order-1)]
        #     midpoints = [0;(tvec[1:end-1]+tvec[2:end])./2]
        #     # midpoints = range(0,μ,length=order)[1:end]
        # end
        # # tvec = range(0,1,length=order+1)
        # # tvec = tvec[1:end-1]
        # #[0;range(0.01,20,length=order-1)]
        # #[0;0.01;20]
        # #[0;range(0.1*μ,2*μ,length=order-1)]
        # #2*((exp(1).^(range(0,1,length=order))).-1)./(exp(1)-1)
        # H = zeros(order,order)
        # # H[1,:] = ME.α
        # for n in 1:order
        #     u = ME.α*exp(ME.Q*tvec[n])
        #     # u = ME.α*(ME.Q^-1)*(exp(ME.Q*tvec[n]*μ)-exp(ME.Q*tvec[n-1]*μ))
        #     u = u./sum(u)
        #     H[n,:] = u
        # end
        # # H[end,:] = ME.α*(ME.Q^-1)*[-exp(ME.Q*tvec[end]*μ)]
        #
        # MEf = (α = [1;zeros(order-1)]', Q = H*ME.Q*H^-1)
        # # MEf2 = (α = [0;1;zeros(order-2)]', Q = H*ME.Q*H^-1)
        #
        # # timesvector = 0:0.01:2*μ
        # # let
        # #     f = zeros(length(timesvector),order)
        # #     g = zeros(length(timesvector),order)
        # #     c = 0
        # #     for t in timesvector
        # #         c = c+1
        # #         g[c,:] = H*exp(ME.Q*t)*ME.q
        # #     end
        # #     plot(timesvector,g)
        # #     scatter!(tvec,zeros(size(tvec)))
        # # end
        #
        # # D = ME.q.*D
        # # D = I(order)[end:-1:1,:]
        # if order>1
        #     idx = sum(tvec.<=μ)+1
        # else
        #     idx = order
        # end
        # uvec = midpoints#tvec#[1:idx]
        # D = zeros(order,order)
        # for n in 2:length(uvec)
        #     dt = uvec[n]-uvec[n-1]
        #     u = (I-exp(ME.Q*dt))*exp(ME.Q*uvec[n-1])*ones(order)
        #
        #     # if n==2
        #     #     dt = (tvec[n]-tvec[n-1])/2
        #     #     u = (I-exp(MEf.Q*μ*dt))*ones(order)
        #     # else
        #     #     dt = tvec[n]-tvec[n-1]
        #     #     u = (I-exp(MEf.Q*μ*dt))*exp(MEf.Q*μ*midpoints[n-2])*ones(order)
        #     # end
        #     D[:, length(uvec)-n+2] = u
        # end
        # u = exp(ME.Q*uvec[end])*ones(order)
        # # u = exp(MEf.Q*μ*midpoints[end])*ones(order)
        # D[:,1] = u
        # D = D[:,end:-1:1]*H
        # D = integrateD(100000,CMEParams[string(order)])
        #
        # U = zeros(order,order)
        # for n in 2:length(uvec)
        #     dt = uvec[n]-uvec[n-1]
        #     u = e(n,order)'*((I-exp(ME.Q*dt))*exp(ME.Q*uvec[n-1]))^-1*ME.Q
        #     u = u./sum(u)
        #
        #     # if n==2
        #     #     dt = (tvec[n]-tvec[n-1])/2
        #     #     u = (I-exp(MEf.Q*μ*dt))*ones(order)
        #     # else
        #     #     dt = tvec[n]-tvec[n-1]
        #     #     u = (I-exp(MEf.Q*μ*dt))*exp(MEf.Q*μ*midpoints[n-2])*ones(order)
        #     # end
        #     U[:, n-1] = u
        # end
        # u = e(1,order)'*exp(ME.Q*uvec[end])^-1*ME.Q
        # u = u./sum(u)
        # # u = exp(MEf.Q*μ*midpoints[end])*ones(order)
        # U[:,end] = u
        # MEU = (α = ME.α*U^-1, Q = U*ME.Q*U^-1)
        # D = I(order)[:,end:-1:1]

        function MakeGlobalApprox(;NCells = 3,up, down,T,C,bkwd=false,jumpMatrixD=I)
            αup,Qup = up
            αdown,Qdown = down
            N₋ = sum(C.<=0)
            N₊ = sum(C.>=0)
            NPhases = length(C)
            order = length(αup)
            qup = -sum(Qup,dims=2)
            qdown = -sum(Qdown,dims=2)
            Q = zeros(NCells*order*NPhases,NCells*order*NPhases)
            for n in 1:NCells
                for i in 1:NPhases
                    idx = (1:order) .+ (n-1)*(order*NPhases) .+ (i-1)*order
                    if C[i]>0
                        Q[idx,idx] = Qup*abs(C[i])
                        if n<NCells
                            Q[idx,idx .+ order*NPhases] = qup*αup*abs(C[i])
                        end
                    elseif C[i]<0
                        Q[idx,idx] = Qdown*abs(C[i])
                        if n>1
                            Q[idx,idx .- order*NPhases] = qdown*αdown*abs(C[i])
                        end
                    end
                end
            end
            T₋₋ = T[C.<=0,C.<=0]
            T₊₋ = T[C.>=0,:].*((C.<0)')
            T₋₊ = T[C.<=0,:].*((C.>0)')
            T₊₊ = T[C.>=0,C.>=0]

            inLower = [kron(diagm(abs.(C).*(C.<0)),qdown)[:,C.<0]; zeros((NCells-1)*NPhases*order,N₋)]
            outLower = [kron(1,kron(T₋₊,αup)) zeros(N₋,N₊+(NCells-1)*NPhases*order)]
            inUpper = [zeros((NCells-1)*NPhases*order,N₊);kron(diagm(abs.(C).*(C.>0)),qup)[:,C.>0]]
            outUpper = [zeros(N₊,N₋+(NCells-1)*NPhases*order) kron(1,kron(T₊₋,αdown))]

            if bkwd
                Tdiag = diagm(diag(T))
                Toff = T-diagm(diag(T))
                T₊₋ = T.*((C.>0)*(C.<0)')
                T₋₊ = T.*((C.<0)*(C.>0)')
                Tchange = T₊₋ + T₋₊
                Tnochange = T-Tchange

                D = zeros(size(Q))
                for c in 1:NCells, i in 1:length(C), j in 1:length(C)
                    idxi = ((i-1)*order) .+ (1:order) .+ ((c-1)*order*NPhases)
                    idxj = ((j-1)*order) .+ (1:order) .+ ((c-1)*order*NPhases)
                    if i!=j
                        if C[i]>0 && C[j]<0
                            D[idxi,idxj] = T[i,j].*jumpMatrixD
                        elseif C[i]<0 && C[j]>0
                            D[idxi,idxj] = T[i,j].*jumpMatrixD
                        else
                            D[idxi,idxj] = T[i,j].*I(order)
                        end
                    else
                        D[idxi,idxj] = T[i,i].*I(order)
                    end
                end
                B = [
                    T₋₋ outLower;
                    inLower D+Q inUpper;
                    outUpper T₊₊;
                ]
            else
                B = [
                    T₋₋ outLower;
                    inLower kron(I(NCells),kron(T,I(order)))+Q inUpper;
                    outUpper T₊₊;
                ]
            end
            return Q, B
        end

        Δ = μ
        Nodes = collect(0:Δ:10)
        NCells = length(Nodes)-1
        Q, B = MakeGlobalApprox(
            NCells = NCells,
            up = ME,
            down = ME,
            T = T,
            C = C,
            bkwd = true,
            jumpMatrixD = D,
        )
        if order<6
            display(B)
        end

        # DGMesh = SFFM.MakeMesh(model=model,NBases=1,Nodes=collect(Nodes[1]:Δ/order:Nodes[end]),Basis="lagrange")
        DGMesh = SFFM.MakeMesh(model=model,NBases=order,Nodes=Nodes,Basis="lagrange")
        All = SFFM.MakeAll(model=model,mesh=DGMesh)

        initDist = zeros(1,size(All.B.B,1))
        initDist[1] = 1

        temp = initDist*exp(Matrix(All.B.B)*tim)#SFFM.EulerDG(D=All.B.B,y=t,x0=initDist)#
        DGdist_t = SFFM.Coeffs2Dist(model=model,mesh=DGMesh,Coeffs=temp,type="probability")

        initDist = zeros(1,size(B,1))
        initDist[1] = 1
        dist_t = initDist*exp(B*tim)#SFFM.EulerDG(D=B,y=t,x0=initDist)#
        display(dist_t)
        display(DGdist_t.distribution)
        pm_t = dist_t[[1:N₋;(end-N₊+1):end]]
        dist_t = reshape(dist_t[N₋+1:end-N₊],order,NPhases,NCells)

        mesh = SFFM.MakeMesh(model=model,NBases=order,Nodes=Nodes,Basis="lagrange")
        simDist = SFFM.Sims2Dist(model=model,mesh=mesh,sims=sims,type="probability")

        localerrME = sum(abs.(pm_t-simDist.pm))
        localerrDG = sum(abs.(DGdist_t.pm-simDist.pm))

        # q = plot(layout = (2,1))
        for n in 1:NCells
            x = simDist.x[n]
            for i in 1:NPhases
                if n==1
                    label1 = "DG"
                    label2 = "ME"
                else
                    label1 = false
                    label2 = false
                end
                yvalsME = [sum(dist_t[:,i,n])]
                # yvalsDG = [sum(DGdist_t.distribution[:,(1:order).+(n-1)*order,i])]
                yvalsDG = [sum(DGdist_t.distribution[:,n,i])]
                # scatter!([x],yvalsDG,label=label1,subplot=i,color=:blue,markershape=:rtriangle)
                # scatter!([x],yvalsME,label=label2,subplot=i,color=:black,markershape=:ltriangle)
                localerrME += abs(sum(yvalsME-simDist.distribution[:,n,i]))
                localerrDG += abs(sum(yvalsDG-simDist.distribution[:,n,i]))
            end
        end
        # display(q)
        # for i in 1:NPhases
        #     scatter!(q,simDist.x,simDist.distribution[:,:,i][:],subplot=i,label="sim")
        # end
        # display(q)
        push!(globalerrME,localerrME)
        push!(globalerrDG,localerrDG)
    end
    p = plot!(orders,log.(globalerrME),label="CME",linestyle=:dash,colour=2)
    plot!(p,orders,log.(globalerrDG),label="DG NBases",linestyle=:solid,colour=3)
    display(p)
end
# push!(globalerrME,localerrME)
# push!(globalerrDG,localerrDG)
# MẼ.q
# #
# timesvector = 0:0.01:2*π./ME.Q[3,2]
# let
#     f = zeros(length(timesvector),order)
#     g = zeros(length(timesvector),order)
#     c = 0
#     for t in timesvector
#         c = c+1
#         g[c,:] = H*exp(ME.Q*t)*ME.q
#     end
#     plot(timesvector,g)
# end
# #
# #
# NB = (α = [1;zeros(order-1)]', Q = MEf.Q)
# plot!(timesvector,density(timesvector,NB))
#
#
# plot(timesvector, orbit(timesvector,MEf)*D[:,end:-1:1], legend=:outertopright)
# plot!(timesvector, density(timesvector,MEf))
# scatter!(tvec, 0*tvec)
#
# MEtemp = (α = orbit(uvec[1]+uvec[2]/2,MEf)*D, Q = MEf.Q)
# plot(timesvector, density(timesvector,MEtemp), legend=:outertopright)
# for i in 2:length(uvec)
#     MEtemp = (α = orbit(uvec[i]+uvec[2]/2,MEf)*D, Q = MEf.Q)
#     plot!(timesvector, density(timesvector,MEtemp), legend=:outertopright)
# end
#
# plot!()

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
module TestModels

  using SearchLight

  import SearchLight: storableFields

  ######## Model from Genie-Searchligth-example-app extracted ############
  export Book, BookWithInterns, Callback, BookWithAuthor, Author
  export addBook, addBooks, seed, seedBook, seedAuthor
  using SearchLight, Dates

  mutable struct Book <: AbstractModel

    ### FIELDS
    id::DbId
    title::String
    author::String
    cover::String

    ### VALIDATION
    # validator::ModelValidator

    ### CALLBACKS
    # before_save::Function
    # after_save::Function
    # on_save::Function
    # on_find::Function
    # after_find::Function

    ### SCOPES
    # scopes::Dict{Symbol,Vector{SearchLight.SQLWhereEntity}}

    ### constructor
    Book(;
      ### FIELDS
      id = DbId(),
      title = "",
      author = "",
      cover = "",

      ### VALIDATION
      # validator = ModelValidator([
      #   ValidationRule(:title, BooksValidator.not_empty)
      # ]),

      ### CALLBACKS
      # before_save = (m::Todo) -> begin
      #   @info "Before save"
      # end,
      # after_save = (m::Todo) -> begin
      #   @info "After save"
      # end,
      # on_save = (m::Todo, field::Symbol, value::Any) -> begin
      #   @info "On save"
      # end,
      # on_find = (m::Todo, field::Symbol, value::Any) -> begin
      #   @info "On find"
      # end,
      # after_find = (m::Todo) -> begin
      #   @info "After find"
      # end,

      ### SCOPES
      # scopes = Dict{Symbol,Vector{SearchLight.SQLWhereEntity}}()

    ) = new(id, title, author, cover                                                             ### FIELDS
            # validator,                                                                  ### VALIDATION
            # before_save, after_save, on_save, on_find, after_find                       ### CALLBACKS
            # scopes                                                                      ### SCOPES
            )
  end

  mutable struct BookWithInterns <: AbstractModel
    ### INTERNALS
    _table_name::String
    _id::String
    _serializable::Vector{Symbol}

    ### FIELDS
    id::DbId
    title::String
    author::String
    cover::String

    ### VALIDATION
    # validator::ModelValidator

    ### CALLBACKS
    # before_save::Function
    # after_save::Function
    # on_save::Function
    # on_find::Function
    # after_find::Function

    ### SCOPES
    # scopes::Dict{Symbol,Vector{SearchLight.SQLWhereEntity}}

    ### constructor
    BookWithInterns(;
      ### FIELDS
      id = DbId(),
      title = "",
      author = "",
      cover = "",

      ### VALIDATION
      # validator = ModelValidator([
      #   ValidationRule(:title, BooksValidator.not_empty)
      # ]),

      ### CALLBACKS
      # before_save = (m::Todo) -> begin
      #   @info "Before save"
      # end,
      # after_save = (m::Todo) -> begin
      #   @info "After save"
      # end,
      # on_save = (m::Todo, field::Symbol, value::Any) -> begin
      #   @info "On save"
      # end,
      # on_find = (m::Todo, field::Symbol, value::Any) -> begin
      #   @info "On find"
      # end,
      # after_find = (m::Todo) -> begin
      #   @info "After find"
      # end,

      ### SCOPES
      # scopes = Dict{Symbol,Vector{SearchLight.SQLWhereEntity}}()

    ) = new("bookwithinterns", "id", Symbol[],                                                ### INTERNALS
            id, title, author, cover                                                             ### FIELDS
            # validator,                                                                  ### VALIDATION
            # before_save, after_save, on_save, on_find, after_find                       ### CALLBACKS
            # scopes                                                                      ### SCOPES
            )
  end

  Base.@kwdef mutable struct Callback <: AbstractModel
    id::DbId = DbId()
    title::String = ""
    indicator::Bool = true
    created_at::String = string(Dates.now())
    # callbacks
    before_save::Function = (m::Callback) -> begin
      @info "Do something before saving"
    end
    after_save::Function = (m::Callback) -> begin
      @info "Do something after saving"
    end
  end

  function storableFields(m::Type{BookWithInterns})
    Dict(["id" => "id",
          "author" => "author",
          "title" => "title",
          "cover" => "cover"])
  end

  mutable struct BookWithAuthor<:AbstractModel
      id::DbId
      id_author::DbId
      title::String
      BookWithAuthor(;title="",id_author= DbId()) = new(DbId(),id_author,title)
  end

  mutable struct Author<:AbstractModel
      id::DbId
      first_name::String
      last_name::String
      books::Array{BookWithAuthor,1}
      Author(;firstname ="", lastname="", books=Array{BookWithAuthor,1}()) = new(DbId(),firstname,lastname,books)
  end

  function storableFields(m::Type{Author})
    Dict(["id" => "id",
     "first_name" => "firstname",
     "last_name" => "lastname",
     "books" => ""])
  end

  function Base.print(io::IO, book::BookWithAuthor)
      print(io, "title: $(book.title)")
  end

  function Base.show(io::IO, author::Author)
      print(io,"Firstname: $(author.first_name) Lastname: $(author.last_name) books: $(length(author.books))")
  end

  """
    function addBook(author::Author, book::BookWithAuthor)
      attach a book to an author
  """
  function addBook(author::Author, book::BookWithAuthor)
      book.id_author = author.id
      push!(author.books, book)
  end

  function addBooks(author::Author, books::Vector{BookWithAuthor})
      map(book->addBook(author, book), books)
  end

  function seed()
    BillGatesBooks = [
      ("The Best We Could Do", "Thi Bui"),
      ("Evicted: Poverty and Profit in the American City", "Matthew Desmond"),
      ("Believe Me: A Memoir of Love, Death, and Jazz Chickens", "Eddie Izzard"),
      ("The Sympathizer!", "Viet Thanh Nguyen"),
      ("Energy and Civilization, A History", "Vaclav Smil")
    ]
  end

  function seedBook()
    [ "A Time to Kill",
    "The Firm",
    "The Pelican Brief",
    "The Client",
    "The Chamber"]
  end

  function seedAuthor()
    ["MAO ZEDONG",
    "AGATHA CHRISTIE",
    "BARBARA CARTLAND",
    "CORÍN TELLADO",
    "CHARLES SCHULZ",
    "EIICHIRO ODA",
    "JAMES PATTERSON",
    "HORATIO ALGER",
    "LEO TOLSTOY"]
  end

end ### End Module


text:
module HalfIntegers

export
    # Abstract types
    HalfInteger,

    # Concrete types
    BigHalfInt,
    Half,
    HalfInt,
    HalfInt8,
    HalfInt16,
    HalfInt32,
    HalfInt64,
    HalfInt128,
    HalfUInt,
    HalfUInt8,
    HalfUInt16,
    HalfUInt32,
    HalfUInt64,
    HalfUInt128,

    # Functions
    half,
    ishalfinteger,
    onehalf,
    twice

"""
    HalfInteger <: Real

Abstract supertype for half-integers. Here, every number ``n/2`` where ``n\\in\\mathbb{Z}``
is considered a half-integer, regardless of whether ``n`` is even or odd.
"""
abstract type HalfInteger <: Real end

HalfInteger(x::Real) = HalfInt(x)
HalfInteger(x::Integer) = Half{typeof(x)}(x)
HalfInteger(x::HalfInteger) = x
HalfInteger(x::Rational{T}) where T = Half{T}(x)
HalfInteger(x::BigFloat) = BigHalfInt(x)

(T::Type{<:AbstractFloat})(x::HalfInteger) = T(float(x))
(T::Type{<:Integer})(x::HalfInteger) =
    isinteger(x) ? T(twice(x) >> 1) : throw(InexactError(Symbol(T), T, x))
(::Type{Bool})(x::HalfInteger) = invoke(Bool, Tuple{Real}, x)
@static if VERSION ≥ v"1.5.0-DEV.820"
    function (::Type{Rational})(x::HalfInteger)
        tx = twice(x)
        if isinteger(x)
            Base.unsafe_rational(tx >> 1, one(tx))
        else
            Base.unsafe_rational(tx, oftype(tx, 2))
        end
    end
    function (::Type{Rational{T}})(x::HalfInteger) where T
        tx = twice(x)
        if isinteger(x)
            Base.unsafe_rational(T, tx >> 1, one(tx))
        else
            Base.unsafe_rational(T, tx, oftype(tx,2))
        end
    end
else
    (T::Type{<:Rational})(x::HalfInteger) = (tx=twice(x); T(tx,oftype(tx,2)))
end

Base.ArithmeticStyle(::Type{<:HalfInteger}) = Base.ArithmeticWraps()

for op in (:<, :≤, :(==))
    @eval Base.$op(x::HalfInteger, y::HalfInteger) = $op(twice(x), twice(y))
    @eval Base.$op(x::HalfInteger, y::Integer) = ((pt,py)=promote(twice(x),y); $op(pt, twice(py)))
    @eval Base.$op(x::Integer, y::HalfInteger) = ((px,pt)=promote(x,twice(y)); $op(twice(px), pt))
    @eval Base.$op(x::HalfInteger, y::Rational) = $op(twice(x)*denominator(y), twice(numerator(y)))
    @eval Base.$op(x::Rational, y::HalfInteger) = $op(twice(numerator(x)), twice(y)*denominator(x))
end

Base.:+(x::T, y::T) where T<:HalfInteger = half(twice(x)+twice(y))
Base.:+(x::HalfInteger) = half(+twice(x))
Base.:-(x::T, y::T) where T<:HalfInteger = half(twice(x)-twice(y))
Base.:-(x::HalfInteger) = half(-twice(x))

Base.:*(x::T, y::T) where T<:HalfInteger = float(x)*float(y)
Base.:*(x::HalfInteger, y::Integer) = half(twice(x)*y)
Base.:*(x::Integer, y::HalfInteger) = y*x

Base.:/(x::T, y::T) where T<:HalfInteger = twice(x)/twice(y)

Base.://(x::HalfInteger, y::HalfInteger) = twice(x)//twice(y)
Base.://(x::HalfInteger, y::Real) = twice(x)//twice(y)
Base.://(x::Real, y::HalfInteger) = twice(x)//twice(y)

Base.:^(x::Real, y::HalfInteger) = x^float(y)
Base.:^(::Irrational{:ℯ}, x::HalfInteger) = exp(x)

@static if VERSION < v"1.4.0-DEV.208"
    Base.div(x::T, y::T) where T<:HalfInteger = div(twice(x), twice(y))
else
    Base.div(x::T, y::T, r::RoundingMode) where T<:HalfInteger = div(twice(x), twice(y), r)
end
Base.fld(x::T, y::T) where T<:HalfInteger = fld(twice(x), twice(y)) # can be removed in Julia 2.0
Base.cld(x::T, y::T) where T<:HalfInteger = cld(twice(x), twice(y)) # can be removed in Julia 2.0
Base.rem(x::T, y::T) where T<:HalfInteger = half(rem(twice(x), twice(y)))
Base.mod(x::T, y::T) where T<:HalfInteger = half(mod(twice(x), twice(y)))

Base.fld1(x::T, y::T) where T<:HalfInteger = fld1(twice(x), twice(y))

const HalfIntegerOrInteger = Union{HalfInteger, Integer}

# Checked arithmetic
Base.checked_abs(x::HalfInteger) = half(Base.checked_abs(twice(x)))
Base.checked_neg(x::HalfInteger) = half(Base.checked_neg(twice(x)))

# These methods are only called if at least one of the arguments is a `HalfInteger`,
# otherwise the `Base.checked_$op(::Integer, ::Integer)` methods from `Base` are used
for op in (:checked_add, :checked_sub, :checked_mod, :checked_rem)
    @eval Base.$op(x::HalfIntegerOrInteger, y::HalfIntegerOrInteger) =
        half(Base.$op(checked_twice(x), checked_twice(y)))
end
for op in (:checked_cld, :checked_div, :checked_fld)
    @eval Base.$op(x::HalfIntegerOrInteger, y::HalfIntegerOrInteger) =
        Base.$op(checked_twice(x), checked_twice(y))
end

Base.checked_mul(x::HalfInteger, y::HalfInteger) = x*y # uses float arithmetic
Base.checked_mul(x::HalfInteger, y::Integer) = half(Base.checked_mul(twice(x), y))
Base.checked_mul(x::Integer, y::HalfInteger) = Base.checked_mul(y, x)

for f in (:add_with_overflow, :sub_with_overflow)
    @eval function Base.Checked.$f(x::HalfInteger, y::HalfInteger)
        r, flag = Base.Checked.$f(twice(x), twice(y))
        half(r), flag
    end
end

Base.Checked.mul_with_overflow(x::HalfInteger, y::HalfInteger) = (x*y, false) # uses float arithmetic
function Base.Checked.mul_with_overflow(x::HalfInteger, y::Integer)
    r, flag = Base.Checked.mul_with_overflow(twice(x), y)
    half(r), flag
end
Base.Checked.mul_with_overflow(x::Integer, y::HalfInteger) = Base.Checked.mul_with_overflow(y, x)

# `lcm`/`gcd`/`gcdx` are only extended to `HalfInteger`s if they are defined for `Rational`s
@static if VERSION ≥ v"1.4.0-DEV.699"
    Base.gcd(x::HalfInteger) = x
    Base.lcm(x::HalfInteger) = x

    Base.gcd(a::T, b::T) where T<:HalfInteger = half(gcd(twice(a), twice(b)))
    Base.lcm(a::T, b::T) where T<:HalfInteger = half(lcm(twice(a), twice(b)))

    function Base.gcdx(x::T, y::T) where T<:HalfInteger
        d, u, v = gcdx(twice(x), twice(y))
        half(d), u, v
    end
end

Base.ceil(T::Type{<:Integer}, x::HalfInteger)  = round(T, x, RoundUp)
Base.floor(T::Type{<:Integer}, x::HalfInteger) = round(T, x, RoundDown)
Base.trunc(T::Type{<:Integer}, x::HalfInteger) = round(T, x, RoundToZero)

Base.round(T::Type{<:Integer}, x::HalfInteger, r::RoundingMode=RoundNearest) =
    T(twice(round(x, r)) >> 1)
Base.round(T::Type{<:Integer}, x::HalfInteger, ::typeof(RoundDown)) = T(twice(x) >> 1)

Base.round(x::HalfInteger, ::typeof(RoundNearest)) =
    isinteger(x)             ? +x           :
    isone(mod1(twice(x), 4)) ? x-onehalf(x) : x+onehalf(x)
Base.round(x::HalfInteger, ::typeof(RoundNearestTiesAway)) =
    isinteger(x) ? +x           :
    x < zero(x)  ? x-onehalf(x) : x+onehalf(x)
Base.round(x::HalfInteger, ::typeof(RoundNearestTiesUp)) =
    isinteger(x) ? +x : x+onehalf(x)
Base.round(x::HalfInteger, ::typeof(RoundToZero)) =
    isinteger(x) ? +x           :
    x < zero(x)  ? x+onehalf(x) : x-onehalf(x)
Base.round(x::HalfInteger, ::typeof(RoundUp)) =
    round(x, RoundNearestTiesUp)
Base.round(x::HalfInteger, ::typeof(RoundDown)) =
    isinteger(x) ? +x : x-onehalf(x)

Base.big(x::HalfInteger) = BigHalfInt(x)

Base.float(x::HalfInteger) = twice(x)/2

Base.isfinite(x::HalfInteger) = isfinite(twice(x))

Base.isinteger(x::HalfInteger) = iseven(twice(x))

@static if VERSION ≥ v"1.6.0-DEV.999"
    Base.ispow2(x::HalfInteger) = ispow2(twice(x))
end

Base.show(io::IO, x::HalfInteger) =
    isinteger(x) ? print(io, twice(x) >> 1) : print(io, twice(x), "/2")

Base.sign(x::HalfInteger) = oftype(x, sign(twice(x)))

Base.denominator(x::HalfInteger) = (t=twice(x); isinteger(x) ? one(t) : oftype(t,2))
Base.numerator(x::HalfInteger)   = (t=twice(x); isinteger(x) ? t >> 1 : t)

Base.sinpi(x::HalfInteger) = sinpihalf(twice(x))
Base.cospi(x::HalfInteger) = cospihalf(twice(x))

Base.sinc(x::HalfInteger) = sinchalf(twice(x))
Base.cosc(x::HalfInteger) = coschalf(twice(x))

# Compute sin(x*π/2)
function sinpihalf(x::Integer)
    xm4 = mod(x,4)
    T = float(typeof(x))
    if (xm4 == 0) | (xm4 == 2)
        x < 0 ? -zero(T) : +zero(T)
    elseif xm4 == 1
        +one(T)
    else
        -one(T)
    end
end
# Compute cos(x*π/2)
function cospihalf(x::Integer)
    xm4 = mod(x,4)
    T = float(typeof(x))
    if (xm4 == 1) | (xm4 == 3)
        +zero(T)
    elseif xm4 == 0
        +one(T)
    else
        -one(T)
    end
end

# Compute sinc(x/2)
function sinchalf(x::Integer)
    iszero(x) && return one(float(typeof(x)))
    iseven(x) && return zero(float(typeof(x)))
    d = inv(π*x)
    isone(mod(x % Int8, Int8(4))) ? 2*d : -2*d
end
# Compute cosc(x/2)
function coschalf(x::Integer)
    T = float(typeof(x))
    iszero(x) && return zero(T)
    xm4 = mod(x % Int8, Int8(4))
    x⁻¹ = inv(x)
    if xm4 == 0
        2*x⁻¹
    elseif xm4 == 1
        (-4/T(π))*x⁻¹*x⁻¹
    elseif xm4 == 2
        -2*x⁻¹
    else
        (4/T(π))*x⁻¹*x⁻¹
    end
end

@static if VERSION ≥ v"1.6.0-DEV.292"
    Base.sincospi(x::HalfInteger) = sincospihalf(twice(x))

    # Compute sincos(x*π/2)
    function sincospihalf(x::Integer)
        xm4 = mod(x,4)
        T = float(typeof(x))
        if xm4 == 0
            (x < 0 ? -zero(T) : +zero(T), +one(T))
        elseif xm4 == 1
            (+one(T), +zero(T))
        elseif xm4 == 2
            (x < 0 ? -zero(T) : +zero(T), -one(T))
        else
            (-one(T), +zero(T))
        end
    end
end

Base.all(::typeof(isinteger), r::AbstractUnitRange{<:HalfInteger}) = !any(!isinteger, r)
Base.all(::typeof(!isinteger), r::AbstractUnitRange{<:HalfInteger}) = !any(isinteger, r)

Base.any(::typeof(isinteger), r::AbstractUnitRange{<:HalfInteger}) =
    !isempty(r) & isinteger(first(r))
Base.any(::typeof(!isinteger), r::AbstractUnitRange{<:HalfInteger}) =
    !isempty(r) & !isinteger(first(r))

# Hashing

Base.decompose(x::HalfInteger) = twice(x), 0, 2

Base.hash(x::HalfInteger, h::UInt) = hashhalf(twice(x), h)

# Compute hash(half(x), h)
hashhalf(x::Integer, h::UInt) = invoke(hash, Tuple{Real,UInt}, half(x), h)

# Version for integers with ≤ 64 bits, adapted from hash(::Rational{<:Base.BitInteger64}, ::UInt)
function hashhalf(x::Base.BitInteger64, h::UInt)
    iseven(x) && return hash(x >> 1, h)
    if abs(x) < 9007199254740992
        return hash(ldexp(Float64(x),-1),h)
    end
    h = Base.hash_integer(1, h)
    h = Base.hash_integer(-1, h)
    h = Base.hash_integer(x, h)
    return h
end

"""
    half([T<:HalfInteger,] x)

For an integer value `x`, return `x/2` as a half-integer of type `T`. If `T` is not given,
return an appropriate `HalfInteger` type. Throw an `InexactError` if `x` is not an integer.

!!! compat "HalfIntegers 1.1"
    `missing` as an argument requires at least HalfIntegers 1.1.

# Examples

```jldoctest
julia> half(3)
3/2

julia> half(-5.0)
-5/2

julia> half(HalfInt16, 8)
4
```
"""
half(x::Real) = half(HalfInteger, x)
half(::Type{HalfInteger}, x::Number) = half(HalfInteger, Integer(x))
half(::Type{HalfInteger}, x::Integer) = half(Half{typeof(x)}, x)

"""
    half([T<:Complex{<:HalfInteger},] x)

For a complex value `x` with integer real and imaginary parts, return `x/2` as a complex
number of type `T`. If `T` is not given, return an appropriate complex number type with
half-integer parts. Throw an `InexactError` if the real or the imaginary part of `x` are not
integer.

!!! compat "HalfIntegers 1.1"
    `missing` as an argument requires at least HalfIntegers 1.1.

# Examples

```jldoctest
julia> half(3 + 2im)
3/2 + 1*im

julia> half(Complex{HalfInt32}, 1.0 + 5.0im)
1/2 + 5/2*im
```
"""
half(x::Complex) = half(Complex, x)
half(::Type{Complex{T}}, x::Number) where T<:HalfInteger = Complex(half(T,real(x)), half(T,imag(x)))
half(::Type{Complex}, x::Number) = Complex(half(real(x)), half(imag(x)))

half(::Missing) = missing
@static if VERSION ≥ v"1.3"
    half(::Type{T}, x) where T>:Missing = half(Base.nonmissingtype_checked(T), x)
else
    half(::Type{T}, x) where T>:Missing = half(Base.nonmissingtype(T), x)
end
half(::Type{>:Missing}, ::Missing) = missing
half(::Type{T}, ::Missing) where T =
    throw(MissingException("cannot convert a missing value to type $T: use Union{$T, Missing} instead"))

"""
    twice(x)

Return `2x`. If `x` is a `HalfInteger`, return an appropriate integer type.

!!! compat "HalfIntegers 1.1"
    `missing` as an argument requires at least HalfIntegers 1.1.

# Examples

```jldoctest
julia> twice(2)
4

julia> twice(1.5)
3.0

julia> twice(HalfInt8(3/2))  # returns an Int8
3

julia> twice(HalfInt32(3.0) + HalfInt32(2.5)*im)  # returns a Complex{Int32}
6 + 5im
```
"""
twice(x) = x + x
twice(x::Complex) = Complex(twice(real(x)), twice(imag(x)))

@static if VERSION ≥ v"1.5.0-DEV.820"
    function twice(x::Rational)
        if iseven(x.den)
            Base.unsafe_rational(oftype(x.num+x.num, x.num), x.den >> 1)
        else
            Base.unsafe_rational(Base.checked_add(x.num, x.num), oftype(x.den >> 1, x.den))
        end
    end
    function twice(x::Rational{BigInt})
        if iseven(x.den)
            Base.unsafe_rational(x.num, x.den >> 1)
        else
            Base.unsafe_rational(twice(x.num), x.den)
        end
    end
else
    function twice(x::Rational)
        if iseven(x.den)
            Rational(oftype(x.num+x.num, x.num), x.den >> 1)
        else
            Rational(Base.checked_add(x.num, x.num), oftype(x.den >> 1, x.den))
        end
    end
    function twice(x::Rational{BigInt})
        if iseven(x.den)
            Rational(x.num, x.den >> 1)
        else
            Rational(twice(x.num), x.den)
        end
    end
end

"""
    twice(T<:Integer, x)
    twice(T<:Complex{<:Integer}, x)

Return `2x` as a number of type `T`.

!!! compat "HalfIntegers 1.1"
    `missing` as an argument requires at least HalfIntegers 1.1.

# Examples

```jldoctest
julia> twice(Int16, HalfInt(5/2))
5

julia> twice(Complex{BigInt}, HalfInt(5/2) + HalfInt(3)*im)
5 + 6im
```
"""
twice(T::Type{<:Integer}, x::Number) = T(twice(x))
twice(T::Type{<:Integer}, x::Integer) = twice(T(x)) # convert to T first to reduce probability of overflow
twice(::Type{Complex{T}}, x::Number) where T<:Integer = Complex(twice(T,real(x)), twice(T,imag(x)))

function twice(T::Type{<:Integer}, x::Rational)
    if x.den == 1
        twice(T, x.num)
    elseif x.den == 2
        convert(T, x.num)
    else
        throw(InexactError(nameof(T), T, twice(x)))
    end
end

@static if VERSION ≥ v"1.3"
    twice(::Type{T}, x) where T>:Missing = twice(Base.nonmissingtype_checked(T), x)
else
    twice(::Type{T}, x) where T>:Missing = twice(Base.nonmissingtype(T), x)
end
twice(::Type{>:Missing}, ::Missing) = missing
twice(::Type{T}, ::Missing) where T =
    throw(MissingException("cannot convert a missing value to type $T: use Union{$T, Missing} instead"))

"""
    checked_twice(x)

Calculate `2x`, checking for overflow. (not exported)

!!! compat "HalfIntegers 1.4"
    This function requires at least HalfIntegers 1.4.

# Examples

```jldoctest
julia> x = typemax(Int64)
9223372036854775807

julia> twice(x)
-2

julia> HalfIntegers.checked_twice(x)
ERROR: OverflowError: 9223372036854775807 + 9223372036854775807 overflowed for type Int64
Stacktrace:
[...]
```
"""
checked_twice(x) = _checked_twice(Base.ArithmeticStyle(x), x)
checked_twice(x::HalfInteger) = twice(x)

_checked_twice(::Base.ArithmeticStyle, x) = twice(x)
_checked_twice(::Base.ArithmeticWraps, x) = Base.checked_add(x, x)

"""
    onehalf(x)

Return the value 1/2 in the type of `x` (`x` can also specify the type itself).

!!! compat "HalfIntegers 1.1"
    `missing` as an argument requires at least HalfIntegers 1.1.

!!! compat "HalfIntegers 1.2"
    Some abstract types as arguments require at least HalfIntegers 1.2.

# Examples

```jldoctest
julia> onehalf(7//3)
1//2

julia> onehalf(HalfInt)
1/2
```
"""
onehalf(x::Union{Number,Missing}) = onehalf(typeof(x))

onehalf(T::Type{<:Number}) = T(onehalf(HalfInteger))
onehalf(T::Type{<:HalfInteger}) = half(T, 1)
onehalf(T::Type{<:AbstractFloat}) = T(0.5)
@static if VERSION ≥ v"1.5.0-DEV.820"
    onehalf(::Type{Rational}) = Base.unsafe_rational(1, 2)
    onehalf(::Type{Rational{T}}) where T = Base.unsafe_rational(T, 1, 2)
else
    onehalf(T::Type{<:Rational}) = T(1, 2)
end
onehalf(::Type{Complex{T}}) where T = Complex{T}(onehalf(T), zero(T))
onehalf(::Type{Missing}) = missing

"""
    ishalfinteger(x)

Test whether `x` is numerically equal to some half-integer.

# Examples

```jldoctest
julia> ishalfinteger(3.5)
true

julia> ishalfinteger(2)
true

julia> ishalfinteger(1//3)
false
```
"""
ishalfinteger(x) = isinteger(twice(x))
ishalfinteger(x::Rational) = (denominator(x) == 1) | (denominator(x) == 2)
ishalfinteger(::HalfIntegerOrInteger) = true
ishalfinteger(::AbstractIrrational) = false
ishalfinteger(::Missing) = missing

@static if VERSION ≥ v"1.7.0-DEV.263"
    Base.range_start_stop_length(start::T, stop::T, len::Integer) where T<:HalfInteger =
        Base._linspace(float(T), start, stop, len)
else
    Base._range(start::T, ::Nothing, stop::T, len::Integer) where T<:HalfInteger =
        Base._linspace(float(T), start, stop, len)
end

Base._linspace(::Type{T}, start::HalfInteger, stop::HalfInteger, len::Integer) where T =
    LinRange{T}(start, stop, len)
Base._linspace(::Type{T}, start::HalfInteger, stop::HalfInteger, len::Integer) where T<:Base.IEEEFloat =
    Base._linspace(T, twice(start), twice(stop), len, 2)

Base.in(x::Real, r::AbstractUnitRange{<:HalfInteger}) =
    ishalfinteger(x) & (first(r) ≤ x ≤ last(r)) & !(isinteger(x) ⊻ isinteger(first(r)))

function Base.intersect(r::AbstractUnitRange{<:HalfInteger},
                        s::AbstractUnitRange{<:HalfIntegerOrInteger})
    fir = max(first(r),first(s))
    las = min(last(r), last(s))
    las = isinteger(first(r)) ⊻ isinteger(first(s)) ? oftype(las, fir-one(fir)) : las
    fir:las
end
Base.intersect(r::AbstractUnitRange{<:Integer}, s::AbstractUnitRange{<:HalfInteger}) = s ∩ r

function Base.intersect(r::Union{AbstractUnitRange{<:HalfInteger}, StepRange{<:HalfInteger}},
                        s::Union{AbstractUnitRange{<:HalfIntegerOrInteger},
                                 StepRange{<:HalfIntegerOrInteger}})
    # All integers are promoted to a common type T to reduce the risk of overflow
    T = promote_type(typeof(twice(first(r))), typeof(twice(step(r))),
                     typeof(twice(first(s))), typeof(twice(step(s))))
    tr = twice(T,first(r)):twice(T,step(r)):twice(T,last(r))
    ts = twice(T,first(s)):twice(T,step(s)):twice(T,last(s))
    trange = tr ∩ ts
    half(first(trange)):half(step(trange)):half(last(trange))
end
Base.intersect(r::Union{AbstractUnitRange{<:Integer}, StepRange{<:Integer}},
               s::Union{AbstractUnitRange{<:HalfInteger}, StepRange{<:HalfInteger}}) = s ∩ r

"""
    Half{T<:Integer} <: HalfInteger

Type for half-integers ``n/2`` where ``n`` is of type `T`.

Aliases for `Half{T}` are defined for all standard `Signed` and `Unsigned` integer types, so
that, e.g., `HalfInt` can be used instead of `Half{Int}`:

| `T`       | Alias for `Half{T}`                |
| :-------- | :--------------------------------- |
| `Int`     | `HalfInt`                          |
| `Int8`    | `HalfInt8`                         |
| `Int16`   | `HalfInt16`                        |
| `Int32`   | `HalfInt32`                        |
| `Int64`   | `HalfInt64`                        |
| `Int128`  | `HalfInt128`                       |
| `BigInt`  | `BigHalfInt` (*not* `HalfBigInt`!) |
| `UInt`    | `HalfUInt`                         |
| `UInt8`   | `HalfUInt8`                        |
| `UInt16`  | `HalfUInt16`                       |
| `UInt32`  | `HalfUInt32`                       |
| `UInt64`  | `HalfUInt64`                       |
| `UInt128` | `HalfUInt128`                      |
"""
struct Half{T<:Integer} <: HalfInteger
    twice::T

    global half(::Type{Half{T}}, x::Number) where T<:Integer = new{T}(x)
end
Half{T}(x::Real) where T<:Integer = half(Half{T}, twice(T,x))
Half{T}(x::Half{T}) where T<:Integer = x

twice(x::Half) = x.twice

Base.promote_rule(::Type{Half{H}}, T::Type{<:Real}) where H<:Integer = promote_type(Rational{H}, T)
Base.promote_rule(::Type{Half{H}}, T::Type{<:Integer}) where H<:Integer = Half{promote_type(H, T)}
Base.promote_rule(::Type{Half{T}}, ::Type{Half{S}}) where {T<:Integer,S<:Integer} = Half{promote_type(T,S)}

function Base.parse(T::Type{<:Half}, s::AbstractString)
    result = tryparse(T, s)
    result === nothing && throw(ArgumentError("cannot parse $(repr(s)) as $T."))
    result
end

function Base.tryparse(::Type{Half{T}}, s::AbstractString) where T<:Integer
    matched = match(r"^\s*((?:\+|-|)\s*\d+)\s*(/\s*2|)\s*$", s)
    matched === nothing && return nothing
    num = tryparse(T, matched.captures[1])
    num === nothing && return nothing
    return isempty(matched.captures[2]) ? Half{T}(num) : half(num)
end

Base.hastypemax(::Type{Half{T}}) where T<:Integer = Base.hastypemax(T)
Base.typemax(::Type{Half{T}}) where T<:Integer = half(typemax(T))
Base.typemin(::Type{Half{T}}) where T<:Integer = half(typemin(T))

Base.widen(::Type{Half{T}}) where T<:Integer = Half{widen(T)}

const HalfInt = Half{Int}
const HalfInt8 = Half{Int8}
const HalfInt16 = Half{Int16}
const HalfInt32 = Half{Int32}
const HalfInt64 = Half{Int64}
const HalfInt128 = Half{Int128}
const BigHalfInt = Half{BigInt}
const HalfUInt = Half{UInt}
const HalfUInt8 = Half{UInt8}
const HalfUInt16 = Half{UInt16}
const HalfUInt32 = Half{UInt32}
const HalfUInt64 = Half{UInt64}
const HalfUInt128 = Half{UInt128}

Base.float(::Type{BigHalfInt}) = BigFloat

end # module


text:
push!(LOAD_PATH, "/home/ian/Grad/Research/Disks/DiskJockey/")
push!(LOAD_PATH, "/n/home07/iczekala/DiskJockey/")

# Using the analytic form of the FT of the Gaussian, and the u,v sampling and noise
# from a real SMA dataset, make a fake dataset.

using gridding
using visibilities
using constants
using image
using gauss_model
using constants

# Read in the dataset
#dv = DataVis("../data/V4046Sgr.hdf5", 12)
dvarr = DataVis("../data/V4046Sgr.hdf5")
nvis = length(dvarr[1].VV)
nlam = 1 #length(dvarr)

# Realistic Gaussian will have scale dimensions
const mu_RA = 2. # [arcsec]
const mu_DEC = -0.5 # [arcsec]
const s_x = 1.2 # [arcsec]
const s_y = 1.0 # [arcsec]
const p0 = [mu_RA, mu_DEC, s_x, s_y] # [arcsec]

# Using the analytic formula for the FT of the Gaussian, compute the true
# visibilities sampled at the u,v points in the SMA dataset
dvarr_fake = Array(DataVis, nlam)

const scale = 0.5
for k=1:nlam
    dv = dvarr[k]
    # Make a new array for each DataVis, otherwise references will be kept and duplicated
    VV_fake = Array(Complex128, nvis)
    invsig = Array(Float64, nvis)
    for i=1:nvis
        model = FTGauss(dv.uu[i], dv.vv[i], p0, k)
        noise = scale * (randn() + randn()*im) # Just adding in some noise here
        VV_fake[i] = model + noise
        invsig[i] = 1./scale
    end
    dvarr_fake[k] = DataVis(dv.lam, dv.uu, dv.vv, VV_fake, invsig)
end

# Make the fake dataset using these Gaussian visibilities and noise terms

# # Instead of using the u,v points in the SMA dataset, actually sample somewhat
# # more uniformly on a grid.
# nfake = 40
# u = linspace(-150, 150, nfake)
# v = linspace(-150, 150, nfake)
# uu = Array(Float64, nfake^2)
# vv = Array(Float64, nfake^2)
# VV_fake = Array(Complex128, nfake^2)
# invsig = Array(Float64, nfake^2)
# const scale = 0.5
# for i=1:nfake
#     for j=1:nfake
#         k = (i -1) * nfake + j
#         uu[k] = u[i]
#         vv[k] = v[j]
#         model = FTGauss(u[i], v[j], p0)
#
#         noise = scale * (randn() + randn()*im) # Just adding in some noise here
#         VV_fake[k] = model + noise
#         invsig[k] = 1./scale
#     end
# end
#
# # Make the fake dataset using these Gaussian visibilities and noise terms
# dv_fake = DataVis(dv.lam, uu, vv, VV_fake, invsig)

visibilities.write(dvarr_fake, "../data/V4046Sgr_fake.hdf5")

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
export iauObl06
"""
Mean obliquity of the ecliptic, IAU 2006 precession model.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  canonical model.

Given:
   date1,date2  double   TT as a 2-part Julian Date (Note 1)

Returned (function value):
               double   obliquity of the ecliptic (radians, Note 2)

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The result is the angle between the ecliptic and mean equator of
   date date1+date2.

Reference:

   Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauObl06(double date1, double date2)

function iauObl06(date1::Real, date2::Real)
   return ccall((:iauObl06, libsofa_c), Cdouble,
        (Cdouble, Cdouble),
        convert(Float64, date1), convert(Float64, date2))
end


text:
function bench()
  @time length(Primes(100)) # warm up JIT
#  println(@time count(x->true, Primes(1000000000))) # about 1.5 seconds slower counting over iteration
  println(@time length(Primes(1000000000)))
end
bench()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module TestExperimental
using Accessors
using PhaseSpaceIO
using Test
using PhaseSpaceIO.Testing

@testset "Collect" begin
    c = Collect(5)
    iter = randn(10)
    @test c(iter) == iter[1:5]
    c = Collect(20)
    @test c(iter) == iter
end

@testset "Accessors" begin
    p = IAEAParticle(photon,2,3,4,5,6,0,0,1,true,(),())
    q = @set p.x = 10
    @test q.x == 10
    @test p.x == 4
end

@testset "histogram" begin
    p = arbitrary(IAEAParticle{0,0})
    p = @set p.E = 1
    p = @set p.x = 1
    p = @set p.weight = 1
    p11_w1 = p
    p11_w2 = @set p11_w1.weight = 2
    p21_w1 = @set p11_w1.E = 2
    p12_w1 = @set p11_w1.x = 2
    edges = ([0,1.5,3],)
    h = @inferred histmap(p->p.E,
                                     [p11_w1, p21_w1],
                                     edges = edges)
    @test h.edges == edges
    @test h.weights == [1,1]

    h = @inferred histmap(p->p.E,
                                     [p11_w2, p11_w1, p21_w1],
                                     edges=edges)
    @test h.weights == [3,1]
    h = @inferred histmap(p->p.E,
                                     [p11_w2, p11_w1, p21_w1],
                                     edges=edges, weight_function=p->1)
    @test h.weights == [2,1]

    edges = ([0,1.5,3], [0.9,1.1,2.1, 2.2])
    h = @inferred histmap(p->p.E, p->p.x,
                                     [p11_w2, p12_w1, p12_w1, p21_w1],
                                     edges=edges)
    @test h.edges == edges
    @test h.weights == [2 2 0; 1 0 0]
end

@testset "download" begin
    key = "VarianClinaciX_6MV_05x05"
    dir = tempdir()
    path = iaea_download(key, dir=dir)
    path2 = iaea_download(key, dir=dir)
    @test path2 == path
    ps = iaea_iterator(collect, path)
    @test length(ps) == 432410
end

@testset "binning 1d" begin
    items = [4, 1.1, 1.3, -3, 3, 3]
    edges=(0:4,)
    b = binning(abs, items, edges=edges)
    @test b.edges == edges
    @test sort(items) == sort(vcat(b.content...))

    @test b.content[1] == Float64[]
    @test b.content[2] == [1.1, 1.3]
    @test b.content[3] == [-3.,3,3]
    @test b.content[4] == [4.]
    @test_throws BoundsError b.content[5]
end

@testset "binning 2d" begin
    items = [
        (x=1.5, y=1.5),
        (x=2.5, y=1.5),
        (x=3.5, y=1.5),
        (x=3.5, y=2.5),
        (x=3.5, y=2.5),
    ]
    edges=(1:4,1:3)
    b = binning(items, edges=edges) do item
        item.x, item.y
    end
    @test b.edges == edges
    @test sort(items) == sort(vcat(b.content...))
    @test b.content[1,1] == [(x = 1.5, y = 1.5)]
    @test b.content[2,1] == [(x = 2.5, y = 1.5)]
    @test b.content[3,1] == [(x = 3.5, y = 1.5)]
    @test b.content[3,2] == [(x = 3.5, y = 2.5), (x=3.5,y=2.5)]
end

end


text:
function get_corrected_measurement(obj::KalmanEstimation)
    return jcall(obj, "getCorrectedMeasurement", EstimatedMeasurement, ())
end

function get_corrected_spacecraft_states(obj::KalmanEstimation)
    return jcall(obj, "getCorrectedSpacecraftStates", Vector{SpacecraftState}, ())
end

function get_current_date(obj::KalmanEstimation)
    return jcall(obj, "getCurrentDate", AbsoluteDate, ())
end

function get_current_measurement_number(obj::KalmanEstimation)
    return jcall(obj, "getCurrentMeasurementNumber", jint, ())
end

function get_estimated_measurements_parameters(obj::KalmanEstimation)
    return jcall(obj, "getEstimatedMeasurementsParameters", ParameterDriversList, ())
end

function get_estimated_orbital_parameters(obj::KalmanEstimation)
    return jcall(obj, "getEstimatedOrbitalParameters", ParameterDriversList, ())
end

function get_estimated_propagation_parameters(obj::KalmanEstimation)
    return jcall(obj, "getEstimatedPropagationParameters", ParameterDriversList, ())
end

function get_physical_estimated_covariance_matrix(obj::KalmanEstimation)
    return jcall(obj, "getPhysicalEstimatedCovarianceMatrix", RealMatrix, ())
end

function get_physical_estimated_state(obj::KalmanEstimation)
    return jcall(obj, "getPhysicalEstimatedState", RealVector, ())
end

function get_physical_innovation_covariance_matrix(obj::KalmanEstimation)
    return jcall(obj, "getPhysicalInnovationCovarianceMatrix", RealMatrix, ())
end

function get_physical_kalman_gain(obj::KalmanEstimation)
    return jcall(obj, "getPhysicalKalmanGain", RealMatrix, ())
end

function get_physical_measurement_jacobian(obj::KalmanEstimation)
    return jcall(obj, "getPhysicalMeasurementJacobian", RealMatrix, ())
end

function get_physical_state_transition_matrix(obj::KalmanEstimation)
    return jcall(obj, "getPhysicalStateTransitionMatrix", RealMatrix, ())
end

function get_predicted_measurement(obj::KalmanEstimation)
    return jcall(obj, "getPredictedMeasurement", EstimatedMeasurement, ())
end

function get_predicted_spacecraft_states(obj::KalmanEstimation)
    return jcall(obj, "getPredictedSpacecraftStates", Vector{SpacecraftState}, ())
end


text:
using AlphaZero
using AlphaZero.Examples: games, experiments
using AlphaZero.Scripts: dummy_run, test_game

using Test

const CI = get(ENV, "CI", nothing) == "true"
const FULL = !CI

@testset "Testing Games" begin
  test_game(games["grid-world"])
  test_game(games["tictactoe"])
  test_game(games["connect-four"])
  @test true
end

@testset "Dummy Runs" begin
  dir = "sessions/test-tictactoe"
  @test dummy_run(experiments["tictactoe"], nostdout=false) == nothing
  if FULL
    @test dummy_run(experiments["connect-four"], nostdout=false) == nothing
  end
end


text:
# A Pythagorean triplet is a set of three natural numbers, a < b < c, for which,
#
# a^2 + b^2 = c^2
# For example, 3^2 + 4^2 = 9 + 16 = 25 = 5^2.
#
# There exists exactly one Pythagorean triplet for which a + b + c = 1000.
# Find the product abc.

using ProjectEulerSolutions

# Cycle through (from n down to n/2) to find sets of a and b that satisfy
# a + b + c = n.
function p009solution(n::Integer=4)::Integer

    maxval = fld(n, 2)
    for a = maxval:-1:2
        for b = maxval:-1:2
            c = sqrt(a*a + b*b)
            if c == floor(c) && a + b + c == n
                return a * b * c
            end
        end
    end

    # Should never happen
    return -1
end

# Note, there are probably faster options which use triplet generators.

p009 = Problems.Problem(p009solution)

Problems.benchmark(p009, 1000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
module SimpleMath
    using ExportPublic

    _secret_pi = 22/7
    my_pi = _secret_pi

    function add_two(a::Int, b::Int)
        _implementation_detail(a) + _implementation_detail(b)
    end

    function _implementation_detail(a::Int)
        a
    end

    function add_squared(a::Int, b::Int)
        _square_it(a) + _square_it(b)
    end

    function _square_it(a::Int)
        a ^ 2
    end

    @exportPublic()
end


text:
# Julia wrapper for header: /usr/include/scip/cons_superindicator.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrSuperindicator(scip)
    ccall((:SCIPincludeConshdlrSuperindicator, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateConsSuperindicator(scip, cons, name, binvar, slackcons, initial, separate, enforce, check, propagate, _local, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsSuperindicator, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_VAR}, Ptr{SCIP_CONS}, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, binvar, slackcons, initial, separate, enforce, check, propagate, _local, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicSuperindicator(scip, cons, name, binvar, slackcons)
    ccall((:SCIPcreateConsBasicSuperindicator, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_VAR}, Ptr{SCIP_CONS}), scip, cons, name, binvar, slackcons)
end

function SCIPgetBinaryVarSuperindicator(cons)
    ccall((:SCIPgetBinaryVarSuperindicator, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_CONS},), cons)
end

function SCIPgetSlackConsSuperindicator(cons)
    ccall((:SCIPgetSlackConsSuperindicator, libscip), Ptr{SCIP_CONS}, (Ptr{SCIP_CONS},), cons)
end

function SCIPtransformMinUC(scip, success)
    ccall((:SCIPtransformMinUC, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{UInt32}), scip, success)
end

function SCIPdialogExecChangeMinUC(scip, dialog, dialoghdlr, nextdialog)
    ccall((:SCIPdialogExecChangeMinUC, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_DIALOG}, Ptr{SCIP_DIALOGHDLR}, Ptr{Ptr{SCIP_DIALOG}}), scip, dialog, dialoghdlr, nextdialog)
end


text:
module Errno
include("errno_h.jl")
export
    E2BIG,
    EACCES,
    EADDRINUSE,
    EADDRNOTAVAIL,
    EADV,
    EAFNOSUPPORT,
    EAGAIN,
    EALREADY,
    EBADE,
    EBADF,
    EBADFD,
    EBADMSG,
    EBADR,
    EBADRQC,
    EBADSLT,
    EBFONT,
    EBUSY,
    ECANCELED,
    ECHILD,
    ECHRNG,
    ECOMM,
    ECONNABORTED,
    ECONNREFUSED,
    ECONNRESET,
    EDEADLK,
    EDESTADDRREQ,
    EDOM,
    EDOTDOT,
    EDQUOT,
    EEXIST,
    EFAULT,
    EFBIG,
    EHOSTDOWN,
    EHOSTUNREACH,
    EHWPOISON,
    EIDRM,
    EILSEQ,
    EINPROGRESS,
    EINTR,
    EINVAL,
    EIO,
    EISCONN,
    EISDIR,
    EISNAM,
    EKEYEXPIRED,
    EKEYREJECTED,
    EKEYREVOKED,
    EL2HLT,
    EL2NSYNC,
    EL3HLT,
    EL3RST,
    ELIBACC,
    ELIBBAD,
    ELIBEXEC,
    ELIBMAX,
    ELIBSCN,
    ELNRNG,
    ELOOP,
    EMEDIUMTYPE,
    EMFILE,
    EMLINK,
    EMSGSIZE,
    EMULTIHOP,
    ENAMETOOLONG,
    ENAVAIL,
    ENETDOWN,
    ENETRESET,
    ENETUNREACH,
    ENFILE,
    ENOANO,
    ENOBUFS,
    ENOCSI,
    ENODATA,
    ENODEV,
    ENOENT,
    ENOEXEC,
    ENOKEY,
    ENOLCK,
    ENOLINK,
    ENOMEDIUM,
    ENOMEM,
    ENOMSG,
    ENONET,
    ENOPKG,
    ENOPROTOOPT,
    ENOSPC,
    ENOSR,
    ENOSTR,
    ENOSYS,
    ENOTBLK,
    ENOTCONN,
    ENOTDIR,
    ENOTEMPTY,
    ENOTNAM,
    ENOTRECOVERABLE,
    ENOTSOCK,
    ENOTTY,
    ENOTUNIQ,
    ENXIO,
    EOPNOTSUPP,
    EOVERFLOW,
    EOWNERDEAD,
    EPERM,
    EPFNOSUPPORT,
    EPIPE,
    EPROTO,
    EPROTONOSUPPORT,
    EPROTOTYPE,
    ERANGE,
    EREMCHG,
    EREMOTE,
    EREMOTEIO,
    ERESTART,
    ERFKILL,
    EROFS,
    ESHUTDOWN,
    ESOCKTNOSUPPORT,
    ESPIPE,
    ESRCH,
    ESRMNT,
    ESTALE,
    ESTRPIPE,
    ETIME,
    ETIMEDOUT,
    ETOOMANYREFS,
    ETXTBSY,
    EUCLEAN,
    EUNATCH,
    EUSERS,
    EXDEV,
    EXFULL
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
import Knet.Ops21: conv, relu

using CUDA.CUDNN:
    cudnnConvolutionForward,
    cudnnConvolutionForward!,
    cudnnConvolutionBackwardFilter,
    cudnnConvolutionBackwardData,
    cudnnGetConvolutionNdForwardOutputDim,
    cudnnSetConvolutionMathType,
    cudnnSetConvolutionReorderType,
    cudnnSetConvolutionGroupCount,
    cudnnFindConvolutionForwardAlgorithmEx,
        cudnnConvolutionFwdAlgoPerf_t,
    cudnnFindConvolutionBackwardFilterAlgorithmEx,
        cudnnConvolutionBwdFilterAlgoPerf_t,
    cudnnFindConvolutionBackwardDataAlgorithmEx,
        cudnnConvolutionBwdDataAlgoPerf_t,
    cudnnConvolutionDescriptor,
        cudnnConvolutionDescriptor_t,
        cudnnCreateConvolutionDescriptor,
        cudnnSetConvolutionNdDescriptor,
        cudnnDestroyConvolutionDescriptor,
    cudnnConvolutionForwardOutput,
    cudnnTensorDescriptor,
    cudnnFilterDescriptor,
    cudnnConvolutionMode_t,
        CUDNN_CONVOLUTION,       # 0
        CUDNN_CROSS_CORRELATION, # 1
    cudnnActivationMode_t,
        CUDNN_ACTIVATION_SIGMOID,      # 0
        CUDNN_ACTIVATION_RELU,         # 1
        CUDNN_ACTIVATION_TANH,         # 2
        CUDNN_ACTIVATION_CLIPPED_RELU, # 3
        CUDNN_ACTIVATION_ELU,          # 4
        CUDNN_ACTIVATION_IDENTITY,     # 5
    cudnnNanPropagation_t,
        CUDNN_NOT_PROPAGATE_NAN, # 0
        CUDNN_PROPAGATE_NAN,     # 1
    cudnnMathType_t,
        CUDNN_DEFAULT_MATH,                    # 0
        CUDNN_TENSOR_OP_MATH,                  # 1
        CUDNN_TENSOR_OP_MATH_ALLOW_CONVERSION, # 2
        CUDNN_FMA_MATH,                        # 3
    cudnnReorderType_t,
        CUDNN_DEFAULT_REORDER, # 0
        CUDNN_NO_REORDER,      # 1
    cudnnConvolutionFwdAlgo_t,
        CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_GEMM,         # 0
        CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_PRECOMP_GEMM, # 1
        CUDNN_CONVOLUTION_FWD_ALGO_GEMM,                  # 2
        CUDNN_CONVOLUTION_FWD_ALGO_DIRECT,                # 3
        CUDNN_CONVOLUTION_FWD_ALGO_FFT,                   # 4
        CUDNN_CONVOLUTION_FWD_ALGO_FFT_TILING,            # 5
        CUDNN_CONVOLUTION_FWD_ALGO_WINOGRAD,              # 6
        CUDNN_CONVOLUTION_FWD_ALGO_WINOGRAD_NONFUSED,     # 7
        CUDNN_CONVOLUTION_FWD_ALGO_COUNT,                 # 8
    cudnnConvolutionBwdFilterAlgo_t,
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_0,                 # 0, /* non-deterministic */
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_1,                 # 1,
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_FFT,               # 2,
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_3,                 # 3, /* non-deterministic */
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_WINOGRAD,          # 4, /* not implemented */
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_WINOGRAD_NONFUSED, # 5,
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_FFT_TILING,        # 6,
        CUDNN_CONVOLUTION_BWD_FILTER_ALGO_COUNT,             # 7
    cudnnConvolutionBwdDataAlgo_t,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_0,                 # 0, /* non-deterministic */
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_1,                 # 1,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_FFT,               # 2,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_FFT_TILING,        # 3,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_WINOGRAD,          # 4,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_WINOGRAD_NONFUSED, # 5,
        CUDNN_CONVOLUTION_BWD_DATA_ALGO_COUNT,             # 6
    cudnnTensorFormat_t,
        CUDNN_TENSOR_NCHW,        # 0, /* row major (wStride = 1, hStride = w) */
        CUDNN_TENSOR_NHWC,        # 1, /* feature maps interleaved ( cStride = 1 )*/
        CUDNN_TENSOR_NCHW_VECT_C, # 2, /* each image point is vector of element of C, vector length in data type */
    cudnnDataType,
    convdims,
    math_mode,
    handle

function conv(
    w::GPUVal, x::GPUVal;
    z = nothing,
    bias = nothing,
    activation = nothing,
    normalization = nothing,

    alpha::Real = 1,
    beta::Real = 0,

    padding::Union{Integer,Vector{<:Integer},Tuple{<:Integer,Vararg{Int}}} = 0,  # >= 0
    stride::Union{Integer,Vector{<:Integer},Tuple{<:Integer,Vararg{Int}}} = 1,   # >= 1
    dilation::Union{Integer,Vector{<:Integer},Tuple{<:Integer,Vararg{Int}}} = 1, # >= 1
    groups::Integer = 1,
    crosscorrelation::Bool = true,
    channelmajor::Bool = false,

    format::cudnnTensorFormat_t = channelmajor ? CUDNN_TENSOR_NHWC : CUDNN_TENSOR_NCHW,
    mode::cudnnConvolutionMode_t = crosscorrelation ? CUDNN_CROSS_CORRELATION : CUDNN_CONVOLUTION,
    mathType::cudnnMathType_t = math_mode(),
    reorderType::cudnnReorderType_t = CUDNN_DEFAULT_REORDER,
    convDesc::cudnnConvolutionDescriptor = cudnnConvolutionDescriptor(convdims(padding,size(x),format), convdims(stride,size(x),format), convdims(dilation,size(x),format), mode, cudnnDataType(eltype(x)), mathType, reorderType, Cint(groups)),
    y = cudnnConvolutionForwardOutput(x, cudnnTensorDescriptor(x;format), cudnnFilterDescriptor(w;format), convDesc, format),

    dw = Ref{Any}(nothing),
    dx = Ref{Any}(nothing),
    dz = Ref{Any}(nothing),
    dbias = Ref{Any}(nothing),
)
    a = (activation === relu && normalization ∈ (nothing, identity) ? CUDNN_ACTIVATION_RELU : CUDNN_ACTIVATION_IDENTITY)
    r = cudnnConvolutionForward!(y, w, x, convDesc; activation=a, bias, z, alpha, beta, format, dw, dx, dz, dbias)
    if normalization ∉ (nothing, identity); r = normalization(r); end
    if a === CUDNN_ACTIVATION_IDENTITY && activation ∉ (nothing, identity); r = activation(r); end
    return r
end


text:
using Knet, Test, Base.Iterators, Printf, LinearAlgebra, Random, CUDA, IterTools
Knet.atype() = KnetArray{Float32}
import Pkg
Pkg.add("ArgParse")
using ArgParse
include("tree.jl")
include("util.jl")
using JLD2

function process_train_data(opt)
    word_manager = SymbolsManager(true)
    init_from_file(word_manager, "$(opt["data_dir"])/vocab.q.txt", opt["min_freq"], opt["max_vocab_size"])
    form_manager = SymbolsManager(true)
    init_from_file(form_manager,"$(opt["data_dir"])/vocab.f.txt", 0, opt["max_vocab_size"])
    data = []
    open("$(opt["data_dir"])/$(opt["train"]).txt", "r") do file
        for line in eachline(file)
            l_list = split(line, "\t")
            w_list = get_symbol_idx_for_list(word_manager, split(strip(l_list[1]), ' '))
            r_list = get_symbol_idx_for_list(form_manager, split(strip(l_list[2]), ' '))
            cur_tree = convert_to_tree(r_list, 1, length(r_list), form_manager)
            push!(data, (w_list, r_list, cur_tree))
        end
    end
    out_mapfile = string(opt["data_dir"]) * "/map.jld2"
    out_datafile = string(opt["data_dir"]) * "/train.jld2"
    managers = [word_manager, form_manager]
    @save out_mapfile managers
    @save out_datafile  data
end

function serialize_data(opt)
    data = []
    #managers = pkl.load( open("{}/map.jld2".format(opt["data_dir"]), "rb" ) )
    @load  (string(opt["data_dir"]) * "/map.jld2") managers
    word_manager, form_manager = managers
    #-------
    open("$(opt["data_dir"])/$(opt["test"]).txt", "r") do file
        for line in eachline(file)
            l_list = split(line, "\t")
            w_list = get_symbol_idx_for_list(word_manager, split(strip(l_list[1]), ' '))
            r_list = get_symbol_idx_for_list(form_manager, split(strip(l_list[2]), ' '))
            cur_tree = convert_to_tree(r_list, 1, length(r_list), form_manager)
            #data.append((w_list, r_list, cur_tree))
            push!(data, (w_list, r_list, cur_tree))
        end
    end
    out_datafile = string(opt["data_dir"]) * "/test.jld2"
    @save out_datafile  data
    #-------
end

s = ArgParseSettings()
@add_arg_table! s begin
    "--data_dir"
        arg_type = String
        help = "data dir"
        default = "../data/"
    "--train"
        help = "train dir"
        arg_type = String
        default = "train"
    "--test"
        help = "minimum word frequency"
        arg_type = String
        default = "test"
    "--min_freq"
        help = "minimum word frequency"
        arg_type = Int
        default = 2
    "--max_vocab_size"
        help = "max vocab size"
        arg_type = Int
        default = 15000
    "--seed"
        help = "random number generator seed"
        arg_type = Int
        default = 123
end

parsed_args = parse_args(s)
Random.seed!(parsed_args["seed"])
process_train_data(parsed_args)
serialize_data(parsed_args)


text:
@testset "Params delete!" begin
  w = rand(2,3)
  b = rand(2)
  ps = Params([w,b])
  delete!(ps, b)
  @test length(ps.order) == length(ps.params) == 1
  @test first(ps.order) == first(ps.params) == w
end

@testset "copyto!" begin
  x = [0,0]
  ps = Params([x])
  copyto!(ps, [1, 2])
  @test x == [1, 2]

  x = [0,0]
  y = [0]
  ps = Params([x, y])
  copyto!(ps, [1, 2, 3])
  @test x == [1, 2]
  @test y == [3]

  ps = Params([[1,2]])
  x = [0, 0]
  copyto!(x, ps)
  @test x == [1, 2]

  ps = Params([[1,2], [3]])
  x = [0, 0, 0]
  copyto!(x, ps)
  @test x == [1, 2, 3]
end


text:
function adjust_axis(x_lim,y_lim)

	# scaling factors
	al_x = [0.05, 0.05];  # x axis (low, high)
	al_y = [0.05, 0.05];  # y axis (low, high)

	# additional axis movement
	if x_lim[1]==0.0; a=-1; else a=0; end
	if x_lim[2]==0.0; b=1; else b=0; end
	if y_lim[1]==0.0; c=-1; else c=0; end
	if y_lim[2]==0.0; d=1; else d=0; end

	xlim = Float64[0,0]; ylim = Float64[0,0];
	xlim[1] = x_lim[1]-abs(x_lim[1]*al_x[1])+a;
	xlim[2] = x_lim[2]+abs(x_lim[2]*al_x[2])+b;
	ylim[1] = y_lim[1]-abs(y_lim[1]*al_y[1])+c;
	ylim[2] = y_lim[2]+abs(y_lim[2]*al_y[2])+d;

	xlims!((xlim[1],xlim[2]))
	ylims!((ylim[1],ylim[2]))
end

"""
allPlots(n;idx=2)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/10/2017, Last Modified: 11/10/2017 \n
--------------------------------------------------------------------------------------\n
"""
function allPlots(n;idx::Int64=1,kwargs...)
  if !isdir(n.r.resultsDir); resultsDir!(n); end
  stp = [statePlot(n,idx,st;kwargs...) for st in 1:n.ocp.state.num]
  ctp = [controlPlot(n,idx,ctr;kwargs...) for ctr in 1:n.ocp.control.num]

  if n.s.ocp.evalCostates && n.s.ocp.integrationMethod == :ps && n.s.ocp.evalConstraints
    csp = [costatePlot(n,idx,st;kwargs...) for st in 1:n.ocp.state.num]
    all = [stp;ctp;csp]
  else
    all = [stp;ctp]
  end

  h = plot(all...,size=_pretty_defaults[:size])
  if !_pretty_defaults[:simulate]; savefig(string(n.r.resultsDir,"main.",_pretty_defaults[:format])) end
  return h
end

"""
stp=statePlot(n,r.eval_num,7);
stp=statePlot(n,idx,st);
stp=statePlot(n,idx,st;(:legend=>"test1"));
stp=statePlot(n,idx,st,stp;(:append=>true));
stp=statePlot(n,idx,st,stp;(:lims=>false));
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/10/2017, Last Modified: 9/21/2017 \n
--------------------------------------------------------------------------------------\n
"""
function statePlot(n,idx::Int64,st::Int64,args...;kwargs...)
  kw = Dict(kwargs);

  # check to see if user would like to add to an existing plot
  if !haskey(kw,:append); append=false;
  else; append = get(kw,:append,0);
  end
  if !append; stp=plot(0,leg=:false); else stp=args[1]; end

  # check to see if user would like to plot limits
  if !haskey(kw,:lims); lims=true;
  else; lims=get(kw,:lims,0);
  end

  # check to see if user would like to label legend
  if !haskey(kw,:legend); legend_string = "";
  else; legend_string = get(kw,:legend,0);
  end

  # check length of n.r.ocp.dfs to see if there is data for the provided idx
  if length(n.r.ocp.dfs) >= idx
    ocpPlot = true
  else
    ocpPlot = false
  end

	if ocpPlot
  	t_vec = linspace(0.0,n.r.ocp.dfs[end][:t][end],_pretty_defaults[:L])
	else
		t_vec = linspace(0.0,n.r.ip.dfsplant[end][:t][end],_pretty_defaults[:L])
	end

  if lims
		# plot the upper limits
		if n.ocp.mXU[st]!=false
			if !isnan(n.ocp.XU[st]);plot!(n.r.ocp.tst,n.ocp.XU_var[st,:],line=_pretty_defaults[:limit_lines][2],label=string(legend_string,"max"));end
		else
			if !isnan(n.ocp.XU[st]);plot!(t_vec,n.ocp.XU[st]*ones(_pretty_defaults[:L],1),line=_pretty_defaults[:limit_lines][2],label=string(legend_string,"max"));end
		end

    # plot the lower limits
    if n.ocp.mXL[st]!=false
      if !isnan(n.ocp.XL[st]);plot!(n.r.ocp.tst,n.ocp.XL_var[st,:],line=_pretty_defaults[:limit_lines][1],label=string(legend_string,"min"));end
    else
      if !isnan(n.ocp.XL[st]);plot!(t_vec,n.ocp.XL[st]*ones(_pretty_defaults[:L],1),line=_pretty_defaults[:limit_lines][1],label=string(legend_string,"min"));end
    end
  end

  # plot the values TODO if there are no lims then you cannot really see the signal
	if !_pretty_defaults[:plantOnly] && ocpPlot
    if n.s.ocp.integrationMethod == :ps && _pretty_defaults[:polyPts] && isequal(n.r.ocp.dfsOpt[:status][idx],:Optimal) && !n.s.ocp.linearInterpolation && n.s.ocp.interpolationOn
      int_color = 1
      for int in 1:n.ocp.Ni
          if int_color > length(_pretty_defaults[:mpc_lines]) # reset colors
            int_color = 1
          end
          plot!(n.r.ocp.AlltpolyPts[idx][int],n.r.ocp.AllXpolyPts[idx][st][int],line=_pretty_defaults[:mpc_lines][int_color], label=string("poly. # ", int))
          int_color = int_color + 1
      end
      scatter!((n.r.ocp.dfs[idx][:t],n.r.ocp.dfs[idx][n.ocp.state.name[st]]),marker=_pretty_defaults[:mpc_markers],label=string(legend_string,"colloc. pts."))
    else
      plot!(n.r.ocp.dfs[idx][:t],n.r.ocp.dfs[idx][n.ocp.state.name[st]],marker=_pretty_defaults[:mpc_markers],line=_pretty_defaults[:mpc_lines][1],label=string(legend_string,"mpc"))
    end
	end
  #  MPC leads plant by an index
  if _pretty_defaults[:plant] #&& idx!=1
		# values
		temp = [n.r.ip.dfsplant[jj][n.ocp.state.name[st]] for jj in 1:idx];
	  vals=[idx for tempM in temp for idx=tempM];

		# time
		temp = [n.r.ip.dfsplant[jj][:t] for jj in 1:idx];
		time=[idx for tempM in temp for idx=tempM];

    plot!(time,vals,line=_pretty_defaults[:plant_lines][1],label=string(legend_string,"plant"));
  end
  if _pretty_defaults[:X0p]
    X0 = n.r.ip.X0p[idx][1]
    t = n.r.ip.X0p[idx][2]
    scatter!((t,X0[st]),marker=_pretty_defaults[:X0p_marker],label="X0p")
  end
  adjust_axis(xlims(),ylims());
	xlims!(t_vec[1],t_vec[end]);
  plot!(size=_pretty_defaults[:size]);
  yaxis!(n.ocp.state.description[st]); xaxis!("time (s)");
  if !_pretty_defaults[:simulate]; savefig(string(n.r.resultsDir,n.ocp.state.name[st],".",_pretty_defaults[:format])); end
  return stp
end

"""
stp=statePlot(n,idx,st1,st2);
stp=statePlot(n,idx,st1,st2;(:legend=>"test1"));
stp=statePlot(n,idx,st1,st2,stp;(:append=>true),(:lims=>false));
# to compare two different states
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/10/2017, Last Modified: 5/28/2017 \n
--------------------------------------------------------------------------------------\n
"""
function statePlot(n,idx::Int64,st1::Int64,st2::Int64,args...;kwargs...)
  kw = Dict(kwargs);

  # check to see if user would like to add to an existing plot
  if !haskey(kw,:append); kw_ = Dict(:append => false); append = get(kw_,:append,0);
  else; append = get(kw,:append,0);
  end
  if !append; stp=plot(0,leg=:false); else stp=args[1]; end

  #if idx > length(n.r.ocp.dfs)
  #  warn("Cannot plot idx = ", idx, " because length(n.r.ocp.dfs) = ", length(n.r.ocp.dfs), ". \n
  #        Skipping idx in statePlot().")
  #  return stp
  #end

  # check length of n.r.ocp.dfs to see if there is data for the provided idx
  if length(n.r.ocp.dfs) >= idx
    ocpPlot = true
  else
    ocpPlot = false
  end

  # check to see if user would like to plot limits
  if !haskey(kw,:lims); kw_ = Dict(:lims => true); lims = get(kw_,:lims,0);
  else; lims = get(kw,:lims,0);
  end

  # check to see if user would like to label legend
  if !haskey(kw,:legend); kw_ = Dict(:legend => ""); legend_string = get(kw_,:legend,0);
  else; legend_string = get(kw,:legend,0);
  end

  # plot the limits
  # TODO check if all constraints are given
	# TODO make it work for linear varying stateTol
  if lims
    if !isnan(n.ocp.XU[st1]);plot!([n.ocp.XU[st1],n.ocp.XU[st1]],[n.ocp.XL[st2],n.ocp.XU[st2]],line=_pretty_defaults[:limit_lines][2],label=string(n.ocp.state.name[st1],"_max"));end
    if !isnan(n.ocp.XL[st1]);plot!([n.ocp.XL[st1],n.ocp.XL[st1]],[n.ocp.XL[st2],n.ocp.XU[st2]],line=_pretty_defaults[:limit_lines][1],label=string(n.ocp.state.name[st1],"_min"));end
    if !isnan(n.ocp.XU[st1]);plot!([n.ocp.XL[st1],n.ocp.XU[st1]],[n.ocp.XU[st2],n.ocp.XU[st2]],line=_pretty_defaults[:limit_lines][4],label=string(n.ocp.state.name[st2],"_max"));end
    if !isnan(n.ocp.XL[st1]);plot!([n.ocp.XL[st1],n.ocp.XU[st1]],[n.ocp.XL[st2],n.ocp.XL[st2]],line=_pretty_defaults[:limit_lines][3],label=string(n.ocp.state.name[st2],"_min"));end
  end

  # plot the values
	if !_pretty_defaults[:plantOnly] && ocpPlot
		plot!(n.r.ocp.dfs[idx][n.ocp.state.name[st1]],n.r.ocp.dfs[idx][n.ocp.state.name[st2]],line=_pretty_defaults[:mpc_lines][1],label=string(legend_string,"mpc"));
	end

  if _pretty_defaults[:plant] #&& idx!=1
		# values
		temp = [n.r.ip.dfsplant[jj][n.ocp.state.name[st1]] for jj in 1:idx]
		vals1 = [idx for tempM in temp for idx=tempM]

		# values
		temp = [n.r.ip.dfsplant[jj][n.ocp.state.name[st2]] for jj in 1:idx]
		vals2 = [idx for tempM in temp for idx=tempM]

		plot!(vals1,vals2,line=_pretty_defaults[:plant_lines][1],label=string(legend_string,"plant"));
  end

  adjust_axis(xlims(),ylims())
  plot!(size=_pretty_defaults[:size])
  xaxis!(n.ocp.state.description[st1])
  yaxis!(n.ocp.state.description[st2])
  if !_pretty_defaults[:simulate] savefig(string(n.r.resultsDir,n.ocp.state.name[st1],"_vs_",n.ocp.state.name[st2],".",_pretty_defaults[:format])); end
  return stp
end

"""
ctrp=controlPlot(n,idx,ctr);
ctrp=controlPlot(n,idx,ctr,ctrp;(:append=>true));
# to plot control signals
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/10/2017, Last Modified: 9/19/2017 \n
--------------------------------------------------------------------------------------\n
"""
function controlPlot(n,idx::Int64,ctr::Int64,args...;kwargs...)
  kw = Dict(kwargs);

  # check to see if user would like to add to an existing plot
  if !haskey(kw,:append); kw_ = Dict(:append => false); append = get(kw_,:append,0);
  else; append = get(kw,:append,0);
  end
  if !append; ctrp=plot(0,leg=:false); else ctrp=args[1]; end

  #if idx > length(n.r.ocp.dfs)
  #  warn("Cannot plot idx = ", idx, " because length(n.r.ocp.dfs) = ", length(n.r.ocp.dfs), ". \n
  #        Skipping idx in controlPlot().")
  #  return ctrp
  #end
  # check length of n.r.ocp.dfs to see if there is data for the provided idx
  if length(n.r.ocp.dfs) >= idx
    ocpPlot = true
  else
    ocpPlot = false
  end

  # check to see if user would like to plot limits
  if !haskey(kw,:lims); kw_ = Dict(:lims => true); lims = get(kw_,:lims,0);
  else; lims = get(kw,:lims,0);
  end

  # check to see if user would like to label legend
  if !haskey(kw,:legend); kw_ = Dict(:legend => ""); legend_string = get(kw_,:legend,0);
  else; legend_string = get(kw,:legend,0)
  end

  if ocpPlot
  	t_vec = linspace(0.0,n.r.ocp.dfs[end][:t][end],_pretty_defaults[:L])
	else
		t_vec = linspace(0.0,n.r.ip.dfsplant[end][:t][end],_pretty_defaults[:L])
	end

  # plot the limits
  if lims
    if !isnan(n.ocp.CU[ctr]);plot!(t_vec,n.ocp.CU[ctr]*ones(_pretty_defaults[:L],1),line=_pretty_defaults[:limit_lines][2],label="max");end
    if !isnan(n.ocp.CL[ctr]); plot!(t_vec,n.ocp.CL[ctr]*ones(_pretty_defaults[:L],1),line=_pretty_defaults[:limit_lines][1],label="min"); end
  end

  if !_pretty_defaults[:plantOnly] && ocpPlot
    if n.s.ocp.integrationMethod == :ps && _pretty_defaults[:polyPts] && isequal(n.r.ocp.dfsOpt[:status][idx],:Optimal) && !n.s.ocp.linearInterpolation && n.s.ocp.interpolationOn
      int_color = 1
      for int in 1:n.ocp.Ni
        if int_color > length(_pretty_defaults[:mpc_lines]) # reset colors
          int_color = 1  # reset colors
        end
        plot!(n.r.ocp.AlltpolyPts[idx][int],n.r.ocp.AllUpolyPts[idx][ctr][int],line=_pretty_defaults[:mpc_lines][int_color], label=string("poly. # ", int))
        int_color = int_color + 1
      end
      scatter!(n.r.ocp.dfs[idx][:t],n.r.ocp.dfs[idx][n.ocp.control.name[ctr]],marker=_pretty_defaults[:mpc_markers],label=string(legend_string,"colloc. pts."))
    else
      plot!(n.r.ocp.dfs[idx][:t],n.r.ocp.dfs[idx][n.ocp.control.name[ctr]],marker=_pretty_defaults[:mpc_markers],line=_pretty_defaults[:mpc_lines][1],label=string(legend_string,"mpc"))
    end
  end

  if _pretty_defaults[:plant] #&& idx!=1
		# values
		temp = [n.r.ip.dfsplant[jj][n.ocp.control.name[ctr]] for jj in 1:idx]
	  vals = [idx for tempM in temp for idx=tempM];

		# time
		temp = [n.r.ip.dfsplant[jj][:t] for jj in 1:idx]
		time = [idx for tempM in temp for idx=tempM]

		plot!(time,vals,line=_pretty_defaults[:plant_lines][1],label=string(legend_string,"plant"));
  end
  adjust_axis(xlims(),ylims());
	xlims!(t_vec[1],t_vec[end]);
  plot!(size=_pretty_defaults[:size]);
  yaxis!(n.ocp.control.description[ctr]);	xaxis!("time (s)");
	if !_pretty_defaults[:simulate] savefig(string(n.r.resultsDir,n.ocp.control.name[ctr],".",_pretty_defaults[:format])) end
  return ctrp
end

"""
costatesPlot(n)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 11/10/2017, Last Modified: 11/10/2017 \n
--------------------------------------------------------------------------------------\n
"""
function costatesPlot(n;idx::Int64=1,kwargs...)
  if !(n.s.ocp.evalCostates && n.s.ocp.integrationMethod == :ps && n.s.ocp.evalConstraints)
    error("   if n.s.ocp.evalCostates && n.s.ocp.integrationMethod == :ps && n.s.ocp.evalConstraints  \n
              must all be true.")
  end

  if !isdir(n.r.resultsDir); resultsDir!(n); end

  csp = [costatePlot(n,idx,st;kwargs...) for st in 1:n.ocp.state.num];
  h = plot(csp...,size=_pretty_defaults[:size])
  if !_pretty_defaults[:simulate]; savefig(string(n.r.resultsDir,"cs_main.",_pretty_defaults[:format])) end
  return h
end

"""
costatePlot(n,idx,st)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 11/10/2017, Last Modified: 11/10/2017 \n
--------------------------------------------------------------------------------------\n
"""
function costatePlot(n,idx::Int64,st::Int64;kwargs...)
  if !(n.s.ocp.evalCostates && n.s.ocp.integrationMethod == :ps && n.s.ocp.evalConstraints)
    error("   if n.s.ocp.evalCostates && n.s.ocp.integrationMethod == :ps && n.s.ocp.evalConstraints  \n
              must all be true.")
  end
  kw = Dict(kwargs);

  # check to see if user would like to add to an existing plot
  if !haskey(kw,:append); append=false;
  else; append = get(kw,:append,0);
  end
  if !append; csp=plot(0,leg=:false); else csp=args[1]; end

  if length(n.r.ocp.dfs) >= idx
    ocpPlot = true
  else
    ocpPlot = false
  end

  # check to see if user would like to label legend
  if !haskey(kw,:legend); legend_string = "";
  else; legend_string = get(kw,:legend,0);
  end

	if ocpPlot
  	t_vec=linspace(0.0,n.r.ocp.dfs[end][:t][end],_pretty_defaults[:L]);
	else
		t_vec=linspace(0.0,n.r.ip.dfsplant[end][:t][end],_pretty_defaults[:L]);
	end

  if n.s.ocp.integrationMethod == :ps
    t_st_int = [n.r.ocp.tst[n.ocp.Nck_cum[int]+1:n.ocp.Nck_cum[int+1]+1] for int in 1:n.ocp.Ni] # this is redundant in interpolateLagrange!()

    int_color = 1
    for int in 1:n.ocp.Ni
        if int_color > length(_pretty_defaults[:mpc_lines]) # reset colors
          int_color = 1
        end
        plot!(n.r.ocp.AlltpolyPts[idx][int],n.r.ocp.AllCSpolyPts[idx][st][int],line=_pretty_defaults[:mpc_lines][int_color], label=string("poly. # ", int))
        int_color = int_color + 1
        scatter!(t_st_int[int][1:end-1],n.r.ocp.CS[st][int],marker=_pretty_defaults[:mpc_markers],label=string(legend_string,"costate pts."))
    end
  #else # NOTE for now :tm methods do not have a costate option
  #  plot!(n.r.ocp.dfs[idx][:t],n.r.ocp.CS[st],marker=_pretty_defaults[:mpc_markers],line=_pretty_defaults[:mpc_lines][1],label=string(legend_string,"mpc"))
  end

  adjust_axis(xlims(),ylims());
  xlims!(t_vec[1],t_vec[end]);
  plot!(size=_pretty_defaults[:size]);
  yaxis!(string(n.ocp.state.description[st]," costate")); xaxis!("time (s)");

  if !_pretty_defaults[:simulate]; savefig(string(n.r.resultsDir,n.ocp.state.name[st]," costate.",_pretty_defaults[:format])); end
  return csp
end


text:
using LinearAlgebra
using FFTW
using FastGaussQuadrature

abstract type AbstractMaxwellSolver end

export Maxwell1DFEM

"""
    maxwell_solver = MaxwellFEM1D( mesh, degree )

1D Maxwell spline finite element solver on a periodic grid

- `Lx`                   : length of Periodic domain
- `delta_x`              : cell size
- `n_dofs`               : number of cells (and grid points)
- `s_deg_0`              : spline degree 0-forms
- `s_deg_1`              : spline degree 1-forms
- `mass_0`               : coefficients of 0-form mass matrix
- `mass_1`               : coefficients of 1-form mass matrix
- `eig_mass0`            : eigenvalues of circulant 0-form mass matrix
- `eig_mass1`            : eigenvalues of circulant 1-form mass matrix
- `eig_weak_ampere`      : eigenvalues of circulant update matrix for Ampere
- `eig_weak_poisson`     : eigenvalues of circulant update matrix for Poisson
- `plan_fw`              : fft plan (forward)
- `plan_bw`              : fft plan (backward)

"""
struct Maxwell1DFEM <: AbstractMaxwellSolver
    xmin::Float64
    Lx::Float64
    delta_x::Float64
    n_dofs::Int32
    s_deg_0::Int32
    s_deg_1::Int32
    mass_0::Vector{Float64}
    mass_1::Vector{Float64}
    eig_mass0::Vector{Float64}
    eig_mass1::Vector{Float64}
    eig_weak_ampere::Vector{Float64}
    eig_weak_poisson::Vector{Float64}

    plan_fw::FFTW.FFTWPlan
    plan_bw::FFTW.FFTWPlan
    work::Vector{Float64}
    wsave::Vector{Float64}
    eigvals::Vector{Float64}

    function Maxwell1DFEM(mesh::OneDGrid, degree::Int)
        xmin = mesh.xmin[1]
        n_dofs = mesh.nx[1]
        Lx = mesh.xmax[1] - mesh.xmin[1]
        delta_x = Lx / n_dofs
        s_deg_0 = degree
        s_deg_1 = degree - 1

        mass_0 = zeros(Float64, s_deg_0 + 1)
        mass_1 = zeros(Float64, s_deg_0)

        if s_deg_0 == 1 # linear and constant splines

            # Upper diagonal coeficients  of linear spline mass matrix (from Eulerian numbers)
            mass_0[1] = 4.0 / 6.0
            mass_0[2] = 1.0 / 6.0
            # Upper diagonal coeficients  of constant spline mass matrix
            mass_1[1] = 1.0

        elseif s_deg_0 == 2 # quadratic and linear splines

            # Upper diagonal coeficients  of quadratic spline mass matrix (from Eulerian numbers)
            mass_0[1] = 66.0 / 120.0
            mass_0[2] = 26.0 / 120.0
            mass_0[3] = 1.0 / 120.0
            # Upper diagonal coeficients  of linear spline mass matrix (from Eulerian numbers)
            mass_1[1] = 4.0 / 6.0
            mass_1[2] = 1.0 / 6.0

        elseif s_deg_0 == 3

            # Upper diagonal coeficients  of cubic spline mass matrix (from Eulerian numbers)
            mass_0[1] = 2416.0 / 5040.0
            mass_0[2] = 1191.0 / 5040.0
            mass_0[3] = 120.0 / 5040.0
            mass_0[4] = 1.0 / 5040.0
            # Upper diagonal coeficients  of quadratic spline mass matrix (from Eulerian numbers)
            mass_1[1] = 66.0 / 120.0
            mass_1[2] = 26.0 / 120.0
            mass_1[3] = 1.0 / 120.0

        else
            throw(ArgumentError("Wrong value of degree = $degree  (1,2 or 3)"))
        end

        eig_mass0 = zeros(Float64, n_dofs)
        eig_mass1 = zeros(Float64, n_dofs)
        eig_weak_ampere = zeros(Float64, n_dofs)
        eig_weak_poisson = zeros(Float64, n_dofs)

        work = zeros(Float64, n_dofs)
        wsave = zeros(Float64, n_dofs)
        plan_fw = FFTW.plan_r2r(work, FFTW.R2HC)
        plan_bw = FFTW.plan_r2r(wsave, FFTW.HC2R)

        # Compute eigenvalues of circulant Ampere update matrix M_0^{-1} D^T M_1
        # and circulant Poisson Matrix (D^T M_1 D)^{-1}
        # zero mode vanishes due to derivative matrix D^T

        eig_weak_ampere[1] = 0.0
        eig_weak_poisson[1] = 0.0  # Matrix is not invertible: 0-mode is set to 0
        eig_mass0[1] = 1.0  # sum of coefficents is one
        eig_mass1[1] = 1.0  # sum of coefficents is one

        for k in 1:(n_dofs ÷ 2 - 1)
            coef0 = mass_0[1]
            coef1 = mass_1[1]
            for j in 1:(s_deg_0 - 1)
                cos_mode = cos(2 * pi * j * k / n_dofs)
                coef0 = coef0 + 2 * mass_0[j + 1] * cos_mode
                coef1 = coef1 + 2 * mass_1[j + 1] * cos_mode
            end
            # add last term for larger matrix
            j = s_deg_0
            coef0 = coef0 + 2 * mass_0[j + 1] * cos(2 * pi * j * k / n_dofs)
            # compute eigenvalues
            eig_mass0[k + 1] = coef0 # real part
            eig_mass0[n_dofs - k + 1] = 0.0 # imaginary part
            eig_mass1[k + 1] = coef1 # real part
            eig_mass1[n_dofs - k + 1] = 0.0 # imaginary part
            cos_mode = cos(2 * pi * k / n_dofs)
            sin_mode = sin(2 * pi * k / n_dofs)
            eig_weak_ampere[k + 1] = (coef1 / coef0) * (1 - cos_mode) # real part
            eig_weak_ampere[n_dofs - k + 1] = -(coef1 / coef0) * sin_mode   # imaginary part
            eig_weak_poisson[k + 1] = 1.0 / (coef1 * ((1 - cos_mode)^2 + sin_mode^2))  # real part
            eig_weak_poisson[n_dofs - k + 1] = 0.0
        end

        # N/2 mode
        coef0 = mass_0[1]
        coef1 = mass_1[1]
        for j in 1:(s_deg_0 - 1)
            coef0 = coef0 + 2 * mass_0[j + 1] * cos(pi * j)
            coef1 = coef1 + 2 * mass_1[j + 1] * cos(pi * j)
        end

        # add last term for larger matrix
        j = s_deg_0
        coef0 = coef0 + 2 * mass_0[j + 1] * cos(pi * j)

        # compute eigenvalues
        eig_mass0[n_dofs ÷ 2 + 1] = coef0
        eig_mass1[n_dofs ÷ 2 + 1] = coef1
        eig_weak_ampere[n_dofs ÷ 2 + 1] = 2.0 * (coef1 / coef0)
        eig_weak_poisson[n_dofs ÷ 2 + 1] = 1.0 / (coef1 * 4.0)

        eigvals = zeros(Float64, n_dofs)

        return new(
            xmin,
            Lx,
            delta_x,
            n_dofs,
            s_deg_0,
            s_deg_1,
            mass_0,
            mass_1,
            eig_mass0,
            eig_mass1,
            eig_weak_ampere,
            eig_weak_poisson,
            plan_fw,
            plan_bw,
            work,
            wsave,
            eigvals,
        )
    end
end

export compute_rhs_from_function!

"""
   compute_rhs_from_function(self, func, degree, coefs_dofs)

Compute the FEM right-hand-side for a given function f and periodic splines of given degree.

Its components are ``\\int f N_i dx`` where ``N_i`` is the B-spline starting at ``x_i``.
"""
function compute_rhs_from_function!(
    coefs_dofs::Vector{Float64}, self::Maxwell1DFEM, func::Function, degree::Int64
)
    bspl = zeros(Float64, (degree + 1, degree + 1))

    # take enough Gauss points so that projection is exact for splines of degree deg
    # rescale on [0,1] for compatibility with B-splines
    x, w = gausslegendre(degree + 1)

    x .= 0.5 .* (x .+ 1.0)
    w .= 0.5 .* w

    # Compute bsplines at gauss_points
    for k in 1:(degree + 1)
        bspl[k, :] .= uniform_bsplines_eval_basis(degree, x[k])
    end

    # Compute coefs_dofs = int f(x)N_i(x)
    for i in 1:(self.n_dofs)
        coef = 0.0
        # loop over support of B spline
        for j in 1:(degree + 1)
            # loop over Gauss points
            for k in 1:(degree + 1)
                coef =
                    coef +
                    w[k] * func(self.delta_x * (x[k] + i + j - 2)) * bspl[k, degree + 2 - j]
            end
        end
        # rescale by cell size
        coefs_dofs[i] = coef * self.delta_x
    end
end

function solve_circulant!(self, eigvals, rhs)
    n = self.n_dofs
    # Compute res from rhs, using eigenvalue of circulant  matrix
    self.work .= rhs
    # Forward FFT
    mul!(self.wsave, self.plan_fw, self.work)
    self.wsave[1] = self.wsave[1] * eigvals[1]
    for k in 2:(n ÷ 2)
        re_p = self.wsave[k] * eigvals[k] - self.wsave[n - k + 2] * eigvals[n - k + 2]
        im_p = self.wsave[k] * eigvals[n - k + 2] + self.wsave[n - k + 2] * eigvals[k]
        self.wsave[k] = re_p
        self.wsave[n - k + 2] = im_p
    end
    self.wsave[n ÷ 2 + 1] = self.wsave[n ÷ 2 + 1] * eigvals[n ÷ 2 + 1]
    # Backward FFT
    mul!(self.work, self.plan_bw, self.wsave)

    return self.work ./= n
end

export compute_e_from_rho!

function compute_e_from_rho!(e::Vector{Float64}, self::Maxwell1DFEM, rho::Vector{Float64})

    # Compute potential phi from rho, using eigenvalue of circulant
    # inverse matrix
    solve_circulant!(self, self.eig_weak_poisson, rho)
    # Compute spline coefficients of Ex from those of phi
    for i in 2:(self.n_dofs)
        e[i] = self.work[i - 1] - self.work[i]
    end
    # treat Periodic point
    return e[1] = self.work[self.n_dofs] - self.work[1]
end

export compute_e_from_j!
"""
    compute_e_from_j!(e, maxwell_solver, current, component)

Compute ``E_i`` from ``j_i`` integrated over the time interval using weak Ampere formulation
"""
function compute_e_from_j!(
    e::Vector{Float64}, self::Maxwell1DFEM, current::Vector{Float64}, component::Int64
)
    n = self.n_dofs
    fill!(self.eigvals, 0.0)

    # Multiply by inverse mass matrix  using the eigenvalues of the circulant
    # inverse matrix

    if (component == 1)
        for i in 1:(n ÷ 2 + 1)
            self.eigvals[i] = 1.0 / self.eig_mass1[i]
        end
        solve_circulant!(self, self.eigvals, current)
    elseif (component == 2)
        for i in 1:(n ÷ 2 + 1)
            self.eigvals[i] = 1.0 / self.eig_mass0[i]
        end
        solve_circulant!(self, self.eigvals, current)
    else
        throw(ArgumentError("Component $component not implemented "))
    end

    # Update the electric field and scale
    self.work ./= self.delta_x
    return e .-= self.work
end

export l2norm_squared

"""
    l2norm_squared(maxwell_solver, coefs_dofs, degree)

Compute square of the L2norm

"""
function l2norm_squared(self, coefs_dofs, degree)

    # Multiply coefficients by mass matrix (use diagonalization FFT and mass matrix eigenvalues)
    if (degree == self.s_deg_0)
        solve_circulant!(self, self.eig_mass0, coefs_dofs)

    elseif (degree == self.s_deg_1)
        solve_circulant!(self, self.eig_mass1, coefs_dofs)
    end

    # Multiply by the coefficients from the left (inner product)
    r = sum(coefs_dofs .* self.work)
    # Scale by delt_x
    return r .* self.delta_x
end

export l2norm_squared2

"""
    l2norm_squared(maxwell_solver, coefs_dofs, degree)

Compute square of the L2norm

"""
function l2norm_squared2(self, coefs_dofs, degree)

    # Multiply coefficients by mass matrix (use diagonalization FFT and mass matrix eigenvalues)
    if (degree == self.s_deg_0)
        solve_circulant!(self, self.eig_mass0, coefs_dofs)

    elseif (degree == self.s_deg_1)
        solve_circulant!(self, self.eig_mass1, coefs_dofs)
    end

    # Multiply by the coefficients from the left (inner product)
    r = sum(self.work .* self.work)
    # Scale by delt_x
    return r .* self.delta_x
end

export l2projection!

"""
    l2projection!(coefs_dofs, maxwell, func, degree)

Compute the L2 projection of a given function f on periodic splines
of given degree
"""
function l2projection!(
    coefs_dofs::Vector{Float64}, self::Maxwell1DFEM, func::Function, degree::Int64
)
    n = self.n_dofs
    eigvals = zeros(Float64, n)

    # Compute right-hand-side
    compute_rhs_from_function!(coefs_dofs, self, func, degree)

    # Multiply by inverse mass matrix
    if (degree == self.s_deg_0)
        for i in 1:(n ÷ 2 + 1)
            eigvals[i] = 1.0 / self.eig_mass0[i]
        end
    elseif (degree == self.s_deg_0 - 1)
        for i in 1:(n ÷ 2 + 1)
            eigvals[i] = 1.0 / self.eig_mass1[i]
        end
    else
        throw(ArgumentError("degree $degree not available"))
    end

    solve_circulant!(self, eigvals, coefs_dofs)

    # Account for scaling in the mass matrix by dx

    return coefs_dofs .= coefs_dofs ./ self.delta_x
end

export compute_e_from_b!

"""
    compute_e_from_b!(field_out, maxwell_solver, delta_t, field_in)

compute Ey from Bz using weak Ampere formulation

"""
function compute_e_from_b!(
    field_out::Vector{Float64},
    self::Maxwell1DFEM,
    delta_t::Float64,
    field_in::Vector{Float64},
)
    coef = delta_t / self.delta_x

    # Compute potential weak curl of Bz using eigenvalue of circulant inverse matrix
    solve_circulant!(self, self.eig_weak_ampere, field_in)
    # Update bz from self value
    return field_out .+= coef .* self.work
end

export compute_b_from_e!
"""
    compute_b_from_e!( field_out, maxwell_solver, delta_t, field_in)

Compute Bz from Ey using strong 1D Faraday equation for spline coefficients
```math
B_z^{new}(x_j) = B_z^{old}(x_j) - \\frac{\\Delta t}{\\Delta x} (E_y(x_j) - E_y(x_{j-1})
```
"""
function compute_b_from_e!(
    field_out::Vector{Float64},
    self::Maxwell1DFEM,
    delta_t::Float64,
    field_in::Vector{Float64},
)
    coef = delta_t / self.delta_x
    # relation betwen spline coefficients for strong Ampere
    for i in 2:(self.n_dofs)
        field_out[i] = field_out[i] + coef * (field_in[i - 1] - field_in[i])
    end
    # treat Periodic point
    return field_out[1] = field_out[1] + coef * (field_in[end] - field_in[1])
end

export compute_rderivatives_from_basis!

function compute_rderivatives_from_basis!(
    field_out::Vector{Float64}, self::Maxwell1DFEM, field_in::Vector{Float64}
)
    coef = 1 / self.delta_x
    # relation betwen spline coefficients for strong Ampere
    for i in 1:(self.n_dofs - 1)
        @inbounds field_out[i] = coef * (field_in[i] - field_in[i + 1])
    end
    # treat Periodic point
    return field_out[end] = coef * (field_in[end] - field_in[1])
end

export compute_lderivatives_from_basis!

function compute_lderivatives_from_basis!(
    field_out::Vector{Float64}, self::Maxwell1DFEM, field_in::Vector{Float64}
)
    coef = 1 / self.delta_x
    # relation betwen spline coefficients for strong Ampere
    for i in 2:(self.n_dofs)
        @inbounds field_out[i] = coef * (field_in[i] - field_in[i - 1])
    end
    # treat Periodic point
    return field_out[1] = coef * (field_in[1] - field_in[end])
end

"""
    inner_product( maxwell_solver, coefs1_dofs, coefs2_dofs, degree )

-  `maxwell_solver` : Maxwell solver object
-  `coefs1_dofs` : Coefficient for each DoF
-  `coefs2_dofs` : Coefficient for each DoF
-  `degree : Specify the degree of the basis functions

return squared L2 norm

"""
function inner_product(self, coefs1_dofs, coefs2_dofs, degree)

    # Multiply coefficients by mass matrix (use diagonalization FFT and mass matrix eigenvalues)
    if degree == self.s_deg_0
        solve_circulant!(self, self.eig_mass0, coefs2_dofs)

    elseif degree == self.s_deg_1
        solve_circulant!(self, self.eig_mass1, coefs2_dofs)
    end

    # Multiply by the coefficients from the left (inner product)
    # Scaled by delt_x

    return sum(coefs1_dofs .* self.work) * self.delta_x
end


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

plot(x=[0.236033, 0.346517, 0.312707, 0.00790928, 0.488613, 0.210968, 0.951916, 0.999905, 0.251662, 0.986666, 0.555751, 0.437108],
     y=[0.424718, 0.773223, 0.28119, 0.209472, 0.251379, 0.0203749, 0.287702, 0.859512, 0.0769509, 0.640396, 0.873544, 0.278582],
     color=repeat(["a","b","c"], outer=[4]),
     Scale.color_discrete_manual("red","purple","green"))


text:
module BinaryHeaps

export
    MutableBinaryHeapWithReduction,
    update!

using TypeParams
using UnPack
using Base.Order

include("walk.jl")
include("binary_heaps.jl")
include("permutations.jl")
include("reduction.jl")

include("MutableBinaryHeapsWithReduction.jl")

end # module


text:
module VoxelPlots

include("MagicaVoxel.jl")
include("Goxel.jl")

end # module


text:
module PkgEval

import Pkg.TOML
using Pkg
using Base: UUID
using Dates
using ProgressMeter
using DataFrames
using Random
using Libdl
using Sockets
using JSON
using Printf

# immutable: in package directory
versions_file() = joinpath(dirname(@__DIR__), "deps", "Versions.toml")
releases_file() = joinpath(dirname(@__DIR__), "deps", "Releases.toml")
registries_file() = joinpath(dirname(@__DIR__), "deps", "Registries.toml")

# mutable: in .cache directory
cache_dir() = joinpath(get(ENV, "XDG_CACHE_HOME", joinpath(homedir(), ".cache")), "PkgEval")
download_dir(name) = joinpath(cache_dir(), "downloads", name)
storage_dir() = joinpath(cache_dir(), "storage")
extra_versions_file() = joinpath(cache_dir(), "Versions.toml")

# fixed locations
registry_dir() = joinpath(first(DEPOT_PATH), "registries")
registry_dir(name) = joinpath(registry_dir(), name)

# utils
isdebug(group) = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, group, PkgEval) !== nothing

include("registry.jl")
include("julia.jl")
include("run.jl")
include("report.jl")

function __init__()
    mkpath(cache_dir())
    touch(extra_versions_file())
end

end # module