text:
# TaylorSeries.jl
#
# Julia module for handling Taylor series of arbitrary but finite order
#
# utils_Taylor1.jl contains the constructors and methods for 1-variable expansions
#
# utils_TaylorN.jl contains the constructors and methods for N-variable expansions
#
# Last modification: 2014.06.07
#
# Luis Benet & David P. Sanders
# UNAM
#

module TaylorSeries

import Base: zero, one
import Base: convert, promote_rule, promote, eltype, length
import Base: real, imag, conj, ctranspose
import Base: rem, mod, mod2pi
import Base: sqrt, exp, log, sin, cos, tan#, square

abstract AbstractSeries{T<:Number,N} <: Number

include("utils_Taylor1.jl")

include("utils_TaylorN.jl")

## The following routines combine Taylor and TaylorN, so they must appear defining
##   Taylor and TaylorN and some of its functionalities

# infostr
infostr{T<:Number}(a::Taylor{T}) =
    string(a.order, "-order Taylor{", T, "}:\n")
infostr{T<:Number}(a::TaylorN{T}) =
    string(a.order, "-order TaylorN{", T, "} in ", a.numVars, " variables:\n")

# pretty_print
function pretty_print{T<:Number}(a::Taylor{T})
    print( infostr(a) )
    z = zero(T)
    space = " "
    a == z && (print(string( space, z)); return)
    strout = space
    ifirst = true
    for i = 0:a.order
        monom = i==0 ? "" : i==1 ? " * x0" : string(" * x0^", i)
        @inbounds c = a.coeffs[i+1]
        c == z && continue
        cadena = numbr2str(c, ifirst)
        strout = string(strout, cadena, monom, space)
        ifirst = false
    end
    println(strout)
end
function pretty_print{T<:Number}(a::TaylorN{T})
    print( infostr(a) )
    z = zero(T)
    space = " "
    a == z && (println(string( space, z)); return)
    varstring = {}
    for ivar=1:a.numVars
        push!(varstring,string(" * x",ivar))
    end
    strout = string(" ")
    ifirst = true
    iIndices = zeros(Int, a.numVars)
    for pos = 1:length(a.coeffs)
        monom = ""
        iIndices = indicesTable[end][pos]
        if pos>1
            for ivar=1:a.numVars
                powivar = iIndices[ivar]
                if powivar == 1
                    monom = string(monom, varstring[ivar])
                elseif powivar > 1
                    monom = string(monom, varstring[ivar], "^", powivar)
                end
            end
        end
        @inbounds c = a.coeffs[pos]
        c == z && continue
        cadena = numbr2str(c, ifirst)
        strout = string(strout, cadena, monom, space)
        ifirst = false
    end
    println(strout)
end
function pretty_print{T<:Number}(a::Union(Array{Taylor{T},1},Array{TaylorN{T},1}))
    for i=1:length(a)
        pretty_print(a[i])
        println("")
    end
end

# make string from a number; for complex numbers, use
function numbr2str{T<:Real}(zz::T, ifirst::Bool=false)
    plusmin = zz > 0 ? "+ " : "- "
    if ifirst
        plusmin = zz > 0 ? " " : "-"
    end
    return string(plusmin, abs(zz))
end
function numbr2str{T}(zz::Complex{T}, ifirst::Bool=false)
    zT = zero(T)
    zz == zero(Complex{T}) && return zT
    zre, zim = reim(zz)
    cadena = ""
    if zre > zT
        if ifirst
            cadena = string("   ( ", abs(zre)," ")
        else
            cadena = string(" + ( ", abs(zre)," ")
        end
        if zim > zT
            cadena = string(cadena, "+ ", abs(zim), " im )")
        elseif zim < zT
            cadena = string(cadena, "- ", abs(zim), " im )")
        else
            cadena = string(cadena, ")")
        end
    elseif zre < zT
        cadena = string(" - ( ", abs(zre), " ")
        if zim > zT
            cadena = string(cadena, "- ", abs(zim), " im )")
        elseif zim < zT
            cadena = string(cadena, "+ ", abs(zim), " im )")
        else
            cadena = string(cadena, ")")
        end
    else
        if zim > zT
            if ifirst
                cadena = string("   ( ", abs(zim), " im )")
            else
                cadena = string(" + ( ", abs(zim), " im )")
            end
        else
            cadena = string(" - ( ", abs(zim), " im )")
        end
    end
    return cadena
end

# evalTaylor(TaylorN, Array{Taylor,1})
function evalTaylor{T<:Number,S<:Number}(a::TaylorN{T}, vT::Array{Taylor{S},1})
    numVars = NUMVARS[end]
    @assert length(vT) == numVars
    R = promote_type(T,S)
    order = length(vT[1])
    sumaT = Taylor(zero(R), order)
    z = zero(sumaT)
    nCoefTot = sizeTable[end]
    iIndices = zeros(Int, numVars)
    for pos = nCoefTot:-1:1
        @inbounds iIndices[1:end] = indicesTable[end][pos]
        a.coeffs[pos] == zero(T) && continue
        val = Taylor(a.coeffs[pos], order)
        for k = 1:numVars
            @inbounds val = val*(vT[k])^iIndices[k]
        end
        sumaT += val
    end
    return sumaT
end

## Exports to Taylor and TaylorN ##
export Taylor, diffTaylor, integTaylor, evalTaylor, deriv, pretty_print
#
export TaylorN
export set_maxOrder, get_maxOrder, set_numVars, get_numVars

end


text:
export TypicalSyndrome, NaturalAllAtOnce, NaturalInterleaved, CompactAllAtOnce,
    CompactInterleaved,
    make_noise_model_for_paper

# The different circuits evaluated
struct TypicalSyndrome <: AbstractFastSyndrome end
struct NaturalAllAtOnce <: AbstractFastSyndrome end
struct NaturalInterleaved <: AbstractFastSyndrome end
struct CompactAllAtOnce <: AbstractFastSyndrome end
struct CompactInterleaved <: AbstractFastSyndrome end

### Make the noise model used in the paper
p0 = 0.1
starting_model = Dict{Symbol, Float64}(
    :t1_t => 100_000,  # ns
    :t1_c => 1_000_000,  # ns
    :dur_t => 50,  # ns
    :dur_tt => 200,  # ns
    :dur_tc => 200,  # ns
    :dur_loadstore => 150,  # ns
    :dur_meas => 200,  # ns
    :p_t => p0 / 10,
    :p_tt => p0,
    :p_tc => p0,
    :p_loadstore => p0,
    :p_meas => p0,
    :cavity_depth => 10,
)
const_keys = [:dur_tt, :dur_t, :dur_tc, :dur_loadstore, :dur_meas,
              :cavity_depth]
coherence_keys = [:t1_t, :t1_c]
error_rate_keys = [:p_t, :p_tt, :p_tc, :p_loadstore, :p_meas]
sensitivity_base_p = 2e-3

function make_noise_model_for_paper(base_error::Float64, override_pairs=())
    error_factor = base_error / p0
    actual_model = Dict{Symbol, Float64}()
    for k in const_keys
        actual_model[k] = starting_model[k]
    end
    for k in coherence_keys
        actual_model[k] = starting_model[k] / error_factor
    end
    for k in error_rate_keys
        actual_model[k] = starting_model[k] * error_factor
    end
    for (k, v) in override_pairs
        actual_model[k] = v
    end
    NoiseModel(actual_model)
end

### Util functions
"""
    combine_flip_probs(p1, p2, ...)

Return the probability that an odd number of events occur
with the given probabilities.
"""
function combine_flip_probs(p_flip::Float64...)
    p_total = 0
    for p in p_flip
        p_total = p_total + p - 2*p_total*p
    end
    p_total
end

"""
    combine_error_probs(p1, p2, ...)

Return the probability that at least one event occurs
with the given probabilities.
"""
function combine_error_probs(p_err::Float64...)
    1 - prod(1 .- p_err)
end

function coherence_error(t1, duration)
    1 - MathConstants.e ^ (-duration / t1)
end

### Calculate edge weights and simulation noise parameters

function calculate_qubit_error_single_pauli(model::NoiseModel;
        t_t::Float64 = 0.0,
        t_c::Float64 = 0.0,
        n_t::Int = 0,
        n_tt::Int = 0,
        n_tc::Int = 0,
        n_loadstore::Int = 0,
        n_meas::Int = 0,
    )
    combine_flip_probs(
        2/3 * coherence_error(model.t1_t, t_t),
        2/3 * coherence_error(model.t1_c, t_c),
        repeat([2/3*model.p_t], n_t)...,
        repeat([8/15*model.p_tt], n_tt)...,
        repeat([8/15*model.p_tc], n_tc)...,
        repeat([2/3*model.p_loadstore], n_loadstore)...,
        repeat([model.p_meas], n_meas)...,
    )
end

# Typical Syndrome
# Z: One CNOT per data, then measure
# X: Hadamard, one CNOT (reverse dir) per data, hadamard, measure
function matching_space_edge(::TypicalSyndrome, model::NoiseModel,
                             is_x::Bool, dist::Int, first_layer::Bool)
    # Mainly 4 CNOTs; one per plaquette
    p = calculate_qubit_error_single_pauli(model,
        t_t = 4*model.dur_tt + 2*model.dur_t + model.dur_meas,
        n_t = 0,
        n_tt = 4,
    )
    -log(p)
end
function matching_time_edge(::TypicalSyndrome, model::NoiseModel,
                            is_x::Bool, dist::Int, first_layer::Bool)
    # Mainly 4 CNOTs
    # Also two Hadamards if is X
    p = calculate_qubit_error_single_pauli(model,
        t_t = 4*model.dur_tt + (is_x ? 2*model.dur_t : 0),
        n_t = (is_x ? 2 : 0),
        n_tt = 4,
        n_meas = 1,
    )
    -log(p)
end
function simulation_noise_parameters(::TypicalSyndrome, model::NoiseModel,
                                     ctx::CodeDistanceSim)
    t_data = 4*model.dur_tt + 2*model.dur_t + model.dur_meas
    t_anc_z = 4*model.dur_tt
    t_anc_x = 4*model.dur_tt + 2*model.dur_t
    Dict{Symbol, Float64}(
        :p_data_layer1 => coherence_error(model.t1_t, t_data),
        :p_data => coherence_error(model.t1_t, t_data),
        :p_anc_z => combine_error_probs(
            calculate_qubit_error_single_pauli(model,
                t_t = t_anc_z,
                n_t = 0),
            model.p_meas),
        :p_anc_x => combine_error_probs(
            calculate_qubit_error_single_pauli(model,
                t_t = t_anc_x,
                n_t = 2),
            model.p_meas),
        :p_cnot1 => model.p_tt,
        :p_cnot => model.p_tt,
    )
end

# NaturalAllAtOnce
# Z: One CNOT per data, then measure
# X: Hadamard, one CNOT (reverse dir) per data, hadamard, measure
# Extra store-load error in first layer and stored for (cavity-1)*d layers
function matching_space_edge(::NaturalAllAtOnce, model::NoiseModel,
                             is_x::Bool, m_dist::Int, first_layer::Bool)
    # Estimate of error on data qubit
    t_round = 4*model.dur_tt + 2*model.dur_t + model.dur_meas
    p = if first_layer
        calculate_qubit_error_single_pauli(model,
            t_t = t_round,
            t_c = ((model.cavity_depth-1)
                   * (m_dist * t_round + 2*model.dur_loadstore)),
            n_t = 0,
            n_tt = 4,
            n_loadstore = 2,
        )
    else
        calculate_qubit_error_single_pauli(model,
            t_t = t_round,
            t_c = 0.0,
            n_t = 0,
            n_tt = 4,
            n_loadstore = 0,
        )
    end
    -log(p)
end
function matching_time_edge(::NaturalAllAtOnce, model::NoiseModel,
                            is_x::Bool, m_dist::Int, first_layer::Bool)
    # Estimate of error on ancilla qubit
    # Same as TypicalSyndrome
    p = calculate_qubit_error_single_pauli(model,
        t_t = 4*model.dur_tt + (is_x ? 2*model.dur_t : 0),
        n_t = (is_x ? 2 : 0),
        n_tt = 4,
        n_meas = 1,
    )
    -log(p)
end
function simulation_noise_parameters(::NaturalAllAtOnce, model::NoiseModel,
                                     ctx::CodeDistanceSim)
    t_round = 4*model.dur_tt + 2*model.dur_t + model.dur_meas
    t_t_data = t_round
    t_t_anc_z = 4*model.dur_tt
    t_t_anc_x = 4*model.dur_tt + 2*model.dur_t
    t_c1_data = ((model.cavity_depth-1)
                 * (ctx.m_dist * t_round + 2*model.dur_loadstore))
    Dict{Symbol, Float64}(
        :p_data_layer1 => combine_error_probs(
            coherence_error(model.t1_c, t_c1_data),
            model.p_loadstore,
            model.p_loadstore,
            coherence_error(model.t1_t, t_t_data)),
        :p_data => coherence_error(model.t1_t, t_t_data),
        :p_anc_z => combine_error_probs(
            calculate_qubit_error_single_pauli(model,
                t_t = t_t_anc_z,
                n_t = 0),
            model.p_meas),
        :p_anc_x => combine_error_probs(
            calculate_qubit_error_single_pauli(model,
                t_t = t_t_anc_x,
                n_t = 2),
            model.p_meas),
        :p_cnot1 => model.p_tt,
        :p_cnot => model.p_tt,
    )
end

# NaturalInterleaved
# Z: One CNOT per data, then measure
# X: Hadamard, one CNOT (reverse dir) per data, hadamard, measure
# Extra store-load error in *every* layer and stored for (cavity-1) layers
function matching_space_edge(::NaturalInterleaved, model::NoiseModel,
                             is_x::Bool, m_dist::Int, first_layer::Bool)
    # Estimate of error on data qubit
    t_round = 4*model.dur_tt + 2*model.dur_t + model.dur_meas
    p = calculate_qubit_error_single_pauli(model,
        t_t = 4*model.dur_tt,
        t_c = ((model.cavity_depth-1) * (t_round + 2*model.dur_loadstore)
               + 2*model.dur_t + model.dur_meas),
        n_t = 0,
        n_tt = 4,
        n_loadstore = 2,
    )
    -log(p)
end
function matching_time_edge(::NaturalInterleaved, model::NoiseModel,
                            is_x::Bool, m_dist::Int, first_layer::Bool)
    # Estimate of error on ancilla qubit
    # Same as TypicalSyndrome
    p = calculate_qubit_error_single_pauli(model,
        t_t = 4*model.dur_tt + (is_x ? 2*model.dur_t : 0),
        n_t = (is_x ? 2 : 0),
        n_tt = 4,
        n_meas = 1,
    )
    -log(p)
end
function simulation_noise_parameters(::NaturalInterleaved, model::NoiseModel,
                                     ctx::CodeDistanceSim)
    t_round = 4*model.dur_tt + 2*model.dur_t + model.dur_meas
    t_t_data = 4*model.dur_tt
    t_t_anc_z = 4*model.dur_tt
    t_t_anc_x = 4*model.dur_tt + 2*model.dur_t
    t_c1_data = ((model.cavity_depth-1) * (t_round + 2*model.dur_loadstore)
                 + 2*model.dur_t + model.dur_meas)
    p_data = combine_error_probs(
        coherence_error(model.t1_c, t_c1_data),
        model.p_loadstore,
        model.p_loadstore,
        coherence_error(model.t1_t, t_t_data))
    Dict{Symbol, Float64}(
        :p_data_layer1 => p_data,
        :p_data => p_data,
        :p_anc_z => combine_error_probs(
            calculate_qubit_error_single_pauli(model,
                t_t = t_t_anc_z,
                n_t = 0),
            model.p_meas),
        :p_anc_x => combine_error_probs(
            calculate_qubit_error_single_pauli(model,
                t_t = t_t_anc_x,
                n_t = 2),
            model.p_meas),
        :p_cnot1 => model.p_tt,
        :p_cnot => model.p_tt,
    )
end

# CompactAllAtOnce
# CNOTs are staggered and interleaved to share qubits between data and ancilla
# so take almost twice as long per layer
# Z: One CNOT per data, then measure
# X: Hadamard, one CNOT (reverse dir) per data, hadamard, measure
# Extra store-load error each layer and stored for (cavity-1)*d layers
function matching_space_edge(::CompactAllAtOnce, model::NoiseModel,
                             is_x::Bool, m_dist::Int, first_layer::Bool)
    # Estimate of error on data qubit
    t_round = (6*model.dur_tt + model.dur_tc + 2*model.dur_t + model.dur_meas
               + 2*model.dur_loadstore)
    p = calculate_qubit_error_single_pauli(model,
            t_t = 3*model.dur_tt,
            t_c = ((first_layer
                        ? (model.cavity_depth-1) * (m_dist * t_round)
                        : 0)
                   + (t_round - 3*model.dur_tt)),
            n_t = 0,
            n_tt = 3,
            n_tc = 1,
            n_loadstore = 2,
        )
    -log(p)
end
function matching_time_edge(::CompactAllAtOnce, model::NoiseModel,
                            is_x::Bool, m_dist::Int, first_layer::Bool)
    # Estimate of error on ancilla qubit
    # One CNOT is transmon-to-cavity
    p = calculate_qubit_error_single_pauli(model,
        t_t = 3*model.dur_tt + model.dur_tc + (is_x ? 2*model.dur_t : 0),
        n_t = (is_x ? 2 : 0),
        n_tt = 3,
        n_tc = 1,
        n_meas = 1,
    )
    -log(p)
end
function simulation_noise_parameters(::CompactAllAtOnce, model::NoiseModel,
                                     ctx::CodeDistanceSim)
    t_round = (6*model.dur_tt + model.dur_tc + 2*model.dur_t + model.dur_meas
               + 2*model.dur_loadstore)
    t_t_data = 3*model.dur_tt
    t_t_anc_z = 3*model.dur_tt + model.dur_tc
    t_t_anc_x = 3*model.dur_tt + model.dur_tc + 2*model.dur_t
    t_c_data = t_round - 3*model.dur_tt
    t_c1_data = t_c_data + (model.cavity_depth-1) * (ctx.m_dist * t_round)
    Dict{Symbol, Float64}(
        :p_data_layer1 => combine_error_probs(
            coherence_error(model.t1_c, t_c1_data),
            model.p_loadstore,
            model.p_loadstore,
            coherence_error(model.t1_t, t_t_data)),
        :p_data => combine_error_probs(
            coherence_error(model.t1_c, t_c_data),
            model.p_loadstore,
            model.p_loadstore,
            coherence_error(model.t1_t, t_t_data)),
        :p_anc_z => combine_error_probs(
            calculate_qubit_error_single_pauli(model,
                t_t = t_t_anc_z,
                n_t = 0),
            model.p_meas),
        :p_anc_x => combine_error_probs(
            calculate_qubit_error_single_pauli(model,
                t_t = t_t_anc_x,
                n_t = 2),
            model.p_meas),
        :p_cnot1 => model.p_tc,
        :p_cnot => model.p_tt,
    )
end

# CompactInterleaved
# CNOTs are staggered and interleaved to share qubits between data and ancilla
# so take almost twice as long per layer
# Z: One CNOT per data, then measure
# X: Hadamard, one CNOT (reverse dir) per data, hadamard, measure
# Extra store-load error each layer and stored for (cavity-1) layers between
function matching_space_edge(::CompactInterleaved, model::NoiseModel,
                             is_x::Bool, m_dist::Int, first_layer::Bool)
    # Estimate of error on data qubit
    t_round = (6*model.dur_tt + model.dur_tc + 2*model.dur_t + model.dur_meas
               + 2*model.dur_loadstore)
    p = calculate_qubit_error_single_pauli(model,
            t_t = 3*model.dur_tt,
            t_c = (model.cavity_depth-1) * t_round + (t_round - 3*model.dur_tt),
            n_t = 0,
            n_tt = 3,
            n_tc = 1,
            n_loadstore = 2,
        )
    -log(p)
end
function matching_time_edge(::CompactInterleaved, model::NoiseModel,
                            is_x::Bool, m_dist::Int, first_layer::Bool)
    # Estimate of error on ancilla qubit
    # Same as CompactAllAtOnce
    p = calculate_qubit_error_single_pauli(model,
        t_t = 3*model.dur_tt + model.dur_tc + (is_x ? 2*model.dur_t : 0),
        n_t = (is_x ? 2 : 0),
        n_tt = 3,
        n_tc = 1,
        n_meas = 1,
    )
    -log(p)
end
function simulation_noise_parameters(::CompactInterleaved, model::NoiseModel,
                                     ctx::CodeDistanceSim)
    t_round = (6*model.dur_tt + model.dur_tc + 2*model.dur_t + model.dur_meas
               + 2*model.dur_loadstore)
    t_t_data = 3*model.dur_tt
    t_t_anc_z = 3*model.dur_tt + model.dur_tc
    t_t_anc_x = 3*model.dur_tt + model.dur_tc + 2*model.dur_t
    t_c_data = (model.cavity_depth-1) * t_round + (t_round - 3*model.dur_tt)
    p_data = combine_error_probs(
        coherence_error(model.t1_c, t_c_data),
        model.p_loadstore,
        model.p_loadstore,
        coherence_error(model.t1_t, t_t_data))
    Dict{Symbol, Float64}(
        :p_data_layer1 => p_data,
        :p_data => p_data,
        :p_anc_z => combine_error_probs(
            calculate_qubit_error_single_pauli(model,
                t_t = t_t_anc_z,
                n_t = 0),
            model.p_meas),
        :p_anc_x => combine_error_probs(
            calculate_qubit_error_single_pauli(model,
                t_t = t_t_anc_x,
                n_t = 2),
            model.p_meas),
        :p_cnot1 => model.p_tc,
        :p_cnot => model.p_tt,
    )
end

ext: jl

max_stars_repo_name: cduck/VLQ

max_issues_repo_name: cduck/VLQ

max_forks_repo_name: cduck/VLQ


text:
using Knet
include("readdata.jl")

# Similar to model02, trying more fully connected layers

@knet function model03(x; T=51, D=4, W=3, o...)
    c1 = cbf73(x; o..., out=64)
    c2 = cbf73(c1; o..., out=128)
    c3 = cbf73(c2; o..., out=256)
    c = pool(c3; o..., window=T)
    f1 = wbf73(c; o..., out=128)
    f2 = wbf73(f1; o..., out=128)
    return wbf73(f2; o..., out=2, f=:soft)
end

function train03(; N=2^13, lr=0.001, adam=true, nbatch=128, o...)
    global f = compile(:model03; o...)
    setp(f, lr=lr, adam=adam)
    global d = GulayseData(nbatch)
    sloss = zloss = -1; n = nextn = 1
    for (x,ygold) in d
        ypred = forw(f, x)
        sl = softloss(ypred,ygold); sloss = (n==1 ? sl : 0.99 * sloss + 0.01 * sl)
        zl = zeroone(ypred,ygold);  zloss = (n==1 ? zl : 0.99 * zloss + 0.01 * zl)
        n==nextn && (println((n,sloss,1-zloss)); nextn*=2)
        back(f, ygold, softloss)
        update!(f)
        (n += 1) > N && break
    end
end

train03()

max_stars_repo_name: JuliaTagBot/BindingPrediction

max_issues_repo_name: JuliaTagBot/BindingPrediction

max_forks_repo_name: JuliaTagBot/BindingPrediction


text:
"""
generate_overlay_code(prn::Int64)

Generate the overlay code for PRNs 1-63 using the S1 register only.
"""
function generate_overlay_code(prn::Int64)
    S1 = overlay_parms[prn]["init"]

end

"""
overlay_parms

Contains the 11-bit shift register taps and initial conditions
for each PRN.
"""
const overlay_parms = Dict( 1 => Dict("taps" => (3, 6, 9, 11),
                                      "init" => (1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0)),
                            2 => Dict("taps" => (4, 8, 9, 11),
                                      "init" => (1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0)),
                            3 => Dict("taps" => (6, 8, 9, 11),
                                      "init" => (0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1)),
                            4 => Dict("taps" => (1, 8, 9, 11),
                                      "init" => (1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1)),
                            5 => Dict("taps" => (1, 2, 3, 8, 10, 11),
                                      "init" => (1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0)),
                            6 => Dict("taps" => (5, 6, 10, 11),
                                      "init" => (1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0)),
                            7 => Dict("taps" => (3, 5, 6, 7, 10, 11),
                                      "init" => (0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0)),
                            8 => Dict("taps" => (6, 8, 10, 11),
                                      "init" => (0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0)),
                            9 => Dict("taps" => (2, 7, 10, 11),
                                      "init" => (1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1)),
                           10 => Dict("taps" => (2, 3, 4, 7, 10, 11),
                                      "init" => (0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1)),
                           11 => Dict("taps" => (3, 5, 6, 7, 8, 9, 10, 11),
                                      "init" => (0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1)),
                           12 => Dict("taps" => (1, 4, 7, 8, 10, 11),
                                      "init" => (0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0)),
                           13 => Dict("taps" => (1, 3, 4, 6, 7, 8, 10, 11),
                                      "init" => (1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0)),
                           14 => Dict("taps" => (1, 2, 4, 7, 8, 9, 10, 11),
                                      "init" => (0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0)),
                           15 => Dict("taps" => (1, 2, 4, 5, 7, 8, 9, 11),
                                      "init" => (1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0)),
                           16 => Dict("taps" => (3, 5, 9, 11),
                                      "init" => (1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0)),
                           17 => Dict("taps" => (2, 4, 5, 7, 8, 9, 10, 11),
                                      "init" => (1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1)),
                           18 => Dict("taps" => (2, 4, 6, 7, 10, 11),
                                      "init" => (1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0)),
                           19 => Dict("taps" => (2, 4, 5, 6, 7, 11),
                                      "init" => (1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1)),
                           20 => Dict("taps" => (2, 5, 6, 7, 8, 11),
                                      "init" => (0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0)),
                           21 => Dict("taps" => (1, 3, 4, 7, 8, 9, 10, 11),
                                      "init" => (0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0)),
                           22 => Dict("taps" => (1, 2, 5, 6, 7, 8, 10, 11),
                                      "init" => (0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0)),
                           23 => Dict("taps" => (2, 3, 4, 5, 8, 11),
                                      "init" => (0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0)),
                           24 => Dict("taps" => (2, 4, 7, 11),
                                      "init" => (1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1)),
                           25 => Dict("taps" => (1, 4, 5, 9, 10, 11),
                                      "init" => (1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1)),
                           26 => Dict("taps" => (1, 4, 8, 11),
                                      "init" => (1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1)),
                           27 => Dict("taps" => (3, 5, 7, 8, 10, 11),
                                      "init" => (1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0)),
                           28 => Dict("taps" => (4, 5, 6, 11),
                                      "init" => (0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)),
                           29 => Dict("taps" => (1, 2, 3, 4, 7, 9, 10, 11),
                                      "init" => (1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1)),
                           30 => Dict("taps" => (1, 3, 4, 5, 8, 11),
                                      "init" => (0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0)),
                           31 => Dict("taps" => (1, 2, 3, 4, 5, 8, 9, 11),
                                      "init" => (1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0)),
                           32 => Dict("taps" => (1, 4, 5, 6, 10, 11),
                                      "init" => (0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1)),
                           33 => Dict("taps" => (1, 4, 7, 9, 10, 11),
                                      "init" => (1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0)),
                           34 => Dict("taps" => (1, 4, 6, 7, 10, 11),
                                      "init" => (1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1)),
                           35 => Dict("taps" => (2, 4, 6, 9, 10, 11),
                                      "init" => (1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0)),
                           36 => Dict("taps" => (2, 3, 4, 9, 10, 11),
                                      "init" => (0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0)),
                           37 => Dict("taps" => (5, 6, 7, 11),
                                      "init" => (0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1)),
                           38 => Dict("taps" => (1, 3, 5, 6, 7, 11),
                                      "init" => (0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1)),
                           39 => Dict("taps" => (1, 2, 6, 11),
                                      "init" => (1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0)),
                           40 => Dict("taps" => (2, 3, 4, 6, 7, 8, 9, 11),
                                      "init" => (0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0)),
                           41 => Dict("taps" => (5, 6, 7, 8, 10, 11),
                                      "init" => (0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1)),
                           42 => Dict("taps" => (3, 4, 5, 9, 10, 11),
                                      "init" => (1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0)),
                           43 => Dict("taps" => (1, 4, 5, 6, 8, 11),
                                      "init" => (0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1)),
                           44 => Dict("taps" => (2, 3, 5, 6, 7, 8, 9, 11),
                                      "init" => (1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1)),
                           45 => Dict("taps" => (1, 2, 4, 6, 10, 11),
                                      "init" => (0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0)),
                           46 => Dict("taps" => (3, 4, 5, 7, 8, 11),
                                      "init" => (0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1)),
                           47 => Dict("taps" => (3, 6, 8, 11),
                                      "init" => (1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0)),
                           48 => Dict("taps" => (1, 3, 4, 6, 8, 11),
                                      "init" => (0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0)),
                           49 => Dict("taps" => (1, 2, 3, 5, 6, 7, 9, 11),
                                      "init" => (0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0)),
                           50 => Dict("taps" => (1, 2, 7, 8, 9, 11),
                                      "init" => (0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)),
                           51 => Dict("taps" => (1, 3, 4, 5, 7, 8, 10, 11),
                                      "init" => (1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1)),
                           52 => Dict("taps" => (1, 3, 5, 6, 10, 11),
                                      "init" => (0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0)),
                           53 => Dict("taps" => (2, 4, 5, 6, 8, 9, 10, 11),
                                      "init" => (1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0)),
                           54 => Dict("taps" => (1, 2, 6, 9, 10, 11),
                                      "init" => (1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1)),
                           55 => Dict("taps" => (3, 7, 10, 11),
                                      "init" => (1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1)),
                           56 => Dict("taps" => (4, 6, 7, 11),
                                      "init" => (1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1)),
                           57 => Dict("taps" => (7, 9, 10, 11),
                                      "init" => (0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0)),
                           58 => Dict("taps" => (3, 5, 8, 11),
                                      "init" => (1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0)),
                           59 => Dict("taps" => (3, 8, 9, 11),
                                      "init" => (0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0)),
                           60 => Dict("taps" => (5, 6, 9, 11),
                                      "init" => (0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0)),
                           61 => Dict("taps" => (5, 9, 10, 11),
                                      "init" => (0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0)),
                           62 => Dict("taps" => (1, 2, 3, 4, 7, 8, 10, 11),
                                      "init" => (1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0)),
                           63 => Dict("taps" => (1, 2, 3, 4, 5, 6, 8, 11),
                                      "init" => (0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1)))


text:
# NOTE: this is mostly just copied from wikipedia
function gmres_cycle(A::SparseMatrixCSC, x0::Vector{Float64}, b::Vector{Float64}, m::Int64; err_thresh = 1e-20)

    # compute initial residual, norm(b) and initial error
    r = b - A * x0
    bnorm = norm(b)
    err = norm(r) / bnorm

    n = length(x0)
    Q = zeros(Float64, (n,m+1))
    H = zeros(Float64, (m+2,m+1))
    sn = zeros(Float64, m)
    cn = zeros(Float64, m)
    e1 = zeros(Float64, m+1)
    e1[1] = 1.0

    beta = norm(r) * e1
    Q[:,1] = r / norm(r)

    # begin iteration
    nk = 0
    for k=1:m
        nk += 1

        # do arnoldi iteration
        q = A * Q[:,k]
        for i=1:k
            H[i,k] = transpose(q) * Q[:,i]
            q = q - H[i,k] * Q[:,i]
        end
        qnorm = norm(q)
        H[k+1,k] = qnorm
        Q[:,k+1] = q / qnorm

        # apply givens rotation
        for i=1:k-1
            temp = cn[i] * H[i,k] + sn[i] * H[i+1,k]
            H[i+1,k] = -sn[i] * H[i,k] + cn[i] * H[i+1,k]
            H[i,k] = temp
        end

        # update sn and cn
        t = sqrt(H[k,k]^2 + H[k+1,k]^2)
        cn[k] = H[k,k] / t
        sn[k] = H[k+1,k] / t

        # eliminate H[k+1, k]
        H[k,k] = cn[k] * H[k,k] + sn[k] * H[k+1,k]
        H[k+1,k] = 0.0

        # update residual
        beta[k+1] = -sn[k]*beta[k]
        beta[k] = cn[k]*beta[k]
        err = abs(beta[k+1]) / bnorm

        if err < err_thresh
            break
        end
    end

    # compute correction to x
    y = H[1:nk,1:nk] \ beta[1:nk]
    x = x0 + Q[:,1:nk] * y
    return x, err
end

function gmres_solve(A::SparseMatrixCSC, x0::Vector{Float64}, b::Vector{Float64}, m::Int64; maxiters=10000, err_thresh = 1e-8)

    x = x0[:]
    err = 1.0
    iteration = 0

    while (err > err_thresh)
        x, err = gmres_cycle(A, x, b, m, err_thresh=err_thresh)
        iteration += 1

        if iteration > maxiters
            @warn "Iteration limit reached, defaulting to direct solver."
            return A \ b
        end
    end

    return x
end


text:
include("pricing_model.jl")

Random.seed!(1)
runPxExperiment(500)

max_stars_repo_name: evanmunro/personalized-policy

max_issues_repo_name: evanmunro/personalized-policy

max_forks_repo_name: evanmunro/personalized-policy


text:
function add_variable!(jm::JuMP.Model, name::Symbol, lower_bound::Union{Float64,Nothing}=nothing)
    if lower_bound===nothing
        jm[name] = JuMP.@variable(jm, base_name=string(name))
    else
        jm[name] = JuMP.@variable(jm, base_name=string(name), lower_bound=lower_bound)
    end
end

function add_parameter_to_jump!(jm, parameter::ScalarParameter)
    jmp_p = @eval(JuMP.@NLparameter($jm, $(parameter.name) == $(parameter.value)))
    jm[parameter.name] = jmp_p
end

function add_parameter_to_jump!(jn, parameter::IndexedParameter)
    jm[parameter.name] = @eval(JuMP.@NLparameter($jm, [$( ( :($(gensym())=$i) for i in parameter.indices)... )], base_name=string($(QuoteNode(parameter.name)))))
end

function add_sector_to_jump!(jm, sector::ScalarSector)
    add_variable!(jm, sector.name, 0.)
end

function add_sector_to_jump!(jm, sector::IndexedSector)
        jm[sector.name] = @eval(JuMP.@variable($jm, [$( ( :($(gensym())=$i) for i in sector.indices)... )], base_name=string($(QuoteNode(sector.name))), lower_bound=0.))
end

function add_commodity_to_jump!(jm, commodity::ScalarCommodity)
    add_variable!(jm, commodity.name, 0.)
end

function add_commodity_to_jump!(jm, commodity::IndexedCommodity)
    jm[commodity.name] = @eval(JuMP.@variable($jm, [$( ( :($(gensym())=$i) for i in commodity.indices)... )], base_name=string($(QuoteNode(commodity.name))), lower_bound=0.))
end

function add_consumer_to_jump!(jm, consumer::ScalarConsumer)
    add_variable!(jm, consumer.name, 0.)
end

function add_consumer_to_jump!(jm, consumer::IndexedConsumer)
    jm[consumer.name] = @eval(JuMP.@variable($jm, [$( ( :($(gensym())=$i) for i in consumer.indices)... )], base_name=string($(QuoteNode(consumer.name))), lower_bound=0.))
end

function build_variables!(m, jm)
    # Add all parameters

    for p in m._parameters
        add_parameter_to_jump!(jm, p)
        # jmp_p = @eval(JuMP.@NLparameter($jm, $(p.name) == $(p.value)))
        # jm[p.name] = jmp_p
    end

    # Add all required variables

    for s in m._sectors
        add_sector_to_jump!(jm, s)
    end

    for c in m._commodities
        add_commodity_to_jump!(jm, c)
    end

    # Add compensated supply variables
    for s in m._productions
        for o in s.outputs
            add_variable!(jm, get_comp_supply_name(o))
        end
    end

    # Add compensated demand variables
    for s in m._productions
        for i in s.inputs
            add_variable!(jm, get_comp_demand_name(i))
        end
    end

    for c in m._consumers
        add_consumer_to_jump!(jm, c)
    end

    # Add final demand variables
    for demand_function in m._demands
        for demand in demand_function.demands
            add_variable!(jm, get_final_demand_name(demand))
        end
    end
end


text:
@testset "Isogenies" begin
  K = GF(7)
  E1 = EllipticCurve(K, [1, 2, 3, 4, 5])
  E2 = EllipticCurve(K, [1, 2, 3, 1, 1])
  phi = @inferred isogeny_from_kernel(E1, division_polynomial_univariate(E1,3)[1])
  @test @inferred domain(phi) == E1
  @test @inferred codomain(phi) == E2
  @test is_isomorphic(E1, codomain(phi))
  @test @inferred image(phi, E1([1,1])) == E2([6,2])

  K = GF(3)
  E1 = EllipticCurve(K, [1, 0, 1, 1, 0])
  E2 = EllipticCurve(K, [1, 0, 1, 0, 2])
  Kx, x = PolynomialRing(K);
  f = x+2
  phi = @inferred isogeny_from_kernel(E1, f)
  @test @inferred domain(phi) == E1
  @test @inferred codomain(phi) == E2
  @test @inferred isinfinite(image(phi, E1([1,2])))

  E1 = EllipticCurve([1, 2, 3, 4, 5])
  E2 = EllipticCurve([1, 2, 3, 979//16, 19067//64])
  phi = @inferred isogeny_from_kernel(E1, division_polynomial_univariate(E1,2)[1])
  @test @inferred domain(phi) == E1
  @test @inferred codomain(phi) == E2
  @test is_isomorphic(E1, codomain(phi))

  K= GF(2,4)
  a = gen(K)
  E1 = EllipticCurve(K,[a^2,1-a,1,0,a])
  E2 = EllipticCurve(K,[a^2,1-a,1,a^8,1])
  f = division_polynomial_univariate(E1,5)[1]
  phi = @inferred isogeny_from_kernel(E1, f)
  @test @inferred domain(phi) == E1
  @test @inferred codomain(phi) == E2
  @test @inferred image(phi, E1([a^2,a^2])) == E2([1,a^14])
end


text:
using Random

"""

    sort_biosystem_compounds(compound_structs,
        biosystem_compounds,
        sortkey,
        zero_mass_behavior)

Return a vector of (compound,mass) sorted by increasing mass.
Always shuffles the order of the zero-mass compounds.

...
# Arguments
- `compound_structs::Compounds`: compounds with metadata
- `biosystem_compounds::IDs`: compounds to sort
- `sortkey::Symbol=:exact_mass`: how to sort compounds. Can also sort by `:mol_weight`
- `zero_mass_behavior::String="end"`: what to do with compounds that have no mass. Can also be `"random"`
...
"""
function sort_biosystem_compounds(
    compound_structs::Compounds,
    biosystem_compounds::IDs,
    sortkey::Symbol=:exact_mass,
    zero_mass_behavior::String="end")

    cpd_masses = [(i,getproperty(compound_structs[i],sortkey)) for i in biosystem_compounds if getproperty(compound_structs[i],sortkey)!=0]
    cpd_zero_masses = [(i,getproperty(compound_structs[i],sortkey)) for i in biosystem_compounds if getproperty(compound_structs[i],sortkey)==0]

    sort!(cpd_masses, by= x->x[2])

    if zero_mass_behavior=="end"
        return vcat(cpd_masses,shuffle(cpd_zero_masses))

    elseif zero_mass_behavior=="random"
        for tup in cpd_zero_masses
            insert!(cpd_masses,rand(1:length(cpd_masses)),tup) # this will permit insertion at longer indices as cpd_masses grows
        end
        return cpd_masses
    end
end

"""
    mix_it_up(tuples,beta,n_swaps)

Swap random tuples within `tuples`.
Based on formula from Handorf et al. 2008.
"""
function mix_it_up!(
    tuples::Vector{<:Tuple{String,Real}},
    beta::Real,
    n_swaps::Int,
    rng::AbstractRNG=Random.default_rng()) #where {R <: Real}

    for _ in 1:n_swaps
        i,j = rand(1:length(tuples),2) ## 2 random indices
        if swap_random(tuples[i],tuples[j],beta,rng) == true
            tuples[i], tuples[j] = tuples[j], tuples[i] ## changes them simultaneously
        end
    end
end

"""
    swap_random(ti,tj,beta)

Return `true` if tuples should be flipped.
"""
function swap_random(
    ti::Tuple{String,<:Real},
    tj::Tuple{String,<:Real},
    beta::Real,
    rng::AbstractRNG=Random.default_rng())

    p = getflipprobability(ti,tj,beta)

    rand(rng,Float64)<p ? true : false
end

"""
    getflipprobability(ti,tj,beta)

Get the probability of flipping tuples ti and tj.
Based on formula from Handorf et al. 2008.
"""
function getflipprobability(
    ti::Tuple{String,Real},
    tj::Tuple{String,Real},
    beta::Real)

    diff = ti[2] - tj[2]

    ## Determine probability of flipping (another way of writing below)
    # ti[1]==0 || tj[1]==0 ? p=.5 : diff>0 ? p = math.e**(-diff/float(beta)) : p=1.0

    if ti[2]==0 || tj[2]==0
        p = .5
    elseif diff > 0
        p = exp(-diff/beta)
    else
        p = 1.0
    end

    return p

end

"""
    randomizecompounds(biosystem_compounds,
        compound_structs,
        n_runs,
        n_swaps,
        beta,
        sortkey,
        zero_mass_behavior)

Return a vector of run results.
Each run result is a randomized list of biosystem_compounds.

...
# Arguments
- `biosystem_compounds::IDs`: compounds to sort
- `compound_structs::Compounds`: compounds with metadata
- `n_runs::Int`: number of runs
- `n_swaps::Int=1000`: number of compounds to within biosystem_compounds to swap when randomizing
- `beta::Real=20`: Mixing coefficient. See Handorf et al. 2008 for details.
- `sortkey::Symbol=:exact_mass`: how to sort compounds. Can also sort by `:mol_weight`
- `zero_mass_behavior::String="end"`: what to do with compounds that have no mass. Can also be `"random"`
...
"""
function randomizecompounds(
    biosystem_compounds::IDs,
    compound_structs::Compounds,
    n_runs::Int,
    n_swaps::Int=1000,
    beta::Real=20,
    sortkey::Symbol=:exact_mass,
    zero_mass_behavior::String="end",
    rng::AbstractRNG=Random.default_rng())

    randomized_cpd_lists = [Vector{String}(undef,length(biosystem_compounds)) for _ in 1:n_runs] # allocate output memory
    ## Generate seeds and RNG objects based on rng
    rng_seeds = [rand(rng,1:typemax(Int)) for i in 1:n_runs]
    rng_list = [MersenneTwister(i) for i in rng_seeds]
    for r in 1:n_runs
        cpds_masses = sort_biosystem_compounds(compound_structs, biosystem_compounds, sortkey, zero_mass_behavior)
        mix_it_up!(cpds_masses,beta,n_swaps,rng_list[r])
        randomized_cpd_lists[r] = [c[1] for c in cpds_masses]
    end
    return randomized_cpd_lists
end

max_stars_repo_name: hbsmith/BioXP

max_issues_repo_name: hbsmith/BioXP

max_forks_repo_name: hbsmith/BioXP


text:
@test_throws DimensionMismatch shuffleobs((X, rand(149)))
@test_throws DimensionMismatch shuffleobs((X, rand(149)), obsdim=:last)

@testset "typestability" begin
    for var in vars
        @test_throws MethodError shuffleobs(var, ObsDim.Undefined())
        @test typeof(@inferred(shuffleobs(var))) <: SubArray
        @test typeof(@inferred(shuffleobs(var, ObsDim.Last()))) <: SubArray
        @test typeof(@inferred(shuffleobs(var, ObsDim.First()))) <: SubArray
        #@test_throws ErrorException @inferred(shuffleobs(var, obsdim=:last))
        @test_throws ErrorException @inferred(shuffleobs(var, obsdim=1))
    end
    for tup in tuples
        @test typeof(@inferred(shuffleobs(tup))) <: Tuple
        @test typeof(@inferred(shuffleobs(tup, ObsDim.Last()))) <: Tuple
        @test_throws ErrorException @inferred(shuffleobs(tup, obsdim=:last))
    end
end

@testset "Array and SubArray" begin
    for var in vars
        @test size(shuffleobs(var)) == size(var)
        @test size(shuffleobs(var, obsdim=1)) == size(var)
    end
    # tests if all obs are still present and none duplicated
    @test vec(sum(shuffleobs(X1),dims=2)) == fill(11325,10)
    @test vec(sum(shuffleobs(X1',obsdim=1),dims=1)) == fill(11325,10)
    @test sum(shuffleobs(Y1)) == 11325
    @test sum(shuffleobs(Y1, obsdim=:first)) == 11325
end

println("<HEARTBEAT>")

@testset "Tuple of Array and SubArray" begin
    for var in ((X,yv), (Xv,y), tuples...)
        @test_throws MethodError shuffleobs(var, ObsDim.Undefined())
        @test_throws MethodError shuffleobs(var...)
        @test typeof(shuffleobs(var)) <: Tuple
        @test all(map(x->(typeof(x)<:SubArray), shuffleobs(var)))
        @test all(map(x->(nobs(x)===150), shuffleobs(var)))
    end
    # tests if all obs are still present and none duplicated
    # also tests that both paramter are shuffled identically
    x1, y1, z1 = shuffleobs((X1,Y1,X1))
    @test vec(sum(x1,dims=2)) == fill(11325,10)
    @test vec(sum(z1,dims=2)) == fill(11325,10)
    @test sum(y1) == 11325
    @test all(x1' .== y1)
    @test all(z1' .== y1)
    x1, y1 = shuffleobs((X1',Y1), obsdim=1)
    @test vec(sum(x1,dims=1)) == fill(11325,10)
    @test sum(y1) == 11325
    @test all(x1 .== y1)
end

@testset "SparseArray" begin
    for var in (Xs, ys)
        @test typeof(shuffleobs(var)) <: DataSubset
        @test nobs(shuffleobs(var)) == nobs(var)
        @test nobs(shuffleobs(var, obsdim=:first)) == nobs(var, obsdim=:first)
    end
    # tests if all obs are still present and none duplicated
    @test vec(sum(getobs(shuffleobs(sparse(X1))),dims=2)) == fill(11325,10)
    @test vec(sum(getobs(shuffleobs(sparse(X1'),obsdim=1)),dims=1)) == fill(11325,10)
    @test sum(getobs(shuffleobs(sparse(Y1)))) == 11325
    @test sum(getobs(shuffleobs(sparse(Y1), obsdim=:first))) == 11325
end

@testset "Tuple of SparseArray" begin
    for var in ((Xs,ys), (X,ys), (Xs,y), (Xs,Xs), (XX,X,ys))
        @test_throws MethodError shuffleobs(var, ObsDim.Undefined())
        @test_throws MethodError shuffleobs(var...)
        @test typeof(shuffleobs(var)) <: Tuple
        @test nobs(shuffleobs(var)) == nobs(var)
    end
    # tests if all obs are still present and none duplicated
    # also tests that both paramter are shuffled identically
    x1, y1 = getobs(shuffleobs((sparse(X1),sparse(Y1))))
    @test vec(sum(x1,dims=2)) == fill(11325,10)
    @test sum(y1) == 11325
    @test all(x1' .== y1)
    x1, y1 = getobs(shuffleobs((sparse(X1'),sparse(Y1)), obsdim=1))
    @test vec(sum(x1,dims=1)) == fill(11325,10)
    @test sum(y1) == 11325
    @test all(x1 .== y1)
end

@testset "ObsView" begin
    # tests if all obs are still present and none duplicated
    ov = @inferred shuffleobs(obsview(X1))
    @test ov isa ObsView
    x1 = getobs(ov)
    @test sum(x1) == fill(11325,10)
    # also tests that both paramter are shuffled identically
    ov = @inferred shuffleobs(obsview((X1,X1)))
    @test ov isa ObsView
    x1 = getobs(ov)
    for i = 1:length(x1)
        @test x1[i][1] == x1[i][2]
    end
    for i = 1:2
        @test sum(getindex.(x1,i)) == fill(11325,10)
    end
end

@testset "BatchView" begin
    # tests if all obs are still present and none duplicated
    bv1 = batchview(X1, 30)
    bv = @inferred shuffleobs(bv1)
    @test length(bv) == length(bv1)
    @test bv isa BatchView{<:SubArray,<:SubArray}
    @test vec(sum(sum(bv),dims=2)) == fill(11325,10)
end


text:
### Exports ###

export PatchClamp

### Types ###

# TODO: Change dir to be a type parameter
@with_kw struct PatchClamp{A<:AbstractArray} <: AbstractNode
  dir::Symbol
  buffer::A
  conns::Vector{SynapticConnection} = SynapticConnection[]
end
PatchClamp(dir::Symbol, sz::Vararg{Int,N} where N) =
  PatchClamp(dir, Float32, sz...)
PatchClamp(dir::Symbol, typ::Type{T}, sz::Vararg{Int,N} where N) where T =
  PatchClamp(dir=dir, buffer=zeros(T, sz...))

### Methods ###

Base.size(pc::PatchClamp) = size(pc.buffer)
Base.getindex(pc::PatchClamp, idx...) =
  getindex(pc.buffer, idx...)

reinit!(pc::PatchClamp) =
  fill!(pc.buffer, zero(eltype(pc.buffer)))

function addedge!(pc::PatchClamp, dstcont, dst, op, conf)
  @assert first(size(dstcont)) == first(size(pc.buffer)) "Input and output sizes must match"
  @assert pc.dir === :output "PatchClamp is not configured for output"
  if op == :input
    push!(pc.conns, SynapticConnection(dst, op, nothing))
  else
    error("PatchClamp cannot handle op: $op")
  end
end
deledge!(pc::PatchClamp, dst, op) =
  filter!(conn->conn.uuid!==dst, pc.conns)
deledges!(pc::PatchClamp) = empty!(pc.conns)
connections(pc::PatchClamp) = pc.conns

function eforward!(cont::CPUContainer{PC}, args) where PC<:PatchClamp
  pc = root(cont)

  if pc.dir === :output
    pc.buffer .= zero(eltype(pc.buffer))
    for conn in connections(pc)
      cont = args[findfirst(arg->arg[2]==conn.uuid, args)][3]
      inputs = cont[:]
      pc.buffer .+= inputs
    end
  end
end


text:
import Base.convert

"""
Converter: Packed_MvNormal -> MvNormal
"""
function convert(::Type{Distributions.MvNormal}, pv::Dict{String, Any})
    len = length(Float64.(pv["mean"]))
    mat = reshape(Float64.(pv["cov"]), len, len)
    return Distributions.MvNormal(Float64.(pv["mean"]), mat)
end

"""
Converter: MvNormal -> Packed_MvNormal
"""
function convert(::Type{Dict{String, Any}}, mvNormal::Distributions.MvNormal)
    v = mvNormal.Σ.mat[:]
    return JSON.parse(JSON.json(Packed_MvNormal(mvNormal.μ, v, "MvNormal")))
end

"""
Converter: Packed_Normal -> Normal
"""
function convert(::Type{Distributions.Normal}, pv::Dict{String, Any})
    return Distributions.Normal(Float64(pv["mean"]), Float64(pv["std"]))
end

"""
Converter: Normal -> Packed_Normal
"""
function convert(::Type{Dict{String, Any}}, normal::Distributions.Normal)
    return JSON.parse(JSON.json(Packed_Normal(normal.μ, normal.σ, "Normal")))
end

"""
Converter: Packed_AliasingScalarSampler -> AliasingScalarSampler
"""
function convert(::Type{IncrementalInference.AliasingScalarSampler}, pv::Dict{String, Any})
    sampler = IncrementalInference.AliasingScalarSampler(Float64.(pv["samples"]), Float64.(pv["weights"]); SNRfloor=pv["quantile"])
    return sampler
end

"""
Converter: AliasingScalarSampler -> Packed_AliasingScalarSampler
"""
function convert(::Type{Packed_AliasingScalarSampler}, sampler::IncrementalInference.AliasingScalarSampler)
    packed = Packed_AliasingScalarSampler(sampler.domain, sampler.weights.values, 0.0, "AliasingScalarSampler")
    return JSON.parse(JSON.json(packed))
end


text:
@testset "fit moments numerical" begin
    DN = LogitNormal(1.3)
    x = rand(DN, 1_000_000);
    m = mean(DN)
    @test abs(m - mean(x))/mean(x) <= 1e-3
    s2 = var(DN)
    @test abs(s2 - var(x))/var(x) <= 1e-2 #1e-4 too strong for random numbers
    s2m = var(DN; mean=m) # specify so that not need to recompute
    @test abs(s2m - var(x))/var(x) <= 1e-2 #1e-4 too strong for random numbers
    sd = std(DN) # specify so that not need to recompute
    @test abs(sd - std(x))/std(x) <= 1e-2 #1e-4 too strong for random numbers
    # sd2 = std(DN, mean=m) # specify so that not need to recompute
    # @test abs(sd2 - std(x))/std(x) <= 1e-2 #1e-4 too strong for random numbers
end;

function test_lognormal_mode(g)
    mo = mode(g)
    @test pdf(g, mo-1e-6) < pdf(g,mo) && pdf(g, mo+1e-6) < pdf(g,mo)
end

###### numerical estimation of moments
@testset "Logitnormal seek mode" begin
    #plot(g); vline!([mode(g)])
    # median
    g = LogitNormal(); test_lognormal_mode(g)
    g = LogitNormal(0,1.4); test_lognormal_mode(g)
    # larger one of two modes, infer the lower one
    g = LogitNormal(0,1.6); test_lognormal_mode(g)
    # two modes, larger is on the right
    g = LogitNormal(0.1,1.6); test_lognormal_mode(g)
    # two modes, larger is on the left
    g = LogitNormal(-0.1,1.6); test_lognormal_mode(g)
end

@testset "fit two quantiles" begin
    qpl = @qp_m(0.3)
    qpu = @qp_u(0.9)
    d = fit(LogitNormal, qpl, qpu);
    @test quantile.(d, [qpl.p, qpu.p]) ≈ [qpl.q, qpu.q]
    d = fit(LogitNormal, qpu, qpl) # sort
    @test quantile.(d, [qpl.p, qpu.p]) ≈ [qpl.q, qpu.q]
end;

@testset "fit to quantilepoint and mean" begin
    # not implemented
end;

@testset "fit to quantilepoint and mode" begin
    # see LogitNormals.jl
    # d = LogitNormal(1,1)
    # m = log(mode(d))
    # qp = @qp(quantile(d,0.95),0.95)
    # dfit = fit_mode_quantile(LogitNormal, mode(d), qp)
    # @test dfit ≈ d
    # dfit = fit(LogitNormal, mode(d), qp, Val(:mode))
    # @test dfit ≈ d
    # # with lower quantile
    # qp = @qp(quantile(d,0.025),0.025)
    # dfit = fit_mode_quantile(LogitNormal, mode(d), qp)
    # @test mode(dfit) ≈ mode(d) && quantile(dfit, qp.p) ≈ qp.q
end;

@testset "fit to quantilepoint and median" begin
    d = LogitNormal(1,1)
    qp = @qp(quantile(d,0.95),0.95)
    dfit = fit(LogitNormal, median(d), qp, Val(:median))
    @test dfit ≈ d
end;

@testset "fit by mode and upper quantile" begin
    #plot(g); vline!([mode(g)])
    # median
    #g = fit_mode_quantile(LogitNormal, 0.68, @qp(0.71,0.99))
    g = fit(LogitNormal, 0.68, @qp_u(0.9), Val(:mode))
    @test mode(g) ≈ 0.68
    @test quantile(g, 0.95) ≈ 0.9
    #plot(g)
    # call with differt subtype of Real, here Rational
    g = fit(LogitNormal, 2//3, @qp_u(0.9), Val(:mode))
    @test mode(g) ≈ 2//3
    @test quantile(g, 0.95) ≈ 0.9
end

function is_logit_slope_monotone(d, upper=0.5, lower=0.0, decreasing = false)
  x = range(lower,upper;length=41)[2:40]	# plotting grid
  dx = pdf.(d, x)	#density function
  if decreasing
    all(diff(dx) .<= 0)
  else
    all(diff(dx) .>= 0)
  end
end

@testset "fit by single mode and flat" begin
    d9 = d = fit_mode_flat(LogitNormal, 0.9)
    @test isapprox( mode(d), 0.9, atol=1e-4)
    @test is_logit_slope_monotone(d, 0.9)
    @test is_logit_slope_monotone(d, 0.0, 0.9, true)
    d1 = d = fit_mode_flat(LogitNormal, 0.1)
    @test isapprox( mode(d), 0.1, atol=1e-4)
    @test is_logit_slope_monotone(d, 0.1)
    @test is_logit_slope_monotone(d, 0.0, 0.1, true)
    @test d1.σ == d9.σ
    @test d1.μ == -d9.μ
    d5 = d = fit_mode_flat(LogitNormal, 0.5)
    @test isapprox( mode(d), 0.5, atol=1e-4)
    @test is_logit_slope_monotone(d, 0.5)
    @test is_logit_slope_monotone(d, 0.0, 0.5, true)
    #plot(d1)
end

@testset "shifloNormal" begin
    d = shifloNormal(1,3)
    @test isapprox(mode(d), 2.0, atol=0.001)
    @test minimum(d) == 1.0
    @test maximum(d) == 3.0
    @test scale(d) == 3-1
    @test location(d) == 1.0
    #plot(d)
end;


text:
# ------------ Utility selection, grads and params ------------ #

abstract type UtilitySelection <: Selection end

function fill_array!(val::T, arr::Vector{K}, f_ind::Int) where {K, T <: UtilitySelection}
    sorted_toplevel_keys = sort(collect(keys(val.utility)))
    sorted_tree_keys  = sort(collect(keys(val.tree)))
    idx = f_ind
    for k in sorted_toplevel_keys
        v = val.utility[k]
        n = fill_array!(v, arr, idx)
        idx += n
    end
    for k in sorted_tree_keys
        n = fill_array!(get_sub(val, k), arr, idx)
        idx += n
    end
    idx - f_ind
end

function from_array(schema::T, arr::Vector{K}, f_ind::Int) where {K, T <: UtilitySelection}
    # Subtypes of UtilitySelection must support an empty constructor.
    sel = T()
    sorted_toplevel_keys = sort(collect(keys(schema.utility)))
    sorted_tree_keys  = sort(collect(keys(schema.tree)))
    idx = f_ind
    for k in sorted_toplevel_keys
        (n, v) = from_array(schema.utility[k], arr, idx)
        idx += n
        sel.utility[k] = v
    end
    for k in sorted_tree_keys
        (n, v) = from_array(get_sub(schema, k), arr, idx)
        idx += n
        sel.tree[k] = v
    end
    (idx - f_ind, sel)
end

# ------------ Gradients ------------ #

struct Gradients <: UtilitySelection
    tree::Dict{Address, Gradients}
    utility::Dict{Address, Any}
    Gradients() = new(Dict{Address, Gradients}(), Dict{Address, Any}())
    Gradients(tree, utility) = new(tree, utility)
end

has_top(ps::Gradients, addr::T) where T <: Address = haskey(ps.utility, addr)
function has_top(ps::Gradients, addr::T) where T <: Tuple
    isempty(addr) && return false
    length(addr) == 1 && return has_top(ps, addr[1])
    return has_top(ps.tree[addr[1]], addr[2 : end])
end

get_top(ps::Gradients, addr::T) where T <: Address = getindex(ps.utility, addr)
function get_top(ps::Gradients, addr::T) where T <: Tuple
    isempty(addr) && error("get_top: index error - tuple address is empty.")
    length(addr) == 1 && return get_top(ps, addr[1])
    return get_top(ps.tree[addr[1]], addr[2 : end])
end

has_sub(ps::Gradients, addr::T) where T <: Address = haskey(ps.tree, addr)
function has_sub(ps::Gradients, addr::T) where T <: Tuple
    isempty(addr) && return false
    length(addr) == 1 && return has_sub(ps, addr[1])
    return has_sub(ps.tree[addr[1]], addr[2 : end])
end

get_sub(ps::Gradients, addr::T) where T <: Address = getindex(ps.tree, addr)
function get_sub(ps::Gradients, addr::T) where T <: Tuple
    isempty(addr) && error("get_sub: index error - tuple address is empty.")
    length(addr) == 1 && return get_sub(ps, addr[1])
    return get_sub(ps.tree[addr[1]], addr[2 : end])
end

getindex(ps::Gradients, addr::T) where T <: Address = getindex(ps.utility, addr)
getindex(ps::Gradients, addr::Tuple{}) = error("Gradients (getindex): empty tuple as index.")
getindex(ps::Gradients, addr::Tuple{T}) where T <: Address = getindex(ps, addr[1])
function getindex(ps::Gradients, addr::T) where T <: Tuple
    has_sub(ps, addr[1]) && return getindex(get_sub(ps, addr[1]), addr[2 : end])
    error("Gradients (getindex): invalid index at $addr.")
end
function getindex(ps::Gradients, addrs...)
    has_top(ps, addrs) && return get_top(ps, addrs)
    has_sub(ps, addrs) && return get_sub(ps, addrs)
    error("Gradients (getindex): invalid index at $addr.")
end

# ------------ Builders ------------ #

function push!(ps::Gradients, addr::T, val) where T <: Address
    has_top(ps, addr) && begin
        ps.utility[addr] += val
        return
    end
    ps.utility[addr] = val
end

function push!(ps::Gradients, addr::T, val) where T <: Tuple
    isempty(addr) && return
    length(addr) == 1 && push!(ps, addr[1], val)
    hd = addr[1]
    tl = addr[2 : end]
    if has_sub(ps, hd)
        push!(get_sub(ps, hd), tl, val)
    else
        sub = Gradients()
        push!(sub, tl, val)
        ps.tree[hd] = sub
    end
end

# ------------ Adjoint ------------ #

Zygote.@adjoint Gradients(tree, utility) = Gradients(tree, utility), s_grad -> (nothing, nothing)

# ------------ Combining two sets of gradients ------------ #

function merge(sel1::Gradients,
               sel2::Gradients)
    utility = merge(sel1.utility, sel2.utility)
    tree = Dict{Address, Gradients}()
    for k in setdiff(keys(sel2.tree), keys(sel1.tree))
        tree[k] = sel2.tree[k]
    end
    for k in setdiff(keys(sel1.tree), keys(sel2.tree))
        tree[k] = sel1.tree[k]
    end
    for k in intersect(keys(sel1.tree), keys(sel2.tree))
        tree[k] = merge(sel1.tree[k], sel2.tree[k])
    end
    return Gradients(tree, utility)
end

+(a_grads::Gradients, b_grads::Gradients) = merge(a_grads, b_grads)

Zygote.@adjoint merge(a_grads, b_grads) = merge(a_grads, b_grads), s_grad -> (nothing, nothing)

# ------------ Parameters, empty and learnable ------------ #

abstract type Parameters <: UtilitySelection end

struct EmptyParameters <: Parameters end

Parameters() = EmptyParameters()
has_top(np::EmptyParameters, addr) = false
get_top(np::EmptyParameters, addr) = error("(get_param) called on instance of EmptyParameters. No parameters!")
has_sub(np::EmptyParameters, addr) = false
get_sub(np::EmptyParameters, addr) = np

struct LearnableParameters <: Parameters
    tree::Dict{Address, LearnableParameters}
    utility::Dict{Address, Any}
    LearnableParameters() = new(Dict{Address, LearnableParameters}(), Dict{Address, Any}())
end

has_sub(ps::LearnableParameters, addr::T) where T <: Address = haskey(ps.tree, addr)
function has_sub(ps::LearnableParameters, addr::T) where T <: Tuple
    isempty(addr) && return false
    length(addr) == 1 && return has_sub(ps, addr[1])
    has_sub(ps, addr[1]) && has_sub(get_sub(ps, addr[1]), addr[2 : end])
end

get_sub(ps::LearnableParameters, addr::T) where T <: Address = getindex(ps.tree, addr)
function get_sub(ps::LearnableParameters, addr::T) where T <: Tuple
    isempty(addr) && return Parameters()
    length(addr) == 1 && return get_sub(ps, addr[1])
    return get_sub(ps.tree[addr[1]], addr[2 : end])
end

set_sub!(ps::LearnableParameters, addr::T, sub::LearnableParameters) where T <: Address = ps.tree[addr] = sub
function set_sub!(ps::LearnableParameters, addr::T, sub::LearnableParameters) where T <: Tuple
    isempty(addr) && return
    length(addr) == 1 && set_sub!(ps, addr[1], sub)
    has_sub(ps, addr[1]) && set_sub!(get_sub(ps, addr[1]), addr[2 : end], sub)
end

has_top(ps::LearnableParameters, addr::T) where T <: Address = haskey(ps.utility, addr)
function has_top(ps::LearnableParameters, addr::T) where T <: Tuple
    isempty(addr) && return false
    length(addr) == 1 && return has_top(ps, addr[1])
    has_sub(ps, addr[1]) && return has_top(get_sub(ps, addr[1]), addr[2 : end])
end

get_top(ps::LearnableParameters, addr::T) where T <: Address = getindex(ps.utility, addr)
function get_top(ps::LearnableParameters, addr::T) where T <: Tuple
    isempty(addr) && error("LearnableParameters (get_top): invalid index at $addr.")
    length(addr) == 1 && return get_top(ps, addr[1])
    get_top(ps.tree[addr[1]], addr[2 : end])
end

getindex(ps::LearnableParameters, addrs...) = get_top(ps, addrs)

# ------------ Builders ------------ #

function push!(ps::LearnableParameters, addr::T, val) where T <: Address
    ps.utility[addr] = val
end

function push!(ps::LearnableParameters, addr::T, val) where T <: Tuple
    isempty(addr) && return
    length(addr) == 1 && push!(ps, addr[1], val)
    hd = addr[1]
    tl = addr[2 : end]
    if has_sub(ps, hd)
        push!(get_sub(ps, hd), tl, val)
    else
        sub = LearnableParameters()
        push!(sub, tl, val)
        ps.tree[hd] = sub
    end
end

function parameters(arr::Array{Pair{T, K}}) where {T <: Tuple, K}
    top = LearnableParameters()
    map(arr) do (k, v)
        push!(top, k, v)
    end
    return top
end

function parameters(p::Pair{T, K}) where {T <: Symbol, K <: Parameters}
    top = LearnableParameters()
    set_sub!(top, p[1], p[2])
    return top
end

# ------------ Adjoint ------------ #

Zygote.@adjoint LearnableParameters(tree, utility) = LearnableParameters(tree, utility), s_grad -> (nothing, nothing)

# ------------ Merging two sets of parameters ------------ #

function merge(sel1::LearnableParameters,
               sel2::LearnableParameters)
    utility = merge(sel1.utility, sel2.utility)
    tree = Dict{Address, Gradients}()
    for k in setdiff(keys(sel2.tree), keys(sel1.tree))
        tree[k] = sel2.tree[k]
    end
    for k in setdiff(keys(sel1.tree), keys(sel2.tree))
        tree[k] = sel1.tree[k]
    end
    for k in intersect(keys(sel1.tree), keys(sel2.tree))
        tree[k] = merge(sel1.tree[k], sel2.tree[k])
    end
    return LearnableParameters(tree, utility)
end

+(a::LearnableParameters, b::LearnableParameters) = merge(a, b)

# ------------ update_parameters links into Flux optimiser APIs ------------ #

function update_parameters(opt, a::LearnableParameters, b::Gradients)
    p_arr = array(a, Float64)
    gs_arr = array(b, Float64)
    update!(opt, p_arr, -gs_arr)
    return selection(a, p_arr)
end

# ------------ Pretty printing utility selections ------------ #

function collect!(par::T, addrs::Vector, chd::Dict, chs::K) where {T <: Tuple, K <: UtilitySelection}
    for (k, v) in chs.utility
        push!(addrs, (par..., k))
        chd[(par..., k)] = v
    end
    for (k, v) in chs.tree
        collect!((par..., k), addrs, chd, v)
    end
end

function collect!(addrs::Vector, chd::Dict, chs::K) where K <: UtilitySelection
    for (k, v) in chs.utility
        push!(addrs, (k, ))
        chd[(k, )] = v
    end
    for (k, v) in chs.tree
        collect!((k, ), addrs, chd, v)
    end
end

import Base.collect
function collect(chs::K) where K <: UtilitySelection
    addrs = []
    chd = Dict()
    collect!(addrs, chd, chs)
    return addrs, chd
end

function Base.display(chs::Gradients; show_values = true)
    println("  __________________________________\n")
    println("             Gradients\n")
    addrs, chd = collect(chs)
    if show_values
        for a in addrs
            println(" $(a) : $(chd[a])")
        end
    else
        for a in addrs
            println(" $(a)")
        end
    end
    println("  __________________________________\n")
end

function Base.display(chs::LearnableParameters; show_values = true)
    println("  __________________________________\n")
    println("             Parameters\n")
    addrs, chd = collect(chs)
    if show_values
        for a in addrs
            println(" $(a) : $(chd[a])")
        end
    else
        for a in addrs
            println(" $(a)")
        end
    end
    println("  __________________________________\n")
end


text:
# This file was generated, do not modify it. # hide
ŷ = predict_mode(mach, Xtrain)
mcr = misclassification_rate(ŷ, ytrain)
println(rpad("MNC mcr:", 10), round(mcr, sigdigits=3))


text:
# three class for each phoneme
function build_phoneme_vocab(entries)
    c = 1
    p2i, i2p = Dict(), Dict()
    for sample in entries
        phns = sample["phonemes"]
        for (_,_,phn) in phns
            phnb = phn*"-b" # begin
            if !haskey(p2i, phnb)
                p2i[phnb] = c
                i2p[c] = phnb
                c += 1
            end

            phnm = phn*"-m" # middle
            if !haskey(p2i, phnm)
                p2i[phnm] = c
                i2p[c] = phnm
                c += 1
            end

            phne = phn*"-e" # end
            if !haskey(p2i, phne)
                p2i[phne] = c
                i2p[c] = phne
                c += 1
            end
        end
    end
    return p2i, i2p
end

# separate vocabulary builder
function build_vocab(entries; kk="phonemes")
    c = 1
    w2i, i2w = Dict(), Dict()
    for sample in entries
        words = sample[kk]
        for (_,_,w) in words
            if !haskey(w2i, w)
                w2i[w] = c
                i2w[c] = w
                c += 1
            end
        end
    end
    return w2i, i2w
end

function make_input(samples)
    return cat(4, map(si->si[1], samples)...)
end

function make_output(samples, p2i)
    return map(s->p2i[s[2]], samples)
end

# feature extractor
function mfsc(arr; sr=16000.0, preprocess=true)
    if preprocess
        nch = 1
        ns = div(length(arr), nch)
        arr = reshape(arr, nch, ns)' / (1<<15)
    end

    feacalc(
        arr; defaults=:rasta, sadtype=:none, normtype=:mvn, augtype=:ddelta,
        dcttype=-1, sr=sr)
end

# read features arrays from hdf5 file
function make_data(
    filepath, entries, split="train";
    remove_sa1=true, nframes=15, sr=16000.0, triphones=true, extra=[])
    samples = filter(ei->ei["split"]==split, entries)
    if remove_sa1
        samples = filter(si->si["sample"] != "sa1", samples)
    end

    features = load(filepath, "features")
    data = []
    i = 0
    for sample in samples
        feature = features[sample["longid"]]'
        # i+=1; info("i=$i")
        ss = []
        ll = div(nframes-1,2)
        isr(x) = Int(round(x/(sr/100)))
        # @show sample["phonemes"]
        pp = map(p->(max(1,isr(p[1])), min(size(feature,2),isr(p[2])), p[3]), sample["phonemes"])
        # @show pp
        for k = max(1+ll,pp[1][1]):size(feature,2)-ll
            while !(k >= pp[1][1] && k <= pp[1][2])
                # @show k, pp[1][1], pp[1][2]
                shift!(pp)
                length(pp) == 0 && break
            end

            length(pp) == 0 && break
            ff = feature[:,k-ll:k+ll]
            ff = reshape(ff, 40, 1, nframes*3, 1)
            push!(ss, (ff, pp[1][3]))
        end

        # for (ti,tf,phn) in sample["phonemes"]
        #     ll = div(nframes-1,2)
        #     ff = zeros(size(feature,1), size(feature,2)+2*ll+1)
        #     ff[:,1+ll:end-1-ll] = feature

        #     ti0 = Int(round(ti/(sr/100)))
        #     tf0 = Int(round(tf/(sr/100)))
        #     tm = div(tf0-ti0,2)+ti0
        #     ll = div((nframes-1),2)

        #     # @show ti0,tf0,tm,ll,size(ff)
        #     push!(ss, (reshape(ff[:,tm:tm+2ll],40,1,nframes*3,1), phn))
        # end
        push!(data, ss)
        push!(extra, sample["longid"])
    end

    return data
end

function get_features(filepath,entries)
    filepath = abspath(filepath)
    features = h5open(filepath, "w") do f
        data = Dict()
        for entry in entries
            data[entry["longid"]] = read(f, entry["longid"])
        end
        data
    end
    return features
end

# transition probabilities, a frequentist approach
function get_probabilities(jsondata, p2i, i2p; sr=16000.0)
    posteriors = zeros(Float32, length(p2i), length(p2i))
    priors = zeros(Float32, length(p2i))

    entries = filter(
        j->j["split"]=="train" && j["sample"] != "sa1", jsondata)
    for entry in entries
        isr(x) = Int(round(x/(sr/100)))
        segments = entry["phonemes"]
        segments = map(s->(isr(s[1]),isr(s[2]),s[3]), segments)
        for k = 1:length(segments)-1
            curr = segments[k]
            next = segments[k+1]

            dif = curr[2]-curr[1]
            priors[p2i[curr[3]]] += dif+1
            posteriors[p2i[curr[3]],p2i[next[3]]] += 1
            posteriors[p2i[curr[3]],p2i[curr[3]]] += dif
        end
        curr = segments[end]
        dif = segments[end][2]-segments[end][1]
        posteriors[p2i[curr[3]],p2i[curr[3]]] += dif
        priors[p2i[segments[end][3]]] += dif+1
    end

    # scale
    priors = priors / sum(priors)
    posteriors = posteriors ./ sum(posteriors,2)

    return posteriors, priors
end

function generate(ysoft, wlen)
    inds = mapslices(indmax, Array(ysoft), 1)
    candidate = []
    for k = 1:wlen:length(inds)
        arr = inds[k:min(length(inds),k+wlen-1)]
        mydict = Dict()
        for i = 1:length(arr)
            mydict[arr[i]] = get(mydict, arr[i], 0) + 1
        end
        myarr = [(k,v) for (k,v) in mydict]
        sort!(myarr, by=x->x[2], rev=true)
        push!(candidate, myarr[1][1])
    end
    res = []
    push!(res, pop!(candidate))
    while !isempty(candidate)
        el = shift!(candidate)
        el != res[end] && push!(res, el)
    end
    return res
end

function make_ygold(fake)
    ygold = [fake[1]]
    for k = 2:length(fake)
        if ygold[end] != fake[k]
            push!(ygold, fake[k])
        end
    end
    return ygold
end

max_stars_repo_name: ilkerkesen/phoneme-convnet

max_issues_repo_name: ilkerkesen/phoneme-convnet

max_forks_repo_name: ilkerkesen/phoneme-convnet


text:
function execute_sts_adjusted_daily_api_call(user_model::AVKUserModel, stock_symbol::String; data_type::Symbol = :json, outputsize::Symbol = :compact, logger::Union{Nothing,AbstractLogger} = nothing)::(Union{T, Nothing} where T<:Any)
    # some error checks -
    # is user_models valid?
    check_result = _check_user_model(user_model)
    if (check_result != nothing && typeof(check_result.value) == AVKError)
        return check_result
    end

    # do we have the API key?
    check_result = _check_missing_api_key(user_model)
    if (check_result != nothing && typeof(check_result.value) == AVKError)
        return check_result
    end

    # do we have a stock_symbol -
    check_result = _check_missing_symbol(stock_symbol)
    if (check_result != nothing && typeof(check_result.value) == AVKError)
        return check_result
    end

    # get key's and id -
    api_key = user_model.alphavantage_api_key

    #call to alpha_vantage_api to get data
    url = "$(alphavantage_api_url_string)?function=TIME_SERIES_DAILY_ADJUSTED&symbol=$(stock_symbol)&apikey=$(api_key)&datatype=$(string(data_type))&outputsize=$(string(outputsize))"
    api_call_result = _http_get_call_with_url(url)
    if (typeof(api_call_result.value) == AVKError)
        return api_call_result
    end#keep the same as daily call

    #check that result value is a string
    #this is also new but doesn't change between time periods
    if (typeof(api_call_result.value) isa String == false)
        println("Call result is not valid type")#error from not correct data output
    else
        api_call_raw_data = api_call_result.value
    end

    #call to logger
    if logger != nothing
        debug_log_api_call(logger, user_model, url)
    end
    #parse json if data called as a .json
    if (data_type == :json)
        #process json
        data_series_key = "Time Series (Daily)"
        return _process_raw_json_api_data_sts_adjusted(api_call_raw_data, data_series_key)

    elseif (data_type == :csv)
        #return process .csv
        return _process_raw_csv_api_data(api_call_raw_data, data_series_key)
    else
        #tell user they requested an unsupported type of data
        error_message = "$(data_type) isn't supported by AlphaVantage. Supported values are {:json, :csv}"
        return Result{AVKError}(AVKError(error_message))
    end
end

max_stars_repo_name: teamalphavantage/AlphaVantageAPI

max_issues_repo_name: teamalphavantage/AlphaVantageAPI

max_forks_repo_name: teamalphavantage/AlphaVantageAPI


text:
module TestDistance
using ParallelKMeans: colwise!
using Test

@testset "naive singlethread colwise" begin
    X = [1.0 3.0 4.0; 2.0 5.0 6.0]
    y = [1.0, 2.0]
    r = Vector{Float64}(undef, 3)

    colwise!(r, X, y, 1)
    @test all(r .≈ [0.0, 13.0, 25.0])
end

@testset "multithread colwise" begin
    X = [1.0 3.0 4.0; 2.0 5.0 6.0]
    y = [1.0, 2.0]
    r = Vector{Float64}(undef, 3)

    colwise!(r, X, y, 2)
    @test all(r .≈ [0.0, 13.0, 25.0])
end

end # module


text:
module RIOPA

import MPI

function __init__()
    # - Initialize here instead of main so that the MPI context can be available
    # for tests.
    # - Conditional allows for case when external users (or tests) have already
    # initialized an MPI context.
    if !MPI.Initialized()
        MPI.Init()
    end
end

include(joinpath("helper", "Ratios.jl"))

include(joinpath("core", "Config.jl"))
include(joinpath("core", "Args.jl"))
include(joinpath("core", "DataStream.jl"))
include(joinpath("core", "DataSet.jl"))
include(joinpath("core", "datagen", "DataGen.jl"))
include(joinpath("core", "io", "IO.jl"))
include(joinpath("core", "Ctrl.jl"))

include(joinpath("core", "io", "HDF5IOBackend.jl"))

# include(joinpath("hello", "adios2.jl"))
include(joinpath("hello", "hdf5.jl"))
include(joinpath("hello", "hello.jl"))

import MLStyle: @match

function main(args)::Int32
    inputs = RIOPA.parse_inputs(args)

    config_filename = inputs["config"]

    command = inputs["%COMMAND%"]
    @match command begin
        "hello" => RIOPA.hello(config_filename)
        "generate-config" => RIOPA.generate_config(config_filename)
        nothing => RIOPA.Ctrl.run(config_filename)
    end

    return 0
end

end


text:
addprocs(20);
@everywhere include("loadMod.jl");
@everywhere const GUROBI_ENV = Gurobi.Env();

caseList = [13,33,123];
NN = 1000;
Δt = 0.25;
N = 5;
iterMax = 20;
pathListDRaw = load("pathHist_1000.jld");
TList = [24,36,48,72,96];

for ci in 1:length(caseList)
    pathDictA = pathListDRaw["pathDict"][ci];
    stochNomOut = Dict();
    for T in TList
        τ = Int64(1/6*T);
        for j in procs()
            remotecall_fetch(readInData,j,ci,caseList,T);
            #remotecall_fetch(readInData_old,j,T,ωSet0,10000,0);
        end
        pathDict = pathDictA[T];
        pathDict1 = Dict();
        noP = 1;
        foundBool = false;
        while !(foundBool)
            if length(pathDict[noP]) == 1
                pathDict1[1] = pathDict[noP];
                foundBool = true;
            else
                noP += 1;
            end
        end

        # train the stochastic programming strategy
        cutDictPG = preGen(τ, T, Δt, N, iterMax);
        # cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(τ, T, Δt, N, true, false, 10, 10, cutDictPG);
        cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(τ, T, Δt, N, false, false,
            max(Int64(round(500/N)),20), max(Int64(round(500/N)),20),cutDictPG);
        for j in procs()
            remotecall_fetch(cutIni,j,cutDict);
        end

        solSDDP, LBSDDP, costSDDP = exeForward(τ, T, Δt, 1, false, pathDict1);
        stochNomOut[T] = [solSDDP,LBSDDP,costSDDP];

        save("stochNomResults_$(ci).jld", "stochNomOut", stochNomOut);
    end
end


text:
@testset "Quadrature" begin
    tests = ["gauss_quadrature",
            ]
    for t in tests
      include("$(t).jl")
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
module Valuation

using HTTP
using JSON
using Dates

using ..CurrenciesBase
using ..Baskets
import ..CurrenciesBase: filltype

export valuate, ExchangeRateTable, ecbrates

# Valuate Monetary or Basket using a table of exchange rates.

"""
A table of exchange rates denominated in some currency. In which currency it is
denominated is unimportant, so long as all exchange rates are denominated in the
same currency. The denomination currency need not even exist in the table.

Optionally, `ExchangeRateTable` objects may contain information about the date
for which they apply. If this isn't provided, then the current date is used.
"""
struct ExchangeRateTable <: AbstractDict{Symbol, Float64}
    date::Date
    table::Dict{Symbol, Float64}
end
ExchangeRateTable(table) = ExchangeRateTable(Date(now()), table)
ExchangeRateTable(entries::Pair{Symbol, Float64}...) =
    ExchangeRateTable(Dict(entries...))
ExchangeRateTable(date::Date, entries::Pair{Symbol, Float64}...) =
    ExchangeRateTable(date, Dict(entries...))

Base.iterate(ert::ExchangeRateTable, s...) = iterate(ert.table, s...)
Base.length(ert::ExchangeRateTable) = length(ert.table)
Base.haskey(ert::ExchangeRateTable, k) = haskey(ert.table, k)
Base.getindex(ert::ExchangeRateTable, k) = ert.table[k]

const ECBCache = Dict{Date, ExchangeRateTable}()

function ecbrates_fresh(datestr::AbstractString)
    error("ECB rates data not supported yet")
    # get fixer.io data
    resp = JSON.parse(String(take!(HTTP.get("https://api.fixer.io/$datestr"))))
    date = Date(resp["date"])
    table = Dict{Symbol, Float64}()
    for (k, v) in resp["rates"]
        table[Symbol(k)] = 1.0/v
    end
    table[Symbol(resp["base"])] = 1.0

    # cache and return
    ECBCache[date] = ExchangeRateTable(date, table)
end
ecbrates_fresh(date::Date) = ecbrates_fresh(string(date))
ecbrates_fresh() = ecbrates_fresh("latest")

"""
    ecbrates()           → ExchangeRateTable
    ecbrates(date::Date) → ExchangeRateTable

Get an `ExchangeRateTable` from European Central Bank data for the specified
date, or for the most recent date available. European Central Bank data is not
available for all currencies, and it is often out of date. This function
requires a connection to the Internet, and reraises whatever exception is thrown
from the `Requests` package if the connection fails for any reason.

Data is retrieved from https://fixer.io/, and then cached in memory to avoid
excessive network traffic. Because of the nature of cached data, an application
running for a long period of time may receive data that is one week or more out
of date.
"""
function ecbrates()
    date = Date(now(Dates.UTC))
    # try last seven days
    for _ in 1:7
        if haskey(ECBCache, date)
            return ECBCache[date]
        end
        date -= Dates.Day(1)
    end

    ecbrates_fresh()
end
ecbrates(date) = haskey(ECBCache, date) ? ECBCache[date] : ecbrates_fresh(date)

"""
    valuate(table, as::Symbol, amount::Monetary)   → Monetary{as}
    valuate(table, as::DataType, amount::Monetary) → as

Reduce the given `Monetary` or `Basket` to a value in a single specified
currency or of a single specified type, using the given exchange rate table. The
exchange rate table can either be an `ExchangeRateTable` or any other
`Associative` mapping `Symbol` to `Real`.

    rates = ExchangeRateTable(:USD => 1.0, :CAD => 0.75)
    valuate(rates, :CAD, 21USD)  # 28CAD
"""
function valuate(table, as::Type{Monetary{U,V,W}}, amount::Monetary{T}) where {T,U,V,W}
    rate = table[T] / table[U]
    amount / majorunit(amount) * rate * majorunit(as)
end

function valuate(table, as::Type{Monetary{U,V,W}}, amount::Basket) where {U,V,W}
    acc = 0.0
    for m in amount
        from = currency(m)
        acc += m / majorunit(m) * table[from]
    end
    acc / table[U] * majorunit(as)
end

function valuate(table, ::Type{U}, amount::AbstractMonetary) where U<:Monetary
    valuate(table, filltype(U), amount)
end

function valuate(table, T::Symbol, amount::AbstractMonetary)
    valuate(table, Monetary{T}, amount)
end

end


text:
import Base: RefValue

struct ReferenceEntry{T}
    ref::RefValue{T}
    count::Int
end
Base.pointer(ref::ReferenceEntry) = _get_pointer(ref.ref)
Base.getindex(ref::ReferenceEntry) = ref.ref[]

# Objects (id) -> Pointers
const ObjectReference = IdDict{Any, RefValue}()
# Pointers -> Reference
const ReferenceCache = Dict{Ptr{Cvoid}, ReferenceEntry}()

const TypedCompatibleArray{T} = Union{DenseArray{T}, Base.ReinterpretArray{T}}
_get_pointer(ref::RefValue{T}) where T = Ptr{Cvoid}(pointer_from_objref(T.mutable ? ref[] : ref))
_get_pointer(ref::RefValue{T}) where T <: TypedCompatibleArray = Ptr{Cvoid}(pointer(ref[]))
get_reference(id::Union{UInt64, Ptr}) = get(ReferenceCache, Ptr{Cvoid}(id), nothing)

delete_reference(id::Union{UInt64, Ptr}) = let ref = pop!(ReferenceCache, Ptr{Cvoid}(id), nothing)
    if isnothing(ref)
        false
    else
        delete!(ObjectReference, ref[])
        true
    end
end

_make_ref(object) = Ref(object)
_make_ref(arr::Base.ReinterpretArray) = _make_ref(parent(arr))

function make_reference(object)
    ref = get!(ObjectReference, object) do
        _make_ref(object)
    end
    ptr = _get_pointer(ref)
    entry = get(ReferenceCache, ptr) do
        ReferenceEntry(ref, 0)
    end
    ReferenceCache[ptr] = ReferenceEntry(entry.ref, entry.count + 1)
end
reference(object) = make_reference(object)
function reference(id::Union{UInt64, Ptr})
    ptr = Ptr{Cvoid}(id)
    entry = pop!(ReferenceCache, ptr, nothing)
    isnothing(entry) && return nothing
    ReferenceCache[ptr] = ReferenceEntry(entry.ref, entry.count + 1)
end
dereference(object::T) where T = initialized() && haskey(ObjectReference, object) ? dereference(
    _get_pointer(ObjectReference[object])
) : nothing
function dereference(id::Union{UInt64, Ptr})
    ptr = Ptr{Cvoid}(id)
    entry = pop!(ReferenceCache, ptr, nothing)
    isnothing(entry) && return nothing
    ref = ReferenceEntry(entry.ref, entry.count - 1)
    if entry.count > 1
        ReferenceCache[ptr] = ref
    else
        delete!(ObjectReference, ref[])
    end
    ref
end


text:
copy_mapslice_helper!(tgt::AbstractArray, src::LabeledArray) = copy!(tgt, peel(src))
copy_mapslice_helper!(tgt::AbstractArray, src::AbstractArray) = copy!(tgt, src)

mapslices_darr_larr(f::Function, arr::AbstractArray, dims::AbstractVector) = begin
  if length(dims) != length(unique(dims))
    throw(ArgumentError("the dims argument should be an array of distinct elements."))
  end
  dimtypes = ntuple(ndims(arr)) do d
    if d in dims
      Colon()
    else
      0
    end
  end
  result = mapslices_darr_larr_inner(f, arr, dimtypes)
  if !isa(result,AbstractArray) || isempty(result)
    result
  else
    peeloff_zero_array_if_necessary(result)
  end
end

peeloff_zero_array_if_necessary(x) = x
peeloff_zero_array_if_necessary{T}(arr::AbstractArray{T,0}) = arr[1]
peeloff_zero_array_if_necessary{K}(arr::DictArray{K,0}) = mapvalues(x->x[1], peel(arr))
peeloff_zero_array_if_necessary{T}(arr::LabeledArray{T,0}) = peeloff_zero_array_if_necessary(peel(arr))

mapslices_lambda_acc_int(acc, x) = acc + (x==Int)
mapslices_lambda_first(it) = it[1]
mapslices_lambda_second_int(it) = it[2] == Int
mapslices_lambda_colon(d) = Colon()
mapslices_lambda_one(d) = 1
mapslices_lambda_axis_tuple_func(testres, ressize, ndimstestres) = d -> begin
  if d<=ndimstestres
    pickaxis(testres, d)
  else
    DefaultAxis(ressize[d-ndimstestres])
  end
end

@generated mapslices_darr_larr_inner{V,N,T}(f::Function, arr::AbstractArray{V,N}, dims::T) = begin
  dimtypes = dims.types
  slice_ndims = foldl(mapslices_lambda_acc_int, 0, dimtypes)
  slice_indices = Int[map(mapslices_lambda_first, Iterators.filter(mapslices_lambda_second_int, enumerate(dimtypes)))...]
  slice_exp = if slice_ndims == N
    :(arr[coords...])
  else
    :(view(arr, coords...))
  end
  quote
    if isempty(arr)
      if $slice_ndims == 0
        return Nullable{Any}()
      else
        return similar(arr, Nullable{Any}, size(arr)[$slice_indices])
      end
    end
    sizearr = size(arr)
    coords = Array{Any}($N)
    # assume that the types are all the same.
    fill!(coords, Colon())
    coords[$slice_indices] = 1
    testslice = $slice_exp
    # testres is the first component.
    testres = nalift(f(testslice))
    reseltype = typeof(testres)

    ressize = sizearr[$slice_indices]
    result = if reseltype<:LabeledArray
      ndimstestres = ndims(testres)
      large_result = LabeledArray(similar(peel(testres), (size(testres)...,ressize...)), ntuple(mapslices_lambda_axis_tuple_func(testres, ressize, ndimstestres), ndimstestres+length(ressize)))
      colons = ntuple(mapslices_lambda_colon, ndims(testres))
      ST = typeof(view(large_result, colons..., ntuple(mapslices_lambda_one, length(ressize))...))
      temp_result = similar(arr, ST, ressize)
      @nloops $slice_ndims i temp_result begin
        @nref($slice_ndims,temp_result,i) = view(large_result, colons..., @ntuple($slice_ndims,i)...)
      end
      # it's important to assign this mutable function's return to the same variable because of fallback.
      newresult = mapslices_darr_larr_inner_typed_wrapper!(temp_result, f, arr, dims, testres)
      if isnull(newresult)
        create_labeled_array_mapslices_inner(large_result, ndimstestres, arr, $slice_indices)
      else
        newresult.value
      end
    elseif reseltype<:AbstractArray
      ndimstestres = ndims(testres)
      large_result = similar(testres, (size(testres)...,ressize...))
      colons = ntuple(mapslices_lambda_colon, ndims(testres))
      ST = typeof(view(large_result, colons..., ntuple(mapslices_lambda_one, length(ressize))...))
      temp_result = similar(arr, ST, ressize)
      @nloops $slice_ndims i temp_result begin
        @nref($slice_ndims,temp_result,i) = view(large_result, colons..., @ntuple($slice_ndims,i)...)
      end
      # it's important to assign this mutable function's return to the same variable because of fallback.
      newresult = mapslices_darr_larr_inner_typed_wrapper!(temp_result, f, arr, dims, testres)
      if isnull(newresult)
        if isa(arr, LabeledArray)
          create_labeled_array_mapslices_inner(LabeledArray(large_result), ndimstestres, arr, $slice_indices)
        else
          large_result
        end
      else
        newresult.value
      end
    else
      res = similar_mapslices_inner(arr, testres, $slice_indices)
      newresult = mapslices_darr_larr_inner_typed_wrapper!(res, f, arr, dims, testres)
      isnull(newresult) ? res : newresult.value
    end
    result
  end
end

mapslices_darr_larr_inner_typed_wrapper!{M,K,N,T,U<:AbstractArray,V}(result::AbstractArray{U,M}, f::Function, arr::AbstractArray{K,N}, dims::T, testres::V) =
  mapslices_darr_larr_inner_typed!(result, f, arr, dims, testres)

@generated mapslices_darr_larr_inner_typed!{M,K,N,T,U<:AbstractArray,V}(result::AbstractArray{U,M}, f::Function, arr::AbstractArray{K,N}, dims::T, testres::V) = begin
  dimtypes = dims.types
  slice_ndims = foldl(mapslices_lambda_acc_int, 0, dimtypes)
  slice_indices = Int[map(mapslices_lambda_first, Iterators.filter(mapslices_lambda_second_int, enumerate(dimtypes)))...]
  slice_exp = if slice_ndims == N
    :(arr[coords...])
  else
    :(view(arr, coords...))
  end
  quote
    sizearr = size(arr)
    coords = Array{Any}($N)
    # assume that the types are all the same.
    #ldict_keys_sofar = Nullable{Vector{KK}}() #testres.keys
    size_sofar = size(testres)
    is_the_first = true
    same_size::Bool = true
    counts_so_far = 0
    @nloops $slice_ndims i j->1:sizearr[$slice_indices[j]] begin
      fill!(coords, Colon())
      @nexprs $slice_ndims j->(coords[$slice_indices[j]] = i_j)
      oneslice = $slice_exp
      oneres = if is_the_first
        is_the_first = false
        testres
      else
        nalift(f(oneslice))
      end
      #same_size &= size_sofar == size(oneres)
      if size_sofar != size(oneres)
        fallback_result = similar(arr, typeof(testres), sizearr[$slice_indices])
        return Nullable(mapslices_darr_larr_inner_typed_fallback_wrapper!(fallback_result, result, f, arr, dims, counts_so_far))
      end
      #Haven't implemented the case when same_size==false yet.
      try
        copy_mapslice_helper!(@nref($slice_ndims, result, i), oneres)
      catch
        # if an exception occurs, just try once again with the fallback version.
        fallback_result = similar(arr, promote_type(typeof(oneres), typeof(testres)), sizearr[$slice_indices])
        return Nullable(mapslices_darr_larr_inner_typed_fallback_wrapper!(fallback_result, result, f, arr, dims, counts_so_far))
      end
      counts_so_far += 1
    end
    Nullable()
  end
end

# TODO: VV and VVV should be actually the same, but what to do with VV=Nullable{Any} and VVV=Nullable?
mapslices_darr_larr_inner_typed_wrapper!{KK,VV,VVV,M,K,N,T}(result::AbstractArray{LDict{KK,VV},M}, f::Function, arr::AbstractArray{K,N}, dims::T, testres::LDict{KK,VVV}) =
  mapslices_darr_larr_inner_typed!(result, f, arr, dims, testres)

@generated mapslices_darr_larr_inner_typed!{KK,VV,VVV,M,K,N,T}(result::AbstractArray{LDict{KK,VV},M}, f::Function, arr::AbstractArray{K,N}, dims::T, testres::LDict{KK,VVV}) = begin
  dimtypes = dims.types
  slice_ndims = foldl(mapslices_lambda_acc_int, 0, dimtypes)
  slice_indices = Int[map(mapslices_lambda_first, Iterators.filter(mapslices_lambda_second_int, enumerate(dimtypes)))...]
  slice_exp = if slice_ndims == N
    :(arr[coords...])
  else
    :(view(arr, coords...))
  end
  quote
    sizearr = size(arr)
    coords = Array{Any}($N)
    # assume that the types are all the same.
    ldict_keys_sofar = Nullable{Vector{KK}}() #testres.keys
    same_ldict::Bool = true
    is_the_first = true
    counts_so_far = 0
    @nloops $slice_ndims i j->1:sizearr[$slice_indices[j]] begin
      fill!(coords, Colon())
      @nexprs $slice_ndims j->(coords[$slice_indices[j]] = i_j)
      oneslice = $slice_exp
      oneres = if is_the_first
        is_the_first = false
        testres
      else
        nalift(f(oneslice))
      end
      same_ldict &= (isnull(ldict_keys_sofar) || ldict_keys_sofar.value == oneres.keys)
      if same_ldict
        if isnull(ldict_keys_sofar)
          ldict_keys_sofar = Nullable(oneres.keys)
        end
      else
        fallback_result = similar(arr, typeof(testres), sizearr[$slice_indices])
        return Nullable(mapslices_darr_larr_inner_typed_fallback_wrapper!(fallback_result, result, f, arr, dims, counts_so_far))
      end
      try
        @nref($slice_ndims, result, i) = oneres
      catch
        # if an exception occurs, just try once again with the fallback version.
        fallback_result = similar(arr, promote_type(typeof(oneres), typeof(testres)), sizearr[$slice_indices])
        return Nullable(mapslices_darr_larr_inner_typed_fallback_wrapper!(fallback_result, result, f, arr, dims, counts_so_far))
      end
      counts_so_far += 1
    end
    Nullable()
  end
end

mapslices_darr_larr_inner_typed_wrapper!{U,M,K,N,T}(result::AbstractArray{U,M}, f::Function, arr::AbstractArray{K,N}, dims::T, testres::U) =
  mapslices_darr_larr_inner_typed!(result, f, arr, dims, testres)

@generated mapslices_darr_larr_inner_typed!{U,M,K,N,T}(result::AbstractArray{U,M}, f::Function, arr::AbstractArray{K,N}, dims::T, testres::U) = begin
  dimtypes = dims.types
  slice_ndims = foldl(mapslices_lambda_acc_int, 0, dimtypes)
  slice_indices = Int[map(mapslices_lambda_first, Iterators.filter(mapslices_lambda_second_int, enumerate(dimtypes)))...]
  slice_exp = if slice_ndims == N
    :(arr[coords...])
  else
    :(view(arr, coords...))
  end
  quote
    sizearr = size(arr)
    coords = Array{Any}($N)
    is_the_first = true
    counts_so_far = 0
    @nloops $slice_ndims i j->1:sizearr[$slice_indices[j]] begin
      fill!(coords, Colon())
      @nexprs $slice_ndims j->(coords[$slice_indices[j]] = i_j)
      oneslice = $slice_exp
      oneres = if is_the_first
        is_the_first = false
        testres
      else
        nalift(f(oneslice))
      end
      try
        @nref($slice_ndims, result, i) = oneres
      catch
        # if an exception occurs, just try once again with the fallback version.
        fallback_result = similar(arr, promote_type(typeof(oneres), typeof(testres)), sizearr[$slice_indices])
        return Nullable(mapslices_darr_larr_inner_typed_fallback_wrapper!(fallback_result, result, f, arr, dims, counts_so_far))
      end
      counts_so_far += 1
    end
    #@rap wrap_array nalift result
    Nullable()
  end
end

mapslices_darr_larr_inner_typed_fallback_wrapper!{V,U,M,K,N,T}(fallback_result::AbstractArray{V,M}, result::AbstractArray{U,M}, f::Function, arr::AbstractArray{K,N}, dims::T, counts_so_far::Int) =
  mapslices_darr_larr_inner_typed_fallback!(fallback_result, result, f, arr, dims, counts_so_far)

@generated mapslices_darr_larr_inner_typed_fallback!{V,U,M,K,N,T}(fallback_result::AbstractArray{V,M}, result::AbstractArray{U,M}, f::Function, arr::AbstractArray{K,N}, dims::T, counts_so_far::Int) = begin
  @debug_stmt @show "mapslices_darr_larr_inner_typed! exception occurs presumably due to type mismatch. a fallback version will run."
  dimtypes = dims.types
  slice_ndims = foldl(mapslices_lambda_acc_int, 0, dimtypes)
  slice_indices = Int[map(mapslices_lambda_first, Iterators.filter(mapslices_lambda_second_int, enumerate(dimtypes)))...]
  slice_exp = if slice_ndims == N
    :(arr[coords...])
  else
    :(view(arr, coords...))
  end
  quote
    sizearr = size(arr)
    coords = Array{Any}($N)
    counts = 0
    @nloops $slice_ndims i j->1:sizearr[$slice_indices[j]] begin
      fill!(coords, Colon())
      @nexprs $slice_ndims j->(coords[$slice_indices[j]] = i_j)
      oneslice = $slice_exp
      icoords = CartesianIndex(@ntuple($slice_ndims, i))
      try
        if counts < counts_so_far
          setindex_mapslices_inner!(fallback_result, @nref($slice_ndims, result, i), icoords)
          #@nref($slice_ndims, fallback_result, i) = @nref($slice_ndims, result, i)
        else
          setindex_mapslices_inner!(fallback_result, nalift(f(oneslice)), icoords)
        end
      catch
        # if an exception occurs, just try once again with the fallback version.
        @debug_stmt @show "mapslices_darr_larr_inner_typed! exception occurs presumably due to type mismatch. a fallback version will run."
        oneres = nalift(f(oneslice))
        if promote_type(typeof(oneres), V) == V
          rethrow()
        end
        fallback_result = similar(arr, promote_type(typeof(oneres), V), sizearr[$slice_indices])
        return mapslices_darr_larr_inner_typed_fallback_wrapper!(fallback_result, result, f, arr, dims, counts)
      end
      counts += 1
    end
    fallback_result
  end
end

similar_mapslices_inner(arr::AbstractArray, testres, slice_indices) = similar(arr, typeof(testres), size(arr)[slice_indices])
similar_mapslices_inner(arr::DictArray, testres::LDict, slice_indices) = begin
  ressize = size(arr)[slice_indices]
  i = 0
  @rap create_dictarray_nocheck mapvalues(testres) do r
    i += 1
    similar(arr, typeof(r), ressize)
  end
end
similar_mapslices_inner(arr::LabeledArray, testres, slice_indices) = begin
  newdata = similar_mapslices_inner(peel(arr), testres, slice_indices)
  newaxes = ntuple(length(slice_indices)) do d
    pickaxis(arr, slice_indices[d])
  end
  LabeledArray(newdata, newaxes)
end

create_labeled_array_mapslices_inner(large_result::LabeledArray, ndimstestres, arr::LabeledArray, slice_indices) = begin
  LabeledArray(peel(large_result), ntuple(ndims(large_result)) do d
    if d <= ndimstestres
      pickaxis(large_result, d)
    else
      pickaxis(arr, slice_indices[d-ndimstestres])
    end
  end)
end

create_labeled_array_mapslices_inner(large_result::LabeledArray, ndimstestres, arr::AbstractArray, slice_indices) = begin
  LabeledArray(peel(large_result), ntuple(ndims(large_result)) do d
    if d <= ndimstestres
      pickaxis(large_result, d)
    else
      DefaultAxis(size(large_result, d))
    end
  end)
end
create_labeled_array_mapslices_inner(large_result::AbstractArray, ndimstestres, arr::LabeledArray, slice_indices) = begin
  LabeledArray(large_result, ntuple(ndims(large_result)) do d
    if d <= ndimstestres
      DefaultAxis(size(large_result, d))
    else
      pickaxis(arr, slice_indices[d-ndimstestres])
    end
  end)
end
create_labeled_array_mapslices_inner(large_result::AbstractArray, ndimstestres, arr::AbstractArray, slice_indices) = large_result

setindex_mapslices_inner!(fallback_result, reselem, coords) = (fallback_result[coords]=reselem)
setindex_mapslices_inner!{V<:AbstractArray}(fallback_result::AbstractArray{V}, reselem::AbstractArray, coords) = begin
  result = similar(fallback_result[1], size(reselem)) #similar(reselem, eltype(fallback_result))
  copy!(result, reselem)
  fallback_result[coords] = result
end
#@generated setindex_mapslices_inner!{V<:AbstractArray,N,T,M}(fallback_result::AbstractArray{V,N}, #reselem::AbstractArray{T,M}, coords) = quote #copy!(fallback_result[coords], reselem)
#  @nloops $M i reselem begin
#    fallback


text:
using ConcreteAbstractions
using Base.Test

@base type AbstractFoo{T}
    a
    b::Int
    c::T
    d::Vector{T}
end

@extend type Foo <: AbstractFoo
    e::T
end

@testset begin
    @test Foo <: AbstractFoo

    foo = Foo(10,10,10,[10],5)
    @test isa(foo, Foo{Int})
    @test fieldnames(foo) == [:a,:b,:c,:d,:e]
    @test isa(foo.a, Int)
    @test isa(foo.b, Int)
    @test isa(foo.c, Int)
    @test isa(foo.d, Vector{Int})
    @test isa(foo.e, Int)
end


text:
# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors    #src
# This Source Code Form is subject to the terms of the Mozilla Public License   #src
# v.2.0. If a copy of the MPL was not distributed with this file, You can       #src
# obtain one at https://mozilla.org/MPL/2.0/.                                   #src

# # The knapsack problem

# Formulate and solve a simple knapsack problem:
#
#     max sum(p_j x_j)
#      st sum(w_j x_j) <= C
#         x binary

using JuMP
import GLPK
import Test

function example_knapsack(; verbose = true)
    profit = [5, 3, 2, 7, 4]
    weight = [2, 8, 4, 2, 5]
    capacity = 10
    model = Model(GLPK.Optimizer)
    @variable(model, x[1:5], Bin)
    ## Objective: maximize profit
    @objective(model, Max, profit' * x)
    ## Constraint: can carry all
    @constraint(model, weight' * x <= capacity)
    ## Solve problem using MIP solver
    optimize!(model)
    if verbose
        println("Objective is: ", objective_value(model))
        println("Solution is:")
        for i in 1:5
            print("x[$i] = ", value(x[i]))
            println(", p[$i]/w[$i] = ", profit[i] / weight[i])
        end
    end
    Test.@test termination_status(model) == MOI.OPTIMAL
    Test.@test primal_status(model) == MOI.FEASIBLE_POINT
    Test.@test objective_value(model) == 16.0
    return
end

example_knapsack()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text: media_impact_fun(x,y) = exp(-x*y)