text:
export Tparam
export PDataIter
export PDataLDLt, PDataK, PDataST

abstract type TPData{T} end  # Ancestor of all PreProcess data

abstract type PDataFact{T} <: TPData{T}  end # Variants using matricial factorization

abstract type PDataIter{T} <: TPData{T}  end # Variants using iterative (Krylov) solvers

abstract type Tparams{T}  end # Ancestor of variant specific parameters

mutable struct Tparam{T} <: Tparams{T} end

mutable struct PDataK{T} <: PDataIter{T}
    d :: Array{T,1}             # (H+λI)\g ; on first call = g
    λ :: T                      # "active" value of λ; on first call = 0
    ζ :: T                      # Inexact Newton order parameter: stop when ||∇q||<||g||^(1+ζ )
    τ :: T                      # temporary testing parameter for decreaseARCqK

    indmin :: Int               # index of best shift value  within "positive". On first call = 0

    positives :: Array{Int,1}   # indices of the shift values yielding (H+λI)⪰0
    xShift :: Array{T,2}        # solutions for each shifted system
    shifts :: Array{T,1}        # values of the shifts
    nshifts :: Int              # number of shifts
    norm_dirs :: Array{T,1}     # norms of xShifts
    OK :: Bool                  # preprocess success
end

mutable struct PDataST{T} <: PDataIter{T}
    H             # Hessian representation, most likely as a linear operator or a sparse matrix
    g             # gradient vector
    ζ :: T        # Inexact Newton order parameter: stop when ||∇q||<||g||^(1+ζ)

    OK :: Bool    # preprocess success
end

mutable struct PDataLDLt{T} <: PDataFact{T}
    L        :: Array{T,2}          # could be sparse
    D        :: Array{T,2}          # block diagonal 1X1 and 2X2
    pp       :: Array{Int,1}        # permutation vector LDL' = H[pp,pp]
                                    # P*v = v[pp]
                                    # P'*v = v[invperm(pp)]
    #s::Array{T,1}                  # Scaling matrix S=diagm(s)
    Δ        :: Array{T,1}          # diagonal, eigenvalues of D
    Q        :: Array{T,2}          # orthogonal matrix, eigenvectors of D:  should be sparse
                                    # QΔQ'  =  D
    g̃        ::Array{T,1}           # transformed gradient
    λ        ::T
    success  ::Bool                 # previous iteration was successfull
    OK       ::Bool                 # preprocess success

    PDataLDLt() = new{Nothing}()
    PDataLDLt(L, D, pp, Δ, Q, g, l, success, OK) = new{eltype(L)}(L, D, pp, Δ, Q, g, l, success, OK)
end

mutable struct PDataSpectral{T} <: PDataFact{T}
    V::Array{T,2}           # orthogonal matrix, eigenvectors of Hessian matrix
    Δ::Array{T,1}           # eigenvalues of Hessian matrix
    g̃::Array{T,1}           # transformed gradient
    λ::T
    success::Bool           # previous iteration was successfull
    OK::Bool                # preprocess success

    PDataSpectral() = new{T}()
    PDataSpectral(V, Λ, g, l, success, OK) = new{eltype(V)}(V, Λ, g, l, success, OK)
end


text:
module Evaluate
    using ..Models: Mixture
    using ..Fitter: prepare_sufficient_statistics, update_cluster_sufficient_statistics!, update_data_sufficient_statistics!
    using ..Fitter: compute_data_posterior_parameters!, get_clusters, compute_existing_cluster_log_predictives, new_cluster_log_predictive
    using ..Fitter: compute_data_log_predictives, prepare_auxillary_variables, gibbs_sample!
    using ..Utils: compute_normalized_weights, generate_all_clusterings
    import ..Fitter: compute_joint_log_likelihood

    export compute_posterior_probability
    export compute_posterior_probability_confidence_interval

    using Statistics

    function update_sufficient_statistics!(sufficient_stats, observation, cluster, model, data, aux_var; update_type="add")
        update_cluster_sufficient_statistics!(sufficient_stats, cluster, update_type=update_type)
        update_data_sufficient_statistics!(sufficient_stats, cluster, observation, model, data, aux_var, update_type=update_type)
        compute_data_posterior_parameters!(cluster, sufficient_stats, model.data_parameters)
    end

    function compute_joint_log_likelihood(clustering::Vector{Int64}, data::Matrix{T}, model::Mixture) where {T <: Real}
        sufficient_stats = prepare_sufficient_statistics(model, data)
        aux_variables = prepare_auxillary_variables(model, data)

        for (observation, cluster) in enumerate(clustering)
            update_sufficient_statistics!(sufficient_stats, observation, cluster, model, data, aux_variables)
        end
        log_likelihood = compute_joint_log_likelihood(sufficient_stats, model, aux_variables)
        return log_likelihood
    end

    function compute_joint_log_likelihood(clustering::Vector{Vector{Int64}}, data::Matrix{T}, model::Mixture) where {T <: Real}
        sufficient_stats = prepare_sufficient_statistics(model, data)
        aux_variables = prepare_auxillary_variables(model, data)

        for cluster in clustering
            cluster_id = minimum(cluster)
            for observation in cluster
                update_sufficient_statistics!(sufficient_stats, observation, cluster_id, model, data, aux_variables)
            end
        end
        log_likelihood = compute_joint_log_likelihood(sufficient_stats, model, aux_variables)
        return log_likelihood
    end

    function compute_all_clustering_log_likelihoods(n, data::Matrix{T}, model::Mixture) where {T <: Real}
        all_clusterings = generate_all_clusterings(n)
        clustering_log_probabilities = Vector{Float64}(undef, length(all_clusterings))
        for (i, clustering) = enumerate(all_clusterings)
            clustering_log_probabilities[i] = compute_joint_log_likelihood(clustering, data, model)
        end
        return clustering_log_probabilities
    end

    function compute_posterior_probability(clustering::Vector{Int64}, data::Matrix{T}, model::Mixture) where {T <: Real}
        n = length(clustering)

        all_clusterings_log_probs = compute_all_clustering_log_likelihoods(n, data, model)
        denominator = sum(exp.(all_clusterings_log_probs))

        log_numerator = compute_joint_log_likelihood(clustering, data, model)

        log_probability = log_numerator - log(denominator)
        return exp(log_probability)
    end

    function compute_posterior_probability(clustering::Vector{Vector{Int64}}, data::Matrix{T}, model::Mixture) where {T <: Real}
        n = size(data)[2]
        clustering_vector = Vector{Int64}(undef, n)
        for cluster in clustering
            cluster_id = minimum(cluster)
            for observation in cluster
                clustering_vector[observation] = cluster_id
            end
        end
        return compute_posterior_probability(clustering_vector, data, model)
    end

    function compute_mcmc_se(clustering::Vector{Int64}, markov_chain::Matrix{Int64})
        chain_length = size(markov_chain)[2]
        batch_length = floor(Int64, sqrt(chain_length))
        num_batches = ceil(Int64, sqrt(chain_length))

        batch_probabilities = Vector{Float64}(undef, num_batches)
        for batch = 1:num_batches
            start_index = batch_length*(batch-1) + 1
            stop_index = batch_length*batch
            batch_probabilities[batch] = compute_posterior_probability(clustering, markov_chain[:, start_index:stop_index])
        end
        posterior_probability = mean(batch_probabilities)
        sigma_squared = (batch_length/(num_batches-1))*sum((batch_probabilities .- posterior_probability).^2)
        mcmc_se = sqrt(sigma_squared)/sqrt(chain_length)
        return mcmc_se
    end

    function compute_mcmc_se(clustering::Vector{Vector{Int64}}, markov_chain::Matrix{Int64})
        n = size(markov_chain)[1]
        clustering_vector = Vector{Int64}(undef, n)
        for cluster in clustering
            cluster_id = minimum(cluster)
            for observation in cluster
                clustering_vector[observation] = cluster_id
            end
        end
        return compute_mcmc_se(clustering_vector, markov_chain)
    end

    function compute_posterior_probability(clustering::Vector{Int64}, markov_chain::Matrix{Int64})
        chain_length = size(markov_chain)[2]
        num_iterations_equal = count(vec(mapslices(z->(clustering == vec(z)), markov_chain, dims=1)))
        posterior_probability = num_iterations_equal/chain_length
        return posterior_probability
    end

    function compute_posterior_probability(clustering::Vector{Vector{Int64}}, markov_chain::Matrix{Int64})
        n = size(markov_chain)[1]
        clustering_vector = Vector{Int64}(undef, n)
        for cluster in clustering
            cluster_id = minimum(cluster)
            for observation in cluster
                clustering_vector[observation] = cluster_id
            end
        end
        return compute_posterior_probability(clustering_vector, markov_chain)
    end

    function compute_posterior_probability_confidence_interval(clustering, markov_chain::Matrix{Int64})
        posterior_probability = compute_posterior_probability(clustering, markov_chain)
        mcmc_se = compute_mcmc_se(clustering, markov_chain)
        return posterior_probability, 1.96*mcmc_se
    end

end

ext: jl

max_stars_repo_name: realseanla/ntl-stick-breaking-julia

max_issues_repo_name: realseanla/ntl-stick-breaking-julia

max_forks_repo_name: realseanla/ntl-stick-breaking-julia


text:
using PyPlot, DelimitedFiles, Random, Revise,
      HiQGA.transD_GP
##
Random.seed!(23)
rMin = 500 #m
rMax = 5200  #m
nRx  = 32
freqs   = [0.25, 0.75, 1.75, 3.25] #Hz
rRx    = collect(LinRange(rMin,rMax,nRx))  # Ranges to receivers   (m)
zRx    = 1000-eps(1.)
zTx    = 975.
RxAzim = 0.
TxDip  = 0.
# Note that the receiver depth needs to be in same model layer as transmitter.
zfixed = [-1e6,    0      ]
ρfixed = [1e13,    0.3    ]
zstart = zRx + eps(1.)
extendfrac, dz = 1.008, 12.
zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac, dz=dz, n=126)
z, ρ, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed)
##
F = transD_GP.CSEM1DEr.RadialErLagged(zTx    = [zTx],
                      rRx    = rRx,
                      freqs  = freqs,
                      zRx    = [zRx],
                      RxAzim = RxAzim,
                      TxDip  = TxDip)
##
ρ[(z.>=zstart) .& (z.<1400)] .= 10^0.0
ρ[(z.>=1400) .&(z.<1700)] .= 10^0.3
ρ[(z.>=1700) .&(z.<1800)] .= 10^0.5
ρ[(z.>=1800) .&(z.<1900)] .= 10^0.2
ρ[(z.>=1900) .&(z.<2000)] .= 10^1.4 # 1.4 or 0.5
ρ[(z.>=2000) .&(z.<2400)] .= 10^0.0
ρ[(z.>=2400) .&(z.<2500)] .= 10^0.3
ρ[(z.>=2500) .&(z.<3500)] .= 10^0.7
ρ[(z.>=3500)]              .= 10^1.0
##
transD_GP.CSEM1DInversion.plotmodelfield!(F, z, ρ)
d, σ = transD_GP.addnoise(F, z, ρ, noisefrac=0.05,
               dz=dz, extendfrac=extendfrac, nfixed=nfixed)
##
csem = transD_GP.CSEMRadialEr(F, d, σ, z=z, ρ=ρ, nfixed=nfixed)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
"""
    abstract type NonlinearOperator <: GridapType end

- [`residual!(b::AbstractVector,op::NonlinearOperator,x::AbstractVector)`](@ref)
- [`jacobian!(A::AbstractMatrix,op::NonlinearOperator,x::AbstractVector)`](@ref)
- [`zero_initial_guess(op::NonlinearOperator)`](@ref)
- [`allocate_residual(op::NonlinearOperator,x::AbstractVector)`](@ref)
- [`allocate_jacobian(op::NonlinearOperator,x::AbstractVector)`](@ref)

"""
abstract type NonlinearOperator <: GridapType end

"""
    residual!(b::AbstractVector,op::NonlinearOperator,x::AbstractVector)
"""
function residual!(b::AbstractVector,op::NonlinearOperator,x::AbstractVector)
  @abstractmethod
end

"""
    residual(op::NonlinearOperator,x::AbstractVector)
"""
function residual(op::NonlinearOperator,x::AbstractVector)
  b = allocate_residual(op,x)
  residual!(b,op,x)
  b
end

"""
    jacobian!(A::AbstractMatrix,op::NonlinearOperator,x::AbstractVector)
"""
function jacobian!(A::AbstractMatrix,op::NonlinearOperator,x::AbstractVector)
  @abstractmethod
end

"""
    jacobian(op::NonlinearOperator,x::AbstractVector)
"""
function jacobian(op::NonlinearOperator,x::AbstractVector)
  A = allocate_jacobian(op,x)
  jacobian!(A,op,x)
  A
end

"""
    residual_and_jacobian!(
      b::AbstractVector, A::AbstractMatrix,
      op::NonlinearOperator, x::AbstractVector)
"""
function residual_and_jacobian!(
  b::AbstractVector, A::AbstractMatrix, op::NonlinearOperator, x::AbstractVector)
  residual!(b,op,x)
  jacobian!(A,op,x)
  (b, A)
end

"""
    residual_and_jacobian(op::NonlinearOperator,x::AbstractVector)
"""
function residual_and_jacobian(op::NonlinearOperator,x::AbstractVector)
  b, A = allocate_residual_and_jacobian(op,x)
  residual_and_jacobian!(b,A,op,x)
  (b, A)
end

"""
    zero_initial_guess(op::NonlinearOperator)
"""
function zero_initial_guess(op::NonlinearOperator)
  @abstractmethod
end

"""
    allocate_residual(op::NonlinearOperator,x::AbstractVector)
"""
function allocate_residual(op::NonlinearOperator,x::AbstractVector)
  @abstractmethod
end

"""
    allocate_jacobian(op::NonlinearOperator,x::AbstractVector)
"""
function allocate_jacobian(op::NonlinearOperator,x::AbstractVector)
  @abstractmethod
end

"""
    allocate_residual_and_jacobian(op::NonlinearOperator,x::AbstractVector)
"""
function allocate_residual_and_jacobian(op::NonlinearOperator,x::AbstractVector)
  b = allocate_residual(op,x)
  A = allocate_jacobian(op,x)
  (b,A)
end

"""
    test_nonlinear_operator(
      op::NonlinearOperator,
      x::AbstractVector,
      b::AbstractVector,
      pred=isapprox;
      jac=nothing)
"""
function test_nonlinear_operator(
  op::NonlinearOperator,
  x::AbstractVector,
  b::AbstractVector,
  pred=isapprox;
  jac=nothing)

  b1 = allocate_residual(op,x)
  residual!(b1,op,x)
  @test pred(b,b1)

  x0 = zero_initial_guess(op)

  if jac != nothing
    nrows, ncols = size(jac)
    A = allocate_jacobian(op,x)
    jacobian!(A,op,x)
    @test pred(A,jac)

    residual_and_jacobian!(b1,A,op,x)
    @test pred(A,jac)
    @test pred(b,b1)

    b1, A = residual_and_jacobian(op,x)
    @test pred(A,jac)
    @test pred(b,b1)
  end

end

# Mock for testing purposes

struct NonlinearOperatorMock <: NonlinearOperator end

function residual!(b::AbstractVector,::NonlinearOperatorMock,x::AbstractVector)
  b[1] = (x[1]-2)*(x[1]-1)
  b[2] = x[2]-3
  b
end

function jacobian!(A::AbstractMatrix,::NonlinearOperatorMock,x::AbstractVector)
  A[1,1] = (x[1]-1) + (x[1]-2)
  A[1,2] = 0
  A[2,1] = 0
  A[2,2] = 1
  A
end

function zero_initial_guess(op::NonlinearOperatorMock)
  x = allocate_residual(op,Float64[])
  fill!(x,zero(eltype(x)))
  x
end

function allocate_residual(op::NonlinearOperatorMock,x::AbstractVector)
  T = eltype(x)
  n = 2
  similar(x,T,n)
end

function allocate_jacobian(op::NonlinearOperatorMock,x::AbstractVector)
  T = eltype(x)
  n = 2
  m = 2
  zeros(T,n,m)
end


text:
#########################
#    SET DIRECTORIES    #
#########################

# load support functions
include("support_functions.jl")

####################
#    STRUCTURES    #
####################

struct AttributeTable

    fp_table::String
    key::Symbol
    fields_to_dict::Array{Symbol, 1}
    clean_table_fields::Bool
    table
    field_maps
    key_vals

    function AttributeTable(fp_table::String, key::Symbol, fields_to_dict::Array{Symbol, 1}, clean_table_fields::Bool)

        # verify table exists and check keys
        table = read_csv(check_path(fp_table, false), true)
        fields_to_dict = [x for x in fields_to_dict if x != key]
        check_fields!(table, [[key]; fields_to_dict])

        # check key
        if length(Set(table[:, key])) < nrow(table)
            error("Invalid key '$(key)' found in $(fp_table): the key is not unique. Check the table and specify a unique key.")
        end

        # clean the keys?
        if clean_table_fields
            clean_field_names!(table)
            fields_to_dict = clean_field_names(fields_to_dict)
            key = clean_field_names([key])[1]
        end

        # clear RST formatting
        table[!, key] = replace.(replace.(table[!, key], "`" => ""), "\$" => "")

        # next, create dict maps
        field_maps = Dict()
        for fld in fields_to_dict
            field_fwd = "$(key)_to_$(fld)"
            field_rev = "$(fld)_to_$(key)"

            field_maps[field_fwd] = build_dict(table[:, [key, fld]])
            # check for 1:1 correspondence before adding reverse
            vals_unique = Set(table[:, fld])
            if (length(vals_unique) == nrow(table))
                field_maps[field_rev] = build_dict(table[:, [fld, key]])
            end
        end

        # get all values associated with they key
        key_vals = sort(collect(Set(table[:, key])))

        return new(fp_table, key, fields_to_dict, clean_table_fields, table, field_maps, key_vals)
    end
end

## high level directory structure
dir_proj = dirname(dirname(@__FILE__))
#dir_data = check_path(joinpath(dirname(dir_proj), "data"), false)
dir_bin = joinpath(dir_proj, "bin")
dir_out = check_path(joinpath(dir_proj, "out"), true)
dir_ref = check_path(joinpath(dir_proj, "ref"), false)


text:
module D3GraphPlot

using LightGraphs, JSON, Blink, Atom, Juno

include("d3graphelectron.jl")
include("d3graphjuno.jl")

function graphtojsonstring(g::Graph)
    data = Dict()

    data["nodes"] = []
    for v in 1:nv(g)
        push!(data["nodes"], Dict("name" => v))
    end

    data["edges"] = []
    for e in edges(g)
        push!(data["edges"], Dict("source" => (src(e)-1), "target" => (dst(e)-1)))
    end

    data_json_string = JSON.json(data)
    # print(data_json_string)
    return data_json_string
end

function d3BlinkPlot(g::Graph)
    #load javascript function drawgraph() that is based on d3
    url = string("file://", dirname(@__FILE__), "/../resources/index.html")
    w = Window()
    loadurl(w, url)

    #prepare the dataset to be passed to drawgraph()
    data_json_string = graphtojsonstring(g)

    #call the javascript function drawgraph()
    @js w (drawgraph($data_json_string);)
end

function d3plot(g::Graph)
    #load javascript function drawgraph() that is based on d3 in the pane Plots
    url = string("file://", dirname(@__FILE__), "/../resources/index.html")

    #prepare the dataset to be passed to drawgraph()
    data_json_string = graphtojsonstring(g)

    #call the javascript function drawgraph()
    plot = Plot(data_json_string, Base.Random.uuid1())

    #display the plot
    display_blink(SyncPlot(plot))

    return plot
end

export d3BlinkPlot, d3plot

end # module


text:
abstract type DelayDiffEqAlgorithm <: AbstractDDEAlgorithm end
abstract type AbstractMethodOfStepsAlgorithm{constrained} <: DelayDiffEqAlgorithm end

struct MethodOfSteps{algType,F,constrained} <: AbstractMethodOfStepsAlgorithm{constrained}
  alg::algType
  fpsolve::F
end

"""
    MethodOfSteps(alg; constrained = false, fpsolve = NLFunctional())

Construct an algorithm that solves delay differential equations by the method of steps,
where `alg` is an ODE algorithm from OrdinaryDiffEq.jl upon which the calculation of
steps is based.

If the algorithm is `constrained` only steps of size at most the minimal delay will be
taken. If it is unconstrained, fixed-point iteration `fpsolve` is applied for step sizes
that exceed the minimal delay.

Citations:

General Approach

ZivariPiran, Hossein, and Wayne H. Enright. "An efficient unified approach for the numerical
solution of delay differential equations." Numerical Algorithms 53.2-3 (2010): 397-417.

State-Dependent Delays

S. P. Corwin, D. Sarafyan and S. Thompson in "DKLAG6: a code based on continuously imbedded
sixth-order Runge-Kutta methods for the solution of state-dependent functional differential
equations", Applied Numerical Mathematics, 1997.
"""
MethodOfSteps(alg; constrained = false, fpsolve = NLFunctional()) =
  MethodOfSteps{typeof(alg),typeof(fpsolve),constrained}(alg, fpsolve)

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
"""
    main_generate_uncertainties_custom
A main file allowing to generate uncertainties for a pscopf instance.
This provides an example to customize the uncertainties generation method
 by overloading the function `PSCOPF.generate_values`.

Parameters:
data_path : the path to network files (pscopf_limits, pscopf_units, pscopf_ptdf) and pscopf_init
out_path : path where pscopf_uncertainties will be written
copy_data_files : If true, input files will be rewritten to out_path
nb_scenarios : Number of uncertain scenarios generated
TS : interest timepoints for which an uncertain value is generated
uncertainties_distribution : The distribution of the injections uncertainties (c.f. uncertainties_distribution.txt)

To define a custom function for generating uncertain values, One can overload the function :
PSCOPF.generate_values(uncertainty_distribution,
                        ech::Dates.DateTime, ts::Dates.DateTime,
                        nb_scenarios::Int64)
"""

using Dates
using Distributions
using DataStructures

root_path = dirname(@__DIR__)
push!(LOAD_PATH, root_path);
cd(root_path)
include(joinpath(root_path, "src", "PSCOPF.jl"));

#########################
# INPUT & PARAMS
#########################

# data_path is the path used to read the network and generators initial state:
# pscopf_units : the description of the units
# pscopf_limits : the flow capacity of each branch
# pscopf_ptdf : the ptdf coefficients per (branch, bus_id)
# pscopf_init : the description of the initial state of units (just before the first interest timepoit T)
data_path = joinpath(@__DIR__, "..", "data", "2buses_small_usecase")

# out_path is the path where uncertainties file will be written:
# pscopf_uncertainties : the randomly generated nodal injections
out_path = joinpath(data_path, "custom_uncertainties")

# If true network and generator initial state files will be written to the out_path
copy_data_files = true

# Number of scenarios to be generated
nb_scenarios = 5

# interest time points (ie T)
ts1 = Dates.DateTime("2015-01-01T11:00:00")
TS = PSCOPF.create_target_timepoints(ts1) #T: 11h, 11h15, 11h30, 11h45

# uncertainties_distribution :
#distribution parameters used to generate injections uncertainties
uncertainties_distribution = SortedDict("wind_1" => [ 0, 10],
                                        "wind_2" => [10, 20],
                                        "bus_1" =>  [20, 50],
                                        "bus_2" =>  [50, 70])

function PSCOPF.generate_values(uncertainty_distribution,
                        ech::Dates.DateTime, ts::Dates.DateTime,
                        nb_scenarios::Int64)
    return rand(Uniform(uncertainty_distribution[1], uncertainty_distribution[2]), nb_scenarios)
end

#########################
# EXECUTION
#########################
function generate_uncerainties(output_path, input_path,
                                uncertainties_distribution, nb_scenarios,
                                TS::Vector{DateTime}, #T
                                ECH::Vector{DateTime}=Vector{DateTime}(); #ECH
                                write_data::Bool=true
                                )
    # load network
    println("load network")
    network = PSCOPF.Data.pscopfdata2network(input_path)

    #define ECH (ie t)
    if isempty(ECH)
        ECHs = []
        for mode in [PSCOPF.PSCOPF_MODE_1, PSCOPF.PSCOPF_MODE_2, PSCOPF.PSCOPF_MODE_3]
            push!(ECHs, PSCOPF.generate_ech(network, TS, mode) )
        end
        horizon_timepoints = sort(unique(Iterators.flatten(ECHs)))
    else
        horizon_timepoints = ECH
    end
    println("Uncertainties will be generated for horizon timepoints (ie. t) : ", horizon_timepoints)

    # generate uncertainties
    uncertainties = PSCOPF.generate_uncertainties(network, TS, horizon_timepoints,
                                                uncertainties_distribution, nb_scenarios)
    println("Uncertainties generated!")

    # write files
    PSCOPF.PSCOPFio.write(output_path, uncertainties)
    if write_data
        PSCOPF.PSCOPFio.write(output_path, network)
        generators_init_state = PSCOPF.PSCOPFio.read_initial_state(input_path)
        PSCOPF.PSCOPFio.write(output_path, generators_init_state)
    end
    println("Uncertainties written to ", output_path)
end

generate_uncerainties(out_path, data_path,
                        uncertainties_distribution, nb_scenarios,
                        TS,
                        write_data=copy_data_files)

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
module SDSS

import DataFrames
import FITSIO
import Grid
import WCS

import ..PSF

const band_letters = ['u', 'g', 'r', 'i', 'z']

const default_mask_planes = Set(
    ["S_MASK_INTERP", "S_MASK_SATUR", "S_MASK_CR", "S_MASK_GHOST"])

function load_stamp_catalog_df(cat_dir, stamp_id, blob; match_blob=false)
    # These files are generated by
    # https://github.com/dstndstn/tractor/blob/master/projects/inference/testblob2.py
    cat_fits = FITSIO.FITS("$cat_dir/cat-$stamp_id.fits")
    num_cols = FITSIO.read_key(cat_fits[2], "TFIELDS")[1]
    ttypes = [FITSIO.read_key(cat_fits[2], "TTYPE$i")[1] for i in 1:num_cols]

    df = DataFrames.DataFrame()
    for i in 1:num_cols
        tmp_data = read(cat_fits[2], ttypes[i])
        df[symbol(ttypes[i])] = tmp_data
    end

    close(cat_fits)

    if match_blob
        camcol_matches = df[:camcol] .== blob[3].camcol_num
        run_matches = df[:run] .== blob[3].run_num
        field_matches = df[:field] .== blob[3].field_num
        df = df[camcol_matches & run_matches & field_matches, :]
    end

    df
end

"""
Load data from a psField file, which contains the point spread function.

Args:
 - field_dir: The directory of the file
 - run_num: The run number
 - camcol_num: The camcol number
 - field_num: The field number
 - b: The filter band (a number from 1 to 5)

Returns:
 - RawPSFComponents.

The point spread function is represented as an rnrow x rncol image
showing what a true point source would look as viewed through the optics.
This image varies across the field, and to parameterize this variation,
the PSF is represented as a linear combination of four "eigenimages",
with weights that vary across the image.  See get_psf_at_point()
for more details.
"""
function load_psf_data(field_dir, run_num, camcol_num, field_num, b)
    @assert 1 <= b <= 5
    psf_filename = "$field_dir/psField-$run_num-$camcol_num-$field_num.fit"
    psf_fits = FITSIO.FITS(psf_filename)
    psf_hdu = psf_fits[b + 1]

    nrows = FITSIO.read_key(psf_hdu, "NAXIS2")[1]
    nrow_b = read(psf_hdu, "nrow_b")[1]
    ncol_b = read(psf_hdu, "ncol_b")[1]
    rnrow = read(psf_hdu, "rnrow")[1]
    rncol = read(psf_hdu, "rncol")[1]
    cmat = convert(Array{Float64, 3}, read(psf_hdu, "c"))
    rrows = convert(Array{Float64, 2}, reduce(hcat, read(psf_hdu, "rrows")))
    close(psf_fits)

    # Only the first (nrow_b, ncol_b) submatrix of cmat is used for reasons obscure
    # to the author.
    PSF.RawPSFComponents(rrows, rnrow, rncol, cmat[1:nrow_b, 1:ncol_b, :])
end

"""
Load relevant data from a photoField file.

Args:
 - field_dir: The directory of the file
 - run_num: The run number
 - camcol_num: The camcol number
 - field_num: The field number
 - b: The filter band id (a number from 1 to 5)

Returns:
 - band_gain: An array of gains for the five bands
 - band_dark_variance: An array of dark variances for the five bands
"""
function load_photo_field(field_dir, run_num, camcol_num, field_num)
    pf_filename = "$field_dir/photoField-$run_num-$camcol_num.fits"
    pf_fits = FITSIO.FITS(pf_filename)
    @assert length(pf_fits) == 2

    field_num_int = parse(Int, field_num)
    field_row = read(pf_fits[2], "field") .== field_num_int
    band_gain = collect(read(pf_fits[2], "gain")[:, field_row])
    band_dark_variance = collect(read(pf_fits[2], "dark_variance")[:, field_row])

    close(pf_fits)

    band_gain, band_dark_variance
end

"""
Load the raw electron counts, calibration vector, and sky background from a field.

Args:
 - field_dir: The directory of the file
 - run_num: The run number
 - camcol_num: The camcol number
 - field_num: The field number
 - b: The filter band (a number from 1 to 5)
 - gain: The gain for this band (e.g. as read from photoField)

Returns:
 - nelec: An image of raw electron counts in nanomaggies
 - calib_col: A column of calibration values (the same for every column of the image)
 - sky_grid: A CoordInterpGrid bilinear interpolation object
 - sky_x: The x coordinates at which to evaluate sky_grid to match nelec.
 - sky_y: The y coordinates at which to evaluate sky_grid to match nelec.
 - sky_image: The sky interpolated to the original image size.
 - wcs: A WCS.WCSTransform object for convert between world and pixel
   coordinates.

The meaing of the frame data structures is thoroughly documented here:
http://data.sdss3.org/datamodel/files/BOSS_PHOTOOBJ/frames/RERUN/RUN/CAMCOL/frame.html
"""
function load_raw_field(field_dir, run_num, camcol_num, field_num, b, gain)
    @assert 1 <= b <= 5
    b_letter = band_letters[b]

    img_filename = "$field_dir/frame-$b_letter-$run_num-$camcol_num-$field_num.fits"
    img_fits = FITSIO.FITS(img_filename)
    @assert length(img_fits) == 4

    # This is the sky-subtracted and calibrated image.
    processed_image = read(img_fits[1])

    # Read in the sky background.
    sky_image_raw = read(img_fits[3], "ALLSKY")
    sky_x = collect(read(img_fits[3], "XINTERP"))
    sky_y = collect(read(img_fits[3], "YINTERP"))

    # Get the WCS coordinates.
    header_str = FITSIO.read_header(img_fits[1], ASCIIString)
    wcs = WCS.from_header(header_str)[1]

    # These are the column types (not currently used).
    ctype = [FITSIO.read_key(img_fits[1], "CTYPE1")[1],
             FITSIO.read_key(img_fits[1], "CTYPE2")[1]]

    # This is the calibration vector:
    calib_col = read(img_fits[2])
    calib_image = [ calib_col[row] for
                    row in 1:size(processed_image)[1],
                    col in 1:size(processed_image)[2] ]

    close(img_fits)

    # Interpolate the sky to the full image.  Combining the example from
    # http://data.sdss3.org/datamodel/files/BOSS_PHOTOOBJ/frames/RERUN/RUN/CAMCOL/frame.html
    # ...with the documentation from the IDL language:
    # http://www.exelisvis.com/docs/INTERPOLATE.html
    # ...we can see that we are supposed to interpret a point (sky_x[i], sky_y[j])
    # with associated row and column (if, jf) = (floor(sky_x[i]), floor(sky_y[j]))
    # as lying in the square spanned by the points
    # (sky_image_raw[if, jf], sky_image_raw[if + 1, jf + 1]).
    # ...keeping in mind that IDL uses zero indexing:
    # http://www.exelisvis.com/docs/Manipulating_Arrays.html
    sky_grid_vals = ((1:1.:size(sky_image_raw)[1]) - 1, (1:1.:size(sky_image_raw)[2]) - 1)
    sky_grid = Grid.CoordInterpGrid(sky_grid_vals, sky_image_raw[:,:,1],
                                    Grid.BCnearest, Grid.InterpLinear)

    # This interpolation is really slow.
    sky_image = [ sky_grid[x, y] for x in sky_x, y in sky_y ]

    # Convert to raw electron counts.  Note that these may not be close to integers
    # due to the analog to digital conversion process in the telescope.
    nelec = gain * convert(Array{Float64, 2}, (processed_image ./ calib_image .+ sky_image))

    nelec, calib_col, sky_grid, sky_x, sky_y, sky_image, wcs
end

"""
Set the pixels in mask_img to NaN in the places specified by the fpM file.

Args:
 - mask_img: The image to be masked (updated in place)
 - field_dir: The directory of the file
 - run_num: The run number
 - camcol_num: The camcol number
 - field_num: The field number

Returns:
 - Updates mask_img in place by setting to NaN all the pixels specified.

 This is based on the function setMaskedPixels in astrometry.net:
 https://github.com/dstndstn/astrometry.net/
"""
function mask_image!(mask_img, field_dir, run_num, camcol_num, field_num, band;
                     python_indexing = false,
                     mask_planes = default_mask_planes)
    # The default mask planes are those used by Dustin's astrometry.net code.
    # See the comments in sdss/dr8.py for fpM.setMaskedPixels
    # and the function sdss/common.py:fpM.setMaskedPixels
    #
    # interp = pixel was bad and interpolated over
    # satur = saturated
    # cr = cosmic ray
    # ghost = artifact from the electronics.

    # http://data.sdss3.org/datamodel/files/PHOTO_REDUX/RERUN/RUN/objcs/CAMCOL/fpM.html
    band_letter = band_letters[band]
    fpm_filename = "$field_dir/fpM-$run_num-$band_letter$camcol_num-$field_num.fit"
    fpm_fits = FITSIO.FITS(fpm_filename)

    # The last header contains the mask.
    fpm_mask = fpm_fits[12]
    fpm_hdu_indices = read(fpm_mask, "value")

    # Only these rows contain masks.
    masktype_rows = find(read(fpm_mask, "defName") .== "S_MASKTYPE")

    # Apparently attributeName lists the meanings of the HDUs in order.
    mask_types = read(fpm_mask, "attributeName")
    plane_rows = findin(mask_types[masktype_rows], mask_planes)

    # Make sure each mask is present.  Is this check appropriate for all mask files?
    @assert length(plane_rows) == length(mask_planes)

    for fpm_i in plane_rows
        # You want the HDU in 2 + fpm_mask.value[i] for i in keep_rows (in a 1-indexed language).
        mask_index = 2 + fpm_hdu_indices[fpm_i]
        cmin = read(fpm_fits[mask_index], "cmin")
        cmax = read(fpm_fits[mask_index], "cmax")
        rmin = read(fpm_fits[mask_index], "rmin")
        rmax = read(fpm_fits[mask_index], "rmax")
        row0 = read(fpm_fits[mask_index], "row0")
        col0 = read(fpm_fits[mask_index], "col0")

        @assert all(col0 .== 0)
        @assert all(row0 .== 0)
        @assert length(rmin) == length(cmin) == length(rmax) == length(cmax)

        for block in 1:length(rmin)
            # The ranges are for a 0-indexed language.
            @assert cmax[block] + 1 <= size(mask_img)[1]
            @assert cmin[block] + 1 >= 1
            @assert rmax[block] + 1 <= size(mask_img)[2]
            @assert rmin[block] + 1 >= 1

            # Some notes:
            # See astrometry.net//sdss/common.py:SetMaskedPixels, which I believe
            # is currently incorrect.
            # - In contrast with  julia, the numpy matrix index range [3:5, 3:5] contains four
            #   pixels, not six.  However, if the numpy is correct, then fpM files contain
            #   many bad rows that don't get masked at all since cmin == cmax
            #   or rmin == rmax.  For this reason, I think the python might be erroneous.
            # - For some reason, the sizes are inconsistent if the rows are read first.
            #   I presume that either these names are strange or I am supposed to read
            #   the image from the field and transpose it.
            # - Julia is 1-indexed, not 0-indexed.

            # Give the option of using Dustin's python indexing or not.
            if python_indexing
                mask_rows = (cmin[block] + 1):(cmax[block])
                mask_cols = (rmin[block] + 1):(rmax[block])
            else
                mask_rows = (cmin[block] + 1):(cmax[block] + 1)
                mask_cols = (rmin[block] + 1):(rmax[block] + 1)
            end

            mask_img[mask_rows, mask_cols] = NaN
        end
    end
    close(fpm_fits)
end

"""
Read a catalog entry.

Args:
 - field_dir: The directory of the file
 - run_num: The run number
 - camcol_num: The camcol number
 - field_num: The field number
 - bandnum: The band number from which to read galaxy properties.  Defaults
    to 3, the r band

Returns:
 - A data frame containing catalog entries in the rows.

 The data format is documented here:
 http://data.sdss3.org/datamodel/files/BOSS_PHOTOOBJ/RERUN/RUN/CAMCOL/photoObj.html

 This is based on the function get_sources in tractor/sdss.py:
 https://github.com/dstndstn/tractor/
"""
function load_catalog_df(field_dir, run_num, camcol_num, field_num; bandnum=3)

    cat_filename = "$field_dir/photoObj-$run_num-$camcol_num-$field_num.fits"
    cat_fits = FITSIO.FITS(cat_filename)

    cat_hdu = cat_fits[2]

    # Eliminate "bright" objects.
    # photo_flags1_map is defined in astrometry.net/sdss/common.py
    # ...where we see that photo_flags1_map['BRIGHT'] = 2
    const bright_bitmask = 2
    is_bright = read(cat_hdu, "objc_flags") & bright_bitmask .> 0
    has_child = read(cat_hdu, "nchild") .> 0

    # This is the position.  In tractor, it is passed into tractor:RaDecPos(ra, dec)
    objid = read(cat_hdu, "objid");
    ra = read(cat_hdu, "ra")
    dec = read(cat_hdu, "dec")

    # 6 = star, 3 = galaxy, others can be ignored.
    objc_type = read(cat_hdu, "objc_type");
    is_star = objc_type .== 6
    is_gal = objc_type .== 3
    is_bad_obj = !(is_star | is_gal)

    # Read in the galaxy types.
    fracdev = read(cat_hdu, "fracdev")[bandnum, :][:];
    has_dev = fracdev .> 0.
    has_exp = fracdev .< 1.
    is_comp = has_dev & has_exp
    is_bad_fracdev = (fracdev .< 0.) | (fracdev .> 1)

    # Read the fluxes.
    psfflux = read(cat_hdu, "psfflux")

    # Record the cmodelflux if the galaxy is composite, otherwise use
    # the flux for the appropriate type.
    cmodelflux = read(cat_hdu, "cmodelflux")
    devflux = read(cat_hdu, "devflux")
    expflux = read(cat_hdu, "expflux")

    # Only collect the following properties for a particular band.
    # NB: the phi quantites in tractor are multiplied by -1.
    phi_dev_deg = read(cat_hdu, "phi_dev_deg")[bandnum, :][:]
    phi_exp_deg = read(cat_hdu, "phi_exp_deg")[bandnum, :][:]

    theta_dev = read(cat_hdu, "theta_dev")[bandnum, :][:]
    theta_exp = read(cat_hdu, "theta_exp")[bandnum, :][:]

    ab_exp = read(cat_hdu, "ab_exp")[bandnum, :][:]
    ab_dev = read(cat_hdu, "ab_dev")[bandnum, :][:]

    close(cat_fits)

    # Match the column names in the stamp-making script.
    cat_df = DataFrames.DataFrame(objid=objid, ra=ra, dec=dec,
                                  is_star=is_star, is_gal=is_gal, frac_dev=fracdev,
                                  ab_exp=ab_exp, theta_exp=theta_exp, phi_exp=phi_exp_deg,
                                  ab_dev=ab_dev, theta_dev=theta_dev, phi_dev=phi_dev_deg)

    cat_df[:run] = run_num
    cat_df[:camcol] = camcol_num
    cat_df[:field] = field_num

    for b=1:length(band_letters)
        band_letter = band_letters[b]
        cat_df[symbol(string("psfflux_", band_letter))] = psfflux[b, :][:]
        cat_df[symbol(string("compflux_", band_letter))] = cmodelflux[b, :][:]
        cat_df[symbol(string("expflux_", band_letter))] = expflux[b, :][:]
        cat_df[symbol(string("devflux_", band_letter))] = devflux[b, :][:]
    end

    is_bad = is_bad_fracdev | is_bad_obj | is_bright | has_child
    bad_frac = sum(is_bad) / length(objid)
    println("Proportion of bad rows: $bad_frac")
    cat_df = cat_df[!is_bad, :]

    cat_df
end

end


text:
using Random

const Series = AbstractVector{Int}

abstract type InfoDist end

mutable struct MIDist <: InfoDist
    joint::Matrix{Int}
    m1::Vector{Int}
    m2::Vector{Int}
    N::Int
    MIDist() = new(zeros(Int, 2, 2), zeros(Int, 2), zeros(Int, 2), 0)
end

MIDist(xs::Series, ys::Series) = accumulate!(MIDist(), xs, ys)

function entropy(dist::MIDist)
    mi = 0.0
    for i in 1:2, j in 1:2
        p = dist.joint[i,j]
        if !iszero(p)
            mi += p * log2(p / (dist.m1[i] * dist.m2[j]))
        end
    end
    log2(dist.N) + mi/dist.N
end

function accumulate!(dist::MIDist, xs::Series, ys::Series)
    dist.N += length(xs)
    @inbounds for i in eachindex(xs)
        x, y = xs[i], ys[i]
        dist.m1[x] += 1
        dist.m2[y] += 1
        dist.joint[x, y] += 1
    end
    dist
end

@inline function clear!(dist::MIDist)
    dist.joint[:] .= 0
    dist.m1[:] .= 0
    dist.m2[:] .= 0
    dist.N = 0
    dist
end

function mutualinfo!(dist::MIDist, xs::Series, ys::Series)
    entropy(accumulate!(dist, xs, ys))
end

mutualinfo(xs::Series, ys::Series) = mutualinfo!(MIDist(), xs, ys)

mutable struct TEDist <: InfoDist
    k::Int
    states::Vector{Int}
    histories::Vector{Int}
    sources::Vector{Int}
    predicates::Vector{Int}
    N::Int

    function TEDist(k::Int)
        q = 2^k
        states = zeros(Int, 4q)
        histories = zeros(Int, q)
        sources = zeros(Int, 2q)
        predicates = zeros(Int, 2q)
        new(k, states, histories, sources, predicates, 0)
    end
end

@inline function clear!(dist::TEDist)
    dist.states[:] .= 0
    dist.histories[:] .= 0
    dist.sources[:] .= 0
    dist.predicates[:] .= 0
    dist.N = 0
    dist
end

function accumulate!(dist::TEDist, xs::Series, ys::Series)
    rng = dist.k:(length(ys)-1)
    dist.N = length(rng)
    history, q = 0, 1
    @inbounds for i in 1:dist.k
        q *= 2
        history = 2history + ys[i] - 1;
    end
    @inbounds for i in rng
        src = xs[i] - 1
        future = ys[i + 1] - 1
        source = 2history + src
        predicate = 2history + future
        state = 2predicate + src

        dist.states[state + 1] += 1
        dist.histories[history + 1] += 1
        dist.sources[source + 1] += 1
        dist.predicates[predicate + 1] += 1

        history = predicate - q*(ys[i - dist.k + 1] - 1)
    end
    dist
end

function entropy(xs::AbstractArray{Int}, N::Int)
    h = N * log2(N)
    for i in eachindex(xs)
        n = xs[i]
        if !iszero(n)
            h -= n * log2(n)
        end
    end
    h / N
end

function entropy(dist::TEDist)
    entropy(dist.sources, dist.N) +
    entropy(dist.predicates, dist.N) -
    entropy(dist.states, dist.N) -
    entropy(dist.histories, dist.N)
end

function transferentropy!(dist::TEDist, xs::Series, ys::Series)
    entropy(accumulate!(dist, xs, ys))
end

function transferentropy(xs::Series, ys::Series, k::Int)
    transferentropy!(TEDist(k), xs, ys)
end

function significance(rng::AbstractRNG, dist::InfoDist, measure::Function,
                      xs::Series, ys::Series;
                      nperms=1000, pvalue=0.05, spotcheck=1000)
    gt = measure(dist, xs, ys)
    clear!(dist)

    count = 0
    xsperm = xs[:]

    if nperms <= spotcheck
        for _ in 1:nperms
            count += (measure(dist, shuffle!(rng, xsperm), ys) ≥ gt)
            clear!(dist)
        end
        p = count / (nperms + 1)
        se = sqrt((p * (1 - p)) / (nperms + 1))

        return gt, p, se
    else
        for _ in 1:spotcheck
            count += (measure(dist, shuffle!(rng, xsperm), ys) ≥ gt)
            clear!(dist)
        end
        p = count / (spotcheck + 1)
        se = sqrt((p * (1 - p)) / (spotcheck + 1))

        if abs(p - pvalue) < 2pvalue
            for _ in 1:(nperms - spotcheck)
                count += (measure(dist, shuffle!(rng, xsperm), ys) ≥ gt)
                clear!(dist)
            end
            p = count / (nperms + 1)
            se = sqrt((p * (1 - p)) / (nperms + 1))
        end

        return gt, p, se
    end
end

function issig(datum::NTuple{3, Float64}; p::Float64=0.05)
    if p < zero(p) || p > one(p)
        throw(DomainError(p, "must be between 0.0 and 1.0"))
    end
    datum[2] < p
end

function signify(data::AbstractArray{NTuple{3,Float64}}; p::Float64=0.05)
    sigdata = first.(data)
    sigdata[.!issig.(data; p=p)] .= 0.0
    sigdata
end

max_stars_repo_name: ELIFE-ASU/dglmoore

max_issues_repo_name: ELIFE-ASU/dglmoore

max_forks_repo_name: ELIFE-ASU/dglmoore


text:
module MCMCSamples

using ArgCheck
using Lazy

export
    MCMCDraws,
    MCMCChain,
    pool

const MCMCDrawsVector = AbstractVector

const MCMCDrawsDict = Dict{Symbol, MCMCDrawsVector}

"""
Structure for storing MCMC draws.

# Interface

Draws are `AbstractVector`s, indexed by `Symbol`s. Uniform length is
enforced. `keys` are ordered, allowing determnistic traversal, but
this is not guaranteed in the internals.

# Usage

Use this structure as a collection of vectors which are ordered so
that corresponding elements match, but the contents are not
necessarily contiguous or come from the same chain — ie a pooled
sample. Useful for posterior inference and predictive checks.

For convergence analysis, use `MCMCChains`.
"""
struct MCMCDraws
    dict::MCMCDrawsDict
    function MCMCDraws(dict::MCMCDrawsDict)
        @argcheck !isempty(dict) "Need at least one vector of draws."
        N = length(first(values(dict)))
        @argcheck all(length(x) == N for x in values(dict)) "Inconsistent vector length."
        new(dict)
    end
end

MCMCDraws(pairs::Pair{Symbol, <:AbstractArray}...) = MCMCDraws(MCMCDrawsDict(pairs))

MCMCDraws(itr) = MCMCDraws(MCMCDrawsDict(itr))

Base.length(draws::MCMCDraws) = length(first(values(draws.dict)))

Base.keys(draws::MCMCDraws) = sort!(collect(keys(draws.dict)))

Base.getindex(draws::MCMCDraws, var::Symbol) = draws.dict[var]

"""
Test if all keys of `itr` are the same.
"""
function all_same_keys(itr)
    @argcheck !isempty(itr)
    keys1 = keys(first(itr))
    all(keys(draws) == keys1 for draws in drop(itr, 1))
end

const MCMCIndex = AbstractVector{Int}

"Test if `index` is valid for MCMC (nonempty, positive, increasing)."
function is_valid_index(index::MCMCIndex)
    !isempty(index) && first(index) > 0 && issorted(index)
end

const MCMCChainMeta = Dict{Symbol, Any}

"""
Markov Chain Monte Carlo draws with additional information.

# Fields

- `index` is an `AbstractVector` of integers, for indexing
  draws. Usually a `UnitRange` like `1000:2000`, or a `StepRange`
  `5000:10:10000` for a thinned sample.

- first `discard` elements of `index` are not meant to be used for
  posterior inference (because they are part of the adaptation, or
  burn-in, etc). These are discarded when pooled.

- `draws`: vectors and variable names, see `MCMCDraws`.

- `meta`: a dictionary indexed by symbols, for other metadata, related
  to sampling (eg mass matrix, adaptation information, information for
  useful continuation of sampling, etc).

# Usage

Use this structure (or a collection of them) for convergence analysis.

Pool a collection of these structures for posterior inference (and
predictive checks, etc).
"""
struct MCMCChain
    index::MCMCIndex
    discard::Int
    draws::MCMCDraws
    meta::MCMCChainMeta
    function MCMCChain(index, discard, draws, meta)
        @argcheck is_valid_index(index) "Invalid index."
        @argcheck length(index) == length(draws) "Inconsistent lengths."
        @argcheck 0 ≤ discard ≤ length(index)
        new(index, discard, draws, meta)
    end
end

@forward MCMCChain.draws Base.keys, Base.length, Base.getindex

# this is the preferred interface
MCMCChain(draws::MCMCDraws; discard = 0, index = (1:length(draws))+discard,
          meta = MCMCChainMeta()) = MCMCChain(index, discard, draws, meta)

"""
Pool MCMC draws into a single `MCMCDraws` structure.
"""
function pool(itr)
    @argcheck all_same_keys(itr) "Incompatible variable names."
    MCMCDraws(key => vcat((@view(c[key][(c.discard+1):end]) for c in itr)...)
              for key in keys(first(itr)))
end

pool(chains::MCMCChain...) = pool(chains)

end # module


text:
using Random, Distributions, Optim, Plots, Statistics

abstract type TypeData end

struct ObservedData
    x::Array{Float64}
    w::Array{Float64}
    y::Array{Float64}
end

struct ExperimentTuner
    T::Int
    n::Int
    ξ::Float64
    α::Array{Float64}
    β₀::Array{Float64}
end

mutable struct IterativeUpdater
    update::Function
    β::Array{Float64, 2}
    π::Array{Float64, 1}
end

function IterativeUpdater(update::Function, tuner::ExperimentTuner)
    k = length(tuner.β₀)
    β = zeros(Float64, (tuner.T, k))
    π = zeros(Float64, tuner.T)
    IterativeUpdater(update, β, π)
end

function IterativeUpdater(βfixed::Array{Float64}, tuner::ExperimentTuner)
    k = length(tuner.β₀)
    β = zeros(Float64, (tuner.T, k))
    π = zeros(Float64, tuner.T)
    function fixedβ(a, types::TypeData, t, tuner::ExperimentTuner)
        return βfixed
    end
    IterativeUpdater(fixedβ, β, π)
end

function broadcastArray(β::Array{Float64}, n)
    if ndims(β)== 1
        β  = zeros(n, 2) .+ transpose(β)
    end
    return β
end

function runExperiment(tuner::ExperimentTuner, methods::Array{IterativeUpdater},
                       typeDGP::Type, objective::Function)
    types = typeDGP(tuner.n)
    for m in methods
        m.β[1, :] = m.update(tuner.β₀, types, 1, tuner)
        m.π[1] = mean(objective(m.β[1, :], types))
    end
    for t in 2:tuner.T
        types = typeDGP(tuner.n)
        for m in methods
            m.β[t, :] = m.update(m.β[t-1, :], types, t, tuner)
            m.π[t] = mean(objective(m.β[t, :], types))
        end
    end
end

function naive_solution(n::Int, typeDGP::Type, naiveObjective::Function)
    types = typeDGP(n)
    res = Optim.optimize(beta -> naiveObjective(beta, types.z, types), [1, 0.0])
    return Optim.minimizer(res)
end

function fk_solution(n::Int, typeDGP::Type, negObjective::Function)
    types = typeDGP(n)
    res = Optim.optimize(beta -> negObjective(beta, types), [1, 0.0])
    bstar = Optim.minimizer(res)
    return bstar
end

#utilities and tests
function linReg(y, x)
    x = [ones(length(y)) x ]
    β = inv(x'*x)*x'*y
    return β
end

#specific types
struct SCTypeData <: TypeData
    γ::Array{Float64}
    z::Array{Float64}
    r::Array{Float64}
end

function SCTypeData(n=1000)
    z = rand(Normal(0, 1), n)
    r = rand(Normal(0, 1), n)
    #used to be 1.5
    γ = rand(Uniform(0, 3), n)
    return SCTypeData(γ, z, r)
end

function ObservedData(types::SCTypeData, β::Array{Float64})
    β = broadcastArray(β, length(types.γ))
    x = types.z .+ types.γ.*(β[:, 2])
    y = types.z .+ types.r
    w = β[:, 1] .+ x.*β[:, 2]
    return ObservedData(x, w, y)
end

function negloss(β::Array{Float64, 1}, types::SCTypeData)
    data = ObservedData(types, β)
    Π = sum((data.y - data.w).^2)
    return Π
end

function loss(β::Array{Float64}, types::SCTypeData)
    data = ObservedData(types, β)
    π = (data.y - data.w).^2
    return -π
end

function SCUpdate(β, types::SCTypeData, t, tuner::ExperimentTuner)
    step = tuner.α[1]/(2)
    ϵ = [ zeros(tuner.n) tuner.ξ.*rand([-1,1], (tuner.n, length(β)-1))]
    βp = β' .+ ϵ
    data = ObservedData(types, βp)
    π = loss(βp, types)

    γ̂ = linReg(π, ϵ[:, 2])[2]
    β_new = zeros(2)
    β_new[2] = β[2] + step*γ̂
    β_new[1] = mean(data.y)  - mean(data.x.*β_new[2])
    return β_new
end
function robustUpdate(β, types::SCTypeData, t, tuner::ExperimentTuner)
    step = tuner.α/(2)
    ϵ = tuner.ξ.*rand([-1,1], (tuner.n, length(β)))
    βp = β' .+ ϵ
    data = ObservedData(types, βp)
    π = loss(βp, types)
    γ̂ = linReg(π, ϵ)[2:(length(β)+1)]
    return β .+ step.*γ̂
end

function fixedXObjective(β, x::Vector, types::SCTypeData)
    w = β[1] .+ β[2] .* x
    y = types.z .+ types.r
    π = (y .- w).^2
    return sum(π)
end

function naiveUpdate(β::Array{Float64}, types::SCTypeData, t::Int, tuner::ExperimentTuner)
    data = ObservedData(types, β)
    return linReg(data.y, data.x)
end

function runSCExperiment(T)
    n = 5000
    β_fk = fk_solution(100000, SCTypeData, negloss)

    β_naive = naive_solution(100000, SCTypeData, fixedXObjective)
    #tuner = ExperimentTuner(T, n, 0.035, [3e-3], [-0.65, 0.75])
    tuner = ExperimentTuner(T, n, 0.075, [3e-3], [-0.65, 0.75])
    methods = [IterativeUpdater(β_fk, tuner),
               IterativeUpdater(SCUpdate, tuner),
               IterativeUpdater(naiveUpdate, tuner),
               IterativeUpdater(β_naive, tuner)]
    runExperiment(tuner, methods, SCTypeData, loss)
    βs = []
    fk_profits = mean(methods[1].π)

    for m in methods
        push!(βs, mean(m.β[:, 2], dims=2))
        #print regret
        println(mean(m.π))
        println(mean(m.π) - fk_profits)
        println(m.β[T, :])
    end
    plot( βs, xlabel = "t", ylabel="Prediction Weight",
    label=["Full Information" "Learning via Experiment" "Repeated Risk Min" "Naive Risk Min"], ylim = [0.5, 1.2])
    savefig("figures/sc_sim.pdf")
end

max_stars_repo_name: evanmunro/personalized-policy

max_issues_repo_name: evanmunro/personalized-policy

max_forks_repo_name: evanmunro/personalized-policy


text:
function update_perceiv_publ_opinion!(
    state::Tuple{AbstractGraph, AbstractArray}, agent_idx::Integer
)
    graph, agent_list = state
    this_agent = agent_list[agent_idx]
    # get neighborhood opinion as baseline
    input = inneighbors(graph, agent_idx)
    if length(input) != 0
        input_opinion_mean = mean([agent_list[input_agent].opinion for input_agent in input])
    else
        input_opinion_mean = this_agent.opinion
    end
    # compute feed opinion
    feed_opinions = [post.opinion for post in this_agent.feed]
    feed_weights = [post.weight for post in this_agent.feed]
    if length(feed_opinions) > 0
        feed_opinion_mean = (
            sum([opinion * weight for (opinion, weight) in zip(feed_opinions, feed_weights)]) /
            sum(feed_weights)
        )
    else
        feed_opinion_mean = this_agent.opinion
    end
    # perceived public opinion is the mean between the feed and neighborhood opinion
    this_agent.perceiv_publ_opinion = mean([input_opinion_mean, feed_opinion_mean])
    return state
end

# suppress output of include()
;


text:
struct Model{F1<:Function,F2<:Function,F3<:Function}
  μ!::F1
  G::F2
  n::Int64
  simulate::F3
  Tx::Type
  Ty::Type
end

function simulateModel(model::Model, n::Int64)
  ys = Vector{model.Ty}(undef, n)
  rng = getRNG()
  x::model.Tx = model.Tx()
  for p in 1:n
    model.μ!(x, rng)
    ys[p] = model.simulate(x, rng)
  end
  return ys
end


text:
"""
    mcse(x::AbstractVector{<:Real}; method::Symbol=:imse, kwargs...)

Compute the Monte Carlo standard error (MCSE) of samples `x`.
The optional argument `method` describes how the errors are estimated. Possible options are:

- `:bm` for batch means [^Glynn1991]
- `:imse` initial monotone sequence estimator [^Geyer1992]
- `:ipse` initial positive sequence estimator [^Geyer1992]

[^Glynn1991]: Glynn, P. W., & Whitt, W. (1991). Estimating the asymptotic variance with batch means. Operations Research Letters, 10(8), 431-435.

[^Geyer1992]: Geyer, C. J. (1992). Practical Markov Chain Monte Carlo. Statistical Science, 473-483.
"""
function mcse(x::AbstractVector{<:Real}; method::Symbol=:imse, kwargs...)
    return if method === :bm
        mcse_bm(x; kwargs...)
    elseif method === :imse
        mcse_imse(x)
    elseif method === :ipse
        mcse_ipse(x)
    else
        throw(ArgumentError("unsupported MCSE method $method"))
    end
end

function mcse_bm(x::AbstractVector{<:Real}; size::Int=floor(Int, sqrt(length(x))))
    n = length(x)
    m = min(div(n, 2), size)
    m == size || @warn "batch size was reduced to $m"
    mcse = StatsBase.sem(Statistics.mean(@view(x[(i + 1):(i + m)])) for i in 0:m:(n - m))
    return mcse
end

function mcse_imse(x::AbstractVector{<:Real})
    n = length(x)
    lags = [0, 1]
    ghat = StatsBase.autocov(x, lags)
    Ghat = sum(ghat)
    @inbounds value = Ghat + ghat[2]
    @inbounds for i in 2:2:(n - 2)
        lags[1] = i
        lags[2] = i + 1
        StatsBase.autocov!(ghat, x, lags)
        Ghat = min(Ghat, sum(ghat))
        Ghat > 0 || break
        value += 2 * Ghat
    end

    mcse = sqrt(value / n)

    return mcse
end

function mcse_ipse(x::AbstractVector{<:Real})
    n = length(x)
    lags = [0, 1]
    ghat = StatsBase.autocov(x, lags)
    @inbounds value = ghat[1] + 2 * ghat[2]
    @inbounds for i in 2:2:(n - 2)
        lags[1] = i
        lags[2] = i + 1
        StatsBase.autocov!(ghat, x, lags)
        Ghat = sum(ghat)
        Ghat > 0 || break
        value += 2 * Ghat
    end

    mcse = sqrt(value / n)

    return mcse
end


text:
using DynamicPolynomials
using MultivariatePolynomials
using Test
using LinearAlgebra

# TODO move to MP
@testset "Issue #70" begin
    @ncpolyvar y0 y1 x0 x1
    p = x1 * x0 * x1
    @test subs(p, x0 => y0, x1 => y1) == y1 * y0 * y1
    @test subs(p, x0 => 1) == x1^2
    @test p(x0 => y0, x1 => y1) == y1 * y0 * y1
end

@testset "Issue #79, Issue #80 and Issue #92" begin
    @polyvar x[1:2]
    p1 = x[1] * 0.0 + x[2] * 0
    p2 = ( x[1] + x[2] ) * 0.0
    @test variables(p1) == x
    @test variables(p1) == variables(p2)

    @polyvar χ[1:4]
    poly_array_int = [1 0 1 0; 0 1 0 1]*χ;
    poly_array_float = Float64.([1 0 1 0; 0 1 0 1])*χ;

    vars_array_int = variables.(poly_array_int);
    vars_array_float = variables.(poly_array_float)
    @test isa(vars_array_int, Vector{<:Vector{<:PolyVar}})
    @test isa(vars_array_float, Vector{<:Vector{<:PolyVar}})
    @test vars_array_int == vars_array_float
    @test vars_array_int[1] == vars_array_float[2] == χ

    # This is for Issue #92
    p3 = subs(p1, x => x)
    @test variables(p3) == x

    p4 = x[1] + x[2]
    @test variables(subs(p4, x[2] => 1)) == [x[1]]
    @test variables(subs(p4, x[1] => 1)) == [x[2]]
end

include("mono.jl")
include("poly.jl")
include("comp.jl")
include("mutable_arithmetics.jl")

module newmodule
    using Test
    import DynamicPolynomials
    @testset "Polyvar macro hygiene" begin
        # Verify that the @polyvar macro works when the package has been activated
        # with `import` instead of `using`.
        DynamicPolynomials.@polyvar x y
        @test isa(x, DynamicPolynomials.PolyVar)
        @test isa(y, DynamicPolynomials.PolyVar)
    end
end

include("mvp.jl")

max_stars_repo_name: blegat/DynamicPolynomials

max_issues_repo_name: blegat/DynamicPolynomials

max_forks_repo_name: blegat/DynamicPolynomials


text:
using DataStructures
using DataFrames
using CSV
using Random

# sample based on count,
# i.e. only keep values in (lower bound, upper bound) from sample
function key_params(df, lower_bound, upper_bound)
    df[lower_bound .< df.count .< upper_bound, :]
end

function post_date_right_open_range_params(sample, lower_bound, upper_bound)
    results = []
    for ix in 1:length(sample)
        start_offset = sample[ix][1]
        count_sum = 0
        for v in sample[ix:length(sample)]
            count_sum += v[2]
        end
        if lower_bound < count_sum < upper_bound
            push!(results, [start_offset, count_sum])
        end
    end
    return results
end

function post_month_params(sample, lower_bound, upper_bound)
    results = []
    for ix in 0:length(sample)÷4 - 1
        start_ix = ix*4 + 1
        end_ix = min(ix*4 + 4, length(sample))
        count_sum = 0
        for v in sample[start_ix:end_ix]
            count_sum += v[2]
        end
        start_day = sample[start_ix][1]
        end_day = sample[min(end_ix + 1, length(sample))][1] # end_ix + 1 to get the first day of the next month
        if lower_bound < count_sum < upper_bound
            push!(results, [[start_day, end_day], count_sum])
        end
    end
    return results
end

function add_months(sourceYear, sourceMonth, sourceDay, nMonths)
    newMonth = sourceMonth - 1 + nMonths
    newYear = sourceYear + newMonth ÷ 12
    newMonth = newMonth % 12 + 1
    newDay = 1
    return (newYear, newMonth, newDay)
end

function convert_posts_histo(histogram, timestamps) # TODO: use timestamps
  week_posts = []
  for month in 1:length(histogram)
    # split total into 4 weeks
    monthTotal = histogram[month]
    (year, month, day) = add_months(2010, 01, 01, month-1)
    push!(week_posts, [(year, month, day   ) , monthTotal ÷ 4])
    push!(week_posts, [(year, month, day+7 ) , monthTotal ÷ 4])
    push!(week_posts, [(year, month, day+14) , monthTotal ÷ 4])
    push!(week_posts, [(year, month, day+21) , monthTotal ÷ 4])
  end

  return week_posts
end

# if length(ARGS) < 2
#     println("arguments: <input dir> <output dir>")
#     exit(1)
# end

# indir = ARGS[1] * "/"
# outdir = ARGS[2] * "/"

cd("../paramgenerator")
indir = "../out/build/"
outdir = "../substitution_out/"

###############################
### load CSVs #################
###############################

files = readdir(indir)

# activity factors
postsPerCountryFactorFiles = filter(f -> endswith(f, "postsPerCountryFactors.csv"), files)
tagClassCountFactorFiles = filter(f -> endswith(f, "tagClassCountFactors.csv"), files)
tagCountFactorFiles = filter(f -> endswith(f, "tagCountFactors.csv"), files)
firstNameCountFactorFiles = filter(f -> endswith(f, "firstNameCountFactors.csv"), files)
miscFactorFiles = filter(f -> endswith(f, "miscFactors.csv"), files)

# person factors
personFactorFiles = filter(x -> occursin(r"personFactors\.csv$", x), files)

# friends
friendsFiles = filter(x -> occursin(r"friendList\d+\.csv$", x), files)

# factors

## activityFactors
countFriends = CSV.read(indir * friendsFiles[1], DataFrame; delim='|', header=["person", "friendCount"])

country_sample = CSV.read(indir * postsPerCountryFactorFiles[1], DataFrame; delim='|', header=["country", "count"])
tagclass_posts = CSV.read(indir * tagClassCountFactorFiles[1], DataFrame; delim='|', header=["tagClass", "count"])
tag_posts = CSV.read(indir * tagCountFactorFiles[1], DataFrame; delim='|', header=["tag", "count"])
# unused in the BI workload
#nameFactors = CSV.read(indir * firstNameCountFactorFiles[1], DataFrame; delim='|', header=["name", "count"])
timestamps = CSV.read(indir * miscFactorFiles[1], DataFrame; delim='|')

## person, friend, and foaf factors
personFactorFile = personFactorFiles[1]
friendsFile = friendsFiles[1]

friends = CSV.read(indir * friendsFile, DataFrame; delim='|', header=["person", "friend"])
personFactors = CSV.read(indir * personFactorFile, DataFrame; delim='|', header=["person", "name", "f", "p", "pl", "pt", "g", "w", "pr", "numMessages", "numForums"])

###############################
### preprocess data frames ####
###############################

# create arrays from nested semicolon-separated fields
personFactors[!, :numMessages] = map.(x -> parse(Int64, x), split.(personFactors[!, :numMessages], ";"))
personFactors[!, :numForums] = map.(x -> parse(Int64, x), split.(personFactors[!, :numForums], ";"))
personFactorsAggregated = combine(
    groupby(personFactors, []),
    :numMessages => sum => :numMessages,
    :numForums => sum => :numForums # unused
)

# determine factors for friends by computing the sums of the factors grouped by person for the join expression:
# friends join_{friends[friend] = personFactors[person]} personFactors
friendsFactors =
    combine(
        groupby(
            innerjoin(friends, personFactors, on = [:friend => :person])
            , :person
        ),
        :f => sum => :f,
        :p => sum => :p,
        :pl => sum => :pl,
        :pt => sum => :pt,
        :g => sum => :g,
        :w => sum => :w,
        :pr => sum => :pr
    )

# determine factors for friends of friends
foafFactors =
    combine(
        groupby(
            innerjoin(friends, friendsFactors, on = [:friend => :person])
            , :person
        ),
        :f => sum => :f,
        :p => sum => :p,
        :pl => sum => :pl,
        :pt => sum => :pt,
        :g => sum => :g,
        :w => sum => :w,
        :pr => sum => :pr
    )

# shuffling/ sorting and sampling data
persons = personFactors[:, :person]
shuffle!(MersenneTwister(1234), persons)
sort!(tag_posts, order(:count, rev=true))
sort!(tagclass_posts, order(:count, rev=true))

# aggregating data
total_posts = sum(tag_posts.count)
person_sum = sum(country_sample.count)

posts_histogram = personFactorsAggregated.numMessages[1]
week_posts = convert_posts_histo(posts_histogram, timestamps)
non_empty_weeks = length(filter(x -> x[2] != 0, week_posts))

# computing posts/month
post_lower_threshold = (total_posts÷(non_empty_weeks÷4))*0.8
post_upper_threshold = (total_posts÷(non_empty_weeks÷4))*1.2
post_months = post_month_params(week_posts, post_lower_threshold, post_upper_threshold)

# some constants
path_bounds = [2, 5] # TODO: path bound for BI Q16 // 10 need a revision anyways
language_codes = [["ar"], ["tk"], ["tk"], ["uz"], ["uz"], ["uz"], ["uz"], ["uz"], ["uz"], ["uz"], ["uz"], ["uz", "tk"], ["uz", "tk"]]
post_lengths = [20,40,113,97,240]

###############################
### compute params for BI #####
###############################

# Q1 // 1
post_date_right_open_range_params(week_posts, 0.3*total_posts, 0.6*total_posts)

# Q3 // 2
bi3 = post_months

# Q4 // 3
bi4 = []
for i in eachrow(key_params(tagclass_posts, total_posts÷20, total_posts÷10))
    for j in eachrow(key_params(country_sample, total_posts÷150, total_posts÷50))
        push!(bi4, [i[1], j[1]])
    end
end
bi4

# Q5 // 4
bi5 = key_params(country_sample, total_posts÷200, total_posts÷100)

# Q6 // 5
bi6 = key_params(tag_posts, total_posts÷1300, total_posts÷900)

# Q7 // 6
bi7 = key_params(tag_posts, total_posts÷900, total_posts÷600)

# Q8 // 7
bi8 = key_params(tag_posts, total_posts÷600, total_posts÷300)

# Q10 // 8
bi10 = []
for i in eachrow(key_params(tag_posts, total_posts÷900, total_posts÷600))
    for j in post_date_right_open_range_params(week_posts, 0.3*total_posts, 0.6*total_posts)
        push!(bi10, [i[1], j[1]])
    end
end
bi10

# Q14 // 9
post_months

# Q16 // 10
bi16 = []
r16 = MersenneTwister(1234)
for i in eachrow(key_params(tagclass_posts, total_posts÷30, total_posts÷10))
    for j in eachrow(key_params(country_sample, total_posts÷80, total_posts÷20))
        person = rand(r16, persons)
        push!(bi16, [i[1], j[1], person])
        # TODO: add bounds
    end
end
bi16

# Q17 // 11
bi17 = key_params(country_sample, total_posts÷200, total_posts÷100)

# Q18 // 12
bi18 = []
for week in post_date_right_open_range_params(week_posts, 0.3*total_posts, 0.6*total_posts)
    for length in post_lengths
        for langs in language_codes
            push!(bi18, [week, length, langs])
        end
    end
end
bi18

# Q21 // 13
bi21 = []
for c in eachrow(key_params(country_sample, total_posts÷200, total_posts÷100))
    push!(bi21, [c[:country], 1356994800000])
end
bi21

# Q22 // 14
bi22 = []
for i in eachrow(key_params(country_sample, total_posts÷120, total_posts÷40))
    for j in eachrow(key_params(country_sample, total_posts÷120, total_posts÷40))
        if i[:country] < j[:country]
            push!(bi22, [i[:country], j[:country]])
        end
    end
end
bi22

# Q25 // 15
bi25 = []
r25 = MersenneTwister(1234)
for day_range in post_months
    for i in 1:min(length(persons), 10)
        person1Id = rand(r25, persons)
        person2Id = rand(r25, filter(p -> p != person1Id, persons))
        push!(bi25, [person1Id, person2Id, day_range[1], day_range[2]])
    end
end
bi25


text:
import Base: mod, real, imag, reim, conj, promote_rule

export GaussMod, AbstractMod

# mod support for Gaussian integers until officially adopted into Base
mod(z::Complex{<:Integer}, n::Integer) = Complex(mod(real(z), n), mod(imag(z), n))

"""
`GaussMod{m}(v)` creates a modular number in mod `m` with value `mod(v,m)`.
`GaussMod{m}()` is equivalent to `GaussMod{m}(0)`.
"""
struct GaussMod{N} <: AbstractMod
    val::Complex{Int}
    function GaussMod(x::Complex{T},N::Int) where T<:Integer
        @assert N>1 "modulus must be at least 2"
        new{N::Int}(mod(x,N))
    end
end

GaussMod{N}(x::Complex{T}=0+0im) where {N,T<:Integer} = GaussMod(x,N)
GaussMod{N}(x::T) where {N,T<:Integer}  = GaussMod{N}(Complex{Int}(x))
GaussMod{N}(x::S,y::T) where {N,S<:Integer,T<:Integer} = GaussMod{N}(x + y*im)
GaussMod(x::Mod{N}) where N = GaussMod{N}(x.val)
GaussMod{N}(x::Mod{N}) where N = GaussMod{N}(x.val)

GaussMod{N}(x::Rational{T}) where {N,T<:Integer} = GaussMod{N}(numerator(x)) / GaussMod{N}(denominator(x))
GaussMod{N}(x::Complex{Rational{T}}) where {N,T<:Integer} = GaussMod{N}(real(x)) + GaussMod{N}(imag(x))*im

Mod{N}(a::S,b::T) where {N,S<:Integer,T<:Integer} = GaussMod{N}(a+b*im)

modulus(x::GaussMod{T}) where T = T

real(x::GaussMod{N}) where N = Mod{N}(real(x.val))
real(x::Mod) = x

imag(x::GaussMod{N}) where N = Mod{N}(imag(x.val))
imag(x::Mod{N}) where N = Mod{N}(0)

reim(x::AbstractMod) = (real(x),imag(x))

conj(x::GaussMod{N}) where N = GaussMod{N}(conj(x.val))

Mod{N}(x::Complex) where N = GaussMod{N}(x)

function Mod{N}(x::GaussMod{N}) where N
    if imag(x) == 0
        return Mod{N}(real(x))
    end
    error("Cannot convert $x to type Mod{$N} (nonzero imaginary part)\nPerhaps use: real(x)")
end

Mod(x::GaussMod{N}) where N = Mod{N}(x)

# ARITHMETIC

function(+)(x::GaussMod{N}, y::GaussMod{N}) where N
    xx = widen(x.val)
    yy = widen(y.val)
    zz = mod(xx+yy,N)
    return Mod{N}(zz)
end

(-)(x::GaussMod{N}) where N = GaussMod{N}(-x.val)
(-)(x::GaussMod{N},y::GaussMod{N}) where N = x + (-y)

function (*)(x::GaussMod{N}, y::GaussMod{N}) where N
    xx = widen(x.val)
    yy = widen(y.val)
    zz = mod(xx*yy,N)
    return Mod{N}(zz)
end

function inv(x::GaussMod{N}) where N
    try
        a,b = reim(x)
        bot = inv(a*a + b*b)
        aa = a*bot
        bb = -b*bot
        return aa + bb*im
    catch
        error("$x is not invertible")
    end
end

is_invertible(x::GaussMod) = is_invertible(real(x*x'))

(/)(x::GaussMod{N}, y::GaussMod{N}) where N = x * inv(y)

(//)(x::GaussMod{N}, y::GaussMod{N}) where N = x/y
(//)(x::Number, y::GaussMod{N}) where N = x/y
(//)(x::GaussMod{N}, y::Number) where N = x/y

show(io::IO, z::GaussMod{N}) where N = print(io,"Mod{$N}($(z.val))")

# Random

rand(::Type{GaussMod{N}}) where N = GaussMod{N}(rand(Int) + im*rand(Int))
rand(::Type{GaussMod{N}},dims::Integer...) where N = GaussMod{N}.(rand(Int,dims...)) + im*GaussMod{N}.(rand(Int,dims...))

promote_rule(::Type{GaussMod{N}}, ::Type{Mod{N}}) where N = GaussMod{N}
promote_rule(::Type{Mod{N}}, ::Type{GaussMod{N}}) where N = GaussMod{N}

promote_rule(::Type{GaussMod{N}}, ::Type{T}) where {N, T<:Integer} = GaussMod{N}
promote_rule(::Type{T}, ::Type{GaussMod{N}}) where {N, T<:Integer} = GaussMod{N}

promote_rule(::Type{GaussMod{N}}, ::Type{Complex{T}}) where {N, T<:Integer} = GaussMod{N}
promote_rule(::Complex{Type{T}}, ::Type{GaussMod{N}}) where {N, T<:Integer} = GaussMod{N}

promote_rule(::Type{Mod{N}}, ::Type{Complex{T}}) where {N, T<:Integer} = GaussMod{N}
promote_rule(::Complex{Type{T}}, ::Type{Mod{N}}) where {N, T<:Integer} = GaussMod{N}

promote_rule(::Type{GaussMod{N}}, ::Type{T}) where {N, T<:Rational} = GaussMod{N}
promote_rule(::Type{T}, ::Type{GaussMod{N}}) where {N, T<:Rational} = GaussMod{N}

promote_rule(::Type{GaussMod{N}}, ::Type{T}) where {N, T<:Complex} = GaussMod{N}
promote_rule(::Type{T}, ::Type{GaussMod{N}}) where {N, T<:Complex} = GaussMod{N}

promote_rule(::Type{Mod{N}}, ::Type{T}) where {N, T<:Complex} = GaussMod{N}
promote_rule(::Type{T}, ::Type{Mod{N}}) where {N, T<:Complex} = GaussMod{N}


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

function fuse_loops(block::Vector{Any})
    first = block[1]
    new_block = Any[first]
    if isa(first, Expr)
        if first.head == :for
            first.args[2].args = fuse_loops(first.args[2].args)
        elseif first.head == :block
            first.args = fuse_loops(first.args)
        end
    end
    for statement in block[2:end]
        if isa(statement, Expr) && statement.head == :for
            # Because we may have update_param calls in between fusable loops
            # we first remove those
            # updates = []
            # is_update_param_call = (node) -> isa(node, Expr) && node.head == :call && node.args[1] == :update_param
            # while length(new_block) > 0 && is_update_param_call(new_block[end])
            #     push!(updates, pop!(new_block))
            # end
            if isa(new_block[end], Expr) && new_block[end].head == :for
                curr_loop = statement
                prev_loop = new_block[end]
                curr_len = curr_loop.args[1].args[2]
                prev_len = prev_loop.args[1].args[2]
                if curr_len == prev_len
                    prev_loopvar = get_loopvar(prev_loop)
                    curr_loopvar = get_loopvar(curr_loop)
                    curr_body = AstWalk(curr_loop.args[2], symbol_replacer, Dict(curr_loopvar => prev_loopvar)).args
                    append!(prev_loop.args[2].args, curr_body)
                    prev_loop.args[2].args = fuse_loops(prev_loop.args[2].args)
                    # pop!(new_block)
                    # append!(new_block, updates)
                    # push!(new_block, prev_loop)
                elseif is_tiled_loop(prev_loop) && is_tiled_loop(curr_loop) && !LATTE_DISABLE_TILE_FUSION
                    forward_factor = get_tile_fusion_factor_forward(curr_len)
                    backward_factor = get_tile_fusion_factor_backward(prev_len)
                    prev_end = prev_len.args[2].args[2]
                    curr_end = curr_len.args[2].args[2]
                    if curr_end * forward_factor == prev_end
                        prev_len.args[2].args[2] = div(prev_end, forward_factor)
                        if forward_factor > 1
                            unshift!(prev_loop.args[2].args, :(TILE_SIZE *= $forward_factor))
                            push!(prev_loop.args[2].args, :(TILE_SIZE /= $forward_factor))
                        end
                        append!(prev_loop.args[2].args, curr_loop.args[2].args)
                        prev_loop.args[2].args = fuse_loops(prev_loop.args[2].args)
                        # pop!(new_block)
                        # append!(new_block, updates)
                        # push!(new_block, prev_loop)
                    elseif curr_end == prev_end * backward_factor
                        if backward_factor > 1
                            unshift!(curr_loop.args[2].args, :(TILE_SIZE *= $backward_factor))
                            push!(curr_loop.args[2].args, :(TILE_SIZE /= $backward_factor))
                        end
                        prev_len.args[4] *= get_tile_fusion_factor_backward(curr_len)
                        append!(prev_loop.args[2].args, curr_loop.args[2].args)
                        prev_loop.args[2].args = fuse_loops(prev_loop.args[2].args)
                        # pop!(new_block)
                        # append!(new_block, updates)
                        # push!(new_block, prev_loop)
                    else
                        # append!(new_block, updates)
                        statement.args[2].args = fuse_loops(statement.args[2].args)
                        push!(new_block, statement)
                    end
                else
                    # append!(new_block, updates)
                    statement.args[2].args = fuse_loops(statement.args[2].args)
                    push!(new_block, statement)
                end
            else
                # append!(new_block, updates)
                statement.args[2].args = fuse_loops(statement.args[2].args)
                push!(new_block, statement)
            end
        else
            push!(new_block, statement)
        end
    end
    filter((x) -> !(isa(x, Symbol) && x == :NOFUSE), new_block)
end


text:
environmentsEndpoint = "api/$curApiVersion/environments"
environmentEndpoint = "$environmentsEndpoint/{1}"

"""
$(SIGNATURES)
Gets all environments, if you specify details=true, you also receive associated sessions.
"""
function getEnvironments(;details=false)::Vector{EnvironmentResponse}
    config = getGraffConfig()
    if config == nothing
        error("Graff config is not set, please call setGraffConfig with a valid configuration.")
    end
    url = "$(config.apiEndpoint)/$(format(environmentsEndpoint))?details=$(details)"
    response = @mock _sendRestRequest(config, HTTP.get, url)
    if(response.status != 200)
        error("Error getting environments, received $(response.status) with body '$(String(response.body))'.")
    end
    environments = JSON2.read(String(response.body), Vector{EnvironmentResponse})
    # #Sort
    sort!(environments, by=(s -> s.id))
    return environments
end

"""
$(SIGNATURES)
Get a specific environment.
"""
function getEnvironment(envId::String)::EnvironmentResponse
    config = getGraffConfig()
    if config == nothing
        error("Graff config is not set, please call setGraffConfig with a valid configuration.")
    end
    url = "$(config.apiEndpoint)/$(format(environmentEndpoint, envId))"
    response = @mock _sendRestRequest(config, HTTP.get, url)
    if(response.status != 200)
        error("Error getting environment, received $(response.status) with body '$(String(response.body))'.")
    end
    return JSON2.read(String(response.body), EnvironmentResponse)
end

"""
$(SIGNATURES)
Delete a specific environment.
Return: Nothing if success, error if failed.
"""
function deleteEnvironment(envId::String)::Nothing
    config = getGraffConfig()
    if config == nothing
        error("Graff config is not set, please call setGraffConfig with a valid configuration.")
    end
    url = "$(config.apiEndpoint)/$(format(environmentEndpoint, envId))"
    response = @mock _sendRestRequest(config, HTTP.delete, url)
    if(response.status != 200)
        error("Error deleting environment, received $(response.status) with body '$(String(response.body))'.")
    end
    return nothing
end

"""
$(SIGNATURES)
Create an environment in Graff.
Return: Returns the created session.
"""
function addEnvironment(environment::EnvironmentRequest)::EnvironmentResponse
    config = getGraffConfig()
    if config == nothing
        error("Graff config is not set, please call setGraffConfig with a valid configuration.")
    end
    url = "$(config.apiEndpoint)/$(format(environmentEndpoint, environment.id))"
    body = JSON.json(environment)
    response = @mock _sendRestRequest(config, HTTP.post, url, data=body, headers=Dict{String, String}("Content-Type" => "application/json"))
    if(response.status != 200)
        error("Error creating session, received $(response.status) with body '$(String(response.body))'.")
    end
    body = String(response.body)
    return JSON2.read(body, SessionDetailsResponse)
end


text:
using FPGAUtils
using Documenter

DocMeta.setdocmeta!(FPGAUtils, :DocTestSetup, :(using FPGAUtils); recursive=true)

makedocs(;
    modules=[FPGAUtils],
    authors="Ian McInerney <ian.s.mcinerney@ieee.org> and contributors",
    repo="https://github.com/JuliaEmbedded/FPGAUtils.jl/blob/{commit}{path}#{line}",
    sitename="FPGAUtils.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaEmbedded.github.io/FPGAUtils.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaEmbedded/FPGAUtils.jl",
)


text:
module TestECDF

using QuantEcon
using Base.Test
using FactCheck

# set up
srand(42)
obs = rand(40)
e = ECDF(obs)

facts("Testing ecdf.jl") do
    # 1.1 is larger than all obs, so ecdf should be 1
    @fact ecdf(e, 1.1) --> roughly(1.0)

    # -1.0 is small than all obs, so ecdf should be 0
    @fact ecdf(e, -1.0) --> roughly(0.0)

    # larger values should have larger values on ecdf
    let x = rand()
        F_1 = ecdf(e, x)
        F_2 = ecdf(e, x*1.1)
        @fact F_1 <= F_2 --> true
    end

end  # facts
end  # module


text:
# This file roughly corresponds to functions documented in the
# Courses API: https://canvas.instructure.com/doc/api/courses

"""
    Canvas.courses(; kwargs...) -> Vector{Course}, page_data

Return a paginated list of courses for the authenticated user.
Return a vector of [`Course`](@ref)s and a dictionary with page data.

**Canvas API documentation**:
[List your courses (`GET /api/v1/courses`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.index)
"""
function courses(; kwargs...)
    json, page_data = paged_request("GET", "/api/v1/courses"; kwargs...)
    return Course.(json), page_data
end

"""
    Canvas.courses(user::User; kwargs...) -> Vector{Course}, page_data

Return a paginated list of courses for the specified user.
Return a vector of [`Course`](@ref)s and a dictionary with page data.

**Canvas API documentation**:
[List courses for a user (`GET /api/v1/users/:user_id/courses`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.user_index)
"""
function courses(user::User; kwargs...)
    json, page_data = paged_request("GET", "/api/v1$(Internals.pid(user))/courses"; kwargs...)
    return Course.(json), page_data
end

"""
    Canvas.create_course(acc::Account; kwargs...) -> Course

Create a new course. Return a [`Course`](@ref).

**Canvas API documentation**:
[Create a new course (`POST /api/v1/accounts/:account_id/courses`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.create)
"""
function create_course(acc::Account; kwargs...)
    json = request("POST", "/api/v1$(Internals.pid(acc))/courses"; kwargs...)
    return Course(json)
end

"""
    Canvas.upload_file(c::Course, file::String; kwargs...) -> File

Upload a file to the course. Return a [`File`](@ref).

**Canvas API documentation**:
[Upload a file to a course (`POST /api/v1/courses/:course_id/files`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.create_file)
"""
function upload_file(c::Course, file::String; kwargs...)
    return Internals.upload_file("/api/v1$(Internals.pid(c))/files", file; kwargs...)
end

"""
    Canvas.users(course::Course; kwargs...) -> Vector{User}, page_data

Return a paginated list of users in the course.
Return a vector of [`User`](@ref)s and a dictionary with page data.

**Canvas API documentation**:
[List users in course (`GET /api/v1/courses/:course_id/users`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.users)
"""
function users(course::Course; kwargs...)
    json, page_data = paged_request("GET", "/api/v1$(Internals.pid(course))/users"; kwargs...)
    return User.(json), page_data
end

"""
    Canvas.recent_students(course::Course; kwargs...) -> Vector{User}, page_data

Return a paginated list of users in the course ordered by how recently they have
logged in. Return a vector of [`User`](@ref)s and a dictionary of page data.

**Canvas API documentation**:
[List recently logged in students (`GET /api/v1/courses/:course_id/recent_students`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.recent_students)
"""
function recent_students(course::Course; kwargs...)
    json, page_data = paged_request("GET", "/api/v1$(Internals.pid(course))/recent_students"; kwargs...)
    return User.(json), page_data
end

"""
    Canvas.user(course::Course, user::User; kwargs...) -> User

Return information about a single user in the course. Return a [`User`](@ref).

**Canvas API documentation**:
[Get single user (`GET /api/v1/courses/:course_id/users/:id`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.user)
"""
function user(course::Course, user::User; kwargs...)
    json = request("GET", "/api/v1$(Internals.pid(course))$(Internals.pid(user))"; kwargs...)
    return User(json)
end

"""
    Canvas.content_share_users(course::Course; kwargs...) -> Vector{User}, page_data

Return a paginated list of users you can share content with.
Return a vector of [`User`](@ref)s and a dictionary of page data.

**Canvas API documentation**:
[Search for content share users (`GET /api/v1/courses/:course_id/content_share_users`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.content_share_users)
"""
function content_share_users(course::Course; kwargs...)
    json, page_data = paged_request("GET", "/api/v1$(Internals.pid(course))/content_share_users"; kwargs...)
    return User.(json), page_data
end

"""
    Canvas.preview_html(course::Course; kwargs...) -> Dict

Preview HTML content processed for this course. Return a dictionary.

**Canvas API documentation**:
[Preview processed html (`POST /api/v1/courses/:course_id/preview_html`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.preview_html)
"""
function preview_html(course::Course; kwargs...)
    json = request("GET", "/api/v1$(Internals.pid(course))/preview_html"; kwargs...)
    return json
end

"""
    Canvas.activity_stream(course::Course; summary::Bool=false, kwargs...) -> Vector{Dict}, page_data

Return the authenticated users course-specific activity stream paginated.

**Request endpoints:**
 - `GET /api/v1/courses/:course_id/activity_stream`
 - `GET /api/v1/courses/:course_id/activity_stream/summary`

**Canvas API documentation**:
 - [Course activity stream (`GET /api/v1/courses/:course_id/activity_stream`)]
   (https://canvas.instructure.com/doc/api/courses#method.courses.activity_stream)
 - [Course activity stream (`GET /api/v1/courses/:course_id/activity_stream/summary`)]
   (https://canvas.instructure.com/doc/api/courses#method.courses.activity_stream_summary)
"""
function activity_stream(course::Course; summary=false, kwargs...)
    endpoint = "/api/v1$(Internals.pid(course))/activity_stream"
    summary && (endpoint *= "/summary")
    json, page_data = paged_request("GET", endpoint; kwargs...)
    return json, page_data
end

"""
    Canvas.todo(course::Course; kwargs...) -> Vector{Dict}, page_data

Return the authenticated users course-specific todo items.

**Canvas API documentation**:
[Course TODO items (`GET /api/v1/courses/:course_id/todo`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.todo_items)
"""
function todo(course::Course; kwargs...)
    json, page_data = paged_request("GET", "/api/v1$(Internals.pid(course))/todo"; kwargs...)
    return json, page_data
end

"""
    Canvas.delete_course(course::Course; kwargs...) -> Dict

Delete or conclude an existing course.

**Canvas API documentation**:
[Delete/Conclude a course (`DELETE /api/v1/courses/:id`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.destroy)
"""
function delete_course(course::Course; kwargs...)
    json = request("DELETE", "/api/v1$(Internals.pid(course))"; kwargs...)
    return json
end

"""
    Canvas.settings(course::Course; kwargs...) -> Dict

Return course settings as a dictionary.

**Canvas API documentation**:
[Get course settings (`GET /api/v1/courses/:course_id/settings`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.api_settings)
"""
function settings(course::Course; kwargs...)
    json = request("GET", "/api/v1$(Internals.pid(course))/settings"; kwargs...)
    return json
end

"""
    Canvas.update_settings(course::Course; kwargs...) -> Dict

Update course settings. Return a dictionary.

**Canvas API documentation**:
[Update course settings (`PUT /api/v1/courses/:course_id/settings`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.update_settings)
"""
function update_settings(course::Course; kwargs...)
    json = request("PUT", "/api/v1$(Internals.pid(course))/settings"; kwargs...)
    return json
end

# TODO: https://canvas.instructure.com/doc/api/courses#method.courses.student_view_student

"""
    Canvas.course(course::Course; kwargs...) -> Course

Return details about a single course specified by id.
Return a [`Course`](@ref).

**Canvas API documentation**:
[Get a single course (`GET /api/v1/courses/:id`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.show)
"""
function course(course::Course; kwargs...)
    json = request("GET", "/api/v1$(Internals.pid(course))"; kwargs...)
    return Course(json)
end

"""
    Canvas.update_course(course::Course; kwargs...) -> Course

Update an existing course and return the new details.
Return a [`Course`](@ref).

**Canvas API documentation**:
[Update a course (`PUT /api/v1/courses/:id`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.update)
"""
function update_course(course::Course; kwargs...)
    json = request("PUT", "/api/v1$(Internals.pid(course))"; kwargs...)
    return Course(json) # Maybe? Or just return json?
end

# TODO: batch update
# """
#     Canvas.update_courses(acc::Account; kwargs...) -> Progress

#  - Request endpoint: `PUT /api/v1/accounts/:account_id/courses`
#  - Canvas API documentation:
#    [*Update courses*](https://canvas.instructure.com/doc/api/courses#method.courses.batch_update)
# """
# function update_courses(acc::Account; kwargs...)
#     json = request("PUT", "/api/v1$(Internals.pid(acc))/courses"; kwargs...)
#     return Progress(json)
# end

"""
    Canvas.reset_content(course::Course; kwargs...) -> Course

Delete the course and create a new equivalent course with no content but all sections
and users moved overs. Return a [`Course`](@ref).

**Canvas API documentation**:
[Reset a course (`POST /api/v1/courses/:course_id/reset_content`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.reset_content)
"""
function reset_content(course::Course; kwargs...)
    json = request("POST", "/api/v1$(Internals.pid(course))/reset_content"; kwargs...)
    return Course(json)
end

"""
    Canvas.effective_due_dates(course::Course; kwargs...) -> Dict

For each assignment in the course, returns each assigned student's ID
and their corresponding due date along with some grading period data.

**Canvas API documentation**:
[Get effective due dates (`GET /api/v1/courses/:course_id/effective_due_dates`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.effective_due_dates)
"""
function effective_due_dates(course::Course; kwargs...)
    json = request("GET", "/api/v1$(Internals.pid(course))/effective_due_dates"; kwargs...)
    return json
end

"""
    Canvas.permissions(course::Course; kwargs...) -> Dict

Return permissions for the authenticated user in the specified course.
Return a dictionary.

**Canvas API documentation**:
[Permissions (`GET /api/v1/courses/:course_id/permissions`)]
(https://canvas.instructure.com/doc/api/courses#method.courses.permissions)
"""
function permissions(course::Course; kwargs...)
    json = request("GET", "/api/v1$(Internals.pid(course))/permissions"; kwargs...)
    return json
end


text:
# -------------------------------------------------------------------------------
# Scores
# -------------------------------------------------------------------------------
# Primal problem
function scores(model::AbstractModel, data::Primal, w::AbstractVector)
    scores(model, data.X, w)
end

function scores(model::AbstractModel, X::AbstractMatrix, w::AbstractVector)
    X * w
end

# Dual problem - train data
function scores(model::PatMat, data::Dual{<:DTrain}, α::AbstractVector, β::AbstractVector)
    .- vec(data.K * vcat(α, β))[data.type.inv_perm]
end

function scores(model::Union{AbstractTopPushK, PatMatNP}, data::Dual{<:DTrain}, α::AbstractVector, β::AbstractVector)
    s = vec(vcat(α, β)'*data.K)
    s[data.ind_β] .*= -1
    return s[data.type.inv_perm]
end

function scores(model::AbstractModel,
                Xtrain::AbstractMatrix,
                ytrain::BitArray{1},
                α::AbstractVector,
                β::AbstractVector;
                kwargs...)

    return scores(model, Dual(model, Xtrain, ytrain; kwargs...), α, β)
end

# Dual problem - validation and test data
function scores(model::AbstractModel, data::Dual{<:Union{DValidation, DTest}}, α::AbstractVector, β::AbstractVector)
    vec(vcat(α, β)'*data.K)
end

function scores(model::AbstractModel,
                Xtrain::AbstractMatrix,
                ytrain::BitArray{1},
                Xvalid::AbstractMatrix,
                yvalid::BitArray{1},
                α::AbstractVector,
                β::AbstractVector;
                kwargs...)

    return scores(model, Dual(model, Xtrain, ytrain, Xvalid, yvalid; kwargs...), α, β)
end

function scores(model::AbstractModel,
                Xtrain::AbstractMatrix,
                ytrain::BitArray{1},
                Xtest::AbstractMatrix,
                α::AbstractVector,
                β::AbstractVector;
                kwargs...)

    return scores(model, Dual(model, Xtrain, ytrain, Xtest; kwargs...), α, β)
end

# Dual problem - load function
function scores(model::AbstractModel, file::AbstractString, α::AbstractVector, β::AbstractVector; kwargs...)

    data = Dual(file; kwargs...)
    s    = scores(model, data, α, β)
    close(data.io)
    return s
end

# -------------------------------------------------------------------------------
# Predict
# -------------------------------------------------------------------------------
# Primal problems
function predict(model::AbstractModel, data::Primal, w::AbstractVector, t::Real)
    scores(model, data, w) .>= t
end

function predict(model::AbstractModel, X::AbstractMatrix, w::AbstractVector, t::Real)
    scores(model, X, w) .>= t
end

# Dual problems
function predict(model::AbstractModel, data::Dual, α::AbstractVector, β::AbstractVector, t::Real)
    scores(model, data, α, β) .>= t
end

function predict(model::AbstractModel,
                 Xtrain::AbstractMatrix,
                 ytrain::BitArray{1},
                 α::AbstractVector,
                 β::AbstractVector,
                 t::Real;
                 kwargs...)

    return scores(model, Dual(model, Xtrain, ytrain; kwargs...), α, β) .>= t
end

function predict(model::AbstractModel,
                 Xtrain::AbstractMatrix,
                 ytrain::BitArray{1},
                 Xvalid::AbstractMatrix,
                 yvalid::BitArray{1},
                 α::AbstractVector,
                 β::AbstractVector,
                 t::Real;
                 kwargs...)

    return scores(model, Dual(model, Xtrain, ytrain, Xvalid, yvalid; kwargs...), α, β) .>= t
end

function predict(model::AbstractModel,
                 Xtrain::AbstractMatrix,
                 ytrain::BitArray{1},
                 Xtest::AbstractMatrix,
                 α::AbstractVector,
                 β::AbstractVector,
                 t::Real;
                 kwargs...)

    return scores(model, Dual(model, Xtrain, ytrain, Xtest; kwargs...), α, β) .>= t
end

function predict(model::AbstractModel, file::AbstractString, α::AbstractVector, β::AbstractVector, t::Real; kwargs...)

    return  scores(model, file, α, β; kwargs...) .>= t
end


text:
immutable PosixPath <: AbstractPath
    parts::Tuple
end

function PosixPath(str::AbstractString)
    str = ASCIIString(str)

    if isempty(str)
        return PosixPath(tuple("."))
    end

    tokenized = split(str, POSIX_PATH_SEPARATOR)
    if isempty(tokenized[1])
        tokenized[1] = POSIX_PATH_SEPARATOR
    end
    return PosixPath(tuple(tokenized...))
end

# The following should be implemented in the concrete types
==(a::PosixPath, b::PosixPath) = parts(a) == parts(b)
Base.string(path::PosixPath) = joinpath(parts(path)...)
parts(path::PosixPath) = path.parts

function isabs(path::PosixPath)
    if parts(path)[1] == POSIX_PATH_SEPARATOR
        return true
    else
        return false
    end
end

drive(path::PosixPath) = ""

function root(path::PosixPath)
    if parts(path)[1] == POSIX_PATH_SEPARATOR
        return POSIX_PATH_SEPARATOR
    else
        return ""
    end
end

function expanduser(path::PosixPath)
    p = parts(path)

    if p[1] == "~"
        if length(p) > 1
            return PosixPath(tuple(homedir(), p[2:end]...))
        else
            return PosixPath(tuple(homedir()))
        end
    end

    return path
end


text:
include("../../src/io/pack.jl")

## simple types (numbers)

type = Int32

vbytes = UInt8[0x1, 0x0, 0x0, 0x0, 0x2, 0x0, 0x0, 0x0]
vector = Int32[1,2]

# tests:
vbytes == pack_vec(vector)
vector == unpack_vec(vbytes, type)

## composite types (tuples)

type = Tuple{Int16, Int32, Int16}
types = fieldtypes(type)

bytes = UInt8[0x1, 0x0, 0x2, 0x0, 0x0, 0x0, 0x3, 0x0]
elem = (Int16(1), Int32(2), Int16(3))

vbytes = repeat(bytes, 10)
vector = Iterators.repeated(elem, 10) |> collect

# tests:

# one element
pack_vec([elem]) == bytes
unpack_vec(bytes, types...) == [elem]

# vector
vbytes == pack_vec(vector)
vector == unpack_vec(vbytes, types...)