text:
#Methods to transform a PEP to a companion linearized form and solve the corresponding linearized pencil
using SparseArrays

export companion
export polyeig #Wrapper around the solver for a linearized PEP pencil

"""
    E,A = companion(nep::PEP);

Linearizes a  polynomial eigenvalue problem (PEP) a to the companion form, as in the paper by Mehrmann and Voss.
More precisely, for a k-th degree PEP with n-by-n coefficient matrices,
this returns matrices E and A, both kn-by-kn, corresponding to the linearized problem
```math
Ax = λEx
```

# Example
```julia-repl
julia> pep = nep_gallery("pep0");
julia> E,A = companion(pep);
julia> λ, V = eigen(A,E);
julia> minimum(svd(compute_Mder(pep,λ[1])).S)
2.4845217786736996e-12
```

# References
* V. Mehrmann and H. Voss, Non-linear eigenvalue problems, a challenge for modern eigenvalue methods, GAMM‐Mitteilungen (2004)
"""
    function companion(pep::PEP)

        n = size(pep,1) #Size of monomial coefficient matrices
        d = size(pep.A,1)-1 #Degree of pep
        T = eltype(pep.A[1]) #Deduce the type of elements in the PEP matrices

        ##### Check sparsity of the problem and allocate memory accordingly #####
        if issparse(pep)
            E = spzeros(T, d*n, d*n)
            A = spzeros(T, n*d, n*d)

            #(d-1)n-by-(d-1)n matrix (Used to construct both E and A)
            Iblock = sparse(kron(Matrix{T}(I, d-1, d-1), Matrix{T}(I, n, n)))
        else
            E = zeros(T, d*n, d*n)
            A = zeros(T, n*d, n*d)

            Iblock = kron(Matrix{T}(I, d-1, d-1), Matrix{T}(I, n, n))
        end

        E[1:n,1:n] = pep.A[d+1];#Fill block (1,1)
        E[n+1:d*n,n+1:d*n] = Iblock;#Fill all blocks on the diagonal with eye(n)

        #####Construct A #####

        #First row block of A
        for i=1:d
           A[1:n,(i-1)*n+1:i*n] = pep.A[d-i+1];
        end
        #Lower part of A
        A[n+1:d*n,1:(d-1)*n] = T(-1.0)*Iblock

        return E,-A

    end

"""
    λ,v = polyeig([eltype],nep::PEP,[eigsolvertype,])

Linearizes a  polynomial eigenvalue problem (PEP) a to the companion form
and solves the corresponding linear eigenvalue problem; see [`companion`](@ref).
The `eigsolvertype` is optinal can be used to specify how the linear problem
is solved; see [`eig_solve`](@ref), and [`EigSolver`](@ref).

# Example
```julia-repl
julia> pep = nep_gallery("pep0");
julia> λ,V = polyeig(pep);
julia> minimum(svd(compute_Mder(pep,λ[1])).S)
1.2050763381899922e-14
julia> norm(compute_Mlincomb(pep,λ[2],vec(V[:,2])))
1.0245470569036458e-12
```
"""
    polyeig(pep::PEP,vargs...)=polyeig(ComplexF64,pep,vargs...)
    function polyeig(::Type{T},pep::PEP,eigsolvertype::Type=DefaultEigSolver) where T

        #Linearize to Ax = λEx
        E,A = companion(pep);

        solver::EigSolver = eigsolvertype(A,E);
        D,V = eig_solve(solver,target=one(T),nev=size(A,1));

        D = Vector{T}(D)
        V = Matrix{T}(V)

        return D,V[1:size(pep,1),:]
    end

"""
    λ,v = polyeig([eltype],nep::ChebPEP)

Computes a companion linearization for the NEP represented
in a Chebyshev basis, and returns eigenpairs.

# Example
```julia
julia> using LinearAlgebra
julia> nep=nep_gallery("dep0");
julia> chebpep=ChebPEP(nep,9,-3,1,cosine_formula_cutoff=5);
julia> (λv,V)=polyeig(chebpep);
julia> ii=argmin(abs.(λv));
julia> λ=λv[ii];
julia> v=V[:,ii];
julia> norm(compute_Mlincomb(chebpep,λ,v))
1.3543968603949142e-14
julia> # Actually, it's not a bad approx to the original NEP either
julia> norm(compute_Mlincomb(nep,λ,v))
4.326355966047557e-6
```

See also [`ChebPEP`](@ref).

# References:

* Amiraslani, A., Corless, R. M. & Lancaster, P. "Linearization of matrix polynomials expressed in poly-nomial bases" IMA J. Numer. Anal.,29 (2009): 141–157.
* Effenberger and Kressner. "Chebyshev interpolation for nonlinear eigenvalue problems." BIT Numerical Mathematics 52.4 (2012): 933-951.

"""
polyeig(pep::ChebPEP,vargs...)=polyeig(ComplexF64,pep,vargs...)
function polyeig(::Type{T}, pep::ChebPEP{TT,Ftype}) where {T,Ftype,TT}
    k=pep.k
    Fk=get_Av(pep);

    n=size(pep,1);
    L0=zeros(T,n*(k-1),n*(k-1));
    L1=zeros(T,n*(k-1),n*(k-1));
    II=Matrix{Ftype}(I,n,n);

    for j=1:(k-2)
        L0[((j-1)*n) .+ (1:n), j*n .+ (1:n)]=II
        L0[j*n .+ (1:n), ((j-1)*n) .+ (1:n)]=II
    end

    for j=1:k-1
        L0[((k-2)*n) .+ (1:n), (n*(j-1)).+ (1:n)]=-Fk[j];
    end

    L0[((k-2)*n) .+ (1:n), (n*(k-3)) .+ (1:n)] += Fk[k];

    for j=1:k-2
        factor=2;
        if (j==1)
            factor=1
        end
        L1[((j-1)*n) .+ (1:n), ((j-1)*n) .+ (1:n)]=factor*II
    end
    L1[((k-2)*n) .+ (1:n),
       ((k-2)*n) .+ (1:n)]=2*Fk[k]

    LL=eigen(L0,L1);

    a=pep.a; b=pep.b;
    return ((b-a)*(LL.values .+ 1)/2 .+ a,
            hcat(map(i -> normalize(LL.vectors[1:n,i]), 1:(n*(k-1)))...))

end


text:
module Bonobo

using DataStructures
using NamedTupleTools

"""
    AbstractNode

The abstract type for a tree node. Your own type for `Node` given to [`initialize`](@ref) needs to subtype it.
The default if you don't provide your own is [`DefaultNode`](@ref).
"""
abstract type AbstractNode end

"""
    AbstractSolution{Node<:AbstractNode, Value}

The abstract type for a `Solution` object. The default is [`DefaultSolution`](@ref).
It is parameterized by `Node` and `Value` where `Value` is the value which describes the full solution i.e the value for every variable.
"""
abstract type AbstractSolution{Node<:AbstractNode, Value} end

"""
    BnBNode

Holds the necessary information of every node.
This needs to be added by every `AbstractNode` as `std::BnBNode`

```julia
id :: Int
lb :: Float64
ub :: Float64
```
"""
mutable struct BnBNode
    id :: Int
    lb :: Float64
    ub :: Float64
end

"""
    DefaultNode <: AbstractNode

The default structure for saving node information.
Currently this includes only the necessary `std::BnBNode` which needs to be part of every [`AbstractNode`](@ref).
"""
mutable struct DefaultNode <: AbstractNode
    std :: BnBNode
end

"""
    DefaultSolution{Node<:AbstractNode,Value} <: AbstractSolution{Node, Value}

The default struct to save a solution of the branch and bound run.
It holds
```julia
objective :: Float64
solution  :: Value
node      :: Node
```
Both the `Value` and the `Node` type are determined by the [`initialize`](@ref) method.

`solution` holds the information to obtain the solution for example the values of all variables.
"""
mutable struct DefaultSolution{Node<:AbstractNode,Value} <: AbstractSolution{Node, Value}
    objective :: Float64
    solution  :: Value
    node      :: Node
end

"""
    AbstractTraverseStrategy

The abstract type for a traverse strategy.
If you implement a new traverse strategy this must be the supertype.

If you want to implement your own strategy the [`get_next_node`](@ref) function needs a new method
which dispatches on the `traverse_strategy` argument.
"""
abstract type AbstractTraverseStrategy end

"""
    AbstractBranchStrategy

The abstract type for a branching strategy.
If you implement a new branching strategy, this must be the supertype.

If you want to implement your own strategy, you must implement a new method for [`get_branching_variable`](@ref)
which dispatches on the `branch_strategy` argument.
"""
abstract type AbstractBranchStrategy end

"""
    BFS <: AbstractTraverseStrategy

The BFS traverse strategy always picks the node with the lowest bound first.
If there is a tie then the smallest node id is used as a tie breaker.
"""
struct BFS <: AbstractTraverseStrategy end

"""
    FIRST <: AbstractBranchStrategy

The `FIRST` strategy always picks the first variable which isn't fixed yet and can be branched on.
"""
struct FIRST <: AbstractBranchStrategy end

"""
    MOST_INFEASIBLE <: AbstractBranchStrategy

The `MOST_INFEASIBLE` strategy always picks the variable which is furthest away from being "fixed" and can be branched on.
"""
struct MOST_INFEASIBLE <: AbstractBranchStrategy end

mutable struct Options
    traverse_strategy   :: AbstractTraverseStrategy
    branch_strategy     :: AbstractBranchStrategy
    atol                :: Float64
    rtol                :: Float64
end

"""
    BnBTree{Node<:AbstractNode,Root,Value,Solution<:AbstractSolution{Node,Value}}

Holds all the information of the branch and bound tree.
"""
mutable struct BnBTree{Node<:AbstractNode,Root,Value,Solution<:AbstractSolution{Node,Value}}
    incumbent::Float64
    lb::Float64
    solutions::Vector{Solution}
    nodes::PriorityQueue{Int,Node}
    root::Root
    branching_indices::Vector{Int}
    num_nodes::Int
    sense::Symbol
    options::Options
end
Base.broadcastable(x::BnBTree) = Ref(x)

include("util.jl")
include("node.jl")
include("branching.jl")

"""
    initialize(; kwargs...)

Initialize the branch and bound framework with the the following arguments.
Later it can be dispatched on `BnBTree{Node, Root, Solution}` for various methods.

# Keyword arguments
- `traverse_strategy` [`BFS`] currently the only supported traverse strategy is [`BFS`](@ref). Should be an [`AbstractTraverseStrategy`](@ref)
- `branch_strategy` [`FIRST`] currently the only supported branching strategies are [`FIRST`](@ref) and [`MOST_INFEASIBLE`](@ref). Should be an [`AbstractBranchStrategy`](@ref)
- `atol` [1e-6] the absolute tolerance to check whether a value is discrete
- `rtol` [1e-6] the relative tolerance to check whether a value is discrete
- `Node` [`DefaultNode`](@ref) can be special structure which is used to store all information about a node.
    - needs to have `AbstractNode` as the super type
    - needs to have `std :: BnBNode` as a field (see [`BnBNode`](@ref))
- `Solution` [`DefaultSolution`](@ref) stores the node and several other information about a solution
- `root` [`nothing`] the information about the root problem. The type can be used for dispatching on types
- `sense` [`:Min`] can be `:Min` or `:Max` depending on the objective sense
- `Value` [`Vector{Float64}`] the type of a solution

Return a [`BnBTree`](@ref) object which is the input for [`optimize!`](@ref).
"""
function initialize(;
    traverse_strategy = BFS,
    branch_strategy = FIRST,
    atol = 1e-6,
    rtol = 1e-6,
    Node = DefaultNode,
    Value = Vector{Float64},
    Solution = DefaultSolution{Node,Value},
    root = nothing,
    sense = :Min,
)
    return BnBTree{Node,typeof(root),Value,Solution}(
        NaN,
        NaN,
        Vector{Solution}(),
        PriorityQueue{Int,Node}(),
        root,
        get_branching_indices(root),
        0,
        sense,
        Options(traverse_strategy(), branch_strategy(), atol, rtol)
    )
end

function get_branching_indices end

"""
    optimize!(tree::BnBTree)

Optimize the problem using a branch and bound approach.

The steps are the following:
```julia
while !terminated(tree) # as long as there are open nodes
    # get the next open node depending on the traverse strategy
    node = get_next_node(tree, tree.options.traverse_strategy)
    # needs to be implemented by you
    # Should evaluate the current node and return the lower and upper bound
    # if the problem is infeasible both values should be set to NaN
    lb, ub = evaluate_node!(tree, node)
    # updates the upper and lower bound of the node struct
    set_node_bound!(tree.sense, node, lb, ub)

    # update the best solution
    updated = update_best_solution!(tree, node)
    updated && bound!(tree, node.id)

    # remove the current node
    close_node!(tree, node)
    # needs to be implemented by you
    # calls get_branching_variable and branch_on_variable! the latter must be implemented by you
    branch!(tree, node)
end
```

every function of the above can be overriden by your own method.
"""
function optimize!(tree::BnBTree)
    while !terminated(tree)
        node = get_next_node(tree, tree.options.traverse_strategy)
        lb, ub = evaluate_node!(tree, node)
        # if the problem was infeasible we simply close the node and continue
        if isnan(lb) && isnan(ub)
            close_node!(tree, node)
            continue
        end

        set_node_bound!(tree.sense, node, lb, ub)

        updated = update_best_solution!(tree, node)
        updated && bound!(tree, node.id)

        close_node!(tree, node)
        branch!(tree, node)
    end
end

"""
    terminated(tree::BnBTree)

Return true when the branch and bound loop in [`optimize!`](@ref) should be terminated.
Default behavior is to terminate the loop only when no nodes exist in the priority queue.
"""
terminated(tree::BnBTree) = isempty(tree.nodes)

"""
    set_node_bound!(objective_sense::Symbol, node::AbstractNode, lb, ub)

Set the bounds of the `node` object to the lower and upper bound given.
Internally everything is stored as a minimization problem. Therefore the objective_sense `:Min`/`:Max` is needed.
"""
function set_node_bound!(objective_sense::Symbol, node::AbstractNode, lb, ub)
    if isnan(ub)
        ub = Inf
    end
    if objective_sense == :Min
        node.lb = lb
        node.ub = ub
    else
        node.lb = -lb
        node.ub = -ub
    end
end

"""
    bound!(tree::BnBTree, current_node_id)

Close all nodes which have a lower bound higher or equal to the incumbent
"""
function bound!(tree::BnBTree, current_node_id)
    for (_,node) in tree.nodes
        if node.id != current_node_id && node.lb >= tree.incumbent
            close_node!(tree, node)
        end
    end
end

"""
    close_node!(tree::BnBTree, node::AbstractNode)

Delete the node from the priority queue from the nodes
"""
close_node!(tree::BnBTree, node::AbstractNode) = delete!(tree.nodes, node.id)

"""
    update_best_solution!(tree::BnBTree, node::AbstractNode)

Update the best solution when we found a better incumbent.
Calls [`add_new_solution!`] if this is the case
"""
function update_best_solution!(tree::BnBTree, node::AbstractNode)
    isinf(node.ub) && return false
    node.ub >= tree.incumbent && return false

    tree.incumbent = node.ub

    add_new_solution!(tree, node)
    return true
end

"""
    add_new_solution!(tree::BnBTree{N,R,V,S}, node::AbstractNode) where {N,R,V,S<:DefaultSolution{N,V}}

Currently it changes the general solution itself by calling [`get_relaxed_values`](@ref) which needs to be implemented by you.
# Todo: Add a possibility to store several solutions based on some options.

This function needs to be implemented by you if you have a different type of Solution object than [`DefaultSolution`](@ref).
"""
function add_new_solution!(tree::BnBTree{N,R,V,S}, node::AbstractNode) where {N,R,V,S<:DefaultSolution{N,V}}
    sol = DefaultSolution(node.ub, get_relaxed_values(tree, node), node)
    if isempty(tree.solutions)
        push!(tree.solutions, sol)
    else
        tree.solutions[1] = sol
    end
end

"""
    get_relaxed_values(tree::BnBTree, node::AbstractNode)

Get the values of the current node. This is always called only after [`evaluate_node!`](@ref) is called.
It is used to store a `Solution` object.
Return the type of `Value` given to the [`initialize`](@ref) method.
"""
function get_relaxed_values end

"""
    get_solution(tree::BnBTree; result=1)

Return the solution values of the problem.
See [`get_objective_value`](@ref) to obtain the objective value instead.
"""
function get_solution(tree::BnBTree{N,R,V,S}; result=1) where {N,R,V,S<:DefaultSolution{N,V}}
    return tree.solutions[result].solution
end

"""
    get_objective_value(tree::BnBTree; result=1)

Return the objective value
"""
function get_objective_value(tree::BnBTree{N,R,V,S}; result=1) where {N,R,V,S<:DefaultSolution{N,V}}
    if tree.sense == :Max
        return -tree.solutions[result].objective
    else
        return tree.solutions[result].objective
    end
end

export BnBTree, BnBNode, AbstractNode, AbstractSolution

export AbstractTraverseStrategy

end

ext: jl


text:
function fetch_aliases(xroot)
    alias_nodes = findall("//@alias", xroot)
    OrderedDict(
        alias.parentelement["name"] => alias.parentelement["alias"] for alias ∈ alias_nodes
    )
end

function build_alias_graph(alias_verts, aliases_dict)
    aliases_g = SimpleDiGraph(length(alias_verts))

    for (j, (src, dst)) ∈ enumerate(aliases_dict)
        i = findfirst(dst .== alias_verts)
        add_edge!(aliases_g, i, j)
    end
    aliases_g
end

const aliases_dict = fetch_aliases(xroot)
const aliases_dict_rev = OrderedDict(v => k for (k, v) ∈ aliases_dict)

"""
Whether this type is an alias for another name.
"""
isalias(name) = name ∈ keys(aliases_dict)

"""
Whether an alias was built from this name.
"""
hasalias(name) = name ∈ values(aliases_dict)

alias_verts = unique(vcat(keys(aliases_dict)..., values(aliases_dict)...))

aliases_g = build_alias_graph(alias_verts, aliases_dict)

aliases(aliases_g::SimpleDiGraph, index) = getindex.(Ref(alias_verts), outneighbors(aliases_g, index))
aliases(name) = (index = findfirst(alias_verts .== name); isnothing(index) ? String[] : aliases(aliases_g, index))

follow_alias(aliases_g::SimpleDiGraph, index) = (indices = inneighbors(aliases_g, index); isempty(indices) ? alias_verts[index] : length(indices) > 1 ? error("More than one indices returned for $(alias_verts[index]) when following alias $(getindex.(Ref(alias_verts), indices))") : alias_verts[first(indices)])
follow_alias(name) = (index = findfirst(alias_verts .== name); isnothing(index) ? name : follow_alias(aliases_g, index))

# for el ∈ topological_sort_by_dfs(aliases_g)
#     elname = alias_verts[el]
#     println("$(lpad("$elname ($el)", 120)) => $(isalias(elname) ? "ALIAS" : "NOALIAS")")
# end

@assert follow_alias("VkPhysicalDeviceMemoryProperties2KHR") == "VkPhysicalDeviceMemoryProperties2"
@assert follow_alias("VkPhysicalDeviceMemoryProperties2") == "VkPhysicalDeviceMemoryProperties2"


text:
module RenderManager

export RenderObject

using CoreExtended

using ..GraphicsManager
using ..StorageManager
using ..ModelManager
using ..CameraManager
using ..EntityManager
using ..ShaderManager

const Entity=EntityManager.Typ
const Model=ModelManager.Typ
const Camera=CameraManager.Typ

render = ()->nothing # render whole szene
OnPreRender = ()->nothing
OnPostRender = ()->nothing

function setOnRender(f1::Function, f2::Function)
	global OnPreRender = f1
	global OnPostRender = f2
end

type RenderObject
	id::Symbol
	order::Symbol
	enabled::Bool
	entity::Entity
	model::Model
	camera::Camera
	shaderProgram::AbstractGraphicsShaderProgram
	modes::Dict{Symbol,Any}
	setModes::Function
	RenderObject(id=:_) = new(id,:_,true,nothing,nothing,nothing,nothing,Dict(),()->nothing)
end

function link(this::RenderObject)
	ModelManager.linkTo(this.model)
	ShaderManager.linkTo(this.shaderProgram)
	this.camera.entity=this.entity
	CameraManager.update(this.camera)
end

function register(this::RenderObject, createEntity=true, createModel=false, createShaderProgram=false, createCamera=false)
	this.entity = createEntity ? EntityManager.create(this.id) : EntityManager.getCurrent()
	this.model = createModel ? ModelManager.create(this.id) : ModelManager.getCurrent()
	this.camera = createCamera ? CameraManager.create(this.id) : CameraManager.getCurrent()
	this.shaderProgram = ShaderManager.getCurrent()
	#RenderManager.createShaderProgram(:FG, "shaders/default.glsl")
end

function update()
	ShaderManager.watchShaderReload(RenderManager.reload)
end

function init(this::RenderObject)
	setMode(this, :POLYGONMODE,[:FRONT_AND_BACK,:FILL])
end

function reload()
	#all=Expr(:block,:())

	sortedList=createSortedDict(RenderObject)

	order=0
	for (k,this) in list
		order+=1
		ex=:()
		for (k,v) in this.modes
			ex=:($(ex.args...); $(GRAPHICSDRIVER().getMode(k,v).args...);)
		end

		this.setModes=()->eval(:(function(); $(ex.args...); end;))

		if this.order == :_ this.order=Symbol(order) end
		if !haskey(sortedList,this.order) sortedList[this.order]=this; end
	end

	global render = stabilize(eval(:(function();
		for (_,this) in $sortedList;
			if !this.enabled; continue; end;
				setCurrent(this);
				link(); # shaderprogram, model, texture, ...
				getShaderProperties(this); 	# set propberties
				this.setModes();
				$OnPreRender();
				renderModel(this);
				$OnPostRender();
			end;
	end;)))
end

function start()
	reset()
end

function resets()
	ModelManager.reset()
	ShaderManager.reset()
end

function unlinks()
	ModelManager.unlink()
	ShaderManager.unlink()
end

# preset for manager: includes code here
presetManager(RenderManager,RenderObject)

function copy(l::RenderObject,r::RenderObject)
	l.model = r.model
	l.camera = r.camera
	l.entity = r.entity
	l.shaderProgram = r.shaderProgram
end

function unlink(typ::DataType)
	#obj = typ == Camera ? EMPTY_CAMERA : (typ == Model ? EMPTY_MODEL : (typ == Model ? nothing : nothing))
	if obj != nothing link(obj) end
end

function link(ent::Entity)
	this = getCurrent()
	this.entity = ent
end

function link(camera::Camera)
	this = getCurrent()
	this.camera = camera
end

function link(model::Model)
	this = getCurrent()
	this.model = model
end

function link(program::AbstractGraphicsShaderProgram)
	this = getCurrent()
	this.shaderProgram = program
end

is(id::Symbol) = RenderManager.getCurrent().id == id
isEntity(id::Symbol) = RenderManager.getCurrent().entity.id == id
isModel(id::Symbol) = RenderManager.getCurrent().model.id == id
isCamera(id::Symbol) = RenderManager.getCurrent().camera.id == id

getMVP() = RenderManager.getCurrent().camera.mvpMat

getData(k::Symbol) = ModelManager.getData(getCurrent().model,k)

function init()
	ShaderManager.setListenerOnShaderPropertyUpdate(:OnShaderPropertyUpdateAttributes, (
		"iVertex" => () ->  getData(:VERTEX)
		,"iUV" => () -> (d=getData(:UV); d != nothing ? d : getData(:VERTEX))
		,"iNormal" => () -> (d=getData(:NORMAL); d != nothing ? d : getData(:VERTEX))
		)
	)
end

function getShaderProperties(this::RenderObject)
	ShaderManager.getShaderProperties(this.shaderProgram, ModelManager.getMainBuffer(this.model).data)
end

setRenderMode(this::RenderObject, mode::Symbol) = GRAPHICSDRIVER().update(ModelManager.getGroup(this.model,:VERTEX).data, mode)
setMode(this::RenderObject, mode::Symbol, value::Any) = this.modes[mode]=value

setRenderMode(mode::Symbol) = setRenderMode(getCurrent(), mode)
setMode(mode::Symbol, value::Any) = GRAPHICSDRIVER().setMode(mode,value)

function renderModel(this::RenderObject)
	model=this.model
	buffer=ModelManager.getMainBuffer(model)
	if model.enabled && ModelManager.isLinked(model) && buffer != nothing && buffer.data != nothing && (data=StorageManager.getData(buffer, 1)) != nothing
		ShaderManager.render(buffer.data, data)
	end
end

end #RenderManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
abstract type VariantIterator end

function offset_first_variant(x::Pgen)::UInt64
    x.header.variant_block_offsets[1]
end

@inline function Base.eltype(vi::VariantIterator)
    Variant
end

"""
    VariantIteratorBase(p::Pgen)
Variant iterator that iterates from the beginning of the Pgen file
"""
struct VariantIteratorBase <: VariantIterator
    p::Pgen
    v::Variant
end

function Base.iterate(vi::VariantIteratorBase,
    state=(1, offset_first_variant(vi.p)))
    if state[1] > vi.p.header.n_variants
        return nothing
    else
        idx = state[1]
        vi.v.index = state[1]
        vi.v.offset = state[2]
        vi.v.record_type = vi.p.header.variant_types[idx]
        vi.v.length = vi.p.header.variant_lengths[idx]
        nextstate = (idx + 1, state[2] + vi.v.length)
        return (vi.v, nextstate)
    end
end

@inline function Base.length(vi::VariantIteratorBase)
    vi.p.header.n_variants
end

@inline function Base.size(vi::VariantIteratorBase)
    (vi.p.header.n_variants, )
end

"""
    iterator(p::Pgen; startidx=1)

Retrieve a variant iterator for `p`.
"""
function iterator(p::Pgen; startidx=1)
    if startidx == 1
        v = Variant(0, 0, 0, 0)
        VariantIteratorBase(p, v)
    else
        @assert false "Not implemented."
    end
end

@inline function set_first_variant!(v::Variant, p::Pgen)
    v.index = 1
    v.offset = offset_first_variant(p)
    v.record_type = p.header.variant_types[1]
    v.length = p.header.variant_lengths[1]
    v
end


text:
# Code to compute the actual sensitivity indices

using OnlineStats

function computeSensSi(S::SobolEvals)

    sensSi = zeros(S.P, size(S.fA,2))

    vy = var(S.fA,1)

    for i = 1:S.P
        vari = sum(S.fA .* S.fC[i],1)./(S.N-1.0)
        mmi = mean(S.fA,1).*mean(S.fC[i],1)
        sensSi[i,:] = (vari.-mmi)./vy
    end

    sensSi
end

function computeSensSiT(S::SobolEvals)

    sensSiT = zeros(S.P, size(S.fA,2))

    vy = var(S.fA,1)

    for i = 1:S.P
        vari = (1/(2*S.N)).*sum((S.fB .- S.fC[i]).^2,1)
        sensSiT[i,:] = vari./vy
    end

    sensSiT
end

function computeSensSiboot!(sensSi::Array{Float64,2},S::SobolEvals,inds)

    vy = var(S.fA[inds,:],1)

    for i = 1:S.P
        vari = sum(S.fA[inds,:] .* S.fC[i][inds,:],1)
        mmi = mean(S.fA[inds,:],1).* mean(S.fC[i][inds,:],1)
        sensSi[i,:] = (vari/(S.N-1.0)-mmi)./vy
    end

    sensSi
end

function computeSensSiTboot!(sensSiT::Array{Float64,2},S::SobolEvals,inds)

    vy = var(S.fA[inds,:],1)

    for i = 1:S.P
        vari = (1/(2*S.N)).*sum((S.fB[inds,:] .- S.fC[i][inds,:]).^2,1)
        sensSiT[i,:] = vari./vy
    end

    sensSiT
end

function computeSensSij(S::SobolEvalsSingleAndJoint)

    sensSi = zeros(S.P, size(S.fA,2))
    sensSij = zeros(length(S.Cpairs), size(S.fA,2))

    vy = var(S.fA,1)

    for i = 1:S.P
        vari = sum(S.fA .* S.fC[i],1)
        mmi = mean(S.fA,1).* mean(S.fC[i],1)
        sensSi[i,:] = (vari/(S.N-1.0)-mmi)./vy
    end

    for i = 1:length(S.Cpairs)
        vari = sum(S.fA .* S.fCij[i],1)
        mmi = mean(S.fA,1).* mean(S.fCij[i],1)
        sensSij[i,:] = (vari/(S.N-1.0)-mmi)./vy
    end
    for i = 1:length(S.Cpairs)
        sensSij[i,:] -= (sensSi[S.Cpairs[i][1],:] .+ sensSi[S.Cpairs[i][2],:])
    end

    sensSi, sensSij
end

function computeSensSijT(S::SobolEvalsSingleAndJoint)

    sensSiT = zeros(S.P, size(S.fA,2))

    vy = var(S.fA,1)

    for i = 1:S.P
        vari = (1/(2*S.N)).*sum((S.fB .- S.fC[i]).^2,1)
        sensSiT[i,:] = vari./vy
    end

    sensSiT
end

function computeSensSijboot!(sensSi::Array{Float64,2},sensSij::Array{Float64,2},S::SobolEvalsSingleAndJoint,inds)

    vy = var(S.fA[inds,:],1)

    for i = 1:S.P
        vari = sum(view(S.fA,inds,:) .* view(S.fC[i],inds,:),1)
        mmi = mean(view(S.fA,inds,:),1).* mean(view(S.fC[i],inds,:),1)
        sensSi[i,:] = (vari/(S.N-1.0)-mmi)./vy
    end

    for i = 1:length(S.Cpairs)
        vari = sum(view(S.fA,inds,:) .* view(S.fCij[i],inds,:),1)
        mmi = mean(view(S.fA,inds,:),1).* mean(view(S.fCij[i],inds,:),1)
        sensSij[i,:] = (vari/(S.N-1.0)-mmi)./vy
    end
    for i = 1:length(S.Cpairs)
        sensSij[i,:] -= (sensSi[S.Cpairs[i][1],:] .+ sensSi[S.Cpairs[i][2],:])
    end

    sensSi, sensSij
end

function computeSensSijTboot!(sensSiT::Array{Float64,2},S::SobolEvalsSingleAndJoint,inds)

    vy = var(S.fA[inds,:],1)

    for i = 1:S.P
        vari = (1/(2*S.N)).*sum((S.fB[inds,:] .- S.fC[i][inds,:]).^2,1)
        sensSiT[i,:] = vari./vy
    end

    sensSiT
end

function sensSijSTD(S::SobolEvalsSingleAndJoint, N::Int)

    sensSij = computeSensSij(S)
    sensSTD = Array{OnlineStats.Series,2}(size(sensSij[1]))
    apxhist = Array{OnlineStats.Hist,2}(size(sensSij[1]))
    sensSijSTD = Array{OnlineStats.Series,2}(size(sensSij[2]))
    apxhistSij = Array{OnlineStats.Hist,2}(size(sensSij[2]))
    for i = 1:size(sensSij[1],1), j = 1:size(sensSij[1],2)
        apxhist[i,j] = Hist(-0.1:0.01:1.1)
        sensSTD[i,j] = Series(sensSij[1][i,j], apxhist[i,j], Mean(), Variance(), Extrema())
    end
    for i = 1:size(sensSij[2],1), j = 1:size(sensSij[2],2)
        apxhistSij[i,j] = Hist(-0.1:0.01:1.1)
        sensSijSTD[i,j] = Series(sensSij[2][i,j], apxhistSij[i,j], Mean(), Variance(), Extrema())
    end

    for i = 2:N
        samp = rand(1:S.N,S.N)
        computeSensSijboot!(sensSij[1],sensSij[2],S,samp)
        for i = 1:size(sensSij[1],1), j = 1:size(sensSij[1],2)
            fit!(sensSTD[i,j],sensSij[1][i,j])
        end
	for i = 1:size(sensSij[2],1), j = 1:size(sensSij[2],2)
            fit!(sensSijSTD[i,j],sensSij[2][i,j])
        end
    end
    sensSTD,sensSijSTD
end

function sensSiSTD(S::SobolEvals, N::Int)

    sensSi = computeSensSi(S)
    sensSTD = Array{OnlineStats.Series,2}(size(sensSi))
    apxhist = Array{OnlineStats.Hist,2}(size(sensSi))
    for i = 1:size(sensSi,1), j = 1:size(sensSi,2)
        apxhist[i,j] = Hist(-0.1:0.001:1.1)
        sensSTD[i,j] = Series(sensSi[i,j], apxhist[i,j], Mean(), Variance(), Extrema())
    end

    for i = 2:N
        samp = rand(1:S.N,S.N)
        computeSensSiboot!(sensSi,S,samp)
        for i = 1:size(sensSi,1), j = 1:size(sensSi,2)
            fit!(sensSTD[i,j],sensSi[i,j])
        end
    end
    sensSTD
end

function sensSiTSTD(S::SobolEvals, N::Int)

    sensSiT = computeSensSiT(S)
    sensSTD = Array{OnlineStats.Series,2}(size(sensSiT))
    apxhist = Array{OnlineStats.Hist,2}(size(sensSiT))
    for i = 1:size(sensSiT,1), j = 1:size(sensSiT,2)
        apxhist[i,j] = Hist(-0.1:0.001:1.1)
        sensSTD[i,j] = Series(sensSiT[i,j], apxhist[i,j], Mean(), Variance(), Extrema())
    end

    for i = 2:N
        samp = rand(1:S.N,S.N)
        computeSensSiTboot!(sensSiT,S,samp)
        for i = 1:size(sensSiT,1), j = 1:size(sensSiT,2)
            fit!(sensSTD[i,j],sensSiT[i,j])
        end
    end
    sensSTD
end

max_stars_repo_name: mstobb/VarianceBasedSensitivity


text:
using SSE
using Test

@testset "binning" begin
    sb = DynamicBin(10, Float64, Int32)
    @test length(sb) == 0
    for i = 1:10
        push!(sb, (2.3, 3))
    end
    @test length(sb) == 10
    @test sb.data[1][4] == 2.3
    @test sb.data[2][4] == 3
    @test fitted(sb)
    push!(sb, (2.3, 3))
    @test !fitted(sb)
end


text:
using Revise
using Geodesy
using MinimalRides
using Dates
using MinimalRides: make_subset, DistanceSmoothing, results, mileage, smooth, load_raw, Ride
using DataFrames
using Underscores
using BenchmarkTools

x = LatLon(-27.468937, 153.023628)
ECEF(LLA(x), wgs84)[1]
ECEF(x, wgs84)[1]

origin = [1, 153.023628, -27.468937] # City Hall, Brisbane, Australia
point = [2, 153.025900, -27.465933]  # Central Station, Brisbane, Australia

v = [1234, 153.023628, -27.468937]
MinimalRides.Location{ECEF}(MinimalRides.Location(v))

p1 = MinimalRides.Location(origin)
p1a = MinimalRides.Location{ECEF}(origin)
p2 = MinimalRides.Location(point)
p2a = MinimalRides.Location{ECEF}(point)

MinimalRides.distance(p1, p2) * 1000  # 402.27127280646783
sqrt(MinimalRides.distance(p1a, p2a)) # 401.5431022017651

using BenchmarkTools
@benchmark MinimalRides.distance($p1, $p2)
# BenchmarkTools.Trial:
#   memory estimate:  0 bytes
#   allocs estimate:  0
#   --------------
#   minimum time:     44.489 ns (0.00% GC)
#   median time:      45.148 ns (0.00% GC)
#   mean time:        48.273 ns (0.00% GC)
#   maximum time:     629.762 ns (0.00% GC)
#   --------------
#   samples:          10000
#   evals/sample:     990

@benchmark MinimalRides.distance(Ref($p1a)[], Ref($p2a)[])
# BenchmarkTools.Trial:
#   memory estimate:  0 bytes
#   allocs estimate:  0
#   --------------
#   minimum time:     1.599 ns (0.00% GC)
#   median time:      1.615 ns (0.00% GC)
#   mean time:        1.680 ns (0.00% GC)
#   maximum time:     23.345 ns (0.00% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1000

########################################
# JSON load test
########################################
url = "https://github.com/epogrebnyak/rides-minimal/blob/master/sample_jsons/sample_jsons.zip?raw=true"
cache = "/tmp/sample_json.zip"
load_raw(url, cache);

########################################
# LatLon to ECEF conversion
########################################
RAW_DATA_URL_UPDATED = ("https://dl.dropboxusercontent.com" *
    "/sh/hibzl6fkzukltk9/AABSMicBJlwMnlmA3ljt1uY5a" *
    "/data_samples-json2.zip")
CACHE_DATA = "/media/win/Data/data_samples-json2.zip"

rides = load_raw(RAW_DATA_URL_UPDATED, CACHE_DATA) .|> Ride;
rides_ecef = Ride{ECEF}.(rides);

@benchmark MinimalRides.Ride{ECEF}.($rides)

DAYS = Date.(["2020-05-21"])
TYPES = [:freight]

subset = make_subset(rides, DAYS, TYPES)

subset_ecef = make_subset(rides_ecef, DAYS, TYPES)
map(x -> x.info.id, subset) == map(x -> x.info.id, subset_ecef) # Sanity check

mileage(rides[10])
mileage(rides_ecef[10])

@benchmark mileage($rides[10])
# BenchmarkTools.Trial:
#   memory estimate:  0 bytes
#   allocs estimate:  0
#   --------------
#   minimum time:     275.165 μs (0.00% GC)
#   median time:      277.340 μs (0.00% GC)
#   mean time:        295.902 μs (0.00% GC)
#   maximum time:     2.431 ms (0.00% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1
@benchmark mileage($rides_ecef[10])
# BenchmarkTools.Trial:
#   memory estimate:  0 bytes
#   allocs estimate:  0
#   --------------
#   minimum time:     6.500 μs (0.00% GC)
#   median time:      6.909 μs (0.00% GC)
#   mean time:        8.071 μs (0.00% GC)
#   maximum time:     22.354 μs (0.00% GC)
#   --------------
#   samples:          10000
#   evals/sample:     5

smooth(DistanceSmoothing(0.5), rides[10]) |> length
smooth(DistanceSmoothing(0.5), rides_ecef[10]) |> length

@benchmark smooth(DistanceSmoothing(0.5), $rides[10])
# BenchmarkTools.Trial:
#   memory estimate:  24.53 KiB
#   allocs estimate:  10
#   --------------
#   minimum time:     290.215 μs (0.00% GC)
#   median time:      295.307 μs (0.00% GC)
#   mean time:        313.093 μs (0.00% GC)
#   maximum time:     2.755 ms (0.00% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1
@benchmark smooth(DistanceSmoothing(0.5), $rides_ecef[10])
# BenchmarkTools.Trial:
#   memory estimate:  32.50 KiB
#   allocs estimate:  10
#   --------------
#   minimum time:     14.488 μs (0.00% GC)
#   median time:      15.557 μs (0.00% GC)
#   mean time:        17.449 μs (3.25% GC)
#   maximum time:     5.704 ms (99.40% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1

df = results(subset; smoothalg = DistanceSmoothing(0.5), search_radius = 2.5 * 1.2) |> DataFrame
pairs_df = @_ sort(df, :op, rev = true) |> filter((_.op > 0.75) & (_.len > 400), __)
@benchmark(results($subset))
# BenchmarkTools.Trial:
#   memory estimate:  3.45 MiB
#   allocs estimate:  1721
#   --------------
#   minimum time:     1.008 s (0.00% GC)
#   median time:      1.036 s (0.00% GC)
#   mean time:        1.036 s (0.00% GC)
#   maximum time:     1.061 s (0.00% GC)
#   --------------
#   samples:          5
#   evals/sample:     1

df_ecef = results(subset_ecef; smoothalg = DistanceSmoothing(0.5), search_radius = 2.5 * 1.2) |> DataFrame
pairs_df_ecef = @_ sort(df_ecef, :op, rev = true) |> filter((_.op > 0.75) & (_.len > 400), __)
@benchmark(results($subset_ecef))
# BenchmarkTools.Trial:
#   memory estimate:  3.65 MiB
#   allocs estimate:  1721
#   --------------
#   minimum time:     41.444 ms (0.00% GC)
#   median time:      43.135 ms (0.00% GC)
#   mean time:        43.243 ms (0.00% GC)
#   maximum time:     50.773 ms (0.00% GC)
#   --------------
#   samples:          116
#   evals/sample:     1


text:
using OrdinaryDiffEq
using DataDrivenDiffEq
using ModelingToolkit
using Distributions
using Random
using Plots

Random.seed!(123);

function sir_ode(u,p,t)
    (s,i,r) = u
    (β,γ) = p
    ds = -β*s*i
    di = β*s*i - γ*i
    dr = γ*i
    [ds,di,dr]
end;

p = [0.5,0.25]
u0 = [0.99, 0.01, 0.0]
tspan = (0.0, 40.0)
δt = 1;

solver = ExplicitRK();

sir_prob = ODEProblem(sir_ode, u0, tspan, p)
sir_sol = solve(sir_prob, solver, saveat = δt);

dd_prob = ContinuousDataDrivenProblem(sir_sol);

@parameters t
@variables u[1:3](t)
Ψ = Basis([u; u[1]*u[2]], u, independent_variable = t)

res_koopman = solve(dd_prob, Ψ, DMDPINV())
sys_koopman = result(res_koopman);

equations(sys_koopman)

parameter_map(res_koopman)

res_sindy = solve(dd_prob, Ψ, STLSQ(),digits=1)
sys_sindy = result(res_sindy);

equations(sys_sindy)

parameter_map(res_sindy)

sir_data = Array(sir_sol);

A = 500.0 # Smaller values of A = noisier data
noisy_data = sir_data
# Note that we can't draw from u0 as R(0)=0
for i in 2:size(sir_data)[2]
    noisy_data[1:3,i] = rand(Dirichlet(A*sir_data[1:3,i]))
end;

scatter(sir_sol.t,noisy_data',title="Noisy data",xlabel="Time",ylabel="Proportion",labels=["S+noise" "I+noise" "R+noise"])
plot!(sir_sol,labels=["S" "I" "R"],legend=:left)

noisy_dd_prob = ContinuousDataDrivenProblem(noisy_data,sir_sol.t,GaussianKernel());

noisy_res_koopman = solve(noisy_dd_prob, Ψ, DMDPINV())
noisy_sys_koopman = result(noisy_res_koopman)
equations(noisy_sys_koopman), parameter_map(noisy_res_koopman)

noisy_res_sindy = solve(noisy_dd_prob, Ψ, STLSQ())
noisy_sys_sindy = result(noisy_res_sindy)
equations(noisy_sys_sindy), parameter_map(noisy_res_sindy)

A = 50.0 # Smaller values of A = noisier data
v_noisy_data = sir_data
# Note that we can't draw from u0 as R(0)=0
for i in 2:size(sir_data)[2]
    v_noisy_data[1:3,i] = rand(Dirichlet(A*sir_data[1:3,i]))
end;

scatter(sir_sol.t,v_noisy_data',title="Very noisy data",xlabel="Time",ylabel="Proportion",labels=["S+noise" "I+noise" "R+noise"])
plot!(sir_sol,labels=["S" "I" "R"],legend=:left)

v_noisy_dd_prob = ContinuousDataDrivenProblem(v_noisy_data,sir_sol.t,GaussianKernel())
v_noisy_res_sindy = solve(v_noisy_dd_prob, Ψ, STLSQ())
v_noisy_sys_sindy = result(v_noisy_res_sindy)
equations(v_noisy_sys_sindy), parameter_map(v_noisy_res_sindy)


text:
@system Planting(Stage) begin
    planting_date ~ hold

    planted(planting_date, t=calendar.time) => (t >= planting_date) ~ flag
end


text:
@testset "Mesh" begin
	@testset "display" begin include("display.jl") end
	@testset "Constraints" begin include("testConstraints.jl") end
	@testset "Regular vs. Tensor" begin include("regularVStensor.jl") end
	@testset "testDiffOps.jl" begin include("testDiffOps.jl") end
	@testset "testInterpolationMatrix.jl" begin include("testInterpolationMatrix.jl") end
	@testset "testAvgOps.jl" begin include("testAvgOps.jl") end
	@testset "testBoundaryNodes.jl" begin include("testBoundaryNodes.jl") end
	@testset "testPersistency.jl" begin include("testPersistency.jl") end
	@testset "testIntPolyChain.jl" begin include("testIntPolyChain.jl") end
	@testset "testMassMatrices.jl" begin include("testMassMatrices.jl") end
	@testset "testForward.jl" begin include("testForward.jl") end
end


text:
module InteractiveFixedEffectModels
##############################################################################
##
## Dependencies
##
##############################################################################
using DataFrames
using FillArrays
using FixedEffects
using LeastSquaresOptim
using LinearAlgebra
using Printf
using Statistics
using StatsBase
using StatsFuns
using StatsModels
using Tables
using Vcov

using Reexport
@reexport using FixedEffectModels

##############################################################################
##
## Exported methods and types
##
##############################################################################

export InteractiveFixedEffectTerm,
InteractiveFixedEffectModel,
ife,
regife
##############################################################################
##
## Load files
##
##############################################################################
include("types.jl")
include("utils/tss.jl")
include("utils/formula.jl")

include("methods/gauss_seidel.jl")
include("methods/ls.jl")
include("fit.jl")
end


text:
############################################################################################
@testset "Batchsize computation: " begin
    Nchains = 10
    Nthreads = 11
    batchsize = BaytesSMC.compute_batchsize(Nchains, Nthreads)
    @test batchsize == 1

    Nchains = 11
    Nthreads = 10
    batchsize = BaytesSMC.compute_batchsize(Nchains, Nthreads)
    @test batchsize == 2

    Nchains = 10
    Nthreads = 10
    batchsize = BaytesSMC.compute_batchsize(Nchains, Nthreads)
    @test batchsize == 1
end


text:
"""
"""
mask

mask(n::NTuple{N}) where {N} =
    mask(ntuple(i -> dirichlet, Val(N)), n)

function mask(bc::NTuple{N}, n::NTuple{N}) where {N}
    opn = _mask.(bc, n)

    ntuple(Val(N)) do i

    end
    center = kron(first.(opn)...)

    center, face, edge, node =
end

function mask(bc::NTuple{1}, n::NTuple{1})
    opn = _mask.(bc, n)
    opn[1]
end

function mask(bc::NTuple{2}, n::NTuple{2})
    opn = _mask.(bc, n)
    ApplyArray(kron, opn[2], opn[1])
end

function mask(bc::NTuple{3}, n::NTuple{3})
    opn = _mask.(bc, n)
    ApplyArray(kron, opn[3], opn[2], opn[1])
end

_mask(::Dirichlet, n::Int) =
    (Diagonal(BitVector(i != n for i in 1:n)),
     Eye{Bool}(n))

_mask(::Periodic, n::Int) =
    (Eye{Bool}(n),
     Eye{Bool}(n))


text:
module FlujoFluidos

# Write your package code here.

export Reynolds

include("accesorios.jl")
include("fluidos.jl")
include("fricción.jl")
include("tuberias.jl")

end


text:
# missing methods in TaylorSeries.jl:

degree(f::Taylor1) = f.order
setindex!(f::Taylor1, x, i) = f.coeffs[i+1] = x

doc"""
A `Taylor1Model` is a Taylor model in 1 variable, providing a bounding box for a function as a polynomial $p$ plus an interval bound $\Delta$ for a function $f(x)$ of a single variable.

Fields:
- `n`: degree of the polynomial
- `x0`: expansion point
- `I`: interval over which the Taylor model is defined / valid
- `p`: the polynomial, represented as `TaylorSeries.Taylor1`
- `Δ`: the interval bound

The approximation is $f(x) \simeq \sum_i p_i (x - x_0)^i$.
"""
struct Taylor1Model{T,S}
    n::Int      # degree
    x0::Interval{T}  # expansion point
    I::Interval{T}   # interval over which the Taylor model is valid
    p::Taylor1{S}  # Taylor Taylor1nomial
    Δ::Interval{T}   # interval remainder bound
end

Taylor1Model(n, x0, I::Interval{T}, p, Δ) where{T} = Taylor1Model(n, interval(x0), I, Taylor1{Interval{T}}(p), interval(Δ))

# convenience constructors with same n, x0, I:
Taylor1Model(f, p, Δ) = Taylor1Model(f.n, f.x0, f.I, p, Δ)
Taylor1Model(f, Δ) = Taylor1Model(f, f.p, Δ)

include("arithmetic.jl")
include("functions.jl")
include("bound.jl")

include("integrate.jl")
include("draw.jl")

import Base.show
function show(io::IO, f::Taylor1Model)
    print(io,
    """Taylor1 model of degree $(f.n):
     - x0:  $(f.x0)
     -  I:  $(f.I)
     -  p: $(f.p)
     -  Δ:  $(f.Δ)
    """
    )
end

zero(f::Taylor1Model{T}) where {T} = Taylor1Model(f, zeros(length(f.p)), 0)

degree(f::Taylor1Model) = degree(f.p)

# Taylor1Model for a variable:
taylor1_var(n::Int, x0, I::T) where {T} = Taylor1Model(n, Interval(x0), I, Taylor1(T[x0, 1], n), T(0))

taylor1_var(n::Int, I::Interval) = taylor1_var(n, mid(I), I)

import Base.copy
copy(f::Taylor1Model) = Taylor1Model(f, copy(f.p), f.Δ)

"""
Evaluate a Taylor1Model at a point
"""
function (f::Taylor1Model)(x)

    x0 = f.x0

    if x ∈ f.I
        return (f.p)(x - x0) + f.Δ
    else
        throw(ArgumentError("Cannot evaluate Taylor1Model at point $x outside interval of definition $(f.I)"))
    end
end

"Note that TaylorSeries indices start at 0."
Base.getindex(f::Taylor1Model, i::Integer) = f.p[i]

# plot recipe for plotting 1D Taylor1Models
@recipe function g(f::Taylor1Model)

    x0, I, n, p, Δ = f.x0, f.I, f.n, f.p, f.Δ

    alpha --> 0.5
    seriestype := :shape

    xs = linspace(I.lo, I.hi, 100)
    evals = f.(xs)

    ylos = [y.lo for y in evals]
    yhis = [y.hi for y in evals]

    xs = [xs; reverse(xs); xs[1]]
    ys = [ylos; reverse(yhis); ylos[1]]

    xs, ys
end


text:
### Basic types and constructors for CompressedStack ###

## Compressed SubBlocks Signatures
type Signature{T}
  first::Int
  last::Int
  context::T
  pos::Int # position of the input
end

## A Partially Compressed Block is composed of the signatures of its SubBlocks
typealias Block{T} Vector{Signature{T}}
## Each level of compressed Blocks (first and second) are stored in Levels
typealias Levels{T} Vector{Block{T}}
## Type Data to store the data and their index
type Data{D}
  data::D
  index::Int
end

## Constructor for Signature
# Constructor for singleton signature
function Signature{T}(index::Int, context::T, pos::Int)
  Signature(index, index, context, pos)
end
# Constructor for a Signature (from a Block)
function Signature{T}(block::Block{T})
  context = block[1].context
  first = block[1].first
  pos = block[1].pos
  last = block[end].last
  Signature(first, last, context, pos)
end

## General CompressedStack (i.e. with parametric context T and data type D)
# D can also be used as an index when the data type is too big
type CompressedStack{T,D}
  # Structure constraints
  size::Int # Size of the input in #elements
  space::Int # Maximum space order of the compressed stack
  depth::Int # Depth (#levels of compression) based on size and space
  # IO
  input::IOStream # Pointer on a file to avoid memory consumption
  pos::Int # Position of the input before the reading
  output::Nullable{IOStream} # Pointer to an optional (Nullable) output file
  # Functions defining the behavior of the stack in case of push/pop
  push_condition::Function
  push_action::Function
  pop_condition::Function
  pop_action::Function
  # First Blocks
  first_partial::Levels{T} # Levels of partially compressed blocks
  first_explicit::Vector{Data{D}}
  first_sign::Nullable{Signature{T}}
  # Second Blocks
  second_partial::Levels{T} # Levels of partially compressed blocks
  second_explicit::Vector{Data{D}}
  second_sign::Nullable{Signature{T}}
  # Fully Compressed Block (only a signature possibly empty [Nullable])
  compressed::Nullable{Signature{T}}
  # Stack's running information
  index::Int # Stock the maximum index read in the input file
  context::Nullable{T} # Current context to use while making block signatures
  # types
  context_type::DataType
  data_type::DataType
end

# General Constructor for CompressedStack
function CompressedStack(size::Int, space::Int, input::IOStream,
  context_type::DataType, data_type::DataType,
  push_action::Function, push_condition::Function,
  pop_action::Function, pop_condition::Function;
  pos = position(input),
  index = 0, context = Nullable{context_type}(), output = Nullable{IOStream}())

  depth = convert(Int,ceil(log(space, size - 0.1))) - 1
  compressed = Nullable{Signature{context_type}}()
  first_explicit = Vector{Data{data_type}}()
  first_partial = Levels{context_type}()
  second_explicit = Vector{Data{data_type}}()
  second_partial = Levels{context_type}()
  sign_explicit = compressed

  # Initialization of each Block at each level
  for i in 1:(depth-1)
    push!(first_partial,Vector{Signature{context_type}}())
    push!(second_partial,Vector{Signature{context_type}}())
  end

  # Call to the basic constructor
  CompressedStack(size, space, depth, input, pos, output, push_condition,
  push_action, pop_condition, pop_action, first_partial, first_explicit,
  sign_explicit, second_partial, second_explicit, sign_explicit,
  compressed, index, context, context_type, data_type)
end

# Constructor for the reconstruction procedure
function CompressedStack(stack::CompressedStack, size::Int, context, index::Int,
  pos::Int)

  input = stack.input
  seek(input, pos)

  CompressedStack(size * stack.space, stack.space, input, stack.context_type,
  stack.data_type, stack.push_action, stack.push_condition,
  stack.pop_action, stack.pop_condition;
  index = index, context = Nullable(context), output = Nullable{IOStream}())
end

### Constructor for CompressedStack from a file input ##
## No given output file, $(name)_out is generated
# Size and space given in the input file
function CompressedStack(name::AbstractString, pop_action::Function,
  push_action::Function, pop_condition::Function, push_condition::Function,
  context_type::DataType, data_type::DataType)

  input = open(name, "r")
  output_name = name * "_out"
  output = open(output_name, "w")
  settings = get_settings(input)

  CompressedStack(settings[1], settings[2], input, context_type, data_type,
  push_action, push_condition, pop_action, pop_condition;
  output= Nullable{IOStream}(output))
end

# Size and space given by the user
function CompressedStack(name::AbstractString, pop_action::Function,
  push_action::Function, pop_condition::Function, push_condition::Function,
  context_type::DataType, data_type::DataType, size::Int, space::Int)

  input = open(name, "r")
  output_name = name * "_out"
  output = open(output_name, "w")
  settings = get_settings(input)

  CompressedStack(size, space, input, context_type, data_type,
  push_action, push_condition, pop_action, pop_condition;
  output= Nullable{IOStream}(output))
end


text:
using Main.JuliaPkgsList: getMetaDataValue, getName
using Main.JuliaPkgsList: getRepo, getStarCount, getLatestVersion
using Main.JuliaPkgsList: loadPkgsData, sortStarred

@testset "JuliaPkgsList.jl :: getMetaDataValue" begin
    FOO = "foo"
    m1 = Dict()
    m2 = Dict(FOO => 5)
    m3 = Dict(FOO => nothing)
    d = map(m -> Dict(JuliaPkgsList.METADATA_KEY => m), [m1, m2, m3])
    @test getMetaDataValue(d[1], FOO, 66) == 66
    @test getMetaDataValue(d[1], FOO, -7) == -7
    @test getMetaDataValue(d[2], FOO, 66) == 5
    @test getMetaDataValue(d[3], FOO, 66) == 66
end

@testset "JuliaPkgsList.jl :: pkg-info extraction" begin
    pkgs = loadPkgsData(SAMPLE_FILE)
    @test length(pkgs) == 148
    juliaPkg = pkgs[105]
    @test getName(juliaPkg) == "julia"
    @test getRepo(juliaPkg) == JULIA_REPO
    @test getStarCount(juliaPkg) == 31228
    @test getLatestVersion(juliaPkg) == "1.5.3"
    @test JuliaPkgsList.getNameAndVersion(juliaPkg) == "julia, 1.5.3"
    @test JuliaPkgsList.getRepoAndVersion(juliaPkg) == "$(JULIA_REPO), 1.5.3"
    yieldMacrosPkg = pkgs[20]
    @test getName(yieldMacrosPkg) == "AbstractYieldMacros"
    @test getStarCount(yieldMacrosPkg) == 0
    @test getLatestVersion(yieldMacrosPkg) == "0.1.0"
    # sorting
    sorted = sortStarred(pkgs)
    @test sorted[1] == juliaPkg
end


text:
using Helloworld
using Documenter

makedocs(;
    modules=[Helloworld],
    authors="Keishi Sando <nobuta05@gmail.com> and contributors",
    repo="https://github.com/nobuta05/Helloworld.jl/blob/{commit}{path}#L{line}",
    sitename="Helloworld.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://nobuta05.github.io/Helloworld.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/nobuta05/Helloworld.jl",
)


text:
# Test cmap.jl

# Tests could be more comprehensive (of course).  In some cases all
# that is done is to execute the function in some way just to make
# sure no exceptions are thrown as a result of something breaking with
# a new version of Julia.

println("testing cmap")

# cmap: Test cmap by calling it with no arguments. This generates a
# list of all the colour maps and in the process generates them all.
cmap()

#=  See if this is what causes Travis to fail...
# equalisecolourmap: Generate a colour map in Lab space with an uneven
# ramp in lightness and check that this is corrected
rgblab = "LAB"
labmap = zeros(256,3)
labmap[1:127,1] = linspace(0,20,127)
labmap[128:256,1] = linspace(20,100,129)
formula = "CIE76"
W = [1,0,0]
sigma = 1
rgbmap = equalisecolourmap(rgblab, labmap, formula, W, sigma)
# Convert to Nx3 array and then back to lab space. Then check that dL
# is roughly constant
labmap2 = srgb2lab(rgbmap)
dL = gradient(labmap2[:,1])
@test maximum(dL[2:end-1]) - minimum(dL[2:end-1]) < 1e-1
=#

# linearrgbmap
rgbmap = linearrgbmap([1,1,1],99)
rgbmap2 = PerceptualColourMaps.convert(Array{Float64,2},rgbmap)
# check middle colour is [.5 .5 .5]
@test maximum(abs(rgbmap2[50:50,:] - [0.5 0.5 0.5])) < 1e-6

# UInt32colormap
ui32 = UInt32colormap(cmap("L3")[1])

# srgb2lab and lab2srgb
rgb = [1 0 0
       0 1 0
       0 0 1]

# lab values of rgb values above
lab = [ 53.2377   80.1029   67.2144
        87.7346  -86.1606   83.1862
        32.2966   79.2009 -107.8615]

clab = srgb2lab(rgb)
crgb = lab2srgb(clab)
@test maximum(abs(rgb-crgb)) < 1e-1
@test maximum(abs(lab-clab)) < 1e-1


text:
module HPInt

push!(LOAD_PATH, ".")
using Error
using Lexer
using Images
using Cairo
export parse, calc, NumVal, ClosureVal, analyze, MatrixVal

#
# =================================================
#

abstract OWL
abstract Environment
abstract RetVal
type Num <: OWL
	n::Real
end
type NumVal <: RetVal
     	n::Real
end
type Binop <: OWL
	op::Function
	lhs::OWL
	rhs::OWL
end
type Plus <: OWL
	operands::Array{OWL}
end
type Unop <: OWL
	op::Function
	operand::OWL
end
type Id <: OWL
     	name::Symbol
end
type If0 <: OWL
     	condition::OWL
	zero_branch::OWL
	nonzero_branch::OWL
end
type Binder <: OWL
     	name::Id
	binding_expr::OWL
end
type With <: OWL
     	binders::Array{Binder}
	body::OWL
end
type FunDef <: OWL
     	formal_parameters::Array{Id}
	fun_body::OWL
end
type FunApp <: OWL
     	fun_expr::OWL
	arg_exprs::Array{OWL}
end
type ClosureVal <: RetVal
     	params::Array{Id}
	body::OWL
	env::Environment
end
type And <: OWL
	operands::Array{OWL}
end
type mtEnv <: Environment
end
type SymVal <: OWL
     	name::Id
	value::Any
end
type CEnvironment <: Environment
     	symvals::Array{SymVal}
	parent::Environment
end
type SimpleLoad <: OWL
     	location::AbstractString
end
type SimpleSave <: OWL
	matrix::OWL
	location::AbstractString
end
type RenderText <: OWL
     	words::AbstractString
	xpos::OWL
	ypos::OWL
end
type Emboss <: OWL
     	matrix::OWL
end
type DropShadow <: OWL
     	matrix::OWL
end
type InnerShadow <: OWL
     	matrix::OWL
end
type MatrixVal <: RetVal
     	n::Array{Float32,2}
end

#
# =================================================
#

# The collatz primitive
function collatz(n::Real)
	return collatz_helper(n, 0)
end
function collatz_helper(n::Real, num_iters::Int)
	if n == 1
	      	return num_iters
	end
	if mod(n,2)==0
		return collatz_helper(n/2, num_iters+1)
	else
		return collatz_helper(3*n+1, num_iters+1)
	end
end
# The simple_load primitive
function simple_load(img_path::AbstractString)
  	im = Images.load(img_path);
  	tmp = Images.permuteddimsview(channelview(im),(2,3,1));
  	d = Images.data(ImageMeta(tmp));
  	d = d[:,:,1];  # just the r channel
  	d = convert(Array{Float32,2}, d);
  	return d
end
# The simple_save primitive
function simple_save(output::Array, img_path::AbstractString)
    	output[output .> 1.0] = 1.0
    	output[output .< 0.0] = 0.0
    	tmpc = convert(Array{UInt32,2}, floor(output*255.0))
    	tmp_output =  tmpc + tmpc*256 + tmpc*65536 + 0xff000000
    	c2 = CairoImageSurface(transpose(tmp_output), Cairo.FORMAT_ARGB32)
    	write_to_png(c2, img_path)
    	return 42
end
# The render_text primitive
function render_text(text_str::AbstractString, xpos, ypos)
  	data = Matrix{UInt32}(256, 256);
  	c = CairoImageSurface(data, Cairo.FORMAT_ARGB32);
  	cr = CairoContext(c);
  	set_source_rgb(cr, 1., 1., 1.);
  	rectangle(cr, 0., 0., 256., 256.);
  	fill(cr);
  	set_source_rgb(cr, 0, 0, 0);
  	select_font_face(cr, "Sans", Cairo.FONT_SLANT_NORMAL, Cairo.FONT_WEIGHT_BOLD);
  	set_font_size(cr, 90.0);
  	move_to(cr, xpos, ypos);
  	show_text(cr, text_str);
  	# tmp is an Array{UInt32,2}
  	tmp = cr.surface.data;
  	# grab just the blue channel, and convert the array to an array of floats
  	tmp2 = convert(Array{Float32,2}, tmp & 0x000000ff) / 255.0;
  	tmp2 = convert(Array{Float32,2}, tmp2);

  	return tmp2
end
# The emboss primitive
function emboss(img::Array)
  	f = [ -2. -1. 0.
              -1.  1. 1.
               0.  1. 1. ];
  	f = convert(Array{Float32,2}, f);

  	es = conv2(f, img);
  	es = es[1:256,1:256];
  	return es
end
# The drop_shadow primitive
function drop_shadow(img::Array)
  	foo = convert(Array{Float32,2}, Kernel.gaussian((5.0,5.0),(25,25))[-12:12,-12:12]);
  	foo = foo / maximum(foo);
  	ds = conv2(foo, img);
  	ds = ds[13:256+12,13:256+12];
  	ds = ds / sum(foo);
  	return ds
end
# The inner_shadow primitive
function inner_shadow( img::Array )
  	foo = convert(Array{Float32,2}, Kernel.gaussian((5.0,5.0),(25,25))[-12:12,-12:12]);
  	foo = foo / maximum(foo);
  	is = conv2(foo, 1.0-img);
  	is = is[8:251+12,8:251+12];
  	is = is / sum(foo);
  	is = max(is, img);
  	return is
end

#
# =================================================
#

# All operators will be stored in this dicitonary
op_table = Dict([
			(:+,+),
			(:-,-),
			(:*,.*),
			(:/,./),
			(:mod,mod),
			(:collatz,collatz),
      			(:min, min),
			(:max, max)
		])
reserved_words=[:+,:-,:*,:/,:.*,:./,:mod,:collatz,:min,:max,:if0,:with,:lambda, :and, :simple_load, :simple_save, :render_text, :emboss, :drop_shadow, :inner_shadow] # added new reserved word

#
# =================================================
#

# This is the 'base case' of the parse tree.
function parse(expr::Real)
	return Num(expr)
end
# Parse things with type Id
function parse(expr::Symbol)
	if !(expr in reserved_words) # Make sure the word is not reserved
		return Id(expr)
	else
	   	throw(LispError("Symbol \'$expr\' is a reserved word"))
	end
end
# When an array is encountered use this
function parse(expr::Array{Any})
	# Check to make sure array is not empty
	if length(expr) >= 1
	   	# The operator should be the first argument
		op_symbol = expr[1]
		if haskey(op_table, op_symbol) # Check to make sure this is a valid op.
		   	if length(expr) == 3 # All Operators have 3 or 2 length => 2 or 1 arguments respectively.
			   	if op_symbol != :collatz # Only the collatz operator needs 1 argument.
				   	return Binop(op_table[op_symbol], parse(expr[2]), parse(expr[3]))
				else
					throw(LispError("Too many parameters given for \'$op_symbol\'")) # Throws if 'collatz' is called with more than one argument.
				end
			elseif length(expr) == 2
			       	if op_symbol == :- || op_symbol == :collatz # These are the only two allowed 1-arg symbols
				   	return Unop(op_table[op_symbol], parse(expr[2]))
				else
					throw(LispError("Not enough parameters given for \'$op_symbol\'")) # Throw this if not uniary operator
				end
			elseif length(expr) >= 3 && op_symbol == :+ # Handle new grammer with multiple sums
			        ops = OWL[]
			        for i in 2:size(expr,1)
			       	     	push!(ops, parse(expr[i]))
			        end
			        return Plus(ops)
			else
				throw(LispError("Wrong amount of parameters given for \'$op_symbol\'")) # Throw this if not correct arg length
			end
	  	elseif op_symbol == :and # new grammer rule for 'and'
			if length(expr) >= 3 # Must have more than two parameters
			   	ops = OWL[]
				for i in 2:size(expr,1)
			      	      	push!(ops, parse(expr[i]))
				end
				return And(ops)
			else
				throw(LispError("Cannot perform 'and' operation without 2 or more operands."))
			end
			# Handle non-binary/uniary operators and AST elements seperately
		elseif op_symbol==:if0
		        if length(expr)==4 # Verify correct input length
				return If0(parse(expr[2]), parse(expr[3]), parse(expr[4])) # Good Input gives good return
				# In case it isn't clear, 'If0' is a node in the AST and has a type defined for it.
			else
				throw(LispError("Wrong amount of parameters given for \'If0\'"))
			end
		elseif op_symbol == :with
	       	        if length(expr) == 3 # Verify correct argument length
		   	   	binding_exprs=expr[2] # This should be an array of 'pairs'
				if no_dups(binding_exprs) # Check to make sure there's no duplicates in the binders
					binders=Binder[]
					if length(binding_exprs) == 0
					   	# Do nothing
					elseif typeof(binding_exprs[1]) == Array{Any,1}
					       	for i in 1:size(binding_exprs,1) # Go through list and bind values
						      	if length(binding_exprs[i]) == 2
				      	       	       	  	push!(binders, Binder(parse(binding_exprs[i][1]), parse(binding_exprs[i][2]))) # Evaluate/parse all bindings
							else
								throw(LispError("Must have one symbol per one expression"))
							end
						end
					else
						throw(LispError("Incorrect syntax on with statement"))
					end
					return With(binders, parse(expr[3])) # The third element is the actual body
					# We'll handle the 'With' node in the AST in the calc statement
				else
					throw(LispError("Duplicate symbols in \'with\' statement"))
				end
			else
				throw(LispError("Wrong amount of parameters given for \'with\'"))
			end
		elseif op_symbol == :lambda
	       	        if length(expr) == 3 # Verify correct argument length
		   	   	if no_dups(expr[2]) # Make sure the arguments aren't duplicated
					return FunDef(map(parse, expr[2]),parse(expr[3])) # elements: 2- argument list. 3- function body
				else
					throw(LispError("Duplicate symbols in \'lambda\' expression"))
				end
	       		else
				throw(LispError("Wrong amount of parameters given for \'lambda\'"))
			end
		elseif op_symbol==:simple_load
			if length(expr) == 2 && typeof(expr[2]) == String # Verify correct arguments
			   	return SimpleLoad(expr[2]) # Easy parse I know
			else
				throw(LispError("Invalid syntax for simple_load."))
			end
		elseif op_symbol == :simple_save
		        if length(expr) == 3 && typeof(expr[3]) == String # Verify correct arguments
				return SimpleSave(parse(expr[2]),expr[3]) # Easy parse I know
			else
			   	throw(LispError("Invalid syntax for simple_save"))
			end
		elseif op_symbol == :render_text
		       	if length(expr) == 4 && typeof(expr[2]) == String # Verify correct arguments
			   	return RenderText(expr[2], parse(expr[3]), parse(expr[4])) # Easy parse I know
			else
				throw(LispError("Invalid syntax for render_text. Length = $(length(expr)) and type = $(typeof(expr[2]))"))
			end
		elseif op_symbol == :emboss
		       	if length(expr) == 2 # Verify correct argument length
			   	return Emboss(parse(expr[2])) # Easy parse I know
			else
				throw(LispError("Invalid syntax for emboss."))
			end
		elseif op_symbol == :drop_shadow
		       	if length(expr) == 2 # Verify correct argument length
			   	return DropShadow(parse(expr[2])) # Another easy parse I know
			else
				throw(LispError("Invalid syntax for drop_shadow."))
			end
		elseif op_symbol == :inner_shadow
		       	if length(expr) == 2 # Verify correct argument length
			   	return InnerShadow(parse(expr[2])) # Easy parse I know
			else
				throw(LispError("Invalid syntax for inner_shadow."))
			end
		else # Everything else must evaluate to a function
			arg_exprs=OWL[]
			for i in 2:size(expr,1)
		      	      	push!(arg_exprs, parse(expr[i]))
			end
			return FunApp(parse(op_symbol), arg_exprs) # The op symbol isn't a symbol. It's a function expression
		end
	else
		throw(LispError("No parameters given: Empty array."));
	end
end
# Catch all other parse types.
function parse(expr::Any)
	throw(LispError("Invalid type \'$expr\'"))
end
# This is the helper function that checks for duplicate symbols in with and lambda statements
function no_dups(exprs::Array{Any})
	if length(exprs) > 0
  	   	syms = Any[]
		if typeof(exprs[1]) == Array{Any,1} # With nodes should all have this
		   	if length(exprs[1]) != 0 # Empty lists are okay. Handle seperately
    	   		   	for i in 1:size(exprs,1)
			       	      	if typeof(exprs[i][1]) != Symbol
				   	   	throw(LispError("Must use symbols for binding expressions"))
					elseif in(exprs[i][1], reserved_words) # Check for reserved words
				               	throw(LispError("Reserved word \'$exprs[i][1]\' cannot be used in binding expression."))
      		      	 		elseif in(exprs[i][1], syms) # Check for duplicates of partial forward-list
        		    	      	       	return false # Duplicate encountered
					else
						push!(syms, exprs[i][1]) # Push good symbols
			 		end
				end
			end
			return true
		elseif typeof(exprs[1]) == Symbol # Used in lambda nodes
    	       	        for i in 1:size(exprs,1)
      		      	      	if typeof(exprs[i]) != Symbol # Make sure we're using symbols
        		    	   	throw(LispError("Must use symbols for single-dim argument lists."))
				elseif in(exprs[i], reserved_words) # Check for reserved words
				        throw(LispError("Reserved word \'$exprs[i]\' cannot be used in argument lists."))
				elseif in(exprs[i], syms) # Check for duplicates
			    	        return false # Duplicate encountered
			 	else
					push!(syms, exprs[i]) # All is good in the hood.
			 	end
			end
			return true
		else
			throw(LispError("Improper format. Probably missing parentheses."))
		end
	else
    	   	return true # No need to throw error or even check. Empty lists are okay.
  	end
end
# Default behavior to catch incorrect format. Throw error if it's not a list
function no_dups(exprs::Any)
	throw(LispError("Improper format. Array expected."))
end

#
# =================================================
#

# Analyze number structure
function analyze(owl::Num)
	return owl
end
# Analyze binomial operator tree
function analyze(owl::Binop)
	return Binop(owl.op, analyze(owl.lhs), analyze(owl.rhs)) # Pass on the recursion
end
# Handle plus nodes
function analyze(owl::Plus)
	if length(owl.operands) == 2 # Is this binomial?
		return Binop(+, analyze(owl.operands[1]), analyze(owl.operands[2]))
	else
		return Binop(+, analyze(owl.operands[1]), analyze(Plus(owl.operands[2:end])))
	end
end
# Handle Uni-operator trees
function analyze(owl::Unop) # Just pass on the recursion
	return Unop(owl.op, analyze(owl.operand))
end
# Handle id nodes
function analyze(owl::Id) # Just pass on the recursion
	return owl
end
# Handle if0 nodes
function analyze(owl::If0) # Just pass on the recursion
	return If0(analyze(owl.condition), analyze(owl.zero_branch), analyze(owl.nonzero_branch))
end
# Handle with nodes
function analyze(owl::With)
	syms = Id[]
	exprs = OWL[]
	for i in 1:size(owl.binders,1)
	      	push!(syms, owl.binders[i].name)
		push!(exprs, owl.binders[i].binding_expr)
	end
	return FunApp(FunDef(syms, analyze(owl.body)), map(analyze, exprs))
end
# Handle function definition nodes
function analyze(owl::FunDef) # Just pass on the recursion
	return FunDef(owl.formal_parameters, analyze(owl.fun_body))
end
# Handle function application nodes
function analyze(owl::FunApp) # Just pass on the recursion
	return FunApp(analyze(owl.fun_expr), map(analyze, owl.arg_exprs))
end
# Handle and nodes
function analyze(owl::And)
	if length(owl.operands) == 1 # It shouldn't ever return one on a root call. It will on recursion though.
		return If0(analyze(owl.operands[1]), Num(0), Num(1))
	else
		return If0(analyze(owl.operands[1]), Num(0), analyze(And(owl.operands[2:end])))
	end
end
# Handle simple load nodes
function analyze(owl::SimpleLoad)
	return owl # Leaf node
end
# Handle simple save nodes
function analyze(owl::SimpleSave)
	return SimpleSave(analyze(owl.matrix), owl.location) # Just pass on the recursion
end
# Handle render text nodes
function analyze(owl::RenderText)
	return RenderText(owl.words, analyze(owl.xpos), analyze(owl.ypos)) # Just pass on the recursion
end
# Handle emboss nodes
function analyze(owl::Emboss)
	return Emboss(analyze(owl.matrix)) # Just pass on the recursion
end
# Handle drop shadow nodes
function analyze(owl::DropShadow)
	return DropShadow(analyze(owl.matrix)) # Just pass on the recursion
end
# Handle inner shadow nodes
function analyze(owl::InnerShadow)
	return InnerShadow(analyze(owl.matrix)) # Just pass on the recursion
end
# Handle default function call
function analyze(owl::OWL) # Should never get here. Throw error.
	throw(LispError("Unkown node!")) # Just pass on the recursion
end

#
# =================================================
#

# Calculate the basic OWL/root abstract syntax tree
function calc(ast::OWL)
	 return calc(ast, mtEnv())
end
# Base case is a number or closure which we will return
function calc(owl::Num, env::Environment)
	 return NumVal(owl.n)
end
# Handle binary operations as in the last lab
function calc(owl::Binop, env::Environment) # The environment will just be used on calculating either side.
	 lhs=@spawn calc(owl.lhs, env)
	 rhs=@spawn calc(owl.rhs, env)
	 left = fetch(lhs) # This is for threading both sides
	 right = fetch(rhs)
	 if ((owl.op != ./) || right.n != 0) && ((owl.op != /) || right.n != 0)
	    	 if typeof(left) == NumVal && typeof(right) == NumVal # Make sure things can evaluated in the calcs
		    	return NumVal(owl.op(left.n, right.n)) # Return correct calculation
		 # Handle new data types
		 elseif (typeof(left) == NumVal && typeof(right) == MatrixVal) || (typeof(left) == MatrixVal && typeof(right) == NumVal)
		 	return MatrixVal(owl.op(left.n, right.n))
		 # Handle new data types
		 elseif typeof(left) == MatrixVal && typeof(right) == MatrixVal
	       	 	if length(left.n) == length(right.n) && length(left.n[1]) == length(right.n[1]) # Check correct sizes
			   	if owl.op == *
					return MatrixVal(left.n .* right.n)
				elseif owl.op == /
					return MatrixVal(left.n ./ right.n)
				else
			   	   	return MatrixVal(owl.op(left.n, right.n))
				end
			else
			   	throw(LispError("Cannot perform binary operations on matrices of different sizes."))
			end
		 else
			throw(LispError("Invalid binary op types did not evaluate correctly"))
		 end
	else
		throw(LispError("Cannot divide by zero."))
	end
end
# Handle bnary operations in a similar fashion to binary operations
function calc(owl::Unop, env::Environment) # As in Binop, the environment is just used to calculate the expression recursively
	 unary = calc(owl.operand, env) # Calculate the expression
	 if typeof(unary) == NumVal # Make sure it evaluated to a value
		 if !(owl.op == collatz && unary.n <= 0) # Check to make sure the collatz function is called with a correct input.
			  return NumVal(owl.op(unary.n))
		 else
			  throw(LispError("Cannot perform collatz on \'$unary\'"))
		 end
	 else
		 throw(LispError("Invalid uninary op types did not evaluate correctly"))
	 end
end
# Handle If0 nodes
function calc(owl::If0, env::Environment) # Pass in the environment as always
	 cond = calc(owl.condition, env)
	 if typeof(cond) == NumVal # Make sure that the condition calculated correctly
	    	 if cond.n == 0
		    	  return calc(owl.zero_branch, env) # This is the 'then' branch
		 else
			  return calc(owl.nonzero_branch, env) # This is the 'else' branch
		 end
	 else
		 throw(LispError("Invalid \'If0\' type with value \'$cond\'."))
	 end
end
# Handle Id/Symbols in grammer
function calc(owl::Id, env::Environment)
	 if env == mtEnv() # Check to see if it's an mpty environment.
	    	 throw(LispError("Could not find symbol \'$owl\'"))
	 elseif hasSymVal(env.symvals, owl) != -1 # See helper function
	 	 return getValue(env.symvals, hasSymVal(env.symvals, owl)) # Return symbol value.
	 else
		 return calc(owl, env.parent) # I believe I need this for recursion
	 end
end
# Handle With nodes
function calc(owl::With, env::Environment)
	 symvals = SymVal[] # The With node has a bunch of sumbols we need to evaluate.
	 for i in 1:size(owl.binders,1)
	       	 push!(symvals, SymVal(owl.binders[i].name, calc(owl.binders[i].binding_expr, env)))
	 end
	 return calc(owl.body, CEnvironment(symvals, env)) # Extend the environment and pass it on.
end
# Define the FunDef Node and pass of the closure
function calc(owl::FunDef, env::Environment)
	 return ClosureVal(owl.formal_parameters, owl.fun_body, env) # Super simple I know.
end
# This is the FunApp node
function calc( owl::FunApp, env::Environment )
	 # The function expression should result in a ClosureVal
	 ret_closure = calc(owl.fun_expr, env)
	 if typeof(ret_closure) == ClosureVal # Check for correct evaluation
	    	 if length(owl.arg_exprs) == length(ret_closure.params) # These two lists should be equal.
	 	    	  # Extend the current environment by binding the actual parameters to the formal parameters
	 	 	  symvals = SymVal[] # Build array
			  for param in zip(ret_closure.params, owl.arg_exprs)
			      	  push!(symvals, SymVal(param[1], @spawn calc(param[2], env)))
	 	 	  end
	      	 	  return calc(ret_closure.body, CEnvironment(symvals, ret_closure.env)) # Calculate new environment and return
		 else
		    	  throw(LispError("Number of parameters and arguments does not match in function call/declaration."))
	 	 end
	 else
		 throw(LispError("Invalid type did not return closure"))
	 end
end
# Calculate SimpleLoad nodes
function calc(owl::SimpleLoad, env::Environment)
	return MatrixVal(simple_load(owl.location)) # Construct primitive off of loader function
end
# Calculate SimpleSave nodes
function calc(owl::SimpleSave, env::Environment)
	matrix = calc(owl.matrix, env)
	if typeof(matrix) == MatrixVal # Check correct input
	   	return simple_save(matrix.n, owl.location)
	else
		throw(LispError("Expected MatrixVal in simple_save calculation. Got \'$(typeof(matrix))\' instead."))
	end
end
# Calculate RenderText nodes
function calc(owl::RenderText, env::Environment)
	xpos=@spawn calc(owl.xpos, env)
	ypos=@spawn calc(owl.ypos, env)
	xpos = fetch(xpos) # Thread this too.
	ypos = fetch(ypos)
	if typeof(xpos) == NumVal && typeof(ypos) == NumVal # Check correct input
	   	return MatrixVal(render_text(owl.words, xpos.n, ypos.n))
	else
		throw(LispError("Expected NumVal in render_text."))
	end
end
# Calculate Emboss nodes
function calc(owl::Emboss, env::Environment)
	matrix = calc(owl.matrix, env)
	if typeof(matrix) == MatrixVal # Check correct input
	   	return MatrixVal(emboss(matrix.n))
	else
		throw(LispError("Expected MatrixVal in emboss."))
	end
end
# Calculate DropShadow nodes
function calc(owl::DropShadow, env::Environment)
	matrix = calc(owl.matrix, env)
	if typeof(matrix) == MatrixVal # Check correct input
	   	return MatrixVal(drop_shadow(matrix.n))
	else
		throw(LispError("Expected MatrixVal in drop_shadow."))
	end
end
# Calculate InnerShadow nodes
function calc(owl::InnerShadow, env::Environment)
	matrix = calc(owl.matrix, env)
	if typeof(matrix) == MatrixVal # Check correct input
	   	return MatrixVal(inner_shadow(matrix.n))
	else
		throw(LispError("Expected MatrixVal in inner_shadow."))
	end
end
# Default calc case
function calc(owl::Any)
	 throw(LispError("Cannot calculate unknown type"))
end
# This could also happen with environments.
function calc(owl::Any, env::Environment)
	 throw(LispError("Cannot calculate unknown type"))
end
# Helper function to check for symbol name in environment array
function hasSymVal(symvals::Array{SymVal}, name::Id)
	for i in 1:size(symvals,1) # Iterate through and find
	       	 if symvals[i].name.name == name.name # I added some parameters to the AST nodes
	       	    	 return i # Found
		 end
	end
	return -1 # Not found
end
# Helper function to get value associated with symbol name
function getValue(symvals::Array{SymVal}, index::Int)
	 return fetch(symvals[index].value)
	 throw(LispError("Symbol with \'$name\' not found - Good job on getting here though"))
end

#
# =================================================
#

# This is for debugging
function calc(expr::AbstractString)
	 return calc(analyze(parse(Lexer.lex(expr))))
end

#
# =================================================
#

end # module

max_stars_repo_name: scottcrossen/JuliaDevelopment

max_issues_repo_name: scottcrossen/JuliaDevelopment

max_forks_repo_name: scottcrossen/JuliaDevelopment


text:
#
println("Fb) Symbolic simplification of Racah expressions by means of sum rules.")
#
using SymEngine

if  true
    # Test for sum rules
    rex = RacahAlgebra.selectRacahExpression(14);        println("rex-original        = $rex")
    rex = RacahAlgebra.equivalentForm(rex);             println("rex-equivalent      = $rex")
    wa  = RacahAlgebra.evaluate(rex);                   println("wa-sum rule eval    = $wa")
    ## ws  = [:j => 2//1, :m => 1, :jb => 2, :jc => 1]
    ## wb  = RacahAlgebra.subs(rex, ws);                   println("\nrex-substitution  = $wb")
    ## wc  = RacahAlgebra.subs(wa, ws);                    println("\nrhs-substitution  = $wc")
    ## wd  = RacahAlgebra.evaluateNumerical(wb);           println("lhs-numerical       = $wd")
    ## we  = RacahAlgebra.evaluateNumerical(wc);           println("rhs-numerical       = $we")
end


text:
module ForecastEval
#-----------------------------------------------------------
#PURPOSE
#	Colin T. Bowers module for forecast evaluation
#NOTES
#-----------------------------------------------------------

using 	Statistics,
		StatsBase,
		Distributions,
		DependentBootstrap

import 	Base.show

export 	DMMethod,
		DMBoot,
		DMHAC,
		DMTest,
		dm,
		RCMethod,
		RCBoot,
		RCTest,
		rc,
		SPAMethod,
		SPABoot,
		SPATest,
		spa,
		MCSMethod,
		MCSBoot,
		MCSBootLowRAM,
		MCSTest,
		mcs

include("hacvariance.jl")
include("pvaluelocal.jl")
include("dm.jl")
include("rc.jl")
include("spa.jl")
include("mcs.jl")
include("mcs_lowRAM.jl")

end #module


text:
using multimode
using QuantumOptics

# Particles
const E0 = 0.2
const Nparticles = 1
const Nlevels = 14
const potential = multimode.BoxPotential
const particles = multimode.Bosons(Nparticles, Nlevels, E0, potential)

# Cavity mode
const index1 = 12
const cavitymodes1 = 14
const delta1 = -4.
# const eta1 = 2.0
const kappa1 = 1.
const U0_1 = -2.

# Multimode system
const scaling = 1./4

function spectrum(system, T)
    Ψ₀ = basis_ket(multimode.quantum.basis(system), 1)
    ρ₀ = Ψ₀⊗dagger(Ψ₀)
    Tss = [0,100]
    tout, ρ_t = multimode.quantum.timeevolution_master(system, Tss, ρ₀)
    ρ_ss = ρ_t[end]

    a = embed(multimode.quantum.basis(system), [2], [destroy(multimode.quantum.basis(system.modes[1]))])
    aρ_ss = a*ρ_ss

    x = Complex128[]
    fout(t, rho) = push!(x, expect(dagger(a), rho))

    multimode.quantum.timeevolution_master(system, T, aρ_ss; fout=fout)
    N = length(x)
    x = [reverse(x), x]
    S = abs2(fft(x))
    S = [S[N+1:end], S[1:N]]
    S /= maximum(S)
    return S
end

T = [0:0.1:250]
#etas = [0.5:0.1:3.5]
etas = [0.1:0.5:3.5]
#etas = [2.0:0.1:2.0]
S = zeros(Float64, length(etas), 2*length(T))
for (i, eta1) in enumerate(etas)
    println("η=", eta1)
    mode1 = multimode.CavityMode(index1, cavitymodes1, delta1, eta1, kappa1, U0_1)
    system = multimode.MultimodeSystem(scaling, particles, [mode1])
    S[i,:] = spectrum(system, T)
end

writecsv("spectrum_mode12_b.csv", S)

# using PyCall
# @pyimport matplotlib.pyplot as plt
# plt.plot(S)
# plt.show()


text:
@test factorial(7) == 5040
@test factorial(7,3) == 7*6*5*4
@test binomial(5,3) == 10
p = shuffle([1:1000])
@test isperm(p)
@test all(invperm(invperm(p)) .== p)
push!(p, 1)
@test !isperm(p)
a = randcycle(10)
@test ipermute!(permute!([1:10], a),a) == [1:10]
@test collect(combinations("abc",2)) == ["ab","ac","bc"]
@test collect(permutations("abc")) == ["abc","acb","bac","bca","cab","cba"]
@test collect(filter(x->(iseven(x[1])),permutations([1,2,3]))) == {[2,1,3],[2,3,1]}
@test collect(filter(x->(iseven(x[3])),permutations([1,2,3]))) == {[1,3,2],[3,1,2]}
@test collect(filter(x->(iseven(x[1])),combinations([1,2,3],2))) == {[2,3]}
@test collect(partitions(4)) ==  {[4], [3,1], [2,2], [2,1,1], [1,1,1,1]}
@test collect(partitions(8,3)) == {[6,1,1], [5,2,1], [4,3,1], [4,2,2], [3,3,2]}
@test collect(partitions([1,2,3])) == {{[1,2,3]}, {[1,2],[3]}, {[1,3],[2]}, {[1],[2,3]}, {[1],[2],[3]}}
@test collect(partitions([1,2,3,4],3)) == {{[1,2],[3],[4]}, {[1,3],[2],[4]}, {[1],[2,3],[4]}, {[1,4],[2],[3]}, {[1],[2,4],[3]},{[1],[2],[3,4]}}

@test length(collect(partitions(30))) == length(partitions(30))
@test length(collect(partitions(90,4))) == length(partitions(90,4))
@test length(collect(partitions('a':'h'))) == length(partitions('a':'h'))
@test length(collect(partitions('a':'h',5))) == length(partitions('a':'h',5))

for n = 0:7, k = 1:factorial(n)
    p = nthperm!([1:n], k)
    @test isperm(p)
    @test nthperm(p) == k
end

#Issue 6154
@test binomial(int32(34), int32(15)) == binomial(BigInt(34), BigInt(15)) == 1855967520
@test binomial(int64(67), int64(29)) == binomial(BigInt(67), BigInt(29)) == 7886597962249166160
@test binomial(int128(131), int128(62)) == binomial(BigInt(131), BigInt(62)) == 157311720980559117816198361912717812000

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
using LinearAlgebra
using Random
using Test
import KroneckerTools: kron_mul_elem!, kron_mul_elem_t!, a_mul_kron_b!,
    a_mul_b_kron_c!, a_mul_b_kron_ct!, at_mul_b_kron_c!,
    a_mul_b_kron_c_d!, kron_a_mul_b!, kron_at_mul_b!,
    kron_at_kron_b_mul_c!, QuasiTriangular.QuasiUpperTriangular,
    A_mul_B!, At_mul_B!, A_mul_Bt!, At_mul_Bt!

Random.seed!(123)
#a = rand(2,3)
#b = rand(3,4)
#c = zeros(8)

#gemm!('N', 'N', 1.0, vec(a), 2, 3, vec(b), 4, 1.0, c)
#@test reshape(c,2,4) ≈ a*b

#=
for m in [1, 3]
    for n in [3, 4]
        local a = randn(n,n)
        t = schur(a).T
        depth = 4
        for p = 0:4
            for q = depth - p
                d = randn(m*n^(p+q+1))
                d_orig = copy(d)
                w = similar(d)
                println("m = $m, p = $p, q = $q")
                d = copy(d_orig)
                @time kron_mul_elem_t!(w, a, d, n^p, n^q*m)
                @test w ≈ kron(kron(Matrix{Float64}(I(n^p)),a'),Matrix{Float64}(I(m*n^q)))*d_orig
                d = copy(d_orig)
                @time kron_mul_elem_t!(w, t, d, n^p, n^q*m)
                @test w ≈ kron(kron(Matrix{Float64}(I(n^p)),t'),Matrix{Float64}(I(m*n^q)))*d_orig
                d = copy(d_orig)
                @time kron_mul_elem!(w, a, d, n^p, n^q*m)
                @test w ≈ kron(kron(Matrix{Float64}(I(n^p)),a),Matrix{Float64}(I(m*n^q)))*d_orig
                d = copy(d_orig)
                @time kron_mul_elem!(w, t, d, n^p, n^q*m)
                @test w ≈ kron(kron(Matrix{Float64}(I(n^p)),t),Matrix{Float64}(I(m*n^q)))*d_orig
            end
        end
    end
end
=#

ma = 2
na = 4
a = randn(ma, na)
mb = 4
nb = 3
b = randn(mb, nb)
c = randn(ma, nb)
A_mul_B!(c, 1, a, 1, ma, na, b, 1, nb)
@test c ≈ a*b
mb = 2
nb = 3
b = randn(mb, nb)
c = randn(na, nb)
At_mul_B!(c, 1, a, 1, na, ma, b, 1, nb)
@test c ≈ transpose(a)*b
nb = 4
mb = 3
b = randn(mb, nb)
c = randn(ma, mb)
A_mul_Bt!(c, 1, a, 1, ma, na, b, 1, mb)
@test c ≈ a*transpose(b)
mb = 4
nb = 2
b = randn(mb, nb)
c = randn(na, mb)
At_mul_Bt!(c, 1, a, 1, na, ma, b, 1, mb)
@test c ≈ transpose(a)*transpose(b)

order=2
ma = 2
na = 4
q = 2
p = na^(order-1)
a = rand(ma, na)
b_orig = rand(q*na^order)
c = rand(ma*p*q)
b = copy(b_orig)
kron_mul_elem!(c, a, b, p, q)
@test c ≈ kron(kron(I(p),a),I(q))*b_orig

ma = 2
na = 2
a = randn(ma,na)
mc = 3
order = 4
c = randn(mc,mc)
b = randn(na,mc^order)
b_orig = copy(b)
d = randn(ma,mc^order)
w1 = Vector{Float64}(undef, ma*mc^order)
w2 = Vector{Float64}(undef, ma*mc^order)
a_mul_b_kron_c!(d, a, b, c, order, w1, w2)
cc = c
for i = 2:order
    global cc = kron(cc,c)
end
@test d ≈ a*b_orig*cc

b = copy(b_orig)
a_mul_kron_b!(d,b,c,order)
cc = c
for i = 2:order
    global cc = kron(cc,c)
end
@test d ≈ b_orig*cc

order = 3
ma = 2
na = 4
a = randn(ma,na)
mb1 = 2
nb1 = 4
b1 = randn(mb1,nb1)
mb2 = 2
nb2 = 2
b2 = randn(mb2,nb2)
b = [b1, b2]
c = randn(ma,nb1*nb2)
work = zeros(16)
a_mul_kron_b!(c,a,b,work)
@test c ≈ a*kron(b[1],b[2])

order = 3
ma = 2
na = 4
a = randn(ma,na)
mb = 4
nb = 8
b = randn(mb,nb)
c = randn(2,2)
d = randn(2,2)
work1 = zeros(mb*nb)
work2 = zeros(mb*nb)
e = randn(ma,8)
a_mul_b_kron_c_d!(e, a, b, c, d, order, work1, work2)
@test e ≈ a*b*kron(c,kron(d,d))

order = 2
ma = 2
na = 4
a = randn(ma,na)
mb = 4
nb = 8
b = randn(mb,nb)
c = randn(ma*ma*nb)
d = randn(na*na*mb)
work1 = rand(na*na*mb)
work2 = rand(na*na*mb)
@time kron_at_kron_b_mul_c!(d, a, order, b, c, work1, work2)
@time kron_at_kron_b_mul_c!(d, a, order, b, c, work1, work2)

@test d ≈ kron(kron(a',a'),b)*c

println("test1")
order = 2
ma = 2
na = 4
q = 2
a = rand(ma, na)
b = rand(q*na^order)
c = rand(q*ma^order)
work1 = rand(q*max(ma, na)^order)
work2 = similar(work1)
q=2
kron_a_mul_b!(c, a, order, b, q, work1, work2)
@time kron_a_mul_b!(c, a, order, b, q, work1, work2)
@test c ≈ kron(kron(a,a),Matrix{Float64}(I(q)))*b

println("test2")
b = rand(q*ma^order)
c = rand(q*na^order)
kron_at_mul_b!(c, a, order, b, q, work1, work2)
@time kron_at_mul_b!(c, a, order, b, q, work1, work2)
@test c ≈ kron(kron(a',a'),Matrix{Float64}(I(q)))*b

println("test3")
order = 2
mc = 2
nc = 3
c = rand(mc,nc)
ma = 2
na = 4
a = rand(ma,na)
nb = nc^order
b = rand(na,nb)
d = rand(ma,mc^order)
work1 = rand(ma*max(mc, nc)^order)
work2 = similar(work1)
a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
@time a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
@test d ≈ a*b*kron(c',c')

println("test4")
nb = mc^order
b = rand(ma,nb)
d = rand(na,nc^order)
work1 = rand(na*max(mc, nc)^order)
work2 = similar(work1)
at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
@time at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
@test d ≈ a'*b*kron(c,c)

println("test1a")
order = 2
ma = 4
na = 4
q = 2
a = QuasiUpperTriangular(triu(rand(ma, na)))
b = rand(q*na^order)
c = rand(q*ma^order)
work1 = rand(q*max(ma, na)^order)
work2 = similar(work1)
kron_a_mul_b!(c, a, order, b, q, work1, work2)
@time kron_a_mul_b!(c, a, order, b, q, work1, work2)
@test c ≈ kron(kron(a,a),Matrix{Float64}(I(q)))*b

println("test2a")
b = rand(q*ma^order)
c = rand(q*na^order)
kron_at_mul_b!(c, a, order, b, q, work1, work2)
@time kron_at_mul_b!(c, a, order, b, q, work1, work2)
@test c ≈ kron(kron(a',a'),Matrix{Float64}(I(q)))*b

println("test3a")
order = 2
mc = 2
nc = 3
c = rand(mc,nc)
ma = 4
na = 4
a = rand(ma,na)
nb = nc^order
b = rand(na,nb)
d = rand(ma,mc^order)
work1 = rand(ma*max(mc, nc)^order)
work2 = similar(work1)
a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
@time a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
@test d ≈ a*b*kron(c',c')

println("test4a")
nb = mc^order
b = rand(ma,nb)
d = rand(na,nc^order)
work1 = rand(na*max(mc, nc)^order)
work2 = similar(work1)
at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
@time at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
@test d ≈ a'*b*kron(c,c)

println("test1b")
order = 2
ma = 4
na = 4
q = 2
a = QuasiUpperTriangular(triu(rand(ma, na)))
b = view(rand(q*na^order),:)
c = view(rand(q*ma^order),:)
work1 = view(rand(q*max(ma, na)^order),:)
work2 = view(rand(q*max(ma, na)^order),:)
kron_a_mul_b!(c, a, order, b, q, work1, work2)
@time kron_a_mul_b!(c, a, order, b, q, work1, work2)
@test c ≈ kron(kron(a,a),Matrix{Float64}(I(q)))*b

println("test2b")
b = view(rand(q*ma^order),:)
c = view(rand(q*na^order),:)
work1 = view(rand(q*max(ma, na)^order),:)
work2 = view(rand(q*max(ma, na)^order),:)
kron_at_mul_b!(c, a, order, b, q, work1, work2)
@time kron_at_mul_b!(c, a, order, b, q, work1, work2)
@test c ≈ kron(kron(a',a'),Matrix{Float64}(I(q)))*b

println("test3b")
order = 2
mc = 2
nc = 3
c = view(rand(mc,nc),:,:)
ma = 4
na = 4
a = rand(ma,na)
nb = nc^order
b = view(rand(na,nb),:,:)
d = view(rand(ma,mc^order),:,:)
work1 = view(rand(ma*max(mc, nc)^order),:)
work2 = view(rand(ma*max(mc, nc)^order),:)
a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
@time a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
@test d ≈ a*b*kron(c',c')

println("test4b")
nb = mc^order
b = view(rand(ma,nb),:,:)
d = view(rand(na,nc^order),:,:)
work1 = rand(na*max(mc, nc)^order)
work2 = similar(work1)
at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
@time at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
@test d ≈ a'*b*kron(c,c)


text:
# Supertype for all Firm Identifiers
abstract type FirmIdentifier end

# Some identifiers are essentially strings, this is a super class for all of those
abstract type FirmIdentifierString <: FirmIdentifier end

# Reamining are similar to integers, this deals with those cases
abstract type FirmIdentifierInt <: FirmIdentifier end

# Technically, there is only 1 Cusip, however, WRDS considers two types (NCusip and Cusip)
# which are basicaly just old and current versions, since other than the actual
# numbers, these are exactly the samne, this abstract type deals with all of those
# so I can have common functions
abstract type CusipAll <: FirmIdentifierString end

# checksum calculator for Cusip
function luhn_checksum(s)
    if length(s) != 8
        error("Length must be 8 digits")
    end
    tot = 0
    for (i, v) in enumerate(s)
        n = if '0' ≤ v ≤ '9'
            parse(Int, v)
        elseif 'A' ≤ v ≤ 'Z'
            Int(v) - 65 + 10
        elseif 'a' ≤ v ≤ 'z'
            Int(v) - 97 + 10
        elseif v == '*'
            36
        elseif v == '@'
            37
        elseif v == '#'
            38
        end
        if iseven(i)
            n *= 2
        end
        tot += floor(n / 10) + n % 10
    end
    (10 - tot % 10) % 10 |> Int
end

"""
GVKey is the primary identifier in the Compustat universe
It is a string, though only has numeric values
"""
struct GVKey <: FirmIdentifierString
    val::String
    function GVKey(s::String)
        if length(s) > 6
            error("GVKey must be 6 characters long")
        end
        if occursin(r"\D", s)
            error("GVKey can only contain numeric characters")
        end
        new(lpad(s, 6, "0"))
    end
end

GVKey(n::Real) = GVKey(string(Int(n)))

"""
CIK is a common identifier outside of WRDS
It is a string, though has numeric values
"""
struct CIK <: FirmIdentifierString
    val::String
    function CIK(s::String)
        if length(s) > 10
            error("CIK must be 10 characters long")
        end
        if occursin(r"\D", s)
            error("CIK can only contain numeric characters")
        end
        new(lpad(s, 10, "0"))
    end
end

CIK(n::Real) = CIK(string(Int(n)))

"""
Cusip is a common identifier within and outside of WRDS
WRDS tracks older Cusips as NCusip
all Cusips are made up of 3 parts, issuer (6 characters)
issue (2 characters), and a checksum
most databases in WRDS only use the 8 characters
"""
struct Cusip <: CusipAll
    issuer::String
    issue::String
    checksum::Int
    function Cusip(issuer::AbstractString, issue::AbstractString, checksum::Int)
        if length(issuer) != 6
            error("Issuer identification must be 6 characters")
        end
        if length(issue) != 2
            error("Issue must be 2 characters")
        end
        if (!)(0 ≤ checksum ≤ 9)
            error("Checksum must be between 0 and 9 (inclusive)")
        end
        if luhn_checksum(issuer * issue) != checksum
            error("Invalid Cusip or Checksum")
        end
        new(String(issuer), String(issue), checksum)
    end
end

function Cusip(s::AbstractString)
    if length(s) == 9
        Cusip(s[1:6], s[7:8], parse(Int, s[9]))
    elseif length(s) == 8
        Cusip(s[1:6], s[7:8], luhn_checksum(s))
    else
        error("Cusip must be 8 or 9 characters")
    end
end

Cusip(issuer::AbstractString, issue::AbstractString) = Cusip(issuer, issue, luhn_checksum(issuer * issue))

"""
See note on Cusip
"""
struct NCusip <: CusipAll
    issuer::String
    issue::String
    checksum::Int
    function NCusip(issuer::AbstractString, issue::AbstractString, checksum::Int)
        if length(issuer) != 6
            error("Issuer identification must be 6 characters")
        end
        if length(issue) != 2
            error("Issue must be 2 characters")
        end
        if (!)(0 ≤ checksum ≤ 9)
            error("Checksum must be between 0 and 9 (inclusive)")
        end
        if luhn_checksum(issuer * issue) != checksum
            error("Invalid NCusip or Checksum")
        end
        new(string(issuer), string(issue), checksum)
    end
end

function NCusip(s::AbstractString)
    if length(s) == 9
        NCusip(s[1:6], s[7:8], parse(Int, s[9]))
    elseif length(s) == 8
        NCusip(s[1:6], s[7:8], luhn_checksum(s))
    else
        error("NCusip must be 8 or 9 characters")
    end
end

NCusip(issuer::AbstractString, issue::AbstractString) = NCusip(issuer, issue, luhn_checksum(issuer * issue))

Cusip(x::NCusip) = Cusip(x.issuer, x.issue, x.checksum)
NCusip(x::Cusip) = NCusip(x.issuer, x.issue, x.checksum)

"""
Ticker here refers to the ticker that is on the NYSE (not to be confused with IbesTicker)
While a common identifier outside of WRDS, this is not recommended since they change often
"""
struct Ticker <: FirmIdentifierString
    val::String
end

Ticker(s::AbstractString) = Ticker(string(s))

"""
IbesTicker is the primary identifier in the IBES universe
"""
struct IbesTicker <: FirmIdentifierString
    val::String
end

IbesTicker(s::AbstractString) = IbesTicker(string(s))

"""
Permno is the primary identifier of securities in the CRSP universe, it is also
one of the most common methods of linking between databases since it is easy to
find links to IBES, Compustat (GVKey), and Cusip.
"""
struct Permno <: FirmIdentifierInt
    val::Int
end

Permno(x::Real) = Permno(Int(x))

Permno(x::Permno) = x
Cusip(x::Cusip) = x
CIK(x::CIK) = x
GVKey(x::GVKey) = x
NCusip(x::NCusip) = x
Ticker(x::Ticker) = x
IbesTicker(x::IbesTicker) = x

Permno(x::Missing) = x
Cusip(x::Missing) = x
CIK(x::Missing) = x
GVKey(x::Missing) = x
NCusip(x::Missing) = x
Ticker(x::Missing) = x
IbesTicker(x::Missing) = x

######################################################################################
# As long as convert is specified for an identifier, all of the other functions here
# should work automatically work (as long as the type fits into string or integer)
######################################################################################

# Used for Permno
Base.convert(::Type{T}, x::K) where {T<:Real,K<:FirmIdentifierInt} = convert(T, x.val)

# Used for Permno, GVKey, and CIK
Base.convert(::Type{K}, x::Real) where {K<:FirmIdentifier} = K(x)

# Used for GVKey, CIK, Ticker, and IbesTicker
Base.convert(::Type{T}, x::K) where {T<:AbstractString,K<:FirmIdentifierString} = convert(T, x.val)

# Special case for GVKey and CIK
Base.convert(::Type{T}, x::GVKey) where {T<:Real} = parse(T, x.val)
Base.convert(::Type{T}, x::CIK) where {T<:Real} = parse(T, x.val)

# Used for GVKey, CIK, Ticker, IbesTicker, Cusip, and NCusip
Base.convert(::Type{K}, x::AbstractString) where {K<:FirmIdentifierString} = K(x)

# Used for Cusip and NCusip
Base.convert(::Type{T}, x::K) where {T<:AbstractString,K<:CusipAll} = convert(T, "$(x.issuer)$(x.issue)")

# Comparisons
# These all convert the non-identifier into an identifier before back to the base value
# This ensures consistency in the comparison, no issues with CIK being an Int in one place and a string in another

# I am not sure these are actually useful, they were made to try and make joins possible but that is handled by hash
Base.:(==)(id::K, s::T) where {T<:AbstractString, K<:FirmIdentifier} = convert(T, id) == convert(T, convert(K, s))
Base.:(==)(s::T, id::K) where {T<:AbstractString, K<:FirmIdentifier} = convert(T, convert(K, s)) == convert(T, id)

# isless is valuable for sorting
Base.isless(id::K, s::T) where {T<:AbstractString, K<:FirmIdentifier} = isless(convert(T, id), convert(T, convert(K, s)))
Base.isless(s::T, id::K) where {T<:AbstractString, K<:FirmIdentifier} = isless(convert(T, convert(K, s)), convert(T, id))
Base.isless(id1::FirmIdentifier, id2::FirmIdentifier) = isless(value(id1), value(id2))

# String based identifiers

Base.show(io::IOContext, id::FirmIdentifierString) = show(io, convert(String, id))
# Base.print(io::IOContext, x::GVKey) = print(io, x.val)
Base.string(id::FirmIdentifierString) = convert(String, id)

# tryparse is useful for csv
Base.tryparse(::Type{T}, s::AbstractString) where {T<:FirmIdentifierString} = try T(s) catch nothing end
value(id::FirmIdentifierString) = convert(String, id)

# Special "value" parameter for CusipAll so it can return 9 characters
value(x::CusipAll, l::Int=8) = "$(x.issuer)$(x.issue)$(x.checksum)"[1:l]

Base.promote_rule(::Type{T}, ::Type{K}) where {T<:FirmIdentifierString,K<:AbstractString} = K

Base.hash(id::FirmIdentifierString) = hash(convert(String, id))

# Integer based identifiers

Base.show(io::IOContext, id::FirmIdentifierInt) = show(io, convert(Int, id))
Base.string(id::FirmIdentifierInt) = string(convert(Int, id))
value(id::FirmIdentifierInt) = convert(Int, id)

Base.promote_rule(::Type{T}, ::Type{K}) where {T<:FirmIdentifierInt,K<:Real} = K
Base.tryparse(::Type{T}, s::AbstractString) where {T<:FirmIdentifierInt} = try T(parse(Int, s)) catch nothing end
Base.hash(id::FirmIdentifierInt) = hash(convert(Int, id))

# to prevent errors for missing data
value(x::Missing) = missing

######################################################################################
# Future notes: One thing that would be nice is if I could specify `isequal` for
# items that are "number like", for example, it would be nice if
# GVKey("000015") == 15 returned true (which is easy) and the has values were the same
# (which I cannot figure out how to do). As far as I can tell, I can specify the hash
# value as either a String or an Integer, but not both

# Another future task is to make the has value of a tuple to be the same as the base
# version, right now this works well for normal joins but not well for groupby,
# which of course is what my join function relies on
######################################################################################


text:
# ---
# title: 126. Word Ladder II
# id: problem126
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Array, String, Backtracking, Breadth-first Search
# link: <https://leetcode.com/problems/word-ladder-ii/description/>
# hidden: true
# ---
#
# Given two words ( _beginWord_ and _endWord_ ), and a dictionary's word list,
# find all shortest transformation sequence(s) from _beginWord_ to _endWord_ ,
# such that:
#
#   1. Only one letter can be changed at a time
#   2. Each transformed word must exist in the word list. Note that _beginWord_ is _not_ a transformed word.
#
# **Note:**
#
#   * Return an empty list if there is no such transformation sequence.
#   * All words have the same length.
#   * All words contain only lowercase alphabetic characters.
#   * You may assume no duplicates in the word list.
#   * You may assume _beginWord_ and _endWord_ are non-empty and are not the same.
#
# **Example 1:**
#
#
#
#     Input:
#     beginWord = "hit",
#     endWord = "cog",
#     wordList = ["hot","dot","dog","lot","log","cog"]
#
#     Output:
#     [
#       ["hit","hot","dot","dog","cog"],
#       ["hit","hot","lot","log","cog"]
#     ]
#
#
# **Example 2:**
#
#
#
#     Input:
#     beginWord = "hit"
#     endWord = "cog"
#     wordList = ["hot","dot","dog","lot","log"]
#
#     Output: []
#
#     Explanation:  The endWord "cog" is not in wordList, therefore no possible ** ** transformation.
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end