text:
using LaplaceApproximations
using LaplaceApproximations.GPML
using Distributions
using LinearAlgebra
using StatsFuns
# using Plots
# using Optim

"Matérn (ν = 3/2) correlation kernel"
function mat32(x1, x2, ρ = 1)
    d = norm(x1 .- x2)
    _mat32(d, ρ)
end
function _mat32(d, ρ = 1)
    (1 + sqrt(3) * d / ρ) * exp(-sqrt(3) * d / ρ)
end
function ddρ_mat32(d, σ, ρ)
    3 * d^2 * σ^2 * exp(-sqrt(3) * d / ρ) / ρ^3
end

function mat32cov_funs(x)
    d = norm.(x .- x')
    function f(p)
        p[1] * _mat32.(d, p[2])
    end
    function gr(p)
        ddσ = 2 * p[1] * _mat32.(d, p[2])
        ddρ = ddρ_mat32.(d, p[1], p[2])
        (ddσ, ddρ)
    end
    (f = f, gr = gr)
end

n = 100
xobs = 5rand(n)
sort!(xobs)
xpred = range(0, 5, length = 101)

mat32cov = mat32cov_funs(xobs)
Σ = Matrix{Float64}(undef, n, n)
for j in 1:n, i in 1:n
    Σ[i, j] = mat32(xobs[i], xobs[j])
    Σ[j, i] = Σ[i, j]
end

utrue = rand(MvNormal(Σ))

# Normal example
σobs = 0.05
normy = utrue .+ σobs * randn(n)
normdatalik = GPML.normal_funs(normy, σobs)
normlogq, normyhat = GPML.gpml_laplace(normdatalik, Σ, zeros(n); maxit = 2)

# normpost = GPML.post_funs(normdatalik, Σ)
# normlogq_la, normyhat_la = laplace_approx(normpost.f, zeros(n); maxit = 10)

# normlogq_la - normlogq
# isapprox(normuhat_la, normyhat, normdatalik, zeros(n))

# plot(xpred, getindex.(normpreds, 2))
# scatter!(xobs, normy)
# scatter!(xobs, yobs)

GPML.gpml_grad([1.0, 1.0], mat32cov, normdatalik, zeros(n))

# normval(covpar) = GPML.gpml_grad(covpar, mat32cov, normdatalik, zeros(n))[1]
# normgr(G, covpar) = G .= GPML.gpml_grad(covpar, mat32cov, normdatalik, zeros(n))[2]
# Not currently working; end up with non-PD B matrix at some point
# maximize(normval, normgr, ones(2), ConjugateGradient())
# maximize(normval, ones(2))

# Logistic example
pobs = logistic.(utrue)
yobs = rand.(Bernoulli.(pobs))
ldatalik = GPML.logistic_funs(yobs)
llogq, luhat = GPML.gpml_laplace(ldatalik, Σ, zeros(n))
lpreds = map(xpred) do x
    GPML.gpml_pred(x, ldatalik, luhat, xobs, Σ, mat32, logistic)
end

# plot(xpred, getindex.(lpreds, 2))
# plot!(xobs, utrue)
# scatter!(xobs, yobs)
# plot!(xpred, getindex.(lpreds, 1))

# lres = maximize(p -> GPML.gpml_laplace(
#     ldatalik,
#     p[1]^2 * mat32.(xobs, xobs', p[2]),
#     ones(n))[1], ones(2))
# lcovpars = Optim.maximizer(lres)
# lq, luhat = GPML.gpml_laplace(ldatalik, lcovpars[1] * mat32.(xobs, xobs', covpars[2]), ones(n))

# Probit example
pyobs = @. rand(Bernoulli(normcdf(utrue)))
pdatalik = GPML.probit_funs(pyobs)
plogq, puhat = GPML.gpml_laplace(pdatalik, Σ, zeros(n))
ppreds = map(xpred) do x
    GPML.gpml_pred(x, pdatalik, puhat, xobs, Σ, mat32, normcdf)
end

# plot(xpred, getindex.(ppreds, 2))
# plot!(xobs, utrue)
# plot!(xpred, getindex.(ppreds, 1))
# scatter!(xpred, pyobs)

ppost = GPML.post_funs(pdatalik, Σ)
u = zeros(n)
for _ in 1:10
    global u -= ppost.he(u) \ ppost.gr(u)
end


text:
export all_partitions
"""
`all_partitions(A::Set)` creates a `Set` containing all possible
partitions of the set `A`.

`all_partitions(n::Int)` creates the `Set` of all partitions of
the set `{1,2,...,n}`.

Both of these take an optional second argument `k` to specify that
only partitions with exactly `k` parts should be returned.
"""
function all_partitions(A::Set{T}) where T
    if length(A) < 2
        P = Partition(A)
        return Set([P])
    end
    # Set aside one element of A and recurse
    x = first(A)
    B = deepcopy(A)
    delete!(B,x)

    PB = all_partitions(B)
    PA = Set{Partition{T}}()  # place to hold partitions we create

    for P in PB
        # case 1: include x as a singleton
        P_parts = parts(P)
        push!(P_parts, Set([x]))
        Q = PartitionBuilder(P_parts)
        push!(PA,Q)

        # case 2: insert x into existing parts
        parts_list = collect(parts(P))
        np = length(parts_list)
        for k=1:np
            push!(parts_list[k],x)  # insert x into k'th part
            Q = PartitionBuilder(Set(parts_list)) # build the partition
            push!(PA,Q)
            delete!(parts_list[k],x) # take it back out
        end

    end
    return PA
end

function all_partitions(n::Int)
    if n < 0
        error("argument must be a nonnegative integer")
    end
    A = Set{Int}(collect(1:n))
    return all_partitions(A)
end

function all_partitions(A::Set{T},k::Int) where T
    if k<0
        error("Number of partitions (k) must be nonnegative")
    end
    n = length(A)
    PA = Set{Partition{T}}()
    if k>n  # no partitions possible
        return PA
    end
    if n==0 # special case when A is empty
        push!(PA, Partition(A))
        return PA
    end
    if k==0 # special case, no partitions
        return PA
    end

    x = first(A)
    B = deepcopy(A)
    delete!(B,x)

    # Step 1: partitions in which x is a singleton
    PB = all_partitions(B,k-1)
    for P in PB
        parts_set = parts(P)
        push!(parts_set, Set([x]))  # {x} as a single part
        Q = PartitionBuilder(parts_set)
        push!(PA,Q)
    end

    # Step 2: partitions in which x is with others
    PB = all_partitions(B,k)
    for P in PB
        for j=1:k  # each P has exactly k parts
            parts_list = collect(parts(P))
            push!(parts_list[j], x)  # add x to part j
            Q = PartitionBuilder(Set(parts_list))
            push!(PA,Q)
        end
    end

    return PA
end

function all_partitions(n::Int,k::Int)
    if n < 0
        error("argument must be a nonnegative integer")
    end
    A = Set{Int}(collect(1:n))
    return all_partitions(A,k)
end

ext: jl


text:
module Inference

using ..Core, ..Utilities
using DynamicPPL: Metadata, _tail, VarInfo, TypedVarInfo,
    islinked, invlink!, getlogp, tonamedtuple, VarName, getsym, vectorize,
    settrans!, _getvns, getdist, split_var_str, CACHERESET, AbstractSampler,
    Model, runmodel!, Sampler, SampleFromPrior, SampleFromUniform,
    Selector, AbstractSamplerState, DefaultContext, PriorContext,
    LikelihoodContext, MiniBatchContext, set_flag!, unset_flag!, NamedDist, NoDist
using Distributions, Libtask, Bijectors
using DistributionsAD: VectorOfMultivariate
using LinearAlgebra
using ..Turing: PROGRESS, Turing
using StatsFuns: logsumexp
using Random: GLOBAL_RNG, AbstractRNG, randexp
using DynamicPPL
using AbstractMCMC: AbstractModel, AbstractSampler
using Bijectors: _debug
using MCMCChains: Chains

import AbstractMCMC
import AdvancedHMC; const AHMC = AdvancedHMC
import AdvancedMH; const AMH = AdvancedMH
import ..Core: getchunksize, getADbackend
import DynamicPPL: tilde, dot_tilde, getspace, get_matching_type,
    VarName, _getranges, _getindex, getval, _getvns
import EllipticalSliceSampling
import Random

export  InferenceAlgorithm,
        Hamiltonian,
        GibbsComponent,
        StaticHamiltonian,
        AdaptiveHamiltonian,
        SampleFromUniform,
        SampleFromPrior,
        MH,
        ESS,
        Gibbs,      # classic sampling
        HMC,
        SGLD,
        SGHMC,
        HMCDA,
        NUTS,       # Hamiltonian-like sampling
        DynamicNUTS,
        IS,
        SMC,
        CSMC,
        PG,
        assume,
        dot_assume,
        observe,
        dot_observe,
        resume

#######################
# Sampler abstraction #
#######################
abstract type AbstractAdapter end
abstract type InferenceAlgorithm end
abstract type ParticleInference <: InferenceAlgorithm end
abstract type Hamiltonian{AD} <: InferenceAlgorithm end
abstract type StaticHamiltonian{AD} <: Hamiltonian{AD} end
abstract type AdaptiveHamiltonian{AD} <: Hamiltonian{AD} end

getchunksize(::Type{<:Hamiltonian{AD}}) where AD = getchunksize(AD)
getADbackend(::Hamiltonian{AD}) where AD = AD()

"""
    mh_accept(logp_current::Real, logp_proposal::Real, log_proposal_ratio::Real)

Decide if a proposal ``x'`` with log probability ``\\log p(x') = logp_proposal`` and
log proposal ratio ``\\log k(x', x) - \\log k(x, x') = log_proposal_ratio`` in a
Metropolis-Hastings algorithm with Markov kernel ``k(x_t, x_{t+1})`` and current state
``x`` with log probability ``\\log p(x) = logp_current`` is accepted by evaluating the
Metropolis-Hastings acceptance criterion
```math
\\log U \\leq \\log p(x') - \\log p(x) + \\log k(x', x) - \\log k(x, x')
```
for a uniform random number ``U \\in [0, 1)``.
"""
function mh_accept(logp_current::Real, logp_proposal::Real, log_proposal_ratio::Real)
    # replacing log(rand()) with -randexp() yields test errors
    return log(rand()) + logp_current ≤ logp_proposal + log_proposal_ratio
end

######################
# Default Transition #
######################

struct Transition{T, F<:AbstractFloat}
    θ  :: T
    lp :: F
end

function Transition(spl::Sampler, nt::NamedTuple=NamedTuple())
    theta = merge(tonamedtuple(spl.state.vi), nt)
    lp = getlogp(spl.state.vi)
    return Transition{typeof(theta), typeof(lp)}(theta, lp)
end

function additional_parameters(::Type{<:Transition})
    return [:lp]
end

##########################################
# Internal variable names for MCMCChains #
##########################################

const TURING_INTERNAL_VARS = (internals = [
    "elapsed",
    "eval_num",
    "lf_eps",
    "lp",
    "weight",
    "le",
    "acceptance_rate",
    "hamiltonian_energy",
    "hamiltonian_energy_error",
    "max_hamiltonian_energy_error",
    "is_accept",
    "log_density",
    "n_steps",
    "numerical_error",
    "step_size",
    "nom_step_size",
    "tree_depth",
    "is_adapt",
],)

#########################################
# Default definitions for the interface #
#########################################

function AbstractMCMC.sample(
    model::AbstractModel,
    alg::InferenceAlgorithm,
    N::Integer;
    kwargs...
)
    return AbstractMCMC.sample(Random.GLOBAL_RNG, model, alg, N; kwargs...)
end

function AbstractMCMC.sample(
    rng::AbstractRNG,
    model::AbstractModel,
    alg::InferenceAlgorithm,
    N::Integer;
    chain_type=Chains,
    resume_from=nothing,
    progress=PROGRESS[],
    kwargs...
)
    if resume_from === nothing
        return AbstractMCMC.sample(rng, model, Sampler(alg, model), N;
                                   chain_type=chain_type, progress=progress, kwargs...)
    else
        return resume(resume_from, N; chain_type=chain_type, progress=progress, kwargs...)
    end
end

function AbstractMCMC.psample(
    model::AbstractModel,
    alg::InferenceAlgorithm,
    N::Integer,
    n_chains::Integer;
    kwargs...
)
    return AbstractMCMC.psample(Random.GLOBAL_RNG, model, alg, N, n_chains; kwargs...)
end

function AbstractMCMC.psample(
    rng::AbstractRNG,
    model::AbstractModel,
    alg::InferenceAlgorithm,
    N::Integer,
    n_chains::Integer;
    chain_type=Chains,
    progress=PROGRESS[],
    kwargs...
)
    return AbstractMCMC.psample(rng, model, Sampler(alg, model), N, n_chains;
                                chain_type=chain_type, progress=progress, kwargs...)
end

function AbstractMCMC.sample_init!(
    ::AbstractRNG,
    model::Model,
    spl::Sampler,
    N::Integer;
    kwargs...
)
    # Resume the sampler.
    set_resume!(spl; kwargs...)

    # Set the parameters to a starting value.
    initialize_parameters!(spl; kwargs...)
end

function AbstractMCMC.sample_end!(
    ::AbstractRNG,
    ::Model,
    ::Sampler,
    ::Integer,
    ::Vector;
    kwargs...
)
    # Silence the default API function.
end

function initialize_parameters!(
    spl::Sampler;
    init_theta::Union{Nothing,Vector}=nothing,
    verbose::Bool=false,
    kwargs...
)
    # Get `init_theta`
    if init_theta !== nothing
        verbose && @info "Using passed-in initial variable values" init_theta
        # Convert individual numbers to length 1 vector; `ismissing(v)` is needed as `size(missing)` is undefined`
        init_theta = [ismissing(v) || size(v) == () ? [v] : v for v in init_theta]
        # Flatten `init_theta`
        init_theta_flat = foldl(vcat, map(vec, init_theta))
        # Create a mask to inidicate which values are not missing
        theta_mask = map(x -> !ismissing(x), init_theta_flat)
        # Get all values
        theta = spl.state.vi[spl]
        @assert length(theta) == length(init_theta_flat) "Provided initial value doesn't match the dimension of the model"
        # Update those which are provided (i.e. not missing)
        theta[theta_mask] .= init_theta_flat[theta_mask]
        # Update in `vi`
        spl.state.vi[spl] = theta
    end
end

##########################
# Chain making utilities #
##########################

function _params_to_array(ts::Vector, spl::Sampler)
    names_set = Set{String}()
    # Extract the parameter names and values from each transition.
    dicts = map(ts) do t
        nms, vs = flatten_namedtuple(t.θ)
        for nm in nms
            push!(names_set, nm)
        end
        # Convert the names and values to a single dictionary.
        return Dict(nms[j] => vs[j] for j in 1:length(vs))
    end
    names = collect(names_set)
    vals = [get(dicts[i], key, missing) for i in eachindex(dicts),
        (j, key) in enumerate(names)]

    return names, vals
end

function flatten_namedtuple(nt::NamedTuple)
    names_vals = mapreduce(vcat, keys(nt)) do k
        v = nt[k]
        if length(v) == 1
            return [(string(k), v)]
        else
            return mapreduce(vcat, zip(v[1], v[2])) do (vnval, vn)
                return collect(FlattenIterator(vn, vnval))
            end
        end
    end
    return [vn[1] for vn in names_vals], [vn[2] for vn in names_vals]
end

function get_transition_extras(ts::Vector)
    # Get the extra field names from the sampler state type.
    # This handles things like :lp or :weight.
    extra_params = additional_parameters(eltype(ts))

    # Get the values of the extra parameters.
    local extra_names
    all_vals = []

    # Iterate through each transition.
    for t in ts
        extra_names = String[]
        vals = []

        # Iterate through each of the additional field names
        # in the struct.
        for p in extra_params
            # Check whether the field contains a NamedTuple,
            # in which case we need to iterate through each
            # key/value pair.
            prop = getproperty(t, p)
            if prop isa NamedTuple
                for (k, v) in pairs(prop)
                    push!(extra_names, string(k))
                    push!(vals, v)
                end
            else
                push!(extra_names, string(p))
                push!(vals, prop)
            end
        end
        push!(all_vals, vals)
    end

    # Convert the vector-of-vectors to a matrix.
    valmat = [all_vals[i][j] for i in 1:length(ts), j in 1:length(all_vals[1])]

    return extra_names, valmat
end

# Default Chains constructor.
function AbstractMCMC.bundle_samples(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler,
    N::Integer,
    ts::Vector,
    chain_type::Type{Chains};
    discard_adapt::Bool=true,
    save_state=false,
    kwargs...
)
    # Check if we have adaptation samples.
    if discard_adapt && :n_adapts in fieldnames(typeof(spl.alg))
        ts = ts[(spl.alg.n_adapts+1):end]
    end

    # Convert transitions to array format.
    # Also retrieve the variable names.
    nms, vals = _params_to_array(ts, spl)

    # Get the values of the extra parameters in each Transition struct.
    extra_params, extra_values = get_transition_extras(ts)

    # Extract names & construct param array.
    nms = [nms; extra_params]
    parray = hcat(vals, extra_values)

    # If the state field has average_logevidence or final_logevidence, grab that.
    le = missing
    if :average_logevidence in fieldnames(typeof(spl.state))
        le = getproperty(spl.state, :average_logevidence)
    elseif :final_logevidence in fieldnames(typeof(spl.state))
        le = getproperty(spl.state, :final_logevidence)
    end

    # Check whether to invlink! the varinfo
    if islinked(spl.state.vi, spl)
        invlink!(spl.state.vi, spl)
    end

    # Set up the info tuple.
    if save_state
        info = (range = rng, model = model, spl = spl, vi = spl.state.vi)
    else
        info = NamedTuple()
    end

    # Chain construction.
    return Chains(
        parray,
        string.(nms),
        deepcopy(TURING_INTERNAL_VARS);
        evidence=le,
        info=info,
        sorted=true
    )
end

function AbstractMCMC.bundle_samples(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler,
    N::Integer,
    ts::Vector,
    chain_type::Type{Vector{NamedTuple}};
    discard_adapt::Bool=true,
    save_state=false,
    kwargs...
)
    nts = Vector{NamedTuple}(undef, N)

    for (i,t) in enumerate(ts)
        k = collect(keys(t.θ))
        vs = []
        for v in values(t.θ)
            push!(vs, v[1])
        end

        push!(k, :lp)

        nts[i] = NamedTuple{tuple(k...)}(tuple(vs..., t.lp))
    end

    return map(identity, nts)
end

function save(c::Chains, spl::Sampler, model, vi, samples)
    nt = NamedTuple{(:spl, :model, :vi, :samples)}((spl, model, deepcopy(vi), samples))
    return setinfo(c, merge(nt, c.info))
end

function resume(c::Chains, n_iter::Int; chain_type=Chains, progress=PROGRESS[], kwargs...)
    @assert !isempty(c.info) "[Turing] cannot resume from a chain without state info"

    # Sample a new chain.
    newchain = AbstractMCMC.sample(
        c.info[:range],
        c.info[:model],
        c.info[:spl],
        n_iter;
        resume_from=c,
        reuse_spl_n=n_iter,
        chain_type=Chains,
        progress=progress,
        kwargs...
    )

    # Stick the new samples at the end of the old chain.
    return vcat(c, newchain)
end

function set_resume!(
    s::Sampler;
    resume_from::Union{Chains, Nothing}=nothing,
    kwargs...
)
    # If we're resuming, grab the sampler info.
    if resume_from !== nothing
        s = resume_from.info[:spl]
    end
end

#########################
# Default sampler state #
#########################

"""
A blank `AbstractSamplerState` that contains only `VarInfo` information.
"""
mutable struct SamplerState{VIType<:VarInfo} <: AbstractSamplerState
    vi :: VIType
end

#######################################
# Concrete algorithm implementations. #
#######################################

include("ess.jl")
include("hmc.jl")
include("mh.jl")
include("is.jl")
include("AdvancedSMC.jl")
include("gibbs.jl")
include("../contrib/inference/sghmc.jl")

################
# Typing tools #
################

for alg in (:SMC, :PG, :MH, :IS, :ESS, :Gibbs)
    @eval getspace(::$alg{space}) where {space} = space
end
for alg in (:HMC, :HMCDA, :NUTS, :SGLD, :SGHMC)
    @eval getspace(::$alg{<:Any, space}) where {space} = space
end

floatof(::Type{T}) where {T <: Real} = typeof(one(T)/one(T))
floatof(::Type) = Real # fallback if type inference failed

function get_matching_type(
    spl::AbstractSampler,
    vi::VarInfo,
    ::Type{T},
) where {T}
    return T
end
function get_matching_type(
    spl::AbstractSampler,
    vi::VarInfo,
    ::Type{<:AbstractFloat},
)
    return floatof(eltype(vi, spl))
end
function get_matching_type(
    spl::Sampler{<:Hamiltonian},
    vi::VarInfo,
    ::Type{<:Union{Missing, AbstractFloat}},
)
    return Union{Missing, floatof(eltype(vi, spl))}
end
function get_matching_type(
    spl::Sampler{<:Hamiltonian},
    vi::VarInfo,
    ::Type{<:AbstractFloat},
)
    return floatof(eltype(vi, spl))
end
function get_matching_type(
    spl::Sampler{<:Hamiltonian},
    vi::VarInfo,
    ::Type{TV},
) where {T, N, TV <: Array{T, N}}
    return Array{get_matching_type(spl, vi, T), N}
end
function get_matching_type(
    spl::Sampler{<:Union{PG, SMC}},
    vi::VarInfo,
    ::Type{TV},
) where {T, N, TV <: Array{T, N}}
    return TArray{T, N}
end

##############
# Utilities  #
##############

getspace(spl::Sampler) = getspace(spl.alg)
function ambiguity_error_msg()
    return "Ambiguous `lhs .~ rhs` or `@. lhs ~ rhs` syntax. The broadcasting can either be
    column-wise following the convention of Distributions.jl or element-wise following
    Julia's general broadcasting semantics. Please make sure that the element type of `lhs`
    is not a supertype of the support type of `AbstractVector` to eliminate ambiguity."
end

end # module


text:
using BinaryBuilder, Pkg

name = "MAGMA"
version = v"2.6.2"
sources = [
    GitSource("https://bitbucket.org/icl/magma", "864fa8b28e8c2905a41fbef32f156695378dc104"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/magma/

cp ../make.inc .

make -j${nproc} sparse-lib
make install prefix=${prefix}
install_license COPYRIGHT
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [Platform("x86_64", "linux"; libc="glibc")]

# The products that we will ensure are always built
products = [
    LibraryProduct("libmagma", :libmagma; dont_dlopen=true),
    LibraryProduct("libmagma_sparse", :libmagma_sparse; dont_dlopen=true)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(name="CUDA_full_jll", version=v"11.6.0")),
    Dependency("libblastrampoline_jll")
]

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

max_stars_repo_name: staticfloat/BinaryBuilderRepos

max_issues_repo_name: staticfloat/BinaryBuilderRepos

max_forks_repo_name: staticfloat/BinaryBuilderRepos


text:
export getAcousticPointSource,loc2cs;

function getElasticRhs(Minv,TYPE)
## Generating the right hand side
if Minv.dim==3
	qx = zeros(TYPE,tuple(Minv.n + [1; 0; 0]...));
	qy = zeros(TYPE,tuple(Minv.n + [0; 1; 0]...));
	qz = zeros(TYPE,tuple(Minv.n + [0; 0; 1]...));
	src = [div(Minv.n[1],2),div(Minv.n[2],2),1];
	qz[src[1],src[2],1] = -2.0./(Minv.h[1]);
	qx[src[1],src[2],1] = -1.0./(Minv.h[3]);
	qx[src[1]+1,src[2],1] = +1.0./(Minv.h[3]);
	qy[src[1],src[2],1] = -1.0./(Minv.h[3]);
	qy[src[1],src[2]+1,1] = +1.0./(Minv.h[3]);
	q = vec([qx[:]; qy[:]; qz[:];]);
else
	qx = zeros(TYPE,tuple(Minv.n + [1; 0]...))
	qy = zeros(TYPE,tuple(Minv.n + [0; 1]...))
	qy[div(Minv.n[1],2),1] = -2.0./(Minv.h[1]);
	qx[div(Minv.n[1],2),1] = -1.0./(Minv.h[2]);
	qx[div(Minv.n[1],2)+1,1] = +1.0./(Minv.h[2]);
	# qy[div(Minv.n[1],2),div(Minv.n[2],2)] = -2./(Minv.h[1]);
	# qx[div(Minv.n[1],2),div(Minv.n[2],2)] = -1./(Minv.h[2]);
	# qx[div(Minv.n[1],2)+1,div(Minv.n[2],2)] = +1./(Minv.h[2]);
	q = vec([qx[:]; qy[:]])
end
return q;
end

function getTopPointSrc(Minv)
if Minv.dim==3
	src = [div(Minv.n[1]+1,2);div(Minv.n[2]+1,2);1];
else
    src = [div(Minv.n[1]+1,2);1]
end
return src;
end

function loc2cs(loc1::Int64,loc2::Int64,n::Array{Int64,1})
@inbounds cs = loc1 + (loc2-1)*n[1];
return cs;
end

function loc2cs3D(loc1::Int64,loc2::Int64,loc3::Int64,n::Array{Int64,1})
@inbounds cs = loc1 + (loc2-1)*n[1] + (loc3-1)*n[1]*n[2];
return cs;
end

function loc2cs(n::Array{Int64},sub::Array{Int64})
if length(sub)==2
	return loc2cs(sub[1],sub[2],n);
else
	return loc2cs3D(sub[1],sub[2],sub[3],n);
end
end

function getAcousticPointSource(Minv,TYPE,src = getTopPointSrc(Minv))
## Generating the right hand side
n_nodes = Minv.n .+ 1;
q = zeros(TYPE,tuple(n_nodes...));
q[loc2cs(n_nodes,src)] = 1.0./(norm(Minv.h)^2);
return q,src;
end


text:
include("init.jl")

@testset "Query federate tests"  begin

    broker = createBroker(2)
    vFed1, fedinfo1 = createValueFederate(1, "fed0")
    vFed2, fedinfo2 = createValueFederate(1, "fed1")

    h.helicsFederateRegisterGlobalTypePublication(vFed1, "pub1", "double", "")
    h.helicsFederateRegisterTypePublication(vFed1, "pub2", "double", "")
    h.helicsFederateRegisterTypePublication(vFed2, "pub3", "double", "")
    h.helicsFederateEnterInitializingModeAsync(vFed1)
    h.helicsFederateEnterInitializingMode(vFed2)
    h.helicsFederateEnterInitializingModeComplete(vFed1)

    core = h.helicsFederateGetCoreObject(vFed1)

    q1 = h.helicsCreateQuery("Testfed0", "publications")
    res = h.helicsQueryCoreExecute(q1, core)
    @test res == "[pub1;Testfed0/pub2]"
    # res = h.helicsQueryExecute(q1, vFed2)
    # @test res == "[pub1;Testfed0/pub2]"
    h.helicsQueryFree(q1)

    # q1 = h.helicsCreateQuery("Testfed1", "isinit")
    # res = h.helicsQueryExecute(q1, vFed1)
    # @test res, "true"
    # h.helicsQueryFree(q1)

    # q1 = h.helicsCreateQuery("Testfed1", "publications")
    # res = h.helicsQueryExecute(q1, vFed1)
    # @test res == "[Testfed1/pub3]"
    # h.helicsQueryFree(q1)

    h.helicsCoreFree(core)
    h.helicsFederateFinalizeAsync(vFed1)
    h.helicsFederateFinalize(vFed2)
    h.helicsFederateFinalizeComplete(vFed1)

    destroyFederate(vFed1, fedinfo1)
    destroyFederate(vFed2, fedinfo2)
    destroyBroker(broker)

end

@testset "Query broker tests" begin

    broker = createBroker(2)
    vFed1, fedinfo1 = createValueFederate(1, "fed0")
    vFed2, fedinfo2 = createValueFederate(1, "fed1")
    core = h.helicsFederateGetCoreObject(vFed1)

    q1 = h.helicsCreateQuery("root", "federates")
    res = h.helicsQueryCoreExecute(q1, core)
    name1 = h.helicsFederateGetName(vFed1)
    name2 = h.helicsFederateGetName(vFed2)

    @test "[$name1;$name2]" == res

    res = h.helicsQueryExecute(q1, vFed1)
    @test "[$name1;$name2]" == res

    h.helicsFederateEnterInitializingModeAsync(vFed1)
    h.helicsFederateEnterInitializingMode(vFed2)
    h.helicsFederateEnterInitializingModeComplete(vFed1)
    h.helicsQueryFree(q1)
    h.helicsCoreFree(core)
    h.helicsFederateFinalizeAsync(vFed1)
    h.helicsFederateFinalize(vFed2)
    h.helicsFederateFinalizeComplete(vFed1)

    destroyFederate(vFed1, fedinfo1)
    destroyFederate(vFed2, fedinfo2)
    destroyBroker(broker)
end


text:
module test_julia_tuples_ntuple

using Test

@test tuple() == ()
@test tuple(1,) == (1,)

@test tuple(1:3...) == ntuple(identity, 3) == (1, 2, 3)
@test tuple([1, 2, 3]) == ([1, 2, 3],)

@test ntuple(M -> 2, 1) == (2,)
@test ntuple(M -> 2, 3) == (2, 2, 2)
@test ntuple(M -> 2, Val(1)) == (2,)
@test ntuple(M -> 2, Val(3)) == (2, 2, 2)

end # module test_julia_tuples

module test_julia_tuples_front_tail

using Test

@test Base.front((1,)) == ()
@test Base.front((1,2,3)) == (1, 2)
@test Base.tail((1,)) == ()
@test Base.tail((1,2,3)) == (2, 3)

end # module test_julia_tuples_front_tail


text:
mutable struct Results <: IS.Results
    variables::Dict
    total_cost::Dict
    optimizer_log::Dict
    time_stamp::DataFrames.DataFrame
end

get_variables(r::Results) = r.variables
get_total_cost(r::Results) = r.total_cost
get_optimizer_log(r::Results) = r.optimizer_log
get_time_stamp(r::Results) = r.time_stamp

struct StackedArea
    time_range::Array
    data_matrix::Matrix
    labels::Array
end

struct BarPlot
    time_range::Array
    bar_data::Matrix
    labels::Array
end

struct StackedGeneration
    time_range::Array
    data_matrix::Matrix
    labels::Array
end

struct BarGeneration
    time_range::Array
    bar_data::Matrix
    labels::Array
end

"""
    variable = get_stacked_plot_data(res::IS.Results, variable::String)

This function takes in results and uses a dataframe from whichever variable name was given and converts it to type StackedArea.
StackedArea is the type of struct that signals the plot() function to use the StackedArea plot recipe method.

# Arguments
- `res::IS.Results`: results
- `variable::String`: the variable to be plotted

#Example
```julia
ThermalStandard = get_stacked_plot_data(res, "P_ThermalStandard")
plot(ThermalStandard)
```

# Accepted Key Words
- `sort::Array`: the array of generators to be plotted, in the order to be plotted
"""

function get_stacked_plot_data(res::IS.Results, variable::String; kwargs...)

    sort = get(kwargs, :sort, nothing)
    time_range = res.time_stamp[!, :Range]
    variable = res.variables[Symbol(variable)]
    alphabetical = sort!(names(variable))

    if isnothing(sort)
        variable = variable[:, alphabetical]
    else
        variable = variable[:, sort]
    end

    data_matrix = convert(Matrix, variable)
    labels = collect(names(variable))
    legend = [names(variable)[1]]
    for name in 2:length(labels)
        legend = hcat(legend, string.(labels[name]))
    end

    return StackedArea(time_range, data_matrix, legend)

end

"""
    variable = get_bar_plot_data(res::IS.Results, variable::String)

This function takes in results and uses a dataframe from whichever variable name was given and converts it to type BarPlot.
StackedGeneration is the type of struct that signals the plot() function to use the StackedGeneration plot recipe method.

# Arguments
- `res::IS.Results`: results
- `variable::String`: the variable to be plotted

#Example
```julia
ThermalStandard = get_bar_plot_data(res, "P_ThermalStandard")
plot(ThermalStandard)
```

# Accepted Key Words
- `sort::Array`: the array of generators to be plotted, in the order to be plotted
"""

function get_bar_plot_data(res::IS.Results, variable::String; kwargs...)

    sort = get(kwargs, :sort, nothing)
    time_range = res.time_stamp[!, :Range]
    variable = res.variables[Symbol(variable)]
    alphabetical = sort!(names(variable))

    if isnothing(sort)
        variable = variable[:, alphabetical]
    else
        variable = variable[:, sort]
    end

    data = convert(Matrix, variable)
    bar_data = sum(data, dims = 1)
    labels = collect(names(variable))
    legend = [names(variable)[1]]
    for name in 2:length(labels)
        legend = hcat(legend, string.(labels[name]))
    end

    return BarPlot(time_range, bar_data, legend)
end

"""
    variable = get_stacked_gen_data(res::IS.Results)

This function takes in results and stacks the variables given.
StackedGeneration is the type of struct that signals the plot() function to use the StackedGeneration plot recipe method.

# Arguments
- `res::IS.Results`: results

#Example
```julia
stack = get_stacked_gen_data(res)
plot(stack)
```

# Accepted Key Words
- `sort::Array`: the array of generators to be plotted, in the order to be plotted
"""

function get_stacked_generation_data(res::IS.Results; kwargs...)

    sort = get(kwargs, :sort, nothing)
    time_range = res.time_stamp[!, :Range]
    key_name = collect(keys(res.variables))
    alphabetical = sort!(key_name)

    if !isnothing(sort)
        labels = sort
    else
        labels = alphabetical
    end

    variable = res.variables[Symbol(labels[1])]
    data_matrix = sum(convert(Matrix, variable), dims = 2)
    legend = [key_name[1]]

    for i in 1:length(labels)
        if i !== 1
            variable = res.variables[Symbol(labels[i])]
            legend = hcat(legend, string.(key_name[i]))
            data_matrix = hcat(data_matrix, sum(convert(Matrix, variable), dims = 2))
        end
    end

    return StackedGeneration(time_range, data_matrix, legend)

end

"""
    variable = get_bar_gen_data(res::IS.Results)

This function takes in results and stacks the variables given.
StackedGeneration is the type of struct that signals the plot() function to use the StackedGeneration plot recipe method.

# Arguments
- `res::IS.Results`: results

#Example
```julia
stack = get_stacked_gen_data(res)
plot(stack)
```

# Accepted Key Words
- `sort::Array`: the array of generators to be plotted, in the order to be plotted
"""

function get_bar_gen_data(res::IS.Results)

    time_range = res.time_stamp[!, :Range]
    key_name = collect(keys(res.variables))
    variable = res.variables[Symbol(key_name[1])]
    data_matrix = sum(convert(Matrix, variable), dims = 2)
    legend = [key_name[1]]
    for i in 1:length(key_name)
        if i !== 1
            variable = res.variables[Symbol(key_name[i])]
            legend = hcat(legend, string.(key_name[i]))
            data_matrix = hcat(data_matrix, sum(convert(Matrix, variable), dims = 2))
        end
    end
    bar_data = sum(data_matrix, dims = 1)
    return BarGeneration(time_range, bar_data, legend)

end

"""
    sort_data(results::IS.Results)

This function takes in struct Results, sorts the generators in each variable, and outputs the sorted
results. The generic function sorts the generators alphabetically.

# Arguments
- `results::Results`: the results of the simulation

# Key Words
- `Variables::Dict{Symbol, Array{Symbol}`: the desired variables and their generator order

#Examples
```julia
Variables = Dict(:ON_ThermalStandard => [:Brighton, :Solitude])
sorted_results = sort_data(res_UC; Variables = Variables)
```
***Note:*** only the generators included in key word 'Variables' will be in the
results, and consequently, only these will be plotted.
"""
function sort_data(res::IS.Results; kwargs...)
    order = get(kwargs, :Variables, Dict())
    if !isempty(order)
        labels = collect(keys(order))
    else
        labels = sort!(collect(keys(res.variables)))
    end
    sorted_variables = Dict()
    for label in labels
        sorted_variables[label] = res.variables[label]
    end
    for (k, variable) in sorted_variables
        if !isempty(order)
            variable = variable[:, order[k]]
        else
            alphabetical = sort!(names(variable))
            variable = variable[:, alphabetical]
        end
        sorted_variables[k] = variable
    end
    return Results(sorted_variables, res.total_cost, res.optimizer_log, res.time_stamp)
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.1.6
#   kernelspec:
#     display_name: Julia 1.1.1
#     language: julia
#     name: julia-1.1
# ---

using LinearAlgebra, Distributions, Statistics, Plots
using NearestNeighbors

#Initializing noisy non linear data
x = range(0,stop=1,length=100) |> collect
y = sin.(x .* 3π/2 )
σ = 0.1
d = MvNormal([0], σ .* Matrix(I,1,1))
n = length(x)
y_noise = y .+ rand(d,n)';

scatter( x, y_noise)
plot!(x, y)

# +
"""
    lowess_knn(x, y, tau = .005) -> yest

Locally weighted regression: fits a nonparametric regression curve
to a scatterplot.
The arrays x and y contain an equal number of elements; each pair
(x[i], y[i]) defines a data point in the scatterplot. The function
returns the estimated (smooth) values of y.
The kernel function is the bell shaped function with parameter tau.
Larger tau will result in a smoother curve.
"""
function lowess(x, k = 3)

    m = last(size(x))
    kdt = KDTree(x; leafsize = 10)
    idxs, dists = knn(kdt, x, k, true)
    yest = zeros(Float64, m)

    w = (1 .- dists.^3) .^ 3
    weights = zeros(Float64, m)

    #Looping through all x-points
    for i in 1:m
        weights .= w[i]
        b = [sum(weights .* y), sum(weights .* y .* x)]
        A = SymTridiagonal(
                [sum(weights), sum(weights .* x .* x)],
                [sum(weights .* x)])
        theta = A \ b
        yest[i] = theta[1] + theta[2] * x[i]
    end

    yest

end
# -


text:
using DrWatson, Test, Dates

ps = DrWatson.PATH_SEPARATOR

d = (a = 0.153456453, b = 5, mode = "double")
@test savename(d; digits = 4) == "a=0.1535_b=5_mode=double"
@test savename("n", d) == "n_a=0.153_b=5_mode=double"
@test savename("n$(ps)", d) == "n$(ps)a=0.153_b=5_mode=double"
@test savename("n$(ps)", d; connector = "-") == joinpath("n", "a=0.153-b=5-mode=double")
@test savename(d, "n") == "a=0.153_b=5_mode=double.n"
@test savename("n", d, "n") == "n_a=0.153_b=5_mode=double.n"
@test savename("n", d, "n"; connector = "-") == "n-a=0.153-b=5-mode=double.n"
@test savename(d, allowedtypes = (String,)) == "mode=double"
@test savename(d, connector=" | ", equals=" = ") == "a = 0.153 | b = 5 | mode = double"

tday = today()
@test savename(@dict(tday)) == "tday=$(string(tday))"

rick = (never = "gonna", give = "you", up = "!");
@test savename(rick) == "give=you_never=gonna_up=!"
@test savename(rick; ignores = ["up"]) == "give=you_never=gonna"
@test savename(rick; sort = false) == "never=gonna_give=you_up=!"

x = 3; y = 5.0;
d = Dict(:x => x, :y => y)
n = (x = x, y = y)

@test d == @dict x y
@test Dict("x" => x, "y" => y) == @strdict x y
@test n == @ntuple x y
@test (@savename x y) == savename(d)

z = "lala"
d2 = Dict(:x => x, :y => y, :z => z)
n2 = (x = x, y = y, z = z)

@test d2 == @dict x y z
@test n2 == @ntuple x y z

@test savename(n2; ignores=(:y,))  == "x=3_z=lala"
@test savename(n2; ignores=("y",)) == "x=3_z=lala"
@test savename(n2; accesses=(:x, :y), ignores=(:y,)) == "x=3"

@test savename(@dict x y) == "x=3_y=5.0"
@test savename(@ntuple x y) == "x=3_y=5.0"
w = rand(50)
@test savename(@dict x y w) == savename(@dict x y)
@test (@savename x y w) == savename(@dict x y)
@test savename(@ntuple x y w) == savename(@dict x y)

@test ntuple2dict(@ntuple x y) == @dict x y
@test sort(collect(keys(dict2ntuple(@dict x y)))) == sort(collect(keys(@ntuple x y)))

@test keytype(tostringdict(d2)) == String
@test keytype(tosymboldict(tostringdict(d2))) == Symbol

a = 3; b = 4
c = @ntuple a b
d = 5; e = @dict c d

@test DrWatson.access(e, :c, :a) == a
ff = dict2ntuple(e)
@test DrWatson.access(ff, :c, :a) == a
@test ff.c.a == a

# Expand tests:
a = 3; b = 4
c = @dict a b
d = 5; e = @dict c d

s = savename(e; allowedtypes = (Any,), expand = ["c"])
@test '(' ∈ s
@test ')' ∈ s
@test occursin("a=3", s)
@test occursin("b=4", s)

# More detailed expand tests
struct A
    a
    p
end
DrWatson.default_allowed(::A) = (Any,)
DrWatson.default_expand(::A) = ["p"]

x = A(5, (b = 3, c = 4))

@test savename(x) == "a=5_p=(b=3,c=4)"

# empty container
x = A(5, NamedTuple())
@test !occursin("p", savename(x))

# container with values that are not by default printed in savename
x = A(5, (m = rand(50,50),))
@test !occursin("p", savename(x))

# Scientific notation for savename
a = 1.2345e-7
c = 1
d = "test"
di = @dict a c d

@test savename(di,sigdigits=6) == "a=1.2345e-7_c=1_d=test"
@test savename(di,sigdigits=5) == "a=1.2345e-7_c=1_d=test"
@test savename(di,sigdigits=4) == "a=1.234e-7_c=1_d=test"
@test savename(di,sigdigits=3) == "a=1.23e-7_c=1_d=test"
@test savename(di,sigdigits=2) == "a=1.2e-7_c=1_d=test"
@test savename(di,sigdigits=1) == "a=1e-7_c=1_d=test"
@test savename(di) == "a=0.0_c=1_d=test"

sn = savename(di,sigdigits=4)
_,parsed,_ = parse_savename(sn)
@test parsed["a"] == 1.234e-7

sn = savename(di,sigdigits=1)
_,parsed,_ = parse_savename(sn)
@test parsed["a"] == 1.0e-7

# Test for NaN and Inf compatibility
let
    a = Inf
    b = NaN
    di = @dict a b
    @test savename(di) == "a=Inf_b=NaN"
end

# Dedicated Macro tests
@testset "Macro Tests" begin
    x = 3; y = 5.0; z = 42;

    @test Dict(:x=>x, :y=>y, :z=>z) == @dict x y z
    @test Dict(:x=>x, :b=>y, :z=>z) == @dict x b=y z
    @test Dict(:a=>x, :b=>y, :c=>z) == @dict a=x b=y c=z
    @test Dict(:z=>x, :x=>y, :y=>z) == @dict z=x x=y y=z

    @test Dict("x"=>x, "y"=>y, "z"=>z) == @strdict x y z
    @test Dict("x"=>x, "b"=>y, "z"=>z) == @strdict x b=y z
    @test Dict("a"=>x, "b"=>y, "c"=>z) == @strdict a=x b=y c=z
    @test Dict("z"=>x, "x"=>y, "y"=>z) == @strdict z=x x=y y=z

    @test (x=x, y=y, z=z) == @ntuple x y z
    @test (x=x, b=y, z=z) == @ntuple x b=y z
    @test (a=x, b=y, c=z) == @ntuple a=x b=y c=z
    @test (z=x, x=y, y=z) == @ntuple z=x x=y y=z

    @test savename((; x, y, z)) == @savename x y z
    @test savename((; x, b=y, z)) == @savename x b=y z
    @test savename((; a=x, b=y, c=z)) == @savename a=x b=y c=z
    @test savename((; z=x, x=y, y=z)) == @savename z=x x=y y=z

end


text:
module MyModel

using Epicenter

# expose variables to users of your model
export  model,
        curr_year

# expose methods to users of your model
export  init,
        advance,
        reset

using InnerModule

# -------

type ModelResultType
    num_customers::Int
    revenue::Float64

    ModelResultType(; num_customers = 0, revenue = 0) = new(num_customers, revenue)
end

# -------

type ModelType
    price::Float64
    results_over_time::Array{ModelResultType, 1}

    ModelType(; price = start_price) = new(price, ModelResultType[])
end

# -------

const global start_year = 2014
const global initial_customers = 50
const global start_price = 15

global curr_year
global model

# -------

function init()
    reset()
end

function advance()
    global curr_year += 1

    current_result = ModelResultType()
    current_result.num_customers = int(cos(model.price / 12) * 64)
    current_result.revenue = model.price * current_result.num_customers

    # add current_result, the results of executing the model for this time step,
    # to the end of the model.results array
    push!(model.results_over_time, current_result)

    # length(model.results_over_time) is the index of the last element of the
    # model.results_over_time array. So adding model.results[length(model.results)]
    # to the persistence queue signals that only the last element of the array has
    # been added or changed. The Epicenter platform will check and process this
    # queue approximately every second.
    record(:model, :results_over_time, length(model.results_over_time))
end

function reset()
    global curr_year = start_year
    global model = ModelType()

    current_result = ModelResultType()
    current_result.num_customers = initial_customers
    current_result.revenue = model.price * current_result.num_customers

    push!(model.results_over_time, current_result)

    # We want to record the whole array here! If it was full of runs before, Epicenter needs to know
    # that they were all cleared and that there's only one result now
    record(:model, :results_over_time)
end

end


text:
abstract type Geometry end

struct IntersectResult{T<:Geometry}
    ray::Ray3
    geometry::T
    distance::Float64
    position::Vector3
    normal::Vector3
end

not_hit(g::T) where T <: Geometry = Nullable{IntersectResult{T}}()
hit(ir::Nullable{IntersectResult{T}}) where T = !isnull(ir)
replace_geometry(result::IntersectResult, geometry::Geometry) =
    IntersectResult(result.ray, geometry, result.distance, result.position, result.normal)

struct GeometryWithMaterial{T <: Geometry, V <: Material} <: Geometry
    geometry::T
    material::V
end

function ∩(ray::Ray3, gwm::GeometryWithMaterial)
    result = ray ∩ gwm.geometry
    if !hit(result)
        return result
    end
    Nullable(replace_geometry(result.value, gwm))
end

struct UnionGeometry <: Geometry
    geometries::Vector{Geometry}
end

function ∩(ray::Ray3, ug::UnionGeometry)
    result = not_hit(ug)
    for geometry in ug.geometries
        tmp = ray ∩ geometry
        if hit(tmp) && (!hit(result) || tmp.value.distance < result.value.distance)
            result = tmp
        end
    end
    result
end

struct Sphere <: Geometry
    center::Vector3
    radius::Float64
end

function ∩(ray::Ray3, s::Sphere)
    𝐯 = ray.origin - s.center
    a0 = norm²(𝐯) - s.radius^2
    ddotv = ray.direction ⋅ 𝐯
    if ddotv <= 0
        discr = ddotv^2 - a0
        if discr >= 0
            distance = -ddotv - √discr
            position = get_point(ray, distance)
            normal = normalize(position - s.center)
            return Nullable(IntersectResult(ray, s, distance, position, normal))
        end
    end
    return not_hit(s)
end

struct Plane <: Geometry
    normal::Vector3
    distance::Float64
end

function ∩(ray::Ray3, s::Plane)
     ldotn = ray.direction ⋅ s.normal
     if ldotn >= 0
        return not_hit(s)
     end
     xdotn = ray.origin ⋅ s.normal
     distance = (s.distance - xdotn) / ldotn
     position = get_point(ray, distance)
     normal = s.normal
     return Nullable(IntersectResult(ray, s, distance, position, normal))
end


text:
module LM

using Merlin
using ..LightNLP

include("data.jl")
include("model.jl")
include("nn.jl")

end


text:
# This file is part of the TaylorSeries.jl Julia package, MIT license

# Printing of TaylorSeries objects

# subscriptify is taken from the ValidatedNumerics.jl package, licensed under MIT "Expat".
# superscriptify is a small variation

const subscript_digits = [c for c in "₀₁₂₃₄₅₆₇₈₉"]
const superscript_digits = [c for c in "⁰¹²³⁴⁵⁶⁷⁸⁹"]

function subscriptify(n::Int)
    dig = reverse(digits(n))
    join([subscript_digits[i+1] for i in dig])
end

function superscriptify(n::Int)
    dig = reverse(digits(n))
    join([superscript_digits[i+1] for i in dig])
end

function pretty_print{T<:Number}(a::Taylor1{T})
    z = zero(T)
    @compat space = string(" ")
    bigO = string("+ 𝒪(t", superscriptify(a.order+1), ")")
    a == zero(a) && return string(space, z, space, bigO)
    @compat strout::String = space
    ifirst = true
    for i in eachindex(a.coeffs)
        @compat monom::String = i==1 ? string("") : i==2 ? string(" t") :
            string(" t", superscriptify(i-1))
        @inbounds c = a.coeffs[i]
        c == z && continue
        cadena = numbr2str(c, ifirst)
        strout = string(strout, cadena, monom, space)
        ifirst = false
    end
    strout = strout * bigO
    strout
end
function pretty_print{T<:Number}(a::HomogeneousPolynomial{T})
    z = zero(T)
    @compat space = string(" ")
    a == zero(a) && return string(space, z)
    @compat strout::String = homogPol2str(a)
    strout
end
function pretty_print{T<:Number}(a::TaylorN{T})
    z = zero(T)
    @compat space = string(" ")
    @compat bigO::String  = string(" + 𝒪(‖x‖", superscriptify(a.order+1), ")")
    a == zero(a) && return string(space, z, bigO)
    @compat strout::String = string("")
    ifirst = true
    for ord in eachindex(a.coeffs)
        pol = a.coeffs[ord]
        pol == zero(a.coeffs[ord]) && continue
        @compat cadena::String = homogPol2str( pol )
        strsgn = (ifirst || ord == 1 || cadena[2] == '-') ? string("") : string(" +")
        strout = string( strout, strsgn, cadena)
        ifirst = false
    end
    strout = strout * bigO
    strout
end

function homogPol2str{T<:Number}(a::HomogeneousPolynomial{T})
    numVars = get_numvars()
    order = a.order
    z = zero(T)
    @compat space = string(" ")
    @compat strout::String = space
    ifirst = true
    iIndices = zeros(Int, numVars)
    for pos = 1:size_table[order+1]
        @compat monom::String = string("")
        @inbounds iIndices[:] = coeff_table[order+1][pos]
        for ivar = 1:numVars
            powivar = iIndices[ivar]
            if powivar == 1
                monom = string(monom, name_taylorNvar(ivar))
            elseif powivar > 1
                monom = string(monom, name_taylorNvar(ivar), superscriptify(powivar))
            end
        end
        @inbounds c = a.coeffs[pos]
        c == z && continue
        cadena = numbr2str(c, ifirst)
        strout = string(strout, cadena, monom, space)
        ifirst = false
    end
    return strout[1:end-1]
end
function numbr2str(zz, ifirst::Bool=false)
    zz == zero(zz) && return string( zz )
    plusmin = ifelse( ifirst, string(""), string("+ ") )
    return string(plusmin, zz)
end
function numbr2str{T<:Union{AbstractFloat,Integer,Rational}}(zz::T, ifirst::Bool=false)
    zz == zero(T) && return string( zz )
    plusmin = ifelse( zz < zero(T), string("- "),
                ifelse( ifirst, string(""), string("+ ")) )
    return string(plusmin, abs(zz))
end
function numbr2str{T<:Complex}(zz::T, ifirst::Bool=false)
    zT = zero(zz.re)
    zz == zero(zz) && return string(zT)
    zre, zim = reim(zz)
    cadena = string("")
    if zre > zT
        if ifirst
            cadena = string(" ( ", abs(zre)," ")
        else
            cadena = string(" + ( ", abs(zre)," ")
        end
        if zim > zT
            cadena = string(cadena, "+ ", abs(zim), " im )")
        elseif zim < zT
            cadena = string(cadena, "- ", abs(zim), " im )")
        else
            cadena = string(cadena, ")")
        end
    elseif zre < zT
        cadena = string(" - ( ", abs(zre), " ")
        if zim > zT
            cadena = string(cadena, "- ", abs(zim), " im )")
        elseif zim < zT
            cadena = string(cadena, "+ ", abs(zim), " im )")
        else
            cadena = string(cadena, ")")
        end
    else
        if zim > zT
            if ifirst
                cadena = string("( ", abs(zim), " im )")
            else
                cadena = string("+ ( ", abs(zim), " im )")
            end
        else
            cadena = string("- ( ", abs(zim), " im )")
        end
    end
    return cadena
end

#name_taylorNvar(n::Int) = string("⋅x", subscriptify(n))
name_taylorNvar(i::Int) = string(" ", get_variable_names()[i])

# summary
summary{T<:Number}(a::Taylor1{T}) = string(a.order, "-order ", typeof(a), ":")
function summary{T<:Number}(a::Union{HomogeneousPolynomial{T}, TaylorN{T}})
    string(a.order, "-order ", typeof(a), " in ", get_numvars(), " variables:")
end

# show
function show(io::IO, a::Union{Taylor1, HomogeneousPolynomial, TaylorN})
    # (isa(a, TaylorN) || isa(a, Taylor1)) && println(io, summary(a))
    print(io, pretty_print(a))
end


text:
using Printf
using LinearAlgebra
using LinearAlgebra: Givens, Rotation

# Auxiliary
function adiagmax(A::StridedMatrix)
    adm = zero(typeof(real(A[1])))
    @inbounds begin
        for i = size(A,1)
            adm = max(adm, abs(A[i,i]))
        end
    end
    return adm
end

# Hessenberg Matrix
struct HessenbergMatrix{T,S<:StridedMatrix} <: AbstractMatrix{T}
    data::S
end

Base.copy(H::HessenbergMatrix{T,S}) where {T,S} = HessenbergMatrix{T,S}(copy(H.data))

Base.getindex(H::HessenbergMatrix{T,S}, i::Integer, j::Integer) where {T,S} = i > j + 1 ? zero(T) : H.data[i,j]

Base.size(H::HessenbergMatrix) = size(H.data)
Base.size(H::HessenbergMatrix, i::Integer) = size(H.data, i)

function LinearAlgebra.ldiv!(H::HessenbergMatrix, B::AbstractVecOrMat)
    n = size(H, 1)
    Hd = H.data
    for i = 1:n-1
        G, _ = givens!(Hd, i, i+1, i)
        lmul!(G, view(Hd, 1:n, i+1:n))
        lmul!(G, B)
    end
    ldiv!(Triangular(Hd, :U), B)
end

# Hessenberg factorization
struct HessenbergFactorization{T, S<:StridedMatrix,U} <: Factorization{T}
    data::S
    τ::Vector{U}
end

function _hessenberg!(A::StridedMatrix{T}) where T
    n = LinearAlgebra.checksquare(A)
    τ = Vector{Householder{T}}(undef, n - 1)
    for i = 1:n - 1
        xi = view(A, i + 1:n, i)
        t  = LinearAlgebra.reflector!(xi)
        H  = Householder{T,typeof(xi)}(view(xi, 2:n - i), t)
        τ[i] = H
        lmul!(H', view(A, i + 1:n, i + 1:n))
        rmul!(view(A, :, i + 1:n), H)
    end
    return HessenbergFactorization{T, typeof(A), eltype(τ)}(A, τ)
end
LinearAlgebra.hessenberg!(A::StridedMatrix) = _hessenberg!(A)

Base.size(H::HessenbergFactorization, args...) = size(H.data, args...)

# Schur
struct Schur{T,S<:StridedMatrix} <: Factorization{T}
    data::S
    R::Rotation
end

function Base.getproperty(F::Schur, s::Symbol)
    if s === :T || s === :Schur
        return triu(F.data, -1)
    else
        return getfield(F, s)
    end
end

function wilkinson(Hmm, t, d)
    λ1 = (t + sqrt(t*t - 4d))/2
    λ2 = (t - sqrt(t*t - 4d))/2
    return ifelse(abs(Hmm - λ1) < abs(Hmm - λ2), λ1, λ2)
end

# We currently absorb extra unsupported keywords in kwargs. These could e.g. be scale and permute. Do we want to check that these are false?
function _schur!(
    H::HessenbergFactorization{T};
    tol = eps(real(T)),
    debug = false,
    shiftmethod = :Francis,
    maxiter = 30*size(H, 1),
    kwargs...) where T

    n = size(H, 1)
    istart = 1
    iend = n
    HH = H.data
    τ = Rotation(Givens{T}[])

    # iteration count
    i = 0

    @inbounds while true
        i += 1
        debug && println("iteration: $i")
        if i > maxiter
            throw(ArgumentError("iteration limit $maxiter reached"))
        end

        # Determine if the matrix splits. Find lowest positioned subdiagonal "zero"
        for _istart in iend - 1:-1:1
            if abs(HH[_istart + 1, _istart]) <= tol*(abs(HH[_istart, _istart]) + abs(HH[_istart + 1, _istart + 1]))
                # Check if subdiagonal element H[i+1,i] is "zero" such that we can split the matrix

                # Set istart to the beginning of the second block
                istart = _istart + 1

                debug && @printf("Split! Subdiagonal element is: %10.3e and istart now %6d\n", HH[istart, istart - 1], istart)

                # Set the subdiagonal element to zero to signal that a split has taken place
                HH[istart, istart - 1] = 0
                break
            end

            # If no splits have taken place then we start from the top
            istart = 1
        end

        # if block size is one we deflate
        if istart >= iend
            debug && @printf("Bottom deflation! Block size is one. New iend is %6d\n", iend - 1)
            iend -= 1

        # and the same for a 2x2 block
        elseif istart + 1 == iend
            debug && @printf("Bottom deflation! Block size is two. New iend is %6d\n", iend - 2)
            iend -= 2

        # run a QR iteration
        # shift method is specified with shiftmethod kw argument
        else
            Hmm = HH[iend, iend]
            Hm1m1 = HH[iend - 1, iend - 1]
            if iszero(i % 10)
                # Use Eispack exceptional shift
                β = abs(HH[iend, iend - 1]) + abs(HH[iend - 1, iend - 2])
                d = (Hmm + β)^2 - Hmm*β/2
                t = 2*Hmm + 3*β/2
            else
                d = Hm1m1*Hmm - HH[iend, iend - 1]*HH[iend - 1, iend]
                t = Hm1m1 + Hmm
            end
            debug && @printf("block start is: %6d, block end is: %6d, d: %10.3e, t: %10.3e\n", istart, iend, d, t)

            if shiftmethod == :Francis
                debug && @printf("Francis double shift! Subdiagonal is: %10.3e, last subdiagonal is: %10.3e\n", HH[iend, iend - 1], HH[iend - 1, iend - 2])
                doubleShiftQR!(HH, τ, t, d, istart, iend)
            elseif shiftmethod == :Rayleigh
                debug && @printf("Single shift with Rayleigh shift! Subdiagonal is: %10.3e\n", HH[iend, iend - 1])

                # Run a bulge chase
                singleShiftQR!(HH, τ, Hmm, istart, iend)
            else
                throw(ArgumentError("only support supported shift methods are :Francis (default) and :Rayleigh. You supplied $shiftmethod"))
            end
        end
        if iend <= 2 break end
    end

    return Schur{T,typeof(HH)}(HH, τ)
end
_schur!(A::StridedMatrix; kwargs...) = _schur!(_hessenberg!(A); kwargs...)
LinearAlgebra.schur!(A::StridedMatrix; kwargs...) = _schur!(A; kwargs...)

function singleShiftQR!(HH::StridedMatrix, τ::Rotation, shift::Number, istart::Integer, iend::Integer)
    m = size(HH, 1)
    H11 = HH[istart, istart]
    H21 = HH[istart + 1, istart]
    if m > istart + 1
        Htmp = HH[istart + 2, istart]
        HH[istart + 2, istart] = 0
    end
    G, _ = givens(H11 - shift, H21, istart, istart + 1)
    lmul!(G, view(HH, :, istart:m))
    rmul!(view(HH, 1:min(istart + 2, iend), :), G')
    lmul!(G, τ)
    for i = istart:iend - 2
        G, _ = givens(HH[i + 1, i], HH[i + 2, i], i + 1, i + 2)
        lmul!(G, view(HH, :, i:m))
        HH[i + 2, i] = Htmp
        if i < iend - 2
            Htmp = HH[i + 3, i + 1]
            HH[i + 3, i + 1] = 0
        end
        rmul!(view(HH, 1:min(i + 3, iend), :), G')
        # mul!(G, τ)
    end
    return HH
end

function doubleShiftQR!(HH::StridedMatrix, τ::Rotation, shiftTrace::Number, shiftDeterminant::Number, istart::Integer, iend::Integer)
    m = size(HH, 1)
    H11 = HH[istart, istart]
    H21 = HH[istart + 1, istart]
    Htmp11 = HH[istart + 2, istart]
    HH[istart + 2, istart] = 0
    if istart + 3 <= m
        Htmp21 = HH[istart + 3, istart]
        HH[istart + 3, istart] = 0
        Htmp22 = HH[istart + 3, istart + 1]
        HH[istart + 3, istart + 1] = 0
    else
        # values doen't matter in this case but variables should be initialized
        Htmp21 = Htmp22 = Htmp11
    end
    G1, r = givens(H11*H11 + HH[istart, istart + 1]*H21 - shiftTrace*H11 + shiftDeterminant, H21*(H11 + HH[istart + 1, istart + 1] - shiftTrace), istart, istart + 1)
    G2, _ = givens(r, H21*HH[istart + 2, istart + 1], istart, istart + 2)
    vHH = view(HH, :, istart:m)
    lmul!(G1, vHH)
    lmul!(G2, vHH)
    vHH = view(HH, 1:min(istart + 3, m), :)
    rmul!(vHH, G1')
    rmul!(vHH, G2')
    lmul!(G1, τ)
    lmul!(G2, τ)
    for i = istart:iend - 2
        for j = 1:2
            if i + j + 1 > iend break end
            # G, _ = givens(H.H,i+1,i+j+1,i)
            G, _ = givens(HH[i + 1, i], HH[i + j + 1, i], i + 1, i + j + 1)
            lmul!(G, view(HH, :, i:m))
            HH[i + j + 1, i] = Htmp11
            Htmp11 = Htmp21
            # if i + j + 2 <= iend
                # Htmp21 = HH[i + j + 2, i + 1]
                # HH[i + j + 2, i + 1] = 0
            # end
            if i + 4 <= iend
                Htmp22 = HH[i + 4, i + j]
                HH[i + 4, i + j] = 0
            end
            rmul!(view(HH, 1:min(i + j + 2, iend), :), G')
            # mul!(G, τ)
        end
    end
    return HH
end

_eigvals!(A::StridedMatrix; kwargs...)           = _eigvals!(_schur!(A; kwargs...))
_eigvals!(H::HessenbergMatrix; kwargs...)        = _eigvals!(_schur!(H; kwargs...))
_eigvals!(H::HessenbergFactorization; kwargs...) = _eigvals!(_schur!(H; kwargs...))

# Overload methods from LinearAlgebra to make them work generically
if VERSION > v"1.2.0-DEV.0"
    function LinearAlgebra.eigvals!(
        A::StridedMatrix;
        sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
        kwargs...)

        if ishermitian(A)
            return LinearAlgebra.sorteig!(eigvals!(Hermitian(A)), sortby)
        end
        LinearAlgebra.sorteig!(_eigvals!(A; kwargs...), sortby)
    end

    LinearAlgebra.eigvals!(
        H::HessenbergMatrix;
        sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
        kwargs...) = LinearAlgebra.sorteig!(_eigvals!(H; kwargs...), sortby)

    LinearAlgebra.eigvals!(
        H::HessenbergFactorization;
        sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
        kwargs...) = LinearAlgebra.sorteig!(_eigvals!(H; kwargs...), sortby)
else
    function LinearAlgebra.eigvals!(A::StridedMatrix; kwargs...)
        if ishermitian(A)
            return eigvals!(Hermitian(A))
        else
            return _eigvals!(A; kwargs...)
        end
    end

    LinearAlgebra.eigvals!(H::HessenbergMatrix; kwargs...) = _eigvals!(H; kwargs...)

    LinearAlgebra.eigvals!(H::HessenbergFactorization; kwargs...) = _eigvals!(H; kwargs...)
end

# To compute the eigenvalue of the pseudo triangular Schur matrix we just return
# the values of the 1x1 diagonal blocks and compute the eigenvalues of the 2x2
# blocks on the fly. We can check for exact zeros in the subdiagonal because
# the QR algorithm sets the elements to zero when a split happens
function _eigvals!(S::Schur{T}; tol = eps(real(T))) where T
    HH = S.data
    n = size(HH, 1)
    vals = Vector{complex(T)}(undef, n)
    i = 1
    while i < n
        Hii = HH[i, i]
        if iszero(HH[i + 1, i])
            # 1x1 block
            vals[i] = Hii
            i += 1
        else
            # 2x2 block
            Hi1i1 = HH[i + 1, i + 1]
            d = Hii*Hi1i1 - HH[i, i + 1]*HH[i + 1, i]
            t = Hii + Hi1i1
            x = 0.5*t
            y = sqrt(complex(x*x - d))
            vals[i    ] = x + y
            vals[i + 1] = x - y
            i += 2
        end
    end
    if i == n
        vals[i] = HH[n, n]
    end
    return vals
end


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#

# - `_show_default[1] = false`: switch pretty printing on,
# - `_show_default[1] = true`: pretty printing off, Julia default on.
const _show_default = [false]

function c_info(c::Clock)
    s1 = "state=$(c.state), "
    s2 = "t=$(round(c.time, sigdigits=4)) $(c.unit), "
    s3 = "Δt=$(round(c.Δt, sigdigits=4)) $(c.unit)"
    return s1*s2*s3
end

function sc_info(c::Clock)
    sc1 = "ev:$(length(c.sc.events)), "
    sc2 = "cev:$(length(c.sc.cevents)), "
    sc3 = "sampl:$(length(c.sc.samples))"
    return sc1*sc2*sc3
end

function pretty_print(c::Clock)
    if c.id > 0
        return pretty_print(pclock(c, c.id))
    else
        s1 = "Clock 0, thrd 1 (+ $(length(c.ac)) ac): "
        s2 = c_info(c)
        s3 = ", prc:$(length(c.processes))"
        return s1*s2*s3*"\n  scheduled "*sc_info(c)*"\n"
    end
end

function pretty_print(ac::ActiveClock)
    s1 = "Active clock $(ac.id) on thrd $(ac.thread): "
    s2 = c_info(ac.clock)
    s3 = ", prc:$(length(ac.clock.processes))"
    return s1*s2*s3*"\n   scheduled "*sc_info(ac.clock)*"\n"
end

function Base.show(io::IO, c::AbstractClock)
    if _show_default[end]
        return Base.show_default(IOContext(io, :compact => false), c)
    else
        return print(io, pretty_print(c))
    end
end


text:
# ---------------------------------------------------------------------------
# This is experimental, and currently unused
# ---------------------------------------------------------------------------

export
    Operator,
    OpGraph,
    InputNode,
    OutputNode,
    Learnable,
    Constant,
    AbstractTransformation,
    op,
    add_op,
    UNARY,
    BINARY,
    ELEMENTWISE,
    AGGREGATOR,
    @op

"""
The idea here is that an Operator wraps an arbitrary function, and there are static query methods
which define how to apply the function to input(s) and how to update the output.

In addition, we should be able to define (as part of the type, and thus at compile time) inversion
functions, deriviatives, etc.

some core traits:
- is_mappable: one input, function applied elementwise, input size matches output size
    examples: sin, exp, log, sigmoid, tanh
- is_operator: two (or more?) inputs, output dimension/size is function of inputs
    examples: * (matrix multiply)
- is_element_operator: two (or more?) inputs, applied elementwise, input sizes match output size
    examples: +, -, .*, ./, max
- is_aggregator: one input, output dimension is less than input dimension
    examples: sum, product
"""

@enum OpType UNARY BINARY ELEMENTWISE AGGREGATOR

const _op_types = Dict{Symbol, Tuple{Function,OpType}}(
    :sigmoid => (sigmoid, UNARY),
    :tanh    => (tanh, UNARY),
    :*       => (*, BINARY),
    :+       => (+, ELEMENTWISE),
    :minimum => (minimum, AGGREGATOR),
)

op_function(s::Symbol) = _op_types[s][1]
op_type(s::Symbol)     = _op_types[s][2]

function add_op(s::Symbol, f::Function, t::OpType)
    _op_types[s] = (f, t)
end

is_op(s::Symbol) = haskey(_op_types, s)

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

# wrap a function call (static transformation) and input/output dimensions in a type signature
immutable Operator{OP,F,I,O} <: AbstractTransformation{I,O}
    f::Function
end

# convenience to instantiate an Operator
op(f::Function, I::Int, O::Int) = Operator{op_type(Symbol(f)), Symbol(f), I, O}(f)
op(s::Symbol, I::Int, O::Int) = Operator{op_type(s), s, I, O}(op_function(s))

(o::Operator)(args...) = (o.f)(args...)

LearnBase.transform!(y, o::Operator{UNARY}, args...) = map!(o.f, y, args...)
LearnBase.transform!(y, o::Operator{ELEMENTWISE}, args...) = map!(o.f, y, args...)
LearnBase.transform!(y, o::Operator{BINARY}, args...) = y[:] = o.f(args...)
LearnBase.transform!(y, o::Operator{AGGREGATOR}, args...) = y[:] = o.f(args...)

function Base.show{OP,F,I,O}(io::IO, o::Operator{OP,F,I,O})
    print(io, "Op{$OP, $F, $I, $O}")
end

is_unary{OPTYPE}(o::Operator{OPTYPE})       = OPTYPE == UNARY
is_binary{OPTYPE}(o::Operator{OPTYPE})      = OPTYPE == BINARY
is_elementwise{OPTYPE}(o::Operator{OPTYPE}) = OPTYPE == ELEMENTWISE
is_aggregator{OPTYPE}(o::Operator{OPTYPE})  = OPTYPE == AGGREGATOR

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

# const _mappables = [:sin, :exp, :log, :sigmoid, :tanh, :sqrt, :^]
# const _operators = [:*, :/, :\]
# const _element_ops = [:+, :-, :.*, :./, :max, :min, :.^]
# const _aggregators = [:sum, :product, :minimum, :maximum]

# function is_op(s::Symbol)
#     for v in (_mappables, _operators, _element_ops, _aggregators)
#         if s in v
#             return true
#         end
#     end
#     false
# end
# const _all_ops = vcat(_mappables, _operators, _element_ops, _aggregators)

# @generated function is_mappable{F,I,O}(o::Operator{F,I,O})
#     if F in _mappables
#         I==O || error("Mappables input and output dimensions should match")
#         :(true)
#     else
#         :(false)
#     end
# end
# @generated function is_operator{F,I,O}(o::Operator{F,I,O})
#     F in _operators ? :(true) : :(false)
# end
# @generated function is_element_operator{F,I,O}(o::Operator{F,I,O})
#     F in _element_ops ? :(true) : :(false)
# end
# @generated function is_aggregator{F,I,O}(o::Operator{F,I,O})
#     F in _aggregators ? :(true) : :(false)
# end

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

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

# input and output nodes allow for easy connectivity
immutable InputNode{I} <: AbstractTransformation{I,I} end
immutable OutputNode{O} <: AbstractTransformation{O,O} end

# a Learnable has no inputs, but produces an output (the parameters)
immutable Learnable{O} <: AbstractTransformation{0,O}
    name::Symbol
end

immutable Constant{O} <: AbstractTransformation{0,O}
    value::Float64
end

# an OpGraph is a lightweight representation of a directed graph of AbstractTransformations
type OpGraph{I,O} <: AbstractTransformation{I,O}
    nodes::Vector{AbstractTransformation}
    edges::Vector{NTuple{2,Int}}
    names::Vector{String}
end

# op_type(g::OpGraph) =

using PlotRecipes
@recipe function f(g::OpGraph)
    arrow --> arrow()
    markersize --> 50
    markeralpha --> 0.2
    linealpha --> 0.4
    linewidth --> 2
    names --> g.names
    func --> :tree
    root --> :left
    PlotRecipes.GraphPlot((Int[e[1] for e=g.edges], Int[e[2] for e=g.edges]))
end

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

"""
We want to build an op graph, where each node has input(s) and output(s)

```
# single input, single output
# w and b are assumed 'learnable parameters'
@op affine(x) = x * w + b

# single input, single output
# no learnable parameters
@op sigmoid(x) = 1 / (1 + exp(-x))

# this shouldn't need a special definiton... use the built-in method
# no learnable parameters
@op tanh(x) = tanh(x)

# 2 inputs, 2 outputs
w is a learnable parameter, and can be fused to avoid recomputing
@op f(x, y) = (w * (x-y), w * (x+y))
```

The point of this design is that we should be able to define simple functions
composed of previously defined components, and the `@op` macro will build
a connected graph of transformations.  Some notes:

- If all components and operations inside an op graph are differentiable, then so is the graph
    - macro should build `deriv` method (and others in a similar way)
- If a variable appears in the body, but not in the function signature, we assume it
    is a `LearnableParameter`, and include it (with gradient state objects) in the internals
    of the op graph (or its sub-components).
"""

const _input_index = 1
const _output_index = 2

# NOTE: in the add_item_to_graph methods, we want to:
#   - increase the numnodes
#   - call add_edge! for each connection

# remember that numnodes is also the index of the last-added node

function add_item_to_graph!(block::Expr, input_idx::Int, variables, numnodes::Base.RefValue{Int}, item, parent_idx)
    @show "EXPR",item input_idx variables numnodes[]
    dump(item, maxdepth=20)

    if isa(item, Symbol)
        _add_symbol_to_graph!(block, numnodes, item, parent_idx, variables)
    elseif isa(item, Number)
        _add_const_to_graph!(block, numnodes, item, parent_idx)

    elseif item.head == :call && is_op(item.args[1])
        # node_indices = Int[]
        func = shift!(item.args)

        add_node!(block, numnodes, func, :operation)
        this_idx = numnodes[]
        for x in item.args
            if isa(x, Symbol)
                _add_symbol_to_graph!(block, numnodes, x, this_idx, variables)
                # if x in variables
                #     # TODO: we should have multiple input nodes, one for each variable
                #     # connect the input node to this op
                #     add_edge!(block, _input_index, this_idx)
                # else
                #     # TODO: this is a Learnable... add a node and connect it to this op
                #     add_node!(block, numnodes, x, :learnable)
                #     add_edge!(block, numnodes[], this_idx)
                # end
            elseif isa(x, Number)
                _add_const_to_graph!(block, numnodes, x, this_idx)

            else
                # this is another op... build it
                add_item_to_graph!(block, input_idx, variables, numnodes, x, this_idx)
            end
        end

        # now connect this op to its parent
        add_edge!(block, this_idx, parent_idx)

    elseif item.head == :line
        # if it's a line number block, pass through to wrapped arg
        # add_item_to_graph!(block, input_idx, variables, numnodes, item.args[1])
        info("line")

    else
        error("OpGraph parse error: $item")
    end
    return
end

function _add_symbol_to_graph!(block, numnodes, x, parent_idx, variables)
    if x in variables
        # TODO: we should have multiple input nodes, one for each variable
        # connect the input node to this op
        add_edge!(block, _input_index, parent_idx)
    else
        # TODO: this is a Learnable... add a node and connect it to this op
        add_node!(block, numnodes, x, :learnable)
        add_edge!(block, numnodes[], parent_idx)
    end
end

function _add_const_to_graph!(block, numnodes, x, parent_idx)
    add_node!(block, numnodes, x, :constant)
    add_edge!(block, numnodes[], parent_idx)
end

# function add_item_to_graph!(block::Expr, input_idx::Int, variables, numnodes::Base.RefValue{Int}, item::Symbol, opidx)
#     @show "SYM",item input_idx variables numnodes[]
#     if item in variables
#         # TODO: we should have multiple input nodes, one for each variable
#         add_edge!(block, _input_index, opidx)
#     else
#         # TODO: this is a Learnable
#         add_node!(block, numnodes, item, :learnable)
#         add_edge!(block, numnodes[], opidx)
#     end
#     return opidx
# end

function add_node!(block::Expr, numnodes, node, nodetype::Symbol)
    nodesym = quot(node)
    nodeexpr = if nodetype == :learnable
        :(Learnable{1}($nodesym))
    elseif nodetype == :constant
        :(Constant{1}($node))
    elseif nodetype == :operation
        :(op($nodesym, 1, 1))
    else
        error("how do I add this node? $node $nodetype")
    end
    push!(block.args, esc(:(push!(g.nodes, $nodeexpr))))
    push!(block.args, esc(:(push!(g.names, string($nodesym)))))
    numnodes[] += 1
end

function add_edge!(block::Expr, i::Int, j::Int)
    push!(block.args, esc(:(push!(g.edges, ($i, $j)))))
end

function _op_macro(funcexpr::Expr, inout::NTuple{2,Int} = (1,1))
    dump(funcexpr; maxdepth=20)
    @show inout typeof(inout)

    func_signature, func_body = funcexpr.args

    if !(funcexpr.head in (:(=), :function))
        error("Must wrap a valid function call!")
    end
    if !(isa(func_signature, Expr) && func_signature.head == :call)
        error("Expected `func_signature = ...` with func_signature as a call Expr... got: $func_signature")
    end

    func_name = func_signature.args[1]
    variables = func_signature.args[2:end]
    @show func_name variables

    # TODO: we want to build an OpGraph as a function of Variable, Learnable, and Op components
    # To simplify, we'll assume that:
    #   - Expr --> Op
    #   - Symbol --> haskey(_ops, k) ? Variable : Learnable

    # Variables are not learnable, and they are NOT part of the graph.  They represent the inputs
    # to the OpGraph.

    I,O = inout

    block = Expr(:block)
    push!(block.args, esc(:(
        g = OpGraph{$I, $O}(
            AbstractTransformation[InputNode{$I}(), OutputNode{$O}()],
            NTuple{2,Int}[],
            String["Input", "Output"]
        )
    )))
    @show block

    # graph_nodes = [InputNode(), OutputNode()]
    # graph_edges = []
    numnodes = Ref(2) # input and output
    for item in func_body.args
        opidx = add_item_to_graph!(block, _input_index, variables, numnodes, item, 2)

        # # connect this node to the output node, since it is returned from the function
        # add_edge!(block, opidx, _output_index)
    end

    # # add this to the operators list
    # push!(_operators, func_name)
    push!(block.args, esc(:(Transformations.add_op($(quot(func_name)), g, Transformations.op_type(g)))))

    push!(block.args, esc(:(g)))
    @show block
    block
end

macro op(args...)
    _op_macro(args...)
end

#####
# below this point is old code:

# Base.call(o::Operator, input) = map(o.f, input)
# transform(o::Operator, input) = map(o.f, input)

# # operator(f::Function, I, O) = Operator{Symbol(f), I, O}(f)

# # same dimensions
# function transform!{T,F,I,O}(output::AbstractArray{T,O}, o::Operator{F,I,O}, input::AbstractArray{T,I})
#     map!(o.f, output, input)
# end
# function transform!{T,F,I,O}(output::AbstractArray{T,O}, o::Operator{F,I,O}, input::AbstractArray{T,I})
#     map!(o.f, output, input)
# end

# # @generated function transform!{T,F,I,O}(output::AbstractArray{T,O}, o::Operator{F,I,O}, input::AbstractArray{T,I})
# #     @show T, F, O, I

# #     :(for (i,j) in zip(eachindex(output), eachindex(input))
# #         output[i] = (o.f)(input[j])
# #     end)
# # end


text:
benchmarkpkg("TextStylometry",
             BenchmarkConfig(id = "master",
                             juliacmd = `/usr/bin/julia -O3`))
benchmarkpkg("TextStylometry",
             BenchmarkConfig(id = "master",
                             juliacmd = `/home/bor/Dependencies/julia/julia -O3`))


text:
using ConstraintSolver

CS = ConstraintSolver
include("../../test/sudoku_fcts.jl")

function from_file(filename, sep='\n')
    s = open(filename) do file
        read(file, String)
    end
    str_sudokus = split(strip(s), sep)
    grids = AbstractArray[]
    for str_sudoku in str_sudokus
        str_sudoku = replace(str_sudoku, "."=>"0")
        one_line_grid = parse.(Int, split(str_sudoku,""))
        grid = reshape(one_line_grid, 9, 9)
        push!(grids, grid)
    end
    return grids
end

function solve_all(grids; benchmark=false, single_times=true)
    ct = time()
    grids = grids
    for (i,grid) in enumerate(grids)
        com = CS.init()

        com_grid = create_sudoku_grid!(com, grid)
        add_sudoku_constr!(com, com_grid)

        if single_times
            GC.enable(false)
            t = time()
            status = solve!(com);
            t = time()-t
            GC.enable(true)
            println(i-1,", ", t)
        else
            GC.enable(false)
            status = solve!(com);
            GC.enable(true)
        end
        if !benchmark
            @show com.info
            println("Status: ", status)
            @assert fulfills_sudoku_constr(com_grid)
        end
    end
    println("")
    tt = time()-ct
    println("total time: ", tt)
    println("avg: ", tt/length(grids))
end

function solve_one(grid)
    com = CS.init()

    CS.build_search_space!(com, grid,[1,2,3,4,5,6,7,8,9],0)
    add_sudoku_constr!(com, grid)

    status = solve!(com; backtrack=false);
    return com
end

function main(; benchmark=false, single_times=true)
    solve_all(from_file("benchmark/sudoku/data/top95.txt"); benchmark=benchmark, single_times=single_times)
    # solve_all(from_file("hardest.txt"), "hardest")
end


text:
# @inline function if_else(condition::Bool, @nospecialize(trueval), @nospecialize(falseval))
#            return ifelse(condition, trueval, falseval)
# end

# ModelingToolkit.@register if_else(x, trueval, falseval)
# # This is needed to bring if_else into MTK's scope.
# @eval ModelingToolkit using ..Main: if_else

# nasty(i) = if_else(i > 1, 1, 0)
# @variables t
# nasty(t)

function MassSpringOscillator(input)

  @parameters t
  @parameters k,c,m
  D = Differential(t)
  @variables pos(t) vel(t)

  eqs = [D(pos) ~ vel,
        D(vel) ~ (-1/m)*(c*vel + k*pos - input(t))
  ]

  ps = [k,c,m] .=>  [2.,1.,4.]
  ics = [pos, vel] .=> [1.,0.]
  od = ODESystem(eqs, t, first.(ics),first.(ps))
  tspan = (0.,100.)
  # prob = ODEProblem(od, ics, tspan, ps)
  return od, ics, tspan, ps
end

# od,ic,ps = create(t->0)
# tspan = (0.,10.)
# prob = ODEProblem(od,ic,tspan,ps)