text:
using FunctionalCollections
using JSON

import FunctionalCollections: append
export Node, node, instanceof, props

"""
    Node(instanceof, children...; props...)
    Node(instanceof, children, props)

The building block of WebIO.
A `Node` is simply a wrapper around an _instance_ (some Julia object) together
with some child nodes and additional properties.

The most common type of `Node` is a DOM node. These can be constructed just by
specifying a symbol as the `instanceof` (they are promoted to an instance of
`WebIO.DOM` under the hood).
```jldoctest
julia> Node(:div, Node(:p, "I am a paragraph!", class="important"))
(div
  (p { class="important" }
    "I am a paragraph!"))
```

Nodes with custom (non-`DOM`) instances should have a corresponding
`WebIO.render` method defined.
"""
struct Node{T}
    instanceof::T
    children::PersistentVector{Any}
    props::Dict{Symbol, Any}

    # This needs to be an inner constructor to enforce the promotion of the
    # instanceof attribute.
    function Node(instanceof, children::AbstractVector, props::AbstractDict)
        instanceof = promote_instanceof(instanceof)
        return new{typeof(instanceof)}(
            promote_instanceof(instanceof),
            _pvec(children),
            kwargs2props(props),
        )
    end
end

function Node(instanceof, children...; props...)
    return Node(instanceof, collect(Any, children), Dict(props...))
end

# Can/should this be deprecated?
# The relationship between this and Widgets.jl doesn't make very much sense
# (WebIO extends `node` from Widgets but Widgets essentially relies on this
# implentation for things like `Widgets.div`).
node(args...; kwargs...) = Node(args...; kwargs...)

"""
    promote_instanceof(x)

Promotes an `instanceof` (when constructing a new `Node`) to the appropriate
type.

# Examples
```jldoctest
julia> repr(WebIO.promote_instanceof(:div))
"WebIO.DOM(:html, :div)"
```
"""
promote_instanceof(x) = x
promote_instanceof(s::Symbol) = DOM(:html, s)
promote_instanceof(s::AbstractString) = promote_instanceof(Symbol(s))

nodetype(n::Node) = typename(n.instanceof)
typename(n::T) where {T} = string(T.name.name)

"""
    kwargs2props(kwargs)

Convert keyword arguments (à la `kwargs...`) to a `props` dict.
"""
kwargs2props(kwargs) = Dict{Symbol, Any}(kwargs)

"""
    DOM(namespace, tag)

An instance (for `Node.instanceof`) that represents a specific DOM node type.
"""
struct DOM
    namespace::Symbol
    tag::Symbol
end

export children, setchildren, instanceof, setinstanceof, props, setprops

children(n::Node) = n.children
setchildren(n::Node, children) = Node(n.instanceof, PersistentVector{Any}(collect(children)), n.props)
instanceof(n::Node) = n.instanceof
setinstanceof(n::Node, instanceof) = Node(instanceof, n.children, n.props)
props(n::Node) = n.props
setprops(n::Node, props) = Node(n.instanceof, n.children, props)

# Modifications
export append, setchild, withchild, withlastchild, mergeprops

append(n::Node, cs) = setchildren(n, append(children(n), cs))
setchild(n::Node, i, c) = setchildren(n, assoc(children(n), i, c))
withchild(f, n::Node, i) = setchild(n, i, f(children(n)[i]))
withlastchild(f, n::Node) = setchild(n, length(children(n)), f(children(n)[i]))
function mergeprops(n::Node, p, ps...)
    out = recmerge(props(n), p)
    for x in ps
        recmerge!(out, x)
    end
    setprops(n, out)
end

# Element extension syntax
(n::Node)(x, args...) = append(n, (x, args...))
(n::Node)(props::AbstractDict...) = mergeprops(n, props...)
(n::Node)(;kwargs...) = mergeprops(n, kwargs)
(n::Node)(args...; kwargs...) = n(args...)(;kwargs...)

# Pretty printing
showindent(io, level) = write(io, repeat("  ", level))

function Base.show(io::IO, ::MIME"text/plain", el::Node)
    _show(io, el)
end

function showprops(io, dict)
    write(io, "{ ")
    for (k,v) in dict
        print(io, k, "=", repr(v), " ")
    end
    write(io, "}")
end

function showchildren(io, elems, indent_level)
    length(elems) == 0 && return
    write(io, "\n")
    l = length(elems)
    for i=1:l
        _show(io, elems[i], indent_level+1)
        i != l && write(io, "\n")
    end
end

function _show(io::IO, el::Node, indent_level=0)
    showindent(io, indent_level)
    write(io, "(")
    if !isa(el.instanceof, DOM)
        write(io, string(el.instanceof))
        write(io, ":")
    elseif el.instanceof.namespace != :html
        write(io, el.instanceof.namespace)
        write(io, ":")
        write(io, el.instanceof.tag)
    else
        write(io, el.instanceof.tag)
    end

    if !isempty(props(el))
        write(io, " ")
        showprops(io, props(el))
    end
    showchildren(io, children(el), indent_level)
    write(io, ")")
end

function _show(io::IO, el, indent_level)
    showindent(io, indent_level)
    show(io, el)
end


text:
function theta2delta(th,dim)
    return beta_inc((dim-1)/2,1/2,sin(th)^2)[1]
end
function delta2theta(delt,dim)
    sintheta = sqrt(beta_inc_inv((dim-1)/2,1/2,delt)[1])
    return asin(sintheta)
end

function epsilon2beta(ep,dim,numMode,N)
    return numMode*(1-theta2delta(delta2theta(ep,dim)/2,dim)/numMode)^N/theta2delta(delta2theta(ep,dim)/4,dim)
end

function looseningratio_quad(ep,dim,kappaP)
    loosenfactor = (1-kappaP*(1-cos(delta2theta(ep,dim))))
    if loosenfactor <= 0
        println("Insufficient samples!")
        return Inf
    else
        return 1/(1-kappaP*(1-cos(delta2theta(ep,dim))))
    end
end

function beta2epsilon(beta,dim,numMode,N)
    ep_u = 1.0
    ep_l = 0.0
    while ep_u - ep_l > 1e-10
        ep = (ep_u+ep_l)/2
        beta_ep = epsilon2beta(ep,dim,numMode,N)
        if beta_ep > beta
            ep_l = ep
        else
            ep_u = ep
        end
    end
    return ep_u
end

function probabilistc_stability_certificate(state0,state;B,numMode,d,batchsize,K0,beta=0.01,tol=1e-4)
    dim = size(B,1)
    N = size(state0,2)
    K,ga,P = soslyap_alternating(state0,state;B=B,d=d,batchsize=batchsize,K0=K0,tol=tol)
    eigP = eigvals(P)
    kappaP = maximum(eigP)/minimum(eigP)
    if d == 1
        ep = beta2epsilon(beta,dim,numMode,N)
        return K,ga*looseningratio_quad(ep,dim,kappaP),ga,P
    else
        epsos = beta2epsilon(beta,dim,numMode,N)
        stateK = zeros(dim,N)
        for i in 1:N
            stateK[:,i] = state[:,i]+B*K*state0[:,i]
        end
        boundnorm = norm_max(state0,stateK)/cos(delta2theta(epsos,dim))

        phid = 0.0
        for i in 1:d
            phid += (2-2*cos(delta2theta(epsos,dim)))^(i/2)*binomial(d,i)
        end
        bound = (ga^d*(1+sqrt(kappaP)*phid)+boundnorm^d*sqrt(kappaP)*phid)^(1/d)

        #=for betasos in 0.1*beta:0.1*beta:0.9*beta
            epsos = beta2epsilon(betasos,dim,numMode,N)

            epnorm = minimum([numMode*(1-(beta-betasos)^(1/N)),1])

            stateK = zeros(dim,N)
            for i in 1:N
                stateK[:,i] = state[:,i]+B*K*state0[:,i]
            end

            boundnorm = norm_max(state0,stateK)/cos(delta2theta(epnorm,dim))

            phid = 0.0
            for i in 1:d
                phid += (2-2*cos(delta2theta(epsos,dim)))^(i/2)*binomial(d,i)
            end
            boundnormclosed = boundnorm #+ opnorm(B*K)
            bound = (ga^d*(1+sqrt(kappaP)*phid)+boundnormclosed^d*sqrt(kappaP)*phid)^(1/d)
            push!(boundsearch,bound)
        end=#
        #println("Probability bisection: $boundsearch")
        return K,bound,ga,P

    end

end

function norm_max(state0,state)
    traj_norm = 0.0
    numTraj = size(state0)[2]

    for i = 1:numTraj
        traj_norm_i = norm(state[:,i])/norm(state0[:,i])
        if traj_norm_i > traj_norm
            traj_norm = traj_norm_i
        end
    end

    return traj_norm
end

function probabilistc_LQR_certificate(state0,state;B,numMode,Q,R,kappaU,K0,beta=0.01,tol=1e-4)
    dim = size(B,1)
    N = size(state0,2)

    ep = beta2epsilon(beta,dim,numMode,N)

    lowbound = 1-kappaU*(1-cos(delta2theta(ep,dim)))
    #println("lowbound: $lowbound")
    if lowbound < 0
        lowbound = 0
    end

    K, xi, P = LQR_alternating(state0,state;B=B,Q=Q,R=R,kappaU=kappaU,lowbound=lowbound,K0=K0,ep=ep,tol=tol)

    Z = P-Q-K'*R*K
    eigZ = eigvals(Z)

    ratio = 1-(maximum(eigZ)/minimum(eigZ))*(1-cos(delta2theta(ep,dim)))
    #println("ratio: $ratio")
    #println("eigZ: $eigZ")

    if ratio < xi
        println("Infeasible!")
    else
        println("A feasible solution is found!")
    end
    return K,xi,P,ratio

end

ext: jl

max_stars_repo_name: zhemingwang/DataDrivenSwitchControl

max_issues_repo_name: zhemingwang/DataDrivenSwitchControl

max_forks_repo_name: zhemingwang/DataDrivenSwitchControl


text:
using TotalLeastSquares, Plots, ThreadTools

N = 500
σvec = exp10.(LinRange(-3, 3, 40))
results1 = tmap(σvec) do σ
      inner = map(1:N) do _
            x   = randn(5)
            A   = randn(500,5)
            An  = A + σ*randn(size(A)) .* (rand(size(A)...) .< 0.01)
            y   = A*x
            yn  = y + σ*randn(size(y)) .* (rand(size(y)...) .< 0.01)

            x̂t = tls(An,yn)
            x̂r = rtls(An,yn)

            norm(x-x̂t), norm(x-x̂r)
      end
      mean(getindex.(inner, 1)), mean(getindex.(inner, 2))
end

res_tls1, res_rtls1 = getindex.(results1, 1), getindex.(results1, 2)

scatter(σvec, [res_tls1 res_rtls1], xscale=:log10, yscale=:log10, lab=["TLS" "RTLS"], xlabel="Noise std", ylabel="Parameter error norm", legend=:topleft, title="Noise Sparsity = 0.01")

##
σ = 5
svec = exp10.(LinRange(-3, 0, 40))
results2 = tmap(svec) do s
      inner = map(1:N) do _
            x   = randn(5)
            A   = randn(500,5)
            An  = A + σ*randn(size(A)) .* (rand(size(A)...) .< s)
            y   = A*x
            yn  = y + σ*randn(size(y)) .* (rand(size(y)...) .< s)

            x̂t = tls(An,yn)
            x̂r = rtls(An,yn)

            norm(x-x̂t), norm(x-x̂r)
      end
      mean(getindex.(inner, 1)), mean(getindex.(inner, 2))
end

res_tls2, res_rtls2 = getindex.(results2, 1), getindex.(results2, 2)

scatter((svec), [res_tls2 res_rtls2], xscale=:log10, yscale=:log10, lab=["TLS" "RTLS"], xlabel="Noise sparsity", ylabel="Parameter error norm", legend=:topleft, title="Noise std = $σ")


text:
using Statistics: mean

function nansum(
    x;
    dims :: Union{Nothing, Tuple} = nothing,
)

    y = copy(x)

    nan_idx    = isnan.(y)
    finite_idx = isfinite.(y)

    y[nan_idx] .= 0.0

    if dims == nothing
        return sum(y)
    else
        return sum(y, dims=dims)
    end
end

function nanmean(
    x;
    dims :: Union{Nothing, Tuple} = nothing,
)

    y = copy(x)

    nan_idx    = isnan.(y)
    finite_idx = isfinite.(y)

    y[nan_idx] .= 0.0
    if dims == nothing
        return sum(y) ./ sum(finite_idx)
    else
        return sum(y, dims=dims) ./ sum(finite_idx, dims=dims)
    end
end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
mutable struct Microdata
    nonmissing::BitVector
    model::Dict{Symbol, SM.TermOrTerms}
    mapping::Dict{Symbol, Vector{Int}}
    matrix::AbstractMatrix{Float64}
    weights::AbstractWeights{Int, Float64, Vector{Float64}}
    corr::CorrStructure
end

# CONSTRUCTORS

function Microdata(
        df::AbstractDataFrame,
        model::Dict{Symbol, SM.TermOrTerms};
        hints::Dict{Symbol} = Dict{Symbol, Any}(),
        subset::AbstractVector{Bool} = [true],
        corr::CorrStructure = Heteroscedastic(),
        weights::AbstractWeights = UnitWeights{Float64}(size(df, 1))
    )

    columns = termvars(model)

    nonmissing = .!iszero.(weights) .& corr.nonmissing .& subset
    nonmissing[nonmissing] .&= completecases(view(df, nonmissing, columns))

    corr    = update_corr(corr, nonmissing)
    weights = update_weights(weights, nonmissing)
    data    = view(df, nonmissing, columns)
    model   = apply_schema(model, schema(model, data, hints))
    mapping = asgn(model)
    matrix  = modelcols(model, data)

    return Microdata(nonmissing, model, mapping, matrix, weights, corr)
end

function Microdata(MD::Microdata)
    Microdata(
        copy(MD.nonmissing),
        copy(MD.model),
        copy(MD.mapping),
        MD.matrix,
        MD.weights,
        MD.corr
    )
end

max_stars_repo_name: lbittarello/Microeconometrics

max_issues_repo_name: lbittarello/Microeconometrics

max_forks_repo_name: lbittarello/Microeconometrics


text:
using Recombinase, Plots
##
args, kwargs = Recombinase.series2D(
    Recombinase.discrete(FlippingModel.survival),
    b,
    error = :MouseID,
    Recombinase.Group(color = :Protocol, linestyle = :Wall),
    select = (:Omissions_plus_one, :Leave),
    ribbon = true
    )
p1 = plot(args...;
    kwargs...,
    legend = false,
    grid = false,
    thickness_scaling = 2.5,
    tick_orientation = :out,
    size = (809,500),
    linewidth=2,
    fillalpha= 0,
    xlims = (0,16),
    xticks = 1:2:16)
##
savefig(p1,"/Volumes/GoogleDrive/My Drive/Reports for Zach/NewProtocols/Fitting/Survival_from_data.pdf")
##
figures = "/Users/dariosarra/Documents/Lab/Mainen/Presentations/Walton_visit/figures/unprocessed_plots"
savefig(p1,joinpath(figures,"New_data_survival.pdf"))
##
Simulated[!, :Leave] .= true
args, kwargs = series2D(
    Recombinase.discrete(FlippingModel.survival),
    Simulated,
    error = :MouseID,
    Recombinase.Group(color = :Protocol, linestyle = :Wall),
    select = (:simulation, :Leave),
    ribbon = true
    )
p2 = plot!(p1,args...;kwargs...,
    legend = false,
    fillalpha = 0.25,linealpha=0)
##
savefig(p1,"/Volumes/GoogleDrive/My Drive/Reports for Zach/NewProtocols/Fitting/Survival_from_fit.pdf")
##
args, kwargs = series2D(
    Recombinase.discrete(FlippingModel.survival),
    b,
    error = :MouseID,
    Recombinase.Group(color = :Protocol, linestyle = :Wall),
    select = (:Omissions_plus_one, :Leave),
    ribbon = true
    )
p2 = plot!(p1,args...;kwargs...,
    legend = false
    ,fillalpha = 0,linewidth=2)
##
savefig(p1,"/Volumes/GoogleDrive/My Drive/Reports for Zach/NewProtocols/Fitting/Survival_fit.pdf")
##
figures = "/Users/dariosarra/Documents/Lab/Mainen/Presentations/Walton_visit/figures/unprocessed_plots"
savefig(p1,joinpath(figures,"New_data_survival_fit.pdf"))
###
s = filter(streaks) do row
    row.Correct_start &&
    row.Streak > 10 &&
    3 < row.Streak_within_Block< 23 &&
    row.ReverseStreak > 20 &&
    (ismissing(row.Pre_Interpoke) || row.Pre_Interpoke < 1) &&
    row.Protocol_Day > 1 &&
    row.Day > 20191215
end
s[!,:AftrePlusOne] = s.AfterLast .+ 1
args, kwargs = series2D(
    Recombinase.discrete(Recombinase.density),
    s,
    error = :MouseID,
    Recombinase.Group(color = :Protocol, linestyle = :Wall),
    select = (:AftrePlusOne),
    ribbon = true
    )

p3 = plot!(p1,args...;kwargs...,legend = false,fillalpha = 0,linewidth=2)

plot!(p3,args...;kwargs...,)


text:
using ProgressMeter

x = zeros(10,10)

@time @showprogress 1 for i in 1:10000
        for j in 1:10000
            for k in 1:10
                for l in 1:10
                    x[k,l] = l
                end
            end
        end
    end

print()


text:
# Define BNFNode macros so we can record source locations in the BNF.

using MacroTools
using MacroTools: postwalk

find_name(s::Symbol) = s

function find_name(exp::Expr)
    @assert length(exp.args) >= 1
    find_name(exp.args[1])
end

function rewrite_struct(def::Expr)::Expr
    structpattern = postwalk(rmlines, Meta.parse(
        "struct nameexp_ fields__ end "))
    constructorpattern = postwalk(rmlines, Meta.parse(
        "function fname_(fargs__) fbody__ end"))
    bindings = MacroTools.match(structpattern, def)
    fields = []
    constructors = []
    for e in bindings[:fields]
        e = longdef(e)
        if isexpr(e, :function)
            push!(constructors, e)
        else
            push!(fields, e)
        end
    end
    if length(constructors) == 0
        # Add default constructor
        name = find_name(bindings[:nameexp])
        args = map(find_name, fields)
        push!(constructors,
              :(function $name($(args...))
                    new($(args...))
                end))
    end
    constructors2 = []
    for e in constructors
        # default source argument
        push!(constructors2,
              postwalk(e) do e
                  if iscall(e, :new)
                      Expr(:call, :new, nothing, e.args[2:end]...)
                  else
                      e
                  end
              end)
        # pass through source argument
        bindings2 = MacroTools.match(constructorpattern, e)
        push!(constructors2,
              :(function $(bindings2[:fname])(source::Union{Nothing, LineNumberNode},
                                              $(bindings2[:fargs]...))
                    $(map(bindings2[:fbody]) do e
                          postwalk(e) do e
                              if iscall(e, :new)
                                  Expr(:call, :new, :source, e.args[2:end]...)
                              else
                                  e
                              end
                          end
                      end...)
                end)        )
    end
    :(Base.@__doc__ struct $(bindings[:nameexp])
          source::Union{Nothing, LineNumberNode}
          $(fields...)

          $(constructors2...)
      end)
end

export @bnfnode

macro bnfnode(exp)
    @capture(exp,
             struct nameexp_
                 fields__
             end)
    name = find_name(nameexp)
    mname = Symbol("@$name")
    source = __source__
    @assert source isa LineNumberNode

    quote
        $(rewrite_struct(exp))

        export $(esc(name)), $(esc(mname))

        macro $(esc(name))(args...)
            Expr(:call, $(esc(name)), $(esc(source)), args...)
        end
    end
end


text:
using Test
using CUDA
using CUDAKernels
using KernelAbstractions

@kernel function copy_shared_kernel!(a,b)
    tid = @index(Global, Linear)

    #N = @uniform prod(@groupsize())

    #shared_tile = @localmem Int32 (N,)
    shared_tile = @localmem Int32 (10,)

    if tid <= length(a)
        @inbounds shared_tile[tid] = a[tid]
    end

    @synchronize()

    if tid <= length(a)
        @inbounds b[tid] = a[tid]
    end
end

function copy_shared!(a, b; numcores = 4, numthreads = 256)

    if isa(a, Array)
        kernel! = copy_shared_kernel!(CPU(), 4)
    else
        kernel! = copy_shared_kernel!(CUDADevice(), 256)
    end

    kernel!(a, b, ndrange=size(a))
end

@testset "CPU / GPU tests" begin
    a = Int64.(rand(1:128,10, 10))
    b = zeros(Int32, 10, 10)

    event = copy_shared!(a, b)
    wait(event)

    if has_cuda_gpu()
        d_a = CuArray(a)
        d_b = similar(d_a)

        event = copy_shared!(d_a, d_b)
        wait(event)

    end
end

max_stars_repo_name: leios/simuleios

max_issues_repo_name: leios/simuleios

max_forks_repo_name: leios/simuleios


text:
function lineplot(;x::AbstractVector = Int[],
                   y::AbstractVector = Int[],
                   group::AbstractVector = Int[],
                   points::Bool = false,
                   pointShape::AbstractString = "circle",
                   pointSize::Real = 10)

    v = VegaVisualization(name = "lineplot")

    add_data!(v, x = x, y = y, group = group)

    #Get unique table name
    table = v.data[1].name

    default_scales!(v)
    default_axes!(v)

    #If non-null array passed, put legend
    if group != Int[]
        legend!(v)
    end

    #Old add_lines! code
      innermarks = VegaMark(_type = "line",
                            properties = VegaMarkProperties(enter = VegaMarkPropertySet(x = VegaValueRef(scale = "x", field = "x"),
                                                                                y = VegaValueRef(scale = "y", field = "y"))
                                                                              )
                              )

      res = VegaMark(_type = "group",
                     from = VegaMarkFrom(data = table,
                                         transform =[VegaTransform(Dict{Any, Any}("type" => "facet", "groupby" => ["group"]))]
                                         ),
                     marks = [innermarks])

      v.marks == nothing? v.marks = [res] : push!(v.marks, res)

    #End add_lines!

    v.marks[1].marks[1].properties.enter.stroke = VegaValueRef(scale = "group", field="group")

    #Make line width a little thicker
    v.marks[1].marks[1].properties.enter.strokeWidth = VegaValueRef(value = 2)

    #Add points if desired
    if points

        copymark = deepcopy(v.marks[1].marks[1])
        push!(v.marks[1].marks, copymark)
        v.marks[1].marks[2]._type = "symbol"
        v.marks[1].marks[2].properties.enter.size = VegaValueRef(value = pointSize)
        v.marks[1].marks[2].properties.enter.shape = VegaValueRef(value = pointShape)

    end

    return v
end

max_stars_repo_name: jpthor/vega


text:
__precompile__(true)

module StatsFuns

using Base: Math.@horner
using Reexport
using SpecialFunctions
import ChainRulesCore
import InverseFunctions

# reexports
@reexport using IrrationalConstants:
    twoπ,       # 2π
    fourπ,      # 4π
    halfπ,      # π / 2
    quartπ,     # π / 4
    invπ,       # 1 / π
    twoinvπ,    # 2 / π
    fourinvπ,   # 4 / π
    inv2π,      # 1 / (2π)
    inv4π,      # 1 / (4π)
    sqrt2,      # √2
    sqrt3,      # √3
    sqrtπ,      # √π
    sqrt2π,     # √2π
    sqrt4π,     # √4π
    sqrthalfπ,  # √(π / 2)
    invsqrt2,   # 1 / √2
    invsqrt2π,  # 1 / √2π
    loghalf,    # log(1 / 2)
    logtwo,     # log(2)
    logπ,       # log(π)
    log2π,      # log(2π)
    log4π       # log(4π)

@reexport using LogExpFunctions:
    xlogx,          # x * log(x) for x > 0, or 0 when x == 0
    xlogy,          # x * log(y) for x > 0, or 0 when x == 0
    xlog1py,        # x * log(1 + y) for x > 0, or 0 when x == 0
    logistic,       # 1 / (1 + exp(-x))
    logit,          # log(x / (1 - x))
    log1psq,        # log(1 + x^2)
    log1pexp,       # log(1 + exp(x))
    log1mexp,       # log(1 - exp(x))
    log2mexp,       # log(2 - exp(x))
    logexpm1,       # log(exp(x) - 1)
    softplus,       # alias of log1pexp
    invsoftplus,    # alias of logexpm1
    log1pmx,        # log(1 + x) - x
    logmxp1,        # log(x) - x + 1
    logaddexp,      # log(exp(x) + exp(y))
    logsubexp,      # log(abs(e^x - e^y))
    logsumexp,      # log(sum(exp(x)))
    softmax,        # exp(x_i) / sum(exp(x)), for i
    softmax!        # inplace softmax

export
    # distrs/beta
    betapdf,            # pdf of beta distribution
    betalogpdf,         # logpdf of beta distribution
    betacdf,            # cdf of beta distribution
    betaccdf,           # ccdf of beta distribution
    betalogcdf,         # logcdf of beta distribution
    betalogccdf,        # logccdf of beta distribution
    betainvcdf,         # inverse-cdf of beta distribution
    betainvccdf,        # inverse-ccdf of beta distribution
    betainvlogcdf,      # inverse-logcdf of beta distribution
    betainvlogccdf,     # inverse-logccdf of beta distribution

    # distrs/binom
    binompdf,           # pdf of binomial distribution
    binomlogpdf,        # logpdf of binomial distribution
    binomcdf,           # cdf of binomial distribution
    binomccdf,          # ccdf of binomial distribution
    binomlogcdf,        # logcdf of binomial distribution
    binomlogccdf,       # logccdf of binomial distribution
    binominvcdf,        # inverse-cdf of binomial distribution
    binominvccdf,       # inverse-ccdf of binomial distribution
    binominvlogcdf,     # inverse-logcdf of binomial distribution
    binominvlogccdf,    # inverse-logccdf of binomial distribution

    # distrs/chisq
    chisqpdf,           # pdf of chi-square distribution
    chisqlogpdf,        # logpdf of chi-square distribution
    chisqcdf,           # cdf of chi-square distribution
    chisqccdf,          # ccdf of chi-square distribution
    chisqlogcdf,        # logcdf of chi-square distribution
    chisqlogccdf,       # logccdf of chi-square distribution
    chisqinvcdf,        # inverse-cdf of chi-square distribution
    chisqinvccdf,       # inverse-ccdf of chi-square distribution
    chisqinvlogcdf,     # inverse-logcdf of chi-square distribution
    chisqinvlogccdf,    # inverse-logccdf of chi-square distribution

    # distrs/fdist
    fdistpdf,           # pdf of F distribution
    fdistlogpdf,        # logpdf of F distribution
    fdistcdf,           # cdf of F distribution
    fdistccdf,          # ccdf of F distribution
    fdistlogcdf,        # logcdf of F distribution
    fdistlogccdf,       # logccdf of F distribution
    fdistinvcdf,        # inverse-cdf of F distribution
    fdistinvccdf,       # inverse-ccdf of F distribution
    fdistinvlogcdf,     # inverse-logcdf of F distribution
    fdistinvlogccdf,    # inverse-logccdf of F distribution

    # distrs/gamma
    gammapdf,           # pdf of gamma distribution
    gammalogpdf,        # logpdf of gamma distribution
    gammacdf,           # cdf of gamma distribution
    gammaccdf,          # ccdf of gamma distribution
    gammalogcdf,        # logcdf of gamma distribution
    gammalogccdf,       # logccdf of gamma distribution
    gammainvcdf,        # inverse-cdf of gamma distribution
    gammainvccdf,       # inverse-ccdf of gamma distribution
    gammainvlogcdf,     # inverse-logcdf of gamma distribution
    gammainvlogccdf,    # inverse-logccdf of gamma distribution

    # distrs/hyper
    hyperpdf,           # pdf of hypergeometric distribution
    hyperlogpdf,        # logpdf of hypergeometric distribution
    hypercdf,           # cdf of hypergeometric distribution
    hyperccdf,          # ccdf of hypergeometric distribution
    hyperlogcdf,        # logcdf of hypergeometric distribution
    hyperlogccdf,       # logccdf of hypergeometric distribution
    hyperinvcdf,        # inverse-cdf of hypergeometric distribution
    hyperinvccdf,       # inverse-ccdf of hypergeometric distribution
    hyperinvlogcdf,     # inverse-logcdf of hypergeometric distribution
    hyperinvlogccdf,    # inverse-logccdf of hypergeometric distribution

    # distrs/nbeta
    nbetapdf,           # pdf of noncentral beta distribution
    nbetalogpdf,        # logpdf of noncentral beta distribution
    nbetacdf,           # cdf of noncentral beta distribution
    nbetaccdf,          # ccdf of noncentral beta distribution
    nbetalogcdf,        # logcdf of noncentral beta distribution
    nbetalogccdf,       # logccdf of noncentral beta distribution
    nbetainvcdf,        # inverse-cdf of noncentral beta distribution
    nbetainvccdf,       # inverse-ccdf of noncentral beta distribution
    nbetainvlogcdf,     # inverse-logcdf of noncentral beta distribution
    nbetainvlogccdf,    # inverse-logccdf of noncentral beta distribution

    # distrs/nbinom
    nbinompdf,           # pdf of negative nbinomial distribution
    nbinomlogpdf,        # logpdf of negative nbinomial distribution
    nbinomcdf,           # cdf of negative nbinomial distribution
    nbinomccdf,          # ccdf of negative nbinomial distribution
    nbinomlogcdf,        # logcdf of negative nbinomial distribution
    nbinomlogccdf,       # logccdf of negative nbinomial distribution
    nbinominvcdf,        # inverse-cdf of negative nbinomial distribution
    nbinominvccdf,       # inverse-ccdf of negative nbinomial distribution
    nbinominvlogcdf,     # inverse-logcdf of negative nbinomial distribution
    nbinominvlogccdf,    # inverse-logccdf of negative nbinomial distribution

    # distrs/nchisq
    nchisqpdf,           # pdf of noncentral chi-square distribution
    nchisqlogpdf,        # logpdf of noncentral chi-square distribution
    nchisqcdf,           # cdf of noncentral chi-square distribution
    nchisqccdf,          # ccdf of noncentral chi-square distribution
    nchisqlogcdf,        # logcdf of noncentral chi-square distribution
    nchisqlogccdf,       # logccdf of noncentral chi-square distribution
    nchisqinvcdf,        # inverse-cdf of noncentral chi-square distribution
    nchisqinvccdf,       # inverse-ccdf of noncentral chi-square distribution
    nchisqinvlogcdf,     # inverse-logcdf of noncentral chi-square distribution
    nchisqinvlogccdf,    # inverse-logccdf of noncentral chi-square distribution

    # distrs/nfdist
    nfdistpdf,           # pdf of noncentral F distribution
    nfdistlogpdf,        # logpdf of noncentral F distribution
    nfdistcdf,           # cdf of noncentral F distribution
    nfdistccdf,          # ccdf of noncentral F distribution
    nfdistlogcdf,        # logcdf of noncentral F distribution
    nfdistlogccdf,       # logccdf of noncentral F distribution
    nfdistinvcdf,        # inverse-cdf of noncentral F distribution
    nfdistinvccdf,       # inverse-ccdf of noncentral F distribution
    nfdistinvlogcdf,     # inverse-logcdf of noncentral F distribution
    nfdistinvlogccdf,    # inverse-logccdf of noncentral F distribution

    # distrs/norm
    normpdf,            # pdf of normal distribution
    normlogpdf,         # logpdf of normal distribution
    normcdf,            # cdf of normal distribution
    normccdf,           # ccdf of normal distribution
    normlogcdf,         # logcdf of normal distribution
    normlogccdf,        # logccdf of normal distribution
    norminvcdf,         # inverse-cdf of normal distribution
    norminvccdf,        # inverse-ccdf of normal distribution
    norminvlogcdf,      # inverse-logcdf of normal distribution
    norminvlogccdf,     # inverse-logccdf of normal distribution

    # distrs/ntdist
    ntdistpdf,           # pdf of noncentral t distribution
    ntdistlogpdf,        # logpdf of noncentral t distribution
    ntdistcdf,           # cdf of noncentral t distribution
    ntdistccdf,          # ccdf of noncentral t distribution
    ntdistlogcdf,        # logcdf of noncentral t distribution
    ntdistlogccdf,       # logccdf of noncentral t distribution
    ntdistinvcdf,        # inverse-cdf of noncentral t distribution
    ntdistinvccdf,       # inverse-ccdf of noncentral t distribution
    ntdistinvlogcdf,     # inverse-logcdf of noncentral t distribution
    ntdistinvlogccdf,    # inverse-logccdf of noncentral t distribution

    # distrs/pois
    poispdf,           # pdf of Poisson distribution
    poislogpdf,        # logpdf of Poisson distribution
    poiscdf,           # cdf of Poisson distribution
    poisccdf,          # ccdf of Poisson distribution
    poislogcdf,        # logcdf of Poisson distribution
    poislogccdf,       # logccdf of Poisson distribution
    poisinvcdf,        # inverse-cdf of Poisson distribution
    poisinvccdf,       # inverse-ccdf of Poisson distribution
    poisinvlogcdf,     # inverse-logcdf of Poisson distribution
    poisinvlogccdf,    # inverse-logccdf of Poisson distribution

    # distrs/tdist
    tdistpdf,           # pdf of student's t distribution
    tdistlogpdf,        # logpdf of student's t distribution
    tdistcdf,           # cdf of student's t distribution
    tdistccdf,          # ccdf of student's t distribution
    tdistlogcdf,        # logcdf of student's t distribution
    tdistlogccdf,       # logccdf of student's t distribution
    tdistinvcdf,        # inverse-cdf of student's t distribution
    tdistinvccdf,       # inverse-ccdf of student's t distribution
    tdistinvlogcdf,     # inverse-logcdf of student's t distribution
    tdistinvlogccdf,    # inverse-logccdf of student's t distribution

    # distrs/srdist
    srdistcdf,           # cdf of studentized range distribution
    srdistccdf,          # ccdf of studentized range distribution
    srdistlogcdf,        # logcdf of studentized range distribution
    srdistlogccdf,       # logccdf of studentized range distribution
    srdistinvcdf,        # inverse-cdf of studentized range distribution
    srdistinvccdf,       # inverse-ccdf of studentized range distribution
    srdistinvlogcdf,     # inverse-logcdf of studentized range distribution
    srdistinvlogccdf,    # inverse-logccdf of studentized range distribution

    # misc
    logmvgamma,         # logarithm of multivariate gamma function
    logmvbeta,          # logarithm of multivariate beta function
    lstirling_asym

## source files
include("misc.jl")
include("rmath.jl")

using .RFunctions

include(joinpath("distrs", "beta.jl"))
include(joinpath("distrs", "binom.jl"))
include(joinpath("distrs", "chisq.jl"))
include(joinpath("distrs", "fdist.jl"))
include(joinpath("distrs", "gamma.jl"))
include(joinpath("distrs", "hyper.jl"))
include(joinpath("distrs", "nbeta.jl"))
include(joinpath("distrs", "nbinom.jl"))
include(joinpath("distrs", "nchisq.jl"))
include(joinpath("distrs", "nfdist.jl"))
include(joinpath("distrs", "norm.jl"))
include(joinpath("distrs", "ntdist.jl"))
include(joinpath("distrs", "pois.jl"))
include(joinpath("distrs", "tdist.jl"))
include(joinpath("distrs", "srdist.jl"))

include("chainrules.jl")
include("inverse.jl")

end # module


text:
"Thrown at @expect failures in step definitions."
struct StepAssertFailure <: Exception
    assertion::String
end

"""
    expect(ex)

Assert that a condition `ex` is true. Throws a StepAssertFailure if false.

# Examples
```
@then "this condition should hold" begin
    @expect 1 == 1
end
```

This will fail and throw a `StepAssertFailure`
```
@then "this condition should not hold" begin
    @expect 1 == 2
end
```
"""
macro expect(ex)
    # Create a human readable form of the expectation.
    if length(ex.args) > 0 && ex.args[1] isa Expr && ex.args[1].head == :globalref
        operator = ex.args[1].args[2]
        readable_ex = copy(ex)
        readable_ex.args[1] = operator
    else
        readable_ex = ex
    end

    quote
        # Check the expectation.
        if !($(esc(ex)))
            throw(StepAssertFailure($(string(readable_ex))))
        end
    end
end


text:
module CTFM

#Use this to print
#=
weave(Pkg.dir("$(pwd())\\FamaMacbeth.jl"),
  informat="script",
  out_path = "$(pwd())\\FamaMacbethOUT.html",
  doctype = "md2html")
=#

if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, StatsBase, GZip, JLD,
  Gadfly, CTModCopy, JuMP, NLopt, HCubature, StaticArrays, CTNumerical,
  ParallelAccelerator, Measures

  #=const WORKERS = 3 #seems like logical cores * (3/4) is the magic value
  const CLEAN_WORKERS = false

  if nworkers() < WORKERS || (nworkers()==1 && WORKERS == 1)
    addprocs(nworkers()==1?WORKERS-nworkers()+1:WORKERS-nworkers())
  end

  @everywhere using DataFrames, Distributions, StatsBase, GZip, JLD,
    Gadfly, CTModCopy, JuMP, NLopt, HCubature, StaticArrays, CTNumerical=#

  const constDef = true
  const DATA_PATH = pwd() * "\\data" #path where we store the data
  const OUTPUT_PATH = pwd() * "\\output"
  const DATE_FORMAT_STR = "yyyymmdd"
  const SP500_NAME = "GSPC"
  const TBILL_NAME = "t-bill"
  const FACTOR_NAME = "FMFactors"
  const ASSET_NAME = "PORTSBM1010"
  #const ASSET_NAME = "PORTSBM23"
  #const ASSET_NAME = "PORTPI55"
  const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
  const HEADER_NAME = "header.tex"

  const PlotContainer = Union{Plot,Gadfly.Compose.Context}
  const CSRegSpec = Tuple{Symbol, Symbol, Int}
  const CPRegSpec = Tuple{Symbol, Int}

  const FACTOR_SYMS = [:Mkt_RF, :SMB, :HML]
  #const FACTOR_SYMS = [:Mkt_RF]
  const RFR_SYM = :RF

  const T_LIMIT = 99999
  #const T_LIMIT = 100

  BLAS.set_num_threads(1) #need this since we are already multi-threading at a higher level

  #pre-process the asset data
  function preProcessFA()::Void

    #read in the data (note the non-standard NA values)
    dateFormat::String = DATE_FORMAT_STR
    assetDF::DataFrame = readtable("$DATA_PATH\\$ASSET_NAME.csv", nastrings=["-999", "-99.99"])
    assetDF[:,:DATE] = ((i::Int)->Date("$i", dateFormat)::Date).(assetDF[:,:DATE])::DataVector{Date}

    factorDF::DataFrame = readtable("$DATA_PATH\\$FACTOR_NAME.csv", nastrings=["-999", "-99.99"])
    factorDF[:,:DATE] = ((i::Int)->Date("$i", dateFormat)::Date).(factorDF[:,:DATE])::DataVector{Date}
    factorDF = factorDF[:,[:DATE; FACTOR_SYMS; RFR_SYM]]

    #merge the dataframes
    DF::DataFrame = join(factorDF, assetDF, on = :DATE, kind=:inner)

    #write out a JLS binary file
    stream::IOStream = open("$DATA_PATH\\$(FACTOR_NAME)_$(ASSET_NAME).jls", "w")
    serialize(stream, DF)
    close(stream)
    return nothing
  end

struct CTFMWindow
  FADF::DataFrame
  FSyms::Vector{Symbol}
  RSyms::Vector{Symbol}

  β::Matrix{Float64} #(NxK)
  α::DataMatrix{Float64} #(NxT)
  λ::Matrix{Float64} #(KxT)

  T::Int #T time slices
  N::Int #N assets
  K::Int #k factors
end

#basic constructor to get coefficients
function CTFMWindow(FADF::DataFrame, FSyms::Vector{Symbol}, RSyms::Vector{Symbol})

  BLAS.set_num_threads(1) #we will do our own multi-threading here
  N::Int = length(RSyms)
  K::Int = length(FSyms)
  T::Int = min(size(FADF,1), T_LIMIT)
  FADF = FADF[(end - T + 1):(end),:]

  #println("size: ",size(FADF,1))

  #prepare the first stage
  #form the RHS
  FExpr::CTExpr = parse(join((String).(FSyms), " + "))
  RHSNames::Vector{Symbol} = [:intercept; FSyms]

  #run the N 1st stage regressions
  β = Matrix{Float64}(N,K)
  @fastmath Threads.@threads for i::Int ∈ 1:N
    β[i,:] .= ((CTLM(FADF, FExpr, RSyms[i], XNames=RHSNames)).β)[2:end]
  end

  #prepare the 2nd stage
  βDF = DataFrame(β)

  names!(βDF, FSyms)
  TSyms::Vector{Symbol} = ((i::Int)->Symbol("T$(i)")).(1:T)

  #transpose the R Matrix and append to the right of the factor results
  FADFRows = eachrow(FADF[:,RSyms])
  vToTranspose = DataVector{Float64}(Vector{Float64}(N))
  for t::Int ∈ 1:T
    βDF[:,TSyms[t]] = DataVector{Float64}(Vector{Float64}(N)) #allocate space
    for i::Int ∈ 1:N # copy the entries 1 by 1
      vToTranspose[i] = FADFRows[t][i]
    end
    βDF[:,TSyms[t]] .= vToTranspose
  end

    #showall(βDF)

  α::DataArray{Float64,2} = DataArray{Float64,2}(Matrix{Float64}(N,T))
  α .= NA

  #run the 2nd stage regressions
  cλ::Matrix{Float64} = Matrix{Float64}(K+1,T)
  @fastmath Threads.@threads for t::Int ∈ 1:T
    cλ[:,t] .= CTLM(βDF,FExpr,TSyms[t],XNames=RHSNames).β
  end

  @fastmath Threads.@threads for t::Int ∈ 1:T
    α[:,t] = βDF[:,TSyms[t]] .- cλ[1,t]
  end

  α -= β * cλ[2:end,:]

#=println("cλ1: ", cλ[1,1])
println("βDF[:,TSyms[1]]: ", βDF[:,TSyms[1]])
println("β * cλ[2:end,t]:", β * cλ[2:end,1])
println("α1: ", α[:,1])=#

  #store the resutls in an object and return it
  return CTFMWindow(FADF, FSyms, RSyms, β, α, cλ[2:end,:], T, N, K)
end

#provides the mean risk premia
meanλ(fm::CTFMWindow)::Vector{Float64} = vec(mean(fm.λ,2))

#get the standard errors of the estimates
function σλ(fm::CTFMWindow)::Vector{Float64}
  λ₀::Vector{Float64} = meanλ(fm)
  return ((1.0/fm.T^2) .* ((k::Int)->sum((fm.λ[k,:]-λ₀[k]).^2)).(1:fm.K)).^0.5
end

#provides the mean risk premia across time
αi(fm::CTFMWindow)::Vector{Float64} = ((i::Int)->mean(dropna(fm.α[i,:]))).(1:fm.N)

#performs the cochrane pricing error test and returns (χ²,p)
function χ²α(fm::CTFMWindow)::Tuple{Float64,Float64}
  α₀::Vector{Float64} = αi(fm)
  #BLAS.set_num_threads(8)

  #first build the covariance matrix

  #pre-allcoate the space
  αtMα₀::DataVector{Float64} = DataVector{Float64}(Vector{Float64}(fm.N))
  covα::Matrix{Float64} = zeros(fm.N, fm.N)
  covαt::DataArray{Float64,2} = DataArray{Float64,2}(Matrix{Float64}(fm.N,fm.N)) #temporary for holding NAs
  TMat::Matrix{Int} = zeros(Int,fm.N, fm.N) #counts the number of non-missing values
  NAMat::BitArray{2} = trues(fm.N, fm.N)

  for t::Int ∈ 1:fm.T
    αtMα₀ .= fm.α[:,t] .- α₀
    covαt .= αtMα₀ * (αtMα₀)'
    NAMat .= (isna).(covαt)
    TMat .+= 1 .- NAMat
    covαt[NAMat] .= 0.0
    covα .+= covαt
  end

  covα .= covα ./ (TMat.^2)

  #calculate the inverse via the cholesky process
  covαInv::Matrix{Float64} = covα \ eye(fm.N)
  χ²::Float64 = α₀' * covαInv * α₀
  p::Float64 = cdf(Chisq(fm.N-1), χ²)
  return χ², p
end

#get the standard errors of the pricing errors
function σα(fm::CTFMWindow)::Vector{Float64}
  αi₀::Vector{Float64} = αi(fm)
  σ::Vector{Float64} = Vector{Float64}(fm.N)
  for i::Int ∈ 1:fm.N
    T::Float64 = length(fm.α[i,:])
    σ[i] = (sum((dropna(fm.α[i,:]).-αi₀[i]).^2) / T^2)^0.5
  end
  return σ #((1.0/fm.T^2) .* ((i::Int)->sum((dropna(fm.α[i,:]).-αi₀[i]).^2)).(1:fm.N)).^0.5
end

function getDF(pathStr::String, refreshData::Bool, preProcessFunc::Function)::DataFrame

  #re-process the data if needed or desired
  if !isfile(pathStr) || refreshData
    preProcessFunc()
  end

  #read the binary file
  stream::IOStream = open(pathStr)
  DF::DataFrame = deserialize(stream)
  close(stream)

  return DF
end

#an encapsulating function for making the graphs
function makeFMScriptGraphs(fm::CTFMWindow, label::Symbol; sorts::Tuple{Symbol, Symbol}=(:ME, :BM),
    sortDims::Tuple{Int,Int} = (10,10))::Void
  pathStr::String = "$OUTPUT_PATH\\$(FACTOR_NAME)_$(ASSET_NAME)_$label"
  FSyms::Vector{Symbol} = FACTOR_SYMS
  RFRSym::Symbol = RFR_SYM

  #this structure allows us to get a vector of symbols corresponding to a single
  #dimension of the sort. So LinkedRSyms[1][4][3] corresponds to ME3
  formToGrid(v::Vector)::Matrix = reshape(v, sortDims[2], sortDims[1])

  RSymMat::Matrix{Symbol} = formToGrid(fm.RSyms)
  RSymDict::Dict = Dict(fm.RSyms[i] => i for i::Int ∈ 1:length(fm.RSyms))

  #turn the betas into a matrix
  βMat::NTuple{fm.K, Matrix{Float64}} = (((k::Int)->Matrix{Float64}(sortDims[2], sortDims[1])).(1:(fm.K))...)
  for k = 1:fm.K
    for j = 1:sortDims[1]
      for i = 1:sortDims[2]
        βMat[k][i,j] = fm.β[RSymDict[RSymMat[i,j]],k]
      end
    end
  end

  Gadfly.push_theme(:default)
  plots::Vector{PlotContainer} = Vector{PlotContainer}()
  plotNames::Vector{String} = Vector{String}()

  push!(plotNames, "SMB1st")
  push!(plots, plot(x=collect(1:10), y=vec(mean(βMat[2], 1)),
    Geom.line,Geom.point,
    Guide.ylabel("β (SMB)"), Guide.xlabel("SMB Portfolio Decile"),
    Coord.Cartesian(ymin=-0.5,ymax=1.5, xmin=0, xmax=10),
    Guide.title("First Stage SMB Beta by SMB Decile")))

  push!(plotNames, "HML1st")
  push!(plots, plot(x=collect(1:10), y=vec(mean(βMat[3], 2)),
    Geom.line,Geom.point,
    Guide.ylabel("β (HML)"), Guide.xlabel("HML Portfolio Decile"),
    Coord.Cartesian(ymin=-0.5,ymax=1.5, xmin=0, xmax=10),
    Guide.title("First Stage HML Beta by HML Decile")))

  push!(plotNames, "SMB2nd")
  push!(plots, plot(x=collect(1:10), y=vec(mean(βMat[2], 2)),
    Geom.line,Geom.point,
    Guide.ylabel("β (SMB)"), Guide.xlabel("HML Portfolio Decile"),
    Coord.Cartesian(ymin=-0.5,ymax=1.5, xmin=0, xmax=10),
    Guide.title("First Stage SMB Beta by HML Decile")))

  push!(plotNames, "HML2nd")
  push!(plots, plot(x=collect(1:10), y=vec(mean(βMat[3], 1)),
    Geom.line,Geom.point,
    Guide.ylabel("β (HML)"), Guide.xlabel("SMB Portfolio Decile"),
    Coord.Cartesian(ymin=-0.5,ymax=1.5, xmin=0, xmax=10),
    Guide.title("First Stage HML Beta by SMB Decile")))

  push!(plotNames, "MRF1")
  push!(plots, plot(x=collect(1:10), y=vec(mean(βMat[1], 1)),
    Geom.line,Geom.point,
    Guide.ylabel("β (HML)"), Guide.xlabel("SMB Portfolio Decile"),
    Coord.Cartesian(ymin=-0.5,ymax=1.5, xmin=0, xmax=10),
    Guide.title("First Stage M-Rf Beta by SMB Decile")))

  push!(plotNames, "MRF2")
  push!(plots, plot(x=collect(1:10), y=vec(mean(βMat[1], 1)),
    Geom.line,Geom.point,
    Guide.ylabel("β (HML)"), Guide.xlabel("HML Portfolio Decile"),
    Coord.Cartesian(ymin=-0.5,ymax=1.5, xmin=0, xmax=10),
    Guide.title("First Stage M-Rf Beta by HML Decile")))

  push!(plotNames, "ErrorHist")
  push!(plots, plot(x=σα(fm)./100*255^.5,
    Geom.histogram(bincount=20,density=true),
    Guide.ylabel("Freq"), Guide.xlabel("σ (ann.)"),
    #Coord.Cartesian(ymin=-0.5,ymax=1.5, xmin=0, xmax=10),
    Guide.title("Frequency of Time Series Avg Pricing Errors")))

    push!(plotNames, "ErrorHistT")
    push!(plots, plot(x=((αi(fm)/100+1.0).^255-1.0) ./ (σα(fm)./100*255^.5),
      Geom.histogram(bincount=20,density=true),
      Guide.ylabel("Freq"), Guide.xlabel("T"),
      #Coord.Cartesian(ymin=-0.5,ymax=1.5, xmin=0, xmax=10),
      Guide.title("Frequency of Avg T Statistics")))

  for i ∈ 1:length(plots)
    draw(SVG("$(pathStr)_$(plotNames[i]).svg", 4inch, 4inch), plots[i])
  end

  #βMat::NTuple{fm.K, Matrix{Float64}} = (((k::Int)->formToGrid(fm.β[:,k])).(1:(fm.K))...)

  return nothing
end

function A7Script(; refreshData::Bool = true, rerunReg::Bool = true,
    label::Symbol=Symbol(""), calcχ²::Bool=true)::Void
  FADF::DataFrame = getDF("$DATA_PATH\\$(FACTOR_NAME)_$(ASSET_NAME).jls", refreshData, preProcessFA)

  #rename non-standard column names
  rename!(FADF, [:SMALL_LoBM, :SMALL_HiBM, :BIG_LoBM, :BIG_HiBM],
    [:ME1_BM1, :ME1_BM10, :ME10_BM1, :ME10_BM10])
    #=rename!(FADF, [:LoOP_LoINV, :LoOP_HiINV, :HiOP_LoINV, :HiOP_HiINV],
      [:OP1_INV1, :OP1_INV10, :OP10_INV1, :OP10_INV10])=#
  #build the vector of return symbols
  #showcols(FADF)
  #setup the factor symbols
  FSyms::Vector{Symbol} = FACTOR_SYMS
  RFRSym::Symbol = RFR_SYM

  RSyms::Vector{Symbol} = names(FADF[:,setdiff(names(FADF), [:DATE; RFR_SYM; FACTOR_SYMS])])

  #net out the risk-free rate
  for s::Symbol ∈ RSyms
    FADF[:,s] = FADF[:,s] .- FADF[:,RFRSym]
  end

  #run the FM regressions
  if rerunReg
    gstream::GZipStream = gzopen("$DATA_PATH\\REG_OUT$label.jls.gz", "w")
    serialize(gstream, CTFMWindow(FADF, FSyms, RSyms))
    close(gstream)
  end

  gstream = gzopen("$DATA_PATH\\REG_OUT$label.jls.gz")
  fm::CTFMWindow = deserialize(gstream)
  close(gstream)

  makeFMScriptGraphs(fm, label)

  λ₀::Vector{Float64} = meanλ(fm)
  σλ₀::Vector{Float64} = σλ(fm)

  println(fm.FSyms, "\n", "λ(Ann.): $((1.0.+λ₀/100).^255-1)", "\n", "λ: $λ₀",
    "\n", "σ(λ): $σλ₀", "\n", "T(λ): $(λ₀./σλ₀)")

  if calcχ²
    println("\nCalculating χ² statistics")
    χ²::Float64, p::Float64 = χ²α(fm)
    println("χ²($(fm.N-1))): $χ², p=$p")
  end

  return nothing
end

#uncomment to run code
#=@time begin
  A7Script(refreshData=false, rerunReg=false, calcχ²=true, label=Symbol("$(Dates.today())"))
  #cleanworkers()
end

gc()=#

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
function getdefaultoptparameters(L)
    tol_converged = 1e-3
    tol_𝓧_spacing = 1e-6
    tol_derivative_zero = 0.2
    tol_no_update = 1e-6
    #tol_h = 1e-17
    N_samples_interval = 10
    N_test_positions = 200
    max_iters = 200
    L_interpolant = 4

    tol_converged = convert(BigFloat, tol_converged )
    tol_𝓧_spacing = convert(BigFloat, tol_𝓧_spacing/(L+2))
    tol_derivative_zero = convert(BigFloat, tol_derivative_zero )
    tol_no_update = convert(BigFloat, tol_no_update )
    #tol_h = convert(BigFloat, tol_h )
    N_samples_interval = convert(Int,N_samples_interval)
    N_test_positions = convert(Int,N_test_positions)
    max_iters = convert(Int,max_iters)
    L_interpolant = convert(Int,L_interpolant)
    verbose_flag = true
    plot_flag = false

    config = MiniMaxConfigType( tol_converged,
                                tol_𝓧_spacing,
                                tol_derivative_zero,
                                tol_no_update,
                                #tol_h,
                                N_samples_interval,
                                N_test_positions,
                                max_iters,
                                L_interpolant,
                                verbose_flag,
                                plot_flag)

    return config
end

function getoptparameters(L,file_name::String)

    tol_converged,
    tol_𝓧_spacing,
    tol_derivative_zero,
    tol_no_update,
    #tol_h,
    N_samples_interval,
    N_test_positions,
    max_iters,
    L_interpolant,
    verbose_flag_input,
    plot_flag_input = DelimitedFiles.readdlm(file_name, ',', Float64)

    tol_converged = convert(BigFloat, tol_converged )
    tol_𝓧_spacing = convert(BigFloat, tol_𝓧_spacing/(L+2))
    tol_derivative_zero = convert(BigFloat, tol_derivative_zero )
    tol_no_update = convert(BigFloat, tol_no_update )
    #tol_h = convert(BigFloat, tol_h )
    N_samples_interval = convert(Int,N_samples_interval)
    N_test_positions = convert(Int,N_test_positions)
    max_iters = convert(Int,max_iters)
    L_interpolant = convert(Int,L_interpolant)
    verbose_flag = verbose_flag_input==1 ? true : false
    plot_flag = plot_flag_input==1 ? true : false

    config = MiniMaxConfigType( tol_converged,
                                tol_𝓧_spacing,
                                tol_derivative_zero,
                                tol_no_update,
                                #tol_h,
                                N_samples_interval,
                                N_test_positions,
                                max_iters,
                                L_interpolant,
                                verbose_flag,
                                plot_flag)
    return config
end

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
#==============================================================================#
# PI.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/pi-2018-02-27.normal.json
#==============================================================================#

__precompile__()

module PI

using AWSCore

"""
    using AWSSDK.PI.describe_dimension_keys
    describe_dimension_keys([::AWSConfig], arguments::Dict)
    describe_dimension_keys([::AWSConfig]; ServiceType=, Identifier=, StartTime=, EndTime=, Metric=, GroupBy=, <keyword arguments>)

    using AWSCore.Services.pi
    pi([::AWSConfig], "DescribeDimensionKeys", arguments::Dict)
    pi([::AWSConfig], "DescribeDimensionKeys", ServiceType=, Identifier=, StartTime=, EndTime=, Metric=, GroupBy=, <keyword arguments>)

# DescribeDimensionKeys Operation

For a specific time period, retrieve the top `N` dimension keys for a metric.

# Arguments

## `ServiceType = "RDS"` -- *Required*
The AWS service for which Performance Insights will return metrics. The only valid value for *ServiceType* is: `RDS`

## `Identifier = ::String` -- *Required*
An immutable, AWS Region-unique identifier for a data source. Performance Insights gathers metrics from this data source.

To use an Amazon RDS instance as a data source, you specify its `DbiResourceId` value - for example: `db-FAIHNTYBKTGAUSUZQYPDS2GW4A`

## `StartTime = timestamp` -- *Required*
The date and time specifying the beginning of the requested time series data. You can't specify a `StartTime` that's earlier than 7 days ago. The value specified is *inclusive* - data points equal to or greater than `StartTime` will be returned.

The value for `StartTime` must be earlier than the value for `EndTime`.

## `EndTime = timestamp` -- *Required*
The date and time specifying the end of the requested time series data. The value specified is *exclusive* - data points less than (but not equal to) `EndTime` will be returned.

The value for `EndTime` must be later than the value for `StartTime`.

## `Metric = ::String` -- *Required*
The name of a Performance Insights metric to be measured.

Valid values for `Metric` are:

*   `db.load.avg` - a scaled representation of the number of active sessions for the database engine.

*   `db.sampledload.avg` - the raw number of active sessions for the database engine.

## `PeriodInSeconds = ::Int`
The granularity, in seconds, of the data points returned from Performance Insights. A period can be as short as one second, or as long as one day (86400 seconds). Valid values are:

*   `1` (one second)

*   `60` (one minute)

*   `300` (five minutes)

*   `3600` (one hour)

*   `86400` (twenty-four hours)

If you don't specify `PeriodInSeconds`, then Performance Insights will choose a value for you, with a goal of returning roughly 100-200 data points in the response.

## `GroupBy = [ ... ]` -- *Required*
A specification for how to aggregate the data points from a query result. You must specify a valid dimension group. Performance Insights will return all of the dimensions within that group, unless you provide the names of specific dimensions within that group. You can also request that Performance Insights return a limited number of values for a dimension.
```
 GroupBy = [
        "Group" => <required> ::String,
        "Dimensions" =>  [::String, ...],
        "Limit" =>  ::Int
    ]
```

## `PartitionBy = [ ... ]`
For each dimension specified in `GroupBy`, specify a secondary dimension to further subdivide the partition keys in the response.
```
 PartitionBy = [
        "Group" => <required> ::String,
        "Dimensions" =>  [::String, ...],
        "Limit" =>  ::Int
    ]
```

## `Filter = ::Dict{String,String}`
One or more filters to apply in the request. Restrictions:

*   Any number of filters by the same dimension, as specified in the `GroupBy` or `Partition` parameters.

*   A single filter for any other dimension in this dimension group.

## `MaxResults = ::Int`
The maximum number of items to return in the response. If more items exist than the specified `MaxRecords` value, a pagination token is included in the response so that the remaining results can be retrieved.

## `NextToken = ::String`
An optional pagination token provided by a previous request. If this parameter is specified, the response includes only records beyond the token, up to the value specified by `MaxRecords`.

# Returns

`DescribeDimensionKeysResponse`

# Exceptions

`InvalidArgumentException`, `InternalServiceError` or `NotAuthorizedException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/pi-2018-02-27/DescribeDimensionKeys)
"""
@inline describe_dimension_keys(aws::AWSConfig=default_aws_config(); args...) = describe_dimension_keys(aws, args)

@inline describe_dimension_keys(aws::AWSConfig, args) = AWSCore.Services.pi(aws, "DescribeDimensionKeys", args)

@inline describe_dimension_keys(args) = describe_dimension_keys(default_aws_config(), args)

"""
    using AWSSDK.PI.get_resource_metrics
    get_resource_metrics([::AWSConfig], arguments::Dict)
    get_resource_metrics([::AWSConfig]; ServiceType=, Identifier=, MetricQueries=, StartTime=, EndTime=, <keyword arguments>)

    using AWSCore.Services.pi
    pi([::AWSConfig], "GetResourceMetrics", arguments::Dict)
    pi([::AWSConfig], "GetResourceMetrics", ServiceType=, Identifier=, MetricQueries=, StartTime=, EndTime=, <keyword arguments>)

# GetResourceMetrics Operation

Retrieve Performance Insights metrics for a set of data sources, over a time period. You can provide specific dimension groups and dimensions, and provide aggregation and filtering criteria for each group.

# Arguments

## `ServiceType = "RDS"` -- *Required*
The AWS service for which Performance Insights will return metrics. The only valid value for *ServiceType* is: `RDS`

## `Identifier = ::String` -- *Required*
An immutable, AWS Region-unique identifier for a data source. Performance Insights gathers metrics from this data source.

To use an Amazon RDS instance as a data source, you specify its `DbiResourceId` value - for example: `db-FAIHNTYBKTGAUSUZQYPDS2GW4A`

## `MetricQueries = [[ ... ], ...]` -- *Required*
An array of one or more queries to perform. Each query must specify a Performance Insights metric, and can optionally specify aggregation and filtering criteria.
```
 MetricQueries = [[
        "Metric" => <required> ::String,
        "GroupBy" =>  [
            "Group" => <required> ::String,
            "Dimensions" =>  [::String, ...],
            "Limit" =>  ::Int
        ],
        "Filter" =>  ::Dict{String,String}
    ], ...]
```

## `StartTime = timestamp` -- *Required*
The date and time specifying the beginning of the requested time series data. You can't specify a `StartTime` that's earlier than 7 days ago. The value specified is *inclusive* - data points equal to or greater than `StartTime` will be returned.

The value for `StartTime` must be earlier than the value for `EndTime`.

## `EndTime = timestamp` -- *Required*
The date and time specifiying the end of the requested time series data. The value specified is *exclusive* - data points less than (but not equal to) `EndTime` will be returned.

The value for `EndTime` must be later than the value for `StartTime`.

## `PeriodInSeconds = ::Int`
The granularity, in seconds, of the data points returned from Performance Insights. A period can be as short as one second, or as long as one day (86400 seconds). Valid values are:

*   `1` (one second)

*   `60` (one minute)

*   `300` (five minutes)

*   `3600` (one hour)

*   `86400` (twenty-four hours)

If you don't specify `PeriodInSeconds`, then Performance Insights will choose a value for you, with a goal of returning roughly 100-200 data points in the response.

## `MaxResults = ::Int`
The maximum number of items to return in the response. If more items exist than the specified `MaxRecords` value, a pagination token is included in the response so that the remaining results can be retrieved.

## `NextToken = ::String`
An optional pagination token provided by a previous request. If this parameter is specified, the response includes only records beyond the token, up to the value specified by `MaxRecords`.

# Returns

`GetResourceMetricsResponse`

# Exceptions

`InvalidArgumentException`, `InternalServiceError` or `NotAuthorizedException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/pi-2018-02-27/GetResourceMetrics)
"""
@inline get_resource_metrics(aws::AWSConfig=default_aws_config(); args...) = get_resource_metrics(aws, args)

@inline get_resource_metrics(aws::AWSConfig, args) = AWSCore.Services.pi(aws, "GetResourceMetrics", args)

@inline get_resource_metrics(args) = get_resource_metrics(default_aws_config(), args)

end # module PI

#==============================================================================#
# End of file
#==============================================================================#


text:
#############
# Management: getting info, deleting containers etc
import Base.stat
export list, delete

stat(serv::PyObject) = serv[:stat]() |> Dict
stat(serv::PyObject, container::AbstractString) = serv[:stat](container) |> Dict

function stat(serv::PyObject, container::AbstractString, object::String)
    serv[:stat](container, [object]) |> first
end

function list(serv::PyObject)
	response = serv[:list]() |> first |> Dict
	Dict(d["name"] => d for d in response["listing"])
end

list(serv::PyObject, container::AbstractString) = serv[:list](container) |> Dict

### Delete

delete(serv::PyObject, container::AbstractString) = serv[:delete](container) |> collect

function delete(serv::PyObject, container::AbstractString, object::String)
    serv[:delete](container, [object]) |> first
end


text:
# using MLJ, EvoTrees, MLJScientificTypes
using MLJ, EvoTrees, CSV, DataFrames

num_cols = [
    "ClientPeriod",
    "MonthlySpending",
    "TotalSpent"
];

cat_cols = [
    "Sex",
    "IsSeniorCitizen",
    "HasPartner",
    "HasChild",
    "HasPhoneService",
    "HasMultiplePhoneNumbers",
    "HasInternetService",
    "HasOnlineSecurityService",
    "HasOnlineBackup",
    "HasDeviceProtection",
    "HasTechSupportAccess",
    "HasOnlineTV",
    "HasMovieSubscription",
    "HasContractPhone",
    "IsBillingPaperless",
    "PaymentMethod"
];

all_feature_cols = [num_cols; cat_cols];
target_col = "Churn";

#,types=Dict("Sex"=>CategoricalValue{String, UInt32})
# + tags=[]
df = DataFrame!(CSV.File("./data/train.csv",pool=0.1, missingstrings=[" "]))
categorical!(df,[cat_cols;target_col]);
describe(df,:eltype,:nunique, :nmissing)
dropmissing!(df);
describe(df,:eltype,:nunique, :nmissing)

X = df[!, all_feature_cols];
y = df[!,target_col];

mach_x = machine(ContinuousEncoder(), X)
fit!(mach_x)
X = MLJ.transform(mach_x, X)

X_mat = Array(X)
train, test = partition(eachindex(y), 0.7, shuffle=true); # 70:30 split
x_train, y_train = X_mat[train, :], y[train]
x_test, y_test = X_mat[test, :], y[test]

tree_model = EvoTreeClassifier(T=Float32, max_depth=9, nrounds=200, η = 0.2, colsample=0.3, metric=:mlogloss)
@time model = fit_evotree(tree_model, x_train, y_train, print_every_n = 10);
pred_test = EvoTrees.predict(model, x_test)
minimum(pred_test)
maximum(pred_test)

tree_model = EvoTreeClassifier(T=Float32, max_depth=6, nrounds=200,colsample=0.3)
mach = machine(tree_model, X, y)
train, test = partition(eachindex(y), 0.7, shuffle=true); # 70:30 split
fit!(mach, rows=train, verbosity=1)
pred_test = MLJ.predict(mach, selectrows(X, test))

y_train_num = Float64.([y.level for y in y_train]) .- 1
tree_model = EvoTreeRegressor(T=Float32, loss=:logistic, metric=:logloss, max_depth=9, nrounds=200, η=0.1, colsample=0.5)
@time model = fit_evotree(tree_model, x_train, y_train_num, print_every_n = 25);
pred_test = EvoTrees.predict(model, x_test)
minimum(pred_test)
maximum(pred_test)

y_train_num = Float64.([y.level for y in y_train]) .- 1
tree_model = EvoTreeRegressor(T=Float32, loss=:logistic, metric=:logloss, max_depth=9, nrounds=200, η=0.1, colsample=0.5, device="gpu")
@time model = fit_evotree(tree_model, x_train, y_train_num, print_every_n = 20);
pred_test = EvoTrees.predict(model, x_test)
minimum(pred_test)
maximum(pred_test)

y_train_num = Float64.([y.level for y in y_train]) .- 1
tree_model = EvoTreeRegressor(T=Float32, loss=:linear, metric=:mse, max_depth=9, nrounds=200, η=0.1, colsample=0.5, device="gpu")
@time model = fit_evotree(tree_model, x_train, y_train_num, print_every_n = 20);
pred_test = EvoTrees.predict(model, x_test)
minimum(pred_test)
maximum(pred_test)


text:
# TODO:
# Shared fitness cannot be used until a better distance calculation scheme has
# been implemented. Fitness should not depend on "Individual". Distance schemes
# should be interchangeable.

type SharedFitness{T}
  base::T
  shared::Float

  SharedFitness(f::T) = new(f)
end

type FitnessSharingScheme <: FitnessScheme
  base::FitnessScheme
  radius::Float
  alpha::Float
  dist::Distance

  FitnessSharingScheme(b::FitnessScheme, r::Float, a::Float, d::Distance) =
    new(b, r, a, d)
end

"""
TODO: Document shared fitness schemes.

"""
shared(s::Dict{Any, Any}) =
  SharedFitnessScheme(s["base"], s["radius"], s["alpha"], s["distance"])

uses(s::FitnessSharingScheme) = SharedFitness{uses(s.base)}
maximise(s::FitnessSharingScheme) = maximise(s.base)
assign(s::FitnessSharingScheme, args...) =
  fitness(s.base, args...)

sh(s::FitnessSharingScheme, d::Float) =
  d < s.radius ? 1.0 - (d/s.radius) ^ s.alpha : 0.0

scale!{F}(s::FitnessSharingScheme, inds::Vector{Tuple{Int, F}}) =
  for (i1_id, i1_f) in inds
    i1_f.shared = score(i1_f.base) /
      sum((i2_id, i2_f) -> sh(s, distance(s.distance, i1, i2)), inds)
  end

compare(s::FitnessSharingScheme, x::SharedFitness, y::SharedFitness) =
  if x.shared > y.shared
    maximise(s) ? -1 : 1
  elseif x.shared < y.shared
    maximise(s) ? 1 : -1
  else
    0
  end


text:
type SubProduct{M <: MarketModelMultiProduct}
  product::M
  multiplier::Float64
  numberOfCashflows::Vector{Int}
  cashFlows::Vector{Vector{MarketModelCashFlow}}
  timeIndices::Vector{Int}
  isDone::Bool
end

# product is already cloned
SubProduct(product::MarketModelMultiProduct, multiplier::Float64, isDone::Bool) =
          SubProduct(product, multiplier, Vector{Int}(), Vector{Vector{MarketModelCashFlow}}(), Vector{Int}(), isDone)

type MarketModelComposite <: MarketModelMultiProduct
  components::Vector{SubProduct}
  rateTimes::Vector{Float64}
  evolutionTimes::Vector{Float64}
  finalized::Bool
  currentIndex::Int
  cashFlowTimes::Vector{Float64}
  allEvolutionTimes::Vector{Vector{Float64}}
  isInSubset::Vector{BitArray{1}}
  evolution::EvolutionDescription
end

MarketModelComposite() = MarketModelComposite(Vector{SubProduct}(), Vector{Float64}(), Vector{Float64}(), false, 1, Vector{Float64}(), Vector{Vector{Float64}}(),
                            Vector{BitArray{1}}(), EvolutionDescription())

function add_product!(mm::MarketModelComposite, product::MarketModelMultiProduct, multiplier::Float64 = 1.0)
  product = clone(product)
  d = get_evolution(product)
  if ~isempty(mm.components)
    # enforce pre-conditions
    d1 = get_evolution(mm.components[1].product)
    rateTimes1 = d1.rateTimes
    rateTimes2 = d.rateTimes
    (length(rateTimes1) == length(rateTimes2) && rateTimes1 == rateTimes2) || error("incompatible rate times")
  end

  push!(mm.components, SubProduct(product, multiplier, false))
  push!(mm.allEvolutionTimes, d.evolutionTimes)

  return mm
end

function finalize!(mm::MarketModelComposite)
  ~mm.finalized || error("already finalized")
  ~isempty(mm.components) || error("no subproducts provided")

  # fetch the rate times from the first subproduct (we already checked that they are the same)
  description = get_evolution(mm.components[1].product)
  mm.rateTimes = description.rateTimes

  mm.evolutionTimes, mm.isInSubset = merge_times(mm.allEvolutionTimes)

  allCashFlowTimes = Vector{Float64}()

  # now for each subproduct
  for i in eachindex(mm.components)
    d = get_evolution(mm.components[i].product)

    # collect all possible cash flow times
    cashFlowTimes = possible_cash_flow_times(mm.components[i].product)
    append!(allCashFlowTimes, cashFlowTimes)

    # allocate working vectors
    n = number_of_products(mm.components[i].product)
    mm.components[i].numberOfCashflows = Vector{Int}(n)
    mm.components[i].cashFlows = Vector{MarketModelCashFlow}[
                                [MarketModelCashFlow() for i=1:max_number_of_cashflows_per_product_per_step(mm.components[i].product)] for j = 1:n]
  end

  # all information collected, we can sort and compact the vector of all cashflow times
  sort!(allCashFlowTimes)
  allCashFlowTimes = unique(allCashFlowTimes)

  append!(mm.cashFlowTimes, allCashFlowTimes)

  # and map each product's cash flow time into the total vector
  for i in eachindex(mm.components)
    productTimes = possible_cash_flow_times(mm.components[i].product)
    mm.components[i].timeIndices = Int[findfirst(mm.cashFlowTimes, productTimes[j]) for j in eachindex(productTimes)]
  end

  mm.evolution = EvolutionDescription(mm.rateTimes, mm.evolutionTimes)

  mm.finalized = true

  return mm
end

function next_time_step!(mm::MarketModelComposite, currentState::CurveState, numberCashFlowsThisStep::Vector{Int}, cashFlowsGenerated::Vector{Vector{MarketModelCashFlow}})
  mm.finalized || error("composite not finalized")

  isDone = true
  n = 1
  offset = 0
  for i in eachindex(mm.components)
    if mm.isInSubset[n][mm.currentIndex] && ~mm.components[i].isDone
      # make it evolve
      thisDone = next_time_step!(mm.components[i].product, currentState, mm.components[i].numberOfCashflows, mm.components[i].cashFlows)

      # and copy the results. Time indices need to be remapped so that they point into all cash-flow times
      # Amounts need to be adjusted by the corresponding multiplier
      for j = 1:number_of_products(mm.components[i].product)
        numberCashFlowsThisStep[j+offset] = mm.components[i].numberOfCashflows[j]
        for k = 1:mm.components[i].numberOfCashflows[j]
          cfFrom = mm.components[i].cashFlows[j][k]
          cfTo = cashFlowsGenerated[j+offset][k]
          cfTo.timeIndex = mm.components[i].timeIndices[cfFrom.timeIndex]
          cfTo.amount = cfFrom.amount * mm.components[i].multiplier
        end
      end

      # finally set isDone to false if this product isn't done
      isDone = isDone && thisDone
    else
      for j = 1:number_of_products(mm.components[i].product)
        numberCashFlowsThisStep[j+offset] = 0
      end
    end

    # the offset is updated whether or not the product was evolved
    offset += number_of_products(mm.components[i].product)
    n += 1
  end
  mm.currentIndex += 1
  return isDone
end

function number_of_products(mm::MarketModelComposite)
  res = 0
  for i in eachindex(mm.components)
    res += number_of_products(mm.components[i].product)
  end

  return res
end

function max_number_of_cashflows_per_product_per_step(mm::MarketModelComposite)
  res = 0
  for i in eachindex(mm.components)
    res  = max(res, max_number_of_cashflows_per_product_per_step(mm.components[i].product))
  end

  return res
end

function possible_cash_flow_times(mm::MarketModelComposite)
  mm.finalized || error("composite not finalized")
  return mm.cashFlowTimes
end

function get_evolution(mm::MarketModelComposite)
  mm.finalized || error("composite not finalized")
  return mm.evolution
end

function reset!(mm::MarketModelComposite)
  for i in eachindex(mm.components)
    reset!(mm.components[i].product)
    mm.components[i].isDone = false
  end
  mm.currentIndex = 1
  return mm
end

clone(mm::MarketModelComposite) = MarketModelComposite(deepcopy(mm.components), copy(mm.rateTimes), copy(mm.evolutionTimes), mm.finalized, mm.currentIndex,
                                    copy(mm.cashFlowTimes), deepcopy(mm.allEvolutionTimes), deepcopy(mm.isInSubset), clone(mm.evolution))


text:
module Bodies

import Base: parent

using AstroTime: SECONDS_PER_DAY, SECONDS_PER_CENTURY, value, j2000, seconds
using ItemGraphs: ItemGraph, SimpleGraph, add_edge!, add_vertex!, items

export CelestialBody,
    NAIFId,
    SolarSystemBarycenter,
    Sun,
    along_orbit_radius,
    declination,
    declination_rate,
    ellipsoid,
    equatorial_radius,
    euler_angles,
    euler_rates,
    from_naifid,
    grav_param,
    mean_radius,
    naifid,
    parent,
    polar_radius,
    right_ascension,
    right_ascension_rate,
    rotation_angle,
    rotation_rate,
    ssb,
    subplanetary_radius,
    sun

const NAIFId = Int
const BODIES = ItemGraph{NAIFId, NAIFId}(SimpleGraph())

abstract type CelestialBody end

function grav_param end
function mean_radius end
function polar_radius end
function equatorial_radius end
function subplanetary_radius end
function along_orbit_radius end

function ellipsoid(body::CelestialBody)
    subplanetary_radius(body), along_orbit_radius(body), polar_radius(body)
end

Base.show(io::IO, body::CelestialBody) = print(io, string(nameof(typeof(body))))

from_naifid(id::NAIFId) = from_naifid(Val(id))

path_ids(from::CelestialBody, to::CelestialBody) = items(BODIES, naifid(from), naifid(to))

abstract type Barycenter <: CelestialBody end

struct SolarSystemBarycenter <: Barycenter end
const ssb = SolarSystemBarycenter()
Base.show(io::IO, ::SolarSystemBarycenter) = print(io, "Solar System Barycenter")
naifid(::SolarSystemBarycenter) = 0
from_naifid(::Val{0}) = ssb
add_vertex!(BODIES, 0)

struct Sun <: CelestialBody end
const sun = Sun()
parent(::Sun) = ssb
naifid(::Sun) = 10
from_naifid(::Val{10}) = sun
add_edge!(BODIES, 0, 10)

alpha0(::CelestialBody) = 0.0
alpha1(::CelestialBody) = 0.0
alpha2(::CelestialBody) = 0.0
delta0(::CelestialBody) = 0.0
delta1(::CelestialBody) = 0.0
delta2(::CelestialBody) = 0.0
omega0(::CelestialBody) = 0.0
omega1(::CelestialBody) = 0.0
omega2(::CelestialBody) = 0.0
alpha(::CelestialBody) = zeros(1)
delta(::CelestialBody) = zeros(1)
omega(::CelestialBody) = zeros(1)
theta0(::CelestialBody) = zeros(1)
theta1(::CelestialBody) = zeros(1)

const Δtc = SECONDS_PER_CENTURY
const Δtd = SECONDS_PER_DAY

theta(b::CelestialBody, t) = theta0(b) .+ theta1(b) .* t / Δtc

function right_ascension(b::CelestialBody, ep)
    t = value(seconds(j2000(ep)))
    α = alpha(b)
    α₀ = alpha0(b)
    α₁ = alpha1(b)
    α₂ = alpha2(b)
    θ = theta(b, t)
    α₀ + α₁ * t / Δtc + α₂ * t^2 / Δtc^2 + sum(α .* sin.(θ))
end

function declination(b::CelestialBody, ep)
    t = value(seconds(j2000(ep)))
    δ = delta(b)
    δ₀ = delta0(b)
    δ₁ = delta1(b)
    δ₂ = delta2(b)
    θ = theta(b, t)
    δ₀ + δ₁ * t / Δtc + δ₂ * t^2 / Δtc^2 + sum(δ .* cos.(θ))
end

function rotation_angle(b::CelestialBody, ep)
    t = value(seconds(j2000(ep)))
    ω = omega(b)
    ω₀ = omega0(b)
    ω₁ = omega1(b)
    ω₂ = omega2(b)
    θ = theta(b, t)
    ω₀ + ω₁ * t / Δtd + ω₂ * t^2 / Δtd^2 + sum(ω .* sin.(θ))
end

function right_ascension_rate(b::CelestialBody, ep)
    t = value(seconds(j2000(ep)))
    α = alpha(b)
    α₁ = alpha1(b)
    α₂ = alpha2(b)
    θ = theta(b, t)
    θ₁ = theta1(b)
    α₁ / Δtc + 2 * α₂ * t / Δtc^2 + sum(α .* θ₁ ./ Δtc .* cos.(θ))
end

function declination_rate(b::CelestialBody, ep)
    t = value(seconds(j2000(ep)))
    δ = delta(b)
    δ₁ = delta1(b)
    δ₂ = delta2(b)
    θ = theta(b, t)
    θ₁ = theta1(b)
    δ₁ / Δtc + 2 * δ₂ * t / Δtc^2 - sum(δ .* θ₁ ./ Δtc .* sin.(θ))
end

function rotation_rate(b::CelestialBody, ep)
    t = value(seconds(j2000(ep)))
    ω = omega(b)
    ω₁ = omega1(b)
    ω₂ = omega2(b)
    θ = theta(b, t)
    θ₁ = theta1(b)
    ω₁ / Δtd + 2 * ω₂ * t / Δtd^2 + sum(ω .* θ₁ ./ Δtc .* cos.(θ))
end

function euler_angles(b::CelestialBody, ep)
    right_ascension(b, ep) + π/2, π/2 - declination(b, ep), mod2pi(rotation_angle(b, ep))
end

function euler_rates(b::CelestialBody, ep)
    right_ascension_rate(b, ep), -declination_rate(b, ep), rotation_rate(b, ep)
end

include("planets.jl")
include("minor.jl")
include("satellites.jl")

const ALL_NAMES = Tuple([["Sun", "SolarSystemBarycenter"]
                         collect(PLANET_NAMES .* "Barycenter");
                         collect(PLANET_NAMES);
                         ["Luna", "Phobos", "Deimos"];
                         collect(JUPITER_SATELLITE_NAMES);
                         collect(SATURN_SATELLITE_NAMES);
                         collect(URANUS_SATELLITE_NAMES);
                         collect(NEPTUNE_SATELLITE_NAMES);
                         collect(MINOR_BODY_NAMES);
                         ["Pluto", "PlutoBarycenter"];
                         collect(PLUTO_SATELLITE_NAMES)])

gm = joinpath(@__DIR__, "..", "..", "gen", "gm.jl")
isfile(gm) || error("`gm.jl` has not been generated, yet. Run `gen/gen_gm.jl`")
include(gm)

pck = joinpath(@__DIR__, "..", "..", "gen", "pck.jl")
isfile(pck) || error("`pck.jl` has not been generated, yet. Run `gen/gen_pck.jl`")
include(pck)

end