text:
using Logging, Random
using JuMP, Gurobi
using DecisionProgramming

Random.seed!(42)

@info("Creating the influence diagram.")
diagram = InfluenceDiagram()

add_node!(diagram, DecisionNode("DP", [], ["0-3 patents", "3-6 patents", "6-9 patents"]))
add_node!(diagram, ChanceNode("CT", ["DP"], ["low", "medium", "high"]))
add_node!(diagram, DecisionNode("DA", ["DP", "CT"], ["0-5 applications", "5-10 applications", "10-15 applications"]))
add_node!(diagram, ChanceNode("CM", ["CT", "DA"], ["low", "medium", "high"]))

generate_arcs!(diagram)

X_CT = ProbabilityMatrix(diagram, "CT")
X_CT[1, :] = [1/2, 1/3, 1/6]
X_CT[2, :] = [1/3, 1/3, 1/3]
X_CT[3, :] = [1/6, 1/3, 1/2]
add_probabilities!(diagram, "CT", X_CT)

X_CM = ProbabilityMatrix(diagram, "CM")
X_CM[1, 1, :] = [2/3, 1/4, 1/12]
X_CM[1, 2, :] = [1/2, 1/3, 1/6]
X_CM[1, 3, :] = [1/3, 1/3, 1/3]
X_CM[2, 1, :] = [1/2, 1/3, 1/6]
X_CM[2, 2, :] = [1/3, 1/3, 1/3]
X_CM[2, 3, :] = [1/6, 1/3, 1/2]
X_CM[3, 1, :] = [1/3, 1/3, 1/3]
X_CM[3, 2, :] = [1/6, 1/3, 1/2]
X_CM[3, 3, :] = [1/12, 1/4, 2/3]
add_probabilities!(diagram, "CM", X_CM)

generate_diagram!(diagram, default_utility=false)

@info("Creating the decision model.")
model = Model()
z = DecisionVariables(model, diagram)

@info("Creating problem specific constraints and expressions")

function variables(model::Model, dims::AbstractVector{Int}; binary::Bool=false)
    v = Array{VariableRef}(undef, dims...)
    for i in eachindex(v)
        v[i] = @variable(model, binary=binary)
    end
    return v
end

# Number of states in each node
n_DP = num_states(diagram, "DP")
n_CT = num_states(diagram, "CT")
n_DA = num_states(diagram, "DA")
n_CM = num_states(diagram, "CM")

n_T = 5                     # number of technology projects
n_A = 5                     # number of application projects
I_t = rand(n_T)*0.1         # costs of technology projects
O_t = rand(1:3,n_T)         # number of patents for each tech project
I_a = rand(n_T)*2           # costs of application projects
O_a = rand(2:4,n_T)         # number of applications for each appl. project

V_A = rand(n_CM, n_A).+0.5 # Value of an application
V_A[1, :] .+= -0.5          # Low market share: less value
V_A[3, :] .+= 0.5           # High market share: more value

x_T = variables(model, [n_DP, n_T]; binary=true)
x_A = variables(model, [n_DP, n_CT, n_DA, n_A]; binary=true)

M = 20                      # a large constant
ε = 0.5*minimum([O_t O_a])  # a helper variable, allows using ≤ instead of < in constraints (28b) and (29b)

q_P = [0, 3, 6, 9]          # limits of the technology intervals
q_A = [0, 5, 10, 15]        # limits of the application intervals
z_dP = z.z[1]
z_dA = z.z[2]

@constraint(model, [i=1:n_DP],
    sum(x_T[i,t] for t in 1:n_T) <= z_dP[i]*n_T)            #(25)
@constraint(model, [i=1:n_DP, j=1:n_CT, k=1:n_DA],
    sum(x_A[i,j,k,a] for a in 1:n_A) <= z_dP[i]*n_A)        #(26)
@constraint(model, [i=1:n_DP, j=1:n_CT, k=1:n_DA],
    sum(x_A[i,j,k,a] for a in 1:n_A) <= z_dA[i,j,k]*n_A)    #(27)

#(28a)
@constraint(model, [i=1:n_DP],
    q_P[i] - (1 - z_dP[i])*M <= sum(x_T[i,t]*O_t[t] for t in 1:n_T))
#(28b)
@constraint(model, [i=1:n_DP],
    sum(x_T[i,t]*O_t[t] for t in 1:n_T) <= q_P[i+1] + (1 - z_dP[i])*M - ε)
#(29a)
@constraint(model, [i=1:n_DP, j=1:n_CT, k=1:n_DA],
    q_A[k] - (1 - z_dA[i,j,k])*M <= sum(x_A[i,j,k,a]*O_a[a] for a in 1:n_A))
#(29b)
@constraint(model, [i=1:n_DP, j=1:n_CT, k=1:n_DA],
    sum(x_A[i,j,k,a]*O_a[a] for a in 1:n_A) <= q_A[k+1] + (1 - z_dA[i,j,k])*M - ε)

@constraint(model, [i=1:n_DP, j=1:n_CT, k=1:n_DA], x_A[i,j,k,1] <= x_T[i,1])
@constraint(model, [i=1:n_DP, j=1:n_CT, k=1:n_DA], x_A[i,j,k,2] <= x_T[i,1])
@constraint(model, [i=1:n_DP, j=1:n_CT, k=1:n_DA], x_A[i,j,k,2] <= x_T[i,2])

@info("Creating model objective.")
patent_investment_cost = @expression(model, [i=1:n_DP], sum(x_T[i, t] * I_t[t] for t in 1:n_T))
application_investment_cost = @expression(model, [i=1:n_DP, j=1:n_CT, k=1:n_DA], sum(x_A[i, j, k, a] * I_a[a] for a in 1:n_A))
application_value = @expression(model, [i=1:n_DP, j=1:n_CT, k=1:n_DA, l=1:n_CM], sum(x_A[i, j, k, a] * V_A[l, a] for a in 1:n_A))
@objective(model, Max, sum( sum( diagram.P(convert.(State, (i,j,k,l))) * (application_value[i,j,k,l] - application_investment_cost[i,j,k]) for j in 1:n_CT, k in 1:n_DA, l in 1:n_CM ) - patent_investment_cost[i] for i in 1:n_DP ))

@info("Starting the optimization process.")
optimizer = optimizer_with_attributes(
    () -> Gurobi.Optimizer(Gurobi.Env()),
    "IntFeasTol"      => 1e-9,
    "LazyConstraints" => 1,
)
set_optimizer(model, optimizer)
optimize!(model)

@info("Extracting results.")
Z = DecisionStrategy(z)
S_probabilities = StateProbabilities(diagram, Z)

@info("Printing decision strategy:")
print_decision_strategy(diagram, Z, S_probabilities)

@info("Extracting path utilities")
struct PathUtility <: AbstractPathUtility
    data::Array{AffExpr}
end
Base.getindex(U::PathUtility, i::State) = getindex(U.data, i)
Base.getindex(U::PathUtility, I::Vararg{State,N}) where N = getindex(U.data, I...)
(U::PathUtility)(s::Path) = value.(U[s...])

path_utility = [@expression(model,
    sum(x_A[s[index_of(diagram, "DP")], s[index_of(diagram, "CT")], s[index_of(diagram, "DA")], a] * (V_A[s[index_of(diagram, "CM")], a] - I_a[a]) for a in 1:n_A) -
    sum(x_T[s[index_of(diagram, "DP")], t] * I_t[t] for t in 1:n_T)) for s in paths(diagram.S)]
diagram.U = PathUtility(path_utility)

@info("Computing utility distribution.")
U_distribution = UtilityDistribution(diagram, Z)

@info("Printing utility distribution.")
print_utility_distribution(U_distribution)

@info("Printing statistics")
print_statistics(U_distribution)


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMIFlux
using Test
using FMIZoo

using FMIFlux.FMIImport: fmi2StringToValueReference, fmi2ValueReference
using FMIFlux: fmi2GetSolutionState, fmi2GetSolutionValue, fmi2GetSolutionTime

exportingToolsWindows = [("Dymola", "2022x")]
exportingToolsLinux = [("Dymola", "2022x")]

function runtests(exportingTool)
    ENV["EXPORTINGTOOL"] = exportingTool[1]
    ENV["EXPORTINGVERSION"] = exportingTool[2]

    @testset "Testing FMUs exported from $(ENV["EXPORTINGTOOL"]) ($(ENV["EXPORTINGVERSION"]))" begin

        @testset "Sensitivities" begin
            include("sens.jl")
        end
        @testset "ME-NeuralFMU (Continuous)" begin
            include("hybrid_ME.jl")
        end
        @testset "ME-NeuralFMU (Discontinuous)" begin
            include("hybrid_ME_dis.jl")
        end
        @testset "CS-NeuralFMU" begin
            include("hybrid_CS.jl")
        end
        @testset "Multiple FMUs" begin
            include("multi.jl")
        end
        @testset "Training modes" begin
            include("train_modes.jl")
        end
    end
end

@testset "FMIFlux.jl" begin
    if Sys.iswindows()
        @info "Automated testing is supported on Windows."
        for exportingTool in exportingToolsWindows
            runtests(exportingTool)
        end
    elseif Sys.islinux()
        @info "Automated testing is supported on Linux."
        for exportingTool in exportingToolsLinux
            runtests(exportingTool)
        end
    elseif Sys.isapple()
        @warn "Test-sets are currrently using Windows- and Linux-FMUs, automated testing for macOS is currently not supported."
    end
end

ext: jl


text:
using Pumas, PumasTutorials, CSV

data7BLQ = CSV.read("./tutorials/nca/data/single_dose_IVbolus_7BLQ.csv",missingstring="NA")
data7BLQ

data18BLQ = CSV.read("./tutorials/nca/data/single_dose_IVbolus_18BLQ.csv",missingstring="NA")
data18BLQ

timeu = u"hr"
concu = u"mg/L"
amtu  = u"mg"

pop7 = read_nca(data7BLQ, id=:ID, time=:time, conc=:DV, amt=:DOSE, ii=24timeu,
    route=:Formulation,timeu=timeu, concu=concu, amtu=amtu)
pop18 = read_nca(data18BLQ, id=:ID, time=:time, conc=:DV, amt=:DOSE, ii=24timeu,
    route=:Formulation, timeu=timeu, concu=concu, amtu=amtu)

pop7 = read_nca(data7BLQ, id=:ID, time=:time, conc=:DV, amt=:DOSE, ii=24timeu,
    route=:Formulation,timeu=timeu, concu=concu, amtu=amtu,llq=0.4concu)
pop18 = read_nca(data18BLQ, id=:ID, time=:time, conc=:DV, amt=:DOSE, ii=24timeu,
    route=:Formulation,timeu=timeu, concu=concu, amtu=amtu,llq=0.4concu)

NCA.auc(pop7,auctype=:last,method=:linear)
NCA.auc(pop18,auctype=:last,method=:linear)

NCA.auc(pop7,auctype=:inf,method=:linuplogdown)

NCA.auc_extrap_percent(pop7)

NCA.auc(pop7[1], interval=(0,12).*timeu)

pop7[1]

NCA.auc(pop7[1], interval=[(0,12).*timeu,(0,6).*timeu])

NCA.lambdaz(pop7)

NCA.lambdaz(pop7[1], threshold=3)

NCA.lambdaz(pop7[1], idxs=[18,19,20])

NCA.lambdaz(pop7[1], slopetimes=[18.5,19,19.5].*timeu)

NCA.interpextrapconc(pop7, 22timeu, method=:linear)

cmax = NCA.cmax(pop7[1])

normalizedose(cmax,pop7[1])

AUClast = NCA.auc(pop7[1],auctype=:last)
normalizedose(AUClast,pop7[1])

NCA.lambdazr2(pop7)
NCA.lambdazadjr2(pop7)
NCA.lambdazintercept(pop7)
NCA.lambdaztimefirst(pop7)
NCA.lambdaznpoints(pop7)

NCA.tmax(pop7)
NCA.cmin(pop7)
NCA.tmin(pop7)

NCA.tlast(pop7)
NCA.clast(pop7)

NCA.aumc(pop7)
NCA.aumclast(pop7)

NCA.thalf(pop7)

NCA.cl(pop7)

NCA.vss(pop7)
NCA.vz(pop7)

NCA.accumulationindex(pop7)

pop7[24]

NCA.cl(pop7[24])

report = NCAReport(pop7)
report = NCA.to_dataframe(report)

report = NCAReport(pop7,pred=true)
report = NCA.to_dataframe(report)

CSV.write("./tutorials/nca/report_SD_IVbolus_7BLQ.csv", report)


text:
# ## Data

import CSV
import Random
import TuringModels

using DataFrames

Random.seed!(1)

file_path = joinpath(TuringModels.project_root, "data", "chimpanzees.csv")
df = CSV.read(file_path, DataFrame; delim=';');

# ## Model

using StatsFuns
using Turing

@model function m10_3(y, x₁, x₂)
    α ~ Normal(0, 10)
    βp ~ Normal(0, 10)
    βpC ~ Normal(0, 10)

    logits = α .+ (βp .+ βpC * x₁) .* x₂
    y .~ BinomialLogit.(1, logits)
end

# ## Output

model = m10_3(df.pulled_left, df.condition, df.prosoc_left)
chns = sample(model, NUTS(), 2000)

# \defaultoutput{}

# ## Original output

"""
      Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
 a    0.05   0.13      -0.15       0.25  3284    1
 bp   0.62   0.22       0.28       0.98  3032    1
 bpC -0.11   0.26      -0.53       0.29  3184    1
""";


text:
include("model.jl")
using Flux, Flux.Tracker
import Random

Random.seed!(0);

# Generate Data
N = 1000
x = randn(N)
b0 = .5
b1 = 2.0
p = sigmoid.(b0 .+ b1 * x)
y = (p .> rand(N)) * 1.0

minibatch_size = 500
vp = VP()
metrics = Vector{Metrics}()
loss(y, x) = -elbo(y, x, vp, N, metrics) / N

params = Tracker.Params([getfield(vp, fn) for fn in fieldnames(typeof(vp))])
# grads = Tracker.gradient(() -> loss(y, x), params)

opt = ADAM(1e-1)
niters = 1000
Random.seed!(3);

@time for i in 1:niters
  idx = sample(1:N, minibatch_size)
  Flux.train!(loss, params, [(y[idx], x[idx])], opt)
  if i % (niters / 10) == 0
    println("Progress: $(i)/$(niters) | ELBO: $(-loss(y, x).data)")
  end
end

println("b0_true: $(b0) | b1_true: $(b1)")
println("b0_mean: $(vp.b0[1].data) | b1_mean: $(vp.b1[1].data)")
println("b0_sd: $(exp(vp.b0[2].data)) | b1_sd: $(exp(vp.b1[2].data))")


text:
module MultiFieldSparseMatrixAssemblersTests

using Gridap.Arrays
using Gridap.Algebra
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.Fields
using Gridap.Integration
using SparseArrays
using Gridap.CellData
using Gridap.MultiField
using Gridap.ReferenceFEs
using Gridap.TensorValues

using Test

order = 2

domain = (0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

degree = order
trian = get_triangulation(model)
quad = CellQuadrature(trian,degree)

trian_Γ = SkeletonTriangulation(model)
quad_Γ = CellQuadrature(trian_Γ,degree)

V = TestFESpace(model,ReferenceFE(lagrangian,Float64,order);conformity=:H1)
Q = TestFESpace(model,ReferenceFE(lagrangian,Float64,order-1),conformity=:L2)

U = TrialFESpace(V)
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V,Q])
X = MultiFieldFESpace([U,P])

free_values = rand(num_free_dofs(X))
xh = FEFunction(X,free_values)
uh, ph = xh

dy = get_cell_shapefuns(Y)
dx = get_cell_shapefuns_trial(X)
dv, dq = dy
du, dp = dx

cellmat = integrate(dv*du,quad)
cellvec = integrate(dv*2,quad)
cellids = get_cell_to_bgcell(trian)
cellmatvec = pair_arrays(cellmat,cellvec)

#cellmat_Γ = integrate(  jump(dv)*dp.⁺ + mean(dq)*jump(dp), quad_Γ)
cellmat_Γ = integrate(  jump(dv)*mean(du) + jump(∇(dq))⋅jump(∇(dp)), quad_Γ)

cellvec_Γ = integrate(  jump(dv) + mean(dq),quad_Γ)
cellmatvec_Γ = pair_arrays(cellmat_Γ,cellvec_Γ)
cellids_Γ = get_cell_to_bgcell(trian_Γ)

assem = SparseMatrixAssembler(SparseMatrixCSR{0,Float64,Int},X,Y)

matvecdata = ([cellmatvec,cellmatvec_Γ],[cellids,cellids_Γ],[cellids,cellids_Γ])
matdata = ([cellmat,cellmat_Γ],[cellids,cellids_Γ],[cellids,cellids_Γ])
vecdata = ([cellvec,cellvec_Γ],[cellids,cellids_Γ])
data = (matvecdata,matdata,vecdata)

test_assembler(assem,matdata,vecdata,data)

A = assemble_matrix(assem,matdata)

A = allocate_matrix(assem,matdata)

assem = SparseMatrixAssembler(X,Y)
test_assembler(assem,matdata,vecdata,data)

end # module


text:
__precompile__()

module Compat

using Base.Meta

"""Get just the function part of a function declaration."""
withincurly(ex) = isexpr(ex, :curly) ? ex.args[1] : ex

if VERSION < v"0.5.0-dev+961"
    export walkdir

    function walkdir(root; topdown=true, follow_symlinks=false, onerror=throw)
        content = nothing
        try
            content = readdir(root)
        catch err
            isa(err, SystemError) || throw(err)
            onerror(err)
            #Need to return an empty task to skip the current root folder
            return Task(()->())
        end
        dirs = Array(eltype(content), 0)
        files = Array(eltype(content), 0)
        for name in content
            if isdir(joinpath(root, name))
                push!(dirs, name)
            else
                push!(files, name)
            end
        end

        function _it()
            if topdown
                produce(root, dirs, files)
            end
            for dir in dirs
                path = joinpath(root,dir)
                if follow_symlinks || !islink(path)
                    for (root_l, dirs_l, files_l) in walkdir(path, topdown=topdown, follow_symlinks=follow_symlinks, onerror=onerror)
                        produce(root_l, dirs_l, files_l)
                    end
                end
            end
            if !topdown
                produce(root, dirs, files)
            end
        end
        Task(_it)
    end
end
if VERSION < v"0.6.0-dev.2043"
    Base.take!(t::Task) = consume(t)
end

function rewrite_show(ex)
    if isexpr(ex, :call)
        Expr(:call, rewrite_show(ex.args[1]), ex.args[2:end]...)
    elseif isexpr(ex, :curly)
        Expr(:curly, rewrite_show(ex.args[1]), ex.args[2:end]...)
    else
        :(Base.writemime)
    end
end

function rewrite_dict(ex)
    length(ex.args) == 1 && return ex

    f = ex.args[1]
    if isexpr(f, :curly)
        newex = Expr(:typed_dict, :($(f.args[2])=>$(f.args[3])))
    else
        newex = Expr(:dict)
    end

    for i = 2:length(ex.args)
        pair = ex.args[i]
        !isexpr(pair, :(=>)) && return ex
        push!(newex.args, pair)
    end
    newex
end

function rewrite_ordereddict(ex)
    length(ex.args) == 1 && return ex

    f = ex.args[1]
    newex = Expr(:call, f, :[])

    for i = 2:length(ex.args)
        pair = ex.args[i]
        !isexpr(pair, :(=>)) && return ex
        push!(newex.args[2].args, Expr(:tuple, pair.args...))
    end

    newex
end

# rewrite Julia 0.4-style split or rsplit (str, splitter; kws...)
# into 0.2/0.3-style positional arguments
function rewrite_split(ex, f)
    limit = nothing
    keep = nothing
    for i in 4:length(ex.args)
        if isexpr(ex.args[i], :kw)
            kw = ex.args[i].args
            if kw[1] == :limit
                limit = kw[2]
            elseif kw[1] == :keep
                keep = kw[2]
            end
        end
    end
    if limit == nothing
        if keep == nothing
            return Expr(:call, f, ex.args[2], ex.args[3])
        else
            return Expr(:call, f, ex.args[2], ex.args[3], keep)
        end
    else
        if keep == nothing
            return Expr(:call, f, ex.args[2], ex.args[3], limit)
        else
            return Expr(:call, f, ex.args[2], ex.args[3], limit, keep)
        end
    end
end

# rewrites all subexpressions of the form `a => b` to `(a, b)`
function rewrite_pairs_to_tuples!(expr::Expr)
    if expr.head == :(=>)
        expr.head = :tuple
    end
    for subexpr in expr.args
        isa(subexpr, Expr) && rewrite_pairs_to_tuples!(subexpr)
    end
    return expr
end

function is_quote_symbol(ex::ANY, val::Symbol)
    if isa(ex, QuoteNode)
        return (ex::QuoteNode).value === val
    elseif isa(ex, Expr)
        ex = ex::Expr
        return ex.head === :quote && length(ex.args) == 1 && ex.args[1] === val
    end
    return false
end

is_index_style(ex::Expr) = ex == :(Compat.IndexStyle) || ex == :(Base.IndexStyle) ||
    (ex.head == :(.) && (ex.args[1] == :Compat || ex.args[1] == :Base) &&
         ex.args[2] == Expr(:quote, :IndexStyle))

is_index_style(arg) = false

# rewrites accesses to IOContext dicts
function rewrite_iocontext!(expr::Expr)
    args = expr.args
    nargs = length(args)
    if nargs == 4 && expr.head === :call && args[1] === :get && args[4] === false
        key = args[3]
        if is_quote_symbol(key, :limit) || is_quote_symbol(key, :compact)
            if VERSION >= v"0.5.0-dev+1936" && VERSION < v"0.5.0-dev+4305"
                args[1] = :(Base.limit_output)
                deleteat!(args, 3:4)
            elseif VERSION < v"0.5.0-dev+1936"
                expr.head = :quote
                args[1] = false
                deleteat!(args, 3:4)
            end
        elseif is_quote_symbol(key, :multiline)
            if VERSION < v"0.5.0-dev+4305"
                expr.head = :quote
                args[1] = false
                deleteat!(args, 3:4)
            end
        end
    end
end

# JuliaLang/julia#10543
if !isdefined(Base, :tryparse)
    function tryparse{T}(::Type{T}, args...)
        try
            Nullable(Base.parse(T, args...))
        catch
            Nullable{T}()
        end
    end
end

import Base.unsafe_convert

if VERSION < v"0.5.0-dev+3831"
    Base.Symbol(args...) = symbol(args...)::Symbol
end

if VERSION < v"0.5.0-dev+2396"
    function new_style_call_overload(ex::Expr)
        # Not a function call
        ((ex.head === :(=) || ex.head === :function) &&
         length(ex.args) == 2 && isexpr(ex.args[1], :call)) || return false
        callee = (ex.args[1]::Expr).args[1]
        # Only Expr function name can be call overload
        isa(callee, Expr) || return false
        callee = callee::Expr
        # (a::A)() = ...
        callee.head === :(::) && return true
        # The other case is with type parameter.
        # Filter out everything without one.
        (callee.head === :curly && length(callee.args) >= 1) || return false
        # Check what the type parameter applies to is a Expr(:(::))
        return isexpr(callee.args[1], :(::))
    end
else
    new_style_call_overload(ex::Expr) = false
end

istopsymbol(ex, mod, sym) = ex in (sym, Expr(:(.), mod, Expr(:quote, sym)))

if VERSION < v"0.5.0-dev+4002"
    @inline broadcast_getindex(arg, idx) = arg[(idx - 1) % length(arg) + 1]
    # Optimize for single element
    @inline broadcast_getindex(arg::Number, idx) = arg
    @inline broadcast_getindex{T}(arg::Array{T,0}, idx) = arg[1]

    # If we know from syntax level that we don't need wrapping
    @inline broadcast_getindex_naive(arg, idx) = arg[idx]
    @inline broadcast_getindex_naive(arg::Number, idx) = arg
    @inline broadcast_getindex_naive{T}(arg::Array{T,0}, idx) = arg[1]

    # For vararg support
    @inline getindex_vararg(idx) = ()
    @inline getindex_vararg(idx, arg1) = (broadcast_getindex(arg1, idx),)
    @inline getindex_vararg(idx, arg1, arg2) =
        (broadcast_getindex(arg1, idx), broadcast_getindex(arg2, idx))
    @inline getindex_vararg(idx, arg1, arg2, arg3, args...) =
        (broadcast_getindex(arg1, idx), broadcast_getindex(arg2, idx),
         broadcast_getindex(arg3, idx), getindex_vararg(idx, args...)...)

    @inline getindex_naive_vararg(idx) = ()
    @inline getindex_naive_vararg(idx, arg1) =
        (broadcast_getindex_naive(arg1, idx),)
    @inline getindex_naive_vararg(idx, arg1, arg2) =
        (broadcast_getindex_naive(arg1, idx),
         broadcast_getindex_naive(arg2, idx))
    @inline getindex_naive_vararg(idx, arg1, arg2, arg3, args...) =
        (broadcast_getindex_naive(arg1, idx),
         broadcast_getindex_naive(arg2, idx),
         broadcast_getindex_naive(arg3, idx),
         getindex_naive_vararg(idx, args...)...)

    # Decide if the result should be scalar or array
    # `size() === ()` is not good enough since broadcasting on
    # a scalar should return a scalar where as broadcasting on a 0-dim
    # array should return a 0-dim array.
    @inline should_return_array(::Val{true}, args...) = Val{true}()
    @inline should_return_array(::Val{false}) = Val{false}()
    @inline should_return_array(::Val{false}, arg1) = Val{false}()
    @inline should_return_array(::Val{false}, arg1::AbstractArray) = Val{true}()
    @inline should_return_array(::Val{false}, arg1::AbstractArray,
                                arg2::AbstractArray) = Val{true}()
    @inline should_return_array(::Val{false}, arg1,
                                arg2::AbstractArray) = Val{true}()
    @inline should_return_array(::Val{false}, arg1::AbstractArray,
                                arg2) = Val{true}()
    @inline should_return_array(::Val{false}, arg1, arg2) = Val{false}()
    @inline should_return_array(::Val{false}, arg1, arg2, args...) =
        should_return_array(should_return_array(Val{false}(), arg1, arg2),
                            args...)

    @inline broadcast_return(res1d, shp, ret_ary::Val{false}) = res1d[1]
    @inline broadcast_return(res1d, shp, ret_ary::Val{true}) = reshape(res1d, shp)

    @inline need_full_getindex(shp) = false
    @inline need_full_getindex(shp, arg1::Number) = false
    @inline need_full_getindex{T}(shp, arg1::Array{T,0}) = false
    @inline need_full_getindex(shp, arg1) = shp != size(arg1)
    @inline need_full_getindex(shp, arg1, arg2) =
        need_full_getindex(shp, arg1) || need_full_getindex(shp, arg2)
    @inline need_full_getindex(shp, arg1, arg2, arg3, args...) =
        need_full_getindex(shp, arg1, arg2) || need_full_getindex(shp, arg3) ||
        need_full_getindex(shp, args...)

    function rewrite_broadcast(f, args)
        nargs = length(args)
        # This actually allows multiple splatting...,
        # which is now allowed on master.
        # The previous version that simply calls broadcast so removing that
        # will be breaking. Oh, well....
        is_vararg = Bool[isexpr(args[i], :...) for i in 1:nargs]
        names = [gensym("broadcast") for i in 1:nargs]
        new_args = [is_vararg[i] ? Expr(:..., names[i]) : names[i]
                    for i in 1:nargs]
        # Optimize for common case where we know the index doesn't need
        # any wrapping
        naive_getidx_for = function (i, idxvar)
            if is_vararg[i]
                Expr(:..., :($Compat.getindex_naive_vararg($idxvar,
                                                           $(names[i])...)))
            else
                :($Compat.broadcast_getindex_naive($(names[i]), $idxvar))
            end
        end
        always_naive = nargs == 1 && !is_vararg[1]
        getidx_for = if always_naive
            naive_getidx_for
        else
            function (i, idxvar)
                if is_vararg[i]
                    Expr(:..., :($Compat.getindex_vararg($idxvar,
                                                         $(names[i])...)))
                else
                    :($Compat.broadcast_getindex($(names[i]), $idxvar))
                end
            end
        end
        @gensym allidx
        @gensym newshape
        @gensym res1d
        @gensym idx
        @gensym ret_ary

        res1d_expr = quote
            $res1d = [$f($([naive_getidx_for(i, idx) for i in 1:nargs]...))
                      for $idx in $allidx]
        end
        if !always_naive
            res1d_expr = quote
                if $Compat.need_full_getindex($newshape, $(new_args...))
                    $res1d = [$f($([getidx_for(i, idx) for i in 1:nargs]...))
                              for $idx in $allidx]
                else
                    $res1d_expr
                end
            end
        end

        return quote
            # The `local` makes sure type inference can infer the type even
            # in global scope as long as the input is type stable
            local $(names...)
            $([:($(names[i]) = $(is_vararg[i] ? args[i].args[1] : args[i]))
               for i in 1:nargs]...)
            local $newshape = $(Base.Broadcast).broadcast_shape($(new_args...))
            # `eachindex` is not generic enough
            local $allidx = 1:prod($newshape)
            local $ret_ary = $Compat.should_return_array(Val{false}(),
                                                         $(new_args...))
            local $res1d
            $res1d_expr
            $Compat.broadcast_return($res1d, $newshape, $ret_ary)
        end
    end
end

if VERSION < v"0.6.0-dev.2782"
    function new_style_typealias(ex::ANY)
        isexpr(ex, :(=)) || return false
        ex = ex::Expr
        return length(ex.args) == 2 && isexpr(ex.args[1], :curly)
    end
else
    new_style_typealias(ex) = false
end

function _compat(ex::Expr)
    if ex.head === :call
        f = ex.args[1]
        if VERSION < v"0.5.0-dev+4340" && length(ex.args) > 3 &&
                istopsymbol(withincurly(ex.args[1]), :Base, :show)
            ex = rewrite_show(ex)
        elseif VERSION < v"0.6.0-dev.826" && length(ex.args) == 3 && # julia#18510
                istopsymbol(withincurly(ex.args[1]), :Base, :Nullable)
            ex = Expr(:call, f, ex.args[2], Expr(:call, :(Compat._Nullable_field2), ex.args[3]))
        end
        if VERSION < v"0.5.0-dev+4305"
            rewrite_iocontext!(ex)
        end
    elseif ex.head === :curly
        f = ex.args[1]
        if ex == :(Ptr{Void})
            # Do not change Ptr{Void} to Ptr{Nothing}: 0.4.0-dev+768
            return ex
        elseif VERSION < v"0.6.0-dev.2575" #20414
            ex = Expr(:curly, map(a -> isexpr(a, :call, 2) && a.args[1] == :(<:) ?
                                  :($TypeVar($(QuoteNode(gensym(:T))), $(a.args[2]), false)) :
                                  isexpr(a, :call, 2) && a.args[1] == :(>:) ?
                                  :($TypeVar($(QuoteNode(gensym(:T))), $(a.args[2]), $Any, false)) : a,
                                  ex.args)...)
        end
    elseif ex.head === :macrocall
        f = ex.args[1]
        if VERSION < v"0.5.0-dev+2129" && f === symbol("@boundscheck") && length(ex.args) == 2
            # Handle 0.5 single argument @boundscheck syntax. (0.4 has a two
            # arguments and a very diffferent meaning).  `nothing` is included
            # so we have a consistent return type.
            ex = :($(ex.args[2]); nothing)
        end
    elseif ex.head === :quote && isa(ex.args[1], Symbol)
        # Passthrough
        return ex
    elseif new_style_call_overload(ex)
        callexpr = ex.args[1]::Expr
        callee = callexpr.args[1]::Expr
        is_kw = (length(callexpr.args) >= 2 &&
                 isexpr(callexpr.args[2], :parameters))
        if callee.head === :(::)
            # (:function, (:call, :(:(::), <1>), <2>), <body>) ->
            # (:function, (:call, :(Base.call), :(:(::), <1>), <2>), <body>)
            unshift!(callexpr.args, :(Base.call))
        else
            # (:function, (:call, :(curly, :(:(::), <1>), <3>), <2>), <body>) ->
            # (:function, (:call, :(curly, :(Base.call), <3>), :(:(::), <1>), <2>), <body>)
            obj = callee.args[1]::Expr
            callee.args[1] = :(Base.call)
            insert!(callexpr.args, 2, obj)
        end
        if is_kw
            # Expr(:parameters) is moved to the 3rd argument
            params = callexpr.args[3]
            @assert isexpr(params, :parameters)
            obj = callexpr.args[2]
            callexpr.args[2] = params
            callexpr.args[3] = obj
        end
    elseif VERSION < v"0.5.0-dev+4002" && ex.head == :. && length(ex.args) == 2 # 15032
        if isexpr(ex.args[2], :quote) && isa(ex.args[2].args[1], QuoteNode)
            # foo.:bar -> foo.(:bar) in older Julia
            return Expr(ex.head, _compat(ex.args[1]), ex.args[2].args[1])
        elseif isexpr(ex.args[2], :quote) && isexpr(ex.args[2].args[1], :quote) &&
               isa(ex.args[2].args[1].args[1], Symbol)
            # foo.:bar -> foo.(:bar) in older Julia
            return Expr(ex.head, _compat(ex.args[1]), QuoteNode(ex.args[2].args[1].args[1]))
        elseif isexpr(ex.args[2], :tuple)
            # f.(arg1, arg2...) -> broadcast(f, arg1, arg2...)
            return rewrite_broadcast(_compat(ex.args[1]),
                                     map(_compat, ex.args[2].args))
        elseif !isa(ex.args[2], QuoteNode) &&
               !(isexpr(ex.args[2], :quote) && isa(ex.args[2].args[1], Symbol))
            # f.(arg) -> broadcast(f, arg)
            return rewrite_broadcast(_compat(ex.args[1]), [_compat(ex.args[2])])
        end
    elseif new_style_typealias(ex)
        ex.head = :typealias
    elseif ex.head === :const && length(ex.args) == 1 && new_style_typealias(ex.args[1])
        ex = ex.args[1]::Expr
        ex.head = :typealias
    elseif ex.head === :import
        if VERSION < v"0.5.0-dev+4340" && length(ex.args) == 2 && ex.args[1] === :Base && ex.args[2] === :show
            return quote
                import Base.show
                import Base.writemime
            end
        end
    elseif VERSION < v"0.5.0-dev+5575" #17510
        if isexpr(ex, :comparison)
            if length(ex.args) == 3 && ex.args[2] == :.=
                return :(Base.broadcast!(Base.identity, $(_compat(todotview(ex.args[1]))), $(_compat(ex.args[3]))))
            elseif length(ex.args) > 3 && ex.args[2] == :.=
                return :(Base.broadcast!(Base.identity, $(_compat(todotview(ex.args[1]))), $(_compat(Expr(:comparison, ex.args[3:end]...)))))
            end
        elseif ex.head == :.= && length(ex.args) == 2 # may arise in macro-constructed expressions
            return :(Base.broadcast!(Base.identity, $(_compat(todotview(ex.args[1]))), $(_compat(ex.args[2]))))
        end
    end
    if VERSION < v"0.5.0-dev+5575" #17510
        # transform e.g. x .+= y into x .= x .+ y:
        shead = string(ex.head)
        if first(shead) == '.' && length(shead) > 2 && last(shead) == '=' && length(ex.args) == 2
            return _compat(Expr(:comparison, ex.args[1], :.=, Expr(:call, Symbol(shead[1:end-1]), ex.args...)))
        end
    end
    if VERSION < v"0.6.0-dev.2840"
        if ex.head == :(=) && isa(ex.args[1], Expr) && ex.args[1].head == :call
            a = ex.args[1].args[1]
            if is_index_style(a)
                ex.args[1].args[1] = :(Base.linearindexing)
            elseif isa(a, Expr) && a.head == :curly
                if is_index_style(a.args[1])
                    ex.args[1].args[1].args[1] = :(Base.linearindexing)
                end
            end
        end
    end
    return Expr(ex.head, map(_compat, ex.args)...)
end
function _compat(ex::Symbol)
    if VERSION < v"0.5.0-dev+1343" && (ex == :Future || ex == :RemoteChannel)
        return :RemoteRef
    end
    return ex
end
_compat(ex) = ex

function _get_typebody(ex::Expr)
    args = ex.args
    if ex.head !== :type || length(args) != 3 || args[1] !== true
        throw(ArgumentError("Invalid usage of @compat: $ex"))
    end
    name = args[2]
    if !isexpr(args[3], :block)
        throw(ArgumentError("Invalid type declaration: $ex"))
    end
    body = (args[3]::Expr).args
    filter!(body) do e
        if isa(e, LineNumberNode) || isexpr(e, :line)
            return false
        end
        return true
    end
    return name, body
end

function _compat_primitive(typedecl)
    name, body = _get_typebody(typedecl)
    if length(body) != 1
        throw(ArgumentError("Invalid primitive type declaration: $typedecl"))
    end
    return Expr(:bitstype, body[1], name)
end

function _compat_abstract(typedecl)
    name, body = _get_typebody(typedecl)
    if length(body) != 0
        throw(ArgumentError("Invalid abstract type declaration: $typedecl"))
    end
    return Expr(:abstract, name)
end

macro compat(ex...)
    if VERSION < v"0.6.0-dev.2746" && length(ex) == 2 && ex[1] === :primitive
        return esc(_compat_primitive(ex[2]))
    elseif length(ex) != 1
        throw(ArgumentError("@compat called with wrong number of arguments: $ex"))
    elseif (VERSION < v"0.6.0-dev.2746" && isexpr(ex[1], :abstract) &&
            length(ex[1].args) == 1 && isexpr(ex[1].args[1], :type))
        # This can in principle be handled in nested case but we do not
        # do that to be consistent with primitive types.
        return esc(_compat_abstract(ex[1].args[1]))
    end
    esc(_compat(ex[1]))
end

export @compat, @inline, @noinline

import Base.@irrational

import Base: remotecall, remotecall_fetch, remotecall_wait, remote_do
if VERSION < v"0.5.0-dev+431"
    for f in (:remotecall, :remotecall_fetch, :remotecall_wait, :remote_do)
        @eval begin
            ($f)(f::Function, w::Base.LocalProcess, args...)   = ($f)(w, f, args...)
            ($f)(f::Function, w::Base.Worker, args...)         = ($f)(w, f, args...)
            ($f)(f::Function, id::Integer, args...)            = ($f)(id, f, args...)
        end
    end
end

if VERSION < v"0.5.0-dev+763"
    export SparseArrays
    const SparseArrays = Base.SparseMatrix
end

if VERSION < v"0.5.0-dev+1229"
    const Filesystem = Base.FS
else
    import Base.Filesystem
end

if VERSION < v"0.5.0-dev+1946"
    const supertype = super
    export supertype
end

if VERSION < v"0.5.0-dev+679"
    import Base: cov, cor

    cov(x::AbstractVector, corrected::Bool) = cov(x, corrected=corrected)
    cov(X::AbstractMatrix, vardim::Integer) = cov(X, vardim=vardim)
    cov(X::AbstractMatrix, vardim::Integer, corrected::Bool) = cov(X, vardim=vardim, corrected=corrected)
    cov(x::AbstractVector, y::AbstractVector, corrected::Bool) = cov(x, y, corrected=corrected)
    cov(X::AbstractMatrix, Y::AbstractMatrix, vardim::Integer) = cov(X, Y, vardim=vardim)
    cov(X::AbstractMatrix, Y::AbstractMatrix, vardim::Integer, corrected::Bool) = cov(X, Y, vardim=vardim, corrected=corrected)

    cor(X::AbstractMatrix, vardim::Integer) = cor(X, vardim=vardim)
    cor(X::AbstractMatrix, Y::AbstractMatrix, vardim::Integer) = cor(X, Y, vardim=vardim)
end

if VERSION < v"0.5.0-dev+2228"
    const readstring = readall
    export readstring

    Base.read(s::IO) = readbytes(s)
    Base.read(s::IO, nb) = readbytes(s, nb)

    Base.write(filename::AbstractString, args...) = open(io->write(io, args...), filename, "w")

    Base.read(filename::AbstractString, args...) = open(io->read(io, args...), filename)
    Base.read!(filename::AbstractString, a) = open(io->read!(io, a), filename)
    Base.readuntil(filename::AbstractString, args...) = open(io->readuntil(io, args...), filename)
    Base.readline(filename::AbstractString) = open(Base.readline, filename)
    Base.readlines(filename::AbstractString) = open(Base.readlines, filename)
    Base.readavailable(s::IOStream) = read!(s, @compat Vector{UInt8}(nb_available(s)))
    Base.readavailable(s::IOBuffer) = read(s)

    function Base.write(to::IO, from::IO)
        while !eof(from)
            write(to, readavailable(from))
        end
    end

    function Base.eachline(filename::AbstractString)
        s = open(filename)
        EachLine(s, ()->close(s))
    end
end

if VERSION < v"0.5.0-dev+2023"
    const displaysize = Base.tty_size
    export displaysize
end

# Pull Request https://github.com/JuliaLang/julia/pull/13232
# Rounding and precision functions:

if VERSION >= v"0.5.0-dev+1182"
    import Base:
        setprecision, setrounding, rounding

else  # if VERSION < v"0.5.0-dev+1182"

    export setprecision
    export setrounding
    export rounding

    setprecision(f, ::Type{BigFloat}, prec) = with_bigfloat_precision(f, prec)
    setprecision(::Type{BigFloat}, prec) = set_bigfloat_precision(prec)

    # assume BigFloat if type not explicit:
    setprecision(prec) = setprecision(BigFloat, prec)
    setprecision(f, prec) = setprecision(f, BigFloat, prec)

    Base.precision(::Type{BigFloat}) = get_bigfloat_precision()

    setrounding(f, T, rounding_mode) =
        with_rounding(f, T, rounding_mode)

    setrounding(T, rounding_mode) = set_rounding(T, rounding_mode)

    rounding(T) = get_rounding(T)

end

module LinAlg
    if VERSION < v"0.5.0-dev+2022"
        const checksquare = Base.LinAlg.chksquare
    else
        import Base.LinAlg.checksquare
    end
end

if VERSION < v"0.5.0-dev+2915"
    const issymmetric = issym
    export issymmetric
end

if VERSION < v"0.5.0-dev+977"
    export foreach

    foreach(f) = (f(); nothing)
    foreach(f, itr) = (for x in itr; f(x); end; nothing)
    foreach(f, itrs...) = (for z in zip(itrs...); f(z...); end; nothing)
end

if !isdefined(Base, :istextmime)
    export istextmime
    istextmime(m::@compat(Union{MIME,AbstractString})) = istext(m)
end

function primarytype(t::ANY)
    tn = t.name
    if isdefined(tn, :primary)
        return tn.primary
    else
        return tn.wrapper
    end
end

export @functorize
macro functorize(f)
    if VERSION >= v"0.5.0-dev+3701"
        f === :scalarmax       ? :(Base.scalarmax) :
        f === :scalarmin       ? :(Base.scalarmin) :
        f === :centralizedabs2fun ? :(primarytype(typeof(Base.centralizedabs2fun(0)))) :
        f
    else
        f = f === :identity        ? :(Base.IdFun()) :
            f === :abs             ? :(Base.AbsFun()) :
            f === :abs2            ? :(Base.Abs2Fun()) :
            f === :exp             ? :(Base.ExpFun()) :
            f === :log             ? :(Base.LogFun()) :
            f === :&               ? :(Base.AndFun()) :
            f === :|               ? :(Base.OrFun()) :
            f === :+               ? :(Base.AddFun()) :
            f === :*               ? :(Base.MulFun()) :
            f === :scalarmax       ? :(Base.MaxFun()) :
            f === :scalarmin       ? :(Base.MinFun()) :
            f === :centralizedabs2fun ? :(Base.CentralizedAbs2Fun) :
            f
        if VERSION >= v"0.4.0-dev+4902"
            f = f === :<           ? :(Base.LessFun()) :
                f === :>           ? :(Base.MoreFun()) :
                f
        end
        if VERSION >= v"0.4.0-dev+4902"
            f = f === :conj        ? :(Base.ConjFun()) :
                f
        end
        if VERSION >= v"0.4.0-dev+6254"
            f = f === :-           ? :(Base.SubFun()) :
                f === :^           ? :(Base.PowFun()) :
                f
        end
        if VERSION >= v"0.4.0-dev+6256"
            f = f === :/           ? :(Base.RDivFun()) :
                f === :\           ? :(Base.LDivFun()) :
                f === :div         ? :(Base.IDivFun()) :
                f
        end
        if VERSION >= v"0.4.0-dev+6353"
            f = f === :$           ? :(Base.XorFun()) :
                f === :.+          ? :(Base.DotAddFun()) :
                f === :.-          ? :(Base.DotSubFun()) :
                f === :.*          ? :(Base.DotMulFun()) :
                f === :mod         ? :(Base.ModFun()) :
                f === :rem         ? :(Base.RemFun()) :
                # DotRemFun is defined, but ::call(::DotRemFun, ...) is not until later
                #f === :.%          ? :(Base.DotRemFun()) :
                f === :.<<         ? :(Base.DotLSFun()) :
                f === :.>>         ? :(Base.DotRSFun()) :
                f
        end
        if VERSION >= v"0.4.0-dev+6359"
            f = f === :./          ? :(Base.DotRDivFun()) :
                f
        end
        if VERSION >= v"0.4.0-rc1+59"
            f = f === :max         ? :(Base.ElementwiseMaxFun()) :
                f === :min         ? :(Base.ElementwiseMinFun()) :
                f
        end
        if VERSION >= v"0.5.0-dev+741"
            f = f === :complex     ? :(Base.SparseArrays.ComplexFun()) :
                f === :dot         ? :(Base.SparseArrays.DotFun()) :
                f
        end
        if VERSION >= v"0.5.0-dev+1472"
            f = f === Symbol(".÷") ? :(Base.DotIDivFun()) :
                f === :.%          ? :(Base.DotRemFun()) :
                f
        end
        f
    end
end

if !isdefined(Base, :Threads)
    @eval module Threads
        macro threads(expr)
            return esc(expr)
        end
        threadid() = 1
        nthreads() = 1
        export @threads, threadid, nthreads
    end
    export Threads
end

if !isdefined(Base, :normalize)
    function normalize!(v::AbstractVector, p::Real=2)
        nrm = norm(v, p)
        __normalize!(v, nrm)
    end

    @inline function __normalize!(v::AbstractVector, nrm::AbstractFloat)
        #The largest positive floating point number whose inverse is less than
        #infinity
        δ = inv(prevfloat(typemax(nrm)))
        if nrm ≥ δ #Safe to multiply with inverse
            invnrm = inv(nrm)
            scale!(v, invnrm)
        else # scale elements to avoid overflow
            εδ = eps(one(nrm))/δ
            scale!(v, εδ)
            scale!(v, inv(nrm*εδ))
        end
        v
    end

    copy_oftype{T,N}(A::AbstractArray{T,N}, ::Type{T}) = copy(A)
    copy_oftype{T,N,S}(A::AbstractArray{T,N}, ::Type{S}) = convert(AbstractArray{S,N}, A)

    function normalize(v::AbstractVector, p::Real = 2)
        nrm = norm(v, p)
        if !isempty(v)
            vv = copy_oftype(v, typeof(v[1]/nrm))
            return __normalize!(vv, nrm)
        else
            T = typeof(zero(eltype(v))/nrm)
            return T[]
        end
    end

    export normalize, normalize!
end

if !isdefined(Base, :AsyncCondition)
    type AsyncCondition
        cond::Condition
        handle::Ptr{Void}

        function AsyncCondition(func=nothing)
            this = new(Condition())
            # the callback is supposed to be called with this AsyncCondition
            # as the argument, so we need to create a wrapper callback
            if func == nothing
                function wrapfunc(data)
                    notify(this.cond)

                    nothing
                end
            else
                function wrapfunc(data)
                    notify(this.cond)
                    func(this)

                    nothing
                end
            end
            work = Base.SingleAsyncWork(wrapfunc)
            this.handle = work.handle

            this
        end
    end

    Base.wait(c::AsyncCondition) = wait(c.cond)
else
    import Base.AsyncCondition
end

# 0.5.0-dev+2301, JuliaLang/julia#14766
if !isdefined(Base, :unsafe_write)
    const unsafe_write = write
    export unsafe_write
end

# JuliaLang/julia#16219
if !isdefined(Base, @compat Symbol("@static"))
     macro static(ex)
        if isa(ex, Expr)
            if ex.head === :if
                cond = eval(current_module(), ex.args[1])
                if cond
                    return esc(ex.args[2])
                elseif length(ex.args) == 3
                    return esc(ex.args[3])
                else
                    return nothing
                end
            end
        end
        throw(ArgumentError("invalid @static macro"))
    end
    export @static
end

# JuliaLang/julia#14082
if VERSION < v"0.5.0-dev+4295"
    function repeat(A::AbstractArray;
                    inner=ntuple(x->1, ndims(A)),
                    outer=ntuple(x->1, ndims(A)))
        ndims_in = ndims(A)
        length_inner = length(inner)
        length_outer = length(outer)

        length_inner >= ndims_in || throw(ArgumentError("number of inner repetitions ($(length(inner))) cannot be less than number of dimensions of input ($(ndims(A)))"))
        length_outer >= ndims_in || throw(ArgumentError("number of outer repetitions ($(length(outer))) cannot be less than number of dimensions of input ($(ndims(A)))"))

        ndims_out = max(ndims_in, length_inner, length_outer)

        inner = vcat(collect(inner), ones(Int,ndims_out-length_inner))
        outer = vcat(collect(outer), ones(Int,ndims_out-length_outer))

        size_in = size(A)
        size_out = ntuple(i->inner[i]*size(A,i)*outer[i],ndims_out)::Dims
        inner_size_out = ntuple(i->inner[i]*size(A,i),ndims_out)::Dims

        indices_in = Array(Int, ndims_in)
        indices_out = Array(Int, ndims_out)

        length_out = prod(size_out)
        R = similar(A, size_out)

        for index_out in 1:length_out
            indices_out = ind2sub(size_out, index_out)
            for t in 1:ndims_in
                # "Project" outer repetitions into inner repetitions
                indices_in[t] = mod1(indices_out[t], inner_size_out[t])
                # Find inner repetitions using flooring division
                indices_in[t] = Base.fld1(indices_in[t], inner[t])
            end
            index_in = sub2ind(size_in, indices_in...)
            R[index_out] = A[index_in]
        end

        return R
    end
else
    const repeat = Base.repeat
end

if VERSION < v"0.5.0-dev+4267"
    if OS_NAME == :Windows
        const KERNEL = :NT
    else
        const KERNEL = OS_NAME
    end

    @eval is_apple()   = $(KERNEL == :Darwin)
    @eval is_linux()   = $(KERNEL == :Linux)
    @eval is_bsd()     = $(KERNEL in (:FreeBSD, :OpenBSD, :NetBSD, :Darwin, :Apple))
    @eval is_unix()    = $(is_linux() || is_bsd())
    @eval is_windows() = $(KERNEL == :NT)
    export is_apple, is_linux, is_bsd, is_unix, is_windows
else
    const KERNEL = Sys.KERNEL
end

import Base: srand, rand, rand!

if isdefined(Core, :String) && isdefined(Core, :AbstractString)
    # Not exported in order to not break code on 0.5
    const UTF8String = Core.String
    const ASCIIString = Core.String
else
    const String = Base.ByteString
    if VERSION >= v"0.4.0-dev+5243"
        @compat (::Type{Base.ByteString})(io::Base.AbstractIOBuffer) = bytestring(io)
    elseif VERSION >= v"0.4.0-dev+1246"
        @compat (::Type{Base.ByteString})(io::IOBuffer) = bytestring(io)
    end
    if VERSION >= v"0.4.0-dev+1246"
        @compat (::Type{Base.ByteString})(s::Cstring) = bytestring(s)
        @compat (::Type{Base.ByteString})(v::Vector{UInt8}) = bytestring(v)
        @compat (::Type{Base.ByteString})(p::Union{Ptr{Int8},Ptr{UInt8}}) = bytestring(p)
        @compat (::Type{Base.ByteString})(p::Union{Ptr{Int8},Ptr{UInt8}}, len::Integer) = bytestring(p, len)
        @compat (::Type{Base.ByteString})(s::AbstractString) = bytestring(s)
    end
end

if VERSION < v"0.5.0-dev+2285"
    fieldoffset(T, i) = @compat UInt(fieldoffsets(T)[i])
    export fieldoffset
end

if !isdefined(Base, :view)
    const view = slice
end

if !isdefined(Base, :pointer_to_string)

    function pointer_to_string(p::Ptr{UInt8}, len::Integer, own::Bool=false)
        a = ccall(:jl_ptr_to_array_1d, Vector{UInt8},
                  (Any, Ptr{UInt8}, Csize_t, Cint), Vector{UInt8}, p, len, own)
        ccall(:jl_array_to_string, String, (Any,), a)
    end

    pointer_to_string(p::Ptr{UInt8}, own::Bool=false) =
        pointer_to_string(p, ccall(:strlen, Csize_t, (Cstring,), p), own)

end

if VERSION < v"0.5.0-dev+4612"
    export unsafe_string, unsafe_wrap
    unsafe_wrap(::Type{Compat.String}, p::@compat(Union{Ptr{Int8},Ptr{UInt8}}), own::Bool=false) = pointer_to_string(convert(Ptr{UInt8}, p), own)
    unsafe_wrap(::Type{Compat.String}, p::@compat(Union{Ptr{Int8},Ptr{UInt8}}), len, own::Bool=false) = pointer_to_string(convert(Ptr{UInt8}, p), len, own)
    unsafe_wrap(::Type{Array}, p::Ptr, dims, own::Bool=false) = pointer_to_array(p, dims, own)
    unsafe_string(p::@compat(Union{Ptr{Int8},Ptr{UInt8}})) = bytestring(p)
    unsafe_string(p::@compat(Union{Ptr{Int8},Ptr{UInt8}}), len) = bytestring(p, len)
    if Cstring != Ptr{UInt8}
        unsafe_string(p::Cstring) = unsafe_string(Ptr{UInt8}(p))
    end
end

if !isdefined(Base, :allunique)
    allunique(itr) = length(itr) == length(unique(itr))
    export allunique
end

if !isdefined(Base, :promote_eltype_op)
    if isdefined(Base, :promote_op)
        promote_eltype_op(::Any) = Union{}
        promote_eltype_op{T}(op, ::AbstractArray{T}) = Base.promote_op(op, T)
        promote_eltype_op{T}(op, ::T               ) = Base.promote_op(op, T)
        promote_eltype_op{R,S}(op, ::AbstractArray{R}, ::AbstractArray{S}) = Base.promote_op(op, R, S)
        promote_eltype_op{R,S}(op, ::AbstractArray{R}, ::S) = Base.promote_op(op, R, S)
        promote_eltype_op{R,S}(op, ::R, ::AbstractArray{S}) = Base.promote_op(op, R, S)
        promote_eltype_op(op, A, B, C, D...) = Base.promote_op(op, eltype(A), promote_eltype_op(op, B, C, D...))
    else
        promote_eltype_op(::Any) = Union{}
        promote_eltype_op{T}(op, ::AbstractArray{T}) = T
        promote_eltype_op{T}(op, ::T               ) = T
        promote_eltype_op{R,S}(op, ::AbstractArray{R}, ::AbstractArray{S}) = Base.promote_type(R, S)
        promote_eltype_op{R,S}(op, ::AbstractArray{R}, ::S) = Base.promote_type(R, S)
        promote_eltype_op{R,S}(op, ::R, ::AbstractArray{S}) = Base.promote_type(R, S)
        promote_eltype_op(op, A, B, C, D...) = Base.promote_type(eltype(A), promote_eltype_op(op, B, C, D...))
    end
else
    import Base.promote_eltype_op
end

if VERSION < v"0.5.0-dev+4351"
    Base.join(io::IO, items) =
        print_joined(io, items)
    Base.join(io::IO, items, delim) =
        print_joined(io, items, delim)
    Base.join(io::IO, items, delim, last) =
        print_joined(io, items, delim, last)
    Base.unescape_string(io, s::AbstractString) =
        print_unescaped(io, s)
    Base.escape_string(io, str::AbstractString, esc::AbstractString) =
        print_escaped(io, str, esc)
end

if VERSION < v"0.5.0-dev+4340"
    Base.show(io::IO, mime, x) = writemime(io, mime, x)
end

if VERSION >= v"0.5.0-dev+4338"
    import Base.LinAlg.BLAS.@blasfunc
elseif VERSION >= v"0.5.0-dev+1871"
    import Base.@blasfunc
else
    macro blasfunc(x)
        Expr(:quote, Base.blasfunc(x))
    end
end

import Base: redirect_stdin, redirect_stdout, redirect_stderr
if VERSION < v"0.6.0-dev.374"
    for (F,S) in ((:redirect_stdin, :STDIN), (:redirect_stdout, :STDOUT), (:redirect_stderr, :STDERR))
        @eval function $F(f::Function, stream)
            STDOLD = $S
            $F(stream)
            try f() finally $F(STDOLD) end
        end
    end
end

if VERSION < v"0.6.0-dev.528"
    macro __DIR__()
        Base.source_dir()
    end
    export @__DIR__
end

# PR #17302
# Provide a non-deprecated version of `@vectorize_(1|2)arg` macro which defines
# deprecated version of the function so that the depwarns can be fixed without
# breaking users.
# Packages are expected to use this to maintain the old API until all users
# of the deprecated vectorized function have migrated.
# These macros should raise a depwarn when the `0.5` support is dropped from
# `Compat` and be dropped when the support for `0.6` is dropped from `Compat`.
# Modified based on the version copied from 0.6 Base.
macro dep_vectorize_1arg(S, f)
    S = esc(S)
    f = esc(f)
    T = esc(:T)
    x = esc(:x)
    AbsArr = esc(:AbstractArray)
    ## Depwarn to be enabled when 0.5 support is dropped.
    # depwarn("Implicit vectorized function is deprecated in favor of compact broadcast syntax.",
    #         Symbol("@dep_vectorize_1arg"))
    :(@deprecate $f{$T<:$S}($x::$AbsArr{$T}) @compat($f.($x)))
end

macro dep_vectorize_2arg(S, f)
    S = esc(S)
    f = esc(f)
    T1 = esc(:T1)
    T2 = esc(:T2)
    x = esc(:x)
    y = esc(:y)
    AbsArr = esc(:AbstractArray)
    ## Depwarn to be enabled when 0.5 support is dropped.
    # depwarn("Implicit vectorized function is deprecated in favor of compact broadcast syntax.",
    #         Symbol("@dep_vectorize_2arg"))
    quote
        @deprecate $f{$T1<:$S}($x::$S, $y::$AbsArr{$T1}) @compat($f.($x,$y))
        @deprecate $f{$T1<:$S}($x::$AbsArr{$T1}, $y::$S) @compat($f.($x,$y))
        @deprecate $f{$T1<:$S,$T2<:$S}($x::$AbsArr{$T1}, $y::$AbsArr{$T2}) @compat($f.($x,$y))
    end
end

if VERSION < v"0.5.0-dev+4677"
    using Base.LinAlg: HermOrSym
    Base.chol(A::HermOrSym) = Base.LinAlg.chol!(A.uplo == 'U' ? copy(A.data) : A.data')
    Base.cholfact(A::HermOrSym) = cholfact(A.data, Symbol(A.uplo))
    Base.cholfact!(A::HermOrSym) = cholfact!(A.data, Symbol(A.uplo))

    Base.cholfact(A::HermOrSym, T::Type) = cholfact(A.data, Symbol(A.uplo), T)
    Base.cholfact!(A::HermOrSym, T::Type) = cholfact!(A.data, Symbol(A.uplo), T)
end

# broadcast over same length tuples, from julia#16986
if v"0.5.0-dev+4002" ≤ VERSION < v"0.6.0-dev.693"
    Base.Broadcast.broadcast{N}(f, t::NTuple{N}, ts::Vararg{NTuple{N}}) = map(f, t, ts...)
end

# julia#18510
if VERSION < v"0.6.0-dev.826"
    _Nullable_field2(x) = !x
else
    _Nullable_field2(x) = x
end

# julia#18484
if VERSION < v"0.6.0-dev.848"
    unsafe_get(x::Nullable) = x.value
    unsafe_get(x) = x
    export unsafe_get
    Base.isnull(x) = false
end

# julia#18977
if !isdefined(Base, :xor)
    const xor = $
    const ⊻ = xor
    export xor, ⊻
end

# julia#19246
if !isdefined(Base, :numerator)
    const numerator = num
    const denominator = den
    export numerator, denominator
end

# julia #19950
if !isdefined(Base, :iszero)
    iszero(x) = x == zero(x)
    iszero(x::Number) = x == 0
    iszero(x::AbstractArray) = all(iszero, x)
    export iszero
end

# julia　#20407
if !isdefined(Base, :(>:))
    const >: = let
        _issupertype(a::ANY, b::ANY) = issubtype(b, a)
    end
    export >:
end

# julia#19088
if VERSION < v"0.6.0-dev.1256"
    Base.take!(io::Base.AbstractIOBuffer) = takebuf_array(io)
end

# julia #17323
if VERSION < v"0.5.0-dev+5380"
    export transcode
    transcode{T<:Compat.String}(::Type{T}, src::Union{Compat.String,Vector{UInt8}}) = utf8(src)
    transcode{T<:Compat.String}(::Type{T}, src::Vector{UInt16}) = utf8(utf16(src))
    transcode{T<:Compat.String}(::Type{T}, src::Vector{UInt32}) = utf8(utf32(src))
    transcode(::Type{UInt8}, src::Vector{UInt8}) = src
    transcode(::Type{UInt8}, src) = transcode(Compat.String, src).data
    function transcode(::Type{UInt16}, src::Union{Compat.String,Vector{UInt8}})
        d = utf16(utf8(src)).data
        return resize!(d, length(d)-1) # strip off trailing NUL codeunit
    end
    function transcode(::Type{UInt32}, src::Union{Compat.String,Vector{UInt8}})
        d = utf32(utf8(src)).data
        return resize!(d, length(d)-1) # strip off trailing NUL codeunit
    end
    transcode(::Type{UInt16}, src::Vector{UInt16}) = src
    transcode(::Type{UInt32}, src::Vector{UInt32}) = src
    if Cwchar_t == Int32
        transcode(::Type{Cwchar_t}, src::Vector{Cwchar_t}) = src
        transcode(::Type{Cwchar_t}, src) = reinterpret(Cwchar_t, transcode(UInt32, src))
        transcode(::Type{UInt8}, src::Vector{Cwchar_t}) = transcode(UInt8, reinterpret(UInt32, src))
        transcode(T, src::Vector{Cwchar_t}) = transcode(T, reinterpret(UInt32, src))
    end
end

# julia #17155 function composition and negation
if VERSION < v"0.6.0-dev.1883"
    export ∘
    ∘(f, g) = (x...)->f(g(x...))
    @compat Base.:!(f::Function) = (x...)->!f(x...)
end

if VERSION < v"0.5.0-dev+1450"
    Base.Diagonal(v::AbstractVector) = Diagonal(collect(v))
end
if VERSION < v"0.5.0-dev+3669"
    using Base: promote_op
    import Base: *
    eval(Expr(:typealias, :(SubMatrix{T}), :(SubArray{T,2})))
    (*)(A::SubMatrix, D::Diagonal) =
        scale!(similar(A, promote_op(*, eltype(A), eltype(D.diag))), A, D.diag)
    (*)(D::Diagonal, A::SubMatrix) =
        scale!(similar(A, promote_op(*, eltype(A), eltype(D.diag))), D.diag, A)
end

if VERSION >= v"0.5.0-dev+5509" && VERSION < v"0.6.0-dev.1632"
    # To work around unsupported syntax on Julia 0.4
    include_string("export .&, .|")
    include_string(".&(xs...) = broadcast(&, xs...)")
    include_string(".|(xs...) = broadcast(|, xs...)")
end

if VERSION < v"0.6.0-dev.2093" # Compat.isapprox to allow for NaNs
    using Base.rtoldefault
    function isapprox(x::Number, y::Number; rtol::Real=rtoldefault(x,y), atol::Real=0, nans::Bool=false)
        x == y || (isfinite(x) && isfinite(y) && abs(x-y) <= atol + rtol*max(abs(x), abs(y))) || (nans && isnan(x) && isnan(y))
    end
else
    import Base.isapprox
end

# julia #13998 single-argument min and max
if VERSION < v"0.5.0-dev+1279"
    Base.min(x::Real) = x
    Base.max(x::Real) = x
    Base.minmax(x::Real) = (x, x)
end

module TypeUtils
    using ..Compat: @static
    @static if isdefined(Base, :isabstract)
        using Base: isabstract, parameter_upper_bound, typename
    else
        isabstract(t::DataType) = t.abstract
        if isdefined(Base, :TypeConstructor)
            isabstract(t::TypeConstructor) = isabstract(t.body)
        end
        isabstract(t::ANY) = false
        parameter_upper_bound(t::DataType, idx) = t.parameters[idx].ub
        typename(t::DataType) = t.name
    end
    export isabstract, parameter_upper_bound, typename
end # module TypeUtils

# @view, @views, @__dot__
include("arraymacros.jl")

# julia #18839
if VERSION < v"0.6.0-dev.1024"
    @eval module Iterators
        export countfrom, cycle, drop, enumerate, flatten, product, repeated,
               rest, take, zip, partition

        import Base: eltype, start, next, done, length, size, ndims
        using Base: tuple_type_cons
        using Base: countfrom, cycle, drop, enumerate, repeated, rest, take,
                    zip
        using Compat

        # julia #14805
        if VERSION < v"0.5.0-dev+3256"
            immutable Flatten{I}
                it::I
            end

            flatten(itr) = Flatten(itr)

            eltype{I}(::Type{Flatten{I}}) = eltype(eltype(I))

            function start(f::Flatten)
                local inner, s2
                s = start(f.it)
                d = done(f.it, s)
                # this is a simple way to make this function type stable
                d && throw(ArgumentError("argument to Flatten must contain at least one iterator"))
                while !d
                    inner, s = next(f.it, s)
                    s2 = start(inner)
                    !done(inner, s2) && break
                    d = done(f.it, s)
                end
                return s, inner, s2
            end

            @inline function next(f::Flatten, state)
                s, inner, s2 = state
                val, s2 = next(inner, s2)
                while done(inner, s2) && !done(f.it, s)
                    inner, s = next(f.it, s)
                    s2 = start(inner)
                end
                return val, (s, inner, s2)
            end

            @inline function done(f::Flatten, state)
                s, inner, s2 = state
                return done(f.it, s) && done(inner, s2)
            end
        else
            using Base: flatten
        end

        # julia #14596
        if VERSION < v"0.5.0-dev+2305"
            # Product -- cartesian product of iterators
            @compat abstract type AbstractProdIterator end

            immutable Prod2{I1, I2} <: AbstractProdIterator
                a::I1
                b::I2
            end

            product(a) = Zip1(a)
            product(a, b) = Prod2(a, b)
            eltype{I1,I2}(::Type{Prod2{I1,I2}}) = Tuple{eltype(I1), eltype(I2)}
            length(p::AbstractProdIterator) = length(p.a)*length(p.b)

            function start(p::AbstractProdIterator)
                s1, s2 = start(p.a), start(p.b)
                s1, s2, Nullable{eltype(p.b)}(), (done(p.a,s1) || done(p.b,s2))
            end

            @inline function prod_next(p, st)
                s1, s2 = st[1], st[2]
                v1, s1 = next(p.a, s1)

                nv2 = st[3]
                if isnull(nv2)
                    v2, s2 = next(p.b, s2)
                else
                    v2 = nv2.value
                end

                if done(p.a, s1)
                    return (v1,v2), (start(p.a), s2, oftype(nv2,nothing), done(p.b,s2))
                end
                return (v1,v2), (s1, s2, Nullable(v2), false)
            end

            @inline next(p::Prod2, st) = prod_next(p, st)
            @inline done(p::AbstractProdIterator, st) = st[4]

            immutable Prod{I1, I2<:AbstractProdIterator} <: AbstractProdIterator
                a::I1
                b::I2
            end

            product(a, b, c...) = Prod(a, product(b, c...))
            eltype{I1,I2}(::Type{Prod{I1,I2}}) = tuple_type_cons(eltype(I1), eltype(I2))

            @inline function next{I1,I2}(p::Prod{I1,I2}, st)
                x = prod_next(p, st)
                ((x[1][1],x[1][2]...), x[2])
            end
        else
            using Base: product
        end

        # julia #15409
        if VERSION < v"0.5.0-dev+3510"
            partition{T}(c::T, n::Int) = PartitionIterator{T}(c, n)

            type PartitionIterator{T}
                c::T
                n::Int
            end

            eltype{T}(::Type{PartitionIterator{T}}) = Vector{eltype(T)}

            function length(itr::PartitionIterator)
                l = length(itr.c)
                return div(l, itr.n) + ((mod(l, itr.n) > 0) ? 1 : 0)
            end

            start(itr::PartitionIterator) = start(itr.c)

            done(itr::PartitionIterator, state) = done(itr.c, state)

            function next{T<:Vector}(itr::PartitionIterator{T}, state)
                l = state
                r = min(state + itr.n-1, length(itr.c))
                return sub(itr.c, l:r), r + 1
            end

            function next(itr::PartitionIterator, state)
                v = Vector{eltype(itr.c)}(itr.n)
                i = 0
                while !done(itr.c, state) && i < itr.n
                    i += 1
                    v[i], state = next(itr.c, state)
                end
                return resize!(v, i), state
            end
        else
            using Base: partition
        end
    end
else
    using Base: Iterators
end

if VERSION < v"0.6.0-dev.2840"
    export IndexStyle, IndexLinear, IndexCartesian
    eval(Expr(:typealias, :IndexStyle, :(Base.LinearIndexing)))
    eval(Expr(:typealias, :IndexLinear, :(Base.LinearFast)))
    eval(Expr(:typealias, :IndexCartesian, :(Base.LinearSlow)))
    IndexStyle{T}(::Type{T}) = Base.linearindexing(T)
    IndexStyle(args...) = Base.linearindexing(args...)
end

if VERSION < v"0.6.0-dev.1653"
    for (fname, felt) in ((:zeros,:zero), (:ones,:one))
        @eval begin
            # allow signature of similar
            Base.$fname(a::AbstractArray, T::Type, dims::Tuple) = fill!(similar(a, T, dims), $felt(T))
            Base.$fname(a::AbstractArray, T::Type, dims...) = fill!(similar(a,T,dims...), $felt(T))
            Base.$fname(a::AbstractArray, T::Type=eltype(a)) = fill!(similar(a,T), $felt(T))
        end
    end
end

# https://github.com/JuliaLang/julia/pull/20203
if VERSION < v"0.6.0-dev.2283"
    # not exported
    function readline(s::IO=STDIN; chomp::Bool=true)
        if chomp
            Base.chomp!(Base.readline(s))
        else
            Base.readline(s)
        end
    end
end

# https://github.com/JuliaLang/julia/pull/18727
if VERSION < v"0.6.0-dev.838"
    Base.convert{T}(::Type{Set{T}}, s::Set{T}) = s
    Base.convert{T}(::Type{Set{T}}, s::Set) = Set{T}(s)
end

include("to-be-deprecated.jl")

end # module Compat


text:
## --------  Unit conversions: courses

# Number of model courses per data course. Setting to 2 allows for half courses in the model, which helps with equally spaced grids.
const modelCoursesPerCourse = ncInt(2);

function model_to_data_courses(mCourses)
    return round.(mCourses ./ modelCoursesPerCourse; digits = 1)
end

# Always returns `ncInt`. Do not use to transform Floats (such as std errors).
function data_to_model_courses(dCourses :: I1) where I1 <: Number
    return round(ncInt, dCourses * modelCoursesPerCourse)
end

## ----------  Dollars

# Base year for prices
const baseYear = 2000
# Conversion factor for data dollars to model dollars
const dollarFactor = 1000.0
# Multiple a data dollar amount by this factor to make it annual
const dTimeFactors = Dict{Symbol,Float64}([
    :perYear => 1.0,
    :perMonth => 12.0,
    :perWeek => weeksPerYear,
    :perDay => daysPerYear,
    :perHour => hoursPerYear
]);

"""
    $(SIGNATURES)

Conversions of data to model dollars.
"""
function dollars_data_to_model(dDollars, timeUnit :: Symbol)
    return dDollars ./ dollarFactor .* dTimeFactors[timeUnit]
end

function dollars_model_to_data(mDollars, timeUnit :: Symbol)
    return mDollars .* dollarFactor ./ dTimeFactors[timeUnit]
end

# ----------------


text:
"""
    julia> cattuples((1,2),(3,4))
    (1, 2, 3, 4)
"""
# cattuples(a::Tuple,b::Tuple...) = (a..., cattuples(b...)...)
# cattuples(a::Tuple) = a

cattuples(x) = _cattuples(x...)
@inline _cattuples(x, y...) = (_cattuples(x)..., _cattuples(y)...)
@inline _cattuples(x, y) = (_cattuples(x)..., _cattuples(y)...)
@inline _cattuples(x::Tuple) = _cattuples(x...)
@inline _cattuples(x) = (x,)

"""
    julia> transposetuple(((1,2),(3,4),(5,6)))
    ((1, 3, 5), (2, 4, 6))
"""
function transposetuple(t)
    return Tuple((getindex.(t,j) for j in 1:length(t[1])))
end

allequal(a) = true
allequal(a,b) = a==b
allequal(a,b,c...) = (a==b) && allequal(b,c...)

export randnt

"""
    randnt(width::Int, depth::Int)

Construct a random nested NamedTuple with `Symbol`s at the leaves, width
`width`, and depth `depth`

EXAMPLE:
```
julia> randnt(3,2)
(a = (y = :y, r = :r, p = :p), v = (f = :f, e = :e, v = :v))
```
"""
function randnt(width, depth)
    k = unique(Symbol.(rand('a':'z', width)))
    if depth ≤ 1
        return namedtuple(k)(k)
    else
        nts = Tuple((randnt(width, depth-1) for _ in 1:length(k)))
        return namedtuple(k)(nts)
    end
end

export unwrap
unwrap(x) = x


text:
struct SemiconductorData <: AbstractMaterial
    material::SemiconductorMaterial
    Ec::Number
    Ev::Number
    φf::Number
    φfₚ::Number
    φfₑ::Number
    ψ::Number
    n::Number
    p::Number
    Na::Number
    Nd::Number
    T::Number
    dist::DistributionFunc
#=     μₙ::AbstractMobilityModel
    μₚ::AbstractMobilityModel =#
end

include("MobilityModels.jl")

SemiconductorData(S::SemiconductorData,ψ,φfₑ=0.0*Eunit,φfₚ=0.0*Eunit) = SemiconductorData(S.T,S.dist,S.material;psi=ψ,φfₑ=φfₑ,φfₚ=φfₚ)

function SemiconductorData(T,f::DistributionFunc,S::SemiconductorMaterial;psi=0.0*Eunit,φfₑ=0.0*Eunit,φfₚ=0.0*Eunit)
    E_f = Ef(T,f,S)
    _Ec = Ec(T,S)-E_f-psi
    _Ev = Ev(T,S)-E_f-psi
    _φf = Ei(T,f,S)-E_f
    _φfₑ = φfₑ
    _φfₚ = φfₚ
    _n = n(T,E_f,f,S)
    _p = p(T,E_f,f,S)
    _Na = Na(S.Dopping,T,E_f,0,Ev(T,S))
    _Nd = Nd(S.Dopping,T,E_f,Ec(T,S),Ev(T,S))
    return SemiconductorData(S,_Ec,_Ev,_φf,_φfₑ,_φfₚ,psi,_n,_p,_Na,_Nd,T,f)
end

ψ(S::SemiconductorData) = S.ψ
Eg(S::SemiconductorData) = Ec(S)-Ev(S)
Ec(S::SemiconductorData) = S.Ec
Ev(S::SemiconductorData) = S.Ev
Ef(S::SemiconductorData) = 0*Eunit
Ei(S::SemiconductorData) = S.φf

Na(S::SemiconductorData) = S.Na
Nd(S::SemiconductorData) = S.Nd
ni(S::SemiconductorData) = sqrt(S.n*S.p)
n(S::SemiconductorData,ψ=0.0*Eunit,φfₑ=0.0*Eunit) = n(S.T,ψ+φfₑ,Ec(S),S.dist,S.material)
p(S::SemiconductorData,ψ=0.0*Eunit,φfₚ=0.0*Eunit) = p(S.T,ψ+φfₚ,Ev(S),S.dist,S.material)
#n(S::SemiconductorData) = n(S,0.0*Eunit,φfₑ(S))
#p(S::SemiconductorData) = p(S,0.0*Eunit,φfₚ(S))
φfₚ(S::SemiconductorData) = S.φfₚ
φfₑ(S::SemiconductorData) = S.φfₑ
ϵᵣ(S::SemiconductorData) = S.material.ϵₘ
C(S::SemiconductorData) = Nd(S)-Na(S)
ρ(S::SemiconductorData) = q/ϵ₀/ϵᵣ(S) * (p(S) - n(S)+C(S))
ρ(S::SemiconductorData,ψ) = q/ϵ₀/ϵᵣ(S) * (p(S,ψ) - n(S,ψ)+C(S))
DistFun(S::SemiconductorData,E) = S.dist(E)
Temperature(S::SemiconductorData) = S.T
χ(S::SemiconductorData) = S.material.χ
Nv(S::SemiconductorData) = Nv(Temperature(S),S.material)
Nc(S::SemiconductorData) = Nc(Temperature(S),S.material)

function scalling(S::SemiconductorData,L)
    n_i = ni(S)
    ϕₜ = kb*Temperature(S)
    Vscal = ϕₜ
    Cscal = max(n_i,abs.(C(S)))
    lscal = L
    λ = ϵ₀*11.7*Vscal/(lscal^2*q*Cscal)
    δ = n_i/Cscal
    μscal = 1000
    Ld0 = sqrt(ϵ₀*11.7* ϕₜ/q/abs(Cscal))
    Rscal = ϕₜ*μscal*Cscal/lscal/lscal
    tscal = lscal^2/ϕₜ/μscal
    Jscal = q*Vscal*μscal*Cscal/lscal
    Dscal = μscal*Vscal
    return Vscal,lscal,Cscal,μscal,Dscal,Jscal,Rscal,tscal,δ,λ,Ld0,n_i

end
function Δχ(S::SemiconductorData,δχ)
    M = Δχ(S.material,δχ)
    SN = SemiconductorData(S.T,S.dist,M)
    S.material = M
    S.Eg = SN.Eg
    S.Ec = SN.Ec
    S.Ev = SN.Ev
    S.Ef = SN.Ef
    S.φf = SN.φf
    S.ψ = SN.ψ
    S.Ei = SN.Ei
    S.n = SN.n
    S.p = SN.p
    S.Na = SN.Na
    S.Nd = SN.Nd
    S.N = SN.N
    S.P = SN.P
    S.T = SN.T
end

function ξ(SD::SemiconductorData,psi)
    return Kingston(SD,psi,SD.dist)
end
function ξ(SD::SemiconductorData,psi,phi)
    return Kingston(SD,psi,phi,SD.dist)
end

function KingstonBKP(SD::SemiconductorData,psi::Number,f::BoltzmanDist)
    #convert to psi definition
    psi = psi +(Ef(SD)-Ei(SD))
    n_i = ni(SD)
    E_d = Ed(SD.material.Dopping)
    if Nd(SD.material.Dopping) != 0
        E_d = Ec(SD)-E_d
    else
        E_d = Ev(SD)+E_d
    end
    ϕt = kb*SD.T
    ωDI = (E_d-Ei(SD))/ϕt
    uc = Ec(SD)/ϕt
    ub = (Ef(SD)-Ei(SD))/ϕt
    us = (psi)/ϕt
    p = (ub-us)*sinh(ub)-(cosh(ub)-cosh(us))
    LD = sqrt(ϵ₀*SD.material.ϵₘ*kb*SD.T/2/q/n_i)
    return ϕt*sqrt(2)/LD*sqrt(abs(p))
end

#= function Kingston(SD::SemiconductorData,psi::Number,f::BoltzmanDist)
    #convert to psi definition
    u = psi
    n_i = ni(SD)
    phi = kb*SD.T
    epsilon = ϵ₀*11.7
    N_a = Na(SD)
    Ld = sqrt(2*phi*q*N_a/epsilon)
    F(x) = sqrt(abs(exp(-x/phi)+x/phi-1 + n_i^2/N_a^2*(exp(x/phi)-x/phi-1)))
    return  sign(psi)*Ld*F(u)#*sqrt(abs(q * (N_a ^ 2 * exp(-u / phi) * phi + exp(u / phi) * n_i ^ 2 * phi + (u - phi) * N_a ^ 2 - n_i ^ 2 * (u + phi)) / epsilon / N_a))
end =#

function KingstonBKP(SD::SemiconductorData,psi::Number,f::DistributionFunc)
    E_d = Ed(SD.material.Dopping)
    ϕt = kb*SD.T
    us = psi/ϕt
    si = 1.0
    if Nd(SD.material.Dopping) != 0
        E_d = Ec(SD)-E_d
        si=-1
    else
        E_d = Ev(SD)+E_d
    end
    Fnd(x) = -ϕt * log1p(2 * exp((Ef(SD) - Ec(SD) + x + E_d) / ϕt)) + ϕt * ((Ef(SD) - Ec(SD) + x + E_d) / ϕt)
    Fna(x) = ϕt * log1p(4 * exp((Ev(SD) - x + E_d - Ef(SD)) / ϕt)) - ϕt * ((Ev(SD) - x + E_d - Ef(SD)) / ϕt)
    p1 = Nd(SD.material.Dopping) * (Fnd(-psi)-Fnd(0))
    p2 = -Na(SD.material.Dopping) * (Fna(psi)-Fna(0))
    ηₚ = (Ev(SD)-Ef(SD))/ϕt
    Pp = ϕt*Nv(SD)*(Primitive(f,ηₚ-us)-Primitive(f,ηₚ))
    ηₙ = (Ef(SD)-Ec(SD))/ϕt
    Pn = -ϕt*Nc(SD)*(Primitive(f,ηₙ+us)-Primitive(f,ηₙ))
    K = q/ϵ₀/ϵᵣ(SD) *(p1+p2+Pp+Pn)
    return si*sign(psi)*sqrt(abs(K))
end

function Kingston(SD::SemiconductorData,psi0::Number,f::DistributionFunc)
    E_d = Ed(SD.material.Dopping)
    psi =psi0 + (Ef(SD)-Ei(SD))
    ϕt = kb*SD.T
    si = 1.0
    n_i = ni(SD)
    if Nd(SD.material.Dopping) != 0
        E_d = Ec(SD)-E_d
        si=-1
    else
        E_d = Ev(SD)+E_d
    end

    ωDI = (E_d-Ei(SD))/ϕt
    uc = Ec(SD)/ϕt
    ub = (Ef(SD)-Ei(SD))/ϕt
    us = (psi)/ϕt

    Fnd(x) = -ϕt * log1p(2 * exp((Ef(SD) - Ec(SD) + x + E_d) / ϕt)) + ϕt * ((Ef(SD) - Ec(SD) + x + E_d) / ϕt)
    Fna(x) = ϕt * log1p(4 * exp((Ev(SD) - x + E_d - Ef(SD)) / ϕt)) - ϕt * ((Ev(SD) - x + E_d - Ef(SD)) / ϕt)
    p1 = Nd(SD.material.Dopping)/n_i * log( (1 + 0.5*exp(ωDI - us))/(1 + 0.5*exp(ωDI - ub)) )
    p2 = Na(SD.material.Dopping)/n_i * log( (1 + 0.25*exp(us - ωDI))/(1 + 0.25*exp(ub - ωDI)) )
    ωVI =  (Ev(SD)-Ei(SD))/ϕt
    p3 = -1/(sqrt(pi)/2*f(ωVI))*(2/3*Primitive(f,ωVI-ub)-2/3*Primitive(f,ωVI-us))

    ωCI =  (Ec(SD)-Ei(SD))/ϕt
    p4 = 1/(sqrt(pi)/2*f(-ωCI))*(2/3*Primitive(f,us-ωCI)-2/3*Primitive(f,ub-ωCI))

    LD = sqrt(ϵ₀*SD.material.ϵₘ*kb*SD.T/2/q/n_i)

    return  sign(psi0)*ϕt/LD*sqrt(abs(p1+p2+p3+p4))

end

function Kingston(SD::SemiconductorData,psi0,phib,f::DistributionFunc)
    E_d = Ed(SD.material.Dopping)
    E_f = Ef(SD)
    psi = psi0 + (E_f-Ei(SD))
    ϕt = kb*SD.T
    si = 1.0
    n_i = ni(SD)
    if Nd(SD.material.Dopping) != 0
        E_d = Ec(SD)-E_d
        si=-1
        phibn = 0.0
        phibp = phib
    else
        E_d = Ev(SD)+E_d
        phibn = phib
        phibp = 0.0
    end

    ωDI = (E_d-Ei(SD))/ϕt
    uc = Ec(SD)/ϕt
    ub = (E_f-Ei(SD))/ϕt
    us = (psi)/ϕt

    Fnd(x) = -ϕt * log1p(2 * exp((Ef(SD) - Ec(SD) + x + E_d) / ϕt)) + ϕt * ((Ef(SD) - Ec(SD) + x + E_d) / ϕt)
    Fna(x) = ϕt * log1p(4 * exp((Ev(SD) - x + E_d - Ef(SD)) / ϕt)) - ϕt * ((Ev(SD) - x + E_d - Ef(SD)) / ϕt)
    p1 = Nd(SD.material.Dopping)/n_i * log( (1 + 0.5*exp(ωDI - us))/(1 + 0.5*exp(ωDI - ub)) )
    p2 = Na(SD.material.Dopping)/n_i * log( (1 + 0.25*exp(us - ωDI))/(1 + 0.25*exp(ub - ωDI)) )
    ωVI =  (Ev(SD)-Ei(SD))/ϕt
    p3 = -1/(sqrt(pi)/2*f(ωVI))*(2/3*Primitive(f,ωVI-ub+phibp/ϕt)-2/3*Primitive(f,ωVI-us+phibp/ϕt))

    ωCI =  (Ec(SD)-Ei(SD))/ϕt
    p4 = 1/(sqrt(pi)/2*f(-ωCI))*(2/3*Primitive(f,us-ωCI-phibn/ϕt)-2/3*Primitive(f,ub-ωCI-phibn/ϕt))

    LD = sqrt(ϵ₀*SD.material.ϵₘ*kb*SD.T/2/q/n_i)

    return  sign(psi0)*ϕt/LD*sqrt(abs(p1+p2+p3+p4))

end


text:
using CoupledDipoles, Revise
using Random;
Random.seed!(1135);
using Plots, PlotThemes, LaTeXStrings
theme(:vibrant)

#= Spherical Cloud and Plane Wave laser at z-direction =#
N, ρ = 2_000, 1.0
atoms = Atom(Sphere(), sphere_inputs(N, ρ)...)

s, Δ = 1e-5, 0.0
laser = Laser(PlaneWave3D([0,0,1]), s, Δ)
simulation = LinearOptics(Scalar(), atoms, laser)

ωₙ, Γₙ = get_spectrum(simulation)
modes = classify_modes(simulation)

### ------ Choose beter looking modes--------------
## I choose these modes because they got clear interpretation
loc_idx = modes.loc[end - 5]
sub_idx = modes.sub[4]
super_idx = modes.super[4]

## Small figures
fig_loc = plot_mode(simulation, loc_idx)
plot!(; title="(a) Loc", guidefont=25, title_position=:left)
fig_sub = plot_mode(simulation, sub_idx)
plot!(; title="(b) Sub", guidefont=25, title_position=:left)
fig_super = plot_mode(simulation, super_idx)
plot!(; title="(c) SR", guidefont=25, title_position=:left)

## Main Figure
fig_big = scatter(ωₙ[modes.sub], Γₙ[modes.sub]; markershape=:square, label="Sub")
scatter!(ωₙ[modes.loc], Γₙ[modes.loc]; markershape=:circle, label="Loc")
scatter!(ωₙ[modes.super], Γₙ[modes.super]; markershape=:utriangle, label="SR")

## Add black markers to indicate the modes shown
scatter!(
    [ωₙ[loc_idx]], [Γₙ[loc_idx]]; label="", color=:black, markershape=:square, markersize=7
)
scatter!(
    [ωₙ[sub_idx]], [Γₙ[sub_idx]]; label="", color=:black, markershape=:circle, markersize=7
)
scatter!(
    [ωₙ[super_idx]],
    [Γₙ[super_idx]];
    label="",
    color=:black,
    markershape=:utriangle,
    markersize=7,
)
plot!(;
    left_margin=5Plots.mm,
    right_margin=5Plots.mm,
    top_margin=5Plots.mm,
    bottom_margin=5Plots.mm,
    yscale=:log10,
    framestyle=:box,
    tickfont=20,
    guidefont=25,
    legendfont=17,
    xlabel=L"\omega_n/\Gamma",
    ylabel=L"\Gamma_n/\Gamma",
    legend=:bottomleft,
)

### --------Edges of Lorentz-Lorenz shift ---------
## Original formula of Lorentz-Lorenz shift cames from Eq 4 in :
## Microscopic and Macroscopic Signatures of 3D Anderson Localization of Light (2019)
## (using α = 0.5)
## α = 0.5
## λ_temp = 2π/k₀
## δ_cP = (λ_temp^3*ρ/(8*π^2)) + 0.5*sqrt(3*α*ρ*λ_temp^3/(4*π^2) -1)
## δ_cM = (λ_temp^3*ρ/(8*π^2)) - 0.5*sqrt(3*α*ρ*λ_temp^3/(4*π^2) -1)

k₀ = CoupledDipoles.k₀
## EPL Paper's formula
δ_cP = (π * ρ / (k₀^3)) + 0.5 * sqrt(3π * ρ / (k₀^3) - 1)
δ_cM = (π * ρ / (k₀^3)) - 0.5 * sqrt(3π * ρ / (k₀^3) - 1)

vline!([-δ_cP, -δ_cP]; label="", linestyle=:dashdot, color=:red, linewidth=5)
vline!([-δ_cM, -δ_cM]; label="", linestyle=:dashdot, color=:red, linewidth=5)

### ------------------ Figure 1--------------------
l = @layout([a{0.75w} [b; c; d]])
plot(fig_big, fig_super, fig_sub, fig_loc; layout=l, size=(1500, 900))
plot!(;
    left_margin=10Plots.mm,
    right_margin=5Plots.mm,
    top_margin=5Plots.mm,
    bottom_margin=5Plots.mm,
)

for fmt in ["png", "svg"]
    savefig("Fig1.$(fmt)")
end


text:
using SparseArrays: SparseMatrixCSC, sparse
using LinearAlgebra: mul!

"""
`struct SymmetricBasis`
"""
struct SymmetricBasis{BOP, PROP} <: ACEBasis
   pibasis::PIBasis{BOP}
   A2Bmap::SparseMatrixCSC{PROP, Int}
end

Base.length(basis::SymmetricBasis{BOP, PROP}) where {BOP, PROP} =
      size(basis.A2Bmap, 1)

fltype(basis::SymmetricBasis{BOP, PROP}) where {BOP, PROP} =  PROP
# rfltype(basis::SymmetricBasis) = rfltype(basis.pibasis)

SymmetricBasis(φ::AbstractProperty, args...; kwargs...) =
      SymmetricBasis(PIBasis(args...; kwargs...), φ)

function SymmetricBasis(pibasis, φ::TP) where {TP}

   # AA index -> AA spec
   AAspec = get_spec(pibasis)

   # construct the reverse mapping AA spec -> iAA and A spec -> iA
   invAAspec = Dict{Any, Int}()
   for (iAA, AA) in enumerate(AAspec)
      invAAspec[AA] = iAA
   end
   invAspec = Dict{Any, Int}()
   for (iA, A) in enumerate(get_spec(pibasis.basis1p))
      invAspec[A] = iA
   end

   # FUTURE: here we could analyze the different symbols and choose
   # which ones are of the "m"-type. e.g. there might be m1, m2 symbols. We
   # could create the convention that any symbol starting with m will be
   # treated like an m from the Ylm basis...

   # allocate the datastructure that computes and caches the
   # coupling coefficients
   # TODO: should this be stored with the basis?
   #       or maybe written to a file on disk? and then flushed every time
   #       we finish with a basis construction???
   #rotc = Rot3DCoeffs(rfltype(pibasis))
   rotc = rot3Dcoeffs(φ, rfltype(pibasis))
   # allocate triplet format
   Irow, Jcol, vals = Int[], Int[], TP[]
   # count the number of PI basis functions = number of rows
   idxB = 0

   # loop through AA basis, but skip most of them ...
   for (iAA, AA) in enumerate(AAspec)
      # AA = [b1, b2, ...], each bi = (n = ..., l = .., m = ...)
      # skip it unless all m are zero, because we want to consider each
      # (nn, ll, ...) block only once.
      # the loop over all possible `mm` must be taken care of inside
      # the `coupling_coeffs` implementation
      if !all(b.m == 0 for b in AA)
         continue
      end
      # get the coupling coefficients
      # here, we could help out a bit and make it easier to find the
      # the relevant basis functions?
      # AAcols will be in spec format i.e. named tuples
      U, AAcols = coupling_coeffs(AA, rotc, φ)

      # loop over the rows of U -> each specifies a basis function
      for irow = 1:size(U, 1)
         idxB += 1
         # loop over the columns of U / over brows
         for (icol, bcol) in enumerate(AAcols)
            # put bcol into the correct order
            bcol_ordered = _get_ordered(bcol, invAspec)
            # this is a subtle step: bcol and bcol_ordered are equivalent
            # permutation-invariant basis functions. This means we will
            # add the same PI basis function several times, but in the call to
            # `sparse` the values will just be added.
            if !haskey(invAAspec, bcol_ordered)
               @show bcol_ordered
               @show degree(bcol_ordered, NaiveTotalDegree(), pibasis.basis1p)
            end
            idxAA = invAAspec[bcol_ordered]
            push!(Irow, idxB)
            push!(Jcol, idxAA)
            push!(vals, TP(U[irow, icol]))
         end
      end
   end
   # TODO: filter and throw out everything that hasn't been used!!
   # create CSC: [   triplet    ]  nrows   ncols
   A2Bmap = sparse(Irow, Jcol, vals, idxB, length(AAspec))
   return SymmetricBasis(pibasis, A2Bmap)
end

function _get_ordered(bb, invAspec)
   iAs = [invAspec[b] for b in bb]
   return bb[ sortperm(iAs, rev = true) ]
end

function coupling_coeffs(bb, rotc::Rot3DCoeffs, φ::Invariant)
   # bb = [ b1, b2, b3, ...)
   # bi = (μ = ..., n = ..., l = ..., m = ...)
   #    (μ, n) -> n; only the l and m are used in the angular basis
   if length(bb) == 0
      return [1.0,], [bb,]
   end
   # convert to a format that the Rotations3D implementation can understand
   # this utility function splits the bb = (b1, b2 ...) with each
   # b1 = (μ = ..., n = ..., l = ..., m = ...) into
   #    l, and a new n = (μ, n)
   ll, nn = _b2llnn(bb)
   # now we can call the coupling coefficient construiction!!
   U, Ms = Rotations3D.rpi_basis(rotc, nn, ll)

   # but now we need to convert the m spec back to complete basis function
   # specifications
   rpibs = [ _nnllmm2b(nn, ll, mm) for mm in Ms ]

   return U, rpibs
end

function coupling_coeffs(bb, rotc::Rot3DCoeffsEquiv, φ::EuclideanVector)
   if length(bb) == 0
      error("an equivariant vector basis function cannot have length 0")
   end
   ll, nn = _b2llnn(bb)
   U, Ms = Rotations3D.vec3_symm_basis(rotc, nn, ll)
   rpibs = [ _nnllmm2b(nn, ll, mm) for mm in Ms ]
   return U, rpibs
end

function coupling_coeffs(bb, rotc::Rot3DCoeffs, φ::SphericalVector)
   if length(bb) == 0
      error("an equivariant vector basis function cannot have length 0")
   end
   ll, nn = _b2llnn(bb)
   # A small modification here - the function yvec_symm_basis shall
   # be φ related which specifies the blocks(the type of orbitals)...
   U, Ms = Rotations3D.yvec_symm_basis(rotc, nn, ll, getL(φ))
   rpibs = [ _nnllmm2b(nn, ll, mm) for mm in Ms ]
   return U, rpibs
end

function coupling_coeffs(bb, rotc::Rot3DCoeffs, φ::SphericalMatrix)
   if length(bb) == 0
      error("an equivariant matrix basis function cannot have length 0")
   end
   ll, nn = _b2llnn(bb)
   U, Ms = Rotations3D.mat_symm_basis(rotc, nn, ll, getL(φ)[1], getL(φ)[2])
   rpibs = [ _nnllmm2b(nn, ll, mm) for mm in Ms ]
   return U, rpibs
end

_nnllmm2b(nn, ll, mm) = [ _nlm2b(n, l, m) for (n, l, m) in zip(nn, ll, mm) ]

@generated function _nlm2b(n::NamedTuple{KEYS}, l, m) where {KEYS}
   code = "b = (" * prod("$(k) = n.$(k), " for k in KEYS) * "l = l, m = m )"
   :( $(Meta.parse(code)) )
end

function _b2llnn(bb)
   @assert all( iszero(b.m) for b in bb )
   ll = SVector( [b.l for b in bb]... )
   nn = SVector( [_all_but_lm(b) for b in bb]... )
   return ll, nn
end

"""
return a names tuple containing all values in b except those corresponding
to l and m keys
"""
@generated function _all_but_lm(b::NamedTuple{NAMES}) where {NAMES}
   code = "n = ("
   for k in NAMES
      if !(k in [:l, :m])
         code *= "$(k) = b.$(k), "
      end
   end
   code *= ")"
   quote
      $(Meta.parse(code))
      n
   end
end

# ---------------- Evaluation code

alloc_temp(basis::SymmetricBasis, args...) =
      (  AA = alloc_B(basis.pibasis),
         tmppi = alloc_temp(basis.pibasis) )

alloc_B(basis::SymmetricBasis, args...) =
      zeros(fltype(basis), length(basis))

function evaluate!(B, tmp, basis::SymmetricBasis,
                   cfg::AbstractConfiguration)
   # compute AA
   evaluate!(tmp.AA, tmp.tmppi, basis.pibasis, cfg)
   evaluate!(B, tmp, basis, tmp.AA)
   return B
end

# this function allows us to attach multiple symmetric bases to a single
#  𝑨 basis
#  TODO: this is extremely inefficient for multiple species But it is simple
#        and clean and will do for now...
function evaluate!(B, tmp, basis::SymmetricBasis,
                   AA::AbstractVector{<: Number})
   mul!(B, basis.A2Bmap, AA)
end

# ---- gradients

function gradtype(basis::SymmetricBasis)
   φ = zero(eltype(basis.A2Bmap))
   dAA = zero(gradtype(basis.pibasis))
   return typeof(φ * dAA)
end

alloc_temp_d(basis::SymmetricBasis, nmax::Integer) =
      (  AA = alloc_B(basis.pibasis),
         dAA = alloc_dB(basis.pibasis, nmax),
         tmppi = alloc_temp(basis.pibasis),
         tmpdpi = alloc_temp_d(basis.pibasis, nmax) )

alloc_dB(basis::SymmetricBasis, nmax::Integer) =
      zeros(gradtype(basis), length(basis), nmax)

function evaluate_d!(dB, tmpd, basis::SymmetricBasis,
                     cfg::AbstractConfiguration)
   # compute AA
   evaluate_ed!(tmpd.AA, tmpd.dAA, tmpd.tmpdpi, basis.pibasis, cfg)
   evaluate_d!(dB, tmpd, basis, tmpd.AA, tmpd.dAA)
   return dB
end

function evaluate_d!(dB, tmpd, basis::SymmetricBasis,
                     AA::AbstractVector{<: Number}, dAA)
   mul!(dB, basis.A2Bmap, dAA)
end


text:
############
# LeafNode #
############

function eachleaf!(f, node::AbstractLeafNode)
    @inbounds @simd for i in eachindex(node)
        if isactive(node, i)
            unsafe_setindex!(node, f(unsafe_getindex(node, i)), i)
        end
    end
end

@inline function eachleaf!(f, node::AbstractLeafNode, indices::CartesianIndices)
    @boundscheck checkbounds(node, indices)
    @inbounds @simd for cartesian in indices
        I = Tuple(cartesian)
        if isactive(node, I...)
            unsafe_setindex!(node, f(unsafe_getindex(node, I...)), I...)
        end
    end
end

#############################################
# Node/HashNode/DynamicNode/DynamicHashNode #
#############################################

@inline function _eachleaf!(f, node, i)
    @_propagate_inbounds_meta
    if isactive(node, i)
        eachleaf!(f, unsafe_getindex(node, i))
    end
end
function eachleaf!(f, node::Union{Node, HashNode, DynamicNode, DynamicHashNode})
    for i in eachindex(node)
        @inbounds _eachleaf!(f, node, i)
    end
end
function eachleaf_threads!(f, node::Union{Node, HashNode, DynamicNode, DynamicHashNode})
    Threads.@threads for i in eachindex(node)
        @inbounds _eachleaf!(f, node, i)
    end
end

@inline function _eachleaf!(f, node, I, indices)
    @_propagate_inbounds_meta
    dims = totalsize(Base.tail(TreeSize(node)))
    if isactive(node, I...)
        child = unsafe_getindex(node, I...)
        offset = @. (I - 1) * dims
        childinds = CartesianIndex(offset.+1):CartesianIndex(offset.+dims)
        eachleaf!(f, child, (indices ∩ childinds) .- CartesianIndex(offset))
    end
end
function eachleaf!(f, node::Union{Node, HashNode, DynamicNode, DynamicHashNode}, indices::CartesianIndices)
    @boundscheck checkbounds(TreeView(node), indices)
    start = offset_cartesian(node, Tuple(first(indices))...)
    stop = offset_cartesian(node, Tuple(last(indices))...)
    for cartesian in start:stop
        I = Tuple(cartesian)
        @inbounds _eachleaf!(f, node, I, indices)
    end
end
function eachleaf_threads!(f, node::Union{Node, HashNode, DynamicNode, DynamicHashNode}, indices::CartesianIndices) # threads version
    @boundscheck checkbounds(TreeView(node), indices)
    S = TreeSize(node)
    start = offset_cartesian(node, Tuple(first(indices))...)
    stop = offset_cartesian(node, Tuple(last(indices))...)
    Threads.@threads for cartesian in start:stop
        I = Tuple(cartesian)
        @inbounds _eachleaf!(f, node, I, indices)
    end
end

##############################
# TreeView/AbstractTreeArray #
##############################

function eachleaf!(f, A::TreeView)
    if length(A) > THREADS_THRESHOLD
        eachleaf_threads!(f, rootnode(A))
    else
        eachleaf!(f, rootnode(A))
    end
    A
end

_torange(i::Int) = i:i
_torange(i::Union{UnitRange, Colon}) = i
function eachleaf!(f, A::TreeView{<: Any, N}, I::Vararg{Union{Int, UnitRange, Colon}}) where {N}
    indices = CartesianIndices(to_indices(A, map(_torange, I)))
    @boundscheck checkbounds(A, indices)
    @inbounds if length(A) > THREADS_THRESHOLD
        eachleaf_threads!(f, rootnode(A), indices)
    else
        eachleaf!(f, rootnode(A), indices)
    end
    A
end

eachleaf!(f, A::AbstractTreeArray) = (eachleaf!(f, A.tree); A)
function eachleaf!(f, A::AbstractTreeArray, I...)
    indices = to_indices(A, I)
    @boundscheck checkbounds(A, indices...)
    @inbounds eachleaf!(f, A.tree, indices...)
    A
end


text:
using ADCME
using PyPlot

x0 = rand(100)
x0 = @. x0*0.4 + 0.3
x1 = collect(LinRange(0, 1, 100))
y0 = sin.(2π*x0)
w = Variable(fc_init([1, 20, 20, 20, 1]))
y = squeeze(fc(x0, [20, 20, 20, 1], w))
loss = sum((y - y0)^2)

sess = Session(); init(sess)
BFGS!(sess, loss)
y1 = run(sess, y)
plot(x0, y0, ".", label="Data")
x_dnn = run(sess,  squeeze(fc(x1, [20, 20, 20, 1], w)))
plot(x1, x_dnn,  "--", label="DNN Estimation")
legend()
w1 = run(sess, w)

##############################

μ = Variable(w1)
ρ = Variable(zeros(length(μ)))
σ = log(1+exp(ρ))

function likelihood(z)
    w = μ + σ * z
    y = squeeze(fc(x0, [20, 20, 20, 1], w))
    sum((y - y0)^2) - sum((w-μ)^2/(2σ^2)) + sum((w-w1)^2)
end

function inference(x)
    z = tf.random_normal((length(σ),), dtype=tf.float64)
    w = μ + σ * z
    y = squeeze(fc(x, [20, 20, 20, 1], w))|>squeeze
end

W = tf.random_normal((10, length(w)), dtype=tf.float64)
L = constant(0.0)
for i = 1:10
    global L += likelihood(W[i])
end

y2 = inference(x1)

opt = AdamOptimizer(0.01).minimize(L)
init(sess)
# run(sess, L)
losses = []
for i = 1:2000
    _, l = run(sess, [opt, L])
    push!(losses, l)
    @info i, l
end

Y = zeros(100, 1000)
for i = 1:1000
    Y[:,i] = run(sess, y2)
end

for i = 1:1000
    plot(x1, Y[:,i], "--", color="gray", alpha=0.5)
end
plot(x1, x_dnn, label="DNN Estimation")
plot(x0, y1, ".", label="Data")
legend()

##############################
# Naive Uncertainty Quantification
function inference_naive(x)
    z = tf.random_normal((length(w1),), dtype=tf.float64)
    w = w1 + log(2)*z
    y = squeeze(fc(x, [20, 20, 20, 1], w))|>squeeze
end
y3 = inference(x1)

Y = zeros(100, 1000)
for i = 1:1000
    Y[:,i] = run(sess, y3)
end

for i = 1:1000
    plot(x1, Y[:,i], "--", color="gray", alpha=0.5)
end
plot(x1, x_dnn, label="DNN Estimation")
plot(x0, y1, ".", label="Data")
legend()


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

name = "DexBuilder"
version = v"2.20.0"

# Collection of sources required to build DexBuilder
sources = [
    "https://golang.org/dl/go1.12.4.linux-amd64.tar.gz" =>
    "d7d1f1f88ddfe55840712dc1747f37a790cbcaa448f6c9cf51bbe10aa65442f5",

    "https://github.com/dexidp/dex/archive/v2.20.0.tar.gz" =>
    "a1da82ff82362ef8db020397a007d63505500546d2892b5469098c853abe486e",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
export GOPATH=$WORKSPACE/gopath
export GOCACHE=$WORKSPACE/gocache
mkdir -p $GOPATH/src/github.com/dexidp
mkdir -p $GOCACHE
ln -fs $WORKSPACE/srcdir/dex-2.20.0 $GOPATH/src/github.com/dexidp/dex
export PATH=$PATH:$WORKSPACE/srcdir/go/bin
go get github.com/bitnami/bcrypt-cli
cd $GOPATH/src/github.com/dexidp/dex
make
mkdir -p $prefix/bin
mkdir -p $prefix/lib
cp $GOPATH/src/github.com/dexidp/dex/bin/* $prefix/bin/
cp $GOPATH/bin/bcrypt-cli $prefix/bin/
tar -czvf $prefix/lib/webtemplates.tar.gz -C $GOPATH/src/github.com/dexidp/dex/web static templates themes
exit

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:musl),
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    ExecutableProduct(prefix, "dex", :dex),
    ExecutableProduct(prefix, "grpc-client", :grpcclient),
    ExecutableProduct(prefix, "bcrypt-cli", :bcryptcli)
]

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

]

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

max_stars_repo_name: tanmaykm/DexBuilder

max_issues_repo_name: tanmaykm/DexBuilder

max_forks_repo_name: tanmaykm/DexBuilder