text:
# Internal functions to create the variable DataFrame

function _result_dataframe_vars(variable::JuMP.Containers.DenseAxisArray)

    if length(axes(variable)) == 1
        result = Vector{Float64}(undef, length(first(variable.axes)))

        for t in variable.axes[1]
            result[t] = JuMP.value(variable[t])
        end

        return DataFrames.DataFrame(var = result)

    elseif length(axes(variable)) == 2

        result = Array{Float64, length(variable.axes)}(undef, length(variable.axes[2]), length(variable.axes[1]))
        names = Array{Symbol, 1}(undef, length(variable.axes[1]))

        for t in variable.axes[2], (ix, name) in enumerate(variable.axes[1])
            result[t, ix] = JuMP.value(variable[name, t])
            names[ix] = Symbol(name)
        end

        return DataFrames.DataFrame(result, names)

    elseif length(axes(variable)) == 3
        extra_dims = sum(length(axes(variable)[2:end-1]))
        extra_vars = [Symbol("S$(s)") for s in 1:extra_dims]
        result_df = DataFrames.DataFrame()
        names = vcat(extra_vars, Symbol.(axes(variable)[1]))

        for i in variable.axes[2]
            third_dim = collect(fill(i,size(variable)[end]))
            result = Array{Float64 ,2}(undef, length(last(variable.axes)),
                                              length(first(variable.axes)))
            for t in last(variable.axes), (ix, name) in enumerate(first(variable.axes))
                result[t, ix] = JuMP.value(variable[name, i, t])
            end
            res = DataFrames.DataFrame(hcat(third_dim, result))
            result_df = vcat(result_df, res)
        end

        return DataFrames.names!(result_df, names)

    else
        error("Dimension Number $(length(axes(variable))) not Supported")

    end

end

function _result_dataframe_duals(constraint::JuMP.Containers.DenseAxisArray)

    if length(axes(constraint)) == 1
        result = Vector{Float64}(undef, length(first(constraint.axes)))

        for t in constraint.axes[1]
            try result[t] = JuMP.dual(constraint[t])
            catch
                result[t] = NaN
            end
        end

        return DataFrames.DataFrame(var = result)

    elseif length(axes(constraint)) == 2

        result = Array{Float64, length(variable.axes)}(undef, length(constraint.axes[2]), length(constraint.axes[1]))
        names = Array{Symbol, 1}(undef, length(constraint.axes[1]))

        for t in constraint.axes[2], (ix, name) in enumerate(constraint.axes[1])
            try result[t, ix] = JuMP.dual(constraint[name, t])
            catch
                result[t, ix] =  NaN
            end
            names[ix] = Symbol(name)
        end

        return DataFrames.DataFrame(result, names)

    else
        error("Dimension Number $(length(axes(constraint))) not Supported")

    end

end

# Function to write results dataframes and variables to a dictionary

function get_model_result(op_m::OperationModel)

    results_dict = Dict{Symbol, DataFrames.DataFrame}()

    for (k, v) in vars(op_m.canonical)
        results_dict[k] = _result_dataframe_vars(v)
    end

    return results_dict

end

function get_model_duals(op_m::OperationModel, cons::Vector{Symbol})

    results_dict = Dict{Symbol, DataFrames.DataFrame}()

    for c in cons
        v = con(op_m.canonical, c)
        results_dict[c] = _result_dataframe_duals(v)
    end

    return results_dict

end

# Function to create a dictionary for the optimizer log of the simulation

function get_optimizer_log(op_m::OperationModel)

    canonical_model = op_m.canonical

    optimizer_log = Dict{Symbol, Any}()

    optimizer_log[:obj_value] = JuMP.objective_value(canonical_model.JuMPmodel)
    optimizer_log[:termination_status] = JuMP.termination_status(canonical_model.JuMPmodel)
    optimizer_log[:primal_status] = JuMP.primal_status(canonical_model.JuMPmodel)
    optimizer_log[:dual_status] = JuMP.dual_status(canonical_model.JuMPmodel)
    optimizer_log[:solver] =  JuMP.solver_name(canonical_model.JuMPmodel)
    try
        optimizer_log[:solve_time] = MOI.get(canonical_model.JuMPmodel, MOI.SolveTime())
    catch
        @warn("SolveTime() property not supported by $(optimizer_log[:solver])")
        optimizer_log[:solve_time] = "Not Supported by $(optimizer_log[:solver])"
    end
    return optimizer_log
end

# Function to create a dictionary for the time series of the simulation

function get_time_stamp(op_model::OperationModel)

    initial_time = PSY.get_forecast_initial_times(op_model.sys)[1]
    interval = PSY.get_forecasts_resolution(op_model.sys)
    horizon = PSY.get_forecasts_horizon(op_model.sys)
    range = collect(initial_time:interval:initial_time+ interval.*horizon)
    time_stamp = DataFrames.DataFrame(Range = range[:,1])

    return time_stamp
end

function get_time_stamp(stage::_Stage, start_time::Dates.DateTime)

    interval = PSY.get_forecasts_resolution(stage.sys)
    horizon = PSY.get_forecasts_horizon(stage.sys)
    range = collect(start_time:interval:start_time+ interval.*horizon)
    time_stamp = DataFrames.DataFrame(Range = range[:,1])

    return time_stamp
end


text:
struct Pricer <: AbstractPricer
    pricing_function::Union{Function, OptimizingIR.BasicBlockInterpreter}
    input_riskfactors::OptimizingIR.LookupTable{RiskFactor}
    currency::Currencies.Currency
end

struct CashflowType
    event::Symbol
    maturity::Date
end

struct CashflowPricer <: AbstractCashflowPricer
    pricing_function::Union{Function, OptimizingIR.BasicBlockInterpreter}
    input_riskfactors::OptimizingIR.LookupTable{RiskFactor}
    currency::Currencies.Currency
    price_output_index::Int # In the resulting tuple from the pricing function, identifies the index of the pricing information
    output_index_to_cashflow_type::Dict{Int, CashflowType} # Key: index in the resulting tuple from the pricing function for the cashflow
end

riskfactors(p::Pricer) = p.input_riskfactors
riskfactors(p::CashflowPricer) = p.input_riskfactors

get_functional_currency(p::Pricer) = p.currency
get_functional_currency(p::CashflowPricer) = p.currency

# maps (observable, state) -> variable
# where state is the state of nature of the pricing routine
# (state argument in lower! method)
const ObservablesBuffer = Dict{Tuple{BufferedObservable, Any}, OptimizingIR.ImmutableValue}

struct CompilerContext{M<:PricingModel, D<:MarketData.AbstractMarketDataProvider, IR<:OptimizingIR.Program, P<:Union{AbstractPricer, AbstractCashflowPricer}}
    model::M
    provider::D
    attr::ContractAttributes
    pricing_date::Date
    program::IR
    observables_buffer::ObservablesBuffer
    input_riskfactors_variable::OptimizingIR.ImmutableVariable # identifies the index of the argument for the pricing function which is a Vector of Risk Factor values
    input_riskfactors::OptimizingIR.LookupTable{RiskFactor}    # identifies the risk factor for each index of input_riskfactors_variable
    target_pricer_type::Type{P}
    output_var_to_cashflowtype::Union{Nothing, Dict{OptimizingIR.ImmutableVariable, CashflowType}}
end

struct CompilerResult{IR<:OptimizingIR.Program, P<:Union{AbstractPricer, AbstractCashflowPricer}}
    program::IR
    input_riskfactors::OptimizingIR.LookupTable{RiskFactor}
    currency::Currencies.Currency
    target_pricer_type::Type{P}
    price_output_index::Int # In the resulting tuple from the pricing function, identifies the index of the pricing information
    output_index_to_cashflow_type::Union{Nothing, Dict{Int, CashflowType}} # Key: index in the resulting tuple from the pricing function for the cashflow
end

get_market_data_provider(ctx::CompilerContext) = ctx.provider
get_pricing_model(ctx::CompilerContext) = ctx.model
get_pricing_date(ctx::CompilerContext) = ctx.pricing_date

@generated function get_target_pricer_type(::Type{CompilerContext{M,D,IR,P}}) where {M,D,IR,P}
    @assert P<:Union{AbstractPricer, AbstractCashflowPricer}
    P
end

"""
Returns `true` if `CompilerContext` was asked to create an `AbstractCashflowPricer`.
"""
@generated function is_cashflow_pricer(::Type{T}) :: Bool where {T<:CompilerContext}
    return get_target_pricer_type(T) == AbstractCashflowPricer
end

"""
Returns the initial state of nature for the current pricing model.

The entry point for contract compilation is
`lower!(ctx, contract)`, which is implemented as:

```julia
lower!(ctx::CompilerContext, c::Contract) = lower!(ctx, c, initial_state(ctx))
```

`Observables` are always valued in a given state by `lower!(ctx, obs, state)`,
therefore there's no need to implement `lower!` for observables with no `state` argument.
"""
initial_state(ctx::CompilerContext) = error("`initial_state` must be implemented for each `PricingModel`.")

ext: jl


text:
# one-column version
function op_H( LF, Vpot, v::Array{Float64,1} )
  # Is this OK?
  Hv = -0.5*op_nabla2( LF, v ) + Vpot .* v
end

# one-column version, sparse
function op_H( LF, ∇2::SparseMatrixCSC{Float64,Int64}, Vpot, v::Array{Float64,1} )
  # Is this OK?
  Hv = -0.5 * ∇2 * v + Vpot .* v
end

# multi-column version, sparse
function op_H( LF, ∇2::SparseMatrixCSC{Float64,Int64}, Vpot, v::Array{Float64,2} )
  Hv = -0.5 * ∇2 * v + Vpot .* v
end

# multicolumn version
function H_psi( LF, Vpot, psi::Array{Float64,2} )
  Ncol = size(psi)[2]
  Hpsi = zeros( Float64, Npoints, Ncol )
  for ic = 1:Ncol
    Hpsi[:,ic] = -0.5*op_nabla2( LF, psi[:,ic] ) + Vpot .* psi[:,ic]
  end
  return H_psi
end


text:
"""
    DiscreteBoundaryFunction(func, parameters)

A wrapper for boundary condition functions with optional parameters.
When `parameters=nothing`, the boundary condition `func` is called with the signature

    `func(i, j, grid, clock, model_fields)`

where `i, j` are the indices along the boundary,
where `grid` is `model.grid`, `clock.time` is the current simulation time and
`clock.iteration` is the current model iteration, and
`model_fields` is a `NamedTuple` with `u, v, w`, the fields in `model.tracers`,
and the fields in `model.diffusivities`, each of which is an `OffsetArray`s (or `NamedTuple`s
of `OffsetArray`s depending on the turbulence closure) of field data.

When `parameters` is not `nothing`, the boundary condition `func` is called with
the signature

    `func(i, j, grid, clock, model_fields)`

*Note* that the index `end` does *not* access the final physical grid point of
a model field in any direction. The final grid point must be explictly specified, as
in `model_fields.u[i, j, grid.Nz]`.
"""
struct DiscreteBoundaryFunction{P, F} <: Function
    func :: F
    parameters :: P
end

# Un-parameterized
@inline (bc::DiscreteBoundaryFunction{<:Nothing})(i, j, grid, clock, model_fields) =
    bc.func(i, j, grid, clock, model_fields)

# Parameterized
@inline (bc::DiscreteBoundaryFunction)(i, j, grid, clock, model_fields) =
    bc.func(i, j, grid, clock, model_fields, bc.parameters)

# Don't re-convert DiscreteBoundaryFunctions passed to BoundaryCondition constructor
BoundaryCondition(TBC, condition::DiscreteBoundaryFunction) =
    BoundaryCondition{TBC, typeof(condition)}(condition)


text:
function push_pair!(i,j,d2,matrix)
  d = sqrt(d2)
  matrix[i,j] = d
  return matrix
end

function add_inf!(matrix)
  for i in eachindex(matrix)
    if matrix[i] == 0
      matrix[i]=+Inf
    end
  end
end

function find_data!(matrix, dist, nframe)
  for i in 1:length(matrix[1,:])            # loop though number of solvents
    data = matrix[:,i]                      # all distances between protein atoms and the ith solvent
    dist[nframe,i] = minimum(data)
  end
end

function autocorr_cell(trajectory::Trajectory)

  nprot = length(trajectory.x_solute)  # number of protein atoms
  nsvt  = length(trajectory.x_solvent) # number of solvent atoms
  nframes = trajectory.nframes         # number of frames

  # vector with time - ns
  delta = 0.01
  time = zeros(nframes)
  t1 = 0.
  for i in 1:nframes
    t1 = t1 + delta
    time[i] = t1
  end

  Dist      = zeros(Float64, nframes, nsvt)        # matrix of distances between the atoms
  domain    = zeros(Float64, nframes, nsvt)        # matrix of correlation obtained from the evaluation of the particle position
  evals     = zeros(Float64, nframes, nsvt)        # Matrix of the number of events observed
  sp        = zeros(Float64, nframes, nsvt)        # Matrix of the sample space - all possible events
  cutoff    = 30 # celllistmap parameter

  #if nprot > nsvt
  #  matrix    = zeros(Float64,nsvt, nprot)
  #else
    matrix    = zeros(Float64,nprot, nsvt)
  #end

  for iframe in 1:nframes
    nextframe!(trajectory) # reading coordinates of next frame
    solute    = trajectory.solute        #
    solvent   = trajectory.solvent       # variables to compute the autocorrelation function
    x_solute  = trajectory.x_solute      #
    x_solvent = trajectory.x_solvent     #

    # CellLisMap parameters
    sides = getsides(trajectory,iframe)   # box
    box = Box(sides,cutoff)

    # Initialize auxiliary linked lists (largest set!)
    cl = CellList(x_solute, x_solvent, box)

    # atoms dist - must be allocated and calculated for each frame
    map_pairwise!((x,y,i,j,d2,output) -> push_pair!(i, j, d2, output), matrix, box, cl)

    add_inf!(matrix)

   # if nprot > nsvt
   #   new = transpose(matrix) # due to the CellListMap package, this must be performed to avoid boundserror
   #   find_data!(new, Dist, iframe)
   # else
      find_data!(matrix, Dist, iframe)
   # end

  end
  closetraj(trajectory)
  return Dist, domain, evals, sp, nprot, nsvt, nframes, time
end


text:
module HumanMotionPrediction

# package code goes here
include("blender.jl")

end # module


text:
@testset "Graphs" begin
    # From MATPOWERM files
    path = "./test/data/case6ww.m"
    network = PowerFlowNetwork(path, "MATPOWERM")
    g = SimpleGraph(network)
    @test nv(g) == nbus(network)
    @test ne(g) == nbranch(network)
    edges = [(1, 2), (1, 4), (1, 5), (2, 3), (2, 4), (2, 5), (2, 6), (3, 5), (3, 6), (4, 5), (5, 6)]
    @test ne(g) == length(edges)
    for (src, dst) in edges
        @test has_edge(g, src, dst)
    end

    # From RAWGO files (with non continuous index)
    path = "./test/data/C2S6N02045/scenario_1/"
    network = PowerFlowNetwork(path, "RAWGO")
    g = SimpleGraph(network)
    @test nv(g) == nbus(network)
    @test ne(g) == nbranch(network; distinct_pair=true)
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
"""
This package is intended to provide an abstract interface for working
with tree structures.
Though the package itself is not particularly sophisticated, it defines
the interface that can be used by other packages to talk about trees.
"""
module AbstractTrees

export print_tree, TreeCharSet, TreeIterator, Leaves, PostOrderDFS, Tree,
    AnnotationNode, StatelessBFS, treemap, treemap!, PreOrderDFS,
    ShadowTree, children

import Base: getindex, setindex!, iterate, nextind, print, show,
    eltype, IteratorSize, IteratorEltype, length, push!, pop!
using Base: SizeUnknown, EltypeUnknown
using Markdown

abstract type AbstractShadowTree end

"""
    children(x)

Return the immediate children of node `x`. You should specialize this method
for custom tree structures. It should return an iterable object for which an
appropriate implementation of `Base.pairs` is available.

The default behavior is to assume that if an object is iterable, iterating over
it gives its children. If an object is not iterable, assume it does not have
children.

# Example

```
struct MyNode{T}
    data::T
    children::Vector{MyNode{T}}
end
AbstractTrees.children(node::MyNode) = node.children
```
"""
children(x) = Base.isiterable(typeof(x)) ? x : ()

has_children(x) = children(x) !== ()

include("traits.jl")
include("implicitstacks.jl")
include("printing.jl")
include("indexing.jl")
include("iteration.jl")

## Special cases

# Types which are iterable but shouldn't be considered tree-iterable
children(x::Number) = ()
children(x::Char) = ()
children(x::Task) = ()
children(x::AbstractString) = ()

# Define this here, there isn't really a good canonical package to define this
# elsewhere
children(x::Expr) = x.args

# For AbstractDicts

printnode(io::IO, kv::Pair{K,V}) where {K,V} = printnode(io,kv[1])
children(kv::Pair{K,V}) where {K,V} = (kv[2],)

# For potentially-large containers, just show the type

printnode(io::IO, ::T) where T <: Union{AbstractArray, AbstractDict} = print(io, T)

end # module


text:
using FewShotAnomalyDetection
using UCI
using Flux
using MLDataPattern

hiddenDim = 32
latentDim = 3
numLayers = 3
nonlinearity = "relu"
layerType = "Dense"
β = 0.1 # ratio between reconstruction error and the distance between p(z) and q(z)
α = 0.1 # threshold in the memory that does not matter to us at the moment!
loss_α = 0.1 # importance ratio between anomalies and normal data in mem_loss
memorySize = 256
k = 256
labelCount = 1

batchSize = 100
numBatches = 1000

dataset = "breast-cancer-wisconsin" # name of the UCI dataset. You can find the names in the e.g. in the LODA paper http://agents.fel.cvut.cz/stegodata/pdfs/Pev15-Loda.pdf - last page
# dataset = "statlog-vehicle"

data, normal_labels, anomaly_labels = UCI.get_umap_data(dataset)
# This method iterates over all subdataset in a multiclass dataset or it only runs once for a single class dataset
subdatasets = UCI.create_multiclass(data, normal_labels, anomaly_labels)
subdata, class_label = subdatasets[1]
println(dataset*" "*class_label)
_X_tr, _y_tr, _X_tst, _y_tst = UCI.split_data(subdata, 0.8) #train/test split

# dirty data normalization
_X_tr .-= minimum(_X_tr)
_X_tr ./= maximum(_X_tr)
_X_tst .-= minimum(_X_tst)
_X_tst ./= maximum(_X_tst)

# no need to have it this way but we are used to it from legacy code
train = (_X_tr, _y_tr)
test = (_X_tst, _y_tst)

##############################################
# SVAE itself
##############################################
inputdim = size(train[1], 1)
svae = SVAEbase(inputdim, hiddenDim, latentDim, numLayers, nonlinearity, layerType)
mem = KNNmemory{Float32}(memorySize, inputdim, k, labelCount, (x) -> zparams(svae, x)[1], α)

# Basic Wasserstein loss to train the svae on unlabelled data
trainRepresentation(data) = wloss(svae, data, β, (x, y) -> FewShotAnomalyDetection.mmd_imq(x, y, 1))
# inserts the data into the memory
remember(data, labels) = trainQuery!(mem, data, labels)
# Expects anomalies in the data with correct label (some of them
trainWithAnomalies(data, labels) = FewShotAnomalyDetection.mem_wloss(svae, mem, data, labels, β, (x, y) -> FewShotAnomalyDetection.mmd_imq(x, y, 1), loss_α)

# Unsupervised learning
opt = Flux.Optimise.ADAM(1e-5)
cb = Flux.throttle(() -> println("SVAE: $(trainRepresentation(train[1]))"), 5)
Flux.train!((x) -> trainRepresentation(getobs(x)), Flux.params(svae), RandomBatches((train[1],), size = batchSize, count = numBatches), opt, cb = cb)
println("Train err: $(trainRepresentation(train[1])) vs test error: $(trainRepresentation(test[1]))")

# Adding stuff into the memory
remember(train[1], train[2])
remember(test[1][:, 230:end], test[2][230:end]) # we don't have anomalies in train for UCI datasets :( and let's add just a couple so we don't have to wait so long

numBatches = 1000 # it will take a looong time

# learn with labels
cb = Flux.throttle(() -> println("SVAE mem loss: $(trainWithAnomalies(test[1], test[2]))"), 60)
Flux.train!(trainWithAnomalies, Flux.params(svae), RandomBatches((test[1], test[2]), size = batchSize, count = numBatches), opt, cb = cb)


text:
include("printf.jl")

const _formatters = Dict{ASCIIStr,FmtSpec}()

function _get_formatter(fmt)
    global _formatters

    chkfmt = get(_formatters, fmt, nothing)
    chkfmt === nothing || return chkfmt
    _formatters[fmt] = FmtSpec(fmt)
end

_cfmt_comma(fspec::FmtSpec, x) = addcommasreal(_cfmt(fspec, x))
_cfmt_comma(fspec::FmtSpec{FmtStr}, x::Rational) =  addcommasrat(_cfmt(fspec, x))
_cfmt_comma(fspec::FmtSpec{<:FmtInts}, x) = checkcommas(_cfmt(fspec, x))

function _cfmt(fspec::FmtSpec, x)
    sv = Base.StringVector(23) # Trust that lower level code will expand if necessary
    pos = _fmt(sv, 1, fspec, x)
    resize!(sv, pos - 1)
    String(sv)
end

cfmt(fspec::FmtSpec, x) = fspec.tsep == 0 ? _cfmt(fspec, x) : _cfmt_comma(fspec, x)
cfmt(fmtstr::ASCIIStr, x) = cfmt(_get_formatter(fmtstr), x)

function generate_formatter(fmt::ASCIIStr)
    fspec = _get_formatter(fmt)
    fspec.tsep == 0 ? x -> _cfmt(fspec, x) : x -> _cfmt_comma(fspec, x)
end

function addcommasreal(s)
    len = length(s)
    dpos = findfirst( isequal('.'), s )
    dpos !== nothing && return addcommas(s, len, dpos-1)
    # find the rightmost digit
    for i in len:-1:1
        isdigit( s[i] ) && return addcommas(s, len, i)
    end
    s
end

# commas are added to only the numerator
addcommasrat(s) = addcommas(s, length(s), findfirst( isequal('/'), s )-1)

function checkcommas(s)
    len = length(s)
    for i in len:-1:1
        isdigit( s[i] ) && return addcommas(s, len, i)
    end
    s
end

function addcommas(s::T, len, lst) where {T<:AbstractString}
    lst < 4 && return s
    beg = 1
    while beg < len
        isdigit(s[beg]) && break
        beg += 1
    end
    dig = lst - beg + 1
    dig < 4 && return s

    commas = div(dig - 1, 3)
    sv = Base.StringVector(len + commas)

    for i = 1:beg-1; sv[i] = s[i]; end
    cnt = dig - commas*3
    pos = beg - 1
    for i = beg:lst-3
        sv[pos += 1] = s[i]
        (cnt -= 1) == 0 && (cnt = 3; sv[pos += 1] = ',')
    end
    for i = lst-2:len; sv[i+commas] = s[i]; end
    T(sv)
end
addcommas(s) = (l = length(s); addcommas(s, l, l))

function generate_format_string(;
                                width::Int=-1,
                                precision::Int= -1,
                                leftjustified::Bool=false,
                                zeropadding::Bool=false,
                                commas::Bool=false,
                                signed::Bool=false,
                                positivespace::Bool=false,
                                alternative::Bool=false,
                                conversion::ASCIIStr="f" #aAdecEfFiosxX
                                )
    s = ['%'%UInt8]
    commas &&
        push!(s, '\'')
    alternative && in( conversion[1], "aAeEfFoxX" ) &&
        push!(s, '#')
    zeropadding && !leftjustified && width != -1 &&
        push!(s, '0')
    if signed
        push!(s, '+')
    elseif positivespace
        push!(s, ' ')
    end
    if width != -1
        leftjustified && push!(s, '-')
        append!(s, _codeunits(string( width )))
    end
    precision != -1 &&
        append!(s, _codeunits(string( '.', precision )))
    String(append!(s, _codeunits(conversion)))
end

function format( x::T;
                 width::Int=-1,
                 precision::Int= -1,
                 leftjustified::Bool=false,
                 zeropadding::Bool=false, # when right-justified, use 0 instead of space to fill
                 commas::Bool=false,
                 signed::Bool=false, # +/- prefix
                 positivespace::Bool=false,
                 stripzeros::Bool=(precision== -1),
                 parens::Bool=false, # use (1.00) instead of -1.00. Used in finance
                 alternative::Bool=false, # usually for hex
                 mixedfraction::Bool=false,
                 mixedfractionsep::UTF8Str="_",
                 fractionsep::UTF8Str="/", # num / den
                 fractionwidth::Int = 0,
                 tryden::Int = 0, # if 2 or higher,
                                  # try to use this denominator, without losing precision
                 suffix::UTF8Str="", # useful for units/%
                 autoscale::Symbol=:none, # :metric, :binary or :finance
                 conversion::ASCIIStr=""
                 ) where {T<:Real}
    checkwidth = commas
    if conversion == ""
        if T <: AbstractFloat || T <: Rational && precision != -1
            actualconv = "f"
        elseif T <: Unsigned
            actualconv = "x"
        elseif T <: Integer
            actualconv = "d"
        else
            conversion = "s"
            actualconv = "s"
        end
    else
        actualconv = conversion
    end
    signed && commas && error( "You cannot use signed (+/-) AND commas at the same time")

    T <: Rational && conversion == "s" && (stripzeros = false)
    if ( T <: AbstractFloat && actualconv == "f" || T <: Integer ) && autoscale != :none
        actualconv = "f"
        if autoscale == :metric
            scales = [
                (1e24, "Y" ),
                (1e21, "Z" ),
                (1e18, "E" ),
                (1e15, "P" ),
                (1e12, "T" ),
                (1e9,  "G"),
                (1e6,  "M"),
                (1e3,  "k") ]
            if abs(x) > 1
                for (mag, sym) in scales
                    if abs(x) >= mag
                        x /= mag
                        suffix = string(sym, suffix)
                        break
                    end
                end
            elseif T <: AbstractFloat
                smallscales = [
                    ( 1e-12, "p" ),
                    ( 1e-9,  "n" ),
                    ( 1e-6,  "μ" ),
                    ( 1e-3,  "m" ) ]
                for (mag,sym) in smallscales
                    if abs(x) < mag*10
                        x /= mag
                        suffix = string(sym, suffix)
                        break
                    end
                end
            end
        else
            if autoscale == :binary
                scales = [
                    (1024.0 ^8,  "Yi" ),
                    (1024.0 ^7,  "Zi" ),
                    (1024.0 ^6,  "Ei" ),
                    (1024.0 ^5,  "Pi" ),
                    (1024.0 ^4,  "Ti" ),
                    (1024.0 ^3,  "Gi"),
                    (1024.0 ^2,  "Mi"),
                    (1024.0,     "Ki")
                ]
            else # :finance
                scales = [
                    (1e12, "t" ),
                    (1e9,  "b"),
                    (1e6,  "m"),
                    (1e3,  "k") ]
            end
            for (mag, sym) in scales
                if abs(x) >= mag
                    x /= mag
                    suffix = string(sym, suffix)
                    break
                end
            end
        end
    end

    nonneg = x >= 0
    fractional = 0
    if T <: Rational && mixedfraction
        actualconv = "d"
        actualx = trunc( Int, x )
        fractional = abs(x) - abs(actualx)
    else
        actualx = (parens && !in( actualconv[1], "xX" )) ? abs(x) : x
    end
    s = cfmt( generate_format_string( width=width,
                                      precision=precision,
                                      leftjustified=leftjustified,
                                      zeropadding=zeropadding,
                                      commas=commas,
                                      signed=signed,
                                      positivespace=positivespace,
                                      alternative=alternative,
                                      conversion=actualconv
                                      ),
              actualx)

    if T <:Rational && conversion == "s"
        if mixedfraction && fractional != 0
            num = fractional.num
            den = fractional.den
            if tryden >= 2 && mod( tryden, den ) == 0
                num *= div(tryden,den)
                den = tryden
            end
            fs = string( num, fractionsep, den)
            length(fs) < fractionwidth &&
                (fs = string(repeat( "0", fractionwidth - length(fs) ), fs))
            s = (actualx != 0
                 ? string(rstrip(s), mixedfractionsep, fs)
                 : (nonneg ? fs : string('-', fs)))
            checkwidth = true
        elseif !mixedfraction
            s = replace( s, "//" => fractionsep )
            checkwidth = true
        end
    elseif stripzeros && in( actualconv[1], "fFeEs" )
        dpos = findfirst( isequal('.'), s )
        dpos === nothing && (dpos = length(s))
        if actualconv[1] in "eEs"
            epos = findfirst(isequal(actualconv[1] == 'E' ? 'E' : 'e'), s)
            rpos = (epos === nothing) ? length( s ) : (epos-1)
        else
            rpos = length(s)
        end
        # rpos at this point is the rightmost possible char to start
        # stripping
        stripfrom = rpos+1
        for i = rpos:-1:dpos+1
            if s[i] == '0'
                stripfrom = i
            elseif s[i] ==' '
                continue
            else
                break
            end
        end
        if stripfrom <= rpos
            # everything after decimal is 0, so strip the decimal too
            s = string(s[1:stripfrom-1-(stripfrom == dpos+1)], s[rpos+1:end])
            checkwidth = true
        end
    end

    s = string(s, suffix)

    if parens && !in( actualconv[1], "xX" )
        # if zero or positive, we still need 1 white space on the right
        s = nonneg ? string(' ', strip(s), ' ') : string('(', strip(s), ')')
        checkwidth = true
    end

    if checkwidth && width != -1
        if (len = length(s) - width) > 0
            s = replace( s, " " => ""; count=len )
            if (len = length(s) - width) > 0
                endswith(s, " ") && (s = reverse(replace(reverse(s), " " => ""; count=len)))
                (len = length(s) - width) > 0 && (s = replace( s, "," => ""; count=len ))
            end
        elseif len < 0
            # Todo: should use lpad or rpad here, can be more efficient
            s = leftjustified ? string(s, repeat( " ", -len )) : string(repeat( " ", -len), s)
        end
    end

    s
end


text:
#! format: off
#[PSSE 240 Bus Case system with Renewables](https://www.nrel.gov/grid/test-case-repository.html)

using PowerSimulationsDynamics
using PowerSystemCaseBuilder
using PowerSystems
using Sundials
using Plots
using OrdinaryDiffEq
PSD = PowerSimulationsDynamics

sys = build_system(PSSETestSystems, "psse_240_case_renewable_sys")

using Logging
sim_ida = PSD.Simulation(
    ResidualModel,
    sys, #system
    pwd(),
    (0.0, 20.0), #time span
    BranchTrip(1.0, Line, "CORONADO    -1101-PALOVRDE    -1401-i_10");
    console_level = Logging.Info,
)

PSD.execute!(sim_ida, IDA(), dtmax = 0.01)

res_ida = read_results(sim_ida)
v1101_ida = get_voltage_magnitude_series(res_ida, 1101);
plot(v1101_ida)

sim_rodas = PSD.Simulation(
    MassMatrixModel,
    sys, #system
    pwd(),
    (0.0, 20.0), #time span
    BranchTrip(1.0, Line, "CORONADO    -1101-PALOVRDE    -1401-i_10");
    console_level = Logging.Info,
)

PSD.execute!(
    sim_rodas,
    Rodas4(),
    saveat = 0.01,
    atol = 1e-10,
    rtol = 1e-10,
    initializealg = NoInit(),
)

res_rodas = read_results(sim_rodas)

v1101 = get_voltage_magnitude_series(res_rodas, 1101);
plot(v1101, label = "RODAS4")
plot!(v1101_ida, label = "IDA")

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
using Pkg

# Activate test environment on older Julia versions
@static if VERSION < v"1.2"
    Pkg.activate(@__DIR__)
    Pkg.develop(PackageSpec(; path=dirname(@__DIR__)))
    Pkg.instantiate()
end

using PythonOT
using Documenter
using Test

@testset "PythonOT.jl" begin
    @testset "doctests" begin
        DocMeta.setdocmeta!(PythonOT, :DocTestSetup, :(using PythonOT); recursive=true)
        doctest(
            PythonOT;
            doctestfilters=[
                r"{([a-zA-Z0-9]+,\s?)+[a-zA-Z0-9]+}",
                r"(Array{[a-zA-Z0-9]+,\s?1}|Vector{[a-zA-Z0-9]+})",
                r"(Array{[a-zA-Z0-9]+,\s?2}|Matrix{[a-zA-Z0-9]+})",
            ],
        )
    end
end


text:
#API for publishing and subscribing to message topics
export Publisher, Subscriber, publish

type Publisher{MsgType<:MsgT}
    o::PyObject

    function Publisher(topic::AbstractString; kwargs...)
        @debug("Creating <$(string(MsgType))> publisher on topic: '$topic'")
        rospycls = _get_rospy_class(MsgType)
        return new(__rospy__[:Publisher](ascii(topic), rospycls; kwargs...))
    end
end
Publisher{MsgType<:MsgT}(topic::AbstractString, ::Type{MsgType}; kwargs...) =
    Publisher{MsgType}(ascii(topic); kwargs...)

function publish{MsgType<:MsgT}(p::Publisher{MsgType}, msg::MsgType)
    pycall(p.o["publish"], PyAny, convert(PyObject, msg))
end

if VERSION >= v"0.5.0-dev+3692" #callbacks are broken

type Subscriber{T}
end
Subscriber(args...) = error(
"""Subscribing to a topic is currently broken on julia v0.5 and above. See
https://github.com/jdlangs/RobotOS.jl/issues/15 for ongoing efforts to fix this.""")

else #callbacks not broken

type Subscriber{MsgType<:MsgT}
    o::PyObject
    callback

    function Subscriber(topic::AbstractString, cb, cb_args::Tuple=(); kwargs...)
        @debug("Creating <$(string(MsgType))> subscriber on topic: '$topic'")
        rospycls = _get_rospy_class(MsgType)
        jl_cb(msg::PyObject) = cb(convert(MsgType, msg), cb_args...)
        return new(
            __rospy__[:Subscriber](ascii(topic), rospycls, jl_cb; kwargs...),
            jl_cb
        )
    end
end

Subscriber{MsgType<:MsgT}(
    topic::AbstractString,
    ::Type{MsgType},
    cb,
    cb_args::Tuple=();
    kwargs...
) = Subscriber{MsgType}(ascii(topic), cb, cb_args; kwargs...)

end #version check


text:
function render(f::Function, cam::C, sampler = SingleSampler(T);
                save_interval::Number=Inf, filename::String="render.png",
                resolution::Integer=1) where {T,C<:Camera{T}}
    w,h = width(cam),height(cam)
    img = zeros(RGB{eltype(cam)}, (w,h))

    W = weight(sampler)
    t₀ = time()
    @showprogress "Rendering: " for i = 1:resolution:w
        for j = 1:resolution:h
            img[i,j] = W*mapreduce(ij -> f(cam[ij...]), +, sampler(i,j))
            if resolution > 1
                n = resolution-1
                for i′ = i:i+n
                    for j′ = j:j+n
                        img[i′,j′] = img[i,j]
                    end
                end
            end
            now = time()
            if now - t₀ > save_interval
                save_clamped(filename, img)
                t₀ = now
            end
        end
    end

    img
end

function render(cam::C, acc::A,
                shade::Function, background::RGB{T},
                args...; kwargs...) where {T,C<:Camera{T},O,A<:Accelerator{T,O}}
    render(cam, args...; kwargs...) do ray::Ray
        hit = Intersection(ray, O)
        intersect!(hit, acc)
        is_hit(hit) ? shade(hit) : background
    end
end

export render


text:
export CuEquation, CuProblem, CuState, CuDualState

function CuProblem(g, v, p, eq, ts)
  st = CuState(cxeltype(eq.LC), size(eq.LC), ts.dt)
  Problem(g, v, p, eq, ts, st)
end

mutable struct CuState{T,dim} <: AbstractState
  t::Float64
  step::Int
  dt::Float64
  sol::CuArray{T,dim}
end

CuState(T::DataType, sz::Tuple, dt) = CuState(0.0, 0, dt, CuArray(zeros(T, sz)))

mutable struct CuDualState{T,dimc,dimr} <: AbstractState
  t::T
  step::Int
  dt::T
  solc::CuArray{T,dimc}
  solr::CuArray{T,dimr}
end

function CuDualState(T::DataType, sizec, sizer, dt)
  solc = CuArray(zeros(T, sizec))
  solr = CuArray(zeros(T, sizer))
  CuDualState(0, 0, dt, solc, solr)
end

# Equation Composite Type
"""
This type defines the linear implicit and explicit components of an equation.
The linear implicit part of an equation is defined by an array of coefficients
which multiply the solution. The explicit part of an equation is calculated
by a function that may define linear and nonlinear parts.
"""
struct CuEquation{T,dim} <: AbstractEquation
  LC::CuArray{T,dim} # Coeffs of the eqn's implicit linear part
  calcN!::Function # Function that calcs linear & nonlinear parts
end

struct DualCuEquation{T,dimc,dimr} <: AbstractEquation
  LCc::CuArray{Complex{T},dimc}
  LCr::CuArray{Complex{T},dimr}
  calcN!::Function
end

CuEquation(eq::AbstractEquation) = CuEquation(CuArray(eq.LC), eq.calcN!)
#CuEquation(eq::AbstractEquation) = CuEquation(eq.LC, eq.calcN!)


text:
module Transducers

export AdHocFoldable,
    Broadcasting,
    Cat,
    Completing,
    Consecutive,
    CopyInit,
    Count,
    Dedupe,
    DistributedEx,
    Drop,
    DropLast,
    DropWhile,
    Empty,
    Enumerate,
    Filter,
    FlagFirst,
    GroupBy,
    Init,
    Interpose,
    Iterated,
    KeepSomething,
    Map,
    MapCat,
    MapSplat,
    NondeterministicThreading,
    NotA,
    OfType,
    OnInit,
    Partition,
    PartitionBy,
    PreferParallel,
    ProductRF,
    ReduceIf,
    ReducePartitionBy,
    Reduced,
    Replace,
    Scan,
    ScanEmit,
    Scanx,
    SequentialEx,
    SplitBy,
    TCat,
    Take,
    TakeLast,
    TakeNth,
    TakeWhile,
    TeeRF,
    ThreadedEx,
    Transducer,
    Unique,
    Zip,
    append!!,
    channel_unordered,
    compose,
    dcollect,
    dcopy,
    dreduce,
    dtransduce,
    eduction,
    foldxd,
    foldxl,
    foldxt,
    ifunreduced,
    opcompose,
    push!!,
    reduced,
    reducingfunction,
    right,
    setinput,
    tcollect,
    tcopy,
    transduce,
    unreduced,
    whencombine,
    whencomplete,
    wheninit,
    whenstart,
    withprogress

using Base.Broadcast: Broadcasted
using Base: tail

import Adapt
import Setfield
import Tables
using ArgCheck
using BangBang.Experimental: modify!!, mergewith!!
using BangBang.NoBang: SingletonVector
using BangBang:
    @!, BangBang, Empty, append!!, collector, empty!!, finish!, push!!, setindex!!, union!!
using Baselet
using CompositionsBase: compose, opcompose
using DefineSingletons: @def_singleton
using Distributed: @everywhere, Distributed
using InitialValues:
    GenericInitialValue,
    InitialValue,
    InitialValues,
    SpecificInitialValue,
    asmonoid,
    hasinitialvalue
using Logging: @logmsg, LogLevel
using MicroCollections: UndefVector
using Requires
using Setfield: @lens, @set, set, setproperties
using SplittablesBase: SplittablesBase, amount, halve

# Dummy `ConstructionBase` module for supporting older `Setfield`:
module ConstructionBase
    using Setfield
    @static if isdefined(Setfield, :constructorof)
        const constructorof = Setfield.constructorof
    else
        const constructorof = Setfield.constructor_of
    end
end

@eval using CompositionsBase: $(Symbol("⨟"))
@eval export $(Symbol("⨟"))

@static if VERSION >= v"1.3-alpha"
    using Base.Threads: @spawn
    function nonsticky!(task)
        task.sticky = false
        return task
    end
else
    # Mock `@spawn` using `@async`:
    @eval const $(Symbol("@spawn")) = $(Symbol("@async"))
    nonsticky!(task) = task
end

# `AbstractArrayOrBroadcasted` is an internal detail of `Base`.  But
# we need this exact type for disambiguation...
if isdefined(Base, :AbstractArrayOrBroadcasted)
    const AbstractArrayOrBroadcasted = Base.AbstractArrayOrBroadcasted
else
    const AbstractArrayOrBroadcasted = Union{AbstractArray, Broadcasted}
end

include("AutoObjectsReStacker.jl")
using .AutoObjectsReStacker: restack

include("showutils.jl")
include("basics.jl")
include("core.jl")
include("library.jl")
include("teezip.jl")
include("groupby.jl")
include("broadcasting.jl")
include("consecutive.jl")
include("partitionby.jl")
include("splitby.jl")
include("combinators.jl")
include("simd.jl")
include("executors.jl")
include("processes.jl")
include("threading_utils.jl")
include("reduce.jl")
include("nondeterministic_threading.jl")
include("dreduce.jl")
include("unordered.jl")
include("air.jl")
include("lister.jl")
include("show.jl")
include("comprehensions.jl")
include("progress.jl")
include("deprecated.jl")

include("interop/blockarrays.jl")
include("interop/lazyarrays.jl")

function __init__()
    @require BlockArrays="8e7c35d0-a365-5155-bbbb-fb81a777f24e" begin
        __foldl__(rf, acc, coll::BlockArrays.BlockArray) =
            _foldl_blockarray(rf, acc, coll)
    end
    @require LazyArrays="5078a376-72f3-5289-bfd5-ec5146d43c02" begin
        __foldl__(rf, acc, coll::LazyArrays.Hcat) =
            _foldl_lazy_hcat(rf, acc, coll)
        __foldl__(rf, acc, coll::LazyArrays.Vcat) =
            _foldl_lazy_vcat(rf, acc, coll)
    end
    @require OnlineStatsBase="925886fa-5bf2-5e8e-b522-a9147a512338" begin
        include("interop/onlinestats.jl")
    end
    @require DataFrames="a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
        include("interop/dataframes.jl")
    end
    @require Referenceables="42d2dcc6-99eb-4e98-b66c-637b7d73030e" begin
        include("interop/referenceables.jl")
    end
end

end # module


text:
using Compat
using Compat.Test

using LinearAlgebra

@testset "Conjugate priors" begin
    tests = ["conjugates",
             "conjugates_normal",
             "conjugates_mvnormal"]

    for t in tests
        fpath = "$t.jl"
        println("running $fpath ...")
        include(fpath)
    end
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# # Maxima.jl のインストールの仕方 (Windows編)
#
# 黒木玄
#
# 2018-01-03, 2019-06-08, 2019-12-27, 2021-11-09
#
# Old Version: https://nbviewer.org/gist/genkuroki/6d32a79a04770dc63d366c2f1a0d61a4
#
# (0) 以下にアクセスして色々読む.
#
# https://github.com/nsmith5/Maxima.jl
#
# https://nsmith5.github.io/Maxima.jl/latest/
#
# (1) Maxima のバイナリをインストールする.
#
# https://sourceforge.net/projects/maxima/files/Maxima-Windows/
#
# 注意 2019-12-27：これを書いている時点の最新版 maxima-clisp-sbcl-5.41.0a-win64.exe をインストールすると gnuplot で作画できなくなる. それより少し古い maxima-clisp-sbcl-5.38.1.exe であれば gnuplot で作画可能である.
#
# 注意 2021-11-09: maxima-5.45.1s-win64.exe をダウンロードしてこのノートブックを正常に実行できている。
#
# (2) Maximaがインストールされているディレクトリ以下を探して maxima.bat を見付け, maxima.bat がおいてあるディレクトリを環境変数 PATH に追加する.
#
# (3) julia で
#
# ```julia
# run('maxima')
# ```
#
# を実行して, julia から maxima.bat が起動できることを確認する.
#
# (4) julia で以下を実行する.
#
# ```julia
# using Pkg
# Pkg.add("Maxima")
# ```
#
# (5) 以下を実行できるかを確認する:
#
# ```julia
# using Maxima
# a = m"sin(%pi/3)"
# b = integrate(m"log(sin(x))", 'x')
# c = integrate(m"x*log(sin(x))", 'x')
# ```
#
# Windows が「アクセスを許可してよいか？」のような意味の警告を発して来た場合には「許可」するように答える.
#
# 追記 2021-11-09: 次の一行を追加しないと、 `display(x)` を手動で実行しないとエラーになってしまう。
#
# ```julia
# Base.convert(::Type{Any}, x::MExpr) = x
# ```

# %%
using Maxima
Base.convert(::Type{Any}, x::MExpr) = x

# %%
a = m"sin(%pi/3)"

# %%
b = integrate(m"log(sin(x))", 'x')

# %%
c = integrate(m"x*log(sin(x))", 'x')

# %%
d = factor(MExpr("a^2 - b^2"))

# %%
A = diff(m"log(f(x))", 'x', 4)

# %%
B = integrate(m"log(sin(x))", 'x', m"0", m"%pi/2")

# %%
C = integrate(m"1/sqrt(1+x+x^2)", 'x')

# %%
# gnuplot が立ち上がる
plot2d(m"sin(x)", x=(-4,4), title="Sine Wave")

# %%
# gnuplot が立ち上がる
plot3d(
m"[5*cos(x)*(cos(x/2)*cos(y)+sin(x/2)*sin(2*y)+3.0)-10.0,
-5*sin(x)*(cos(x/2)*cos(y)+sin(x/2)*sin(2*y)+3.0),
5*(-sin(x/2)*cos(y)+cos(x/2)*sin(2*y))], [x, -%pi, %pi], [y, -%pi, %pi]")

# %%
using Base64
showimg(mime, fn; tag="img") = open(fn) do f
    base64 = base64encode(f)
    display("text/html", """<$tag src="data:$mime;base64,$base64" />""")
end

showimg("image/png", "2021-11-09_00.png")
showimg("image/png", "2021-11-09_01.png")

# %%

max_stars_repo_path: 0023/Maxima.jl example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0023/Maxima.jl example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0023/Maxima.jl example.jl

max_forks_repo_name: genkuroki/public


text:
# CAUTION there are several problems about these two packages:
#
# 1. The Cairo and Fontconfig packages must be loaded before ImageMagick.
#
# 2. And I must import Cairo first, then Fontconfig.
#
# 3. Also import Cairo is needed for Compose to draw an image.
import Cairo
import Fontconfig

# FIXME It looks like ImageMagick is not useful.
using ImageMagick

using ProgressMeter

using LightGraphs
# using LightGraphs.SimpleGraphs
using MetaGraphs
using Base.Iterators

using Statistics: I
# cov, etc
using Statistics

using GraphPlot: gplot, circular_layout
using Compose: PNG, draw

using Random

using SparseArrays: sparse

import Base.show
import Base.display
using Distributions

using HDF5
using Match
using RecursiveArrayTools
import CSV

using LinearAlgebra: Diagonal

using BSON: @save, @load

# export gen_graphs_hard, DataSpec, myplot, dataspec_to_id

include("display.jl")
include("data.jl")

function myplot(g, labels=nothing)
#     gplot(g, nodelabel=1:nv(g), layout=circular_layout)
    if isnothing(labels)
        if has_prop(g, :names)
            labels = get_prop(g, :names)
        else
            labels = 1:nv(g)
        end
    end
    gplot(g,
          layout=circular_layout,
          # TODO pass optional label
          nodelabel=labels,
          NODELABELSIZE = 4.0 * 1,
          # nodelabeldist=2,
          # nodelabelc="darkred",
          arrowlengthfrac = is_directed(g) ? 0.15 : 0.0)
end

"""Display a PNG by writing it to tmp file and show the filename. The filename
would be replaced by an Emacs plugin.

"""
function Base.display(d::EmacsDisplay, mime::MIME"image/png", x::AbstractGraph)
    # path, io = Base.Filesystem.mktemp()
    # path * ".png"
    path = tempname() * ".png"

    # Using GraphLayout plotting backend
    #
    # [random_layout, circular_layout, spring_layout, stressmajorize_layout, shell_layout, spectral_layout]
    # default: spring_layout, which is random
    draw(PNG(path), gplot(x,
                          layout=circular_layout,
                          # TODO pass optional label
                          nodelabel=1:nv(x),
                          NODELABELSIZE = 4.0 * 2,
                          # nodelabeldist=2,
                          # nodelabelc="darkred",
                          arrowlengthfrac = is_directed(x) ? 0.15 : 0.0))

    # loc_x, loc_y = layout_spring_adj(x)
    # draw_layout_adj(x, loc_x, loc_y, filename=path)

    println("$(path)")
    println("#<Image: $(path)>")
end

function display_graph_with_label(x::AbstractGraph, labels)
    path = tempname() * ".png"
    draw(PNG(path), gplot(x,
                          layout=circular_layout,
                          nodelabel=labels,
                          NODELABELSIZE = 4.0 * 1,
                          arrowlengthfrac = is_directed(x) ? 0.15 : 0.0))
    println("$(path)")
    println("#<Image: $(path)>")
end

##############################
## Generating DAG
##############################

function gen_ER(d, e)
    erdos_renyi(d, e, is_directed=true)
end

function gen_SF(d, e)
    # static_scale_free(n, m, 2)
    barabasi_albert(d, convert(Int, round(e/d)), is_directed=true)
end

function gen_ER_dag(d, e=d)
    gen_ER(d, e)  |> Graph |> random_orientation_dag |> MetaDiGraph
end

function gen_SF_dag(d, e=d)
    gen_SF(d, e)  |> Graph |> random_orientation_dag |> MetaDiGraph
end

function test()
    g = gen_ER(10, 10)
    g = gen_SF(10, 10)

    is_directed(g)
    is_cyclic(g)

    g2 = LightGraphs.random_orientation_dag(Graph(g))
    is_cyclic(g)
    is_cyclic(g2)

end

""" Get the weight adj matrix.
"""
function dag_W(dag::AbstractGraph{U}) where U
    # adjacency_matrix(g)
    n_v = nv(dag)
    nz = ne(dag)
    # colpt = ones(U, n_v + 1)

    W = zeros(n_v, n_v)

    # rowval = sizehint!(Vector{U}(), nz)
    for j in 1:n_v  # this is by column, not by row.
        dsts = sort(inneighbors(dag, j))
        for i in dsts
            W[i,j] = get_prop(dag, i, j, :weight)
        end
        # colpt[j + 1] = colpt[j] + length(dsts)
        # append!(rowval, dsts)
    end

    # FIXME why init to 1?
    # spmx = SparseMatrixCSC(n_v, n_v, colpt, rowval, ones(Int, nz))
    # return spmx
    return sparse(W)
end

""" Generate model parameters according to the DAG structure and causal mechanism

"""
function gen_weights(dag, zfunc=()->1)
    # dag is a meta graph, I'll need to attach parameters to each edge
    # 1. go through all edges
    # 2. attach edges
    # FIXME to ensure pureness, copy dag
    for e in edges(dag)
        set_prop!(dag, e, :weight, zfunc())
    end
    W = dag_W(dag)
    for e in edges(dag)
        rem_prop!(dag, e, :weight)
    end
    W
end

function gen_data2(W, noise, n)
    d = size(W, 1)
    X = zeros(n, d)
    g = DiGraph(W)

    if noise == :Gaussian
        noise_fn = randn
    elseif noise == :Poisson
        # FIXME fixed hyper-parameter \lambda
        d = Poisson(1)
        noise_fn = n->rand(d, n)
    elseif noise == :Exp
        d = Exponential(1)
        noise_fn = n->rand(d, n)
    elseif noise == :Gumbel
        d = Gumbel(0, 1)
        noise_fn = n->rand(d, n)
    else
        error("Noise model $noise not supported.")
    end

    # topological sort
    # for vertices in order
    for v in topological_sort_by_dfs(g)
        parents = inneighbors(g, v)
        # FIXME
        X[:, v] = X[:, parents] * W[parents, v] + noise_fn(n)
    end
    X
end

function my_uniform_init(din, dout)
    (rand(din, dout) .* 1.5 .+ 0.5) .* rand([1,-1], din, dout)
end

function test()
    # maybe this MLP weights are too small
    Dense(10, 10, initW=my_uniform_init).W
end

function gen_MLP(g, noise, n)
    adj = adjacency_matrix(g)
    d = size(adj, 1)
    # CAUTION MLP requires channel last
    X = zeros(d, n)

    dist = @match noise begin
        :Gaussian => Normal(0,1)
        # FIXME fixed hyper-parameter \lambda
        :Poisson => Poisson(1)
        :Exp => Exponential(1)
        :Gumbel => Gumbel(0, 1)
        _ => error("Noise model $noise not supported.")
    end
    noise_fn = n->rand(dist, n)

    # topological sort
    # for vertices in order
    # 1. get causal order
    for v in topological_sort_by_dfs(g)
        # 2. for causal order nodes
        #    2.1 get parents
        parents = inneighbors(g, v)
        if length(parents) == 0
            X[v, :] = noise_fn(n)
        else
            # 2.2 construct a random Dense Layer, with sigmoid activation. The
            # weights should be in some range, say [0.5,2], and TODO this might
            # be configurable by the spec.k as well
            mlp = Chain(Dense(length(parents), 10, σ,
                              # FIXME maybe this default is better
                              initW = my_uniform_init),
                        Dense(10, 1),
                        x->dropdims(x, dims=1))
            # 3.3 get result by applying the MLP
            X[v, :] = mlp(X[parents, :]) .+ noise_fn(n)
        end
    end
    # CAUTION and this is transposed to the final result
    X'
end

function test()
    # construct a G
    g = gen_ER_dag(13)
    # from g to binary adj matrix
    adj = adjacency_matrix(g)
    # genereate linear X
    W = gen_weights(g, ()->((rand() * 1.5 + 0.5) * rand([1,-1])))
    # generate linear X
    X = gen_data2(W, :Gaussian, 1000)
    # generate quadratic X
    size(gen_MLP(g, :Gaussian, 100))
end

function getch2(X)
    ch2 = Diagonal(dropdims(var(X, dims=1), dims=1))
    return cat(cor(X), ch2, dims=3)
end

function getch3(X)
    vars = dropdims(var(X, dims=1), dims=1)
    vars = (vars .- mean(vars)) ./ std(vars)
    ch3 = Diagonal(vars)
    return cat(cor(X), ch3, dims=3)
end

# TODO more variables
function gen_sup_data_internal(g, spec)
    d = nv(g)
    ds = map(1:spec.N) do i
        # DEBUG different weights
        # W = gen_weights(g)
        # W = gen_weights(g, ()->((rand() + 0.5) * rand([1,-1])))
        # W = gen_weights(g, ()->((rand() * 1.5 + 0.5) * rand([1,-1])))
        # W = gen_weights(g, ()->((rand() * (0.5+spec.k) + 0.5) * rand([1,-1])))
        # W = gen_weights(g, ()->((rand() * 1.5 + 0.5)))

        X = @match spec.mechanism begin
            :Linear => begin
                # [0.5, 0.5+k]
                # k=4, [+-0.5, +-4.5], mean(W) = 2.5
                #
                # X
                # cov(X)
                # var(X, dims=1)
                W = gen_weights(g, ()->((rand() * (0.5+spec.k) + 0.5) * rand([1,-1])))
                # FIXME the number of data points generated
                gen_data2(W, spec.noise, 1000)
            end
            :MLP =>
                # FIXME performance
                gen_MLP(g, spec.noise, 1000)
            _ => error("Unsupported mechanism.")
        end

        # cor(X), W
        # DEBUG one-hot encoding
        # cor(X), Flux.onehotbatch(W_bin, [0,1])

        # FIXME I should record both COR and COV
        # or, maybe I should just record W. g can be recovered from W, and
        @match spec.mat begin
            # DEBUG using adjacency_matrix(g) instead of W
            :COR => (cor(X), adjacency_matrix(g))
            :COV => (cov(X), adjacency_matrix(g))
            :medCOV => ((cov(X) ./ median(var(X, dims=1))), adjacency_matrix(g))
            :maxCOV => ((cov(X) ./ maximum(var(X, dims=1))), adjacency_matrix(g))
            :CH2 => (getch2(X), adjacency_matrix(g))
            :CH3 => (getch3(X), adjacency_matrix(g))
            _ => error("Unsupported matrix type")
        end
    end
    input = map(ds) do x x[1] end
    output = map(ds) do x x[2] end
    cat(input..., dims=4), cat(output..., dims=3)
end

function mycat(aoa; dims)
    # array of array
    arr = similar(aoa[1], size(aoa[1])..., length(aoa))
    for i in 1:length(aoa)
        if dims==3
            arr[:,:,:,i] = aoa[i]
        elseif dims==4
            arr[:,:,:,:,i] = aoa[i]
        else
            error("dims should be 3 or 4")
        end
    end
    # combine the last two dims FIXME the order unchanged?
    reshape(arr, size(arr)[1:end-2]..., :)
end

function gen_sup_data_with_graphs(spec, gs)
    # train data
    ds = @showprogress 0.1 "Generating.." map(gs) do g
        x, y = gen_sup_data_internal(g, spec)
    end
    input = map(ds) do x x[1] end
    output = map(ds) do x x[2] end
    @show size(input)
    @show size(input[1])
    @show size(output)

    # FIXME this splat is too much and cause stack overflow. But looks like I
    # have only the way of cat() ..
    #
    # cat(input..., dims=3), cat(output..., dims=3)

    # instead I'm going to pre-allocate an array and fill
    mycat(input, dims=4), mycat(output, dims=3)
    # This repo might be an alternative:
    # https://github.com/JuliaDiffEq/RecursiveArrayTools.jl
end

function gen_raw_data_with_graphs(spec, gs)
    ds = @showprogress 0.1 "Generating raw XY .." map(gs) do g
        d = nv(g)
        W = gen_weights(g, ()->((rand() * (0.5+spec.k) + 0.5) * rand([1,-1])))
        X = gen_data2(W, spec.noise, 1000)
        W[W .> 0] .= 1
        W[W .< 0] .= 1
        X, W
    end
    input = map(ds) do x x[1] end
    output = map(ds) do x x[2] end
    # mycat(input), mycat(output)
    # input, output
    # FIXME mycat is not robust enough for this
    cat(input..., dims=3), cat(output..., dims=3)
end

struct DataSpec
    d
    # weight range: [0.5, 0.5+k]
    # k: 1, 2, 4, 8
    k
    # gtype: :ER :SF
    gtype
    # noise: :Gaussian :Poisson
    noise

    # :COV or :COR
    mat
    # :Linear, :GP, :Quad, :MLP
    mechanism

    # I'll be training on different d, k, gtype, mat, mechanism
    #
    # for noise, I'll be only testing on different noise

    ng
    N
    bsize
    seed
end

# FIXME previous 10000, 10
function DataSpec(;d, k, gtype, noise, mat=:COV, mechanism=:Linear, ng=nothing, N=nothing, bsize=nothing, seed=1234)
    # FIXME maybe check error here
    #
    # UPDATE set ng and N based on d
    if isnothing(ng) || isnothing(N)
        if d <= 20
            ng = 3000
            N = 3
        elseif d <= 40
            ng = 2000
            N = 2
        elseif d <= 80
            # FIXME this might be too small. But should be more than enough if we
            # just use it as testing data
            ng = 1000
            N = 1
        elseif d <= 150
            ng = 1000
            N = 1
            # DEBUG testing memory limit during training (pullback of EQ layer seems
            # to consume lots of memory)
        elseif d <= 200
            ng = 1000
            N = 1
        elseif d <= 300
            ng = 300
            N = 1
        elseif d <= 400
            ng = 300
            N = 1
        end
    end
    # set bsize; batch size is not relevant to the data saving. Only relevant for data loading
    if isnothing(bsize)
        if d <= 30
            bsize=64
        elseif d<=50
            bsize=32
        elseif d<=70
            bsize=16
        elseif d<=100
            bsize=8
        elseif d<=400
            # FIXME still too large?
            bsize=8
        else
            error("too large d")
        end
    end
    DataSpec(d, k, gtype, noise, mat, mechanism, ng, N, bsize, seed)
end

function dataspec_to_id(spec::DataSpec)
    a = join(["d=$(spec.d)",
              "k=$(spec.k)",
              "gtype=$(spec.gtype)",
              "noise=$(spec.noise)",
              "mat=$(spec.mat)",
            "mec=$(spec.mechanism)",
            "ng=$(spec.ng)",
            "N=$(spec.N)"
              ],
             "_")
end

function dataspec_to_id(specs::Array{DataSpec, N} where N)
    # "[" * join(dataspec_to_id.(specs), "+") * "]"

    # FIXME however, for ensemble model, I actually just want a single model
    "ensemble"
end

function test()
    # this will always connect to the last dimension
    convert(Array, VectorOfArray([randn(10,5,3), randn(10,5,3)]))
end

function gs2hdf5mat(graphs)
    # from array of graphs to matrix suitable for saving as hdf5
    mats = Matrix.(adjacency_matrix.(graphs))
    # this seems to be slow
    @show size(mats)
    @show size(mats[1])
    @info "concating .."
    # onemat = cat(mats..., dims=3)
    onemat = convert(Array, VectorOfArray(mats))
    @info "done"
    onemat
end

function hdf5mat2gs(onemat)
    # and transfer back
    MetaDiGraph.([onemat[:,:,i] for i in 1:size(onemat, 3)])
end

function test()
    graphs = gen_graphs_hard(DataSpec(d=10, k=1, gtype=:ER,
                                      noise=:Gaussian))[1:11]
    recs = hdf5mat2gs(gs2hdf5mat(graphs))
    for (g,r) in zip(graphs, recs)
        g == r || error("Error")
    end
end

function tmp_convert_bson_hdf5()
    # loop through all the g.bson, and generate a g.hdf5
    # walk through the file system
    for dir in readdir("data")
        fname = joinpath("data", dir, "g.bson")
        fnew = joinpath(dirname(fname), "g.hdf5")
        if isfile(fname) && !isfile(fnew)
            @info "converting $fname to $fnew"
            @load fname train_gs test_gs
            h5open(fnew, "w") do file
                write(file, "train_gs", gs2hdf5mat(train_gs))
                write(file, "test_gs", gs2hdf5mat(test_gs))
            end
        end
    end
end

function load_sup_ds(spec, batch_size=100; use_raw=false, suffix="")
    # create "data/" folder is not already there
    if !isdir("data") mkdir("data") end
    # 1. generate graph
    gdir = "data/$(spec.gtype)-$(spec.d)-ng=$(spec.ng)-$(spec.seed)$suffix"
    if !isdir(gdir) mkdir(gdir) end
    gfile = "$gdir/g.hdf5"
    if isfile(gfile)
        # DEBUG loading
        # FIXME I don't need to load the graphs
        train_gs = hdf5mat2gs(h5read(gfile, "train_gs"))
        test_gs = hdf5mat2gs(h5read(gfile, "test_gs"))
    else
        @info "Generating graphs for " spec
        # generate graphs first
        # CAUTION this seed! might overwrite existing seeding for training and experiments.
        # I probably want to generate graphs before-hand?
        Random.seed!(spec.seed)
        graphs = gen_graphs_hard(spec)

        # ratio
        index = convert(Int, round(length(graphs) * 4 / 5))
        train_gs = graphs[1:index]
        test_gs = graphs[index:end]
        # DEBUG writing
        h5open(gfile, "w") do file
            write(file, "train_gs", gs2hdf5mat(train_gs))
            write(file, "test_gs", gs2hdf5mat(test_gs))
        end
    end

    fname = "$gdir/$(dataspec_to_id(spec)).hdf5"
    if ispath(fname)
        train_x = h5read(fname, "train_x")
        train_y = h5read(fname, "train_y")
        test_x = h5read(fname, "test_x")
        test_y = h5read(fname, "test_y")
        raw_x = h5read(fname, "raw_x")
        raw_y = h5read(fname, "raw_y")
    else
        # generate according to train_gs and test_gs
        # TODO use different graphs for ds and test
        @info "generating training ds .."
        train_x, train_y = gen_sup_data_with_graphs(spec, train_gs)
        @info  "generating testing ds .."
        test_x, test_y = gen_sup_data_with_graphs(spec, test_gs)

        # generate raw data for use with other methods
        #
        # FIXME the raw_x should not be just for test graphs. When training
        # NCC/RCC, I need them as training data
        raw_x, raw_y = gen_raw_data_with_graphs(spec, test_gs)

        h5open(fname, "w") do file
            write(file, "train_x", train_x)
            write(file, "train_y", train_y)
            write(file, "test_x", test_x)
            write(file, "test_y", test_y)
            write(file, "raw_x", raw_x)
            write(file, "raw_y", raw_y)
        end
        @info "Saved to $fname"
    end
    if use_raw
#         @info "raw info: " size(raw_x), size(raw_y)
        return DataSetIterator(raw_x, raw_y, batch_size)
    else
        return (DataSetIterator(train_x, train_y, batch_size),
                DataSetIterator(test_x, test_y, batch_size))
    end
end

function gen_graphs(spec)
    # generate n
    f = @match spec.gtype begin
        :ER => gen_ER_dag
        :SF => gen_SF_dag
        :ER2 => (d)->gen_ER_dag(d, 2*d)
        :ER4 => (d)->gen_ER_dag(d, 4*d)
        :SF2 => (d)->gen_SF_dag(d, 2*d)
        :SF4 => (d)->gen_SF_dag(d, 4*d)
        :Bern => error("Not implemented")
        _ => error("Not supported graph type.")
    end
    all_graphs = map(1:spec.ng) do i
        f(spec.d)
    end
    # @show length(all_graphs)
    # filter
    unique_graphs = unique(all_graphs) do g
        # FIXME can this be compared directly?
        adjacency_matrix(g)
    end
    # @show length(unique_graphs)
    unique_graphs
end

function gen_graphs_hard(spec)
    # CAUTION this function will loop forever to generate enough unique graph
    all_gs = []
    @showprogress 0.1 "generating graph.." for i in 1:10
        # @info "graph iter $i .."
        gs = gen_graphs(spec)
        append!(all_gs, gs)
        all_gs = unique(all_gs) do g
            # FIXME can this be compared directly?
            adjacency_matrix(g)
        end

        if length(all_gs) >= spec.ng
            all_gs = all_gs[1:spec.ng]
            break
        end
    end
    spec.ng == length(all_gs) || @warn "Not enougth unique graphs, $(length(all_gs)) < $n"
    all_gs
end

function named_graph(names)
    g = MetaDiGraph(length(names))
    set_prop!(g, :names, names)
    g
end

function named_graph_add_edge!(g, from, to)
    names = get_prop(g, :names)
    from_idx = findall(x->String(x)==String(from), names)[1]
    to_idx = findall(x->String(x)==String(to), names)[1]
    add_edge!(g, from_idx, to_idx)
end

function display_named_graph(g)
    names = get_prop(g, :names)
    path = tempname() * ".png"
    draw(PNG(path), gplot(g,
                          layout=circular_layout,
                          nodelabel=names,
                          NODELABELSIZE = 4.0 * 1,
                          arrowlengthfrac = is_directed(g) ? 0.15 : 0.0))
    println("$(path)")
    println("#<Image: $(path)>")
end

function Sachs_ground_truth()
    # load 2005-Science-Sachs-Causal
    df = CSV.read("Sachs/1.cd3cd28.csv")
    greal = named_graph(names(df))
    display_named_graph(greal)
#     named_graph_add_edge!(greal, :PKC, :PKA)
    # ground truth
    named_graph_add_edge!(greal, :PKC, :PKA)           # green
    named_graph_add_edge!(greal, :PKC, :pjnk)
    named_graph_add_edge!(greal, :PKC, :P38)
    named_graph_add_edge!(greal, :PKC, :praf)
    named_graph_add_edge!(greal, :PKC, :pmek)
    named_graph_add_edge!(greal, :PKA, :pjnk)
    named_graph_add_edge!(greal, :PKA, :P38)
    named_graph_add_edge!(greal, :PKA, :praf)
    named_graph_add_edge!(greal, :PKA, :pmek)
    named_graph_add_edge!(greal, :PKA, :pakts473)
    # FIXME :erk?
    named_graph_add_edge!(greal, :PKA, Symbol("p44/42"))
    named_graph_add_edge!(greal, :praf, :pmek)
    named_graph_add_edge!(greal, :pmek, Symbol("p44/42"))
    named_graph_add_edge!(greal, Symbol("p44/42"), :pakts473)           # green
    named_graph_add_edge!(greal, :plcg, :PIP2)
    named_graph_add_edge!(greal, :plcg, :PIP3)
    named_graph_add_edge!(greal, :PIP3, :PIP2)
    return greal
end

function read_syntren(fname)
    df = CSV.read(fname)
    X = convert(Matrix, df[:, Not(:GENE)])
    newdf = DataFrame(X', df[:, :GENE])
    return newdf
end

# read ground truth
function syntren_ground_truth(fname, node_names)
    lines = readlines(fname)
    # 1. get all the nodes
#     words = map(lines) do line
#         split(line)
#     end
#     unique_words = setdiff(unique(cat(words..., dims=1)), ["ac", "du", "re"])
    # 2. for each line, add edge
    # use node_names so that it is consistent with the prediction
    greal = named_graph(node_names)
    for line in lines
        from, rel, to = split(line)
        # there seems to be many self edges
        if from == to continue end
        # FIXME ac, re
        if rel == "ac"
            named_graph_add_edge!(greal, from, to)
        elseif rel == "re"
            named_graph_add_edge!(greal, from, to)
        elseif rel == "du"
            # FIXME???
            named_graph_add_edge!(greal, from, to)
            named_graph_add_edge!(greal, to, from)

        else
            error("Unsupported relation.")
        end
    end
    return greal
end

# TODO bnlearn dataset and groundtruth

# donwload data in notebook
# Assume data downloaded in data/bnlearn/xxx

function bnlearn_ground_truth()
    greal = named_graph([:A, :B, :C, :D, :E, :F, :G])
    named_graph_add_edge!(greal, :B, :C)
    named_graph_add_edge!(greal, :A, :C)
    named_graph_add_edge!(greal, :D, :B)
    named_graph_add_edge!(greal, :D, :F)
    named_graph_add_edge!(greal, :A, :F)
    named_graph_add_edge!(greal, :G, :F)
    named_graph_add_edge!(greal, :E, :F)
    greal
end

max_stars_repo_name: lihebi/DAG-EQ

max_issues_repo_name: lihebi/DAG-EQ

max_forks_repo_name: lihebi/DAG-EQ


text:
A = rand(4,3)   # simulate a matrix with elements selected at random from (0,1)

B = rand(3, 4)

A'   # "lazy" transpose

A*B

using LinearAlgebra, RCall, StatsModels, Tables

Formaldehyde = rcopy(R"Formaldehyde")

R"""
library(ggplot2)
qplot(x=carb, y=optden, data=Formaldehyde, geom="point")
"""

X = hcat(ones(size(Formaldehyde, 1)), Formaldehyde.carb)

y = Formaldehyde.optden

β = X\y     # least squares estimate

r = y - X*β   #residual

X'r   # not exactly zero but very small entries

using GLM

m1 = fit(LinearModel, @formula(optden ~ 1 + carb), Formaldehyde)

f1 = @formula(optden ~ 1 + carb)
y, X = modelcols(apply_schema(f1, schema(Formaldehyde)), Formaldehyde)

X

X

xpx = X'X

ch = cholesky(xpx)

ch.U'ch.U

ch.U'ch.U ≈ xpx

v = ldiv!(ch.U', X'y)

βc = ldiv!(ch.U, copy(v))     # solution from the Cholesky factorization

βc ≈ β

ldiv!(ch, X'y)

Xy = hcat(X, y)              # augmented model matrix

cha = cholesky(Xy'Xy)        # augmented Cholesky factor

RXX = UpperTriangular(view(cha.U, 1:2, 1:2))

RXX ≈ ch.U

rXy = cha.U[1:2, end]     # creates a copy ("view" doesn't copy)

rXy ≈ v

βac = ldiv!(RXX, copy(rXy))     # least squares solution from the augmented Cholesky

βac ≈ β

abs2(cha.U[end,end]) ≈ sum(abs2, y - X*β)   # check on residual sum of squares

Formrows = Tables.rowtable(Formaldehyde)

chr = cholesky(zeros(3, 3) + I)  # initialize

fill!(chr.factors, 0);   # zero out the contents
chr

fill!(chr.factors, 0)
for r in Formrows
    lowrankupdate!(chr, [1.0, r.carb, r.optden])
end
chr

qrX = qr(X)

qrX.R'qrX.R

qrX.R'qrX.R ≈ xpx


text:
module FESpaceConstructors

using Gridap
using Gridap.Helpers
import Gridap: FESpace

function FESpace(;kwargs...)

  reffe = _get_kwarg(:reffe,kwargs)
  @assert isa(reffe,Symbol) "For the moment, reffe can only be a symbol"

  model = _get_kwarg(:model,kwargs)
  labels = _get_kwarg(:labels,kwargs,nothing)
  conformity = _get_kwarg(:conformity,kwargs,true)

  diritags = _get_kwarg(:diritags,kwargs,Int[])
  dirimasks = _get_kwarg(:dirimasks,kwargs,nothing)

  order = _get_kwarg(:order,kwargs)

  polytope = _get_polytope(model)

  constraint = _get_kwarg(:constraint,kwargs,nothing)

  dim = celldim(model)

  fespace = nothing

  if reffe in [:Lagrangian,:QLagrangian,:PLagrangian]

    T = _get_kwarg(:valuetype,kwargs,nothing)
    if T == nothing
      error("valuetype is a mandatory keyword argument in FESpace constructor for Lagrangian reference FEs")
    end

    if _is_reffe_lagrangian_compatible_with_polytope(reffe,polytope)

      if conformity in [false, :L2]
        if labels == nothing
          fespace =  DLagrangianFESpace(T,model,order,diritags,dirimasks)
        else
          fespace =  DLagrangianFESpace(T,model,labels,order,diritags,dirimasks)
        end

      elseif conformity in [true, :default, :H1]
        if labels == nothing
          fespace = CLagrangianFESpace(T,model,order,diritags,dirimasks)
        else
          fespace = CLagrangianFESpace(T,model,labels,order,diritags,dirimasks)
        end

      else

        s = "Conformity $conformity not implemented for $reffe reference FE on a $(_poly_type(polytope))"
        error(s)

      end

    elseif reffe == :PLagrangian

      if conformity in [false, :L2]

        _reffe = PDiscRefFE(T,polytope,order)
        fespace = DiscFESpace(_reffe,model)

      else

        s = "Conformity $conformity not possible for $reffe reference FE on a $(_poly_type(polytope))"
        error(s)

      end

    else
      error("Reference element $reffe nor implemented on a $(_poly_type(polytope))")
    end

  elseif reffe == :RaviartThomas

    if ! _is_cube(polytope)
      error("RaviartThomas reference FEs can only be constructed on top of a ncube. Check your model.")
    end

    _reffe = RTRefFE(polytope,Float64,order)

    if conformity in [false, :L2]

      fespace = DiscFESpace(_reffe,model)

    elseif conformity in [true, :default, :HDiv]

      if labels == nothing
        _labels = FaceLabels(model)
      else
        _labels = labels
      end

      grid = Grid(model)
      trian = Triangulation(grid)
      graph = GridGraph(model)

      # fespace = ConformingFESpace(_reffe,trian,graph,_labels,diritags)
      fespace = DivConformingFESpace(_reffe,trian,graph,_labels,diritags)

    else

      s = "Conformity $conformity not possible for a $reffe reference FE"
      error(s)

    end

  elseif reffe == :Nedelec

    if ! _is_cube(polytope)
      error("Nedelec reference FEs can only be constructed on top of a ncube. Check your model.")
    end

    _reffe = NedelecRefFE(polytope,Float64,order)

    if conformity in [false, :L2]

      fespace = DiscFESpace(_reffe,model)

    elseif conformity in [true, :default, :HCurl]

      if labels == nothing
        _labels = FaceLabels(model)
      else
        _labels = labels
      end

      grid = Grid(model)
      trian = Triangulation(grid)
      graph = GridGraph(model)

      fespace = ConformingFESpace(_reffe,trian,graph,_labels,diritags)

    else

      s = "Conformity $conformity not possible for a $reffe reference FE"
      error(s)

    end

  else
    error("Reference element $reffe not implemented")
  end

  @assert fespace != nothing

  if constraint == nothing
    return fespace

  elseif constraint == :zeromean
    return ZeroMeanFESpace(fespace,order)

  else
    error("Unknown constraint value $constraint")

  end

end

function _get_kwarg(kwarg,kwargs)

  try
    return kwargs[kwarg]
  catch
    s = "The key-word argument $(kwarg) is mandatory in the FESpace constructor"
    error(s)
  end

end

function _get_kwarg(kwarg,kwargs,value)

  try
    return kwargs[kwarg]
  catch
    return value
  end

end

function _get_polytope(model)

  polys = CellPolytopes(Grid(model))

  @notimplementedif ! isa(polys,ConstantCellValue)

  polys.value

end

function _is_reffe_lagrangian_compatible_with_polytope(reffe,polytope)
  v = (reffe == :Lagrangian)
  v = v || (reffe == :QLagrangian && _is_cube(polytope))
  v = v || (reffe == :PLagrangian && _is_simplex(polytope))
  v
end

_is_cube(p) = all(p.extrusion.array .== HEX_AXIS)

_is_simplex(p) = all(p.extrusion.array .== TET_AXIS)

function _poly_type(p)

  if _is_cube(p)
    t = :ncube
  elseif _is_simplex(p)
    t = :simplex
  else
    @notimplemented
  end

  return t

end

end # module


text:
#### Select column subset from the EBD ####
# Because ebd_warblers_cut.csv cannot be loaded in memory anymore

include("../required.jl")

# Set path to files
cut_file = "./data/raw/ebd_warblers_cut.csv"
head_file = "./data/proc/ebd_warblers_cut_head.csv"

# Load first rows only
@time head1 = CSV.read(cut_file, DataFrame; limit=9)

# Select subset of names
show(stdout, "text/plain", names(head1))
names_subset = [
    # "GLOBAL UNIQUE IDENTIFIER",
    # "CATEGORY",
    "COMMON NAME",
    "SCIENTIFIC NAME",
    # "SUBSPECIES SCIENTIFIC NAME",
    # "OBSERVATION COUNT", # needed in tests
    # "COUNTRY CODE",
    "LATITUDE",
    "LONGITUDE",
    "OBSERVATION DATE",
    # "OBSERVER ID", # needed in tests
    # "SAMPLING EVENT IDENTIFIER", # needed in tests
    # "PROTOCOL TYPE",
    # "DURATION MINUTES",
    # "EFFORT DISTANCE KM",
    # "NUMBER OBSERVERS",
    # "ALL SPECIES REPORTED",
    "GROUP IDENTIFIER",
    "APPROVED",
]
findall(in(names_subset), names(head1))

# Load names subset on head lines first (to test)
head2 = CSV.read(cut_file, DataFrame; limit=9, select=names_subset)

# Test prepare_ebd_data! function on selected columns
prepare_ebd_data!(head2)
select!(head2, [:species, :commonName, :year, :latitude, :longitude, :groupIdentifier])

# Load names subset on full data set (to see if it fit in memory)
@time df = CSV.read(cut_file, DataFrame; select=names_subset)
# Fits in memory 🎉

# Test prepare_ebd_data! function on selected columns
prepare_ebd_data!(df)
select!(df, [:species, :commonName, :year, :latitude, :longitude, :groupIdentifier])
# Works 🎉

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
# ------------------------------ matrix updates -------------------------

function _βΠψdθρ!(tmp, ddm, θ)
    P = Πψtmp(tmp)
    y = psispace(ddm)
    β = discount(ddm)
    n = checksquare(P)
    n == length(y) || throw(DimensionMismatch("size _Πψtmp = $n versus size(y) = $(size(y))"))

    rwrd = reward(ddm)
    σ = theta_ρ(rwrd,θ)
    ρ = _ρ(revenue(rwrd), σ)
    dρdθρ = _dρdθρ(σ)

    Δ = 0.5 * step(y)

    @inbounds for (j,yj) in enumerate(y)
        if j == 1
            @. P[:,j] = β * _dcond_probdρ(yj,y,Δ,ρ) * dρdθρ
        elseif j == n
            @. P[:,j] = -β * _dcond_probdρ(yj,y,-Δ,ρ) * dρdθρ
        else
            @. P[:,j] = β * ( _dcond_probdρ(yj,y,Δ,ρ) - _dcond_probdρ(yj,y,-Δ,ρ) ) * dρdθρ
        end
    end
    return P
end

function _βΠψ!(tmp, ddm, θ)
    P = Πψtmp(tmp)
    y = psispace(ddm)
    β = discount(ddm)
    n = checksquare(P)
    n == length(y) || throw(DimensionMismatch("size _Πψtmp = $n versus size(y) = $(size(y))"))

    rwrd = reward(ddm)
    σ = theta_ρ(rwrd,θ)
    ρ = _ρ(revenue(rwrd), σ)

    Δ = 0.5 * step(y)

    @inbounds for (j,yj) in enumerate(y)
        if j == 1
            @. P[ :,j] = β * normcdf(cond_z(yj, y, Δ, ρ))
        elseif j == n
            @. P[ :,j] = β * normccdf(cond_z(yj, y, -Δ, ρ))
        else
            @. P[ :,j] = β * ( normcdf(cond_z(yj, y, Δ, ρ)) - normcdf(cond_z(yj, y, -Δ, ρ) ))
        end
    end
    return P
end

# ------------------------------ derivative check -------------------------

function check_dΠψ(ddm, θ)

    ψspace = psispace(ddm)

    rwrd = reward(ddm)
    θρ = theta_ρ(rwrd,θ)
    ρ = _ρ(revenue(rwrd), θρ)

    Δ = 0.5 * step(ψspace)

    for yj in ψspace
        for y in ψspace
            fdσ = Calculus.derivative((sig) -> normcdf(cond_z(yj, y, Δ, _ρ(sig))), θρ)
            dσ  = _dρdθρ(θρ) * _dcond_probdρ(yj, y, Δ, ρ)
            abs(fdσ - dσ ) < 1e-7 || throw(error("bad σ grad at σ = $σ, ψ2 = $yj, ψ1 = $y"))
        end
    end
    return true
end


text:
# Superoperator Tests
using Schrodinger
using Test, LinearAlgebra
using Schrodinger: apply_process, operator_to_choi, trace_norm,
    gate_fidelity_choi, kraus_to_natural, natural_to_kraus, natural_to_choi,
    choi_to_natural, kraus_to_choi, choi_to_kraus, choi_to_chi, chi_to_choi

println("Testing Superoperator Functions...")

@testset "Representation Changes" begin
d = 4
for i = 1:20
    U = rand_unitary(d)
    U2 = qeye(d) ⊗ U
    @test operator_to_choi(U) ≈ sum(U2*(ket((i,i),(d,d))*ket((j,j),(d,d))')*U2' for i=0:d-1,j=0:d-1)
end
end


text:
struct SingleVariant
    a::Union{Int32}
end

struct DoubleVariant
    a::Union{Int16, Int32}
end

struct SizeAlignMismatch
    a::Union{Tuple{Int16, Int16, Int16}, Int32}
end

struct UnionInTuple
    a::Tuple{Union{Int16, Int32}}
end

struct GenericInUnion{T}
    a::Union{T, Int32}
end

@testset "Structs with bits unions" begin
    @test begin
        b = JlrsReflect.reflect([SingleVariant])
        sb = JlrsReflect.StringWrappers(b)

        sb[JlrsReflect.basetype(SingleVariant)] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.SingleVariant")]
        pub struct SingleVariant {
            pub a: i32,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([DoubleVariant])
        sb = JlrsReflect.StringWrappers(b)

        sb[JlrsReflect.basetype(DoubleVariant)] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.DoubleVariant")]
        pub struct DoubleVariant {
            #[jlrs(bits_union_align)]
            _a_align: ::jlrs::wrappers::inline::union::Align4,
            #[jlrs(bits_union)]
            pub a: ::jlrs::wrappers::inline::union::BitsUnion<4>,
            #[jlrs(bits_union_flag)]
            pub a_flag: u8,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([SizeAlignMismatch])
        sb = JlrsReflect.StringWrappers(b)

        sb[JlrsReflect.basetype(SizeAlignMismatch)] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.SizeAlignMismatch")]
        pub struct SizeAlignMismatch {
            #[jlrs(bits_union_align)]
            _a_align: ::jlrs::wrappers::inline::union::Align4,
            #[jlrs(bits_union)]
            pub a: ::jlrs::wrappers::inline::union::BitsUnion<6>,
            #[jlrs(bits_union_flag)]
            pub a_flag: u8,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([UnionInTuple])
        sb = JlrsReflect.StringWrappers(b)

        sb[JlrsReflect.basetype(UnionInTuple)] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.UnionInTuple")]
        pub struct UnionInTuple<'frame, 'data> {
            pub a: ::jlrs::wrappers::ptr::ValueRef<'frame, 'data>,
        }"""
    end

    @test_throws ErrorException begin
        b = JlrsReflect.reflect([GenericInUnion])
    end
end


text:
# ---
# title: 1606. Find Servers That Handled Most Number of Requests
# id: problem1606
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Ordered Map
# link: <https://leetcode.com/problems/find-servers-that-handled-most-number-of-requests/description/>
# hidden: true
# ---
#
# You have `k` servers numbered from `0` to `k-1` that are being used to handle
# multiple requests simultaneously. Each server has infinite computational
# capacity but **cannot handle more than one request at a time**. The requests
# are assigned to servers according to a specific algorithm:
#
#   * The `ith` (0-indexed) request arrives.
#   * If all servers are busy, the request is dropped (not handled at all).
#   * If the `(i % k)th` server is available, assign the request to that server.
#   * Otherwise, assign the request to the next available server (wrapping around the list of servers and starting from 0 if necessary). For example, if the `ith` server is busy, try to assign the request to the `(i+1)th` server, then the `(i+2)th` server, and so on.
#
# You are given a **strictly increasing** array `arrival` of positive integers,
# where `arrival[i]` represents the arrival time of the `ith` request, and
# another array `load`, where `load[i]` represents the load of the `ith` request
# (the time it takes to complete). Your goal is to find the **busiest
# server(s)**. A server is considered **busiest** if it handled the most number
# of requests successfully among all the servers.
#
# Return _a list containing the IDs (0-indexed) of the **busiest server(s)**_.
# You may return the IDs in any order.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/08/load-1.png)
#
#
#
#     Input: k = 3, arrival = [1,2,3,4,5], load = [5,2,3,3,3]
#     Output: [1]
#     Explanation:
#     All of the servers start out available.
#     The first 3 requests are handled by the first 3 servers in order.
#     Request 3 comes in. Server 0 is busy, so it's assigned to the next available server, which is 1.
#     Request 4 comes in. It cannot be handled since all servers are busy, so it is dropped.
#     Servers 0 and 2 handled one request each, while server 1 handled two requests. Hence server 1 is the busiest server.
#
#
# **Example 2:**
#
#
#
#     Input: k = 3, arrival = [1,2,3,4], load = [1,2,1,2]
#     Output: [0]
#     Explanation:
#     The first 3 requests are handled by first 3 servers.
#     Request 3 comes in. It is handled by server 0 since the server is available.
#     Server 0 handled two requests, while servers 1 and 2 handled one request each. Hence server 0 is the busiest server.
#
#
# **Example 3:**
#
#
#
#     Input: k = 3, arrival = [1,2,3], load = [10,12,11]
#     Output: [0,1,2]
#     Explanation: Each server handles a single request, so they are all considered the busiest.
#
#
# **Example 4:**
#
#
#
#     Input: k = 3, arrival = [1,2,3,4,8,9,10], load = [5,2,10,3,1,2,2]
#     Output: [1]
#
#
# **Example 5:**
#
#
#
#     Input: k = 1, arrival = [1], load = [1]
#     Output: [0]
#
#
#
#
# **Constraints:**
#
#   * `1 <= k <= 105`
#   * `1 <= arrival.length, load.length <= 105`
#   * `arrival.length == load.length`
#   * `1 <= arrival[i], load[i] <= 109`
#   * `arrival` is **strictly increasing**.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module AttachmapBenchs

using Test
using InplaceArrays.Arrays
using InplaceArrays.Fields
using InplaceArrays.Fields: OtherMockBasis, MockBasis
using FillArrays
using InplaceArrays.TensorValues

@inline function loop(a,cache)
  for i in eachindex(a)
    ai = getindex!(cache,a,i)
  end
end

function bench1(n)
  p1 = Point(2,2)
  p2 = Point(4,2)
  p3 = Point(1,3)
  p4 = Point(5,2)
  x = [p1,p2,p3,p4]
  np = length(x)
  v = 2.0
  d = 2
  ndof = 8
  wi = 3.0
  w = fill(wi,ndof)
  f = OtherMockBasis{d}(ndof)
  r = MockBasis{d}(v,ndof)
  l = n
  ax = fill(x,l)
  af = Fill(f,l)
  ar = Fill(r,l)
  aw = fill(w,l)
  aϕ = lincomb(af,aw)
  ab = attachmap(ar,aϕ)
  abx = evaluate(ab,ax)
  cabx = array_cache(abx)
  @time loop(abx,cabx)
  a∇bx = evaluate(∇(ab),ax)
  ca∇bx = array_cache(a∇bx)
  @time loop(a∇bx,ca∇bx)
end

for n in (1,1,10,1000,100000)
  @eval begin
    println("+++ runing suite for n = $($n) +++")
    bench1($n)
  end
end

end # module


text:
using WCS
using Test

wcs = WCSTransform(2;
                   alt     = 'B',
                   cdelt   = [-0.066667, 0.066667],
                   colnum  = 2,
                   ctype   = ["RA---AIR", "DEC--AIR"],
                   crpix   = [-234.75, 8.3393],
                   crval   = [0.0, -90.0],
                   obsgeo  = [1.0, 2.0, 3.0],
                   pc      = [1.0 0.0; 0.0 1.0],
                   pv      = [(2, 1, 45.0)],
                   velangl = 3,
                   wcsname = "WCSLIB.jl")

@testset "transforms" begin
    pixcoords = [0.0 24.0 45.0;  # x coords
                 0.0 38.0 98.0]  # y coords

    worldcoords = pix_to_world(wcs, pixcoords)

    expected_world = [267.96547027 276.53931377 287.77080792;
                      -73.73660749 -71.97412809 -69.67813884]

    @test maximum(abs.(worldcoords .- expected_world)) < 5e-9
    pixcoords_out = world_to_pix(wcs, worldcoords)
    @test maximum(abs.(pixcoords_out .- pixcoords)) < 1e-9

    # Test Array{Float64, 1} methods of above
    worldcoords = pix_to_world(wcs, pixcoords[:, 1])
    @test maximum(abs.(worldcoords .- expected_world[:, 1])) < 5e-9
    pixcoords_out = world_to_pix(wcs, worldcoords)
    @test maximum(abs.(pixcoords_out .- pixcoords[:, 1])) < 1e-9
end

@testset "getproperty" begin
    @test wcs.ctype == ["RA---AIR", "DEC--AIR"]
    @test wcs.dateavg == ""
    @test wcs.alt == 'B'
    @test wcs.cdelt == [-0.066667, 0.066667]
    @test wcs.cd == [0.0 0.0; 0.0 0.0]
    @test wcs.obsgeo == (1.0, 2.0, 3.0,
                         63.43494882292201, -0.004025899381340196, -6.37813776403742e6)
    # Test that all public properties can be obtained
    for k in propertynames(wcs)
        @test getproperty(wcs, k) == getproperty(wcs, k)
    end
    # Test that accessing all private properties results in an error
    for k in setdiff(propertynames(wcs, true), propertynames(wcs))
        @test_throws ErrorException getproperty(wcs, k)
    end
end

@testset "header methods" begin
    header = "SIMPLE  =                    T / file does conform to FITS standard             BITPIX  =                  -64 / number of bits per data pixel                  NAXIS   =                    2 / number of data axes                            NAXIS1  =                 3636 / length of data axis 1                          NAXIS2  =                 1939 / length of data axis 2                          EXTEND  =                    T / FITS dataset may contain extensions            COMMENT   FITS (Flexible Image Transport System) format is defined in 'AstronomyCOMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H WCSAXES =                    2 / Number of coordinate axes                      CRPIX1  =               1818.0 / Pixel coordinate of reference point            CRPIX2  =                969.5 / Pixel coordinate of reference point            CDELT1  =             -0.00825 / [deg] Coordinate increment at reference point  CDELT2  =              0.00825 / [deg] Coordinate increment at reference point  CUNIT1  = 'deg'                / Units of coordinate increment and value        CUNIT2  = 'deg'                / Units of coordinate increment and value        CTYPE1  = 'RA---CEA'           / Right ascension, cylindrical equal area projectCTYPE2  = 'DEC--CEA'           / Declination, cylindrical equal area projection CRVAL1  =            3.575E+02 / [deg] Coordinate value at reference point      CRVAL2  =            0.000E+00 / [deg] Coordinate value at reference point      PV2_1   =                  1.0 / CEA projection parameter                       LONPOLE =                  0.0 / [deg] Native longitude of celestial pole       LATPOLE =                 90.0 / [deg] Native latitude of celestial pole        RADESYS = 'ICRS'               / Equatorial coordinate system                                                                                COMMENT  WCS header keyrecords produced by WCSLIB 5.9                           "

    ws = WCS.from_header(header)
    @test length(ws) == 1
    @test eltype(ws) === WCSTransform
    w = ws[1]

    @test w.naxis == 2

    # Check that `to_header` succeeds and has the right number
    # of records. (Checking exact equality of input and output
    # strings is too demanding, as the exact formatting differs
    # between minor WCSLIB versions. For example,
    # `0.0E+00` vs `0.000E+00`.)
    header_out = WCS.to_header(w)
    @test length(header_out) == 17 * 80

    # test relax keyword
    faulty_header = replace(header, "RADESYS = 'ICRS'" => "RADECSYS= 'FK5' ")
    wcs_relaxed = WCS.from_header(faulty_header)[1]  # default is maximum relaxation
    @test wcs_relaxed.radesys == "FK5"
    wcs_strict = WCS.from_header(faulty_header; relax=WCS.HDR_NONE)[1]  # strict
    @test wcs_strict.radesys == "ICRS"  # defaults to ICRS since EQUINOX not defined
end

@testset "utilities" begin
    @test_throws ErrorException wcs.obsgeo = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
    for ctrl in 0:2; @test_throws ErrorException WCS.obsfix(ctrl, wcs); end
    @test_throws ErrorException wcs.obsgeo = [1.0, 2.0, 3.0, 4.0, 5.0, WCS.UNDEFINED]
    WCS.obsfix(1, wcs)
    @test wcs.obsgeo == (1.0, 2.0, 3.0,
                         63.43494882292201, -0.004025899381340196, -6.37813776403742e6)
    @test_throws ErrorException wcs.obsgeo = [1.0, 2.0, 3.0, 4.0, 5.0, WCS.UNDEFINED]
    for ctrl in (0, 2); @test_throws ErrorException WCS.obsfix(ctrl, wcs); end
    # `obsfix(0, wcs)` errored but also reset the last element
    @test wcs.obsgeo == (1.0, 2.0, 3.0, 4.0, 5.0, -6.37813776403742e6)

    @test WCS.wcslib_version() > v"5.0"  # just check that it works

    # Test propagating errors from wcslib.
    # TODO: replace this with a public WCS routine that throws an error.
    @test_throws ErrorException WCS.assert_ok(Cint(1))

    # test rejecting non-ascii strings
    # TODO: replace this with a public WCS routine that throws an error.
    @test_throws AssertionError WCS.convert_string(NTuple{72, UInt8}, "àèìoù")
end


text:
using LatexSVG
using Test

include("test_latexdocument.jl")
include("test_latexengines.jl")
include("test_svgprocessing.jl")


text:
module CartesianDomains

export RectangularDomain

using Printf

abstract type AbstractDomain end

include("rectangular_domain.jl")

end


text:
using Base.Test
S = 1000
# Rows of Tr should sum to one
Tr = rand(S, S, 3);
Tr = Tr ./ sum(Tr, 2);
@test_approx_eq sum(Tr, 2) ones(size(Tr)[2:3])
# Columns of Tc should sum to one
Tc = similar(Tr);
for k = 1:3
  Tc[:, :, k] = Tr[:, :, k]';
end
@test_approx_eq sum(Tc, 1) ones(size(Tc)[2:3])

R = rand(S, 3);
Q = zeros(R);
V = zeros(S);

function bellman1(v::Vector{Float64}, t::Array{Float64,3},
                  r::Matrix{Float64}, q::Matrix{Float64}, d::Float64)
  for a = 1:size(t)[3]
    q[:, a] = r[:, a] + d * (t[:, :, a] * v)
  end
  [maximum(q[r, :]) for r = 1:size(t)[1]]
end

function bellman2(v::Vector{Float64}, t::Array{Float64,3},
                  r::Matrix{Float64}, q::Matrix{Float64}, d::Float64)
  for a = 1:size(t)[3]
    for s = 1:size(t)[1]
      q[s, a] = (r[s, a] + d * (t[s, :, a] * v))[1]
    end
  end
  [maximum(q[r, :]) for r = 1:size(t)[1]]
end

function bellman3(v::Vector{Float64}, t::Array{Float64,3},
                  r::Matrix{Float64}, q::Matrix{Float64}, d::Float64)
  for a = 1:size(t)[3]
    for s = 1:size(t)[1]
      q[s, a] = (r[s, a] + d * (v' * t[:, s, a]))[1]
    end
  end
  [maximum(q[r, :]) for r = 1:size(t)[1]]
end

function bellman4(v::Vector{Float64}, t::Array{Float64,3},
                  r::Matrix{Float64}, q::Matrix{Float64}, d::Float64)
  s1, s2, a = size(t)
  @assert length(v) == s1
  @assert size(r) == (s2, a)
  @assert size(q) == (s2, a)
  @inbounds for k = 1:a
    for j = 1:s2
      S = 0.0
      @simd for i = 1:s1
        S += t[i, j, k] * v[i]
      end
      q[j, k] = r[j, k] + d * S
    end
  end
  [maximum(q[r, :]) for r = 1:size(t)[1]]
end

function bellman5(v::Vector{Float64}, t::Array{Float64,3},
                  r::Matrix{Float64}, q::Matrix{Float64}, d::Float64)
  for a = 1:size(t)[3]
    for s1 = 1:size(t)[1]
      S = 0.0
      for s2 = 1:size(t)[1]
        S += t[s2, s1, a] * v[s2]
      end
      q[s1, a] = r[s1, a] + d * S
    end
  end
  [maximum(q[r, :]) for r = 1:size(t)[1]]
end

function bellman6(v::Vector{Float64}, t::Array{Float64,3},
                  r::Matrix{Float64}, q::Matrix{Float64}, d::Float64)
  s1, s2, a = size(t)
  @assert length(v) == s1
  @assert size(r) == (s2, a)
  @assert size(q) == (s2, a)
  @inbounds for k = 1:a
    for i = 1:s2
      S = 0.0
      @simd for j = 1:s1
        S += t[i, j, k] * v[j]
      end
      q[i, k] = r[i, k] + d * S
    end
  end
  [maximum(q[r, :]) for r = 1:size(t)[1]]
end

function bellman7(newv::Vector, v::Vector{Float64}, t::Array{Float64,3},
                  r::Matrix{Float64}, d::Float64, policy::Vector{Int})
  s1, s2, a = size(t)
  @assert length(v) == length(newv) == s1
  @assert size(r) == (s2, a)
  copy!(newv, v)
  @inbounds for k = 1:a
    for j = 1:s2
      S = 0.0
      @simd for i = 1:s1
        S += t[i, j, k] * v[i]
      end
      q = r[j, k] + d * S
      if newv[j] < q
        newv[j] = q
        policy[j] = k
      end
    end
  end
end

function bellman7b(newv::Vector, v::Vector{Float64}, t::Array{Float64,3},
                  r::Matrix{Float64}, d::Float64, policy::Vector{Int})
  s1, s2, a = size(t)
  @assert length(v) == length(newv) == s1
  @assert size(r) == (s2, a)
  @inbounds begin
    for i = 1:s2
      newv[i] = v[i]
    end
    for k = 1:a
      for j = 1:s2
        S = 0.0
        @simd for i = 1:s1
          S += t[i, j, k] * v[i]
        end
        q = r[j, k] + d * S
        if newv[j] < q
          newv[j] = q
          policy[j] = k
        end
      end
    end
  end
end

function bellman8(v::Vector{Float64}, t::Array{Float64,3},
                  r::Matrix{Float64}, q::Matrix{Float64}, d::Float64)
  for a = 1:size(t)[3]
    tv = t[:, :, a] * v
    q[:, a] = r[:, a] + d * tv
  end
  [maximum(q[r, :]) for r = 1:size(t)[1]]
end

function bellman9(v::Vector{Float64}, t::Array{Float64,3},
                  r::Matrix{Float64}, q::Matrix{Float64}, d::Float64)
  s1, s2, a = size(t)
  @assert length(v) == s1
  @assert size(r) == (s1, a)
  @assert size(q) == (s1, a)
  @inbounds for a = 1:size(t)[3]
    tv = t[:, :, a] * v
    @simd for i = 1:s1
      q[i, a] = r[i, a] + d * tv[i]
    end
  end
  [maximum(q[r, :]) for r = 1:size(t)[1]]
end

V1 = bellman1(V, Tr, R, Q, 0.9);
P1 = [indmax(Q[s, :]) for s = 1:size(Tr)[1]];

V2 = bellman2(V, Tr, R, Q, 0.9);
P2 = [indmax(Q[s, :]) for s = 1:size(Tr)[1]];
@test_approx_eq V1  V2
@assert P1 == P2

V3 = bellman3(V, Tc, R, Q, 0.9);
P3 = [indmax(Q[s, :]) for s = 1:size(Tr)[1]];
@test_approx_eq V1  V3
@assert P1 == P3

V4 = bellman4(V, Tc, R, Q, 0.9);
P4 = [indmax(Q[s, :]) for s = 1:size(Tc)[1]];
@test_approx_eq V1  V4
@assert P1 == P4

V5 = bellman5(V, Tc, R, Q, 0.9);
P5 = [indmax(Q[s, :]) for s = 1:size(Tc)[1]];
@test_approx_eq V1  V5
@assert P1 == P5

V6 = bellman6(V, Tr, R, Q, 0.9);
P6 = [indmax(Q[s, :]) for s = 1:size(Tr)[1]];
@test_approx_eq V1  V6
@assert P1 == P6

V7new = similar(V)
V7old = copy(V)
P7 = ones(Int, S)
bellman7(V7new, V7old, Tc, R, 0.9, P7);
@test_approx_eq V1 V7new
@assert P1 == P7

V7bnew = similar(V)
V7bold = copy(V)
P7b = ones(Int, S)
bellman7b(V7bnew, V7bold, Tc, R, 0.9, P7b);
@test_approx_eq V1 V7bnew
@assert P1 == P7b

V8 = bellman8(V, Tr, R, Q, 0.9);
P8 = [indmax(Q[s, :]) for s = 1:size(Tr)[1]];
@test_approx_eq V1  V8
@assert P1 == P8

# timings

n = 100

print("Test 1: ")
@time for x = 1:n
  V1 = bellman1(V1, Tr, R, Q, 0.9)
end

print("Test 2: ")
@time for x = 1:n
  V2 = bellman2(V2, Tr, R, Q, 0.9)
end

print("Test 3: ")
@time for x = 1:n
  V3 = bellman3(V3, Tc, R, Q, 0.9)
end

print("Test 4: ")
@time for x = 1:n
  V4 = bellman4(V4, Tc, R, Q, 0.9)
end

print("Test 5: ")
@time for x = 1:n
  V5 = bellman5(V5, Tc, R, Q, 0.9)
end

print("Test 6: ")
@time for x = 1:n
  V6 = bellman6(V6, Tr, R, Q, 0.9)
end

print("Test 7: ")
@time for x = 1:n
  bellman7(V7new, V7old, Tc, R, 0.9, P7)
end

print("Test 7b: ")
@time for x = 1:n
  bellman7b(V7bnew, V7bold, Tc, R, 0.9, P7b)
end

print("Test 8: ")
@time for x = 1:n
  V8 = bellman6(V8, Tr, R, Q, 0.9)
end