text:
"""
    DecisionRef <: AbstractVariableRef

Holds a reference to the model and the corresponding MOI.VariableIndex.
"""
struct DecisionRef <: JuMP.AbstractVariableRef
    model::JuMP.Model
    index::MOI.VariableIndex
end
"""
    KnownRef <: AbstractVariableRef

Holds a reference to the model and the corresponding MOI.VariableIndex.
"""
struct KnownRef <: JuMP.AbstractVariableRef
    model::JuMP.Model
    index::MOI.VariableIndex
end

# Getters (model) #
# ========================== #
function get_decisions(model::JuMP.Model)
    !haskey(model.ext, :decisions) && return IgnoreDecisions()
    return model.ext[:decisions]
end

function all_decisions(model::JuMP.Model)
    decisions = get_decisions(model)::Decisions
    return all_decisions(decisions)
end

function all_known_decisions(model::JuMP.Model)
    decisions = get_decisions(model)::Decisions
    return all_known_decisions(decisions)
end
"""
    all_decision_variables(model::JuMP.Model)

Returns a list of all decisions currently in the `model`. The decisions are
ordered by creation time.
"""
function all_decision_variables(model::JuMP.Model)
    decisions = get_decisions(model)::Decisions
    return map(decisions.undecided) do index
        DecisionRef(model, index)
    end
end
"""
    all_known_decision_variables(model::JuMP.Model)

Returns a list of all known decisions currently in the `model`. The decisions are
ordered by creation time.
"""
function all_known_decision_variables(model::JuMP.Model)
    decisions = get_decisions(model)::Decisions
    return map(decisions.knowns) do index
        KnownRef(model, index)
    end
end
"""
    num_decisions(model::JuMP.Model)

Return the number of decisions in `model`.
"""
function num_decisions(model::JuMP.Model)
    decisions = get_decisions(model)::Decisions
    return num_decisions(decisions)
end
"""
    num_decisions(model::JuMP.Model)

Return the number of known decisions in `model`.
"""
function num_known_decisions(model::JuMP.Model)
    decisions = get_decisions(model)::Decisions
    return num_known_decisions(decisions)
end

# Getters (refs) #
# ========================== #
function get_decisions(dref::DecisionRef)
    return get_decisions(dref.model)
end

function get_decisions(kref::KnownRef)
    return get_decisions(kref.model)
end
"""
    decision(dref::Union{DecisionRef, KnownRef})

Return the internal `Decision` associated with `dref`.
"""
function decision(dref::Union{DecisionRef, KnownRef})
    decisions = get_decisions(dref)::Decisions
    return decision(decisions, index(dref))
end
"""
    state(dref::DecisionRef)

Return the `DecisionState` of `dref`.
"""
function state(dref::DecisionRef)
    return decision(dref).state
end

# Setters #
# ========================== #
function take_decisions!(model::JuMP.Model, drefs::Vector{DecisionRef}, vals::AbstractVector)
    # Check that all given decisions are in model
    map(dref -> check_belongs_to_model(dref, model), drefs)
    # Check decision length
    length(drefs) == length(vals) || error("Given decision of length $(length(vals)) not compatible with number of decision variables $(length(drefs)).")
    # Update decisions
    for (dref, val) in zip(drefs, vals)
        d = decision(dref)
        # Update state
        d.state = Taken
        # Update value
        d.value = val
    end
    # Update objective and constraints in model
    update_decisions!(model, DecisionsStateChange())
    return nothing
end

function untake_decisions!(model::JuMP.Model, drefs::Vector{DecisionRef})
    # Check that all given decisions are in model
    map(dref -> check_belongs_to_model(dref, model), drefs)
    # Update decisions
    need_update = false
    for dref in drefs
        d = decision(dref)
        if state(d) == Taken
            need_update |= true
            # Update state
            d.state = NotTaken
        end
    end
    # Update objective and constraints in model (if needed)
    need_update && update_decisions!(model, DecisionsStateChange())
    return nothing
end

function update_decisions!(model::JuMP.Model)
    update_decisions!(model, DecisionsStateChange())
end

function update_known_decisions!(model::JuMP.Model)
    decisions = get_decisions(model)
    if decisions isa IgnoreDecisions
        # Nothing to do if decisions are ignored
        # @warn ?
        return nothing
    end
    # Update states
    update_decisions!(model, KnownValuesChange())
    return nothing
end

function update_known_decisions!(model::JuMP.Model, krefs::Vector{KnownRef}, vals::AbstractVector)
    decisions = get_decisions(model)
    if decisions isa IgnoreDecisions
        # Nothing to do if decisions are ignored
        # @warn ?
        return nothing
    end
    # Check that all given decisions are in model
    map(kref -> check_belongs_to_model(kref, model), krefs)
    # Check decision length
    length(krefs) == length(vals) || error("Given decision of length $(length(vals)) not compatible with number of defined known decision variables $(length(krefs)).")
    # Update decisions
    for (kref, val) in zip(krefs, vals)
        d = decision(kref)
        # Update value
        d.value = val
    end
    # Update states
    update_decisions!(model, KnownValuesChange())
    return nothing
end

function update_known_decisions!(model::JuMP.Model, vals::AbstractVector)
    decisions = get_decisions(model)
    if decisions isa IgnoreDecisions
        # Nothing to do if decisions are ignored
        # @warn ?
        return nothing
    end
    # Update values
    update_known_decisions!(decisions, vals)
    # Update states
    update_decisions!(model, KnownValuesChange())
    return nothing
end

# JuMP variable interface #
# ========================== #
function MOI.get(model::JuMP.Model, attr::MOI.AbstractVariableAttribute,
                 dref::Union{DecisionRef, KnownRef})
    check_belongs_to_model(dref, model)
    if MOI.is_set_by_optimize(attr)
        return JuMP._moi_get_result(backend(model), attr, index(dref))
    else
        return MOI.get(backend(model), attr, index(dref))
    end
end

function MOI.set(model::Model, attr::MOI.AbstractVariableAttribute,
                 dref::Union{DecisionRef, KnownRef}, value)
    check_belongs_to_model(dref, model)
    MOI.set(backend(model), attr, index(dref), value)
    return nothing
end

JuMP.name(dref::DecisionRef) = MOI.get(owner_model(dref), MOI.VariableName(), dref)::String
JuMP.name(kref::KnownRef) = MOI.get(owner_model(kref), MOI.VariableName(), kref)::String

function JuMP.set_name(dref::DecisionRef, name::String)
    return MOI.set(owner_model(dref), MOI.VariableName(), dref, name)
end
function JuMP.set_name(kref::KnownRef, name::String)
    return MOI.set(owner_model(dref), MOI.VariableName(), kref, name)
end

function decision_by_name(model::Model, name::String)
    index = MOI.get(backend(model), MOI.VariableIndex, name)
    if index isa Nothing
        return nothing
    else
        return DecisionRef(model, index)
    end
end

function known_decision_by_name(model::Model, name::String)
    index = MOI.get(backend(model), MOI.VariableIndex, name)
    if index isa Nothing
        return nothing
    else
        return KnownRef(model, index)
    end
end

JuMP.index(dref::DecisionRef) = dref.index
JuMP.index(kref::KnownRef) = kref.index

function JuMP.value(dref::DecisionRef; result::Int = 1)::Float64
    if state(dref) == Taken
        # If decision has been fixed the value can be fetched
        # directly
        decision(dref).value
    end
    return MOI.get(owner_model(dref), MOI.VariablePrimal(result), dref)
end

function JuMP.value(kref::KnownRef)::Float64
    decisions = get_decisions(kref)::Decisions
    return decision_value(decisions, index(kref))
end

function JuMP.is_fixed(dref::DecisionRef)
    if state(dref) == Taken
        return true
    end
    return false
end

JuMP.is_fixed(kref::KnownRef) = true

function JuMP.unfix(dref::DecisionRef)
    if state(dref) == NotTaken
        # Nothing to do, just return
        return nothing
    end
    d = decision(dref)
    # Update state
    d.state = NotTaken
    # Prepare modification
    change = DecisionStateChange(index(dref), NotTaken, -value(dref))
    # Update objective and constraints
    update_decisions!(JuMP.owner_model(dref), change)
    return nothing
end

function JuMP.unfix(kref::KnownRef)
    error("Decision with known value cannot be unfixed.")
end
"""
    fix(dref::DecisionRef, val::Number)

Fix the decision associated with `dref` to `val`.

See also [`unfix`](@ref).
"""
function JuMP.fix(dref::DecisionRef, val::Number)
    d = decision(dref)
    if state(dref) == NotTaken
        # Prepare modification
        change = DecisionStateChange(index(dref), Taken, val)
        # Update state
        d.state = Taken
        # Update value
        d.value = val
    else
        # Prepare modification
        change = DecisionStateChange(index(dref), Taken, val - d.value)
        # Just update value
        d.value = val
    end
    # Update objective and constraints
    update_decisions!(JuMP.owner_model(dref), change)
    return nothing
end
"""
    fix(kref::KnownRef, val::Number)

Update the known decision value of `kref` to `val`.
"""
function JuMP.fix(kref::KnownRef, val::Number)
    d = decision(kref)
    # Prepare modification
    change = KnownValueChange(index(kref), val - decision_value(d))
    # Update known value
    d.value = val
    # Update objective and constraints
    update_decisions!(JuMP.owner_model(kref), change)
    return nothing
end

function JuMP.delete(model::JuMP.Model, dref::DecisionRef)
    if model !== owner_model(dref)
        error("The decision you are trying to delete does not " *
              "belong to the model.")
    end
    # TODO What to do here?
end

function JuMP.delete(model::JuMP.Model, kref::KnownRef)
    if model !== owner_model(kref)
        error("The known decision you are trying to delete does not " *
              "belong to the model.")
    end
    # TODO What to do here?
end

JuMP.owner_model(dref::DecisionRef) = dref.model
JuMP.owner_model(kref::KnownRef) = kref.model

struct DecisionNotOwned <: Exception
    dref::DecisionRef
end

struct KnownDecisionNotOwned <: Exception
    kref::KnownRef
end

function JuMP.check_belongs_to_model(dref::DecisionRef, model::AbstractModel)
    if owner_model(dref) !== model
        throw(DecisionNotOwned(dref))
    end
end

function JuMP.check_belongs_to_model(kref::KnownRef, model::AbstractModel)
    if owner_model(kref) !== model
        throw(KnownDecisionNotOwned(kref))
    end
end

function JuMP.is_valid(model::Model, dref::DecisionRef)
    return model === owner_model(dref)
end

function JuMP.is_valid(model::Model, kref::KnownRef)
    return model === owner_model(kref)
end

function JuMP.has_lower_bound(dref::DecisionRef)
    index = MOI.ConstraintIndex{SingleDecision, MOI.GreaterThan{Float64}}(index(dref).value)
    return MOI.is_valid(backend(owner_model(dref)), index)
end
function JuMP.LowerBoundRef(dref::DecisionRef)
    moi_lb =  MOI.ConstraintIndex{SingleDecision, MOI.GreaterThan{Float64}}
    index = moi_lb(index(dref).value)
    return ConstraintRef{JuMP.Model, moi_lb, ScalarShape}(owner_model(dref),
                                                          index,
                                                          ScalarShape())
end
function JuMP.set_lower_bound(dref::Decision, lower::Number)
    new_set = MOI.GreaterThan(convert(Float64, lower))
    if has_lower_bound(dref)
        cindex = MOI.ConstraintIndex{SingleDecision, MOI.GreaterThan{Float64}}(index(dref).value)
        MOI.set(backend(owner_model(dref)), MOI.ConstraintSet(), cindex, new_set)
    else
        MOI.add_constraint(backend(owner_model(dref)), SingleDecision(index(dref)), new_set)
    end
    return nothing
end
function JuMP.delete_lower_bound(dref::DecisionRef)
    JuMP.delete(owner_model(dref), LowerBoundRef(dref))
end
function JuMP.lower_bound(dref::DecisionRef)
    if !has_lower_bound(dref)
        error("Decision $(dref) does not have a lower bound.")
    end
    cset = MOI.get(owner_model(dref), MOI.ConstraintSet(),
                   LowerBoundRef(dref))::MOI.GreaterThan{Float64}
    return cset.lower
end

function JuMP.has_upper_bound(dref::DecisionRef)
    index = MOI.ConstraintIndex{SingleDecision, MOI.LessThan{Float64}}(index(dref).value)
    return MOI.is_valid(backend(owner_model(dref)), index)
end
function JuMP.UpperBoundRef(dref::DecisionRef)
    moi_ub =  MOI.ConstraintIndex{SingleDecision, MOI.LessThan{Float64}}
    index = moi_ub(index(dref).value)
    return ConstraintRef{JuMP.Model, moi_ub, ScalarShape}(owner_model(dref),
                                                          index,
                                                          ScalarShape())
end
function JuMP.set_upper_bound(dref::DecisionRef, lower::Number)
    new_set = MOI.LessThan(convert(Float64, lower))
    if has_upper_bound(dref)
        cindex = MOI.ConstraintIndex{SingleDecision, MOI.LessThan{Float64}}(index(dref).value)
        MOI.set(backend(owner_model(dref)), MOI.ConstraintSet(), cindex, new_set)
    else
        MOI.add_constraint(backend(owner_model(dref)), SingleDecision(index(dref)), new_set)
    end
    return nothing
end
function JuMP.delete_upper_bound(dref::DecisionRef)
    JuMP.delete(owner_model(dref), UpperBoundRef(dref))
end
function JuMP.upper_bound(dref::DecisionRef)
    if !has_upper_bound(dref)
        error("Decision $(dref) does not have a upper bound.")
    end
    cset = MOI.get(owner_model(dref), MOI.ConstraintSet(),
                   UpperBoundRef(dref))::MOI.LessThan{Float64}
    return cset.upper
end

function JuMP.is_integer(dref::DecisionRef)
    index = MOI.ConstraintIndex{SingleDecision, MOI.Integer}(index(dref).value)
    return MOI.is_valid(backend(owner_model(dref)), index)
end
function JuMP.set_integer(dref::DecisionRef)
    if is_integer(dref)
        return nothing
    elseif is_binary(dref)
        error("Cannot set the decision $(dref) to integer as it is already binary.")
    else
        MOI.add_constraint(backend(owner_model(dref)), SingleDecision(index(dref)), MOI.Integer())
    end
end
function unset_integer(dref::DecisionRef)
    JuMP.delete(owner_model(dref), IntegerRef(dref))
    return nothing
end
function JuMP.IntegerRef(dref::DecisionRef)
    moi_int =  MOI.ConstraintIndex{SingleDecision, MOI.Integer}
    index = moi_int(index(dref).value)
    return ConstraintRef{JuMP.Model, moi_int, ScalarShape}(owner_model(dref),
                                                           index,
                                                           ScalarShape())
end

function JuMP.is_binary(dref::DecisionRef)
    index = MOI.ConstraintIndex{SingleDecision, MOI.ZeroOne}(index(dref).value)
    return MOI.is_valid(backend(owner_model(dref)), index)
end
function JuMP.set_binary(dref::DecisionRef)
    if is_binary(dref)
        return nothing
    elseif is_integer(dref)
        error("Cannot set the decision $(dref) to binary as it is already integer.")
    else
        MOI.add_constraint(backend(owner_model(dref)), SingleDecision(index(dref)), MOI.ZeroOne())
    end
end
function unset_binary(dref::DecisionRef)
    JuMP.delete(owner_model(dref), BinaryRef(dref))
    return nothing
end
function JuMP.BinaryRef(dref::DecisionRef)
    moi_bin =  MOI.ConstraintIndex{SingleDecision, MOI.ZeroOne}
    index = moi_bin(index(dref).value)
    return ConstraintRef{JuMP.Model, moi_bin, ScalarShape}(owner_model(dref),
                                                           index,
                                                           ScalarShape())
end

function JuMP.start_value(dref::DecisionRef)
    return MOI.get(owner_model(dref), MOI.VariablePrimalStart(), dref)
end
function JuMP.set_start_value(dref::DecisionRef, value::Number)
    MOI.set(owner_model(dref), MOI.VariablePrimalStart(), dref, Float64(value))
end

JuMP.has_lower_bound(kref::KnownRef) = false
JuMP.LowerBoundRef(kref::KnownRef) = error("Known decision does not have bounds.")
JuMP.set_lower_bound(kref::KnownRef, lower::Number) = error("Known decision does not have bounds.")
JuMP.delete_lower_bound(kref::KnownRef) = error("Known decision does not have bounds.")
JuMP.lower_bound(kref::KnownRef) = error("Known decision does not have bounds.")

JuMP.has_upper_bound(kref::KnownRef) = false
JuMP.UpperBoundRef(kref::KnownRef) = error("Known decision does not have bounds.")
JuMP.set_upper_bound(kref::KnownRef, lower::Number) = error("Known decision does not have bounds.")
JuMP.delete_upper_bound(kref::KnownRef) = error("Known decision does not have bounds.")
JuMP.upper_bound(kref::KnownRef) = error("Known decision does not have bounds.")

JuMP.is_integer(kref::KnownRef) = false
JuMP.set_integer(kref::KnownRef) = error("Known decision does not have integrality constraints.")
JuMP.unset_integer(kref::KnownRef) = error("Known decision does not have integrality constraints.")
JuMP.IntegerRef(kref::KnownRef) = error("Known decision does not have integrality constraints.")

JuMP.is_binary(kref::KnownRef) = false
JuMP.set_binary(kref::KnownRef) = error("Known decision does not have binary constraints.")
JuMP.unset_binary(kref::KnownRef) = error("Known decision does not have binary constraints.")
JuMP.BinaryRef(kref::KnownRef) = error("Known decision does not have binary constraints.")

JuMP.start_value(kref::KnownRef) = error("Known decision does not have start values.")
JuMP.set_start_value(kref::KnownRef, value::Number) = error("Known decision does not have start values.")

function Base.hash(dref::DecisionRef, h::UInt)
    return hash(objectid(owner_model(dref)), hash(dref.index, h))
end
function Base.hash(kref::KnownRef, h::UInt)
    return hash(objectid(owner_model(kref)), hash(kref.index, h))
end

JuMP.isequal_canonical(d::DecisionRef, other::DecisionRef) = isequal(d, other)
function Base.isequal(dref::DecisionRef, other::DecisionRef)
    return owner_model(dref) === owner_model(other) && dref.index == other.index
end
JuMP.isequal_canonical(k::KnownRef, other::KnownRef) = isequal(k, other)
function Base.isequal(kref::KnownRef, other::KnownRef)
    return owner_model(kref) === owner_model(other) && kref.index == other.index
end

Base.iszero(::DecisionRef) = false
Base.copy(dref::DecisionRef) = DecisionRef(dref.model, dref.index)
Base.broadcastable(dref::DecisionRef) = Ref(dref)

Base.iszero(::KnownRef) = false
Base.copy(kref::KnownRef) = KnownRef(dref.model, dref.index)
Base.broadcastable(kref::KnownRef) = Ref(kref)

# JuMP copy interface #
# ========================== #
function JuMP.copy_extension_data(decisions::Decisions, dest::Model, src::Model)
    new_decisions = Decisions()
    for dref in all_decision_variables(src)
        set_decision!(new_decisions, index(dref), decision(dref))
    end
    for kref in all_known_decision_variables(src)
        set_known_decision!(new_decisions, index(kref), decision(kref))
    end
    return new_decisions
end

function Base.getindex(reference_map::JuMP.ReferenceMap, dref::DecisionRef)
    return DecisionRef(reference_map.model,
                       reference_map.index_map[index(dref)])
end

function Base.getindex(reference_map::JuMP.ReferenceMap, kref::KnownRef)
    return DecisionRef(reference_map.model,
                       reference_map.index_map[index(kref)])
end


text:
using LocalMemoize
using Base.Test

@memoize function test_simple(@memo(x), @memo(y), z=8)
    println("First run: $((x, y, z))")
    c = 8x
    @memo a = c + 3
    @memo b = x * 2
    @memo d = a * b * y
    return a + b + d + z
end
@test test_simple(3, 5, 8) == test_simple(3, 5)
@test test_simple(4, 5, 9) != test_simple(3, 5, 9)
cache = get_memos(test_simple)[1]
cache_struct = typeof(cache).parameters[2]
@test fieldnames(cache_struct) == [:a, :b, :d]
@test sort([keys(cache)...]) == [(3, 5), (4, 5)]
@test cache[(3, 5)] == cache_struct(27, 6, 810)
@test cache[(4, 5)] == cache_struct(35, 8, 1400)
@test [fieldtype(cache_struct, i) for i in 1:3] == [Any, Any, Any]
clear_memos(test_simple)
cache = get_memos(test_simple)[1]
@test length(keys(cache)) == 0

@memoize function test_static_type(@memo(x::Float64), y::Float64)
    println("First run: $((x, y))")
    c = x / 2
    @memo a = c + 3
    b = x * 2
    @memo b
    return a * b * y
end
@test test_static_type(1.0, 2.0) == test_static_type(1.0, 2.0)
@test test_static_type(1.0, 2.0) != test_static_type(3.0, 2.0)
cache = get_memos(test_static_type)[1]
cache_struct = typeof(cache).parameters[2]
@test fieldnames(cache_struct) == [:a, :b]
@test sort([keys(cache)...]) == [(1.0,), (3.0,)]
@test cache[(1.0,)] == cache_struct(3.5, 2.0)
@test cache[(3.0,)] == cache_struct(4.5, 6.0)
@test [fieldtype(cache_struct, i) for i in 1:2] == [Float64, Float64]
clear_memos(test_static_type)
cache = get_memos(test_static_type)[1]
@test length(keys(cache)) == 0

@memoize function test_dispatch(@memo(x::Float64), y::Float64)
    @memo a = 2x
    a * y
end
@memoize function test_dispatch(@memo(x::Int64), y::Float64)::Float64
    @memo a = 3x
    a * y
end
@test test_dispatch(1.0, 2.0) == test_dispatch(1.0, 2.0)
@test test_dispatch(Int64(1), 2.0) == test_dispatch(Int64(1), 2.0)

@memoize test_return(x::Int64)::Int64 = x < 3 ? 1 : (test_return(x - 1) + test_return(x - 2))
@test test_return(Int64(5)) == 5
@test test_return(Int64(100)) == 3736710778780434371

ext: jl


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ d0106330-c5a6-11eb-17b7-77a6ed6e1efe
md"""
# Decision Tree

[Reference 1](https://en.wikipedia.org/wiki/Decision_tree)\
[Reference 2](https://scikit-learn.org/stable/modules/tree.html)
"""

# ╔═╡ 88897152-4b3d-4a62-bc67-32970342d24f
md"""
### C4.5 Approach

Algorithm described [here](https://en.wikipedia.org/wiki/C4.5_algorithm)\
My implementation is inspired by the description on [this page](https://cis.temple.edu/~giorgio/cis587/readings/id3-c45.html)

Use __Entropy__ to determine purity:

$$\text{Entropy}(S)=\sum^c_{i=1}-p_i\log_2(p_i)$$

For a given segment of data $S$, $c$ is the number of different class levels, and $p_i$ is the proportion of values falling into class level $i$\
Entropy of 0 means that data is homogenous, and 1 indicates the maximum amount of disorder

Then compute __information gain__ for each feature to decide which to split on:

$$\text{InfoGain}(F)=\text{Entropy}(S_1)-\text{Entropy}(S_2)$$

where $F$ is a feature, $S_1$ is the segment before split, and $S_2$ is the partitions after split\
Entropy on $S_2$ is the weighed sum of entropy on all its partitions ($P_i$), $w_i$ is the proportion of samples falled into $P_i$

$$\text{Entropy}(S_2)=\sum^n_{i=1}w_i\text{Entropy}(P_i)$$

The larger the information gain, the better a feature is splitted to be homogenous

An improvement on information gain is __gain ratios__:

$$\text{GainRatio}=\frac{\text{InfoGain}(F)}{\text{SplitInfo}(F)}$$

where $\text{SplitInfo}$ is the entropy due to the split of $F$

$$\text{SplitInfo}=\sum^c_{j=1}-P_j\log_2(P_j)$$

where each $P_j$ is the proportion of values after partitioning $F$, without considering target classes like entropy function

__Pre-pruning__ (early stopping): As the tree will grow large by dividing features and overfit, we can stop the growing tree after it reaches certain amount of decisions

__Post-pruning__: After the tree grows to its maximum size, cut out nodes and branches that have little effect on accuracy
"""

# ╔═╡ 43e44b5a-62d7-4f95-8685-62470ee55453
module tools include("../tools.jl") end

# ╔═╡ 2a6db899-de2d-497f-8290-f805d59deec6
JuTools = tools.JuTools

# ╔═╡ d5361f35-8091-4f26-83bd-40bd3cd033d3
begin
X_data, Y_data = JuTools.data_generate_cluster_2d(pos1=(30.0, 80.0), pos2=(80.0, 30.0),
    radius1=5.0, radius2=10.0, random_scale=8.0, data_size=1000)
size(X_data), size(Y_data)
end

# ╔═╡ d76345d3-759e-4cde-a371-7c59fdee0b3e
mutable struct DecisionTree
    col_id::Integer
    children::Union{AbstractDict{Function,DecisionTree},Nothing}
    target::Union{AbstractFloat,Nothing}
end

# ╔═╡ a5c720d7-aacc-4e91-a045-37aecfcef305
function create_decision_tree(X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number; max_depth::Integer=10)::DecisionTree
    @assert ndims(X_data) == ndims(Y_data) + 1 == 2
    @assert size(X_data)[1] == size(Y_data)[1]
    @assert max_depth >= 1

    function majority(Y_vec::Array)::Number
        unique_vals = Dict{Number, Integer}()
        for Y_val in Y_vec
            if !haskey(unique_vals, Y_val)
                unique_vals[Y_val] = 1
            else
                unique_vals[Y_val] += 1
            end
        end
        result = sort(collect(unique_vals), by=m->m[2])
        return result[end][1]
    end

    function is_same_class(Y_vec::Array)::Bool
        val = Y_vec[1]
        for m in Y_vec[2:end]
            if m != val
                return false
            end
        end
        return true
    end

    function get_unique_array(X_vec::Array)::Array
        X_max = maximum(X_vec)
        X_min = minimum(X_vec)
        round_digits = -Integer(trunc(log10(X_max-X_min)-2))
        result = []
        for m in X_vec
            if !(trunc(m, digits=round_digits) in result)
                push!(result, m)
            end
        end
        return result
    end

    function info(Y_data::Array)::AbstractFloat
        portion = Dict{Number,Integer}()
        for m in Y_data
            if !haskey(portion, m)
                portion[m] = 1
            else
                portion[m] += 1
            end
        end
        T = size(Y_data)[1]
        result = 0.0
        for m in keys(portion)
            p = portion[m] / T
            result += -p * log2(p)
        end
        return result
    end

    function create_decision_tree_recursive(X_data::Array, Y_data::Array, X_index_visited::Array, max_depth::Integer)::DecisionTree
        if is_same_class(Y_data)
            return DecisionTree(-1, nothing, Y_data[1])
        end
        if max_depth <= 0
            return DecisionTree(-1, nothing, majority(Y_data))
        end
        max_index = -1
        max_ratio = 0.0
        max_fn = nothing
        for i in 1:size(X_data)[2]
            if i in X_index_visited
                continue
            end
            X_vec_unique = get_unique_array(X_data[:, i])
            if length(X_vec_unique) > 10
                # continous data, then choose binary threshold
                fn_max = nothing
                ratio_max = 0
                for j in 1:(length(X_vec_unique)-1)
                    # choose average as threshold
                    threshold = (X_vec_unique[j] + X_vec_unique[j+1]) / 2.0
                    fn = (m) -> m <= threshold
                    # compute information gain
                    Ii = info(Y_data)
                    Y_1 = Y_data[fn.(X_data[:, i])]
                    Y_2 = Y_data[(!fn).(X_data[:, i])]
                    I1 = info(Y_1)
                    I2 = info(Y_2)
                    gain = Ii - (length(Y_1)/length(Y_data)*I1 + length(Y_2)/length(Y_data)*I2)
                    splitInfo = -(length(Y_1)/length(Y_data))*log2(length(Y_1)/length(Y_data)) - (length(Y_2)/length(Y_data))*log2(length(Y_2)/length(Y_data))
                    ratio = gain / splitInfo
                    if fn_max === nothing
                        fn_max = fn
                        ratio_max = ratio
                    else
                        if ratio > ratio_max
                            fn_max = fn
                            ratio_max = ratio
                        end
                    end
                end
                if ratio_max > max_ratio
                    max_index = i
                    max_fn = fn_max
                    max_ratio = ratio_max
                end
            else
                # categorical data, then compute by splitting on each category
                Ii = info(Y_data)
                gain = Ii
                splitInfo = 0.0
                for val in X_vec_unique
                    fn = (m) -> m == val
                    Y0 = Y_data[fn.(X_data[:, i])]
                    I0 = info(Y0)
                    gain -= (length(Y0)/length(Y_data))*I0
                    splitInfo -= (length(Y0)/length(Y_data))*log2(length(Y0)/length(Y_data))
                end
                ratio = gain / splitInfo
                if ratio > max_ratio
                    max_index = i
                    max_fn = nothing
                    max_ratio = ratio
                end
            end
        end
        if max_index < 0
            return DecisionTree(-1, nothing, majority(Y_data))
        end
        X_index_visited = copy(X_index_visited)
        push!(X_index_visited, max_index)
        children = Dict{Function,DecisionTree}()
        if max_fn === nothing
            # categorical
            X_vec_unique = get_unique_array(X_data[:, max_index])
            for val in X_vec_unique
                fn = (m) -> m == val
                identify = fn.(X_data[:, max_index])
                X_part = X_data[identify, :]
                Y_part = Y_data[identity]
                children[fn] = create_decision_tree_recursive(X_part, Y_part, X_index_visited, max_depth-1)
            end
        else
            # continuous
            identify_1 = max_fn.(X_data[:, max_index])
            X_part_1 = X_data[identify_1, :]
            Y_part_1 = Y_data[identify_1]
            children[max_fn] = create_decision_tree_recursive(X_part_1, Y_part_1, X_index_visited, max_depth-1)
            identify_2 = (!max_fn).(X_data[:, max_index])
            X_part_2 = X_data[identify_2, :]
            Y_part_2 = Y_data[identify_2]
            children[(!max_fn)] = create_decision_tree_recursive(X_part_2, Y_part_2, X_index_visited, max_depth-1)
        end
        return DecisionTree(max_index, children, nothing)
    end

    return create_decision_tree_recursive(X_data, Y_data, [], max_depth)
end

# ╔═╡ fe83da85-c0fd-4955-98f5-270a8e6caa77
function predict(X_data::Array{T} where T<:Number, tree::DecisionTree)::Array
    if ndims(X_data) == 1
        X_data = reshape(X_data, (1, size(X_data)[1]))
    end
    @assert ndims(X_data) == 2
    prediction = []
    for i in 1:size(X_data)[1]
        X_vec = X_data[i, :]
        Y_pred = nothing
        tree_copy = tree
        while true
            if tree_copy.target !== nothing
                Y_pred = tree_copy.target
                break
            end
            @assert tree_copy.children !== nothing
            updated = false
            for key_fn in keys(tree_copy.children)
                if key_fn(X_vec[tree_copy.col_id])
                    tree_copy = tree_copy.children[key_fn]
                    updated = true
                    break
                end
            end
            if !updated
                println("Error occured for $X_vec")
                break
            end
        end
        push!(prediction, Y_pred)
    end
    return prediction
end

# ╔═╡ 28c70ec7-ebc5-4826-a0d0-b80cc8603dc7
tree = create_decision_tree(X_data, Y_data, max_depth=20)

# ╔═╡ 2301126e-6659-4162-bb5f-4eadc0601f88
JuTools.compute_accuracy(predict(X_data, tree), Y_data)

# ╔═╡ Cell order:
# ╟─d0106330-c5a6-11eb-17b7-77a6ed6e1efe
# ╟─88897152-4b3d-4a62-bc67-32970342d24f
# ╠═43e44b5a-62d7-4f95-8685-62470ee55453
# ╠═2a6db899-de2d-497f-8290-f805d59deec6
# ╠═d5361f35-8091-4f26-83bd-40bd3cd033d3
# ╠═d76345d3-759e-4cde-a371-7c59fdee0b3e
# ╠═a5c720d7-aacc-4e91-a045-37aecfcef305
# ╠═fe83da85-c0fd-4955-98f5-270a8e6caa77
# ╠═28c70ec7-ebc5-4826-a0d0-b80cc8603dc7
# ╠═2301126e-6659-4162-bb5f-4eadc0601f88

max_stars_repo_name: teamclouday/JuliaLearn

max_issues_repo_name: teamclouday/JuliaLearn

max_forks_repo_name: teamclouday/JuliaLearn


text:
using GLMakie: Theme, Axis3, RGBAf

flo_theme = Theme(
    Axis3 = (
        backgroundcolor = :black,
        protrusions = (0,0,0,0),
        titlealign = :left,
        viewmode = :fit,
        perspectiveness = 0.25
    )
)

set_theme!(flo_theme)


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: jl:light,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.1
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# +
"""Tutorial: CEPA0 and CCD"""

__author__    = ["D. Menendez", "Adam S. Abbott"]
__credit__    = ["D. Menendez", "Adam S. Abbott", "Justin M. Turney"]

__copyright__ = "(c) 2014-2020, The Psi4Julia Developers"
__license__   = "BSD-3-Clause"
__date__      = "2020-08-02"
# -

# # Introduction
# In this tutorial, we will implement the coupled-electron pair approximation (CEPA0) and coupled-cluster doubles (CCD) methods using our spin orbital framework covered in the [previous tutorial](8a_Intro_to_spin_orbital_postHF.ipynb).
#
#
# ### I. Coupled Cluster Theory
#
# In single reference coupled cluster theory, dynamic correlation is acquired by operating an exponential operator on some reference determinant, such as a Hartree-Fock wavefunction, to obtain the coupled cluster wavefunction given by:
#
# \begin{equation}
# \mid \mathrm{\Psi_{CC}} \rangle = \exp(\hat{T}) \mid \mathrm{\Phi} \rangle
# \end{equation}
#
# where $\hat{T} = T_1 + T_2 + ... + T_n$ is the sum of "cluster operators" which act on our reference wavefunction to excite electrons from occupied ($i, j, k$...) to virtual ($a, b, c$...) orbitals. In second quantization, these cluster operators are expressed as:
#
# \begin{equation}
# T_k = \left(\frac{1}{k!}\right)^2 \sum_{\substack{i_1 \ldots i_k \\ a_1 \ldots a_k }} t_{i_1 \ldots i_k}^{a_1 \ldots a_k}   a_{a_1}^{\dagger} \ldots a_{a_k}^{\dagger} a_{i_k} \ldots a_{i_1}
# \end{equation}
#
# where $t$ is the $t$-amplitude, and $a^{\dagger}$ and $a$ are creation and annihilation operators.
#
# ### II. Coupled Cluster Doubles
# For CCD, we only include the doubles cluster operator:
#
# \begin{equation}
# \mid \mathrm{\Psi_{CCD}} \rangle = \exp(T_2) \mid \mathrm{\Phi} \rangle
# \end{equation}
#
# The CCD Schrödinger equation is
#
# \begin{equation}
# \hat{H} \mid \mathrm{\Psi_{CCD}} \rangle = E \mid \mathrm{\Psi_{CCD}}\rangle
# \end{equation}
#
# The details will not be covered here, but if we project the CCD Schrödinger equation on the left by our Hartree-Fock reference determinant $ \langle \mathrm{\Phi}\mid $, assuming intermediate normalization $\langle \Phi \mid \mathrm{\Psi_{CCD}} \rangle = 1$, we obtain:
#
# \begin{equation}
#  \langle \Phi \mid \hat{H} \space \exp(T_2) \mid \Phi \rangle = E
# \end{equation}
#
# which is most easily evaluated with a diagrammatic application of Wick's theorem. Assuming Brillouin's theorem applies (that is, our reference is a Hartree-Fock wavefunction) we obtain:
#
# \begin{equation}
# E_{\mathrm{CCD}} = \tfrac{1}{4} \bar{g}_{ij}^{ab} t_{ab}^{ij}
# \end{equation}
#
# A somewhat more involved derivation is that of the $t$-amplitudes. These are obtained in a similar fashion to the energy expression, this time projecting the CCD Schrödinger equation on the left by a doubly-excited reference determinant $ \langle\Phi_{ij}^{ab}\mid $:
#
# \begin{equation}
# \langle\Phi_{ij}^{ab}\mid \hat{H} \space \exp(T_2) \mid \Phi \rangle
# \end{equation}
#
# I will spare you the details of solving this expectation value as well. But, if one evaluates the diagrams via Wick's theorem and simplifies, the $t$-amplitudes are given by:
#
# \begin{equation}
# t_{ab}^{ij} = (\mathcal{E}_{ab}^{ij})^{-1} \left( \bar{g}_{ab}^{ij} + \tfrac{1}{2} \bar{g}_{ab}^{cd} t_{cd}^{ij} + \tfrac{1}{2} \bar{g}_{kl}^{ij} t_{ab}^{kl}  + \hat{P}_{(a \space / \space b)}^{(i \space / \space j)} \bar{g}_{ak}^{ic} t_{bc}^{jk} - \tfrac{1}{2}\hat{P}_{(a \space / \space b)} \bar{g}_{kl}^{cd} t_{ac}^{ij} t_{bd}^{kl} - \tfrac{1}{2}  \hat{P}^{(i \space / \space j)} \bar{g}_{kl}^{cd} t_{ab}^{ik} t_{cd}^{jl} + \tfrac{1}{4} \bar{g}_{kl}^{cd} t_{cd}^{ij} t_{ab}^{kl} +  \hat{P}^{(i \space / \space j)} \bar{g}_{kl}^{cd} t_{ac}^{ik} t_{bd}^{jl} \right)
# \end{equation}
#
# where $(\mathcal{E}_{ab}^{ij})^{-1}$ is the orbital energy denominator, more familiarly known as
#
# \begin{equation}
# (\mathcal{E}_{ab}^{ij})^{-1} = \frac{1}{\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b}
# \end{equation}
#
# and $\bar{g}_{pq}^{rs}$ is the antisymmetrized two-electron integral in physicist's notation $\langle pq \mid\mid rs \rangle$. $\hat{P}$ is the *antisymmetric permutation operator*. This operator acts on a term to produce the sum of the permutations of the indicated indices, with an appropriate sign factor. Its effect is best illustrated by an example. Consider the fourth term, which is really four terms in one.
#
# $\hat{P}_{(a \space / \space b)}^{(i \space / \space j)} \bar{g}_{ak}^{ic} t_{bc}^{jk}$ produces:
#
# 1. The original: $ \quad \bar{g}_{ak}^{ic} t_{bc}^{jk} \\ $
#
# 2. Permuation of $a$ and $b$: $ \quad  \textrm{-} \bar{g}_{bk}^{ic} t_{ac}^{jk} \\ $
#
# 3. Permuation of $i$ and $j$: $ \quad \, \, \textrm{-} \bar{g}_{ak}^{jc} t_{bc}^{ik} \\ $
#
# 4. Permuation of $a$ and $b$, $i$ and $j$: $ \quad \bar{g}_{bk}^{jc} t_{ac}^{ik} \\ $
#
#
# Note that each permutation adds a sign change. This shorthand notation keeps the equation in a more manageable form.
#
# Since the $t$-amplitudes and the energy depend on $t$-amplitudes, we must iteratively solve these equations until they reach self consistency, and the energy converges to some threshold.
#
# ### III. Retrieving MP2 and CEPA0 from the CCD equations
# It is interesting to note that if we only consider the first term of the expression for the doubles amplitude $t_{ab}^{ij}$ and plug it into the energy expression, we obtain the MP2 energy expression:
#
# \begin{equation}
# t_{ab}^{ij} = (\mathcal{E}_{ab}^{ij})^{-1} \bar{g}_{ab}^{ij}
# \end{equation}
#
# \begin{equation}
# E_{\mathrm{MP2}} = \tfrac{1}{4}  \bar{g}_{ij}^{ab} t_{ab}^{ij}  = \tfrac{1}{4} \bar{g}_{ij}^{ab} \bar{g}_{ab}^{ij}   (\mathcal{E}_{ab}^{ij})^{-1}
# \end{equation}
#
# Furthermore, if we leave out the quadratic terms in the CCD amplitude equation (terms containing two $t$-amplitudes), we obtain the coupled electron-pair approximation (CEPA0):
# \begin{equation}
# t_{ab}^{ij} = (\mathcal{E}_{ab}^{ij})^{-1} \left( \bar{g}_{ab}^{ij} + \tfrac{1}{2} \bar{g}_{ab}^{cd} t_{cd}^{ij} + \tfrac{1}{2} \bar{g}_{kl}^{ij} t_{ab}^{kl}  + \hat{P}_{(a \space / \space b)}^{(i \space / \space j)} \bar{g}_{ak}^{ic} t_{bc}^{jk} \right)
# \end{equation}
#
# The CEPA0 energy expression is identical:
#
# \begin{equation}
# E_{\mathrm{CEPA0}} = \tfrac{1}{4} \bar{g}_{ij}^{ab} t_{ab}^{ij}
# \end{equation}
#
# Using our spin orbital setup for the MO coefficients, orbital energies, and two-electron integrals used in the [previous tutorial](8a_Intro_to_spin_orbital_postHF.ipynb), we are equipped to program the expressions for the CEPA0 and CCD correlation energy.

# ### Implementation: CEPA0 and CCD
# As usual, we import Psi4, NumPy, and TensorOperations, and set the appropriate options.

# +
# ==> Import statements & Global Options <==
using PyCall: pyimport
psi4 = pyimport("psi4")
np   = pyimport("numpy")
using TensorOperations: @tensor
using Formatting: printfmt

psi4.set_memory(Int(2e9))
numpy_memory = 2
psi4.core.set_output_file("output.dat", false)

# +
# ==> Molecule & Psi4 Options Definitions <==
mol = psi4.geometry("""
0 1
O
H 1 1.1
H 1 1.1 2 104
symmetry c1
""")

psi4.set_options(Dict("basis"         => "6-31g",
                      "scf_type"      => "pk",
                      "reference"     => "rhf",
                      "mp2_type"      => "conv",
                      "e_convergence" => 1e-8,
                      "d_convergence" => 1e-8))
# -

# Note that since we are using a spin orbital setup, we are free to use any Hartree-Fock reference we want. Here we choose RHF. For convenience, we let Psi4 take care of the Hartree-Fock procedure, and return the wavefunction object.

# Get the SCF wavefunction & energies
scf_e, scf_wfn = psi4.energy("scf", return_wfn=true)

# Load in information about the basis set and orbitals using MintsHelper and the wavefunction:

mints = psi4.core.MintsHelper(scf_wfn.basisset())
nbf = mints.nbf()           # number of basis functions
nso = 2nbf                  # number of spin orbitals
nalpha = scf_wfn.nalpha()   # number of alpha electrons
nbeta = scf_wfn.nbeta()     # number of beta electrons
nocc = nalpha + nbeta       # number of occupied orbitals
nvirt = 2nbf - nocc         # number of virtual orbitals

# Spin-block our MO coefficients and two-electron integrals, just like in the spin orbital MP2 code:

# +
Ca = np.asarray(scf_wfn.Ca())
Cb = np.asarray(scf_wfn.Cb())
C = [Ca zero(Ca); zero(Cb) Cb]; # direct sum

# Result: | Ca  0 |
#         | 0   Cb|

# +
# Get the two electron integrals using MintsHelper
I = np.asarray(mints.ao_eri())

"""
Function that spin blocks two-electron integrals
Using `np.kron`, we project I into the space of the 2x2 identity, tranpose the result
and project into the space of the 2x2 identity again. This doubles the size of each axis.
The result is our two electron integral tensor in the spin orbital form.
"""
function spin_block_tei(I)
    identity = [ 1.0 0.0; 0.0 1.0]
    I = np.kron(identity, I)
    np.kron(identity, permutedims(I, reverse(1:4)))
end

# Spin-block the two electron integral array
I_spinblock = spin_block_tei(I);
# -

# Convert two-electron integrals to antisymmetrized physicist's notation:

# Converts chemist's notation to physicist's notation, and antisymmetrize
# (pq|rs) ↦ ⟨pr|qs⟩
# Physicist's notation
tmp = permutedims(I_spinblock, (1, 3, 2, 4))
# Antisymmetrize:
# ⟨pr||qs⟩ = ⟨pr|qs⟩ - ⟨pr|sq⟩
gao = tmp - permutedims(tmp, (1, 2, 4, 3));

# Obtain the orbital energies, append them, and sort the columns of our MO coefficient matrix according to the increasing order of orbital energies.

# +
# Get orbital energies
eps_a = np.asarray(scf_wfn.epsilon_a())
eps_b = np.asarray(scf_wfn.epsilon_b())
eps = vcat(eps_a, eps_b)

# Before sorting the orbital energies, we can use their current arrangement to sort the columns
# of C. Currently, each element i of eps corresponds to the column i of C, but we want both
# eps and columns of C to be in increasing order of orbital energies

# Sort the columns of C according to the order of increasing orbital energies
C = C[:, sortperm(eps)]

# Sort orbital energies in increasing order
sort!(eps);
# -

# Finally, we transform our two-electron integrals to the MO basis. Here, we denote the integrals as `gmo` to differentiate from the chemist's notation integrals `I_mo`.

# Transform gao, which is the spin-blocked 4d array of physicist's notation,
# antisymmetric two-electron integrals, into the MO basis using MO coefficients
gmo = @tensor begin
   gmo[P,Q,R,S] := gao[p,Q,R,S] * C[p,P]
   gmo[p,Q,R,S] := gmo[p,q,R,S] * C[q,Q]
   gmo[p,q,R,S] := gmo[p,q,r,S] * C[r,R]
   gmo[p,q,r,S] := gmo[p,q,r,s] * C[s,S]
end
nothing

# Construct the 4-dimensional array of orbital energy denominators:

# Define slices, create 4 dimensional orbital energy denominator tensor
n = [CartesianIndex()]
o = [p ≤ nocc for p in 1:nso]
v = [p > nocc for p in 1:nso]
e_denom = @. inv(-eps[v, n, n, n] - eps[n, v, n, n] + eps[n, n, o, n] + eps[n, n, n, o]);

# We now have everything we need to construct our $t$-amplitudes and iteratively solve for our CEPA0 and CCD energy. To build the $t$-amplitudes, we first construct an empty 4-dimensional array to store them.

# Create space to store t amplitudes
t_amp = zeros(nvirt, nvirt, nocc, nocc);

# # Implementation: CEPA0
# First we will program CEPA0. Recall the expression for the $t$-amplitudes:
#
# \begin{equation}
# t_{ab}^{ij} = (\mathcal{E}_{ab}^{ij})^{-1} \left( \bar{g}_{ab}^{ij} + \tfrac{1}{2} \bar{g}_{ab}^{cd} t_{cd}^{ij} + \tfrac{1}{2} \bar{g}_{kl}^{ij} t_{ab}^{kl}  + \hat{P}_{(a \space / \space b)}^{(i \space / \space j)} \bar{g}_{ak}^{ic} t_{bc}^{jk} \right)
# \end{equation}
#
# These terms translate naturally into code using Julia's `@tensor` function. To access only the occupied and virtual indices of `gmo` we use our slices defined above. The permutation operator terms can be easily obtained by transposing the original result accordingly. To construct each iteration's $t$-amplitude:
#
# ~~~julia
# mp2    = @view gmo[v, v, o, o]
# @tensor cepa1[ a,b,i,j] := 0.5(gmo[v,v,v,v])[a,b,c,d] * t_amp[c,d,i,j]
# @tensor cepa2[ a,b,i,j] := 0.5(gmo[o,o,o,o])[k,l,i,j] * t_amp[a,b,k,l]
# @tensor cepa3a[a,b,i,j] :=    (gmo[v,o,o,v])[a,k,i,c] * t_amp[b,c,j,k]
# cepa3b = -permutedims(cepa3a, (2, 1, 3, 4)) # a <-> b
# cepa3c = -permutedims(cepa3a, (1, 2, 4, 3)) #          i <-> j
# cepa3d =  permutedims(cepa3a, (2, 1, 4, 3)) # a <-> b, i <-> j
# cepa3  =  cepa3a + cepa3b + cepa3c + cepa3d
#
# t_amp_new = @. e_denom * (mp2 + cepa1 + cepa2 + cepa3)
# ~~~
#
# To evaluate the energy, $E_{\mathrm{CEPA0}} = \tfrac{1}{4} \bar{g}_{ij}^{ab} t_{ab}^{ij}$,
#
# ~~~julia
# E_CEPA0 = 1/4 * @tensor scalar((gmo[o,o,v,v])[i,j,a,b] * t_amp_new[a,b,i,j])
# ~~~

# Putting it all together, we initialize the energy, set the max iterations, and iterate the energy until it converges to our convergence criterion:

# +
# Initialize energy
E_CEPA0 = let E_CEPA0 = 0.0, gmo=gmo, o=o,v=v, e_denom = e_denom, t_amp = t_amp

   MAXITER = 50

   for cc_iter in 1:MAXITER
       E_old = E_CEPA0

       # Collect terms
       mp2      = @view gmo[v,v,o,o]
       @tensor cepa1[ a,b,i,j] := 0.5(gmo[v,v,v,v])[a,b,c,d] * t_amp[c,d,i,j]
       @tensor cepa2[ a,b,i,j] := 0.5(gmo[o,o,o,o])[k,l,i,j] * t_amp[a,b,k,l]
       @tensor cepa3a[a,b,i,j] :=    (gmo[v,o,o,v])[a,k,i,c] * t_amp[b,c,j,k]
       cepa3b = -permutedims(cepa3a, (2, 1, 3, 4))
       cepa3c = -permutedims(cepa3a, (1, 2, 4, 3))
       cepa3d =  permutedims(cepa3a, (2, 1, 4, 3))
       cepa3  =  cepa3a + cepa3b + cepa3c + cepa3d

       # Update t amplitude
       t_amp_new = @. e_denom * (mp2 + cepa1 + cepa2 + cepa3)

       # Evaluate Energy
       E_CEPA0 = 1/4 * @tensor scalar((gmo[o,o,v,v])[i,j,a,b] * t_amp_new[a,b,i,j])
       t_amp = t_amp_new
       dE = E_CEPA0 - E_old
       printfmt("CEPA0 Iteration {1:3d}: Energy = {2:4.12f} dE = {3:1.5e}\n", cc_iter, E_CEPA0, dE)

       if abs(dE) < 1.e-8
           @info "CEPA0 Iterations have converged!"
           break
       end

       if cc_iter == MAXITER
           psi4.core.clean()
           error("Maximum number of iterations exceeded.")
       end
   end
   E_CEPA0
end

printfmt("\nCEPA0 Correlation Energy: {:5.15f}\n", E_CEPA0)
printfmt("CEPA0 Total Energy: {:5.15f}\n", E_CEPA0 + scf_e)
# -

# Since `t_amp` is initialized to zero, the very first iteration should be the MP2 correlation energy. We can check the final CEPA0 energy with Psi4. The method is called `lccd`, or linear CCD, since CEPA0 omits the terms with two cluster amplitudes.

psi4.compare_values(psi4.energy("lccd"), E_CEPA0 + scf_e, 6, "CEPA0 Energy")

# # Implementation: CCD
#
# To code CCD, we only have to add in the last four terms in our expression for the $t$-amplitudes:
#
# \begin{equation}
# t_{ab}^{ij} = (\mathcal{E}_{ab}^{ij})^{-1} \left( \bar{g}_{ab}^{ij} + \tfrac{1}{2} \bar{g}_{ab}^{cd} t_{cd}^{ij} + \tfrac{1}{2} \bar{g}_{kl}^{ij} t_{ab}^{kl}  + \hat{P}_{(a \space / \space b)}^{(i \space / \space j)} \bar{g}_{ak}^{ic} t_{bc}^{jk} - \underline{\tfrac{1}{2}\hat{P}_{(a \space / \space b)} \bar{g}_{kl}^{cd} t_{ac}^{ij} t_{bd}^{kl} - \tfrac{1}{2}  \hat{P}^{(i \space / \space j)} \bar{g}_{kl}^{cd} t_{ab}^{ik} t_{cd}^{jl} + \tfrac{1}{4} \bar{g}_{kl}^{cd} t_{cd}^{ij} t_{ab}^{kl} +  \hat{P}^{(i \space / \space j)} \bar{g}_{kl}^{cd} t_{ac}^{ik} t_{bd}^{jl}} \right)
# \end{equation}
#
# which we readily translate into `@tensor`'s:
#
# ~~~julia
# @tensor ccd1a[a,b,i,j] := (gmo[o,o,v,v])[k,l,c,d] * t_amp[a,c,i,j] * t_amp[b,d,k,l]
# ccd1b  = -permutedims(ccd1a, (2, 1, 3, 4))
# ccd1   = -0.5(ccd1a + ccd1b)
#
# @tensor ccd2a[a,b,i,j] := (gmo[o,o,v,v])[k,l,c,d] * t_amp[a,b,i,k] * t_amp[c,d,j,l]
# ccd2b  = -permutedims(ccd2a, (1, 2, 4, 3))
# ccd2   = -0.5(ccd2a + ccd2b)
#
# @tensor ccd3[a,b,i,j] := 1/4 * (gmo[o,o,v,v])[k,l,c,d] * t_amp[c,d,i,j] * t_amp[a,b,k,l]
#
# @tensor ccd4a[a,b,i,j] := (gmo[o,o,v,v])[ k,l,c,d] * t_amp[a,c,i,k] * t_amp[b,d,j,l]
# ccd4b  = -permutedims(ccd4a, (1, 2, 4, 3))
# ccd4   = (ccd4a + ccd4b)
# ~~~
#
# and the energy expression is identical to CEPA0:
# \begin{equation}
# E_{CCD } = \tfrac{1}{4} \bar{g}_{ij}^{ab} t_{ab}^{ij}
# \end{equation}
#
# Adding the above terms to our CEPA0 code will compute the CCD correlation energy (may take a minute or two to run):

# +
# Initialize energy
E_CCD = let E_CCD = 0.0, o=o,v=v, e_denom=e_denom, t_amp=t_amp

   MAXITER = 50

   # Create space to store t amplitudes
   t_amp = zeros(nvirt, nvirt, nocc, nocc)
   for cc_iter in 1:MAXITER
       E_old = E_CCD

       # Collect terms
       mp2      = @view gmo[v,v,o,o]
       @tensor cepa1[ a,b,i,j] := 0.5(gmo[v,v,v,v])[a,b,c,d] * t_amp[c,d,i,j]
       @tensor cepa2[ a,b,i,j] := 0.5(gmo[o,o,o,o])[k,l,i,j] * t_amp[a,b,k,l]
       @tensor cepa3a[a,b,i,j] :=    (gmo[v,o,o,v])[a,k,i,c] * t_amp[b,c,j,k]
       cepa3b = -permutedims(cepa3a, (2, 1, 3, 4))
       cepa3c = -permutedims(cepa3a, (1, 2, 4, 3))
       cepa3d =  permutedims(cepa3a, (2, 1, 4, 3))
       cepa3  =  cepa3a + cepa3b + cepa3c + cepa3d

       @tensor ccd1a_ref[a,b,i,j] := (gmo[o,o,v,v])[k,l,c,d] * t_amp[a,c,i,j] * t_amp[b,d,k,l]
       @tensor ccd1a_tmp[c,b]     := (gmo[o,o,v,v])[k,l,c,d] * t_amp[b,d,k,l]
       @tensor ccd1a[a,b,i,j]     := ccd1a_tmp[c,b]    * t_amp[a,c,i,j]
       println(isapprox(ccd1a_ref, ccd1a))

       ccd1b  = -permutedims(ccd1a, (2, 1, 3, 4))
       ccd1   = -0.5(ccd1a + ccd1b)

       @tensor ccd2a_ref[a,b,i,j] := (gmo[o,o,v,v])[k,l,c,d] * t_amp[a,b,i,k] * t_amp[c,d,j,l]
       @tensor ccd2a_tmp[j,k]     := (gmo[o,o,v,v])[k,l,c,d] * t_amp[c,d,j,l]
       @tensor ccd2a[a,b,i,j]     := ccd2a_tmp[j,k]    * t_amp[a,b,i,k]
       println(isapprox(ccd2a_ref, ccd2a))

       ccd2b  = -permutedims(ccd2a, (1, 2, 4, 3))
       ccd2   = -0.5(ccd2a + ccd2b)

       @tensor ccd3_ref[a,b,i,j] := 1/4 * (gmo[o,o,v,v])[k,l,c,d] * t_amp[c,d,i,j] * t_amp[a,b,k,l]
       @tensor ccd3_tmp[k,l,i,j] :=       (gmo[o,o,v,v])[k,l,c,d] * t_amp[c,d,i,j]
       @tensor ccd3[a,b,i,j]     := 1/4 * ccd3_tmp[k,l,i,j] * t_amp[a,b,k,l]
       println(isapprox(ccd3_ref, ccd3))

       @tensor ccd4a_ref[a,b,i,j] := (gmo[o,o,v,v])[ k,l,c,d] * t_amp[a,c,i,k] * t_amp[b,d,j,l]
       @tensor ccd4a_tmp[l,a,i,d] := (gmo[o,o,v,v])[ k,l,c,d] * t_amp[a,c,i,k]
       @tensor ccd4a[a,b,i,j]     := ccd4a_tmp[l,a,i,d] * t_amp[b,d,j,l]
       println(isapprox(ccd4a_ref, ccd4a))

       ccd4b  = -permutedims(ccd4a, (1, 2, 4, 3))
       ccd4   = ccd4a + ccd4b

       # Update Amplitude
       t_amp_new = @. e_denom * (mp2 + cepa1 + cepa2 + cepa3 + ccd1 + ccd2 + ccd3 + ccd4)

       # Evaluate Energy
       E_CCD = 1/4 * @tensor scalar((gmo[o,o,v,v])[i,j,a,b] * t_amp_new[a,b,i,j])
       t_amp = t_amp_new
       dE = E_CCD - E_old
       printfmt("CCD Iteration {1:3d}: Energy = {2:4.12f} dE = {3:1.5e}\n", cc_iter, E_CCD, dE)

       if abs(dE) < 1.e-8
           @info "CCD Iterations have converged!"
           break
       end

       if cc_iter == MAXITER
           psi4.core.clean()
           error("Maximum number of iterations exceeded.")
       end
   end
   E_CCD
end

printfmt("\nCCD Correlation Energy:    {:15.12f}\n", E_CCD)
printfmt("CCD Total Energy:         {:15.12f}\n", E_CCD + scf_e)
# -

# Unfortunately, Psi4 does not have a CCD code to compare this to. However, Psi4 does have Bruekner CCD, an orbital-optimized variant of CCD. We can qualitatively compare our energies to this energy. The Bruekner-CCD energy should be a little lower than our CCD energy due to the orbital optimization procedure.

psi4_bccd = psi4.energy("bccd", ref_wfn = scf_wfn)
printfmt("\nPsi4 BCCD Correlation Energy:    {:15.12f}\n", psi4_bccd - scf_e)
printfmt("Psi4 BCCD Total Energy:         {:15.12f}\n", psi4_bccd)

# ## References
#
# 1. Modern review of coupled-cluster theory, included diagrammatic derivations of the CCD equations:
# 	> [[Bartlett and Musial:2007](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.79.291)] Rodney J. Bartlett and Monika Musial,  "Coupled-cluster theory in quantum chemistry" *Rev. Mod. Phys.* **79**, 291 (2007)
#
# 2. Background on CEPA:
#     >Kutzelnigg, Werner 1977 *Methods of Electronic Structure Theory* ed. H. F. Schaefer III (Plenum, New York), p 129
#
# 3. More CEPA:
#     > [Koch and Kutzelnigg:1981](https://link.springer.com/article/10.1007/BF00553396) S. Koch and W. Kutzelnigg, *Theor. Chim. Acta* **59**, 387 (1981).
#
# 4. Original CCD Paper:
#     > [Čížek:1966](http://aip.scitation.org/doi/abs/10.1063/1.1727484) Jiří Čížek, "On the Correlation Problem in Atomic and Molecular Systems. Calculation of Wavefunction Components in Ursell‐Type Expansion Using Quantum‐Field Theoretical Methods" *J. Chem. Phys* **45**, 4256 (1966)
#
# 5. Useful notes on diagrams applied to post-HF methods:
#     > A. V. Copan, "Diagram notation" accessed with https://github.com/CCQC/chem-8950/tree/master/2017
#


text:
# # [Shor's Algorithm](@id Shor)

# ## References
# * [Neilsen](https://aapt.scitation.org/doi/abs/10.1119/1.1463744?journalCode=ajp)
# * [An Insightful Blog](https://algassert.com/post/1718)

# The main program of a Shor's algorithm can be summrized in several lines of code.
# For the theory part, please refer the reference materials above.
# It factorize an integer `L`, and returns one of the factors.
using Yao, BitBasis
using Yao.EasyBuild: qft_circuit

# ## Number theory basic
# Before entering the main program, let us defined some useful functions in number theory.

module NumberTheory

export Z_star, Eulerφ, continued_fraction, mod_inverse, rand_primeto, factor_a_power_b
export is_order, order_from_float, find_order

"""
    Z_star(N::Int) -> Vector

returns the Z* group elements of `N`, i.e. {x | gcd(x, N) == 1}
"""
Z_star(N::Int) = filter(i->gcd(i, N)==1, 0:N-1)
Eulerφ(N) = length(Z_star(N))

"""
    continued_fraction(ϕ, niter::Int) -> Rational

obtain `s` and `r` from `ϕ` that satisfies `|s/r - ϕ| ≦ 1/2r²`
"""
continued_fraction(ϕ, niter::Int) = niter==0 ? floor(Int, ϕ) : floor(Int, ϕ) + 1//continued_fraction(1/mod(ϕ, 1), niter-1)
continued_fraction(ϕ::Rational, niter::Int) = niter==0 || ϕ.den==1 ? floor(Int, ϕ) : floor(Int, ϕ) + 1//continued_fraction(1/mod(ϕ, 1), niter-1)

"""
    mod_inverse(x::Int, N::Int) -> Int

Return `y` that `(x*y)%N == 1`, notice the `(x*y)%N` operations in Z* forms a group and this is the definition of inverse.
"""
function mod_inverse(x::Int, N::Int)
    for i=1:N
        (x*i)%N == 1 && return i
    end
    throw(ArgumentError("Can not find the inverse, $x is probably not in Z*($N)!"))
end

"""
    is_order(r, x, N) -> Bool

Returns true if `r` is the order of `x`, i.e. `r` satisfies `x^r % N == 1`.
"""
is_order(r, x, N) = powermod(x, r, N) == 1

"""
    find_order(x::Int, N::Int) -> Int

Find the order of `x` by brute force search.
"""
function find_order(x::Int, N::Int)
    findfirst(r->is_order(r, x, N), 1:N)
end

"""
    rand_primeto(N::Int) -> Int

Returns a random number `2 ≦ x < N` that is prime to `N`.
"""
function rand_primeto(N::Int)
    while true
        x = rand(2:N-1)
        d = gcd(x, N)
        if d == 1
            return x
        end
    end
end

"""
    order_from_float(ϕ, x, L) -> Int

Estimate the order of `x` to `L`, `r`, from a floating point number `ϕ ∼ s/r` using the continued fraction method.
"""
function order_from_float(ϕ, x, L)
    k = 1
    rnum = continued_fraction(ϕ, k)
    while rnum.den < L
        r = rnum.den
        if is_order(r, x, L)
            return r
        end
        k += 1
        rnum = continued_fraction(ϕ, k)
    end
    return nothing
end

"""
    factor_a_power_b(N::Int) -> (Int, Int) or nothing

Factorize `N` into the power form `a^b`.
"""
function factor_a_power_b(N::Int)
    y = log2(N)
    for b = 2:ceil(Int, y)
        x = 2^(y/b)
        u1 = floor(Int, x)
        u1^b == N && return (u1, b)
        (u1+1)^b == N && return (u1+1, b)
    end
end

end

# ## A quantum function to compute `mod`
# Before introducing the main program, let us customize a block for computing the classical function `mod`.
# In a more practical setup, it should be compiled to basic quantum gates. Here we just hack this function for simplicity.

"""
    KMod <: PrimitiveBlock{2}

The first `k` qubits are exponent, and the rest `n-k` are base `a`,
it calculates `mod(a^k*x, L)`, notice `gcd(a, L)` should be 1.
"""
struct KMod <: PrimitiveBlock{2}
    n::Int
    k::Int
    a::Int
    L::Int
    function KMod(n, k, a, L)
        @assert gcd(a, L) == 1 && L<=1<<(n-k)
        new(n, k, a, L)
    end
end

Yao.nqudits(m::KMod) = m.n

function bint2_reader(T, k::Int)
    mask = bmask(T, 1:k)
    return b -> (b&mask, b>>k)
end

function _apply!(reg::ArrayReg, m::KMod)
    nstate = zero(reg.state)

    reader = bint2_reader(Int, m.k)
    for b in basis(reg)
        k, i = reader(b)
        _i = i >= m.L ? i : mod(i*powermod(m.a, m.k, m.L), m.L)
        _b = k + _i<<m.k + 1
        for j in 1:YaoArrayRegister._asint(nbatch(reg))
            @inbounds nstate[_b,j] = reg.state[b+1,j]
        end
    end
    reg.state = nstate
    reg
end

function Yao.mat(::Type{T}, m::KMod) where {T}
    perm = Vector{Int}(undef, 1<<m.n)
    reader = bint2_reader(Int, m.k)
    for b in basis(m.n)
        k, i = reader(b)
        _i = i >= m.L ? i : mod(i*powermod(m.a, k, m.L), m.L)
        _b = k + _i<<m.k + 1
        @inbounds perm[_b] = b+1
    end
    YaoBlocks.LuxurySparse.PermMatrix(perm, ones(T, 1<<m.n))
end

Base.adjoint(m::KMod) = KMod(m.n, m.k, mod_inverse(m.a, m.L), m.L)
Yao.print_block(io::IO, m::KMod) = print(io, "Mod: $(m.a)^k*x % $(m.L) (nqubits = $(nqudits(m)), number of control bits = $(m.k))")

# ## Main Program
# Here, the input `ver` can be either `Val(:quantum)` or `Val(:classical)`,
# where the classical version is for comparison.
using .NumberTheory

function shor(L::Int, ver=Val(:quantum); maxtry=100)
    L%2 == 0 && return 2

    ## find short cut solutions like `a^b`
    res = NumberTheory.factor_a_power_b(L)
    res !== nothing && return res[1]

    for i in 1:maxtry
        ## step 1
        x = NumberTheory.rand_primeto(L)

        ## step 2
        r = get_order(ver, x, L; )
        if r%2 == 0 && powermod(x, r÷2, L) != L-1
            ## step 3
            f1, f2 = gcd(powermod(x, r÷2, L)-1, L), gcd(powermod(x, r÷2, L)+1, L)
            if f1!=1
                return f1
            elseif f2!=1
                return f2
            else
                error("Algorithm Fail!")
            end
        end
    end
end

# Except some shortcuts, in each try, the main program can be summarized in several steps
# 1. randomly pick a number that prime to the input numebr `L`, i.e. `gcd(x, L) = 1`.
# The complexity of this algorithm is polynomial.
# 2. get the order `x`, i.e. finding a number `r` that satisfies `mod(x^r, L) = 1`.
# If `r` is even and `x^(r÷2)` is non-trivial, go on, otherwise start another try.
# Here, trivial means equal to `L-1 (mod L)`.
# 3. According to Theorem 5.2 in Neilsen book,
# one of `gcd(x^(r÷2)-1, L)` and `gcd(x^(r÷2)+1, L)` must be a non-trivial (`!=1`) factor of `L`.
# Notice `powermod(x, r÷2, L)` must be `-1` rather than `1`,
# otherwise the order should be `r/2` according to definition.

# The only difference between classical and quantum version is the order finding algorithm.

# ## Order Finding
# We provided a classical order finding algorithm in `NumberTheory`,
# here we focus on the quantum version.
# The algorithm is consisted
# 1. run the circuit to get a bitstring,
# 2. interpret this bitstring in output register as a rational number `s/r`.
# To achieve this, we first interpret it as a floating point number,
# then the continued fraction algorithm can find the best match for us.
#
# When using the quantum version, we have the flexibility to set key word arguments `nshot`,
# `nbit` (size of input data register) and `ncbit` (size of control register, or output register).
# `nbit` can be simply chosen as the minimum register size to store input,
# while `ncbit` can be estimated with the following function
"""estimate the required size of the output register."""
estimate_ncbit(nbit::Int, ϵ::Real) = 2*nbit + 1 + ceil(Int,log2(2+1/2ϵ))

get_order(::Val{:classical}, x::Int, L::Int; kwargs...) = NumberTheory.find_order(x, L)
function get_order(::Val{:quantum}, x::Int, L::Int; nshots::Int=10,
            nbit::Int=bit_length(L-1), ncbit::Int=estimate_ncbit(nbit, 0.25))
    c = order_finding_circuit(x, L; nbit=nbit, ncbit=ncbit)
    reg = join(product_state(nbit, 1), zero_state(ncbit))

    res = measure(copy(reg) |> c; nshots=nshots)
    for r in res
        ## split bit string b into lower bits `k` and higher bits `r`.
        mask = bmask(1:ncbit)
        k,i = r&mask, r>>ncbit
        ## get s/r
        ϕ = bfloat(k)  #
        ϕ == 0 && continue

        ## order_from_float: given a floating point number,
        ## return the closest rational number with bounded number of continued fraction steps.
        order = NumberTheory.order_from_float(ϕ, x, L)
        if order === nothing
            continue
        else
            return order
        end
    end
    return nothing
end

# #### The circuit used for finding order
"""
    order_finding_circuit(x::Int, L::Int; nbit::Int=bit_length(L-1), ncbit::Int=estimate_ncbit(nbit, 0.25)) -> AbstractBlock

Returns the circuit for finding the order of `x` to `L`,
feeding input `|1>⊗|0>` will get the resulting quantum register with the desired "phase" information.
"""
function order_finding_circuit(x::Int, L::Int; nbit::Int, ncbit::Int)
    N = nbit+ncbit
    chain(N, repeat(N, H, 1:ncbit), KMod(N, ncbit, x, L),
        subroutine(N, qft_circuit(ncbit)', 1:ncbit))
end

# The circuit for order finding is consisted of three parts
# 1. Hadamard gates,
# 2. `KMod` that computes a classical function `mod(a^k*x, L)`.
# `k` is the integer stored in first `K` (or `ncbit`) qubits and the rest `N-K` qubits stores `a`.
# Notice it is not a basic gate, it should have been compiled to multiple gates, which is not implemented in `Yao` for the moment.
# To learn more about implementing arithmatics on a quantum circuit, please read [this paper](https://arxiv.org/abs/1805.12445).
# 3. Inverse quantum fourier transformation.

# ## Run
# Factorizing `15`, you should see `3` or `5`, please report a bug if it is not...
using Random; Random.seed!(129) #src
shor(15, Val(:quantum))


text:
#!/usr/bin/env julia
#=
 The MIT License (MIT)
 Copyright (c) 2018 Michael Hirsch, Ph.D.
=#
using Test

function cohensutherland(xmin, ymax, xmax, ymin, x1, y1, x2, y2)
    #=
    Clips a line to a rectangular area.

    This implements the Cohen-Sutherland line clipping algorithm.  xmin,
    ymax, xmax and ymin denote the clipping area, into which the line
    defined by x1, y1 (start point) and x2, y2 (end point) will be
    clipped.

    If the line intersects with the rectangular clipping area,
    a tuple of the clipped line points will be returned in the form (cx1, cy1, cx2, cy2).
    =#
    INSIDE,LEFT, RIGHT, LOWER, UPPER = 0,1, 2, 4, 8

    function _getclip(xa, ya)
        p = INSIDE  #default is inside

        # consider x
        if xa < xmin
            p |= LEFT
        elseif xa > xmax
            p |= RIGHT
        end

        # consider y
        if ya < ymin
            p |= LOWER # bitwise OR
        elseif ya > ymax
            p |= UPPER #bitwise OR
        end

        return p
     end

# check for trivially outside lines
    k1 = _getclip(x1, y1)
    k2 = _getclip(x2, y2)

#%% examine non-trivially outside points
    #bitwise OR |
    while (k1 | k2) != 0 # if both points are inside box (0000) , ACCEPT trivial whole line in box

        # if line trivially outside window, REJECT
        if (k1 & k2) != 0 #bitwise AND &
            return nothing
        end

        # non-trivial case, at least one point outside window

        # this is NOT a bitwise or
        if k1 != 0
            opt = k1
        elseif k2 != 0
            opt = k2
        else
            jl_error("unexpected k1,k2  $k1  $k2")
        end

        if (opt & UPPER) != 0 #these are bitwise ANDS
            x = x1 + (x2 - x1) * (ymax - y1) / (y2 - y1)
            y = ymax
        elseif (opt & LOWER) != 0
            x = x1 + (x2 - x1) * (ymin - y1) / (y2 - y1)
            y = ymin
        elseif (opt & RIGHT) != 0
            y = y1 + (y2 - y1) * (xmax - x1) / (x2 - x1)
            x = xmax
        elseif (opt & LEFT) != 0
            y = y1 + (y2 - y1) * (xmin - x1) / (x2 - x1)
            x = xmin
        else
            jl_error("Undefined clipping state")
        end

        if opt == k1
            x1, y1 = x, y
            k1 = _getclip(x1, y1)
            #println("checking k1: $x , $y   $k1")
        elseif opt == k2
            x2, y2 = x, y
            k2 = _getclip(x2, y2)
            #println("checking 2: $x , $y   $k2")
        end

    end

    return x1, y1, x2, y2
end

if basename(PROGRAM_FILE) == basename(@__FILE__)
    xmin, ymax, xmax, ymin, x1, y1, x2, y2 = 1,  5, 4, 3, 0,  0, 4, 6

    x1, y1, x2, y2 = cohensutherland(xmin, ymax, xmax, ymin, x1, y1, x2, y2)

    @test x1≈2
    @test y1≈3
    @test x2≈3.3333333
    @test y2≈5

    println("x1,y2,x2,y2 $x1 $y1 $x2 $y2")
end

max_stars_repo_name: scivision/lineclipping-julia

max_issues_repo_name: scivision/lineclipping-julia

max_forks_repo_name: scivision/lineclipping-julia


text:
#######################################################################
# SystemInfo
#######################################################################

struct SystemInfo{Tv,S}
    sampledterms::S
    namesdict::Dict{NameType,Int}
    names::Vector{NameType}
    nsites::Vector{Int}
    norbitals::Vector{Int}
    dims::Vector{Int}       # Hamiltonian block dimensions for each sublattice
    offsets::Vector{Int}    # Hamiltonian block offset for each sublattice
end

SystemInfo{Tv}(terms::S, args...) where {Tv,S} = SystemInfo{Tv,S}(terms, args...)

function SystemInfo(lat::Lattice{E,L,T}, model::Model{Tv}, prevsamples...) where {E,L,T,Tv}
    ns = length(lat.sublats)
    nsites = Vector{Int}(undef, ns)
    namesdict = Dict{NameType,Int}()
    names = Vector{NameType}(undef, ns)
    zeropos = zero(SVector{E,T})
    for (i, sublat) in enumerate(lat.sublats)
        nsites[i] = length(sublat.sites)
        namesdict[sublat.name] = i
        names[i] = sublat.name
    end
    sampledterms = getsamples(namesdict, zeropos, model.terms...)
    allterms = (prevsamples..., sampledterms...)
    norbitals = zeros(Int, ns)
    _fillorbitals!(norbitals, allterms...)
    dims = zeros(Int, ns)
    offsets = zeros(Int, ns + 1)
    for i in eachindex(lat.sublats)
        dims[i] = norbitals[i] * nsites[i]
        offsets[i + 1] = offsets[i] + dims[i]
    end
    return SystemInfo{Tv}(allterms, namesdict, names, nsites, norbitals, dims, offsets)
end

getsamples(namesdict, zeropos) = ()
getsamples(namesdict, zeropos, term, terms...) =
    (_getsamples(namesdict, zeropos, term), getsamples(namesdict, zeropos, terms...)...)

_getsamples(namesdict, zeropos, term) =
    (term, term(zeropos, zeropos), _getsublats(namesdict, term.sublats))

_getsublats(namesdict, s::Missing) = s
_getsublats(namesdict, s::NTuple{N,Any}) where {N} =
    ntuple(i -> (sublatindex(namesdict, first(s[i])),
                 sublatindex(namesdict, last(s[i]))), Val(N))

_fillorbitals!(norbitals) = nothing

function _fillorbitals!(norbitals, sample::Tuple, samples...)
    _fillorbitals!(norbitals, sample...)
    _fillorbitals!(norbitals, samples...)
end

function _fillorbitals!(norbitals, term, ::SMatrix{N,M}, ::Missing) where {N,M}
    N == M ? fill!(norbitals, N) : throw(DimensionMismatch("Inconsitent model orbital dimensions"))
    return nothing
end

_fillorbitals!(norbitals, term, sm::SMatrix, ss::NTuple{N,Tuple{Int,Int}}) where {N} =
    foreach(s -> _fillorbitals!(norbitals, sm, s), ss)
_fillorbitals!(norbitals, ::SMatrix{N,M}, (s1, s2)::Tuple{Int,Int}) where {N,M} =
    (_fillorbitals!(norbitals, Val(N), s1); _fillorbitals!(norbitals, Val(M), s2))

function _fillorbitals!(norbitals, ::Val{N}, s) where {N}
    0 < s <= length(norbitals) || return nothing
    if norbitals[s] == 0
        norbitals[s] = N
    elseif norbitals[s] != N
        throw(DimensionMismatch("Inconsitent model orbital dimensions"))
    end
    return nothing
end

function tosite(row, sysinfo)
    s = findsublat(row, sysinfo.offsets)
    offset = sysinfo.offsets[s]
    norbs = sysinfo.norbitals[s]
    deltaoffset = row - offset - 1
    site = div(deltaoffset, norbs) + 1
    orboffset = rem(deltaoffset, norbs) + 1
    return site, orboffset, s
end

torow(siteindex, sublat, sysinfo) =
    sysinfo.offsets[sublat] + (siteindex - 1) * sysinfo.norbitals[sublat] + 1

function findsublat(row, offsets)
    for n in eachindex(offsets)
        offsets[n] >= row && return n - 1
    end
    return 0
end

#######################################################################
# Operator
#######################################################################

mutable struct Block{Tv,L}
    ndist::SVector{L,Int}
    matrix::SparseMatrixCSC{Tv,Int}
    nlinks::Int
end

function Block(ndist, matrix, sysinfo::SystemInfo)
    b = Block(ndist, matrix, 0)
    isempty(matrix) || updatenlinks!(b, sysinfo)
    return b
end

Base.isempty(b::Block) = isempty(b.matrix)

Base.zero(b::Block{Tv,L}) where {Tv,L} =
    Block(zero(SVector{L,Int}), spzeros(Tv, size(b.matrix)...), 0)

function Base.show(io::IO, b::Block{Tv,L}) where {Tv,L}
    print(io, "Block{$Tv,$L}: Bloch harmonic $(b.ndist) of dimensions $(size(b.matrix)) with $(nnz(b.matrix)) elements")
end

struct Operator{Tv,L}
    matrix::SparseMatrixCSC{Tv,Int}
    intra::Block{Tv,L}
    inters::Vector{Block{Tv,L}}
    boundary::Vector{Tuple{Int,Int}} # ptrs to matrix and intra that contain terms in any inters
end

function Operator{Tv,L}(n::Int) where {Tv,L}
    return Operator(
        spzeros(n,n),
        Block{Tv,L}(zero(SVector{L,Int}), n),
        Block{Tv,L}[],
        Tuple{Int,Int}[])
end

function Base.show(io::IO, op::Operator{L,Tv}) where {L,Tv}
    print(io, "Operator{$L,$Tv}: Bloch operator of dimensions $(size(op.matrix)) with $(nnz(op.matrix)) elements in main matrix")
end

#######################################################################
# Operator API
#######################################################################

nlinks(o::Operator) = nlinks(o.intra) + (isempty(o.inters) ? 0 : sum(nlinks, o.inters))
nlinks(b::Block) = b.nlinks

nsublats(b::Block) = nsublats(b.sysinfo)

function updatenlinks!(b::Block, sysinfo)
    n = 0
    zeron = iszero(b.ndist)
    for (_, (target, source), (row, col), _) in BlockIterator(b, sysinfo)
        (!zeron || row != col) && (n += 1)
    end
    b.nlinks = n
    return nothing
end

insertblochphases!(o::Operator{Tv,L}, kn) where {Tv<:AbstractFloat,L} =
    throw(DomainError(Tv, "Cannot apply Bloch phases to a real Hamiltonian."))

function insertblochphases!(op::Operator{Tv,L}, ϕn, dϕaxis = missing) where {Tv,L}
    length(ϕn) == L || throw(DimensionMismatch(
        "The dimension of the normalized Bloch phases should match the lattice dimension $L"))
    rows = rowvals(op.matrix)
    vals = nonzeros(op.matrix)
    valsintra = nonzeros(op.intra.matrix)
    for (ptrmatrix, ptrintra) in op.boundary
        vals[ptrmatrix] = iszero(ptrintra) ? zero(Tv) : valsintra[ptrintra]
    end
    for inter in op.inters
        # dϕaxis != missing produces derivatives of Bloch phases respect a given axis
        dexp = dϕaxis === missing ? 1.0 : 2pi * im * dot(inter.ndist, dϕaxis)
        phase = dexp * exp(2pi * im * dot(inter.ndist, ϕn))
        for (i,j,v,ptr) in SparseMatrixReader(inter.matrix)
            for ptr in nzrange(op.matrix,j)
                rows[ptr] == i && (vals[ptr] += phase * v; break)
            end
        end
    end
    return op
end

function boundaryoperator(op::Operator{Tv}) where {Tv}
    n = length(op.boundary)
    sb = SparseMatrixBuilder{Tv}(size(op.matrix)...)
    rows = rowvals(op.matrix)
    vals = nonzeros(op.matrix)
    sofar = 1
    for col in 1:size(op.matrix, 2)
        colrange = nzrange(op.matrix, col)
        for bidx in sofar:n
            ptr = first(op.boundary[bidx])
            if ptr in colrange
                sofar = bidx + 1
                pushtocolumn!(sb, rows[ptr], vals[ptr])
            else
                break
            end
        end
        finalizecolumn!(sb)
    end
    matrix = sparse(sb)
    intra = zero(op.intra)
    inters = op.inters
    boundary = extractboundary(matrix, intra, inters)
    return Operator(matrix, intra, inters, boundary)
end


text:
using FNVHash
using Base.Test

# taken from http://www.isthe.com/chongo/tech/comp/fnv/#zero-hash challenge 7
@test fnv32([0xcc,0x24,0x31,0xc4]) == 0
@test fnv32([0xe0, 0x4d, 0x9f, 0xcb]) == 0

TEST(s) = s

function TEST0(s)
    a = Vector{UInt8}(s)
    push!(a, 0x00)
end

R10(s::AbstractString) = s^10
R500(s::AbstractString) = s^500

R10(s) = repeat(s, outer=[10])
R500(s) = repeat(s, outer=[500])

input = Any[
  TEST(""),
  TEST("a"),
  TEST("b"),
  TEST("c"),
  TEST("d"),
  TEST("e"),
  TEST("f"),
  TEST("fo"),
  TEST("foo"),
  TEST("foob"),
  TEST("fooba"),
  TEST("foobar"),
  TEST0(""),
  TEST0("a"),
  TEST0("b"),
  TEST0("c"),
  TEST0("d"),
  TEST0("e"),
  TEST0("f"),
  TEST0("fo"),
  TEST0("foo"),
  TEST0("foob"),
  TEST0("fooba"),
  TEST0("foobar"),
  TEST("ch"),
  TEST("cho"),
  TEST("chon"),
  TEST("chong"),
  TEST("chongo"),
  TEST("chongo "),
  TEST("chongo w"),
  TEST("chongo wa"),
  TEST("chongo was"),
  TEST("chongo was "),
  TEST("chongo was h"),
  TEST("chongo was he"),
  TEST("chongo was her"),
  TEST("chongo was here"),
  TEST("chongo was here!"),
  TEST("chongo was here!\n"),
  TEST0("ch"),
  TEST0("cho"),
  TEST0("chon"),
  TEST0("chong"),
  TEST0("chongo"),
  TEST0("chongo "),
  TEST0("chongo w"),
  TEST0("chongo wa"),
  TEST0("chongo was"),
  TEST0("chongo was "),
  TEST0("chongo was h"),
  TEST0("chongo was he"),
  TEST0("chongo was her"),
  TEST0("chongo was here"),
  TEST0("chongo was here!"),
  TEST0("chongo was here!\n"),
  TEST("cu"),
  TEST("cur"),
  TEST("curd"),
  TEST("curds"),
  TEST("curds "),
  TEST("curds a"),
  TEST("curds an"),
  TEST("curds and"),
  TEST("curds and "),
  TEST("curds and w"),
  TEST("curds and wh"),
  TEST("curds and whe"),
  TEST("curds and whey"),
  TEST("curds and whey\n"),
  TEST0("cu"),
  TEST0("cur"),
  TEST0("curd"),
  TEST0("curds"),
  TEST0("curds "),
  TEST0("curds a"),
  TEST0("curds an"),
  TEST0("curds and"),
  TEST0("curds and "),
  TEST0("curds and w"),
  TEST0("curds and wh"),
  TEST0("curds and whe"),
  TEST0("curds and whey"),
  TEST0("curds and whey\n"),
  TEST("hi"), TEST0("hi"),
  TEST("hello"), TEST0("hello"),
  TEST(b"\xff\x00\x00\x01"), TEST(b"\x01\x00\x00\xff"),
  TEST(b"\xff\x00\x00\x02"), TEST(b"\x02\x00\x00\xff"),
  TEST(b"\xff\x00\x00\x03"), TEST(b"\x03\x00\x00\xff"),
  TEST(b"\xff\x00\x00\x04"), TEST(b"\x04\x00\x00\xff"),
  TEST(b"\x40\x51\x4e\x44"), TEST(b"\x44\x4e\x51\x40"),
  TEST(b"\x40\x51\x4e\x4a"), TEST(b"\x4a\x4e\x51\x40"),
  TEST(b"\x40\x51\x4e\x54"), TEST(b"\x54\x4e\x51\x40"),
  TEST("127.0.0.1"), TEST0("127.0.0.1"),
  TEST("127.0.0.2"), TEST0("127.0.0.2"),
  TEST("127.0.0.3"), TEST0("127.0.0.3"),
  TEST("64.81.78.68"), TEST0("64.81.78.68"),
  TEST("64.81.78.74"), TEST0("64.81.78.74"),
  TEST("64.81.78.84"), TEST0("64.81.78.84"),
  TEST("feedface"), TEST0("feedface"),
  TEST("feedfacedaffdeed"), TEST0("feedfacedaffdeed"),
  TEST("feedfacedeadbeef"), TEST0("feedfacedeadbeef"),
  TEST("line 1\nline 2\nline 3"),
  TEST("chongo <Landon Curt Noll> /\\../\\"),
  TEST0("chongo <Landon Curt Noll> /\\../\\"),
  TEST("chongo (Landon Curt Noll) /\\../\\"),
  TEST0("chongo (Landon Curt Noll) /\\../\\"),
  TEST("http://antwrp.gsfc.nasa.gov/apod/astropix.html"),
  TEST("http://en.wikipedia.org/wiki/Fowler_Noll_Vo_hash"),
  TEST("http://epod.usra.edu/"),
  TEST("http://exoplanet.eu/"),
  TEST("http://hvo.wr.usgs.gov/cam3/"),
  TEST("http://hvo.wr.usgs.gov/cams/HMcam/"),
  TEST("http://hvo.wr.usgs.gov/kilauea/update/deformation.html"),
  TEST("http://hvo.wr.usgs.gov/kilauea/update/images.html"),
  TEST("http://hvo.wr.usgs.gov/kilauea/update/maps.html"),
  TEST("http://hvo.wr.usgs.gov/volcanowatch/current_issue.html"),
  TEST("http://neo.jpl.nasa.gov/risk/"),
  TEST("http://norvig.com/21-days.html"),
  TEST("http://primes.utm.edu/curios/home.php"),
  TEST("http://slashdot.org/"),
  TEST("http://tux.wr.usgs.gov/Maps/155.25-19.5.html"),
  TEST("http://volcano.wr.usgs.gov/kilaueastatus.php"),
  TEST("http://www.avo.alaska.edu/activity/Redoubt.php"),
  TEST("http://www.dilbert.com/fast/"),
  TEST("http://www.fourmilab.ch/gravitation/orbits/"),
  TEST("http://www.fpoa.net/"),
  TEST("http://www.ioccc.org/index.html"),
  TEST("http://www.isthe.com/cgi-bin/number.cgi"),
  TEST("http://www.isthe.com/chongo/bio.html"),
  TEST("http://www.isthe.com/chongo/index.html"),
  TEST("http://www.isthe.com/chongo/src/calc/lucas-calc"),
  TEST("http://www.isthe.com/chongo/tech/astro/venus2004.html"),
  TEST("http://www.isthe.com/chongo/tech/astro/vita.html"),
  TEST("http://www.isthe.com/chongo/tech/comp/c/expert.html"),
  TEST("http://www.isthe.com/chongo/tech/comp/calc/index.html"),
  TEST("http://www.isthe.com/chongo/tech/comp/fnv/index.html"),
  TEST("http://www.isthe.com/chongo/tech/math/number/howhigh.html"),
  TEST("http://www.isthe.com/chongo/tech/math/number/number.html"),
  TEST("http://www.isthe.com/chongo/tech/math/prime/mersenne.html"),
  TEST("http://www.isthe.com/chongo/tech/math/prime/mersenne.html#largest"),
  TEST("http://www.lavarnd.org/cgi-bin/corpspeak.cgi"),
  TEST("http://www.lavarnd.org/cgi-bin/haiku.cgi"),
  TEST("http://www.lavarnd.org/cgi-bin/rand-none.cgi"),
  TEST("http://www.lavarnd.org/cgi-bin/randdist.cgi"),
  TEST("http://www.lavarnd.org/index.html"),
  TEST("http://www.lavarnd.org/what/nist-test.html"),
  TEST("http://www.macosxhints.com/"),
  TEST("http://www.mellis.com/"),
  TEST("http://www.nature.nps.gov/air/webcams/parks/havoso2alert/havoalert.cfm"),
  TEST("http://www.nature.nps.gov/air/webcams/parks/havoso2alert/timelines_24.cfm"),
  TEST("http://www.paulnoll.com/"),
  TEST("http://www.pepysdiary.com/"),
  TEST("http://www.sciencenews.org/index/home/activity/view"),
  TEST("http://www.skyandtelescope.com/"),
  TEST("http://www.sput.nl/~rob/sirius.html"),
  TEST("http://www.systemexperts.com/"),
  TEST("http://www.tq-international.com/phpBB3/index.php"),
  TEST("http://www.travelquesttours.com/index.htm"),
  TEST("http://www.wunderground.com/global/stations/89606.html"),
  TEST(R10("21701")),
  TEST(R10("M21701")),
  TEST(R10("2^21701-1")),
  TEST(R10(b"\x54\xc5")),
  TEST(R10(b"\xc5\x54")),
  TEST(R10("23209")),
  TEST(R10("M23209")),
  TEST(R10("2^23209-1")),
  TEST(R10(b"\x5a\xa9")),
  TEST(R10(b"\xa9\x5a")),
  TEST(R10("391581216093")),
  TEST(R10("391581*2^216093-1")),
  TEST(R10(b"\x05\xf9\x9d\x03\x4c\x81")),
  TEST(R10("FEDCBA9876543210")),
  TEST(R10(b"\xfe\xdc\xba\x98\x76\x54\x32\x10")),
  TEST(R10("EFCDAB8967452301")),
  TEST(R10(b"\xef\xcd\xab\x89\x67\x45\x23\x01")),
  TEST(R10("0123456789ABCDEF")),
  TEST(R10(b"\x01\x23\x45\x67\x89\xab\xcd\xef")),
  TEST(R10("1032547698BADCFE")),
  TEST(R10(b"\x10\x32\x54\x76\x98\xba\xdc\xfe")),
  TEST(R500(b"\x00")),
  TEST(R500(b"\x07")),
  TEST(R500("~")),
  TEST(R500(b"\x7f"))
]

# Taken from http://www.isthe.com/chongo/src/fnv/test_fnv.c fnv1_32_test_vector
results = [
    0x811c9dc5,
    0x050c5d7e,
    0x050c5d7d,
    0x050c5d7c,
    0x050c5d7b,
    0x050c5d7a,
    0x050c5d79,
    0x6b772514,
    0x408f5e13,
    0xb4b1178b,
    0xfdc80fb0,
    0x31f0b262,
    0x050c5d1f,
    0x70772d5a,
    0x6f772bc7,
    0x6e772a34,
    0x6d7728a1,
    0x6c77270e,
    0x6b77257b,
    0x408f5e7c,
    0xb4b117e9,
    0xfdc80fd1,
    0x31f0b210,
    0xffe8d046,
    0x6e772a5c,
    0x4197aebb,
    0xfcc8100f,
    0xfdf147fa,
    0xbcd44ee1,
    0x23382c13,
    0x846d619e,
    0x1630abdb,
    0xc99e89b2,
    0x1692c316,
    0x9f091bca,
    0x2556be9b,
    0x628e0e73,
    0x98a0bf6c,
    0xb10d5725,
    0xdd002f35,
    0x4197aed4,
    0xfcc81061,
    0xfdf1479d,
    0xbcd44e8e,
    0x23382c33,
    0x846d61e9,
    0x1630abba,
    0xc99e89c1,
    0x1692c336,
    0x9f091ba2,
    0x2556befe,
    0x628e0e01,
    0x98a0bf09,
    0xb10d5704,
    0xdd002f3f,
    0x1c4a506f,
    0x6e772a41,
    0x26978421,
    0xe184ff97,
    0x9b5e5ac6,
    0x5b88e592,
    0xaa8164b7,
    0x20b18c7b,
    0xf28025c5,
    0x84bb753f,
    0x3219925a,
    0x384163c6,
    0x54f010d7,
    0x8cea820c,
    0xe12ab8ee,
    0x26978453,
    0xe184fff3,
    0x9b5e5ab5,
    0x5b88e5b2,
    0xaa8164d6,
    0x20b18c15,
    0xf28025a1,
    0x84bb751f,
    0x3219922d,
    0x384163ae,
    0x54f010b2,
    0x8cea8275,
    0xe12ab8e4,
    0x64411eaa,
    0x6977223c,
    0x428ae474,
    0xb6fa7167,
    0x73408525,
    0xb78320a1,
    0x0caf4135,
    0xb78320a2,
    0xcdc88e80,
    0xb78320a3,
    0x8ee1dbcb,
    0xb78320a4,
    0x4ffb2716,
    0x860632aa,
    0xcc2c5c64,
    0x860632a4,
    0x2a7ec4a6,
    0x860632ba,
    0xfefe8e14,
    0x0a3cffd8,
    0xf606c108,
    0x0a3cffdb,
    0xf906c5c1,
    0x0a3cffda,
    0xf806c42e,
    0xc07167d7,
    0xc9867775,
    0xbf716668,
    0xc78435b8,
    0xc6717155,
    0xb99568cf,
    0x7662e0d6,
    0x33a7f0e2,
    0xc2732f95,
    0xb053e78f,
    0x3a19c02a,
    0xa089821e,
    0x31ae8f83,
    0x995fa9c4,
    0x35983f8c,
    0x5036a251,
    0x97018583,
    0xb4448d60,
    0x025dfe59,
    0xc5eab3af,
    0x7d21ba1e,
    0x7704cddb,
    0xd0071bfe,
    0x0ff3774c,
    0xb0fea0ea,
    0x58177303,
    0x4f599cda,
    0x3e590a47,
    0x965595f8,
    0xc37f178d,
    0x9711dd26,
    0x23c99b7f,
    0x6e568b17,
    0x43f0245b,
    0xbcb7a001,
    0x12e6dffe,
    0x0792f2d6,
    0xb966936b,
    0x46439ac5,
    0x728d49af,
    0xd33745c9,
    0xbc382a57,
    0x4bda1d31,
    0xce35ccae,
    0x3b6eed94,
    0x445c9c58,
    0x3db8bf9d,
    0x2dee116d,
    0xc18738da,
    0x5b156176,
    0x2aa7d593,
    0xb2409658,
    0xe1489528,
    0xfe1ee07e,
    0xe8842315,
    0x3a6a63a2,
    0x06d2c18c,
    0xf8ef7225,
    0x843d3300,
    0xbb24f7ae,
    0x878c0ec9,
    0xb557810f,
    0x57423246,
    0x87f7505e,
    0xbb809f20,
    0x8932abb5,
    0x0a9b3aa0,
    0xb8682a24,
    0xa7ac1c56,
    0x11409252,
    0xa987f517,
    0xf309e7ed,
    0xc9e8f417,
    0x7f447bdd,
    0xb929adc5,
    0x57022879,
    0xdcfd2c49,
    0x6edafff5,
    0xf04fb1f1,
    0xfb7de8b9,
    0xc5f1d7e9,
    0x32c1f439,
    0x7fd3eb7d,
    0x81597da5,
    0x05eb7a25,
    0x9c0fa1b5,
    0x53ccb1c5,
    0xfabece15,
    0x4ad745a5,
    0xe5bdc495,
    0x23b3c0a5,
    0xfa823dd5,
    0x0c6c58b9,
    0xe2dbccd5,
    0xdb7f50f9
]

# Taken from http://www.isthe.com/chongo/src/fnv/test_fnv.c fnv1a_32_test_vector
results1a = [
    0x811c9dc5,
    0xe40c292c,
    0xe70c2de5,
    0xe60c2c52,
    0xe10c2473,
    0xe00c22e0,
    0xe30c2799,
    0x6222e842,
    0xa9f37ed7,
    0x3f5076ef,
    0x39aaa18a,
    0xbf9cf968,
    0x050c5d1f,
    0x2b24d044,
    0x9d2c3f7f,
    0x7729c516,
    0xb91d6109,
    0x931ae6a0,
    0x052255db,
    0xbef39fe6,
    0x6150ac75,
    0x9aab3a3d,
    0x519c4c3e,
    0x0c1c9eb8,
    0x5f299f4e,
    0xef8580f3,
    0xac297727,
    0x4546b9c0,
    0xbd564e7d,
    0x6bdd5c67,
    0xdd77ed30,
    0xf4ca9683,
    0x4aeb9bd0,
    0xe0e67ad0,
    0xc2d32fa8,
    0x7f743fb7,
    0x6900631f,
    0xc59c990e,
    0x448524fd,
    0xd49930d5,
    0x1c85c7ca,
    0x0229fe89,
    0x2c469265,
    0xce566940,
    0x8bdd8ec7,
    0x34787625,
    0xd3ca6290,
    0xddeaf039,
    0xc0e64870,
    0xdad35570,
    0x5a740578,
    0x5b004d15,
    0x6a9c09cd,
    0x2384f10a,
    0xda993a47,
    0x8227df4f,
    0x4c298165,
    0xfc563735,
    0x8cb91483,
    0x775bf5d0,
    0xd5c428d0,
    0x34cc0ea3,
    0xea3b4cb7,
    0x8e59f029,
    0x2094de2b,
    0xa65a0ad4,
    0x9bbee5f4,
    0xbe836343,
    0x22d5344e,
    0x19a1470c,
    0x4a56b1ff,
    0x70b8e86f,
    0x0a5b4a39,
    0xb5c3f670,
    0x53cc3f70,
    0xc03b0a99,
    0x7259c415,
    0x4095108b,
    0x7559bdb1,
    0xb3bf0bbc,
    0x2183ff1c,
    0x2bd54279,
    0x23a156ca,
    0x64e2d7e4,
    0x683af69a,
    0xaed2346e,
    0x4f9f2cab,
    0x02935131,
    0xc48fb86d,
    0x2269f369,
    0xc18fb3b4,
    0x50ef1236,
    0xc28fb547,
    0x96c3bf47,
    0xbf8fb08e,
    0xf3e4d49c,
    0x32179058,
    0x280bfee6,
    0x30178d32,
    0x21addaf8,
    0x4217a988,
    0x772633d6,
    0x08a3d11e,
    0xb7e2323a,
    0x07a3cf8b,
    0x91dfb7d1,
    0x06a3cdf8,
    0x6bdd3d68,
    0x1d5636a7,
    0xd5b808e5,
    0x1353e852,
    0xbf16b916,
    0xa55b89ed,
    0x3c1a2017,
    0x0588b13c,
    0xf22f0174,
    0xe83641e1,
    0x6e69b533,
    0xf1760448,
    0x64c8bd58,
    0x97b4ea23,
    0x9a4e92e6,
    0xcfb14012,
    0xf01b2511,
    0x0bbb59c3,
    0xce524afa,
    0xdd16ef45,
    0x60648bb3,
    0x7fa4bcfc,
    0x5053ae17,
    0xc9302890,
    0x956ded32,
    0x9136db84,
    0xdf9d3323,
    0x32bb6cd0,
    0xc8f8385b,
    0xeb08bfba,
    0x62cc8e3d,
    0xc3e20f5c,
    0x39e97f17,
    0x7837b203,
    0x319e877b,
    0xd3e63f89,
    0x29b50b38,
    0x5ed678b8,
    0xb0d5b793,
    0x52450be5,
    0xfa72d767,
    0x95066709,
    0x7f52e123,
    0x76966481,
    0x063258b0,
    0x2ded6e8a,
    0xb07d7c52,
    0xd0c71b71,
    0xf684f1bd,
    0x868ecfa8,
    0xf794f684,
    0xd19701c3,
    0x346e171e,
    0x91f8f676,
    0x0bf58848,
    0x6317b6d1,
    0xafad4c54,
    0x0f25681e,
    0x91b18d49,
    0x7d61c12e,
    0x5147d25c,
    0x9a8b6805,
    0x4cd2a447,
    0x1e549b14,
    0x2fe1b574,
    0xcf0cd31e,
    0x6c471669,
    0x0e5eef1e,
    0x2bed3602,
    0xb26249e0,
    0x2c9b86a4,
    0xe415e2bb,
    0x18a98d1d,
    0xb7df8b7b,
    0x241e9075,
    0x063f70dd,
    0x0295aed9,
    0x56a7f781,
    0x253bc645,
    0x46610921,
    0x7c1577f9,
    0x512b2851,
    0x76823999,
    0xc0586935,
    0xf3415c85,
    0x0ae4ff65,
    0x58b79725,
    0xdea43aa5,
    0x2bb3be35,
    0xea777a45,
    0x8f21c305,
    0x5c9d0865,
    0xfa823dd5,
    0x21a27271,
    0x83c5c6d5,
    0x813b0881
]

@testset fnv1 begin
    for (i, r) in zip(input, results)
        @test fnv32(i) == r
    end
end

@testset fnv1a begin
    for (i, r) in zip(input, results1a)
        @test fnv1a32(i) == r
    end
end


text:
@testset "devices" begin

name(dev)
totalmem(dev)
attribute(dev, CUDAdrv.MAX_THREADS_PER_BLOCK)
@test warpsize(dev) == Int32(32)
capability(dev)
@grab_output show(stdout, "text/plain", dev)

end


text:
@afterscenario() do context, scenario
    delete_all_databases(context[:client])
end


text:
module PyPlotSupport

using Stochy
import Stochy: ERP, support

export histogram, bar

pyplotloaded() = isdefined(Main, :PyPlot)
del(names...) = [Main.PyPlot.(n) for n in names]

function bar(erp::ERP; kwargs...)
    @assert pyplotloaded()
    pybar, xticks, xlim, figure, xlabel, ylabel = del(:bar, :xticks, :xlim, :figure, :xlabel, :ylabel)
    xs = sort(collect(support(erp)))
    ps = map(x->exp(score(erp,x)), xs)
    N = length(xs)
    spacing = 0.05
    width = 1 - 2*spacing
    figure(figsize=(6,4))
    pybar([1:N]+spacing, ps, width; kwargs...)
    xticks([1:N]+0.5, xs)
    xlim(1,N+1)
    xlabel("x")
    ylabel("p(x)")
end

function histogram(erp::ERP; kwargs...)
    @assert pyplotloaded()
    hist, figure, xlabel, ylabel = del(:hist, :figure, :xlabel, :ylabel)
    figure(figsize=(6,4))
    xs = sort(collect(support(erp)))
    ps = map(x->exp(score(erp,x)), xs)
    hist(;x=xs, weights=ps, kwargs...)
    xlabel("x")
    ylabel("count")
end

end # module


text:
using JuliaDB, Compat, GroupedErrors, SputnikUtilities, Images, NamedTuples
using StatPlots
using Compat.Test
school = loadtable(GroupedErrors.exampletable("school.csv"))

function compare_plots(plt1, plt2; sigma = [1,1], eps = 0.02)
    ref1 = joinpath(@__DIR__, "plots", "plot1.png")
    savefig(plt1, ref1)
    ref2 = joinpath(@__DIR__, "plots", "plot2.png")
    savefig(plt2, ref2)
    outcome = Images.test_approx_eq_sigma_eps(Images.load(ref1), Images.load(ref2), sigma, eps)
    rm(ref1)
    rm(ref2)
    outcome
end

@testset "select" begin
    selectdiscrete = (SelectValues(:Minrty, ["Yes"], true), SelectValues(:Sx, ["Male"], false))
    selectcontinuous = (SelectPredicate(:MAch, x -> 6 <= x <= 10),)
    expected = filter(i -> i.Minrty == "Yes" && i.Sx == "Male" && 6 <= i.MAch <= 10, school)
    @test SputnikUtilities.selectdata(school, selectdiscrete, selectcontinuous) == expected
    d2s  = Data2Select(school, selectdiscrete, selectcontinuous)
    @test SelectedData(d2s) == SelectedData(expected, (:Minrty,))
end

@testset "statplots" begin
    selectdiscrete = (SelectValues(:Minrty, String[], true, false), SelectValues(:Sx, ["Male"], false))
    selectcontinuous = (SelectPredicate(:MAch, x -> 6 <= x <= 10),)
    expected = filter(i -> i.Sx == "Male" && 6 <= i.MAch <= 10, school)
    d2s  = Data2Select(school, selectdiscrete, selectcontinuous)
    sd = SelectedData(d2s)

    a = Analysis(data = sd, x = :MAch, plot = density)
    plt1 = process(a)
    plt2 = @df expected density(:MAch, group = {:Minrty})
    @test compare_plots(plt1, plt2) < 0.001

    a = Analysis(data = sd, x = :MAch, y = :SSS, plot = scatter, plot_kwargs = [(:legend, :topleft)])
    plt1 = process(a)
    plt2 = @df expected scatter(:MAch, :SSS, legend = :topleft, group ={:Minrty})
    @test compare_plots(plt1, plt2) < 0.001

    sd = SelectedData(Data2Select(school, (), ()))
    a = Analysis(data = sd, x = :MAch, y = :SSS, plot = scatter)
    plt1 = process(a)
    plt2 = @df school scatter(:MAch, :SSS, legend = false)
    @test compare_plots(plt1, plt2) < 0.001
end

@testset "groupederrors" begin
    school = GroupedErrors.exampletable("school.csv") |> loadtable
    selectdiscrete = (SelectValues(:Minrty, ["Yes", "No"], true), SelectValues(:Sx, ["Male"], false))
    selectcontinuous = (SelectPredicate(:MAch, x -> 6 <= x <= 10),)
    expected = filter(i -> i.Sx == "Male" && 6 <= i.MAch <= 10, school)
    d2s  = Data2Select(school, selectdiscrete, selectcontinuous)
    sd = SelectedData(d2s)
    a = Analysis(data = sd,
                 x = :MAch,
                 y = :density,
                 plot = plot,
                 axis_type = :continuous,
                 smoother = 50.0)
    plt1 = process(a)
    plt2 = @> expected begin
        @splitby _.Minrty
        @x _.MAch :continuous
        @y :density bandwidth = (51.0)*std(column(expected, :MAch))/200
        @plot plot()
    end
    @test compare_plots(plt1, plt2) < 0.001

    a = Analysis(data = sd,
                 compute_error = (:across, :School),
                 x = :MAch,
                 y = :density,
                 plot = plot,
                 axis_type = :continuous,
                 smoother = 50.0)
    plt1 = process(a)
    plt2 = @> expected begin
        @across _.School
        @splitby _.Minrty
        @x _.MAch :continuous
        @y :density bandwidth = (51.0)*std(column(expected, :MAch))/200
        @plot plot()
    end
    @test compare_plots(plt1, plt2) < 0.001

    a = Analysis(data = sd,
                 compute_error = (:across, :School),
                 x = :MAch,
                 y = :SSS,
                 xfunc = mean,
                 yfunc = mean,
                 plot = scatter,
                 axis_type = :pointbypoint)
    plt1 = process(a)
    plt2 = @> expected begin
        @across _.School
        @splitby _.Minrty
        @x _.MAch
        @y _.SSS
        @plot scatter()
    end
    @test compare_plots(plt1, plt2) < 0.001
end


text:
@benchgroup "insertions" begin
  n = 10000
  @bench "ER Generation" g = SimpleGraph($n, 16 * $n)
end


text:
using Test
import AMRVW
const A = AMRVW
using LinearAlgebra
using Random

# test polynomials
p4 = [24.0, -50.0, 35.0, -10.0, 1.0]
p5 = [-120.0, 274.0, -225.0, 85.0, -15.0, 1.0]
p6 = [720.0, -1764.0, 1624.0, -735.0, 175.0, -21.0, 1.0]
pc = [0.0 + 1.0im, -1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 - 1.0im, 1.0 + 0.0im]

# could compare many ways
compare_eigenvalues(M1, M2) = norm(prod(eigvals(M1)) - prod(eigvals(M2))) <= sqrt(eps())

@testset "twisted-bulge-step-any-m" begin
    # bulge step should not change eigenvalues

    N = 10

    ## Real first
    T = Float64
    S = T;
    D = A.SparseDiagonal(T, 0)

    for m in 1:N-1
        @assert 0 < m < N

        ps = randperm(N)
        Ms = A.TwistedChain(A.random_rotator.(T, ps))
        Asc = A.random_rotator.(T, m:-1:1)

        RF = A.RFactorizationIdentity{T, S}()
        M = diagm(0 => ones(S, N+1))
        M1 = Ms * M
        A.bulge_step!(N,m, Ms, D, RF, Asc)
        M2 = Ms * (D * M)
        @test compare_eigenvalues(M1, M2)
    end

    ## complex
    S = Complex{T}
    D =  A.SparseDiagonal(S,N+1)

    for m in 1:N-1
        ps = randperm(N)
        Ms = A.TwistedChain(A.random_rotator.(S, ps))
        Asc = A.random_rotator.(S, m:-1:1)

        RF = A.RFactorizationIdentity{T, S}()
        M = diagm(0 => ones(S, N+1))
        M1 = Ms * M

        A.bulge_step!(N, m, Ms, D, RF, Asc)
        M2 = Ms * (D * M)
        @test compare_eigenvalues(M1, M2)
    end

end

@testset "eigenvalues twisted Q" begin

    # Basic algorithm
    ## real
    T = Float64
    S = Complex{T}
    Qs = A.random_rotator.(T,  [6, 10, 5, 20, 12, 14, 7, 3, 17, 19, 1, 13, 8, 18, 16, 4, 2, 9, 15, 11])
    QF = A.q_factorization(A.TwistedChain(Qs))
    state = A.QRFactorization(QF)
    e1 = eigvals(Matrix(state))
    @test norm(e1 - eigvals(state)) <= sqrt(eps(T))

    ## Complex
    Qs = A.random_rotator.(S,  [6, 10, 5, 20, 12, 14, 7, 3, 17, 19, 1, 13, 8, 18, 16, 4, 2, 9, 15, 11])
    QF = A.q_factorization(A.TwistedChain(Qs))
    state = A.QRFactorization(QF)
    e1 = eigvals(Matrix(state))
    es = eigvals(state)
    @test norm(e1 - es) <= sqrt(eps(T))

    # Different RF
    ## Real
    state = A.amrvw(p6)
    RF1 = state.RF
    RF2 = A.amrvw(rand(6), rand(6)).RF
    R = Matrix(RF1)[1:end-1, 1:end-1]
    RF3 = A.RFactorizationUpperTriangular(R)
    RF4 = A.RFactorizationUpperTriangular(LinearAlgebra.UpperTriangular(Matrix(R)))
    RF5 = A.RFactorizationUnitaryDiagonal(rand([-1.0, 1.0],6))
    RF6 = A.RFactorizationIdentity{T, S}()

    pv = [:right, :left, :right, :left]
    QF = A.QFactorizationTwisted(A.TwistedChain(state.QF.Q.x, pv))
    for RF in (RF1, RF2, RF3, RF4, RF5, RF6)
        state = A.QRFactorization(QF, RF)
        e1 = eigvals(Matrix(state))
        es = eigvals(state)
        @test norm(prod(e1) - prod(es)) < sqrt(eps(T))
    end

    ## Complex
    state = A.amrvw(pc)
    Q = state.QF.Q
    RF1 = state.RF
    RF2 = A.amrvw(rand(S,5), rand(S,5)).RF
    R = Matrix(RF1)#[1:end-1, 1:end-1]
    RF3 = A.RFactorizationUpperTriangular(R)
    RF4 = A.RFactorizationUpperTriangular(LinearAlgebra.UpperTriangular(Matrix(R)))
    ps = [complex(sincos(x)...) for x in rand(T,5)]
    RF5 = A.RFactorizationUnitaryDiagonal(ps)

    RF6 = A.RFactorizationIdentity{T, S}()

    pv = [:right, :left, :right]
    QF = A.QFactorizationTwisted(A.TwistedChain(state.QF.Q.x, pv))
    for RF in (RF1, RF2, RF3, RF4, RF5)
        state = A.QRFactorization(QF, RF)
        e1 = eigvals(Matrix(state))
        es = eigvals(state)
        @test norm(prod(e1) - prod(es)) < sqrt(eps(T))
    end

end