text:
export sum_index
function sum_index(a, b)
    r = a[1]
    @inbounds for i=1:length(a)
        r += a[i] + b[i]
    end
    r
end
export splat
@generated function splat{N}(a, ::Val{N})
    names = [Symbol("arg_$i") for i=1:N]
    tup = Expr(:tuple, names...)
    expr=Expr(:(=), tup, :a)
    expr=quote
        $expr
        $tup
    end
    expr
end

function rand_vector{X<:Vector}(A::Type{X}, N, T)
    rand(T, N)
end
function rand_vector{X<:FixedSizeArrays.FixedVector}(A::Type{X}, N, T)
    rand(A{N, T})
end
function rand_vector{X<:Tuple}(A::Type{X}, N, T)
    tuple(rand(T, N)...)
end
function rand_vector{X<:VecElement}(A::Type{X}, N, T)
    map(VecElement, tuple(rand(T, N)...))
end
function rand_vector{X<:SIMD.Vec}(A::Type{X}, N, T)
    SIMD.Vec(tuple(rand(T, N)...))
end
function rand_vector{X<:StaticArrays.SVector}(A::Type{X}, N, T)
    StaticArrays.rand(StaticArrays.SVector{N,T})
end

function vector_bench()
    constructors = [
        ("Vector", Vector),
        ("FSA", FixedSizeArrays.Vec),
        ("SIMD", SIMD.Vec),
        #("VecTuple", VecElement),
        ("Tuple", Tuple),
        ("StaticArrays", StaticArrays.SVector)
    ]

    element_types = [Float32, Float64, Int64, Int32]
    binary_funs = [sum_index, +, ./]
    mixed_funs = [+, /, *]
    unary_funs = [-, sum, prod]
    vector_of_vector_funs = [-, sum, prod, mean]

    # test if everything works (quicker than benchmarking)
    for (name, c) in constructors

        for N in (3, 10), T in element_types
            a = rand_vector(c, N, T)
            b = rand_vector(c, N, T)
            splat(a, Val{N}())
            for fun in unary_funs
                (fun)(a)
            end
            for fun in binary_funs
                (fun)(a, b)
            end
            for fun in mixed_funs
                x = rand(T)
                fun(a, x)
            end
            for fun in vector_of_vector_funs
                ET = typeof(rand_vector(c, N, T))
                a = ET[rand_vector(c, N, T) for _=1:1000]
                @assert eltype(a) == ET
                (fun)(a)
            end
        end
    end

    println("Okay gang, lets start benchmarkin'")
    for (name, c) in constructors
        suite[name] = BenchmarkGroup()
        for N in (3, 10), T in element_types
            a = rand_vector(c, N, T)
            b = rand_vector(c, N, T)

            #suite[name][splat, N, T] = @benchmarkable splat($a, $(Val{N}()))

            for fun in unary_funs
                suite[name][fun, N, T] = @benchmarkable $(fun)($a)
            end
            for fun in binary_funs
                suite[name][fun, N, T] = @benchmarkable $(fun)($a, $b)
            end
            for fun in mixed_funs
                x = rand(T)
                suite[name]["mixed", fun, N, T] = @benchmarkable $(fun)($a, $x)
            end
            for fun in vector_of_vector_funs
                ET = typeof(rand_vector(c, N, T))
                a = ET[rand_vector(c, N, T) for _=1:50]
                @assert eltype(a) == ET
                suite[name]["vecvec", fun, N, T] = @benchmarkable $(fun)($a)
            end
        end
    end
end


text:
module CounterfactualRegret

include(joinpath("Matrix", "MatrixGames.jl"))
export MatrixGame, MatrixPlayer, clear!, train_both!, train_one!
export evaluate

include(joinpath("Matrix", "CFR.jl"))
export SimpleIIGame, SimpleIIPlayer

include(joinpath("Extensive", "ExtensiveGames.jl"))
export Game, initialhist, isterminal, utility, player, next_hist, infokey, actions
export chance_action, chance_actions, other_player
export infokeytype, histtype

include(joinpath("Extensive", "CFR.jl"))
include(joinpath("Extensive", "CSMCCFR.jl"))
include(joinpath("Extensive", "ESMCCFR.jl"))
include(joinpath("Extensive", "DCFR.jl"))
export CFRSolver, CSCFRSolver, DCFRSolver, ESCFRSolver
export train!

include(joinpath("Extensive", "evaluation.jl"))
export FullEvaluate, MonteCarloEvaluate

include(joinpath("Extensive", "exploitability.jl"))
export ExploitabilitySolver, exploitability

include(joinpath("Extensive", "callback.jl"))

# Miscellaneous Games
include(joinpath("games", "Extensive2Matrix.jl"))
include(joinpath("games", "IIEMatrix.jl"))
include(joinpath("games", "Kuhn.jl"))

end # module

ext: jl

max_stars_repo_name: WhiffleFish/HelloCFR

max_issues_repo_name: WhiffleFish/HelloCFR

max_forks_repo_name: WhiffleFish/HelloCFR


text:
## Internal functions for getting and setting variable info
# Get info
function _variable_info(vref::InfOptVariableRef)::JuMP.VariableInfo
    return JuMP.owner_model(vref).vars[JuMP.index(vref)].info
end

# Set info for infinite variables
function _update_variable_info(vref::InfiniteVariableRef,
                               info::JuMP.VariableInfo)
    parameter_refs = JuMP.owner_model(vref).vars[JuMP.index(vref)].parameter_refs
    JuMP.owner_model(vref).vars[JuMP.index(vref)] = InfiniteVariable(info,
                                                                 parameter_refs)
    return
end

# Set info for point variables
function _update_variable_info(vref::PointVariableRef, info::JuMP.VariableInfo)
    infinite_variable_ref = JuMP.owner_model(vref).vars[JuMP.index(vref)].infinite_variable_ref
    parameter_values = JuMP.owner_model(vref).vars[JuMP.index(vref)].parameter_values
    JuMP.owner_model(vref).vars[JuMP.index(vref)] = PointVariable(info,
                                                          infinite_variable_ref,
                                                          parameter_values)
    return
end

# Set info for hold variables
function _update_variable_info(vref::HoldVariableRef, info::JuMP.VariableInfo)
    parameter_bounds = JuMP.owner_model(vref).vars[JuMP.index(vref)].parameter_bounds
    JuMP.owner_model(vref).vars[JuMP.index(vref)] = HoldVariable(info,
                                                               parameter_bounds)
    return
end

"""
    JuMP.has_lower_bound(vref::InfOptVariableRef)::Bool

Extend [`JuMP.has_lower_bound`](@ref JuMP.has_lower_bound(::JuMP.VariableRef)) to
return a `Bool` whether an `InfiniteOpt` variable has a lower bound.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, 0 <= vref))
julia> has_lower_bound(vref)
true
```
"""
JuMP.has_lower_bound(vref::InfOptVariableRef)::Bool = _variable_info(vref).has_lb

"""
    JuMP.lower_bound(vref::InfOptVariableRef)::Float64

Extend [`JuMP.lower_bound`](@ref JuMP.lower_bound(::JuMP.VariableRef)) to return the
lower bound of an `InfiniteOpt` variable. Errors if `vref` doesn't have a lower
bound.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, 0 <= vref))
julia> lower_bound(vref)
0.0
```
"""
function JuMP.lower_bound(vref::InfOptVariableRef)::Float64
    if !JuMP.has_lower_bound(vref)
        error("Variable $(vref) does not have a lower bound.")
    end
    return _variable_info(vref).lower_bound
end

# TODO modify to eliminate dictionaries if possible like SingleVariable in MOI
# Extend to return the index of the lower bound constraint associated with `vref`.
function JuMP._lower_bound_index(vref::InfOptVariableRef)::Int
    if !JuMP.has_lower_bound(vref)
        error("Variable $(vref) does not have a lower bound.")
    end
    return JuMP.owner_model(vref).var_to_lower_bound[JuMP.index(vref)]
end

# Extend to specify the index `cindex` of a lower bound constraint for `vref`.
function _set_lower_bound_index(vref::InfOptVariableRef, cindex::Int)
    JuMP.owner_model(vref).var_to_lower_bound[JuMP.index(vref)] = cindex
    return
end

"""
    JuMP.set_lower_bound(vref::InfOptVariableRef, lower::Number)

Extend [`JuMP.set_lower_bound`](@ref JuMP.set_lower_bound(::JuMP.VariableRef, ::Number))
to specify the lower bound of an `InfiniteOpt` variable `vref`. Errors if `vref`
is fixed.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, 0 <= vref))
julia> set_lower_bound(vref, -1)

julia> lower_bound(vref)
-1.0
```
"""
function JuMP.set_lower_bound(vref::InfOptVariableRef, lower::Number)
    newset = MOI.GreaterThan(convert(Float64, lower))
    if JuMP.has_lower_bound(vref)
        cindex = JuMP._lower_bound_index(vref)
        JuMP.owner_model(vref).constrs[cindex] = JuMP.ScalarConstraint(vref,
                                                                       newset)
        set_optimizer_model_ready(JuMP.owner_model(vref), false)
    else
        @assert !JuMP.is_fixed(vref) "$vref is fixed, cannot set lower bound."
        cref = JuMP.add_constraint(JuMP.owner_model(vref),
                                   JuMP.ScalarConstraint(vref, newset))
        _set_lower_bound_index(vref, JuMP.index(cref))
        JuMP.owner_model(vref).constr_in_var_info[JuMP.index(cref)] = true
    end
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(true, convert(Float64, lower),
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

"""
    JuMP.LowerBoundRef(vref::InfOptVariableRef)::GeneralConstraintRef

Extend [`JuMP.LowerBoundRef`](@ref JuMP.LowerBoundRef(::JuMP.VariableRef)) to extract
a constraint reference for the lower bound of `vref`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); vref = @hold_variable(m, 0 <= var))
julia> cref = LowerBoundRef(vref)
var ≥ 0.0
```
"""
function JuMP.LowerBoundRef(vref::InfOptVariableRef)::GeneralConstraintRef
    index = JuMP._lower_bound_index(vref)
    model = JuMP.owner_model(vref)
    if model.constrs[index].func isa InfiniteExpr
        return InfiniteConstraintRef(model, index,
                                     JuMP.shape(model.constrs[index]))
    else
        return FiniteConstraintRef(model, index,
                                   JuMP.shape(model.constrs[index]))
    end
end

"""
    JuMP.delete_lower_bound(vref::InfOptVariableRef)

Extend [`JuMP.delete_lower_bound`](@ref JuMP.delete_lower_bound(::JuMP.VariableRef))
to delete lower bound of `vref`. Errors if it doesn't have a lower bound.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, 0 <= vref))
julia> delete_lower_bound(vref)

julia> has_lower_bound(vref)
false
```
"""
function JuMP.delete_lower_bound(vref::InfOptVariableRef)
    JuMP.delete(JuMP.owner_model(vref), JuMP.LowerBoundRef(vref))
    delete!(JuMP.owner_model(vref).var_to_lower_bound, JuMP.index(vref))
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(false, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

"""
    JuMP.has_upper_bound(vref::InfOptVariableRef)::Bool

Extend [`JuMP.has_upper_bound`](@ref JuMP.has_upper_bound(::JuMP.VariableRef)) to
return a `Bool` whether an `InfiniteOpt` variable has an upper bound.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, 0 >= vref))
julia> has_upper_bound(vref)
true
```
"""
JuMP.has_upper_bound(vref::InfOptVariableRef)::Bool = _variable_info(vref).has_ub

"""
    JuMP.upper_bound(vref::InfOptVariableRef)::Float64

Extend [`JuMP.upper_bound`](@ref JuMP.upper_bound(::JuMP.VariableRef)) to return the
upper bound of an `InfiniteOpt` variable. Errors if `vref` doesn't have a upper bound.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, 0 >= vref))
julia> upper_bound(vref)
0.0
```
"""
function JuMP.upper_bound(vref::InfOptVariableRef)::Float64
    if !JuMP.has_upper_bound(vref)
        error("Variable $(vref) does not have a upper bound.")
    end
    return _variable_info(vref).upper_bound
end

# Extend to return the index of the upper bound constraint associated with `vref`.
function JuMP._upper_bound_index(vref::InfOptVariableRef)::Int
    if !JuMP.has_upper_bound(vref)
        error("Variable $(vref) does not have a upper bound.")
    end
    return JuMP.owner_model(vref).var_to_upper_bound[JuMP.index(vref)]
end

# Extend to specify the index `cindex` of a upper bound constraint for `vref`
function _set_upper_bound_index(vref::InfOptVariableRef, cindex::Int)
    JuMP.owner_model(vref).var_to_upper_bound[JuMP.index(vref)] = cindex
    return
end

"""
    JuMP.set_upper_bound(vref::InfOptVariableRef, upper::Number)

Extend [`JuMP.set_upper_bound`](@ref JuMP.set_upper_bound(::JuMP.VariableRef, ::Number))
to specify the upper bound of an `InfiniteOpt` variable `vref`. Errors if `vref`
is fixed.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, 0 >= vref))
julia> set_upper_bound(vref, 1)

julia> upper_bound(vref)
1.0
```
"""
function JuMP.set_upper_bound(vref::InfOptVariableRef, upper::Number)
    newset = MOI.LessThan(convert(Float64, upper))
    if JuMP.has_upper_bound(vref)
        cindex = JuMP._upper_bound_index(vref)
        JuMP.owner_model(vref).constrs[cindex] = JuMP.ScalarConstraint(vref,
                                                                       newset)
        set_optimizer_model_ready(JuMP.owner_model(vref), false)
    else
        @assert !JuMP.is_fixed(vref)
        cref = JuMP.add_constraint(JuMP.owner_model(vref),
                                   JuMP.ScalarConstraint(vref, newset))
        _set_upper_bound_index(vref, JuMP.index(cref))
        JuMP.owner_model(vref).constr_in_var_info[JuMP.index(cref)] = true
    end
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  true, convert(Float64, upper),
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

"""
    JuMP.UpperBoundRef(vref::InfOptVariableRef)::GeneralConstraintRef

Extend [`JuMP.UpperBoundRef`](@ref JuMP.UpperBoundRef(::JuMP.VariableRef)) to extract
a constraint reference for the upper bound of `vref`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); vref = @hold_variable(m, 1 >= var))
julia> cref = UpperBoundRef(vref)
var ≤ 1.0
```
"""
function JuMP.UpperBoundRef(vref::InfOptVariableRef)::GeneralConstraintRef
    index = JuMP._upper_bound_index(vref)
    model = JuMP.owner_model(vref)
    if model.constrs[index].func isa InfiniteExpr
        return InfiniteConstraintRef(model, index,
                                     JuMP.shape(model.constrs[index]))
    else
        return FiniteConstraintRef(model, index,
                                   JuMP.shape(model.constrs[index]))
    end
end

"""
    JuMP.delete_upper_bound(vref::InfOptVariableRef)

Extend [`JuMP.delete_upper_bound`](@ref JuMP.delete_upper_bound(::JuMP.VariableRef))
to delete the upper bound of `vref`. Errors if it doesn't have an upper bound.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, 0 >= vref))
julia> delete_upper_bound(vref)

julia> has_upper_bound(vref)
false
```
"""
function JuMP.delete_upper_bound(vref::InfOptVariableRef)
    JuMP.delete(JuMP.owner_model(vref), JuMP.UpperBoundRef(vref))
    delete!(JuMP.owner_model(vref).var_to_upper_bound, JuMP.index(vref))
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  false, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

"""
    JuMP.is_fixed(vref::InfOptVariableRef)::Bool

Extend [`JuMP.is_fixed`](@ref JuMP.is_fixed(::JuMP.VariableRef)) to return `Bool`
whether an `InfiniteOpt` variable is fixed.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref == 0))
julia> is_fixed(vref)
true
```
"""
JuMP.is_fixed(vref::InfOptVariableRef)::Bool = _variable_info(vref).has_fix

"""
    JuMP.fix_value(vref::InfOptVariableRef)::Float64

Extend [`JuMP.fix_value`](@ref JuMP.fix_value(::JuMP.VariableRef)) to return the fix
value of an `InfiniteOpt` variable. Errors if variable is not fixed.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref == 0))
julia> fix_value(vref)
0.0
```
"""
function JuMP.fix_value(vref::InfOptVariableRef)::Float64
    if !JuMP.is_fixed(vref)
        error("Variable $(vref) is not fixed.")
    end
    return _variable_info(vref).fixed_value
end

# Extend to return the index of the fix constraint associated with `vref`.
function JuMP._fix_index(vref::InfOptVariableRef)::Int
    if !JuMP.is_fixed(vref)
        error("Variable $(vref) is not fixed.")
    end
    return JuMP.owner_model(vref).var_to_fix[JuMP.index(vref)]
end

# Extend to set the index of the fix constraintassociated with `vref`.
function _set_fix_index(vref::InfOptVariableRef, cindex::Int)
    JuMP.owner_model(vref).var_to_fix[JuMP.index(vref)] = cindex
    return
end

"""
    JuMP.fix(vref::InfOptVariableRef, value::Number; force::Bool = false)

Extend [`JuMP.fix`](@ref JuMP.fix(::JuMP.VariableRef, ::Number)) to fix the value of
an `InfiniteOpt` variable. Errors if variable has a lower and/or an upper
bound(s) unless `force = true`.

**Examples**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref == 0); @hold_variable(m, vref2 <= 0))
julia> fix(vref, 3)

julia> fix_value(vref)
3.0

julia> fix(vref2, 2, force = true)

julia> fix_value(vref2)
2.0
```
"""
function JuMP.fix(vref::InfOptVariableRef, value::Number; force::Bool = false)
    new_set = MOI.EqualTo(convert(Float64, value))
    model = JuMP.owner_model(vref)
    if JuMP.is_fixed(vref)  # Update existing fixing constraint.
        cindex = JuMP._fix_index(vref)
        JuMP.owner_model(vref).constrs[cindex] = JuMP.ScalarConstraint(vref,
                                                                       new_set)
        set_optimizer_model_ready(JuMP.owner_model(vref), false)
    else  # Add a new fixing constraint.
        if  JuMP.has_upper_bound(vref) ||  JuMP.has_lower_bound(vref)
            if !force
                error("Unable to fix $(vref) to $(value) because it has " *
                      "existing variable bounds. Consider calling " *
                      "`JuMP.fix(variable, value; force=true)` which will " *
                      "delete existing bounds before fixing the variable.")
            end
            if  JuMP.has_upper_bound(vref)
                 JuMP.delete_upper_bound(vref)
            end
            if  JuMP.has_lower_bound(vref)
                 JuMP.delete_lower_bound(vref)
            end
        end
        cref = JuMP.add_constraint(model, JuMP.ScalarConstraint(vref, new_set))
        _set_fix_index(vref, JuMP.index(cref))
        JuMP.owner_model(vref).constr_in_var_info[JuMP.index(cref)] = true
    end
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(false, info.lower_bound,
                                                  false, info.upper_bound,
                                                  true, convert(Float64, value),
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

"""
    JuMP.FixRef(vref::InfOptVariableRef)::GeneralConstraintRef

Extend [`JuMP.FixRef`](@ref JuMP.FixRef(::JuMP.VariableRef)) to return the constraint
reference of the fix constraint associated with `vref`. Errors `vref` is not
fixed.

**Examples**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); vref = @hold_variable(m, var == 1))
julia> cref = FixRef(vref)
var = 1.0
```
"""
function JuMP.FixRef(vref::InfOptVariableRef)::GeneralConstraintRef
    index = JuMP._fix_index(vref)
    model = JuMP.owner_model(vref)
    if model.constrs[index].func isa InfiniteExpr
        return InfiniteConstraintRef(model, index,
                                     JuMP.shape(model.constrs[index]))
    else
        return FiniteConstraintRef(model, index,
                                   JuMP.shape(model.constrs[index]))
    end
end

"""
    JuMP.unfix(vref::InfOptVariableRef)

Extend [`JuMP.unfix`](@ref JuMP.unfix(::JuMP.VariableRef)) to unfix `vref`. Errors
if it is not fixed.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref == 0))
julia> unfix(vref)

julia> is_fixed(vref)
false
```
"""
function JuMP.unfix(vref::InfOptVariableRef)
    JuMP.delete(JuMP.owner_model(vref), JuMP.FixRef(vref))
    delete!(JuMP.owner_model(vref).var_to_fix, JuMP.index(vref))
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  false, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

"""
    JuMP.start_value(vref::InfOptVariableRef)::Union{Nothing, Float64}

Extend [`JuMP.start_value`](@ref JuMP.start_value(::JuMP.VariableRef)) to return
starting value of `InfiniteOpt` variable if it has one. Returns `nothing` otherwise.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref, start = 0))
julia> start_value(vref)
0.0
```
"""
function JuMP.start_value(vref::InfOptVariableRef)::Union{Nothing, Float64}
    return _variable_info(vref).start
end

"""
    JuMP.set_start_value(vref::InfOptVariableRef, value::Number)

Extend [`JuMP.set_start_value`](@ref JuMP.set_start_value(::JuMP.VariableRef, ::Number))
to specify the start value of `InfiniteOpt` variables.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref, start = 0))
julia> set_start_value(vref, 1)

julia> start_value(vref)
1.0
```
"""
function JuMP.set_start_value(vref::InfOptVariableRef, value::Number)
    info = _variable_info(vref)
    set_optimizer_model_ready(JuMP.owner_model(vref), false)
    _update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           true, Float64(value),
                                           info.binary, info.integer))
    return
end

"""
    JuMP.is_binary(vref::InfOptVariableRef)::Bool

Extend [`JuMP.is_binary`](@ref JuMP.is_binary(::JuMP.VariableRef)) to return `Bool`
whether an `InfiniteOpt` variable is binary.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref, Bin))
julia> is_binary(vref)
true
```
"""
JuMP.is_binary(vref::InfOptVariableRef)::Bool = _variable_info(vref).binary

# Extend to return the index of the binary constraint associated with `vref`.
function JuMP._binary_index(vref::InfOptVariableRef)::Int
    if !JuMP.is_binary(vref)
        error("Variable $(vref) is not binary.")
    end
    return JuMP.owner_model(vref).var_to_zero_one[JuMP.index(vref)]
end

# Extend to specify the index of the binary constraint associated with `vref`.
function _set_binary_index(vref::InfOptVariableRef, cindex::Int)
    JuMP.owner_model(vref).var_to_zero_one[JuMP.index(vref)] = cindex
    return
end

"""
    JuMP.set_binary(vref::InfOptVariableRef)

Extend [`JuMP.set_binary`](@ref JuMP.set_binary(::JuMP.VariableRef)) to specify an
`InfiniteOpt` variable as a binary variable. Errors if `vref` is an integer
variable.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref, Bin))
julia> set_binary(vref)

julia> is_binary(vref)
true
```
"""
function JuMP.set_binary(vref::InfOptVariableRef)
    if JuMP.is_binary(vref)
        return
    elseif JuMP.is_integer(vref)
        error("Cannot set the variable_ref $(vref) to binary as it " *
              "is already integer.")
    end
    cref = JuMP.add_constraint(JuMP.owner_model(vref), JuMP.ScalarConstraint(vref, MOI.ZeroOne()))
    _set_binary_index(vref, JuMP.index(cref))
    JuMP.owner_model(vref).constr_in_var_info[JuMP.index(cref)] = true
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  true, info.integer))
    return
end

"""
    JuMP.BinaryRef(vref::InfOptVariableRef)::GeneralConstraintRef

Extend [`JuMP.BinaryRef`](@ref JuMP.BinaryRef(::JuMP.VariableRef)) to return a
constraint reference to the constraint constrainting `vref` to be binary. Errors
if one does not exist.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); vref = @hold_variable(m, var, Bin))
julia> cref = BinaryRef(vref)
var binary
```
"""
function JuMP.BinaryRef(vref::InfOptVariableRef)::GeneralConstraintRef
    index = JuMP._binary_index(vref)
    model = JuMP.owner_model(vref)
    if model.constrs[index].func isa InfiniteExpr
        return InfiniteConstraintRef(model, index,
                                     JuMP.shape(model.constrs[index]))
    else
        return FiniteConstraintRef(model, index,
                                   JuMP.shape(model.constrs[index]))
    end
end

"""
    JuMP.unset_binary(vref::InfOptVariableRef)

Extend [`JuMP.unset_binary`](@ref JuMP.unset_binary(::JuMP.VariableRef)) to unset
`vref` as a binary variable. Errors if it is not binary.

```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref, Bin))
julia> unset_binary(vref)

julia> is_binary(vref)
false
```
"""
function JuMP.unset_binary(vref::InfOptVariableRef)
    JuMP.delete(JuMP.owner_model(vref), JuMP.BinaryRef(vref))
    delete!(JuMP.owner_model(vref).var_to_zero_one, JuMP.index(vref))
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  false, info.integer))
    return
end

"""
    JuMP.is_integer(vref::InfOptVariableRef)::Bool

Extend [`JuMP.is_integer`](@ref JuMP.is_integer(::JuMP.VariableRef)) to return `Bool`
whether an `InfiniteOpt` variable is integer.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref, Int))
julia> is_integer(vref)
true
```
"""
JuMP.is_integer(vref::InfOptVariableRef)::Bool = _variable_info(vref).integer

# Extend to return the index of the integer constraintassociated with `vref`.
function JuMP._integer_index(vref::InfOptVariableRef)::Int
    if !JuMP.is_integer(vref)
        error("Variable $(vref) is not an integer.")
    end
    return JuMP.owner_model(vref).var_to_integrality[JuMP.index(vref)]
end

# Extend to specify the index of the integer constraint associated with `vref`.
function _set_integer_index(vref::InfOptVariableRef, cindex::Int)
    JuMP.owner_model(vref).var_to_integrality[JuMP.index(vref)] = cindex
    return
end

"""
    JuMP.set_integer(vref::InfOptVariableRef)

Extend [`JuMP.set_integer`](@ref JuMP.set_integer(::JuMP.VariableRef)) to specify an
`InfiniteOpt` variable as a integer variable. Errors if `vref` is an binary
variable.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref, Int))
julia> set_integer(vref)

julia> is_integer(vref)
true
```
"""
function JuMP.set_integer(vref::InfOptVariableRef)
    if JuMP.is_integer(vref)
        return
    elseif JuMP.is_binary(vref)
        error("Cannot set the variable_ref $(vref) to integer as it " *
              "is already binary.")
    end
    cref = JuMP.add_constraint(JuMP.owner_model(vref), JuMP.ScalarConstraint(vref,
                               MOI.Integer()))
    _set_integer_index(vref, JuMP.index(cref))
    JuMP.owner_model(vref).constr_in_var_info[JuMP.index(cref)] = true
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, true))
    return
end

"""
    JuMP.IntegerRef(vref::InfOptVariableRef)::GeneralConstraintRef

Extend [`JuMP.IntegerRef`](@ref JuMP.IntegerRef(::JuMP.VariableRef)) to return a
constraint reference to the constraint constrainting `vref` to be integer.
Errors if one does not exist.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); vref = @hold_variable(m, var, Int))
julia> cref = IntegerRef(vref)
var integer
```
"""
function JuMP.IntegerRef(vref::InfOptVariableRef)::GeneralConstraintRef
    index = JuMP._integer_index(vref)
    model = JuMP.owner_model(vref)
    if model.constrs[index].func isa InfiniteExpr
        return InfiniteConstraintRef(model, index,
                                     JuMP.shape(model.constrs[index]))
    else
        return FiniteConstraintRef(model, index,
                                   JuMP.shape(model.constrs[index]))
    end
end

"""
    JuMP.unset_integer(vref::InfOptVariableRef)

Extend [`JuMP.unset_integer`](@ref JuMP.unset_integer(::JuMP.VariableRef)) to unset
`vref` as an integer variable. Errors if it is not an integer variable.

```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref, Int))
julia> unset_integer(vref)

julia> is_integer(vref)
false
```
"""
function JuMP.unset_integer(vref::InfOptVariableRef)
    JuMP.delete(JuMP.owner_model(vref), JuMP.IntegerRef(vref))
    delete!(JuMP.owner_model(vref).var_to_integrality, JuMP.index(vref))
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, false))
    return
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
name = "H3"
version = v"3.4.2"

# function url2hash(url)
#     path = download(url)
#     open(io-> bytes2hex(BinaryProvider.sha256(io)), path)
# end
# url2hash("https://github.com/uber/h3/archive/v3.4.2.zip") |> println

sources = [
    "https://github.com/uber/h3/archive/v$version.zip" => "c5d024af8f7a852349ffce69fe33b456f96d7c940c4ffbeb872f98318b21e03c",
]

# Bash recipe for building across all platforms
script = """
cd \$WORKSPACE/srcdir/h3-$version/
cat <<EOF > CMakeLists.txt.patch
diff -uNr h3-3.4.2-original/CMakeLists.txt h3-3.4.2/CMakeLists.txt
--- h3-3.4.2-original/CMakeLists.txt    2019-02-23 18:30:30.000000000 +0900
+++ h3-3.4.2/CMakeLists.txt 2019-02-23 17:58:02.000000000 +0900
@@ -53,6 +53,8 @@

 project(h3 LANGUAGES C VERSION \\\${H3_VERSION})

+set(CMAKE_C_FLAGS "-std=c99 \\\${CMAKE_C_FLAGS}")
+
 set(H3_COMPILE_FLAGS "")
 set(H3_LINK_FLAGS "")
 if(NOT WIN32)
@@ -495,21 +497,6 @@
     add_h3_executable(mkRandGeo src/apps/testapps/mkRandGeo.c \\\${APP_SOURCE_FILES})
     add_h3_executable(mkRandGeoBoundary src/apps/testapps/mkRandGeoBoundary.c \\\${APP_SOURCE_FILES})

-    # Benchmarks
-    add_custom_target(benchmarks)
-
-    macro(add_h3_benchmark name srcfile)
-        add_h3_executable(\\\${name} \\\${srcfile} \\\${APP_SOURCE_FILES})
-        add_custom_target(bench_\\\${name} COMMAND \\\${TEST_WRAPPER} \\\$<TARGET_FILE:\\\${name}>)
-        add_dependencies(benchmarks bench_\\\${name})
-    endmacro()
-
-    add_h3_benchmark(benchmarkH3Api src/apps/benchmarks/benchmarkH3Api.c)
-    add_h3_benchmark(benchmarkKRing src/apps/benchmarks/benchmarkKRing.c)
-    add_h3_benchmark(benchmarkH3Line src/apps/benchmarks/benchmarkH3Line.c)
-    add_h3_benchmark(benchmarkH3SetToLinkedGeo src/apps/benchmarks/benchmarkH3SetToLinkedGeo.c)
-    add_h3_benchmark(benchmarkPolyfill src/apps/benchmarks/benchmarkPolyfill.c)
-    add_h3_benchmark(benchmarkPolygon src/apps/benchmarks/benchmarkPolygon.c)
 endif()

 # Installation (https://github.com/forexample/package-example)
EOF
patch -p1 -i CMakeLists.txt.patch
mkdir build
cd build
cmake -DBUILD_SHARED_LIBS=1 -DCMAKE_INSTALL_PREFIX=\$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/\$target/\$target.toolchain ..
make
make install
rm -rf \$WORKSPACE/destdir/bin/{g,h,k}* \$WORKSPACE/destdir/lib/cmake \$WORKSPACE/destdir/logs
ls \$WORKSPACE/destdir/lib
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libh3", :libh3),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
# add the julia packages used
Pkg.add("JSON")
Pkg.add("PyCall")
Pkg.add("ArgParse")
Pkg.add("StatsBase")

max_stars_repo_name: Condimentum/spacetyckiting

max_issues_repo_name: Condimentum/spacetyckiting

max_forks_repo_name: Condimentum/spacetyckiting


text:
#

# wrap in individual modules to avoid name conflicts.
module TestAutomaticDifferentiation
include("automatic_differentiation_test.jl")
end

module TestDual
include("dual_n.jl")
end


text:
srand(5)
include("setup.jl")
simid = 2
simname = "dwj"
hyper = false
data = readcsv("output/volatility/dwj.csv", header=true)[1]
Xd = Vector{Float64}(log.(data[:, 2]/data[1, 2]))
quc = 90
qu = quc/100

td = map(d -> Dates.value(Dates.Date(d)), data[:,1])/1
dates = map(d -> Dates.value(Dates.Date(d)), data[:,1])/365
dates += -dates[1] + 1971 + (Dates.value(Date("1971-07-02")) - Dates.value(Date("1971-01-01")))/365

td -= td[1]
dt = td[2] - td[1]

#assert(all(diff(td) .== 7))
X = data[:,2][:]
sub = 10
lim = (740, 1060)
display(plot(dates, X,  xaxis = ("y", (dates[1],dates[end]), 1972:1974), linewidth = 0.4, ylim = lim, size = (400, 400/1.4), color=colorant"#0044FF",  label=L"X"))
savefig("output/volatility/fig$(simname)pathX.pdf")
lim = (-0.2, 0.2)
display(plot(dates, Xd,  xaxis = ("y", (dates[1],dates[end]), 1972:1974), linewidth = 0.4, ylim = lim, size = (400, 400/1.4), color=colorant"#0044FF",  label=L"Z"))
savefig("output/volatility/fig$(simname)pathZ.pdf")
lim = (-0.1, 0.1)
display(plot(dates[1:end-1], diff(X)./X[1:end-1], xaxis = ("y", (dates[1],dates[end]), 1972:1974), linewidth = 0.4, ylim = lim, size = (400, 400/1.4), color=colorant"#0044FF",  label=L"Y"))
savefig("output/volatility/fig$(simname)pathY.pdf")

N = length(Xd)

if simid == 1
    a1 = 0.0
    a = 5.0
    aζ = a

    M = 6
    hyper = true

    σa = 10.0
elseif simid == 2
    a1 = 0.0
    a = 5.0
    aζ = a

    M = 12
    hyper = true

    σa = 10.0
elseif simid == 3
    a1 = 0.0
    a = 5.0
    aζ = a

    M = 6
    hyper = true

    σa = 12.0
    Xd = X
elseif simid == 4
    a1 = 0.0
    a = 5.0
    aζ = a

    M = 12
    hyper = true

    σa = 12.0
    Xd = X
else
    error("")
end

pri = InverseGamma(0.1, 0.1)
priind= InverseGamma(0.001,0.001)
pri2 = InverseGamma(0.3, 0.3)
#pri2 = InverseGamma(2.5, 2.5) #try

if hyper
    iterations = 200000
else
    iterations = 200000
end
#
#M = (length(td)-1) ÷ H # observed intervals per bin
#r = length(td)-1 - H*M
#assert(H*M == length(td)-1)
H = (N-1)÷M
r = length(td)-1 - H*M

# compute independent posterior

post = Vector(H)
Z = Vector(H)
x = Vector(H)
for i in 1:H
    if i == H
        ii = 1+(i-1)*M:N
    else
        ii = 1+(i-1)*M:(i)*M+1
    end
    tx = td[ii]

    x[i] = diff(Xd[ii])./sqrt.(diff(tx)) # normalize obs
    Z[i] = sum(diff(Xd[ii]).^2)
    post[i] = ConjugatePriors.posterior((0., priind), Normal, x[i])
 #   post[i] = InverseGamma(shape(priind) + (M + (i==H)*r)/2, scale(priind) + Z[i]/2/dt)

end

if simid <= 2
    lim = (-0.01, .042)
else
    lim = (-2, 22)
end
bands = map(p->sqrt(quantile(p, (1 - qu)/2)), post), map(p->sqrt(quantile(p, 1 - (1 - qu)/2)), post)

plot(piecewise( Bridge.SamplePathBand(dates[1:M:end], bands...) ),
    ylim = lim, xaxis = ("y", (dates[1],dates[end]), 1972:1974), size = (400, 400/1.4), color=colorant"#0044FF", label = "marginal $quc\\% cred.~band",);
display(Plots.current())
savefig("output/volatility/fig$(simname)postind$simid.pdf")

v = zeros(H)
ζ = zeros(H)

# initialize MCMC with a draw from the prior
if a1 > 0.01
    v[1] = rand(InverseGamma(a1, a1))
else
    v[1] = 1.0
end
for i in 2:H
    ζ[i] = rand(InverseGamma(aζ, inv(v[i-1]/aζ)))
    v[i] = rand(InverseGamma(a, inv(ζ[i]/a)))
end

subsampl = 1
burnin = 1000
subinds = (1 + burnin):subsampl:iterations
samples = zeros(H, length(subinds))

state = Bridge.mcstart(v)
si = 1
printiter = 10
as = []
acc = 0
for iter in 1:iterations
    # move: parameters
    v[1] = rand(InverseGamma(a1 + aζ + M/2, a1 + aζ/ζ[2] + Z[1]/2/dt))
    for h in 2:H-1
        v[h] = rand(InverseGamma(aζ + a + M/2, (a/ζ[h] + aζ/ζ[h+1] + Z[h]/2/dt)))
    end
    v[H] = rand(InverseGamma(aζ + (M+r)/2, a/ζ[H] + Z[H]/2/dt))

    # move: latent
    if iter in subinds
        samples[:, si] = v
        si += 1
    end
    for h in 2:H
        ζ[h] = rand(InverseGamma(aζ + a, (aζ/v[h-1] + a/v[h])))
    end

    state = Bridge.mcnext(state, v)

    # move: hyper parameters
    if hyper
        a˚ = a + σa*randn()
        while a˚ < eps()
            a˚ = a + σa*randn()
        end

        lq = logpdf(pri2, a)
        lq += (2*(H-1))*(a*log(a) - lgamma(a))
        s = sum(log(v[k-1]*v[k]*ζ[k]*ζ[k]) + (1/v[k-1] + 1/v[k])/ζ[k]  for k = 2:H)
        lq += -a*s

        lq˚ = logpdf(pri2, a˚)
        lq˚ += (2*(H-1))*(a˚*log(a˚) - lgamma(a˚))
        lq˚ += -a˚*s
        mod(iter, printiter) == 0 && print("$iter \t α ", a˚, (lq˚,lq))
        if rand() < exp(lq˚ - lq)*cdf(Normal(0, σa), a)/cdf(Normal(0, σa), a˚)
            acc = acc + 1
            aζ = a = a˚
            mod(iter, printiter) == 0 && print("✓")
        end
        mod(iter, printiter) == 0 && println()
        push!(as, a)
    end

end
hyper && println("acc prob hyper ", acc/iterations)
writecsv("output/volatility/samples$simname$simid.csv", samples')
if hyper
    writecsv("output/volatility/samplesas$simname$simid.csv", as)
end
if simid <= 2
    lim = (-0.01, .027)
else
    lim = (-2, 22)
end
quc = 90
qu = quc/100
bands = mapslices(v->sqrt(quantile(v, (1-qu)/2)), samples, (2,))[:],  mapslices(v->sqrt(quantile(v, 1 - (1 - qu)/2)), samples, (2,))[:]

plot(piecewise( Bridge.SamplePathBand(dates[1:M:end], bands...) ),
    ylim = lim, xaxis = ("y", (dates[1],dates[end]), 1972:1974), size = (400, 400/1.4),color=colorant"#0044FF", label = "marginal $quc\\% cred.~band",);
display(Plots.current())
savefig("output/volatility/fig$(simname)post$simid.pdf")

open(joinpath("output", "volatility", "$(simname)$(simid)info.txt"), "w") do f
    println(f, "α1 ", a1)
    !hyper && println(f, "α ", a)
    println(f, "M ", M)
    println(f, "iterations ", iterations)
    println(f, "samples at ", subinds)
    println(f, "$quc % credible bands")

    println(f, "hyper ", hyper)
    if hyper
        println(f, "acc ", round(acc/iterations, 3))
        println(f,  "σa ", σa)
    end
end

# TODO
# use quantiles
# add samplepath


text:
#!/usr/bin/env julia
# vim:set ft=julia ts=4 sw=4 sts=4 autoindent:

# Performance testing for the parsing modules.
#
# Author:   Pontus Stenetorp    <pontus stenetorp se>
# Version:  2014-05-09

include("../hdr.jl")

require("conllx.jl")
require("dep.jl")
require("hybrid.jl")
require("parse.jl")

using CoNLLX
using DepGraph
using Hybrid
using Parse

open(talpath, "r") do data_f
    sents = collect(Sentence, conllxparse(data_f, useproj=true))
    trainitr = train(sents, epochs=1)
    for pass in 1:2
        before = time()
        next(trainitr)
        duration = time() - before
        throughput = length(sents) / duration
        print("Unlabelled training throughput ")
        @printf("(pass: %d): %.1f ", pass, throughput)
        println("sentence(s) per second")
    end
end

max_stars_repo_name: ninjin/allen

max_issues_repo_name: ninjin/allen

max_forks_repo_name: ninjin/allen


text:
module Kernels
using JuliaPBF.StructPBF
export
    Wendland, WendlandGrad, WendlandGradAbs,
    Poly6, Poly6Grad, Poly6GradAbs

    const oneOverPi = 1 / π

    function Wendland(absDelPos::Float64, inAnsDataStruct::AnalysisDataStruct)
        oneOverH= 1/inAnsDataStruct.kernelRadius
        qq      = absDelPos * oneOverH
        alpha   = 1.75*oneOverH*oneOverH*oneOverPi # For 2-Dim
        return alpha * (1- qq*0.5)^4 * (1 + 2.0*qq)
    end

    function WendlandGrad(delPos::Vec2, inAnsDataStruct::AnalysisDataStruct)
        oneOverH= 1/inAnsDataStruct.kernelRadius
        absDelPos = sqrt(abs2(delPos))
        qq      = absDelPos * oneOverH
        alpha   = 1.75*oneOverH*oneOverH*oneOverPi # For 2-Dim
        if (absDelPos < 0.000000000000001)
            return 0.0
        end
        if (2.0 < qq)
            return 0.0
        end
        return (-5 * alpha * (1- qq*0.5)^3 * qq / absDelPos) * delPos
    end

    function WendlandGradAbs(delPos::Vec2, inAnsDataStruct::AnalysisDataStruct)
        oneOverH= 1/inAnsDataStruct.kernelRadius
        absDelPos = sqrt(abs2(delPos))
        qq      = absDelPos * oneOverH
        alpha   = 1.75*oneOverH*oneOverH*oneOverPi # For 2-Dim
        if (absDelPos < 0.000000000000001)
            return 0.0
        end
        if (2.0 < qq)
            return 0.0
        end
        return (-5 * alpha * (1- qq*0.5)^3 * qq)
    end

    function Poly6(absDelPos::Float64, inAnsDataStruct::AnalysisDataStruct)
        oneOverH  = 1/inAnsDataStruct.kernelRadius
        qq        = absDelPos * oneOverH
        alpha     = 4*(oneOverH^2)*oneOverPi # For 2-Dim
        return alpha * (1- qq*qq)^3
    end

    function Poly6Grad(delPos::Vec2, inAnsDataStruct::AnalysisDataStruct)
        oneOverH  = 1/inAnsDataStruct.kernelRadius
        absDelPos = sqrt(abs2(delPos))
        qq        = absDelPos * oneOverH
        alphaGrad = 4*(oneOverH^4)*oneOverPi # For 2-Dim
        if (absDelPos < 0.000000000000001)
            return 0.0
        end
        if (1.0 < qq)
            return 0.0
        end
        return - 6 * alphaGrad * (1- qq*qq)^2 * delPos
    end
end

max_stars_repo_name: Hipgineer/JuliaPBF

max_issues_repo_name: Hipgineer/JuliaPBF

max_forks_repo_name: Hipgineer/JuliaPBF


text:
#TODO: Inflation-index the wealth distributions in figuresByCategory

##Graph types and constants
const COLORS = ["DodgerBlue", "OrangeRed", "Green",
  "Brown", "Black", "BlueViolet"]

const PlotContainer = Union{Plot,Gadfly.Compose.Context,Vector{Plot}} #holds graphs
const ASSET_YEARS = [2015, 2010, 2005, 2000, 1997]
const ASSET_YEARS_FULL = collect(1997:2015)
const WEALTH_FOCAL_FIELD = :adjnetassets

const RETURN_FOCAL_FIELD = :return
const RETURN_YEARS = ASSET_YEARS
const RETURN_YEARS_FULL = collect(2000:2015)

const CATEGORY_FIGURE_FIELD = :categorynameconsolidated
const BIN_ON_FIELDS = [:fisyr, CATEGORY_FIGURE_FIELD]

#NOTE: Create a new more consolidated category field

getNonmissingCategories(data::NCCSData; categoryField::Symbol=CATEGORY_FIGURE_FIELD) =
  unique(data.df[(!ismissing).(data.df[categoryField]), categoryField])

#const FIGURES_MIN_ASSETS = 1_500_000.

#contains the functions for writing the grpah by year
function wealthFiguresByYear(data::NCCSData; assetYears::Vector{Int} = ASSET_YEARS,
    outputPath::String = GRAPH_PATH, focalField::Symbol=WEALTH_FOCAL_FIELD,
    categoryFigureField::Symbol = CATEGORY_FIGURE_FIELD,
    categoriesUsed::Vector = getNonmissingCategories(data, categoryField = categoryFigureField))

  local plotNames::Vector{String} = Vector{String}()
  local plots::Vector{PlotContainer} = Vector{PlotContainer}()
  local subdf::SubDataFrame
  #local presubdf::SubDataFrame =  view(data.df, data.df[focalField] .≥ figuresMinAssets,:)

  for y::Int ∈ assetYears #main wealth distribution
    subdf = view(data.df, data.df[:fisyr] .== y,:) #make a view for performance reasons

    push!(plotNames, "Wealth_Distribution_of_Non-Profits_($y)-dodge")
    push!(plots, plot(subdf, x=focalField, color=:charitytype,
      Guide.title("Wealth Distribution of Non-Profits ($y)-dodge"),
      Guide.xlabel("Wealth"),Guide.ylabel("Density"),
      Scale.x_log10,
      Geom.histogram(bincount=50, density=false, position=:dodge),
      Coord.cartesian(ymin=0.0, ymax = 2500.0, xmin=5.0)))

    push!(plotNames, "Wealth_Distribution_of_Non-Profits_($y)")
    push!(plots, plot(subdf,
      x=focalField, color=categoryFigureField,
      Guide.title("Wealth Distribution of Non-Profits ($y)"),
      Guide.xlabel("Wealth"), Guide.ylabel("Amount"),  Scale.x_log10,
      Geom.line, Stat.histogram(bincount=50, density=false),
      Coord.cartesian(ymin=0.0, ymax = 3000.0, xmin=5.0)
      ))

  #=  push!(plotNames, "Wealth_Distribution_of_Private_Charities_($y)")
    push!(plots, plot(subdf[subdf[:charitytype] .== :pc,:],
      x=focalField, color=categoryFigureField,
      Guide.title("Wealth Distribution of Non-Profits ($y)"),
      Guide.xlabel("Wealth"), Guide.ylabel("Amount"),
      Scale.x_log10, Geom.line,
      Stat.histogram(bincount=50, density=false),
      Coord.cartesian(ymin=0.0, ymax = 2500.0, xmin=5.0)
      ))

    push!(plotNames, "Wealth_Distribution_of_Private_Foundations_($y)")
    push!(plots, plot(subdf[subdf[:charitytype] .== :pf,:],
      x=focalField, color=categoryFigureField,
      Guide.title("Wealth Distribution of Non-Profits ($y)"),
      Guide.xlabel("Wealth"), Guide.ylabel("Amount"),
      Scale.x_log10, Geom.line,
      Stat.histogram(bincount=50, density=false),
      Coord.cartesian(ymin=0.0, ymax = 1000.0, xmin=5.0)
      ))=#

    #=push!(plotNames, "Log_Wealth_Distribution_of_Non-Profits_($y)")
    push!(plots, plot(subdf,
    x=Symbol('l',focalfield), color=:charitytype,
      Guide.title("Log Wealth Distribution of Non-Profits ($y)"),
      Guide.xlabel("Wealth"),Guide.ylabel("Density"),
      Scale.x_log10,
      Geom.histogram(bincount=50, density=true)))=#
  end

  for i ∈ 1:length(plots) #write the graphs
    try
      draw(PDF("$outputPath\\$(plotNames[i]).pdf", 9inch, 7inch), plots[i])
    catch err
      @error "ERROR writing $(plotNames[i]). Error message $err"
      println("Size: $(size(subdf))")
    end
    #println("Graph $(plotNames[i]) written.")
  end

  println("Cross-sectional Wealth Graphs by Year Drawn")
  return nothing
end

#contains the functions for writing the grpah by year
function returnFiguresByYear(data::NCCSData; returnYears::Vector{Int} = RETURN_YEARS,
    outputPath::String = GRAPH_PATH, focalField::Symbol=RETURN_FOCAL_FIELD,
    categoryFigureField::Symbol = CATEGORY_FIGURE_FIELD,
    categoriesUsed::Vector = getNonmissingCategories(data, categoryField = categoryFigureField))

  local plotNames::Vector{String} = Vector{String}()
  local plots::Vector{PlotContainer} = Vector{PlotContainer}()
  local subdf::SubDataFrame
  local presubdf::SubDataFrame = view(data.df, (!ismissing).(data.df[focalField]),:)
  presubdf = view(presubdf, (!ismissing).(presubdf[categoryFigureField]),:)

  for y::Int ∈ returnYears #main wealth distribution
    subdf = view(presubdf, presubdf[:fisyr] .== y,:) #make a view for performance reasons

    #println("\n####\nYear: $y; rows: $(size(subdf,1)),
    #  rowspf: rows: $(sum(subdf[:charitytype] .== :pf)),
    #  rowspc: rows: $(sum(subdf[:charitytype] .== :pc)))")

    push!(plotNames, "Return_Distribution_of_Non-Profits_($y)_by_Type")
    push!(plots, plot(subdf, x=focalField, color=:charitytype,
      Guide.title("Return Distribution of Non-Profits ($y) by Type"),
      Guide.xlabel("Return"),Guide.ylabel("Density"),
      Geom.histogram(bincount=50, density=true#=, position=:dodge=#),
      Coord.cartesian(xmin=-0.5, xmax=0.5)))

    push!(plotNames, "Return_Distribution_of_Non-Profits_($y)_by_Category")
    push!(plots, plot(subdf,
      x=focalField, color=categoryFigureField,
      Guide.title("Return_Distribution_of_Non-Profits_($y)_by_Category"),
      Guide.xlabel("Return"), Guide.ylabel("Amount"),
      Geom.line, Stat.histogram(bincount=50, density=false),
      Coord.cartesian(ymin=0.0, ymax = 1000.0, xmin=-0.5, xmax=0.5)
      ))

    push!(plotNames, "Return_Distribution_of_Private_Charities_($y)")
    push!(plots, plot(subdf[subdf[:charitytype] .== :pc,:],
      x=focalField, color=categoryFigureField,
      Guide.title("Return Distribution of Non-Profits ($y)"),
      Guide.xlabel("Return"), Guide.ylabel("Amount"),
      Geom.line,
      Stat.histogram(bincount=50, density=false),
      Coord.cartesian(ymin=0.0, ymax = 1000.0, xmin=-0.5, xmax=0.5)
      ))

    push!(plotNames, "Return_Distribution_of_Private_Foundations_($y)")
    push!(plots, plot(subdf[subdf[:charitytype] .== :pf,:],
      x=focalField, color=categoryFigureField,
      Guide.title("Return Distribution of Non-Profits ($y)"),
      Guide.xlabel("Return"), Guide.ylabel("Amount"),
      Geom.line,
      Stat.histogram(bincount=50, density=true),
      Coord.cartesian(ymin=0.0, ymax = 1000.0, xmin=-0.5, xmax=0.5)
      ))

    #=push!(plotNames, "Log_Wealth_Distribution_of_Non-Profits_($y)")
    push!(plots, plot(subdf,
    x=Symbol('l',focalfield), color=:charitytype,
      Guide.title("Log Wealth Distribution of Non-Profits ($y)"),
      Guide.xlabel("Wealth"),Guide.ylabel("Density"),
      Scale.x_log10,
      Geom.histogram(bincount=50, density=true)))=#
  end

  for i ∈ 1:length(plots) #write the graphs
    draw(PDF("$outputPath\\$(plotNames[i]).pdf", 9inch, 7inch), plots[i])
    #println("Graph $(plotNames[i]) written.")
  end

  println("Cross-sectional Return Graphs by Category Drawn")
  return nothing

end

#contains the functions for writing the graphs by category
function wealthFiguresByCategory(data::NCCSData; assetYears::Vector{Int} = ASSET_YEARS,
    outputPath::String = GRAPH_PATH,
    minPointsPerCategory::Int = MIN_POINTS_PER_CATEGORY,
    focalField::Symbol=WEALTH_FOCAL_FIELD,
    categoryFigureField::Symbol = CATEGORY_FIGURE_FIELD,
    categoriesUsed::Vector = getNonmissingCategories(data, categoryField = categoryFigureField))

  local plotNames::Vector{String} = Vector{String}()
  local plots::Vector{PlotContainer} = Vector{PlotContainer}()
  local subdf::SubDataFrame
  local presubdf::SubDataFrame

  #get only the categories with valid values
  presubdf = view(data.df, ((y::Int)-> y ∈ assetYears).(data.df[:fisyr]),:)
  #presubdf = view(presubdf, presubdf[focalField] .≥ figuresMinAssets)

  push!(plotNames, "Wealth_Distribution_of_Non-Profits_by_Year")
  push!(plots, plot(presubdf,
    x=focalField, color=:fisyr,
    Guide.title("Wealth of Non-Profits by Year"),
    Guide.xlabel("Wealth"), Guide.ylabel("Density"),
    Scale.x_log10, Geom.line,
    Stat.histogram(bincount=50, density=true),
    Coord.cartesian(ymin=0.0, ymax = 0.8, xmin=5.)),)

  push!(plotNames, "Wealth_Distribution_of_Public_Charities_by_Year")
  push!(plots, plot(presubdf[presubdf[:charitytype] .== :pc, :],
    x=focalField, color=:fisyr,
    Guide.title("Wealth of Public Charities by Year"),
    Guide.xlabel("Wealth"), Guide.ylabel("Density"),
    Scale.x_log10, Geom.line,
    Stat.histogram(bincount=50, density=true),
    Coord.cartesian(ymin=0.0, ymax = 1., xmin=5.)))

  push!(plotNames, "Wealth_Distribution_of_Private_Foundations_by_Year")
  push!(plots, plot(presubdf[presubdf[:charitytype] .== :pf,:],
    x=focalField, color=:fisyr,
    Guide.title("Wealth of Private Foundations by Year"),
    Guide.xlabel("Wealth"), Guide.ylabel("Density"),
    Scale.x_log10, Geom.line,
    Stat.histogram(bincount=50, density=true),
    Coord.cartesian(ymin=0.0, ymax = 1., xmin=5.)))

  #filter out funds with no category
  presubdf = view(presubdf, (!ismissing).(presubdf[categoryFigureField]),:)

  #=for c::MSymbol ∈ categoriesUsed

    #filter to the category and exclude non-studied asset years
    subdf = view(presubdf, Vector{Bool}(presubdf[categoryFigureField] .== c), :)

    #we need to make sure we have a viable sample
    if size(subdf,1) ≥ minPointsPerCategory
      push!(plotNames, "categories\\Wealth_Distribution_of Non-Profit_($c)")
      push!(plots, plot(subdf,
        x=focalField, color=:fisyr,
        Guide.title("Wealth Distribution of Non-Profits ($c)"),
        Guide.xlabel("Wealth"), Guide.ylabel("Amount"),
        Scale.x_log10, Geom.line,
        Stat.histogram(bincount=50, density=false),
        Coord.cartesian(ymin=0.0, ymax = 1000.0)
        ))
    end

    if size(subdf[subdf[:charitytype] .== :pc,:],1) ≥ minPointsPerCategory
      push!(plotNames, "categories\\Wealth_Distribution_of_Public Charities_($c)")
      push!(plots, plot(subdf[subdf[:charitytype] .== :pc,:],
        x=focalField, color=:fisyr,
        Guide.title("Wealth Distribution of Non-Profits ($c)"),
        Guide.xlabel("Wealth"), Guide.ylabel("Amount"),
        Scale.x_log10, Geom.line,
        Stat.histogram(bincount=50, density=false),
        Coord.cartesian(ymin=0.0, ymax = 1000.0)))
    end

    if size(subdf[subdf[:charitytype] .== :pf,:],1) ≥ minPointsPerCategory
      push!(plotNames, "categories\\Wealth_Distribution_of_Private_Foundations_($c)")
      push!(plots, plot(subdf[subdf[:charitytype] .== :pf,:],
        x=:adjnetassets, color=:fisyr,
        Guide.title("Wealth Distribution of Non-Profits ($c)"),
        Guide.xlabel("Wealth"), Guide.ylabel("Amount"),
        Scale.x_log10, Geom.line,
        Stat.histogram(bincount=50, density=false),
        Coord.cartesian(ymin=0.0, ymax = 1000.0, xmin=5.0)))
    end
  end=#

  for i ∈ 1:length(plots) #write the graphs
    draw(PDF("$outputPath\\$(plotNames[i]).pdf", 9inch, 7inch), plots[i])
    #println("Graph $(plotNames[i]) written.")
  end

  println("Cross-sectional Wealth Graphs by Category Drawn.")

  return nothing
end

#contains the functions for writing the graphs by category
function returnFiguresByCategory(data::NCCSData;
    returnYears::Vector{Int} = RETURN_YEARS,
    outputPath::String = GRAPH_PATH,
    minPointsPerCategory::Int = MIN_POINTS_PER_CATEGORY,
    focalField::Symbol=RETURN_FOCAL_FIELD,
    categoryFigureField::Symbol = CATEGORY_FIGURE_FIELD,
    categoriesUsed::Vector = getNonmissingCategories(data, categoryField = categoryFigureField))

  local plotNames::Vector{String} = Vector{String}()
  local plots::Vector{PlotContainer} = Vector{PlotContainer}()
  local subdf::SubDataFrame
  local presubdf::SubDataFrame

  #get only the categories with valid values
  presubdf = view(data.df, ((y::Int)-> y ∈ returnYears).(data.df[:fisyr]),:)
  presubdf = view(presubdf, (!ismissing).(presubdf[:return]), :)

  push!(plotNames, "Return_Distribution_of_Non-Profits_by_Year")
  push!(plots, plot(presubdf,
    x=focalField, color=:fisyr,
    Guide.title("Return of Non-Profits by Year"),
    Guide.xlabel("Return"), Guide.ylabel("Amount"),
     Geom.line,
    Stat.histogram(bincount=50, density=false),
    Coord.cartesian(xmin=-0.5, xmax=0.5, ymin=0.0, ymax=1000.0)))

  #=push!(plotNames, "Return_Distribution_of_Public Charities_by_Year")
  push!(plots, plot(presubdf[presubdf[:charitytype] .== :pc, :],
    x=focalField, color=:fisyr,
    Guide.title("Return of Public Charities by Year"),
    Guide.xlabel("Return"), Guide.ylabel("Amount"),
     Geom.line,
    Stat.histogram(bincount=50, density=false),
    Coord.cartesian(xmin=-0.5, xmax=0.5, ymin=0.0, ymax=1000.0)))

  push!(plotNames, "Return_Distribution_of_Private Foundations_by_Year")
  push!(plots, plot(presubdf[presubdf[:charitytype] .== :pf,:],
    x=focalField, color=:fisyr,
    Guide.title("Return of Private Foundations by Year"),
    Guide.xlabel("Return"), Guide.ylabel("Amount"),
     Geom.line,
    Stat.histogram(bincount=50, density=false),
    Coord.cartesian(xmin=-0.5, xmax=0.5, ymin=0.0, ymax=1000.0)))=#

  #filter out funds with no category
  presubdf = view(presubdf, (!ismissing).(presubdf[categoryFigureField]),:)

  #=for c::MSymbol ∈ categoriesUsed

    #filter to the category and exclude non-studied asset years
    subdf = view(presubdf, Vector{Bool}(presubdf[categoryFigureField] .== c), :)

    #we need to make sure we have a viable sample
    if size(subdf,1) ≥ minPointsPerCategory
      push!(plotNames, "categories\\Return_Distribution_of_Non-Profit_($c)")
      push!(plots, plot(subdf,
        x=focalField, color=:fisyr,
        Guide.title("Return Distribution of Non-Profits ($c)"),
        Guide.xlabel("Return"), Guide.ylabel("Amount"),
         Geom.line,
        Stat.histogram(bincount=50, density=false),
        Coord.cartesian(xmin=-0.5, xmax=0.5, ymin=0.0, ymax=1000.0)))
    end

    if size(subdf[subdf[:charitytype] .== :pc,:],1) ≥ minPointsPerCategory
      push!(plotNames, "categories\\Return_Distribution_of_Public_Charities_($c)")
      push!(plots, plot(subdf[subdf[:charitytype] .== :pc,:],
        x=focalField, color=:fisyr,
        Guide.title("Return Distribution of Non-Profits ($c)"),
        Guide.xlabel("Return"), Guide.ylabel("Amount"),
         Geom.line,
        Stat.histogram(bincount=50, density=false),
        Coord.cartesian(xmin=-0.5, xmax=0.5, ymin=0.0, ymax=1000.0)))
    end

    if size(subdf[subdf[:charitytype] .== :pf,:],1) ≥ minPointsPerCategory
      push!(plotNames, "categories\\Return_Distribution_of_Private Foundations_($c)")
      push!(plots, plot(subdf[subdf[:charitytype] .== :pf,:],
        x=:adjnetassets, color=:fisyr,
        Guide.title("Return Distribution of Non-Profits ($c)"),
        Guide.xlabel("Return"), Guide.ylabel("Amount"),
         Geom.line,
        Stat.histogram(bincount=50, density=false),
        Coord.cartesian(xmin=-0.5, xmax=0.5, ymin=0.0, ymax=1000.0)))
    end
  end=#

  for i ∈ 1:length(plots) #write the graphs
    draw(PDF("$outputPath\\$(plotNames[i]).pdf", 9inch, 7inch), plots[i])
    #println("Graph $(plotNames[i]) written.")
  end

  println("Cross-sectional Return Graphs by Category Drawn.")

  return nothing
end

function binForFigures(data::NCCSData;
  returnYears::Vector{Int} = RETURN_YEARS_FULL,
  binOnFields::Vector{Symbol}=BIN_ON_FIELDS,
  focalField::Symbol=RETURN_FOCAL_FIELD,
  assetsField::Symbol = :adjnetassets)::DataFrame

  #pre-filter the dataframe
  local subdf::SubDataFrame
  local presubdf::SubDataFrame = view(data.df, (y::Int -> y ∈ returnYears).(data.df[:fisyr]), :)
  presubdf = view(presubdf, (!ismissing).(presubdf[focalField]), :)

  for f ∈ binOnFields #get rid of missing values
    presubdf = view(presubdf, (!ismissing).(presubdf[f]), :)
  end

  outdf::DataFrame = unique(presubdf[binOnFields])
  nrows::Int = size(outdf, 1)

  local N::Vector{MInt} = Vector{MInt}(undef, nrows)
  local equalmean::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local valuemean::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local stddev::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  #local meanm1::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  #local meanp1::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local p5::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local p25::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local p50::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local p75::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local p95::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)

  #make an index for the rows- do it this way for performance
  outIndex::Dict = Dict(Tuple(outdf[i, binOnFields])=>i for i::Int ∈ 1:nrows)
  presubdfBybinOnFields::GroupedDataFrame = groupby(presubdf, binOnFields)
  for subdf::SubDataFrame ∈  presubdfBybinOnFields
    r::Int = outIndex[Tuple(subdf[1, binOnFields])]
    focal::Vector = subdf[focalField]

    #compute the appropriate statistics
    N[r] = length(focal)

    if N[r] .≥ 100 #NOTE: this is a hardcoded cutoff
      equalmean[r] = StatsBase.mean(focal)
      valuemean[r] =
        sum(processmissing(subdf[focalField] .* subdf[assetsField])) /
        sum(processmissing(subdf[(!ismissing).(subdf[focalField]), assetsField]))
      stddev[r] = StatsBase.std(focal)
      p5[r] = StatsBase.quantile(focal, .05)
      p25[r] = StatsBase.quantile(focal, .25)
      p50[r] = StatsBase.quantile(focal, 0.5)
      p75[r] = StatsBase.quantile(focal, 0.75)
      p95[r] = StatsBase.quantile(focal, 0.95)
    end

  end

  #make the dataframe
  outdf[:N] = N
  outdf[:equalmean] = equalmean
  outdf[:valuemean] = valuemean
  outdf[:stddev] = stddev #NOTE: Should this be weighted standard deviation???
  outdf[:valuemeanm1] = outdf[:valuemean] .- outdf[:stddev]
  outdf[:valuemeanp1] = outdf[:valuemean] .+ outdf[:stddev]
  outdf[:p5] = p5
  outdf[:p25] = p25
  outdf[:p50] = p50
  outdf[:p75] = p75
  outdf[:p95] = p95

  for f ∈ binOnFields # we want to get some sub-aggregates
    outdf[Symbol(:p5, f)] = Vector{MFloat64}(undef, nrows)
    outdf[Symbol(:p25, f)] = Vector{MFloat64}(undef, nrows)
    outdf[Symbol(:p75, f)] = Vector{MFloat64}(undef, nrows)
    outdf[Symbol(:p95, f)] = Vector{MFloat64}(undef, nrows)

    for subdf::SubDataFrame ∈ groupby(presubdf, f)
      focal = subdf[focalField]
      rows::Vector{Bool} = (outdf[f] .== subdf[1:1, f])
      outdf[rows, Symbol(:p5, f)] = StatsBase.quantile(focal, .05)
      outdf[rows, Symbol(:p25, f)] = StatsBase.quantile(focal, .25)
      outdf[rows, Symbol(:p75, f)] = StatsBase.quantile(focal, .75)
      outdf[rows, Symbol(:p95, f)] = StatsBase.quantile(focal, .95)
    end
  end
  return outdf

end

#contains the functions for writing the graphs by category
function returnFiguresBinned(data::NCCSData;
    returnYears::Vector{Int} = RETURN_YEARS_FULL,
    outputPath::String = GRAPH_PATH,
    focalField::Symbol=RETURN_FOCAL_FIELD,
    categoryFigureField::Symbol = CATEGORY_FIGURE_FIELD)

  local aggdf::DataFrame = binForFigures(data, returnYears=returnYears,
    focalField=focalField, binOnFields = [:fisyr, categoryFigureField])

  local plotNames::Vector{String} = Vector{String}()
  local plots::Vector{PlotContainer} = Vector{PlotContainer}()

  aggdf = aggdf[((!ismissing).(aggdf[categoryFigureField])) .& ((!ismissing).(aggdf[:valuemean])), :]

  sort!(aggdf, [:fisyr, categoryFigureField])
  push!(plotNames, "Mean_Returns_by_Category")
  push!(plots, plot(aggdf,
    x=categoryFigureField, y=:valuemean, ymin = Symbol(:p5,categoryFigureField),
      ymax = Symbol(:p95,categoryFigureField),
    Stat.x_jitter(range=0.5), color=:fisyr,
    Guide.title("Mean Returns by Category"),
    Guide.xlabel("Category"), Guide.ylabel("Return"),
    Geom.point, Geom.errorbar
    #Stat.histogram(bincount=50, density=true),
    #Coord.cartesian(ymin=0.0, ymax = 3.)
  ))

  push!(plotNames, "Mean_Returns_by_Year")
  push!(plots, plot(aggdf,
    x=:fisyr, y=:valuemean,
    ymin = :p5fisyr, ymax = :p95fisyr,
    #Stat.x_jitter(range=0.5),
    color=categoryFigureField,
    Guide.title("Mean Returns by Year"),
    Guide.xlabel("Year"), Guide.ylabel("Return"),
    Geom.errorbar, Geom.line,
    #Stat.histogram(bincount=50, density=true),
    Coord.cartesian(xmin=1999, xmax=2016)))

  sort!(aggdf, [:fisyr, categoryFigureField])
  push!(plotNames, "Median_Returns_by_Category")
  push!(plots, plot(aggdf,
    x=categoryFigureField, y=:p50, ymin = :p25, ymax = :p75,
    Stat.x_jitter(range=0.5), color=:fisyr,
    Guide.title("Median Returns by Category"),
    Guide.xlabel("Category"), Guide.ylabel("Return"),
    Geom.point, Geom.errorbar
    #Stat.histogram(bincount=50, density=true),
    #Coord.cartesian(ymin=0.0, ymax = 3.)
  ))

  push!(plotNames, "Median_Returns_by_Year")
  push!(plots, plot(aggdf,
    x=:fisyr, y=:p50, ymin = :p25, ymax = :p75,
    #Stat.x_jitter(range=0.5),
    color=categoryFigureField,
    Guide.title("Median Returns by Year"),
    Guide.xlabel("Year"), Guide.ylabel("Return"),
    Geom.errorbar, Geom.line,
    #Stat.histogram(bincount=50, density=true),
    Coord.cartesian(xmin=1999, xmax = 2016)
  ))

  for i ∈ 1:length(plots) #write the graphs
    draw(PDF("$outputPath\\$(plotNames[i]).pdf", 9inch, 7inch), plots[i])
    #println("Graph $(plotNames[i]) written.")
  end

  println("Binned return graphs drawn.")

  return nothing
end

function makeDescriptiveFigures(data::NCCSData; assetYears::Vector{Int} = ASSET_YEARS)::Nothing

  Gadfly.push_theme(:default)

  #futures::Vector{Future} = Vector{Future}()

  #push!(futures, @spawn wealthFiguresByYear(data))
  #push!(futures, @spawn wealthFiguresByCategory(data))
  #push!(futures, @spawn returnFiguresByYear(data))
  #push!(futures, @spawn returnFiguresByCategory(data))
  #push!(futures, @spawn returnFiguresBinned(data))

  #(fetch).(futures)

  wealthFiguresByYear(data)
  wealthFiguresByCategory(data)
  #returnFiguresByYear(data)
  #returnFiguresByCategory(data)
  #returnFiguresBinned(data)
  industryCross(data)

  return nothing
end

#oneoff tasks
function oneOffCategories(data::NCCSData)
  metadf::DataFrame = deepcopy(data.df[data.df[:fisyr] .== 2013, [:category, :categoryname, :naics2name, :naics2]])
  dropmissing!(metadf)
  outdf::DataFrame = by(metadf, [:category, :categoryname, :naics2name, :naics2],
    [:category, :categoryname, :naics2name, :naics2] =>
      x -> (n=length(x.category)))

  CSV.write("$WORKING_PATH\\categoryMapping.csv",outdf)

  return nothing
end

function oneOffTasks(data::NCCSData)
  println(data.df[1:10, :datefiscal])

  return nothing
end

#launching point for the descriptive output
function makeDescriptive(;refreshFilter::Bool = true,
    refreshDescriptiveFigures::Bool = true,
    refreshDescriptiveTables::Bool = true, writePrimaryCSV::Bool = false,
    doOneOffTasks::Bool= false,
    refreshβ::Bool = true,
    refreshSP500::Bool = false)

    data::NCCSData = constructPrimaryOutput(refreshFilter=refreshFilter,
      writePrimaryCSV=writePrimaryCSV, refreshSP500=refreshSP500)

    refreshDescriptiveFigures && makeDescriptiveFigures(data)
    refreshDescriptiveTables && makeDescriptiveTables(data)
    refreshβ && processβ(data)

    doOneOffTasks && oneOffTasks(data)

    return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
@testset "EDAData tests" begin show_testset_section() #Scope for test data

using MDDatasets
using CMDimCircuits.CircuitAnalysis #Calculations
using CMDimCircuits.NetwAnalysis #For manipulating matrices
using CMDimCircuits.EDAData
include("netwanalysis_cmptools.jl")

#==Helper functions
===============================================================================#

@testset ".sNp files: Simple write/read/compare" begin show_testset_description()
	#Input data
	Z0 = 50 #Ohms (reference impedance)
	ℓ = 5e-3 #Meters
	Zc = 200 #Ohms (line impedance)
	μ0 = 4pi*1e-7 #F/m
	ϵ0 = 8.854e-12 #H/m
	f = collect(1:.1:100)*1e9

	#Computations
	len_mm = ℓ/1e-3
	filepath_S = "lineS_$(len_mm)m.s2p"
	filepath_Z = "lineZ_$(len_mm)m.s2p"

	ω = 2pi*f
	α = 0
	β = ω*sqrt(μ0*ϵ0)
	γ = α .+ im*β

	#Convert γ type to DataF1 (function of 1 argument):
	γ = DataF1(f, γ)

	#ABCD matrix:
	A = cosh(γ*ℓ)      ;B = sinh(γ*ℓ)*Zc;
	C = sinh(γ*ℓ)/Zc   ;D = cosh(γ*ℓ);
	T = Network(:ABCD, [A B; C D])

	#Convert to S-paramters:
	S = Network(:S, T)
	(s11, s12, s21, s22) = mx2elem(S)
	@test maxerror(s11, S[1,1]) < SPARAM_THRESH
	@test maxerror(s12, S[1,2]) < SPARAM_THRESH
	@test maxerror(s21, S[2,1]) < SPARAM_THRESH
	@test maxerror(s22, S[2,2]) < SPARAM_THRESH

	ts = typeof(S)
	msg = "typeof(S) = $ts\n" *
	      "TODO: Should be LESS SPECIFIC: NetworkParametersRef{:S,2,DataF1}"
	@warn msg

	#Save data to .s2p file:
	EDAData.write_snp(filepath_S, S)
	EDAData.write_snp(filepath_Z, Network(:Z, S))

	#Re-load data:
	Sread = EDAData.read_snp(filepath_S, numports=2)
		@test :S == Symbol(NPType(Sread))
		@test maxerror(Sread, S) < SPARAM_THRESH

	Zdata = EDAData.read_snp(filepath_Z, numports=2)
		@test :Z == Symbol(NPType(Zdata))
		@test maxerror(Network(:S, Zdata), S) < SPARAM_THRESH
end

@testset ".psf files: Simple write/read/compare" begin show_testset_description()
	@info "      TODO: Add PSFWrite module & tests"
end

@testset ".tr0 files: Simple write/read/compare" begin show_testset_description()
	@info "      TODO: Add tests"
end

end


text:
"""
    abstract type AbstractASTConsumer <: Any
Supertype for ASTConsumers.
"""
abstract type AbstractASTConsumer end

"""
    struct ASTConsumer <: AbstractASTConsumer
"""
struct ASTConsumer <: AbstractASTConsumer
    ptr::CXASTConsumer
end

Base.unsafe_convert(::Type{CXASTConsumer}, x::ASTConsumer) = x.ptr
Base.cconvert(::Type{CXASTConsumer}, x::ASTConsumer) = x


text:
solve(i::BudgetedElementaryPathInstance{T}, ::DynamicProgramming; kwargs...) where T = budgeted_lp_dp(i; kwargs...)

function budgeted_lp_dp(i::BudgetedElementaryPathInstance{T}) where T
  V = Dict{Tuple{T, Int}, Float64}()
  S = Dict{Tuple{T, Int}, Vector{Edge{T}}}()

  # Initialise. For β = 0, this is exactly Bellman-Ford algorithm with costs
  # (instead of rewards). Otherwise, use the same initialisation as
  # Bellman-Ford.
  β0 = lp_dp(ElementaryPathInstance(graph(i), rewards(i), src(i), dst(i)))
  for v in vertices(graph(i))
    S[v, 0] = β0.solutions[v]
    V[v, 0] = length(S[v, 0]) == 0 ? 0 : sum(rewards(i)[e] for e in S[v, 0])
  end

  for β in 1:budget(i)
    for v in vertices(graph(i))
      V[v, β] = -Inf
      S[v, β] = Edge{T}[]
    end

    V[src(i), β] = 0.0
  end

  # Dynamic part.
  for β in 1:budget(i)
    # Loop needed when at least a weight is equal to zero. TODO: remove it when all weights are nonzero?
    while true
      changed = false

      for v in vertices(graph(i))
        for w in inneighbors(graph(i), v)
          # Compute the remaining part of the budget still to use.
          remaining_budget = max(0, β - weight(i, w, v))

          # If the explored subproblem has no solution, skip it.
          if V[w, remaining_budget] == -Inf
            continue
          end

          # If using the solution to the currently explored subproblem would
          # lead to a cycle, skip it.
          if any(src(e) == v for e in S[w, remaining_budget])
            continue
          end

          # Compute the amount of budget already used by the solution to the currently explored subproblem.
          used_budget = 0
          if length(S[w, remaining_budget]) > 0
            used_budget = sum(weight(i, src(e), dst(e)) for e in S[w, remaining_budget])
          end
          if used_budget < remaining_budget
            continue
          end

          # Compute the maximum: is passing through w advantageous?
          if V[w, remaining_budget] + reward(i, w, v) > V[v, β] && used_budget >= remaining_budget
            changed = true

            V[v, β] = V[w, remaining_budget] + reward(i, w, v)
            S[v, β] = vcat(S[w, remaining_budget], Edge(w, v))
          end
        end
      end

      if ! changed
        break
      end
    end
  end

  return BudgetedElementaryPathSolution(i, S[dst(i), budget(i)], V, S)
end


text:
mutable struct FastNode{G<:Group} <: OnlineStat{XY}
    stats::Vector{G}
    ids::Vector{Int}  # self, left, right
    j::Int
    at::Float64
    ig::Float64
end
function FastNode(p, nkeys=2; stat = FitNormal(), id=1)
    FastNode([p * stat for k in 1:nkeys], [id], 0, 0.0, -Inf)
end
Base.show(io::IO, o::FastNode) =
    print(io, "FastNode | $(nkeys(o)) keys × $(nvars(o)) vars | at=$((o.j, o.at))")

function _fit!(o::FastNode, xy)
    x, y = xy
    y in 1:nkeys(o) || error("y must be an integer in 1:nkeys")
    _fit!(o.stats[y], x)
end
nobs(o::FastNode) = sum(nobs, o.stats)
probs(o::FastNode) = nobs.(o.stats) ./ nobs(o)
nkeys(o::FastNode) = length(o.stats)
nvars(o::FastNode) = length(o.stats[1])
Base.getindex(o::FastNode, i) = [stat[i] for stat in o.stats]

_merge!(o::FastNode, o2::FastNode) = _merge!.(o.stats, o2.stats)

function fakedata(::Type{FastNode}, n, p)
    x = randn(n, p)
    y = [(rand() > 1 /(1 + exp(xb))) + 1 for xb in x * (1:p)]
    x, y
end
function classify(o::FastNode)
    out = 1
    n = nobs(o.stats[1])
    for j in 2:nkeys(o)
        n2 = nobs(o.stats[j])
        if n2 > n
            out = j
            n = n2
        end
    end
    out
end

function whichchild(o::FastNode, x::VectorOb)
    o.j == 0 && return rand(o.ids[2:3])
    x[o.j] < o.at ? o.ids[2] : o.ids[3]
end

impurity(p) = entropy(p, 2)

# tree needs `length` method
function split!(o::FastNode, tree)
    n = nobs(o)
    nl = zeros(nkeys(o))  # "prob" left
    nr = zeros(nkeys(o))  # "prob" right
    imp_before = impurity(probs(o))
    ig = -Inf
    ind = 0
    at = -Inf
    sc = zeros(9 * nkeys(o))
    for j in 1:nvars(o)
        stats_j = o[j]
        k = 0
        for stat in stats_j # split candidates
            μ = mean(stat)
            σ = std(stat)
            sc[k+1] = μ - 2σ
            sc[k+2] = μ - 1.5σ
            sc[k+3] = μ - σ
            sc[k+4] = μ - .5σ
            sc[k+5] = μ
            sc[k+6] = μ + .5σ
            sc[k+7] = μ + σ
            sc[k+8] = μ + 1.5σ
            sc[k+9] = μ + 2σ
            k += 9
        end
        for loc in sc
            for k in 1:nkeys(o)
                stat_kj = stats_j[k]  # summarizes variable j for label k
                # @show cdf(stat_kj, loc), std(stat_kj), nobs(stat_kj)
                nl[k] = cdf(stat_kj, loc) * nobs(stat_kj)
                nr[k] = nobs(stat_kj) - nl[k]
            end
            imp_l = impurity(nl ./ sum(nl))
            imp_r = impurity(nr ./ sum(nr))
            imp_after = smooth(imp_l, imp_r, sum(nr) / n)
            new_ig = imp_before - imp_after
            if new_ig > ig
                ig = new_ig
                ind = j
                at = loc
            end
        end
    end
    o.j = ind
    o.at = at
    o.ig = ig
    d = length(tree)
    push!(o.ids, d + 1)
    push!(o.ids, d + 2)
    FastNode(nvars(o), nkeys(o); id = d + 1), FastNode(nvars(o), nkeys(o); id = d + 2)
end

#-----------------------------------------------------------------------# FastTree
"""
    FastTree(p::Int, nclasses=2; stat=FitNormal(), maxsize=5000, splitsize=1000)

Calculate a decision tree of `p` predictors variables and classes `1, 2, …, nclasses`.
Nodes split when they reach `splitsize` observations until `maxsize` nodes are in the tree.
Each variable is summarized by `stat`, which can be `FitNormal()` or `Hist(nbins)`.

# Example

    x = randn(10^5, 10)
    y = rand([1,2], 10^5)

    o = fit!(FastTree(10), zip(eachrow(x),y))

    xi = randn(10)
    classify(o, xi)
"""
struct FastTree{T<:FastNode} <: OnlineStat{XY}
    tree::Vector{T}
    maxsize::Int
    splitsize::Int
end
function FastTree(p, nkeys=2; stat = FitNormal(), maxsize=5000, splitsize=1000)
    tree = [FastNode(p, nkeys; stat=stat)]
    FastTree(tree, maxsize, splitsize)
end
function Base.show(io::IO, o::FastTree)
    print(io, "FastTree(n=", nobs(o))
    print(io, ", size=", length(o.tree))
    print(io, ", maxsize=", o.maxsize)
    print(io, ", splitsize=", o.splitsize)
    print(io, ")")
end
nobs(o::FastTree) = nobs(o.tree[1])
nkeys(o::FastTree) = nkeys(o.tree[1])
nvars(o::FastTree) = nvars(o.tree[1])
Base.length(o::FastTree) = length(o.tree)
function _fit!(o::FastTree, xy)
    x, y = xy
    node = whichleaf(o, x)
    _fit!(node, xy)
    if length(o.tree) < o.maxsize && nobs(node) > o.splitsize
        left, right = split!(node, o.tree)
        push!(o.tree, left)
        push!(o.tree, right)
    end
end

function whichleaf(o, x::VectorOb)
    node = o.tree[1]
    while length(node.ids) > 1
        node = o.tree[whichchild(node, x)]
    end
    node
end

classify(o::FastTree, x::VectorOb) = classify(whichleaf(o, x))
function classify(o::FastTree, x::AbstractMatrix)
    [classify(o, xi) for xi in OnlineStatsBase.eachrow(x)]
end

#-----------------------------------------------------------------------# FastForest
"""
    FastForest(p, nkeys=2; stat=FitNormal(), kw...)

Calculate a random forest where each variable is summarized by `stat`.

# Keyword Arguments

- `nt=100)`: Number of trees in the forest
- `b=floor(Int, sqrt(p))`: Number of random features for each tree to receive
- `maxsize=1000`: Maximum size for any tree in the forest
- `splitsize=5000`: Number of observations in any given node before splitting
- `λ = .05`: Probability that each tree is updated on a new observation

# Example

    x, y = randn(10^5, 10), rand(1:2, 10^5)

    o = fit!(FastForest(10), zip(eachrow(x),y))

    classify(o, x[1,:])
"""
mutable struct FastForest{T<:FastTree} <: OnlineStat{XY}
    forest::Vector{T}
    subsets::Matrix{Int}
    p::Int
    λ::Float64
    n::Int
end
function FastForest(p, nkeys=2; stat = FitNormal(), maxsize=1000, splitsize = 5000,
                    nt=100, b=floor(Int, sqrt(p)), λ = 5 * (1 / nt))
    forest = [FastTree(b, nkeys; stat=stat, maxsize=maxsize, splitsize=splitsize) for i in 1:nt]
    subsets = fill(0, b, nt)
    for j in 1:size(subsets, 2)
        subsets[:, j] = sample(1:p, b; replace=false)
    end
    FastForest(forest, subsets, p, λ, 0)
end
nkeys(o::FastForest) = maximum(nkeys(tree) for tree in o.forest)
function Base.show(io::IO, o::FastForest)
    print(io, "FastForest(")
    print(io, "n=", nobs(o))
    print(io, ", nt=", length(o.forest))
    print(io, ", b=", size(o.subsets, 1))
    print(io, ", avg_size=", mean(length, o.forest))
    print(io, ")")
end
function _fit!(o::FastForest, xy)
    x, y = xy
    o.n += 1
    for (tree, subset) in zip(o.forest, OnlineStatsBase.eachcol(o.subsets))
        rand() < o.λ && fit!(tree, (x[subset], y))
    end
end

function _classify(o::FastForest, x::VectorOb, buffer::Vector{Int})
    buffer[:] .= 0
    for (i, tree) in enumerate(o.forest)
        buffer[classify(tree, x[o.subsets[:, i]])] += 1
    end
    findmax(buffer)[2]
end

classify(o::FastForest, x::VectorOb) = _classify(o, x, zeros(Int, nkeys(o)))
function classify(o::FastForest, x::AbstractMatrix)
    buffer = zeros(Int, nkeys(o))
    [_classify(o, xi, buffer) for xi in OnlineStatsBase.eachrow(x)]
end


text:
#=
Styles for

Head
Subhead
DataEl
=#

abstract AbstractStyleElement
abstract AbstractPresentationStyleElement <: AbstractStyleElement

# const stylenames = [:EmptyStyleElement, :HeadStyleElement, :SubheadStyleElement, :DataStyleElement]
const stylenames = [:ExpandStyleElement, :CartesianStyleElement, :UniformStyleElement]

for r in stylenames
  @eval begin
    immutable $r <: AbstractPresentationStyleElement
      d::Dict{Symbol, Any}
    end
    $r() = $r(Dict{Symbol, Any}())
  end
end

style(::Type{AbstractStyleElement}, data,
      styleel::AbstractPresentationStyleElement) = error("styleel must be concrete")


text:
# Test of view usage

const Rho_       = 1
const Ux_        = 2
const Uy_        = 3
const Uz_        = 4
const Bx_        = 5
const By_        = 6
const Bz_        = 7
const P_         = 8

function flux_test()

   nI, nJ, nK, nVar = 40, 40, 40, 8
   nG = 1

   state_GV = zeros(nI+2*nG,nJ+2*nG,nK+2*nG,nVar)
   Flux_XV = zeros(nI+1,nJ,nK,nVar)

   iMin,iMax = 1+nG, nI+nG
   jMin,jMax = 1+nG, nJ+nG
   kMin,kMax = 1+nG, nK+nG

   LS_XV = @view state_GV[iMin-1:iMax,jMin:jMax,kMin:kMax,:]
   RS_XV = @view state_GV[iMin:iMax+1,jMin:jMax,kMin:kMax,:]

   @time for k = 1:nK, j = 1:nJ, i = 1:nI+1
      bL = LS_XV[i,j,k,Bx_]^2 + LS_XV[i,j,k,By_]^2 + LS_XV[i,j,k,Bz_]^2
      Flux_XV[i,j,k,Ux_] = LS_XV[i,j,k,Ux_]^2 / LS_XV[i,j,k,Rho_] +
         LS_XV[i,j,k,P_] + 0.5*bL - LS_XV[i,j,k,Bx_]^2
   end

   @time @. Flux_XV[:,:,:,Ux_] = LS_XV[:,:,:,Ux_]^2 / LS_XV[:,:,:,Rho_] + LS_XV[:,:,:,P_] +
      0.5*(LS_XV[:,:,:,Bx_]^2 + LS_XV[:,:,:,By_]^2 + LS_XV[:,:,:,Bz_]^2) -
      LS_XV[:,:,:,Bx_]^2

   LS_XV = state_GV[iMin-1:iMax,jMin:jMax,kMin:kMax,:]
   RS_XV = state_GV[iMin:iMax+1,jMin:jMax,kMin:kMax,:]

   @time for k = 1:nK, j = 1:nJ, i = 1:nI+1
      bL = LS_XV[i,j,k,Bx_]^2 + LS_XV[i,j,k,By_]^2 + LS_XV[i,j,k,Bz_]^2
      Flux_XV[i,j,k,Ux_] = LS_XV[i,j,k,Ux_]^2 / LS_XV[i,j,k,Rho_] +
         LS_XV[i,j,k,P_] + 0.5*bL - LS_XV[i,j,k,Bx_]^2
   end

   @time @. Flux_XV[:,:,:,Ux_] = LS_XV[:,:,:,Ux_]^2 / LS_XV[:,:,:,Rho_] + LS_XV[:,:,:,P_] +
      0.5*(LS_XV[:,:,:,Bx_]^2 + LS_XV[:,:,:,By_]^2 + LS_XV[:,:,:,Bz_]^2) -
      LS_XV[:,:,:,Bx_]^2

   return
end


text:
export
  QWModel,
  QWModelCont,
  QWModelDiscr,
  QWDynamics,
  graph,
  model,
  parameters

"""
    QWModel

Abstract supertype of all quantum walk models.
"""
abstract type QWModel end

"""
    graph(model)

Returns `graph` element of `model`.
"""
graph(model::QWModel) = model.graph

"""
    QWModelCont

Abstract supertype of all continuous quantum walk models.
"""
abstract type QWModelCont <: QWModel end
"""
    QWModelDiscr

Abstract supertype of all discrete quantum walk models.
"""
abstract type QWModelDiscr <: QWModel end

"""
    QWDynamics

Abstract supertype of all dynamics on quantum walk models.
"""
abstract type QWDynamics{T<:QWModel} end

"""
    graph(qwd)

Returns `graph` element of `model` from dynamics `qwd`. Equivalent to
`graph(model(qwd))`.
"""
graph(qwd::QWDynamics) = graph(qwd.model)

"""
    model(qwd)

Returns `model` element of dynamics `qwd`.
"""
model(qwd::QWDynamics) = qwd.model

"""
    parameters(qwd)

Returns `parameters` element of dynamics `qwd`.
"""
parameters(qwd::QWDynamics) = qwd.parameters

"""
    measure(qwd::QWDynamics, state[, vertices::Vector{Int}])

Measure `state` according to `qwd`. If `vertices` is provided, probabilities of
given vertices is returned. Otherwise full probability distribution is returned.
Output is of type `Vector{Float64}`.
"""
measure

"""
    evolve(qwd::QWDynamics{<:QWModelDiscr}, state)
    evolve(qwd::QWDynamics{<:QWModelCont}, state, time::Real)

Evolve `state` according to `qwd`. For discrete model single-step evolution
is implemented. Type returned is the same as type of `state`.
"""
evolve

"""
    check_qwdynamics(qwdtype::Type{<:QWDynamics}, model::QWModel, parameters::Dict{Symbol}, ...)
	check_qwdynamics(qwd)

Checks whetver combination of the arguments creates valid quantum walk dynamics
`qwdtype`. Check whether `qwd` is properly parametrized dynamic.
"""
check_qwdynamics


text:
function per_modstokes(F1, F2, L, lambda, xnu=[], ynu=[]; zerok0::Bool=false, ifgrad::Bool=false)
    # Check input
    M = size(F1)
    @assert length(M)>1 && M[1]==M[2] "Grid must be square"
    @assert size(F1)==size(F2)
    N = M[1]
    # Do FFTs in-place
    Fhat1, Fhat2 = fft(F1), fft(F2)
    # Convolve with truncated greens function
    Uhat1, Uhat2 = Fhat1, Fhat2 # Operate in-place
    k1, k2 = k_vectors([N, N], [L, L])
    k1 = ifftshift(k1)
    k2 = ifftshift(k2)
    # Store k==0 values
    Fhat1k0 = Fhat1[1,1]
    Fhat2k0 = Fhat2[1,1]
    for j=1:N
        for i=1:N
            K1 = k1[i]
            K2 = k2[j]
            Ksq = K1*K1 + K2*K2
            G = 1/ (Ksq*(Ksq + lambda^2))
            KdotFhat = K1*Fhat1[i,j] + K2*Fhat2[i,j]
            Uhat1[i,j] = G*(Ksq*Fhat1[i,j] - K1*KdotFhat)
            Uhat2[i,j] = G*(Ksq*Fhat2[i,j] - K2*KdotFhat)
        end
    end
    # Correct at k==0 (where just just got div by zero)
    Uhat1[1,1] = Fhat1k0 / lambda^2
    Uhat2[1,1] = Fhat2k0 / lambda^2

    if zerok0
        # This improves conditioning, but gets answer wrong by a constant
        # (which we can let integral equation solver pick up)
        # NOTE: It's possible that this is true only because we are using an
        # ill-conditioned test case
        info("Zeroing k=0 mode")
        Uhat1[1,1], Uhat2[1,1] = 0.0, 0.0
    end

    Nnu = length(xnu)
    if Nnu > 0
        # NUFFT for off-grid points
        @assert length(xnu)==length(ynu)
        scale = 2*pi/L
        # Rescale and shift points from [-L/2,L/2] to [0,2*pi]
        xnu *= scale
        ynu *= scale
        xnu += pi
        ynu += pi
        # Compute
        unu1 = Array{Complex128}(Nnu)
        unu2 = Array{Complex128}(Nnu)
        opts = finufft_default_opts()
        opts.modeord = 1
        nufft2d2!(xnu, ynu, unu1, 1, 1e-15, Uhat1, opts)
        nufft2d2!(xnu, ynu, unu2, 1, 1e-15, Uhat2, opts)
        # Get real parts and rescale values
        unu1 = real(unu1)/N^2
        unu2 = real(unu2)/N^2
    else
        unu1, unu2 = [], []
    end

    # Gradient
    if ifgrad
        U1x = zero(Uhat1)
        U1y = zero(Uhat1)
        U2x = zero(Uhat2)
        U2y = zero(Uhat2)
        for j=1:N
            for i=1:N
                U1x[i,j] = 1im*k1[i]*Uhat1[i,j]
                U1y[i,j] = 1im*k2[j]*Uhat1[i,j]
                U2x[i,j] = 1im*k1[i]*Uhat2[i,j]
                U2y[i,j] = 1im*k2[j]*Uhat2[i,j]
            end
        end
        ifft!(U1x)
        ifft!(U1y)
        ifft!(U2x)
        ifft!(U2y)
        U1x, U1y, U2x, U2y = real(U1x), real(U1y), real(U2x), real(U2y)
    else
        U1x, U1y, U2x, U2y = [], [], [], []
    end

    # Transform back in-place
    ifft!(Uhat1)
    ifft!(Uhat2)
    # Return real parts
    U1 = real(Uhat1)
    U2 = real(Uhat2)
    return U1, U2, unu1, unu2, U1x, U1y, U2x, U2y
end

# Compute gradient of real periodic function using Fourier transform
function per_gradient(F, L)
    # Check input
    M = size(F)
    @assert length(M)>1 && M[1]==M[2] "Grid must be square"
    N = M[1]
    Fhat = fft(F)
    k1, k2 = k_vectors([N, N], [L, L])
    k1 = ifftshift(k1)
    k2 = ifftshift(k2)
    Fx = zeros(Fhat)
    Fy = zeros(Fhat)
    for j=1:N
        for i=1:N
            K1 = k1[i]
            K2 = k2[j]
            Fx[i,j] = 1im*K1*Fhat[i,j]
            Fy[i,j] = 1im*K2*Fhat[i,j]
        end
    end
    # Transform back in-place
    ifft!(Fx)
    ifft!(Fy)
    # Return real parts
    return real(Fx), real(Fy)
end


text:
using EDFPlus
using Test

edfh = loadfile("EDFPlusTestFile.edf")
sz = size(edfh.EDFsignals)
@test sz == (20010,601)

eann = Annotation()
@test eann.onset == 0.0
ann = Annotation(61.04, "5.25", "They said schöner")
addannotation!(edfh, ann.onset, ann.duration, ann.annotation)
EDFPlus.latintoascii("")
edfh.gender = "Male"
edfh.signalparam[1].label = "Signal"
@test EDFPlus.readBiosemiStatus(edfh) isa Dict

newedfh = writefile!(edfh, "NEWedfplustestfile.edf")
@test size(newedfh.EDFsignals) == sz
@test EDFPlus.bytesperdatapoint(newedfh) == 2
closefile!(newedfh)

bdfh = loadfile("samplefrombiosemicom.bdf")
bsz = size(bdfh.BDFsignals)
@test bsz == (60, 34816)

newbdfh = writefile!(bdfh, "NEWsamplefrombiosemicom.bdf")
@test size(newbdfh.BDFsignals) == bsz
closefile!(newbdfh)
closefile!(bdfh)

newbdfh = writefile!(edfh, "NEWbdfplusfromedfplus.bdf", sigformat=EDFPlus.bdfplus)
@test size(newbdfh.BDFsignals) == sz
@test EDFPlus.bytesperdatapoint(newbdfh) == 3

newedfh2 = writefile!(newbdfh, "NEWedfplusfrombdfplus.edf", sigformat=EDFPlus.edfplus)
@test size(newedfh2.EDFsignals) == sz
@test EDFPlus.bytesperdatapoint(newedfh2) == 2

true


text:
using MPIMeasurements
using Plots
using Unitful
using UnitfulRecipes

pyplot()
default(show = true)

ENV["JULIA_DEBUG"] = "all"

# Add test configurations to path
testConfigDir = normpath(string(@__DIR__), "../test/Scanner/TestConfigs")
addConfigurationPath(testConfigDir)

scannerName_ = "TestSimpleSimulatedScanner"
scanner = MPIScanner(scannerName_)

path = normpath(testConfigDir, "TestSimpleSimulatedScanner/Sequences/Sequence.toml")
sequence = sequenceFromTOML(path)
setupSequence(scanner, sequence)

uMeas, uRef, t = readData(getDAQ(scanner), 1, 5)
showFrames = 1:5
plot(t[:, 1, showFrames], uRef[:,1,1,showFrames])
plot!(t[:, 1, showFrames], uMeas[:,1,1,showFrames])


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

using Markdown
using InteractiveUtils

# ╔═╡ c9aa684a-43a8-11eb-3400-b383901a7f05
begin
	using Pkg
	Pkg.add(url = "https://github.com/heltonmc/LightPropagation.git")
	using LightPropagation
end

# ╔═╡ 4399d234-43a5-11eb-2389-21beb59e41fd
md"# _When is semi-infinite semi-infinite?_

### Learning outcomes:
- Simulate the time-domain reflectance under the diffusion approximation for **semi-finite, slab, and parralelepiped** geometries
- Utilize the non-linear curve fitting functions to extract optical properties from simulated experimental data
- Gain physical intuition of approximations and how to apply them
"

# ╔═╡ 31c8aab0-43a7-11eb-3fea-417ceb3ba0b5
md"
!!! info \"Disclaimer\"
	The following is not a recommendation on when to apply certain approximations as       these will depend on experimental setup and application. Please use this as an         introduction to the software.
"

# ╔═╡ 117d72ca-43a9-11eb-1394-e584e55e76e9
md"
##### Loading LightPropagation.jl
Let's begin by loading the package into the current environment. This can be done as below or as shown on the main page.
"

# ╔═╡ f7a02712-43a8-11eb-1cff-eb83d99a0989
md"
## Forward Models

The functions we are using to simulate the TPSF for the semi-infinite, slab, and parallelpiped geometries are:
- **TPSF\_DA\_semiinf\_refl**
- **TPSF\_DA\_slab_refl**
- **TPSF\_DA\_paralpip_refl**

"

# ╔═╡ f70fa3ee-43a9-11eb-1599-c9d55ca04787
md"
!!! tip \"Available Functions\"
	To view all the available functions it is often useful to check the exported 		functions in the main file. For example if you go to                                   ~/LightPropagation/src/LightPropagation.jl in the repository you will find             `export TPSF_DA_semiinf_refl` lines that indicate available functions.
"

# ╔═╡ acd4eee0-43ab-11eb-14f7-f95c2291d580
md"
Let's simulate the TPSF for a semi-infinite medium under the following conditions:
- t = 0:0.01:5 (ns)
- β = [0.1, 10.0] -- β[1] and β[2] are the absorption and reduced scattering coefficient (1/cm)
- ρ = 1.0 -- source-detector separation (cm)
- ndet = 1.0 -- index of refraction of detector
- nmed = 1.0 -- index of refraction of medium
"

# ╔═╡ 27dee8a8-43ab-11eb-33a2-81af1f08c66d
begin
	t = 0.01:0.01:5
	TPSF_semiinf = TPSF_DA_semiinf_refl(t, [0.1,10.0], 1.0, 1.0, 1.0)
end

# ╔═╡ a32eb024-43ab-11eb-09fa-9dc34018f56f
begin
	using Plots
	plot(t, TPSF_semiinf, yscale=:log10, lw = 2, ylabel="Counts", 							xlabel="time (ns)", label="TPSF_semiinf")
end

# ╔═╡ 5864c6da-43ad-11eb-117a-b7270071791d
md"
##### Comparison of semi-infinite and slab
Now that we know how to run a function and plot it, let's compare the semi-infinite approximation to a slab (bounded in the z-direction) for two thickness (2cm and 8cm)
"

# ╔═╡ 2a29c586-43ad-11eb-007f-45911e0d0ca5
begin
	TPSF_slab_2cm = TPSF_DA_slab_refl(t, [0.1,10.0], 1.0, 1.0, 1.0, 2.0)
	TPSF_slab_8cm = TPSF_DA_slab_refl(t, [0.1,10.0], 1.0, 1.0, 1.0, 8.0)
end

# ╔═╡ d5282152-43ae-11eb-3381-5bb1db31fece
begin
	plot(t, TPSF_semiinf, yscale=:log10, lw = 6, ylabel="Counts", 			 			xlabel="time (ns)", label="TPSF_semiinf") #make lw thicker for clarity
	plot!(t, TPSF_slab_2cm, lw = 2, label = "TPSF_slab_2cm")
	plot!(t, TPSF_slab_8cm, lw = 2, label = "TPSF_slab_8cm", c="black")
end

# ╔═╡ 5bd4961c-43b0-11eb-0d78-c3850f2f643d
md"
Great. Now change the input parameters for each curve (optical properties and SDS) to investigate different domains. Here, for this set of optical properties and SDS the  semi-infinite approximation matches exactly the slab solution when the slab thickness is 8cm (for t<5ns), though it would be a poor approximation for thin samples (2cm).
"

# ╔═╡ e98e6dae-43b3-11eb-2e8a-01012c7befe2
md"
#### Compare all three geometries
Let's introduce the parralelepiped geometry where it is bounded in the x, y, and z direction. Remember the inputs look like this:

TPSF\_DA\_paralpip\_refl(t, β::Array{Float64,1}, ndet::Float64, nmed::Float64, rd::Array{Float64,1}, rs::Array{Float64,1}, L::Array{Float64,1})

You can also check arguments by typing ? in the REPL and typing in the function name as shown in main page. Here, rd is the x, y location of the source on the top surface, and rs is the location of the detector on the top source. These are measured from the corner. Let's have two cubes one with L = 8cm and L = 4cm with the source located in the center so that rs = [4.0, 4.0] and [2.0, 2.0] respectively. We want the detector to be 1cm from the source to match SDS in previous examples. Below they are defined as rd = [3.0, 4.0] and [1.0, 2.0].
"

# ╔═╡ d40bc0b0-43b0-11eb-1f81-597d8d2a017b
begin
	TPSF_paralpip_8cm = TPSF_DA_paralpip_refl(t, [0.1,10.0], 1.0, 1.0, [3.0,4.0], 			[4.0,4.0], [8.0,8.0,8.0])
	TPSF_paralpip_4cm = TPSF_DA_paralpip_refl(t, [0.1,10.0], 1.0, 1.0, [1.0,2.0], 			[2.0,2.0], [4.0,4.0,4.0])
end

# ╔═╡ 0cb1c518-43b6-11eb-3bac-052d89f819ae
begin
	plot(t, TPSF_semiinf, yscale=:log10, lw = 10, ylabel="Counts", 			 			xlabel="time (ns)", label="TPSF_semiinf") #make lw thicker for clarity
	plot!(t, TPSF_slab_2cm, lw = 2, label = "TPSF_slab_2cm")
	plot!(t, TPSF_slab_8cm, lw = 4, label = "TPSF_slab_8cm", c="black")
	plot!(t, TPSF_paralpip_8cm, lw = 2, label = "TPSF_paralpip_8cm", alpha = 0.9)
	plot!(t, TPSF_paralpip_4cm, lw = 2, label = "TPSF_paralpip_4cm", alpha = 0.9)

end

# ╔═╡ 9a7d26ee-4479-11eb-1eb0-f79e3c07fb28
pwd()

# ╔═╡ 46710c9e-43b8-11eb-0172-65d7f2171f56
md"
#### Key Takeaways
- Good agreement at the earliest times (<1ns) where the boundary effects are minimized.
- For confined geometries, (e.g. TPSF\_paralpip\_4cm and TPSF\_slab\_2cm) we see there are significant differences compared to the semi-infinite approximations.
- If the boundaries are significantly away from the source the semi-infinite approximation appears to be a good approximation.
"

# ╔═╡ 3cb5c3e2-43b9-11eb-2299-7503e415f13c
md"
## Inverse Problem

Now that we have visually observed the differences between the three geometries, let's investigate and try to quantify the error in the extracted optical properties when using different approximations.

The first thing we must do is simulate some experimental data. In practice, samples are not infinite so let's use the parralelepiped geometry as the ground 'truth' and let's fit all three models to that data.
"

# ╔═╡ d2ee9b40-43eb-11eb-1f68-7f94589ba3f4
ydata = TPSF_DA_paralpip_refl(t, [0.1,10.0], 1.0, 1.0, [0.5,1.5], [1.5,1.5], [3.0,3.0,3.0])

# ╔═╡ 2d46a2d4-43b7-11eb-10fa-03d55dba74e8
md"
##### Setup Inverse problem

The first thing we must do is set up the input data and model parameters. These use two separate structures. The input data has the following structure:
```julia
@with_kw struct fitDTOF

    #required inputs
    t::Union{AbstractRange{Float64}, AbstractVector{Float64}}
    DTOF::Array{Float64,1}
    IRF::Array{Float64,1}

    nmed::Float64 = 1.5
	ndet::Float64 = 1.51
	lambda::Float64 = 750.0

    #semi-infinite and slab
    ρ::Union{Float64, Missing} = missing

    #slab
    s::Union{Float64, Missing} = missing

    #parallelepiped
    rd::Union{Array{Float64,1}, Missing} = missing
    rs::Union{Array{Float64,1}, Missing} = missing
    L::Union{Array{Float64,1}, Missing} = missing
end
```
"

# ╔═╡ 27b1a654-43ec-11eb-1b4f-9f165debb017
md"
We can see that we must also describe and IRF. For now, let's assume the IRF is as defined below.
"

# ╔═╡ 6f0f914c-43c3-11eb-2561-f7c8cae30e99
IRF = [1; zeros(length(ydata) - 1)]

# ╔═╡ 369062ea-43c3-11eb-3143-a562cfe3afe5
data = fitDTOF(t = t, DTOF = ydata./maximum(ydata), IRF = IRF, ρ = 1.0, nmed = 1.0, ndet = 1.0, s = 3.0, rd = [0.5,1.5], rs = [1.5,1.5], L = [3.0,3.0,3.0])

# ╔═╡ 495b4f02-43c5-11eb-3ea6-6d35ebc81514
m = DTOF_fitparams(model = TPSF_DA_paralpip_refl, risefactor = 0.01, tailfactor = 1e-10,initparams = [0.05, 20.0])

# ╔═╡ 6508edfe-43c5-11eb-2d67-4fcc34883150
f, tfit, yraw, yfit  = getfit(data, m)

# ╔═╡ 50348fb8-43c8-11eb-2bfc-6102f19bcf0e
begin
	m1 = DTOF_fitparams(model = TPSF_DA_semiinf_refl, risefactor = 0.01, tailfactor = 1e-10, initparams = [0.05, 20.0])
	m2 = DTOF_fitparams(model = TPSF_DA_slab_refl, risefactor = 0.01, tailfactor = 		1e-10, initparams = [0.05, 20.0])
end

# ╔═╡ 6fa2a252-43c8-11eb-2ffd-73dfebb5454e
begin
	f1, tfit1, yraw1, yfit1  = getfit(data, m1)
	f2, tfit2, yraw2, yfit2  = getfit(data, m2)
end

# ╔═╡ 8790e19e-43c8-11eb-11d4-e5b36da13513
begin
	scatter(tfit, yraw, m="*", ms = 8, label = "ydata", c="black")
	plot!(tfit, yfit, lw = 3, alpha = 1, label = "paralpip model", c="blue")
	plot!(tfit1, yfit1, lw = 3, alpha = 1,label = "semiinf model", c="red")
	plot!(tfit2, yfit2,lw = 3, alpha = 1, label = "slab model", c="orange")
	plot!(ylabel="Counts", xlabel="time (ns)")
	savefig("fitcomparison.png")
end

# ╔═╡ a8e788d4-43c8-11eb-3c20-b94cd8884937
f.param

# ╔═╡ af161984-43c8-11eb-0715-079cc75433ec
f1.param

# ╔═╡ b231ab40-43c8-11eb-0819-1d16af7aafb4
f2.param

# ╔═╡ c298b8e6-4488-11eb-2b48-b9a7fcc20740
abs.(f.param .- f2.param)./f.param

# ╔═╡ Cell order:
# ╟─4399d234-43a5-11eb-2389-21beb59e41fd
# ╠═31c8aab0-43a7-11eb-3fea-417ceb3ba0b5
# ╟─117d72ca-43a9-11eb-1394-e584e55e76e9
# ╠═c9aa684a-43a8-11eb-3400-b383901a7f05
# ╟─f7a02712-43a8-11eb-1cff-eb83d99a0989
# ╟─f70fa3ee-43a9-11eb-1599-c9d55ca04787
# ╟─acd4eee0-43ab-11eb-14f7-f95c2291d580
# ╠═27dee8a8-43ab-11eb-33a2-81af1f08c66d
# ╠═a32eb024-43ab-11eb-09fa-9dc34018f56f
# ╟─5864c6da-43ad-11eb-117a-b7270071791d
# ╠═2a29c586-43ad-11eb-007f-45911e0d0ca5
# ╠═d5282152-43ae-11eb-3381-5bb1db31fece
# ╟─5bd4961c-43b0-11eb-0d78-c3850f2f643d
# ╟─e98e6dae-43b3-11eb-2e8a-01012c7befe2
# ╠═d40bc0b0-43b0-11eb-1f81-597d8d2a017b
# ╠═0cb1c518-43b6-11eb-3bac-052d89f819ae
# ╠═9a7d26ee-4479-11eb-1eb0-f79e3c07fb28
# ╟─46710c9e-43b8-11eb-0172-65d7f2171f56
# ╟─3cb5c3e2-43b9-11eb-2299-7503e415f13c
# ╠═d2ee9b40-43eb-11eb-1f68-7f94589ba3f4
# ╟─2d46a2d4-43b7-11eb-10fa-03d55dba74e8
# ╟─27b1a654-43ec-11eb-1b4f-9f165debb017
# ╠═6f0f914c-43c3-11eb-2561-f7c8cae30e99
# ╠═369062ea-43c3-11eb-3143-a562cfe3afe5
# ╠═495b4f02-43c5-11eb-3ea6-6d35ebc81514
# ╠═6508edfe-43c5-11eb-2d67-4fcc34883150
# ╠═50348fb8-43c8-11eb-2bfc-6102f19bcf0e
# ╠═6fa2a252-43c8-11eb-2ffd-73dfebb5454e
# ╠═8790e19e-43c8-11eb-11d4-e5b36da13513
# ╠═a8e788d4-43c8-11eb-3c20-b94cd8884937
# ╠═af161984-43c8-11eb-0715-079cc75433ec
# ╠═b231ab40-43c8-11eb-0819-1d16af7aafb4
# ╠═c298b8e6-4488-11eb-2b48-b9a7fcc20740

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation