text:
#===============================================================================

        Defines the main workhorse functions of this package that perform
        MCMC sampling for stochastic processes observed discretely in time.

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

"""
    run!(
        mcmc::MCMC, num_mcmc_steps, data, θinit, callbacks=Callback[];
        kwargs...
    )

The main calling function of this package that initializes, runs and outputs the
results of an MCMC sampler for discretely observed stochastic processes.

# Arguments
- `mcmc::MCMC`: specifies a sequence of updates that constitute a single MCMC step
                and provides additional info about the MCMC backend algorithm to
                be used
- `num_mcmc_steps`: the total number of MCMC steps
- `data`: completely characterizes everything that is known about the
          observations (including everything known about the underlying dynamics
          and the way it was collected),
- `θinit`: an initial guess for the unknown parameters
- `callbacks`: a collection of extra utility functions that do extra work around
               MCMC sampling.
- `exclude_updates`: a standard named argument which lists the update indices
                     and corresponding ranges of MCMC iterations at which given
                     updates are supposed to be omitted.
- `kwargs...`: additional named arguments, passed onto initializers of global
               workspace; they depend on the chosen `MCMCBackend`
"""
function run!(
        mcmc::MCMC, num_mcmc_steps, data, θinit, callbacks=Callback[];
        kwargs...
    )
    init!(
        mcmc,
        num_mcmc_steps,
        data,
        θinit,
        get(kwargs, :exclude_updates, []);
        kwargs...
    )
    local_wss = create_workspaces(mcmc.backend, mcmc)

    init!(callbacks, mcmc.workspace)
    __run!(mcmc.workspace, local_wss, mcmc.updates, mcmc.schedule, callbacks)
    for callback in callbacks
        cleanup!(callback, mcmc.workspace, local_wss, (mcmciter=num_mcmc_steps,))
    end
    mcmc.workspace, local_wss
end

"""
    __run!(global_ws, local_wss, updates, schedule, callbacks)

Internal loops that run the (already initialized) MCMC sampler. `global_ws` and
`local_wss` are the already initialized global workspace and local workspaces
(one for each update) respectively. `schedule` is the iterator over MCMC steps
and `callbacks` is a list of callbacks.
"""
function __run!(global_ws, local_wss, updates, schedule, callbacks)
    #TODO find out whether try-catch is the best way to register user's
    # interruptions for initiating an early exit from the loop.
    #TODO implement functionality that registers user interruption,
    # prompts for confirmation of early exit and whether to save the most recent
    # state.
    for step in schedule
        local_ws, local_update = update_workspaces!(
            updates[step.pidx],
            global_ws,
            local_wss[step.pidx],
            step,
            (step.prev_pidx===nothing ? nothing : local_wss[step.prev_pidx])
        )
        update_callbacks!(callbacks, global_ws, local_wss, step, __PRESTEP)
        update!(local_update, global_ws, local_ws, step)
        update_adaptation!(updates, global_ws, local_ws, step)
        update_callbacks!(callbacks, global_ws, local_wss, step, __POSTSTEP)
    end
end

#            RELATED TO TRANSFER OF INFORMATION BETWEEN WORKSPACES
#-------------------------------------------------------------------------------

"""
    update_workspaces!(
        local_updt::MCMCParamUpdate,
        global_ws::GlobalWorkspace,
        local_ws::LocalWorkspace,
        step,
        prev_ws,
    )
Transfer all the information that is needed at the time of starting the MCMC
update step from the global (or a previously used local) to the currently used
local workspace.
"""
function update_workspaces!(
        local_updt::MCMCParamUpdate,
        global_ws::GlobalWorkspace,
        local_ws::LocalWorkspace,
        step,
        prev_ws,
    )
    state(local_ws) .= state(global_ws, local_updt) #local_ws.sub_ws.state on LHS
    prev_ws == nothing || ( ll(local_ws) .= ll(prev_ws, step.prev_mcmciter) )
    compute_gradients_and_momenta!(local_updt, local_ws, __PREVIOUS)
    local_ws, local_updt
end

#=
function transfer_local_knowledge!(
        local_updt::MCMCParamUpdate,
        local_ws::LocalWorkspace,
        step,
        ::Nothing,
    )
end

function transfer_local_knowledge!(
        local_updt::MCMCParamUpdate,
        local_ws::LocalWorkspace,
        step,
        prev_ws,
    )
    ll(local_ws) .= ll(prev_ws, step.prev_mcmciter)
end
=#

#                          RELATED TO ADAPTATION
#-------------------------------------------------------------------------------

function update_adaptation!(
        updates::AbstractArray{<:MCMCUpdate},
        global_ws::GlobalWorkspace,
        local_ws::LocalWorkspace,
        step
    )
    # there may be a list, you must make sure that the first entry only is used
    # for param updates.
    _accepted = accepted(local_ws, step.mcmciter)

    for (i, updt) in enumerate(updates)
        update_adaptation!(_accepted, updt, global_ws, step, i)
    end
end

#=
function update_adaptation!(
        accepted::Bool, updt::MCMCUpdate, global_ws, step, i
    )
    nothing
end
=#
function update_adaptation!(
        accepted::T, updt::MCMCUpdate, global_ws, step, i
    ) where T <: Union{Bool, AbstractArray{<:Bool}}
    typeof(updt.adpt) <: NoAdaptation && return
    _my_turn = Val(step.pidx==i)
    register_only_on_my_turn(_my_turn, updt.adpt) && return
    register!(updt, updt.adpt, accepted, state(global_ws, updt))
    ttu = time_to_update(_my_turn, updt.adpt)
    #ttu && println("acceptance rate: ", acceptance_rate(updt.adpt))
    #ttu && println("old ϵ: ", updt.rw.ϵ)
    #ttu && println("acceptance rate: ", acceptance_rate(updt.adpt.v[1]))
    #ttu && println("old ρ: ", updt.ρs[1][1])
    ttu && readjust!(updt, updt.adpt, step.mcmciter)
    #ttu && println("new ρ: ", updt.ρs[1][1])
    #ttu && println()
end

# overwrite for your custom adaptations if need be
register_only_on_my_turn(::Val{true}, args...) = false
register_only_on_my_turn(::Val{false}, args...) = true
time_to_update(::Val{false}, args...) = false

#           RELATED TO UPDATE (see the remaining methods in updates.jl)
#-------------------------------------------------------------------------------
#NOTE This is likely to change once conjugate and gradient-based updates are
# implemented and some missing features are discovered
"""
    update!(
        updt::MCMCParamUpdate,
        global_ws::GlobalWorkspace,
        ws::LocalWorkspace,
        step,
    )

General recipe for performing parameter update (doing an update of the main
`state` of the MCMC sampler).
"""
function update!(
        updt::MCMCParamUpdate,
        global_ws::GlobalWorkspace,
        ws::LocalWorkspace,
        step,
    )
    proposal!(updt, global_ws, ws, step) # specific to updt, see `updates.jl`
    set_proposal!(updt, global_ws, ws, step)
    compute_ll!(updt, global_ws, ws, step)
    accept_reject!(updt, global_ws, ws, step)
    update_stats!(global_ws, ws, step)
end

"""
    set_proposal!(
        updt::MCMCParamUpdate,
        global_ws::GlobalWorkspace,
        ws::LocalWorkspace,
        step
    )

The proposal parameter θ° has been sampled by this point. This function
propagates it through workspaces and sets the new parameter in proposal law `P°`
"""
function set_proposal!(
        updt::MCMCParamUpdate,
        global_ws::GlobalWorkspace,
        local_ws::LocalWorkspace,
        step
    )
    _set_local_state°(updt, global_ws, local_ws, step)
    set_parameters!(__PROPOSAL, updt, global_ws, local_ws)
end

function _set_local_state°(
        updt::MCMCParamUpdate,
        global_ws::GlobalWorkspace,
        local_ws::LocalWorkspace,
        step
    )
    θ° = state°(global_ws, step)
    θ° .= state(global_ws)
    subidx(θ°, updt) .= state°(local_ws)
end

"""
    compute_ll!(
        global_ws::GlobalWorkspace,
        ws::LocalWorkspace,
        step,
    )

Evaluate the proposal log-likelihood at the observations.
"""
function compute_ll!(
        updt::MCMCParamUpdate,
        global_ws::GlobalWorkspace,
        local_ws::LocalWorkspace,
        step,
    )
    ll°(local_ws) .= loglikelihood(global_ws, __PROPOSAL)
    ll°(local_ws, step.mcmciter) .= ll°(local_ws)
    compute_gradients_and_momenta!(updt, local_ws, __PROPOSAL)
end

"""
    accept_reject!(updt, global_ws, ws, step, i=1)

Finish computations of the log-likelihood ratio between the target and proposal
in MH acceptance probability and accept/reject respectively.
"""
function accept_reject!(
        updt, global_ws, ws, step, _ll=sum(ll(ws)), _ll°=sum(ll°(ws)), i=1
    )
    llr = (
        _ll° - _ll
        + log_transition_density(__PROPOSAL, updt, ws, i)
        - log_transition_density(__PREVIOUS, updt, ws, i)
        + log_prior(__PROPOSAL, updt, ws, i)
        - log_prior(__PREVIOUS, updt, ws, i)
    )
    accepted = rand(Exponential(1.0)) > -llr
    #println("i: $(step.mcmciter), ll: $(_ll), ll°: $(_ll°), llr: $(llr), a/r: ", (accepted ? "✔" : "✘"))
    register_accept_reject_results!(accepted, updt, global_ws, ws, step, i)
end

#NOTE by default accept conjugate updates?
function accept_reject!(updt::MCMCConjugateParamUpdate, global_ws, ws, step)
    register_accept_reject_results!(true, updt, global_ws, ws, step, 1)
end

"""
    register_accept_reject_results!(
        accepted::Bool,
        updt,
        global_ws::GlobalWorkspace,
        local_ws::LocalWorkspace,
        step,
        i=1
    )
Register the result of accept/reject step.
"""
function register_accept_reject_results!(
        accepted::Bool,
        updt,
        global_ws::GlobalWorkspace,
        local_ws::LocalWorkspace,
        step,
        i=1
    )
    register_accept_reject_results!(accepted, updt, local_ws, step, i)
    set_chain_param!(accepted, updt, global_ws, local_ws, step)
    set_parameters!(__PREVIOUS, updt, global_ws, local_ws)
end

function set_chain_param!(accepted, updt, global_ws, local_ws, step)
    θ = state(global_ws)
    accepted && (
        # update current local and global states
        state(local_ws) .= state°(local_ws);
        subidx(θ, updt) .= state°(local_ws)
    )
    state(global_ws, step) .= θ
end

"""
    register_accept_reject_results!(
        accepted::Bool, updt, ws::LocalWorkspace, step, i=1
    )

Register the results of accept/reject step relevant to a local workspace.
"""
function register_accept_reject_results!(
        accepted::Bool, updt, ws::LocalWorkspace, step, i=1
    )
    ll°(ws, step.mcmciter)[i] = ll°(ws)[i]
    ll(ws, step.mcmciter)[i] = ( accepted ? ll°(ws)[i] : ll(ws)[i] )
    set_accepted!(ws, step.mcmciter, accepted)
end

"""
    log_transition_density(
        ::Previous, updt::MCMCParamUpdate, ws::LocalWorkspace, i
    )

Evaluate the log-density for a transition θ → θ°.
"""
function log_transition_density(
        ::Previous,
        updt::MCMCParamUpdate,
        ws::LocalWorkspace,
        i
    )
    log_transition_density(updt, state(ws), state°(ws))
end

"""
    log_transition_density(
        ::Proposal, updt::MCMCParamUpdate, ws::LocalWorkspace, i
    )

Evaluate the log-density for a transition θ° → θ.
"""
function log_transition_density(
        ::Proposal,
        updt::MCMCParamUpdate,
        ws::LocalWorkspace,
        i
    )
    log_transition_density(updt, state°(ws), state(ws))
end

"""
    log_prior(::Previous, updt::MCMCParamUpdate, ws::LocalWorkspace, i)

Evaluate the log-prior at θ.
"""
function log_prior(::Previous, updt::MCMCParamUpdate, ws::LocalWorkspace, i)
    log_prior(updt, state(ws))
end

"""
    log_prior(::Proposal, updt::MCMCParamUpdate, ws::LocalWorkspace, i)

Evaluate the log-prior at θ°.
"""
function log_prior(::Proposal, updt::MCMCParamUpdate, ws::LocalWorkspace, i)
    log_prior(updt, state°(ws))
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct ScenarioBasedInternal <: ForecastInternal
        label::String
        resolution::Dates.Period
        initial_time::Dates.DateTime
        scenario_count::Int64
        time_series_uuid::UUIDs.UUID
        horizon::Int
        internal::InfrastructureSystemsInternal
    end

A Discrete Scenario Based forecast for a particular data field in a Component.

# Arguments
- `label::String`: label of component parameter forecasted
- `resolution::Dates.Period`
- `initial_time::Dates.DateTime`: forecast availability time
- `scenario_count::Int64`: Number of scenarios
- `time_series_uuid::UUIDs.UUID`: reference to time series data; timestamp - scalingfactor
- `horizon::Int`: length of this forecast
- `internal::InfrastructureSystemsInternal`
"""
mutable struct ScenarioBasedInternal <: ForecastInternal
    "label of component parameter forecasted"
    label::String
    resolution::Dates.Period
    "forecast availability time"
    initial_time::Dates.DateTime
    "Number of scenarios"
    scenario_count::Int64
    "reference to time series data; timestamp - scalingfactor"
    time_series_uuid::UUIDs.UUID
    "length of this forecast"
    horizon::Int
    internal::InfrastructureSystemsInternal
end

function ScenarioBasedInternal(label, resolution, initial_time, scenario_count, time_series_uuid, horizon, )
    ScenarioBasedInternal(label, resolution, initial_time, scenario_count, time_series_uuid, horizon, InfrastructureSystemsInternal(), )
end

function ScenarioBasedInternal(; label, resolution, initial_time, scenario_count, time_series_uuid, horizon, )
    ScenarioBasedInternal(label, resolution, initial_time, scenario_count, time_series_uuid, horizon, )
end

"""Get ScenarioBasedInternal label."""
get_label(value::ScenarioBasedInternal) = value.label
"""Get ScenarioBasedInternal resolution."""
get_resolution(value::ScenarioBasedInternal) = value.resolution
"""Get ScenarioBasedInternal initial_time."""
get_initial_time(value::ScenarioBasedInternal) = value.initial_time
"""Get ScenarioBasedInternal scenario_count."""
get_scenario_count(value::ScenarioBasedInternal) = value.scenario_count
"""Get ScenarioBasedInternal time_series_uuid."""
get_time_series_uuid(value::ScenarioBasedInternal) = value.time_series_uuid
"""Get ScenarioBasedInternal horizon."""
get_horizon(value::ScenarioBasedInternal) = value.horizon
"""Get ScenarioBasedInternal internal."""
get_internal(value::ScenarioBasedInternal) = value.internal

"""Set ScenarioBasedInternal label."""
set_label!(value::ScenarioBasedInternal, val::String) = value.label = val
"""Set ScenarioBasedInternal resolution."""
set_resolution!(value::ScenarioBasedInternal, val::Dates.Period) = value.resolution = val
"""Set ScenarioBasedInternal initial_time."""
set_initial_time!(value::ScenarioBasedInternal, val::Dates.DateTime) = value.initial_time = val
"""Set ScenarioBasedInternal scenario_count."""
set_scenario_count!(value::ScenarioBasedInternal, val::Int64) = value.scenario_count = val
"""Set ScenarioBasedInternal time_series_uuid."""
set_time_series_uuid!(value::ScenarioBasedInternal, val::UUIDs.UUID) = value.time_series_uuid = val
"""Set ScenarioBasedInternal horizon."""
set_horizon!(value::ScenarioBasedInternal, val::Int) = value.horizon = val
"""Set ScenarioBasedInternal internal."""
set_internal!(value::ScenarioBasedInternal, val::InfrastructureSystemsInternal) = value.internal = val

ext: jl


text:
module FlowNodeModelsFlowModelTest

using Test
using ReactiveMP
using ReactiveMP: getforward, getbackward, getjacobian, getinv_jacobian
using ReactiveMP: forward, forward!, backward, backward!, jacobian, jacobian!, inv_jacobian, inv_jacobian!, forward_jacobian, backward_inv_jacobian
using ReactiveMP: det_jacobian, absdet_jacobian, logdet_jacobian, logabsdet_jacobian, detinv_jacobian, absdetinv_jacobian, logabsdetinv_jacobian
@testset "Flow Model" begin

    @testset "Constructor" begin

        # check for single layer
        f = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        model = FlowModel( 2, (layer, ) )
        compiled_model = compile(model)
        @test length(compiled_model.layers)      == 1
        @test typeof(compiled_model.layers[1])   <: AdditiveCouplingLayer
        @test typeof(compiled_model.layers[1].f[1]) <: PlanarFlow

        # check for two layers
        f1 = PlanarFlow()
        layer1 = AdditiveCouplingLayer(f1)
        f2 = PlanarFlow()
        layer2 = AdditiveCouplingLayer(f2; permute=false)
        model = FlowModel( 2, (layer1, layer2) )
        compiled_model = compile(model)
        @test length(compiled_model.layers)      == 3
        @test typeof(compiled_model.layers[1])   <: AdditiveCouplingLayer
        @test typeof(compiled_model.layers[1].f[1]) <: PlanarFlow
        @test typeof(compiled_model.layers[2])   <: PermutationLayer
        @test typeof(compiled_model.layers[3])   <: AdditiveCouplingLayer
        @test typeof(compiled_model.layers[3].f[1]) <: PlanarFlow

    end

    @testset "Get" begin

        # check get functions for single layer model
        f = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        model = FlowModel(2, (layer, ) )
        compiled_model = compile(model)
        @test getlayers(compiled_model)               == compiled_model.layers
        @test length(getlayers(compiled_model))       == 1
        @test typeof(getforward(compiled_model))      <: Function
        @test typeof(getbackward(compiled_model))     <: Function
        @test typeof(getjacobian(compiled_model))     <: Function
        @test typeof(getinv_jacobian(compiled_model)) <: Function

        # check get functions for multi layer model
        f1 = PlanarFlow()
        layer1 = AdditiveCouplingLayer(f1)
        f2 = PlanarFlow()
        layer2 = AdditiveCouplingLayer(f2; permute=false)
        model = FlowModel( 2, (layer1, layer2) )
        compiled_model = compile(model)
        @test getlayers(compiled_model)               == compiled_model.layers
        @test length(getlayers(compiled_model))       == 3
        @test typeof(getforward(compiled_model))      <: Function
        @test typeof(getbackward(compiled_model))     <: Function
        @test typeof(getjacobian(compiled_model))     <: Function
        @test typeof(getinv_jacobian(compiled_model)) <: Function

    end

    @testset "Base" begin

        # check base functions (single layer)
        f = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        model = FlowModel(2,  (layer, ) )
        compiled_model = compile(model)
        @test eltype(compiled_model) == Float64
        @test length(compiled_model) == 1

        # check base functions (multi layer)
        f1 = PlanarFlow()
        layer1 = AdditiveCouplingLayer(f1;)
        f2 = PlanarFlow()
        layer2 = AdditiveCouplingLayer(f2; permute=false)
        model = FlowModel( 2, (layer1, layer2) )
        compiled_model = compile(model)
        @test eltype(compiled_model) == Float64
        @test length(compiled_model) == 3

    end

    @testset "Forward-Backward" begin

        # check forward function (single layer)
        params = [1.0, 2.0, -3.0]
        f = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        model = FlowModel( 2, (layer, ) )
        compiled_model = compile(model, params)
        @test forward(compiled_model, [5.0, 1.5]) == [5.0, 7.4999983369439445]
        @test forward(compiled_model, [4.0, 2.5]) == [4.0, 7.499909204262595]
        @test forward.(compiled_model, [[5.0, 1.5], [4.0, 2.5]]) == [[5.0, 7.4999983369439445], [4.0, 7.499909204262595] ]

        # check forward function (multiple layers)
        params = [1.0, 2.0, -3.0, 1.0, 2.0, -3.0]
        f1 = PlanarFlow()
        layer1 = AdditiveCouplingLayer(f1)
        f2 = PlanarFlow()
        layer2 = AdditiveCouplingLayer(f2; permute=false)
        model = FlowModel( 2, (layer1, layer2) )
        compiled_model = compile(model, params)
        @test forward(compiled_model, [5.0, 1.5]) == [7.4999983369439445, 13.49999833686844]
        @test forward(compiled_model, [4.0, 2.5]) == [7.499909204262595, 12.499909204187064]
        @test forward.(compiled_model, [[5.0, 1.5], [4.0, 2.5]]) == [[7.4999983369439445, 13.49999833686844], [7.499909204262595, 12.499909204187064]]

        # check forward! function (single layer)
        params = [1.0, 2.0, -3.0]
        f = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        model = FlowModel( 2, (layer, ) )
        compiled_model = compile(model, params)
        output = zeros(2)
        forward!(output, compiled_model, [5.0, 1.5])
        @test output == [5.0, 7.4999983369439445]
        forward!(output, compiled_model, [4.0, 2.5])
        @test output == [4.0, 7.499909204262595]

        # check forward! function (multiple layers)
        params = [1.0, 2.0, -3.0, 1.0, 2.0, -3.0]
        f1 = PlanarFlow()
        layer1 = AdditiveCouplingLayer(f1)
        f2 = PlanarFlow()
        layer2 = AdditiveCouplingLayer(f2; permute=false)
        model = FlowModel( 2, (layer1, layer2) )
        compiled_model = compile(model, params)
        output = zeros(2)
        forward!(output, compiled_model, [5.0, 1.5])
        @test output == [7.4999983369439445, 13.49999833686844]
        forward!(output, compiled_model, [4.0, 2.5])
        @test output == [7.499909204262595, 12.499909204187064]

        # check backward function (single layer)
        params = [1.0, 2.0, -3.0]
        f = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        model = FlowModel( 2, (layer, ) )
        compiled_model = compile(model, params)
        @test backward(compiled_model, [5.0, 7.4999983369439445]) == [5.0, 1.5]
        @test backward(compiled_model, [4.0, 7.499909204262595]) == [4.0, 2.5]
        @test backward.(compiled_model, [[5.0, 7.4999983369439445], [4.0, 7.499909204262595]]) == [[5.0, 1.5], [4.0, 2.5]]

        # check backward function (multiple layers)
        params = [1.0, 2.0, -3.0, 1.0, 2.0, -3.0]
        f1 = PlanarFlow()
        layer1 = AdditiveCouplingLayer(f1)
        f2 = PlanarFlow()
        layer2 = AdditiveCouplingLayer(f2; permute=false)
        model = FlowModel( 2, (layer1, layer2) )
        compiled_model = compile(model, params)
        @test backward(compiled_model, [7.4999983369439445, 13.49999833686844]) == [5.0, 1.5]
        @test backward(compiled_model, [7.499909204262595, 12.499909204187064]) == [4.0, 2.5]
        @test backward.(compiled_model, [[7.4999983369439445, 13.49999833686844], [7.499909204262595, 12.499909204187064]]) == [[5.0, 1.5], [4.0, 2.5]]

        # check backward! function (single layer)
        params = [1.0, 2.0, -3.0]
        f = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        model = FlowModel( 2, (layer, ) )
        compiled_model = compile(model, params)
        output = zeros(2)
        backward!(output, compiled_model, [5.0, 7.4999983369439445])
        @test output == [5.0, 1.5]
        backward!(output, compiled_model, [4.0, 7.499909204262595])
        @test output == [4.0, 2.5]

        # check backward! function (multiple layers)
        params = [1.0, 2.0, -3.0, 1.0, 2.0, -3.0]
        f1 = PlanarFlow()
        layer1 = AdditiveCouplingLayer(f1)
        f2 = PlanarFlow()
        layer2 = AdditiveCouplingLayer(f2; permute=false)
        model = FlowModel( 2, (layer1, layer2) )
        compiled_model = compile(model, params)
        output = zeros(2)
        backward!(output, compiled_model, [7.4999983369439445, 13.49999833686844])
        @test output == [5.0, 1.5]
        backward!(output, compiled_model, [7.499909204262595, 12.499909204187064])
        @test output == [4.0, 2.5]

    end

    @testset "Jacobian" begin

        # check jacobian function (single layer)
        params = [1.0, 2.0, -3.0]
        f = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        model = FlowModel( 2, (layer, ) )
        compiled_model = compile(model, params)
        @test jacobian(compiled_model, [3.0, 1.5]) == [1.0 0.0; 1.0197320743308804 1.0]
        @test jacobian(compiled_model, [2.5, 5.0]) == [1.0 0.0; 1.1413016497063289 1.0]
        @test jacobian.(compiled_model, [[3.0, 1.5], [2.5, 5.0]]) == [[1.0 0.0; 1.0197320743308804 1.0], [1.0 0.0; 1.1413016497063289 1.0]]

        # check jacobian! function (single layer)
        params = [1.0, 2.0, -3.0]
        f = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        model = FlowModel( 2, (layer, ) )
        compiled_model = compile(model, params)
        output = zeros(2,2)
        jacobian!(output, compiled_model, [3.0, 1.5])
        @test output == [1.0 0.0; 1.0197320743308804 1.0]
        jacobian!(output, compiled_model, [2.5, 5.0])
        @test output == [1.0 0.0; 1.1413016497063289 1.0]

        # check jacobian function (multiple layers)
        params = [1.0, 2.0, -3.0, 1.0, 2.0, -3.0]
        f1 = PlanarFlow()
        layer1 = AdditiveCouplingLayer(f1)
        f2 = PlanarFlow()
        layer2 = AdditiveCouplingLayer(f2; permute=false)
        model = FlowModel( 2, (layer1, layer2) )
        compiled_model = compile(model, params)
        @test jacobian(compiled_model, [3.0, 1.5]) == [1.0197320743308804 1.0; 2.0197330107171334 1.0000009182669414]
        @test jacobian(compiled_model, [2.5, 5.0]) == [1.1413016497063289 1.0; 2.1413016497136192 1.0000000000063878]
        @test jacobian.(compiled_model, [[3.0, 1.5], [2.5, 5.0]]) == [[1.0197320743308804 1.0; 2.0197330107171334 1.0000009182669414], [1.1413016497063289 1.0; 2.1413016497136192 1.0000000000063878]]

        # check jacobian! function (multiple layers)
        params = [1.0, 2.0, -3.0, 1.0, 2.0, -3.0]
        f1 = PlanarFlow()
        layer1 = AdditiveCouplingLayer(f1)
        f2 = PlanarFlow()
        layer2 = AdditiveCouplingLayer(f2; permute=false)
        model = FlowModel( 2, (layer1, layer2) )
        compiled_model = compile(model, params)
        output = zeros(2,2)
        jacobian!(output, compiled_model, [3.0, 1.5])
        @test output == [1.0197320743308804 1.0; 2.0197330107171334 1.0000009182669414]
        jacobian!(output, compiled_model, [2.5, 5.0])
        @test output == [1.1413016497063289 1.0; 2.1413016497136192 1.0000000000063878]

        # check inv_jacobian function (single layer)
        params = [1.0, 2.0, -3.0]
        f = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        model = FlowModel( 2, (layer, ) )
        compiled_model = compile(model, params)
        @test inv_jacobian(compiled_model, [3.0, 1.5]) == [1.0 0.0; -1.0197320743308804 1.0]
        @test inv_jacobian(compiled_model, [2.5, 5.0]) == [1.0 0.0; -1.1413016497063289 1.0]
        @test inv_jacobian.(compiled_model, [[3.0, 1.5], [2.5, 5.0]]) == [[1.0 0.0; -1.0197320743308804 1.0], [1.0 0.0; -1.1413016497063289 1.0]]

        # check inv_jacobian! function (single layer)
        params = [1.0, 2.0, -3.0]
        f = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        model = FlowModel( 2, (layer, ) )
        compiled_model = compile(model, params)
        output = zeros(2,2)
        inv_jacobian!(output, compiled_model, [3.0, 1.5])
        @test output == [1.0 0.0; -1.0197320743308804 1.0]
        inv_jacobian!(output, compiled_model, [2.5, 5.0])
        @test output == [1.0 0.0; -1.1413016497063289 1.0]

        # check inv_jacobian! function (multiple layers)
        params = [1.0, 2.0, -3.0, 1.0, 2.0, -3.0]
        f1 = PlanarFlow()
        layer1 = AdditiveCouplingLayer(f1)
        f2 = PlanarFlow()
        layer2 = AdditiveCouplingLayer(f2; permute=false)
        model = FlowModel( 2, (layer1, layer2) )
        compiled_model = compile(model, params)
        output = zeros(2,2)
        inv_jacobian!(output, compiled_model, [3.0, 1.5])
        @test output == [-1.0197320743308804 1.0; 2.0197330107171334 -1.0000009182669414]
        inv_jacobian!(output, compiled_model, [2.5, 5.0])
        @test output == [-1.1413016497063289 1.0; 4.412130707993816 -2.9896834976728544]

        # check inv_jacobian function (multiple layers)
        params = [1.0, 2.0, -3.0, 1.0, 2.0, -3.0]
        f1 = PlanarFlow()
        layer1 = AdditiveCouplingLayer(f1)
        f2 = PlanarFlow()
        layer2 = AdditiveCouplingLayer(f2; permute=false)
        model = FlowModel( 2, (layer1, layer2) )
        compiled_model = compile(model, params)
        @test inv_jacobian(compiled_model, [3.0, 1.5]) == [-1.0197320743308804 1.0; 2.0197330107171334 -1.0000009182669414]
        @test inv_jacobian(compiled_model, [2.5, 5.0]) == [-1.1413016497063289 1.0; 4.412130707993816 -2.9896834976728544]
        @test inv_jacobian.(compiled_model, [[3.0, 1.5], [2.5, 5.0]]) == [[-1.0197320743308804 1.0; 2.0197330107171334 -1.0000009182669414], [-1.1413016497063289 1.0; 4.412130707993816 -2.9896834976728544]]

    end

    @testset "Joint processing functions" begin
        model = FlowModel(
            (
                InputLayer(8),
                AdditiveCouplingLayer(PlanarFlow()),
                AdditiveCouplingLayer(PlanarFlow(); permute=false)
            )
        )
        compiled_model = compile(model)
        x = randn(8)
        @test forward_jacobian(compiled_model, x) == (forward(compiled_model, x), jacobian(compiled_model, x))
        @test backward_inv_jacobian(compiled_model, x) == (backward(compiled_model, x), inv_jacobian(compiled_model, x))
        @test inv(jacobian(compiled_model, x)) ≈ inv_jacobian(compiled_model, forward(compiled_model, x))

    end

    @testset "Utility Jacobian" begin

        # check utility functions jacobian (single layer)
        params = [1.0, 2.0, -3.0]
        f = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        model = FlowModel( 2, (layer, ) )
        compiled_model = compile(model, params)
        @test det_jacobian(compiled_model, [1.5, 6.9]) == 1.0
        @test det_jacobian(compiled_model, [2.5, 6.4]) == 1.0
        @test absdet_jacobian(compiled_model, [1.5, 6.9]) == 1.0
        @test absdet_jacobian(compiled_model, [2.5, 6.4]) == 1.0
        @test isapprox(logdet_jacobian(compiled_model, [1.5, 6.9]), 0.0; atol=1e-10)
        @test isapprox(logdet_jacobian(compiled_model, [2.5, 6.4]), 0.0; atol=1e-10)
        @test sum(isapprox.(logabsdet_jacobian(compiled_model, [1.5, 6.9]), (0.0, 1.0); atol=1e-10)) == 2
        @test sum(isapprox.(logabsdet_jacobian(compiled_model, [2.5, 6.4]), (0.0, 1.0); atol=1e-10)) == 2

        # check utility functions jacobian (multiple layers)
        params = [1.0, 2.0, -3.0, 1.0, 2.0, -3.0]
        f1 = PlanarFlow()
        layer1 = AdditiveCouplingLayer(f1)
        f2 = PlanarFlow()
        layer2 = AdditiveCouplingLayer(f2; permute=false)
        model = FlowModel( 2, (layer1, layer2) )
        compiled_model = compile(model, params)
        @test detinv_jacobian(compiled_model, [1.5, 6.9]) ≈ -1.0
        @test detinv_jacobian(compiled_model, [2.5, 6.4]) ≈ -1.0
        @test absdetinv_jacobian(compiled_model, [1.5, 6.9]) ≈ 1.0
        @test absdetinv_jacobian(compiled_model, [2.5, 6.4]) ≈ 1.0
        @test sum(isapprox.(logabsdetinv_jacobian(compiled_model, [1.5, 6.9]), (0.0, -1.0); atol=1e-10)) == 2
        @test sum(isapprox.(logabsdetinv_jacobian(compiled_model, [2.5, 6.4]), (0.0, -1.0); atol=1e-10)) == 2

    end

end

end


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#goldschmidt-log.jl

include("../unum.jl")
using Unums

#implements an goldschmidt-like logarithm algorithm for Float64 as a testground for
#doing it in unums.

#returns a Int16 "1" if it's negative, "0" if it's positive.
function signof(x::Float64)
  ((reinterpret(UInt64, x) & 0x8000_0000_0000_0000) != 0) ? 1 : 0
end

function exponentof(x::Float64)
  Int64((reinterpret(Int64, x) & 0x7FF0_0000_0000_0000) >> 52 - 1023)
end

function castfrac(x::Float64)
    (reinterpret(UInt64, x) & 0x000F_FFFF_FFFF_FFFF) << 12
end

function maskfor(T::Type)
  if (T == Float64)
    return 0xFFFF_FFFF_FFFF_FF00  #56 digits suffices for Float64 (52 digits)
  elseif (T == Float32)
    return 0xFFFF_FFE0_0000_0000  #27 digits of precision seems to suffice for a Float32 (23 digits)
  end
end

o16 = one(UInt16)
z16 = zero(UInt16)
__clz_array=[0x0004,0x0003,0x0002,0x0002, o16, o16, o16, o16, z16,z16,z16,z16,z16,z16,z16,z16]
function leading_zeros(n)
  (n == 0) && return 64
  res::UInt16 = 0
  #use the binary search method
  (n & 0xFFFF_FFFF_0000_0000 == 0) && (n <<= 32; res += 0x0020)
  (n & 0xFFFF_0000_0000_0000 == 0) && (n <<= 16; res += 0x0010)
  (n & 0xFF00_0000_0000_0000 == 0) && (n <<= 8;  res += 0x0008)
  (n & 0xF000_0000_0000_0000 == 0) && (n <<= 4;  res += 0x0004)
  res + __clz_array[(n >> 60) + 1]
end

#simple fused - multiply - add.  Assumes num1 has hidden bits "carry" and num2
#has no hidden bits.
rm = 0x0000_0000_FFFF_FFFF
function sfma(carry, num1, num2)
  (fracprod, _) = Unums.__chunk_mult(num1, num2)
  (_carry, fracprod) = Unums.__carried_add(carry, num1, fracprod)
  ((carry & 0x1) != 0) && ((_carry, fracprod) = Unums.__carried_add(_carry, num2, fracprod))
  ((carry & 0x2) != 0) && ((_carry, fracprod) = Unums.__carried_add(_carry, lsh(num2, 1), fracprod))
  (_carry, fracprod)
end

#performs a simple multiply, Assumes that number 1 has a hidden bit of exactly one
#and number 2 has a hidden bit of exactly zero
#(1 + a)(0 + b) = b + ab
function smult(a::UInt64, b::UInt64)

  (fraction, _) = Unums.__chunk_mult(a, b)
  carry = one(UInt64)

  #only perform the respective adds if the *opposing* thing is not subnormal.
  ((carry, fraction) = Unums.__carried_add(carry, fraction, b))

  #carry may be as high as three!  So we must shift as necessary.
  (fraction, shift, is_ubit) = Unums.__shift_after_add(carry, fraction, _)
  fraction << 1
end

function reassemble(sign::UInt64, ev::UInt64, fv::UInt64)
  number = (fv >> 12) | ((ev + 1023) << 52) | (sign << 63)

  println(bits(number))

  reinterpret(Float64, number)
end

include("logtable.jl")
#ultimately, we may need to have more digits on the end of this value for logarithm.
const log2e = 0x71547652b82fe_000
const __logarithm_magicnumber = 0xb8aa3b27e0000000

function exlg(x::FloatingPoint)
  #exact floating point with the goldschmidt algorithm.
  T = typeof(x)

  isnan(x) && return (nan(T), false)

  x <= 0 && return nan(T)

  #calculate the exponent.
  exp_f::Int64 = exponentof(x) + (issubnormal(x) ? 1 : 0)

  #figure the decimals.
  fraction::UInt64 = castfrac(x)

  if (issubnormal(x))
<<<<<<< HEAD:oddsandends/goldschmidt-log.jl
    shift::UInt64 = clz(fraction) + 1
=======
    shift::UInt64 = leading_zeros(fraction) + 1
>>>>>>> directory-reorg:src/oddsandends/goldschmidt-log.jl
    fraction = fraction << shift
    exp_f -= shift
  end

  sign::UInt64 = 0
  if (exp_f < 0)
    exp_f = -exp_f - 1
    sign = 1
  end
<<<<<<< HEAD:oddsandends/goldschmidt-log.jl
  lz = clz(UInt64(exp_f))
=======
  lz = leading_zeros(UInt64(exp_f))
>>>>>>> directory-reorg:src/oddsandends/goldschmidt-log.jl
  resexp = 63 - lz
  #add the exponent part onto the result fraction.
  resfrac = UInt64(exp_f << (lz + 1))

  #do the goldschmidt-type algorithm.
  #first, "divide by two" by doing a virtual shift left, appending the implied
  #one onto the most significant end.
  fraction = (fraction >> 1) | 0x8000_0000_0000_0000
  diff::UInt64 = 0
  d::Int64 = 0
  m::UInt64 = 0
  intsumdelta::UInt64 = 0
  intsumsofar::UInt64 = 0

  #goldschmidt-like algorithm.
  for idx = 1:32
    #figure the difference between the fraction we have and what's left.  This
    #is equivalent to the operation "2 - x".  For 0.5 < x < 1.0 x(2 - x) is
    #bounded between 0.5 and 1, so we get closer by iterative multiplication.
    #However, because calculating the log of this value is nontrivial, we must
    #use a lookup table which has precalculated logs for fractions f of the form
    # f = 1.0...010....0, indexed by the place of the one.
    diff = 0 - fraction
    #find the place of the top bit
    d = leading_zeros(diff)
    #generate the value (we will do a multiply with this to keep track of our progress.)
    m = 0x8000_0000_0000_0000 >> d

    #look up the value to be added, or, if it's far enough along, just do a bitshift.
    intsumdelta = d > 29 ? (__logarithm_magicnumber >> d) : lt64[d]
    intsumsofar += intsumdelta

    #update frac to be the product of our cumulative fraction and the 1.0...01
    (_, fraction) = Unums.__sfma(zero(UInt64), fraction, m)

    #terminating condition.
    if m < 0x0000_0000_0000_0400
      break
    end
  end

  #now do a flip (if the log is positive)
  iint::UInt64 = (sign == 0) ? -intsumsofar : intsumsofar
  resfrac |= iint >> resexp
  #reassemble the value into the requisite floating point
  reassemble(sign, resexp, resfrac)
end

v = rand(Int64)
x = abs(reinterpret(Float64, v))
z = exlg(x)
a = log2(x)

println("input :    $x")
println("answer:    $(log2(x))")
println("calculate: $z")

println("xbits:     ", bits(x))
println("abits:     ", bits(a))
println("zbits:     ", bits(z))

println("diff:      ", reinterpret(UInt64, a) - reinterpret(UInt64, z))

#######################################################
## testing fun

i2f(i) = reinterpret(Float64, (i >> 12) | 0x3FF0_0000_0000_0000)

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using wgregseq, DataFrames, CSV, BioSequences, CairoMakie, Statistics, StatsBase

# Set path
dir = @__DIR__
home_dir = joinpath(split(dir, "/")[1:end-2])

##
filename = "2022-02-15_twist_order.csv"

df = CSV.read(
    "/$home_dir/data/twist_orders/$filename",
    types=Dict(
        "fwd_primer" => String,
        "rev_primer1" => String,
        "rev_primer2" => String,
        "rev_primer3" => String,
    ),
    DataFrame
)
df.fwd_primer = parse.(Tuple{Int64, Tuple{Int64, Int64}}, df.fwd_primer)
df.rev_primer1 = parse.(Tuple{Int64, Tuple{Int64, Int64}}, df.rev_primer1)
df.rev_primer2 = parse.(Tuple{Int64, Tuple{Int64, Int64}}, df.rev_primer2)
df.rev_primer3 = parse.(Tuple{Int64, Tuple{Int64, Int64}}, df.rev_primer3)

##
primer_df_IDT = DataFrame()
primer_df = DataFrame()

overhang = lowercase(wgregseq.design.import_primer(300, "rev"))

ind0 = 142
i = 0

enzyme_df = wgregseq.enzyme_list

SbfI_site = enzyme_df[enzyme_df.enzyme .== "SbfI", "site"][1]
SalI_site = enzyme_df[enzyme_df.enzyme .== "SalI", "site"][1]

barcode_IDT =  "(N:25252525)" * join((fill("(N)", 19)))
barcode =  join((fill("N", 20)))

for primer in unique(df.fwd_primer)
    primer_ind = primer[1]
    primer_seq = wgregseq.design.import_primer(primer_ind, "fwd")
    name = "SC$(ind0+i)_$(primer_ind)_fwd"
    append!(primer_df, DataFrame(name=name, sequence=primer_seq))
    i += 1
end

for primer in unique(df.rev_primer1)
    primer_ind = primer[1]
    primer_seq = overhang * reverse_complement(wgregseq.design.import_primer(primer_ind, "rev"))
    name = "SC$(ind0+i)_$(primer_ind)_rev"
    append!(primer_df, DataFrame(name=name, sequence=primer_seq))
    i += 1
    #=
    bar_primer_seq = string(reverse_complement(LongDNA{4}(SalI_site))) *
                        barcode_ *
                        string(reverse_complement(LongDNA{4}(SbfI_site))) *
                        string(primer_seq)
    bar_primer_name = "SC$(ind0+i)_$(primer_ind)_SbfI_BC_SalI_rev"
    append!(primer_df, DataFrame(name=bar_primer_name, sequence=bar_primer_seq))
    i += 1
    println(bar_primer_seq)
    =#
end

for primer in unique(df.rev_primer2)
    primer_ind = primer[1]
    primer_seq = overhang * reverse_complement(wgregseq.design.import_primer(primer_ind, "rev")[1:19])
    name = "SC$(ind0+i)_$(primer_ind)_rev"
    append!(primer_df, DataFrame(name=name, sequence=primer_seq))
    i += 1
    #=
    bar_primer_seq = string(reverse_complement(wgregseq.design.import_primer(primer_ind, "rev")))
    bar_primer_seq = reverse_complement(LongDNA{4}(bar_primer_seq * SbfI_site * barcode * SalI_site))
    bar_primer_name = "SC$(ind0+i)_$(primer_ind)_SbfI_BC_SalI_rev"
    append!(primer_df, DataFrame(name=bar_primer_name, sequence=string(primer_seq)))
    i += 1
    println(bar_primer_seq)
    =#
end

for primer in unique(df.rev_primer3)
    primer_ind = primer[1]
    primer_seq = overhang * reverse_complement(wgregseq.design.import_primer(primer_ind, "rev")[1:19])
    name = "SC$(ind0+i)_$(primer_ind)_rev"
    append!(primer_df, DataFrame(name=name, sequence=primer_seq))
    i += 1
    #=
    bar_primer_seq = string(reverse_complement(wgregseq.design.import_primer(primer_ind, "rev")))
    bar_primer_seq = reverse_complement(LongDNA{4}(bar_primer_seq * SbfI_site * barcode * SalI_site))
    bar_primer_name = "SC$(ind0+i)_$(primer_ind)_SbfI_BC_SalI_rev"
    append!(primer_df, DataFrame(name=bar_primer_name, sequence=string(primer_seq)))
    i += 1
    println(bar_primer_seq)
    =#
end
spacer = "CGATAAAC"
barcoding_primer_IDT = spacer * string(reverse_complement(LongDNA{4}(SalI_site))) *
                        barcode_IDT *
                        string(reverse_complement(LongDNA{4}(SbfI_site))) *
                        string(overhang)

barcoding_primer = LongDNA{4}(spacer) * reverse_complement(LongDNA{4}(SalI_site)) *
                        LongDNA{4}(barcode) *
                        reverse_complement(LongDNA{4}(SbfI_site)) *
                        overhang
append!(primer_df, DataFrame(name="SC$(ind0+i)_300_SbfI_BC_SalI_rev", sequence=barcoding_primer))
println("Barcoding Primer for IDT: ", barcoding_primer_IDT)
primer_df
CSV.write("/$home_dir/data/twist_orders/primer_table.csv", primer_df)


text:
using Dates

const pagesizes = Dict( "lpr" => [132, 66], "tn3270" => [80, 43])

pagefit(prn) = haskey(pagesizes, prn) ?
    [div(pagesizes[prn][1], 22), div(pagesizes[prn][2], 12)] : [1, 1]
pagecols(prn) = haskey(pagesizes, prn) ? pagesizes[prn][1] : 20

function centerobject(x, cols)
    content = string(x)
    rpad(lpad(content, div(cols + length(content), 2)), cols)
end

function ljustlines(x, cols)
    arr = Vector{String}()
    for s in split(x, "\n")
        push!(arr, rpad(s, cols)[1:cols])
    end
    join(arr, "\n")
end

function formatmonth(yr, mo)
    dt = Date("$yr-$mo-01")
    dayofweekfirst = dayofweek(dt)
    numweeklines = 1
    str = centerobject(monthname(dt), 20) * "\nMo Tu We Th Fr Sa Su\n"
    str *= " " ^ (3 * (dayofweekfirst - 1)) * lpad(string(1), 2)
    for i = 2:daysinmonth(dt)
        if (i + dayofweekfirst + 5) % 7 == 0
            str *= "\n" * lpad(i, 2)
            numweeklines += 1
        else
            str *= lpad(string(i), 3)
        end
    end
    str *= numweeklines < 6 ? "\n\n\n" : "\n\n"
    ljustlines(str, 20)
end

function formatyear(displayyear, printertype)
    calmonths = [formatmonth(displayyear, mo) for mo in 1:12]
    columns = pagecols(printertype)
    monthsperline = pagefit(printertype)[1]
    joinspaces = max( (monthsperline > 1) ?
        div(columns - monthsperline * 20, monthsperline - 1) : 1, 1)
    str = "\n" * centerobject(displayyear, columns) * "\n"
    monthcal = [split(formatmonth(displayyear, i), "\n") for i in 1:12]
    for i in 1:monthsperline:length(calmonths) - 1
        for j in 1:length(monthcal[1])
            monthlines = map(x->monthcal[x][j], i:i + monthsperline - 1)
            str *= rpad(join(monthlines, " " ^ joinspaces), columns) * "\n"
        end
        str *= "\n"
    end
    str
end

function lineprintcalendar(years)
    for year in years, printer in keys(pagesizes)
        println(formatyear(year, printer))
    end
end

lineprintcalendar(1969)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    BHMM([a::AbstractVector{T}, ]A::AbstractMatrix{T}, B::AbstractVector{<:Categorical}) where F where T
Build an BHMM with transition matrix `A` and observations distributions `B`.
If the initial state distribution `a` is not specified, a uniform distribution is assumed.
Observations distributions can be of different types (for example `Normal` and `Exponential`).
However they must be of the same dimension (all scalars or all multivariates).
# Example
```julia
hmm = BHMM([0.9 0.1; 0.1 0.9], [Normal(0,1), Normal(10,1)])
```
"""
struct BHMM{T} <: AbstractHMM{Univariate}
    a::AbstractVector{T}
    A::AbstractMatrix{T}
    B::AbstractVector{<:Categorical}
    partition::String
    BHMM{T}(a, A, B, partition="") where {T} = assert_BHMM(a, A, B) && new(a, A, B, partition)
end

BHMM(a::AbstractVector{T}, A::AbstractMatrix{T}, B::AbstractVector{<:Categorical}, partition ...) where {T} = BHMM{T}(a, A, B, partition ...)
BHMM(A::AbstractMatrix{T}, B::AbstractVector{<:Categorical}, partition...) where {T} = BHMM{T}(ones(size(A)[1])/size(A)[1], A, B, partition...)

copy(hmm::BHMM) = BHMM(copy(hmm.a), copy(hmm.A), copy(hmm.B), hmm.partition)
size(hmm::BHMM) = (length(hmm.B), length(hmm.B[1].p))

"""
    assert_BHMM(a, A, B)
Throw an `ArgumentError` if the initial state distribution and the transition matrix rows does not sum to 1,
and if the observations distributions does not have the same dimensions.
"""
function assert_BHMM(a::AbstractVector,
                    A::AbstractMatrix,
                    B::AbstractVector{<:Categorical})
    !isprobvec(a) && throw(ArgumentError("Initial state vector a is not a valid probability vector!"))
    !istransmat(A) && throw(ArgumentError("Transition matrix A is not valid!"))

    !all([length(d.p) for d in B] .== length(B[1].p)) && throw(ArgumentError("All distributions must have the same dimensions"))
    !(length(a) == size(A,1) == length(B)) && throw(ArgumentError("Length of initial state vector a, dimension of transition matrix A, and number of distributions B are not the same!"))
    return true
end

function Base.show(io::IO, hmm::BHMM)
    println(io, "Background BHMM")
    println(io, "State Initial and Transition Probabilities")
    print(io, "a: ")
    show(io, hmm.a)
    println(io)
    print(io, "A: ")
    display(hmm.A)
    println(io)
    println(io, "INFORMATIVE SYMBOLS BY STATE")
    for (n,d) in enumerate(hmm.B)
        print(io, "K$n ")
        print_emitters(d)
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

if !isdefined(@__MODULE__, :DefsChannelId)
    const DefsChannelId = String
else
    @warn("Skipping redefinition of DefsChannelId to String")
end


text:
# dynamic_constraints.jl
# Dynamic constraint generation helper functions
# Authored by Arthur J Delarue on 9/27/15

MANHATTAN_NODES = 6134

function chooseStartingConstraints(
	travelTimes::Array{Float64, 2},
	numRides::Array{Int, 2},
	pairs::Array{Array{Int,1},1};
	sample_size::Int = 500,
	num_nodes::Int = MANHATTAN_NODES)
	"""
	Chooses constraints that algorithm will start with. Takes in matrix array of travel times and number of rides (and starting sample size) and returns arrays of travel times/number of rides of same decisions but only relevant info, plus 1D arrays with sources and destinations.
	"""
	# Randomly pick starting constraints
	newTravelTimes, newNumRides = chooseConstraints(travelTimes, numRides, sample_size = sample_size, num_nodes=num_nodes)
	# Create vectors for sources and destinations
	srcs = Int[]
	dsts = Int[]
	sizehint!(srcs, sample_size)
	sizehint!(dsts, sample_size)
	for i in 1:size(travelTimes)[1], j in 1:size(travelTimes)[2]
		if newTravelTimes[i,j] > 0
			# Load sources and destinations: will be useful when adding constraints
			push!(srcs, i)
			push!(dsts, j)
			push!(pairs[i], j)
		end
	end
	return srcs, dsts, newTravelTimes, newNumRides, pairs
end

function updateConstraints(
	travelTimes::Array{Float64,2},
	numRides::Array{Int,2},
	calculatedTravelTimes::Array{Float64,2},
	totalPaths::Array{Any, 1},
	totalNumExpensiveTurns::Array{Any, 1},
	numPaths::Array{Int, 1},
	srcs::Array{Int, 1},
	dsts::Array{Int, 1},
	pairs::Array{Array{Int, 1}, 1};
	numNodePairsToAdd::Int = 0,
	numNodePairsToRemove::Int = 0,
	addOnlyIfAbovePercentile::Float64 = 0.5,
	removeOnlyIfBelowPercentile::Float64 = 0.5)

	"""
	Add and remove node pairs based on selection criteria.
	"""

	# Create new return structures
	newSrcs = deepcopy(srcs)
	newDsts = deepcopy(dsts)
	newPairs = deepcopy(pairs)
	newTotalPaths = deepcopy(totalPaths)
	newTotalNumExpensiveTurns = deepcopy(totalNumExpensiveTurns)
	newNumPaths = deepcopy(numPaths)

	# Create lookup set for constraints already under consideration
	nodePairsAlreadyIn = Set([(srcs[i], dsts[i]) for i = 1:length(srcs)])
	# Select new constraints - first, put everything in vector form to prepare for sorting
	indicesVector = Tuple{Int,Int}[]
	sizehint!(indicesVector, 12000000)
	for i = 1:size(travelTimes)[1], j = 1:size(travelTimes)[1]
		if travelTimes[i,j] > 0
			push!(indicesVector, (i,j))
		end
	end
	# then, compute distance between calculated times and data for all pairs
	newErrorVector = zeros(length(indicesVector))
	for i = 1:length(indicesVector)
		newErrorVector[i] = sqrt(numRides[indicesVector[i][1], indicesVector[i][2]]/travelTimes[indicesVector[i][1], indicesVector[i][2]]) * abs(travelTimes[indicesVector[i][1], indicesVector[i][2]] - calculatedTravelTimes[indicesVector[i][1], indicesVector[i][2]])
	end
	# compute error for nodepairs already in LP
	errorVector = zeros(length(srcs))
	for i = 1:length(srcs)
		errorVector[i] = sqrt(numRides[srcs[i], dsts[i]]/travelTimes[srcs[i], dsts[i]]) * abs(travelTimes[srcs[i], dsts[i]] - calculatedTravelTimes[srcs[i], dsts[i]])
	end
	# sort the ones already in the data by error
	p2 = sortperm(errorVector)
	# sort everything by error
	p = sortperm(newErrorVector)
	# Find ranges where it is acceptable to add/remove constraints
	minIndexToAdd = round(Int, addOnlyIfAbovePercentile * length(indicesVector))
	maxIndexToRemove = round(Int, removeOnlyIfBelowPercentile * length(newSrcs))
	# Now, remove best node pairs - create vector of indices to keep, then splice out the indices we don't want to keep
	indicesToKeep = collect(1:length(newSrcs))
	numNodePairsRemoved = 0
	while numNodePairsRemoved < min(numNodePairsToRemove, maxIndexToRemove)
		i = rand(1:maxIndexToRemove)
		if (newSrcs[p2[i]], newDsts[p2[i]]) in nodePairsAlreadyIn
			idx = findfirst(indicesToKeep, p2[i])
			splice!(indicesToKeep, idx)
			idx = findfirst(newPairs[newSrcs[p2[i]]], newDsts[p2[i]])
			splice!(newPairs[newSrcs[p2[i]]], idx)
			delete!(nodePairsAlreadyIn, (newSrcs[p2[i]], newDsts[p2[i]]))
			numNodePairsRemoved += 1
		end
	end
	# Since we spliced out the unwanted indices, this performs the necessary reduction
	newSrcs = newSrcs[indicesToKeep]
	newDsts = newDsts[indicesToKeep]
	newTotalNumExpensiveTurns = newTotalNumExpensiveTurns[indicesToKeep]
	newNumPaths = newNumPaths[indicesToKeep]
	newTotalPaths = newTotalPaths[indicesToKeep]
	# Add poor node pairs - this is easier, just push the new nodepairs to the end :)
	numNodePairsAdded = 0
	while numNodePairsAdded < numNodePairsToAdd && length(nodePairsAlreadyIn) < (length(indicesVector) - minIndexToAdd)
		i = rand(minIndexToAdd:length(indicesVector))
		if !(indicesVector[p[i]] in nodePairsAlreadyIn)
			push!(newSrcs, indicesVector[p[i]][1])
			push!(newDsts, indicesVector[p[i]][2])
			push!(newTotalPaths, Array{Int}[])
			push!(newTotalNumExpensiveTurns, Int[])
			push!(newNumPaths, 0)
			push!(newPairs[indicesVector[p[i]][1]], indicesVector[p[i]][2])
			sort!(newPairs[indicesVector[p[i]][1]])
			numNodePairsAdded += 1
			push!(nodePairsAlreadyIn, indicesVector[p[i]])
		end
	end
	return newSrcs, newDsts, newTotalPaths, newTotalNumExpensiveTurns, newNumPaths, newPairs
end

function updatePaths(
	paths::Array{Array{Int64,1}, 1},
	numExpensiveTurns::Array{Int,1},
	totalPaths::Array{Any, 1},
	totalNumExpensiveTurns::Array{Any,1},
	numPaths::Array{Int, 1};
	maxNumPathsPerOD::Int=3,
	times::AbstractArray{Float64,2}=zeros(1,1),
	turnCost::Float64=0.0,
	travelTimes::AbstractArray{Float64,2}=zeros(1,1),
	numRides::AbstractArray{Int,2}=zeros(Int,(1,1)),
	srcs::Array{Int,1}=zeros(1),
	dsts::Array{Int,1}=zeros(1),
	dynamicConstraints::Bool=true,
	globalConstraintUpdate::Bool=true
	)
	"""
	Update path array so as not to exceed max number of paths.
	"""
	newTotalPaths = [totalPaths[i] for i = 1:length(totalPaths)]
	newTotalNumExpensiveTurns = [totalNumExpensiveTurns[i] for i = 1:length(totalNumExpensiveTurns)]
	totalPaths = newTotalPaths
	totalNumExpensiveTurns = newTotalNumExpensiveTurns
	if dynamicConstraints && globalConstraintUpdate
		# Deal with single path per O,D case separately because of the bugs it's been causing
		if maxNumPathsPerOD > 1
			for i=1:length(paths)
				index = findfirst(totalPaths[i], paths[i])
				# If path already in paths
				if index != 0
					totalPaths[i][1], totalPaths[i][index] = totalPaths[i][index], totalPaths[i][1]
					totalNumExpensiveTurns[i][1], totalNumExpensiveTurns[i][index] = totalNumExpensiveTurns[i][index], totalNumExpensiveTurns[i][1]			# New path is equality constraint
				# If we still have space to add the path
				else
					numPaths[i] += 1
					if numPaths[i] == 1
						push!(totalPaths[i], paths[i])
						push!(totalNumExpensiveTurns[i], numExpensiveTurns[i])
					else
						push!(totalPaths[i], paths[i])
						push!(totalNumExpensiveTurns[i], numExpensiveTurns[i])
						assert(numPaths[i] == length(totalPaths[i]))
						totalPaths[i][numPaths[i]], totalPaths[i][1] = totalPaths[i][1], totalPaths[i][numPaths[i]]
						totalNumExpensiveTurns[i][numPaths[i]], totalNumExpensiveTurns[i][1] = totalNumExpensiveTurns[i][1], totalNumExpensiveTurns[i][numPaths[i]]
					end
				end
			end
			# Check if paths need to be removed, and if so, remove them
			nPaths = sum(numPaths)
			# Check if we have too many paths
			if nPaths > maxNumPathsPerOD * length(totalPaths)
				# Put all the paths that could be removed (i.e. not equality paths)
				# into a big array (indicesVector)
				# so that we will be able to find them after we sort
				indicesVector = Tuple{Int,Int,Int,Int}[]
				sizehint!(indicesVector, nPaths)
				for i = 1:length(totalPaths), j=1:numPaths[i]
					if j > 1
						push!(indicesVector, (i,j,srcs[i],dsts[i]))
					end
				end
				errorVector = zeros(length(indicesVector))
				# Check if there are any non-equality paths
				if length(indicesVector) > 0
					# Compute error on each path (difference between it and the equality path,
					# 						weighted by path time)
					for i=1:length(indicesVector)
						errorVector[i] = abs(sum([times[totalPaths[indicesVector[i][1]][indicesVector[i][2]][a], totalPaths[indicesVector[i][1]][indicesVector[i][2]][a+1]] for a = 1:(length(totalPaths[indicesVector[i][1]][indicesVector[i][2]])-1)]) + totalNumExpensiveTurns[indicesVector[i][1]][indicesVector[i][2]] * turnCost - sum([times[totalPaths[indicesVector[i][1]][1][a], totalPaths[indicesVector[i][1]][1][a+1]] for a = 1:(length(totalPaths[indicesVector[i][1]][1])-1)]) - totalNumExpensiveTurns[indicesVector[i][1]][1] * turnCost)/(travelTimes[indicesVector[i][3], indicesVector[i][4]])
					end
					# Sort paths by "badness"
					p = sortperm(errorVector)
					# Figure out how many paths need to be removed
					numPathsToRemove = nPaths - maxNumPathsPerOD * length(totalPaths)
					# Initialize an array of indices of the paths that need to be kept (so far all paths are in there)
					pathsToKeep = [collect(1:numPaths[i]) for i=1:length(numPaths)]
					for i = 1:numPathsToRemove
						# Remove appropriate path from index list
						idx = findfirst(pathsToKeep[indicesVector[p[end+1-i]][1]], indicesVector[p[end+1-i]][2])
						splice!(pathsToKeep[indicesVector[p[end+1-i]][1]], idx)
					end
					# Rebuild path arrays using pathsToKeep to filter out removed paths
					totalPaths = [totalPaths[i][pathsToKeep[i]] for i = 1:length(totalPaths)]
					totalNumExpensiveTurns = [totalNumExpensiveTurns[i][pathsToKeep[i]] for i = 1:length(totalNumExpensiveTurns)]
					numPaths = [length(pathsToKeep[i]) for i = 1:length(numPaths)]
				end
			end
		else
			# Recreate path and turn arrays from scratch (since we will overwrite all their contents anyways)
			numDataPoints = length(totalPaths)
			assert(numDataPoints == length(paths))
			totalPaths = Array(Any, numDataPoints)
			totalNumExpensiveTurns = Array(Any, numDataPoints)
			for i = 1:numDataPoints
				totalPaths[i] = Array{Int}[]
				totalNumExpensiveTurns[i] = Int[]
			end
			# If one path per constraint, just update the path... easy-peasy
			for i=1:length(paths)
				push!(totalPaths[i], paths[i])
				push!(totalNumExpensiveTurns[i], numExpensiveTurns[i])
			end
		end
	# Non-dynamic constraint case
	else
		for i=1:length(paths)
			index = findfirst(totalPaths[i], paths[i])
			# If path already in paths
			if index != 0
				totalPaths[i][1], totalPaths[i][index] = totalPaths[i][index], totalPaths[i][1]
				totalNumExpensiveTurns[i][1], totalNumExpensiveTurns[i][index] = totalNumExpensiveTurns[i][index], totalNumExpensiveTurns[i][1]			# New path is equality constraint
			# If we still have space to add the path
			elseif numPaths[i] < maxNumPathsPerOD
				numPaths[i] += 1
				if numPaths[i] == 1
					push!(totalPaths[i], paths[i])
					push!(totalNumExpensiveTurns[i], numExpensiveTurns[i])
				else
					push!(totalPaths[i], paths[i])
					push!(totalNumExpensiveTurns[i], numExpensiveTurns[i])
					assert(numPaths[i] == length(totalPaths[i]))
					totalPaths[i][numPaths[i]], totalPaths[i][1] = totalPaths[i][1], totalPaths[i][numPaths[i]]
					totalNumExpensiveTurns[i][numPaths[i]], totalNumExpensiveTurns[i][1] = totalNumExpensiveTurns[i][1], totalNumExpensiveTurns[i][numPaths[i]]
				end
			# If we need to remove a path
			else
				worstIndex = findWorstPathIndex(totalPaths[i], totalNumExpensiveTurns[i], turnCost, times)
				if worstIndex == 1
					totalPaths[i][1] = paths[i]
					totalNumExpensiveTurns[i][1] = numExpensiveTurns[i]
				else
					totalPaths[i][1], totalPaths[i][worstIndex] = paths[i], totalPaths[i][1]
					totalNumExpensiveTurns[i][1], totalNumExpensiveTurns[i][worstIndex] = numExpensiveTurns[i], totalNumExpensiveTurns[i][1]
				end
			end
		end
	end
	return totalPaths, totalNumExpensiveTurns, numPaths
end

function loadTrainingResults(directoryName::AbstractString, max_iterations::Int, manhattan::Manhattan, turnCost::Float64)
	"""
	Given the name of the directory where the algorithm output was saved, run shortest paths to determine travel times between all pairs of nodes.
	"""
	println("-- Loading results from $(directoryName)...")
	link_times = load("Outputs/$(directoryName)/manhattan-times-$(max_iterations).jld", "times")
	println("-- Computing shortest paths...")
	new_graph, new_edge_dists, new_nodes = modifyGraphForDijkstra(manhattan.network, link_times, manhattan.positions, turn_cost=turnCost)
	old_nodes = getInverseMapping(new_nodes, nv(new_graph))
	new_sp = parallelShortestPathsWithTurnsAuto(manhattan.network, new_graph, new_edge_dists, new_nodes)
	return new_sp.traveltime
end

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing


text:
type Clock
	t::Int
	final_t::Int
end

function gettimestep(c::Clock)
	return c.t
end

function move_forward(c::Clock)
	c.t = c.t + 1
	nothing
end

function finished(c::Clock)
	if c.t>c.final_t
		return true
	else
		return false
	end
end


text:
function create_data_for_test()
    data = zeros(Float32, 2, 1000)
    data[:, 1:250] .= [-1, -1]
    data[:, 251:500] .= [-1, 1]
    data[:, 501:750] .= [1, -1]
    data[:, 751:1000] .= [1, 1]
    return data
end

@testset "Testing Module (Determinstic)" begin
    data = create_data_for_test()
    labels, clusters, weights = fit(data, 100.0, iters = 200, seed = 12345, burnout = 5)
    @test all(data[:, 1:250] .== data[:, 1])
    @test all(data[:, 251:500] .== data[:, 251])
    @test all(data[:, 501:750] .== data[:, 501])
    @test all(data[:, 751:1000] .== data[:, 751])
    @test length(clusters) == 4
    println(weights)
    @test all(weights .>= 0.15)
    labels_histogram = get_labels_histogram(labels)
    for (k, v) in labels_histogram
        @test v == 250
    end
end

@testset "Testing Module (Random mess)" begin
    @everywhere Random.seed!(12345)
    x, labels, clusters = generate_gaussian_data(10^5, 3, 10, 100.0)

    hyper_params = DPMMSubClusters.niw_hyperparams(
        Float32(1.0),
        zeros(Float32, 3),
        Float32(5),
        Matrix{Float32}(I, 3, 3) * 1,
    )

    dp = dp_parallel(
        x,
        hyper_params,
        Float32(1000000000000000000000.0),
        100,
        1,
        nothing,
        true,
        false,
        15,
        labels,
    )
    @test length(dp[1].group.local_clusters) > 1
end

@testset "Multinomial Module And save load" begin
    @everywhere Random.seed!(12345)
    x, labels, clusters = generate_mnmm_data(10^3, 100, 20, 50)
    @test size(x, 1) == 100
    @test size(x, 2) == 10^3
    npzwrite("save_load_test/mnm_data.npy", x')

    dp = dp_parallel("save_load_test/multinomial_params.jl")
    @test length(dp[1].group.local_clusters) > 1
    dp = run_model_from_checkpoint("save_load_test/checkpoint_20.jld2")
    @test length(dp[1].group.local_clusters) > 1
end


text:
PARAM[:capacity] = Dict(
  :futuresfile=>"BBFutures",
  :futuresfileext=>"csv",
  :futurescolumns => [:fullticker, :description, :exchange,
    :volume10d,	:openinterest10d,	:subcategory],
  :periodicity => ["DAILY", "WEEKLY", "MONTHLY", "QUARTERLY", "YEARLY"],
  :pullfields => ["FUT_AGGTE_VOL", "PX_VOLUME", "PX_LAST", "PX_SETTLE",
    "CONTRACT_VALUE", "OPEN_INT"])

function futuresdata(;
  inputpath::String = PARAM[:inputpath],
  futuresfile::String = PARAM[:capacity][:futuresfile],
  futuresfileext::String = PARAM[:capacity][:futuresfileext],
  )

  futures = CSV.read("$inputpath\\$futuresfile.$futuresfileext") |> DataFrame!
  select!(futures, PARAM[:capacity][:futurescolumns])

  parsemfloat(f)=something(tryparse(Float64, f), missing)
  (:test ∈ propertynames(futures)) && (futures.test = parsemfloat.(futures.test))

  #display(describe(futures))
  return futures
end

#a function for pulling lots of tickers and automating common parameter choices
function multipullold(S;
  livepull::Bool="live pull is required",
  tickers::Vector{String}=error("tickers kw are required"),
  fields::Vector{String}=error("fields kw are required"),
  startdate::Date=Date(1900,1,1),
  enddate::Date=today(),
  periodicity::Vector{String}=error("periodicity kw are required"),
  CshAdjAbnormal::Bool =	true,
  CshAdjNormal::Bool =	true,
  currency::String=	"USD",
  DateFormat::String = "D",
  pullname::String=error("pullname is required"),
  inbinstream=IN_BIN_STREAM,
  outbinstream=OUT_BIN_STREAM,
  binext = BIN_EXTENSION)

  local results
  pullpath::String = "$(PARAM[:inputpath])\\$pullname.$binext"

  #WARNING- what to do about dates?

  #set up a parameter dict

  #this will hold each of the result dictionaries
  #allunique(tickers) || error("Tickers are not unique!")

  #batch the tickers for performance
  Ntickers::Int = length(tickers)
  Nperiods::Int = length(periodicity)
  tickerparts=Vector.(Iterators.partition(tickers, Threads.nthreads()) |> collect)
  Nparts::Int = length(tickerparts)

  sessions = (_->startbbsession()).(1:Threads.nthreads())

  if livepull
    try
      resultparts::Vector{Vector{Dict}} = [Vector{Dict}(undef, Nparts) for _ ∈ 1:Nperiods]
      @sync for (i, p) ∈ enumerate(periodicity)
        options::Dict{String, Any} = Dict(
          "CshAdjAbnormal"=>CshAdjAbnormal,
          "CshAdjNormal"=>CshAdjNormal,
          "currency"=>currency,
          #"DateFormat"=>DateFormat,
          )
        if p=="DAILY"
          #options["CDR"] = "7D"
          options["overrideOption"] = "OVERRIDE_OPTION_CLOSE"
          options["nonTradingDayFillOption"] = "ACTIVE_DAYS_ONLY"
          options["nonTradingDayFillMethod"] = "NIL_VALUE"
        end

        #query BB for each partition
        Threads.@threads for j ∈ 1:Nparts
          @async resultparts[i][j] = BLPData.bdh(sessions[Threads.threadid()], tickerparts[j], fields, startdate, enddate,
            periodicity=p, options=options)
        end
      end

    finally
      endbbsession.(sessions)
    end
    #combine the results
    @assert sum([(resultperiod->length.(resultperiod)).(resultparts)...;]) == Ntickers*Nperiods
    results = (resultperiod->reduce(merge, resultperiod)).(resultparts)
    @assert length(results) == Nperiods

    #write out the results
    outbinstream(pullpath, results)
  else
    results = inbinstream(pullpath)
  end

  #println(results)
  #println(results)
  results
end

function multipull(S;
  livepull::Bool="live pull is required",
  tickers::Vector{String}=error("tickers kw are required"),
  fields::Vector{String}=error("fields kw are required"),
  startdate::Date=Date(1900,1,1),
  enddate::Date=today(),
  periodicity::Vector{String}=error("periodicity kw are required"),
  CshAdjAbnormal::Bool =	true,
  CshAdjNormal::Bool =	true,
  currency::String=	"USD",
  DateFormat::String = "D",
  pullname::String=error("pullname is required"),
  inbinstream=IN_BIN_STREAM,
  outbinstream=OUT_BIN_STREAM,
  binext = BIN_EXTENSION)

  local results
  pullpath::String = "$(PARAM[:inputpath])\\$pullname.$binext"

  #WARNING- what to do about dates?

  #set up a parameter dict

  #this will hold each of the result dictionaries
  #allunique(tickers) || error("Tickers are not unique!")

  #batch the tickers for performance
  Ntickers::Int = length(tickers)
  Nperiods::Int = length(periodicity)

  resultparts::Vector{Vector{DataFrame}} = [Vector{DataFrame}(undef, Ntickers) for _ ∈ 1:Nperiods]
  if livepull
    @sync for (i, p) ∈ enumerate(periodicity)
      options::Dict{String, Any} = Dict(
        "CshAdjAbnormal"=>CshAdjAbnormal,
        "CshAdjNormal"=>CshAdjNormal,
        "currency"=>currency,
        #"DateFormat"=>DateFormat,
        )
      if p=="DAILY"
        #options["CDR"] = "7D"
        options["overrideOption"] = "OVERRIDE_OPTION_CLOSE"
        options["nonTradingDayFillOption"] = "ACTIVE_DAYS_ONLY"
        options["nonTradingDayFillMethod"] = "NIL_VALUE"
      end

      #query BB for each partition
      for j ∈ 1:Ntickers
        @async resultparts[i][j] = BLPData.bdh(S, tickers[j], fields, startdate, enddate,
          periodicity=p, options=options) |> DataFrame!
        #=for n ∈ names(resultparts[i][j])
          oldtype = eltype(resultparts[i][j][!,n])
          if !(Missing <: oldtype)
            resultparts[i][j][!,n] = Vector{Union{oldtype,Missing}}=#
      end
    end

    #label each df
    badtickers::Vector{String} = Vector{String}()
    for i ∈ 1:Nperiods
      for j ∈ 1:Ntickers
        if ncol(resultparts[i][j]) > 1#length(setdiff(fields, names(resultparts[i][j]))) == 0
          resultparts[i][j][!,:ticker] .= tickers[j]
        else
          push!(badtickers, "$(tickers[j]):$(periodicity[i])")
          resultparts[i][j] = DataFrame()
        end
      end
    end

    @info "Individual dataframes formed. $(length(badtickers)) bad tickers:\n$(badtickers)"

    #combine the results
    @assert sum([(resultperiod->length(resultperiod)).(resultparts)...;]) == Ntickers*Nperiods
    vcatunion(df1,df2) = vcat(df1,df2, cols=:union)
    results = (resultperiod->reduce(vcatunion, resultperiod)).(resultparts)
    @assert length(results) == Nperiods

    #write out the results
    outbinstream(pullpath, results)
  else
    results = inbinstream(pullpath)
  end

  #println(results)
  #println(results)
  results
end

function writefuturespull(df::AbstractDataFrame, outname::String;
  compress::Bool = size(df,1) > 150_000, #simple heuristic for compressing the CSV
  outputpath::String = PARAM[:outputpath],
  futuresfile::String = PARAM[:capacity][:futuresfile],
  outcsvstream = compress ? OUT_CSV_STREAM : OUT_SMALL_CSV_STREAM,
  csvext = compress ? CSV_EXTENSION : SMALL_CSV_EXTENSION
  )

  p::String = "$outputpath\\$outname.$csvext"
  outcsvstream(p, df)
end

#main entry point
function bbscript(S, ::Val{:futurespull};
  livepull=error("Live pull is required"),
  periodicity=PARAM[:capacity][:periodicity],
  fields=PARAM[:capacity][:pullfields])

  futures::DataFrame = futuresdata()
  tickers::Vector{String}=futures.fullticker

  #=tickers = ["AC1 Comdty", "C 1 Comdty", "DCS1 Comdty", "CZT1 Comdty", "WZ1 Comdty",
  "YW1 Comdty", "CRD1 Comdty", "IBW1 Comdty", "YC1 Comdty", "EP1 Comdty",
  "VB1 Comdty", "VBE1 Comdty", "MIE1 Comdty", "JC1 Comdty", "CBS1 Comdty",
  "JMS1 Comdty", "CXY1 Comdty", "VV1 Comdty", "C71 Comdty", "CT1 Comdty",
  #="QCY1 Comdty", "CTT1 Comdty", "CCL1 Comdty", "SHA1 Comdty", "P71 Comdty",
  "PFC1 Comdty", "PFR1 Comdty", "PAL1 Comdty", "CB1 Comdty", "DCE1 Comdty",
  "SB1 Comdty", "APW1 Comdty", "KO1 Comdty", "CC1 Comdty", "KC1 Comdty",
  "CFD1 Comdty", "DF1 Comdty", "REE1 Comdty", "QC1 Comdty", "QW1 Comdty",
  "SBT1 Comdty", "JCI1 Comdty", "UDR1 Comdty", "C21 Comdty", "KCT1 Comdty",
  "S91 Comdty", "CCT1 Comdty", "POC1 Comdty", "DA1 Comdty", "PMA1 Comdty",
  "JO1 Comdty", "MMR1 Comdty", "CHE1 Comdty", "LE1 Comdty", "V61 Comdty",
  "SF1 Comdty", "AX1 Comdty", "KV1 Comdty", "FEP1 Comdty", "MKP1 Comdty",
  "QSR1 Comdty", "OMP1 Comdty", "CPI1 Comdty", "POI1 Comdty", "CSO1 Comdty",
  "FSP1 Comdty", "BSC1 Comdty", "CSE1 Comdty", "QCR1 Comdty", "P81 Comdty",
  "MDS1 Comdty", "BUT1 Comdty", "JOT1 Comdty", "DRW1 Comdty", "QCE1 Comdty",
  "WWC1 Comdty", "WQS1 Comdty", "SAA1 Comdty", "LC1 Comdty", "LH1 Comdty",
  "FC1 Comdty", "JHT1 Comdty", "LS1 Comdty", "LET1 Comdty", "TFE1 Comdty",
  "OVO1 Comdty", "ZRR1 Comdty", "ZRO1 Comdty", "ANA1 Comdty", "M11 Comdty",
  "RS1 Comdty", "C01 Comdty", "DPR1 Comdty", "Q61 Comdty", "IJ1 Comdty",
  "RR1 Comdty", "O 1 Comdty", "SU1 Comdty", "QRI1 Comdty", "G01 Comdty",
  "FU1 Comdty", "MYDJ1 Comdty", "FY1 Comdty", "AKOJ1 Comdty", "IRL1 Comdty",=#
  #"ZVB1 Comdty", "IY1 Comdty", "IRI1 Comdty", "ZRC1 Comdty", "AE1 Comdty",
  #"SH1 Comdty", "AK1 Comdty", "S 1 Comdty", "SM1 Comdty", "BO1 Comdty",
  #"Q81 Comdty", "BP1 Comdty", "M71 Comdty", "SAR1 Comdty", "SAY1 Comdty",
  #"SAW1 Comdty", "SY1 Comdty", "YK1 Comdty", "FNS1 Comdty", "MSI1 Comdty",
  "VT1 Comdty", "SXY1 Comdty", "LIS1 Comdty", ]=#

  #=results = @btime multipull($S,
    livepull=$livepull,
    tickers=$tickers,
    fields=$fields,
    periodicity=$periodicity,
    pullname="ES1test")=#
  results = multipull(S,
    livepull=livepull,
    tickers=tickers,
    fields=fields,
    periodicity=periodicity,
    pullname="futures")

  results .= (result->leftjoin(result, futures, on=:ticker=>:fullticker)).(results)

  for (i,p) ∈ enumerate(periodicity)
    writefuturespull(results[i], "$(PARAM[:capacity][:futuresfile])_$p")
  end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
function n_tuple!(bs::BlueStyle, fst::FST, s::State)
    style = getstyle(bs)
    line_margin = s.line_offset + length(fst) + fst.extra_margin
    lidx = findlast(n -> n.typ === PLACEHOLDER, fst.nodes)
    fidx = findfirst(n -> n.typ === PLACEHOLDER, fst.nodes)
    multiline_arg = findfirst(is_block, fst.nodes) !== nothing
    multiline_arg && (fst.nest_behavior = AlwaysNest)
    has_closer = is_closer(fst[end])

    # "foo(a, b, c)" is true if "foo" and "c" are on different lines
    src_diff_line = if s.opts.join_lines_based_on_source
        last_arg_idx = findlast(is_iterable_arg, fst.nodes)
        last_arg = last_arg_idx === nothing ? fst[end] : fst[last_arg_idx]
        fst[1].endline != last_arg.startline
    else
        false
    end

    if lidx !== nothing && (line_margin > s.opts.margin || must_nest(fst) || src_diff_line)
        args_range = fidx+1:lidx-1
        args_margin = sum(length.(fst[args_range]))

        nest_to_oneline =
            if can_nest(fst) && (fst.indent + s.opts.indent + args_margin <= s.opts.margin)
                !contains_comment(fst.nodes[args_range])
            else
                false
            end

        # nesting is initiated by source code
        src_nest = src_diff_line && line_margin <= s.opts.margin

        args_diff_line = if src_nest
            # first arg
            first_arg_idx = if fst.typ === TupleN && !has_closer
                1
            elseif is_opener(fst[1])
                # (...)
                findfirst(is_iterable_arg, fst.nodes)
            else
                # f(...)
                findnext(is_iterable_arg, fst.nodes, 2)
            end
            first_arg = first_arg_idx === nothing ? fst[1] : fst[first_arg_idx]
            # last arg
            last_arg_idx = findlast(is_iterable_arg, fst.nodes)
            last_arg = last_arg_idx === nothing ? fst[end] : fst[last_arg_idx]

            first_arg.endline != last_arg.startline
        else
            false
        end

        line_offset = s.line_offset
        if has_closer
            fst[end].indent = fst.indent
        end
        if fst.typ !== TupleN || has_closer
            fst.indent += s.opts.indent
        end

        add_trailing_comma = (src_nest && args_diff_line) || (!src_nest && !nest_to_oneline)

        for (i, n) in enumerate(fst.nodes)
            if n.typ === NEWLINE
                s.line_offset = fst.indent
            elseif (i == fidx || i == lidx) && !src_diff_line && nest_to_oneline
                fst[i] = Newline(length = n.len)
                s.line_offset = fst.indent
            elseif n.typ === PLACEHOLDER
                if src_diff_line
                    si = findnext(
                        n -> n.typ === PLACEHOLDER || n.typ === NEWLINE,
                        fst.nodes,
                        i + 1,
                    )
                    nested = nest_if_over_margin!(style, fst, s, i; stop_idx = si)
                    if has_closer && !nested && n.startline == fst[end].startline
                        # trailing types are automatically converted, undo this if
                        # there is no nest and the closer is on the same in the
                        # original source.
                        if fst[i-1].typ === TRAILINGCOMMA ||
                           fst[i-1].typ === TRAILINGSEMICOLON
                            fst[i-1].val = ""
                            fst[i-1].len = 0
                        end
                    end
                elseif !nest_to_oneline
                    fst[i] = Newline(length = n.len)
                    s.line_offset = fst.indent
                else
                    nest!(style, n, s)
                end
            elseif n.typ === TRAILINGCOMMA
                if add_trailing_comma
                    n.val = ","
                    n.len = 1
                end
                nest!(style, n, s)
            elseif n.typ === TRAILINGSEMICOLON
                n.val = ";"
                n.len = 1
                nest!(style, n, s)
            elseif n.typ === INVERSETRAILINGSEMICOLON
                n.val = ""
                n.len = 0
                nest!(style, n, s)
            elseif has_closer && (i == 1 || i == length(fst.nodes))
                nest!(style, n, s)
            else
                diff = fst.indent - fst[i].indent
                add_indent!(n, s, diff)
                n.extra_margin =
                    !nest_to_oneline ? 1 : i < length(fst.nodes) ? length(fst[i+1]) : 0
                nest!(style, n, s)
            end
        end

        if has_closer
            s.line_offset = fst[end].indent + 1
        end

    else
        extra_margin = fst.extra_margin
        has_closer && (extra_margin += 1)
        nest!(style, fst.nodes, s, fst.indent, extra_margin = extra_margin)
    end
end
@inline n_call!(bs::BlueStyle, fst::FST, s::State) = n_tuple!(bs, fst, s)
@inline n_curly!(bs::BlueStyle, fst::FST, s::State) = n_tuple!(bs, fst, s)
@inline n_macrocall!(bs::BlueStyle, fst::FST, s::State) = n_tuple!(bs, fst, s)
@inline n_ref!(bs::BlueStyle, fst::FST, s::State) = n_tuple!(bs, fst, s)
@inline n_braces!(bs::BlueStyle, fst::FST, s::State) = n_tuple!(bs, fst, s)
@inline n_vect!(bs::BlueStyle, fst::FST, s::State) = n_tuple!(bs, fst, s)
@inline n_parameters!(bs::BlueStyle, fst::FST, s::State) = n_tuple!(bs, fst, s)
@inline n_invisbrackets!(bs::BlueStyle, fst::FST, s::State) = n_tuple!(bs, fst, s)
@inline n_bracescat!(bs::BlueStyle, fst::FST, s::State) = n_tuple!(bs, fst, s)

function n_conditionalopcall!(bs::BlueStyle, fst::FST, s::State)
    style = getstyle(bs)
    if fst[end].typ === Conditional
        conditional_to_if_block!(fst, s)
        nest!(style, fst, s)
    else
        n_conditionalopcall!(DefaultStyle(style), fst, s)
    end
    return
end

function n_binaryopcall!(bs::BlueStyle, fst::FST, s::State)
    style = getstyle(bs)
    if fst.ref !== nothing && parent_is(fst.ref[], n -> is_if(n) || n.head === :macrocall)
        n_binaryopcall!(YASStyle(style), fst, s; indent = fst.indent + s.opts.indent)
    else
        n_binaryopcall!(YASStyle(style), fst, s)
    end
    return
end

function n_chainopcall!(bs::BlueStyle, fst::FST, s::State)
    style = getstyle(bs)
    if fst.ref !== nothing && parent_is(fst.ref[], n -> is_if(n) || n.head === :macrocall)
        n_block!(DefaultStyle(style), fst, s; indent = fst.indent + s.opts.indent)
    else
        n_block!(DefaultStyle(style), fst, s; indent = s.line_offset)
    end
    return
end

n_comparison!(bs::BlueStyle, fst::FST, s::State) = n_chainopcall!(bs, fst, s)


text:
export skipinf, skipnan, skipinfnan
export partitionarray
export pairwisepairs, simpairs

"""
    partitionarray(array::AbstractArray, steps::AbstractVector{<:Integer})
Like Base.Iterators.Partition, except you can apply arbitrary sizes to
partition the array by. The `steps` must add up to the total row length
of the array.

**Example**
```
julia> partitionarray(rand(20,5), [10,3,4,3]) .|> size
[(10, 5), (3, 5), (4, 5), (3, 5)]
```
"""
function partitionarray(array::AbstractArray, steps::AbstractVector{<:Integer})
    # solution brilliantly provided by @stevengj and @mcabbott on Slack and Discourse
    # (https://discourse.julialang.org/t/is-there-a-simple-intuitive-way-to-partition-a-matrix-by-arbitrary-strides-like-i/55863)
    v = axes(array,1)
    v == 1:sum(steps) || error("Steps provided do not sum to length of the first dimension")
    i = firstindex(v)
    tmp = (view(v, i:(i+=s)-1) for s in steps)
    [view(array,r,:) for r in tmp]
end

"""
    pairwisepairs(smp_names::Vector{T}) where T
Given a vector, returns a lazy iterator of tuples of unique all x
all combinations of element pairs, excluding self-comparisons.

**Example**
```
julia> colors = ["red_1", "red_2", "blue_1", "blue_2"] ;

julia> pairwisepairs(colors) |> collect
6-element Array{Tuple{String,String},1}:
 ("red_1", "red_2")
 ("red_1", "blue_1")
 ("red_1", "blue_2")
 ("red_2", "blue_1")
 ("red_2", "blue_2")
 ("blue_1", "blue_2")
```
"""
@inline function pairwisepairs(smp_names::AbstractVector{T}) where T
    len = length(smp_names)
    (tuple(smp_names[i], smp_names[j]) for i in 1:len-1 for j in i+1:len)
end

"""
    _simpairs(data::Vector{String})
Takes a Vector of sample names and returns a Tuple of sample pairs, grouped by simulation
number. This is an internal function used for isolating sibship pairs from simulated shipship
pairs (via `PopGenSims.jl`) to perform `relatedness` estimates only on those pairs.

**Example**
```julia
julia> a = ["sim1_off1", "sim1_off2", "sim2_off1", "sim2_off2"] ;

julia> _simpairs(a)
("sim1_off1", "sim1_off2")
("sim2_off1", "sim2_off2")
```
"""
function simpairs(data::Vector{String})
    n = length(data)
    isodd(n) && throw(ArgumentError("Expected an even number of samples, but got $n"))
    Tuple.(Base.Iterators.partition(sort(data), 2))
end

"""
    skipinf(itr)
Return an iterator over the elements in `itr` skipping `Inf` and `-Inf` values. The returned
object can be indexed using indices of itr if the latter is indexable. Indices
corresponding to `Inf` values are not valid: they are skipped by keys and eachindex,
and a MissingException is thrown when trying to use them. This is effectively `skipmissing`
for `Inf` and `-Inf` values.

Use collect to obtain an `Array` containing the non-`Inf` values in `itr`. Note that even
if `itr` is a multidimensional array, the result will always be a `Vector` since it is not
possible to remove `Inf`s while preserving dimensions of the input.
"""
skipinf(itr) = Iterators.filter(isfinite, itr)

"""
    skipnan(itr)
Return an iterator over the elements in `itr` skipping `NaN` values. The returned
object can be indexed using indices of itr if the latter is indexable. Indices
corresponding to `NaN` values are not valid: they are skipped by keys and eachindex,
and a MissingException is thrown when trying to use them. This is effectively `skipmissing`
for `NaN` values.

Use collect to obtain an `Array` containing the non-`NaN` values in `itr`. Note that even
if `itr` is a multidimensional array, the result will always be a `Vector` since it is not
possible to remove `NaN`s while preserving dimensions of the input.
"""
skipnan(itr) = Iterators.filter(!isnan, itr)

"""
    skipinfnan(itr)
Return an iterator over the elements in `itr` skipping `NaN`, `Inf` and `-Inf` values.
See the docstrings of `skipinf` and `skipnan` more details.
"""
skipinfnan(itr) = Iterators.filter(x -> (isfinite(x) & !isnan(x)), itr)


text:
using Test, PolyChaos

A = [ 0  0  0  0;
 1  0  0  0;
 0  1  0  0;
 0  0  1  0;
 0  0  0  1;
 2  0  0  0;
 1  1  0  0;
 1  0  1  0;
 1  0  0  1;
 0  2  0  0;
 0  1  1  0;
 0  1  0  1;
 0  0  2  0;
 0  0  1  1;
 0  0  0  2]

@test calculateMultiIndices(4, 2) == A
@test calculateMultiIndices(3, 0) == zeros(Int64, 1, 3)
@test_throws DomainError calculateMultiIndices( 3, -1)
@test_throws DomainError calculateMultiIndices( 0,  1)
@test_throws DomainError calculateMultiIndices(-1,  1)

@test PolyChaos.numberPolynomials(3,4) == factorial(3+4) ÷ ( factorial(3)*factorial(4) )
@test findUnivariateIndices(2,A) == [1,3,10]