text:
function Documented(def::StructDefinition{false})
    (; spec) = def
    doc = string(
        "Intermediate wrapper for $(spec.name).",
        extension_doc(spec),
        api_doc(spec),
    )
    Documented(def, doc)
end

function Documented(def::StructDefinition{true})
    (; spec) = def
    doc = string(
        "High-level wrapper for $(spec.name).",
        extension_doc(spec),
        api_doc(spec),
    )
    Documented(def, doc)
end

Documented(def::Constructor{<:StructDefinition{HL, SpecStruct}}) where {HL} = document_constructor(def, def.to.spec, HL)
Documented(def::Constructor{HandleDefinition,APIFunction{APIFunction{CreateFunc}}}) = document_constructor(def, def.from.spec.spec, true)
Documented(def::Constructor{HandleDefinition,APIFunction{CreateFunc}}) = document_constructor(def, def.from.spec, false)

function document_constructor(def::Constructor, spec, is_hl)
    (; p) = def
    doc = string(
        extension_doc(spec),
        args_summary(document_arguments(p, spec, is_hl)),
        api_doc(spec),
    )
    Documented(def, doc)
end

backquoted(arg) = string('`', arg, '`')

function document_return_codes(spec::SpecFunc)
    res = ""
    if !isempty(spec.success_codes) || !isempty(spec.error_codes)
        res *= "Return codes:"
    end
    if !isempty(spec.success_codes) && must_return_success_code(spec)
        res *= string("\n- ", join(backquoted.(remove_vk_prefix.(spec.success_codes)), "\n- "))
    end
    if !isempty(spec.error_codes)
        res *= string("\n- ", join(backquoted.(remove_vk_prefix.(spec.error_codes)), "\n- "))
    end

    isempty(res) && return ""
    '\n'^2 * res
end
document_return_codes(spec::CreateFunc) = document_return_codes(spec.func)
document_return_codes(spec) = ""

function append_to_argdoc!(argdocs, spec, str)
    identifier = wrap_identifier(spec)
    doc = argdocs[identifier]
    argdocs[identifier] = doc * str
end

api_doc(spec::Spec) = string("\n\n", "[API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/$(spec.name).html)")
api_doc(spec::CreateFunc) = api_doc(spec.func)
api_doc(::Nothing) = ""

extension_doc(::Nothing) = ""
function extension_doc(spec)
    ext = extension(spec)
    isnothing(ext) && return ""
    string("\n\n", "Extension: ", replace(string(ext.name), '_' => "\\\\_"))
end

function args_summary(argdocs)
    isempty(argdocs) && return ""
    string("\n\n", "Arguments:\n- ", join(argdocs, "\n- "))
end

Documented(def::APIFunction{SpecFunc}) = Documented(def, def.spec, def.p, false)
Documented(def::APIFunction{APIFunction{SpecFunc}}) = Documented(def, def.spec.spec, def.p, true)
function Documented(def::APIFunction{T}) where {T<:Union{CreateFunc,APIFunction{CreateFunc}}}
    is_hl = T === APIFunction{CreateFunc}
    create_func = is_hl ? def.spec.spec : def.spec
    spec = create_func.batch ? nothing : create_func
    Documented(def, spec, def.p, is_hl)
end

function extra_doc(spec::SpecFunc)
    params = parameters(spec)
    if any(is_data_with_retrievable_size, params)
        """

        !!! warning
            The pointer returned by this function holds memory owned by Julia. It is therefore **your** responsibility to free it after use (e.g. with `Libc.free`)."""
    else
        ""
    end
end
extra_doc(spec::CreateFunc) = extra_doc(spec.func)
extra_doc(spec) = ""

function update_argdocs!(argdocs, spec::SpecFunc)
    params = parameters(spec)
    if any(is_opaque_data, params)
        opaque_params = params[findall(is_opaque_data, params)]
        for param in opaque_params
            append_to_argdoc!(
                argdocs,
                param,
                " (must be a valid pointer with `$(wrap_identifier(len(param)))` bytes)",
            )
        end
    end
end
update_argdocs!(argdocs, spec::CreateFunc) = update_argdocs!(argdocs, spec.func)
update_argdocs!(argdocs, spec) = nothing

function Documented(def::APIFunction, spec, p, is_hl)
    argdocs = document_arguments(p, spec, is_hl)
    update_argdocs!(argdocs, spec)

    doc = string(
        extension_doc(spec),
        document_return_codes(spec),
        args_summary(argdocs),
        extra_doc(spec),
        api_doc(spec),
    )
    Documented(def, doc)
end

parameters(spec::Spec) = children(spec)
parameters(spec::CreateFunc) = [collect(children(spec.func)); collect(children(spec.create_info_struct))]

function document_arguments(p, spec::Spec, is_hl)
    params = parameters(spec)
    args, kwargs = get(Vector, p, :args), get(Vector, p, :kwargs)
    d = Dictionary{Symbol,String}()
    for arg in args
        identifier = name(arg)
        str = @match identifier begin
            :fptr => string(backquoted(arg), ": function pointer used for the API call")
            :fptr_create => string(backquoted(arg), ": function pointer used for creating the handle(s)")
            :fptr_destroy => string(backquoted(arg), ": function pointer used for destroying the handle(s) upon finalization")
            :next_types => string(backquoted(arg), ": types of members to initialize and include as part of the `next` chain")
            _ => nothing
        end
        if !isnothing(str)
            insert!(d, identifier, str)
            continue
        end
        i = findfirst(==(identifier) ∘ wrap_identifier, params)
        str = if !isnothing(i)
            param = params[i]
            t = is_hl ? hl_type(param) : idiomatic_julia_type(param)
            str = backquoted(:($identifier::$t))
            param.is_externsync ? str * " (externsync)" : str
        else
            backquoted(arg)
        end
        insert!(d, identifier, str)
    end
    for kwarg in kwargs
        identifier, value = kwarg.args
        i = findfirst(==(identifier) ∘ wrap_identifier, params)
        str = if !isnothing(i)
            param = params[i]
            t = is_hl ? hl_type(param) : idiomatic_julia_type(param)
            str = string(backquoted(:($identifier::$t)), ": defaults to ", backquoted(value))
            param.is_externsync ? str * " (externsync)" : str
        else
            string(backquoted(identifier), ": defaults to ", backquoted(value))
        end
        insert!(d, identifier, str)
    end
    d
end

function docstring(ex, docstring)
    endswith(docstring, '\n') || (docstring *= '\n')
    Dict(
        :category => :doc,
        :ex => ex,
        :docstring => docstring,
    )
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text: abstract type AbstractPotential end

ext: jl

max_stars_repo_name: heindelj/MolecularPotentials

max_issues_repo_name: heindelj/MolecularPotentials

max_forks_repo_name: heindelj/MolecularPotentials


text:
#
# Computes the total potential energy of the system using r2
#

function utotal2(p,data)
    first_atom = zeros(Int64, data.nc, data.nc)
    next_atom = zeros(Int64, data.N)
    makell!(first_atom, next_atom, p, data)
    ut = 0.
    for ip = 1:data.N
        icell = trunc(Int64, p[ip][1]/data.cutoff) + 1
        jcell = trunc(Int64, p[ip][2]/data.cutoff) + 1
        for ic = icell-1:icell+1
            iw = wrap_cell(data.nc, ic)
            for jc = jcell-1:jcell+1
                jw = wrap_cell(data.nc, jc)
                jp = first_atom[iw,jw]
                while jp > 0
                    if ip < jp
                        r2 = r2pbc(p[ip], p[jp], data.bdim)
                        if r2 ≤ data.cutoff2
                            ut += upair2(r2, data)
                        end
                    end
                    jp = next_atom[jp]
                end
            end
        end
    end
    return ut
end


text:
using OnlineStats

export
    LinearGaussianBayesianNetBehavior,
    LB_TrainParams,
    LB_PreallocatedData

"""
The ordering of the parental instantiations in discrete networks follows the convention
defined in Decision Making Under Uncertainty.

Suppose a variable has three discrete parents. The first parental instantiation
assigns all parents to their first bin. The second will assign the first
parent (as defined in `parents`) to its second bin and the other parents
to their first bin. The sequence continues until all parents are instantiated
to their last bins.

This is a directly copy from Base.sub2ind but allows for passing a vector instead of separate items

Note that this does NOT check bounds
"""
function sub2ind_vec{T<:Integer, S<:Integer}(dims::AbstractVector{T}, I::AbstractVector{S})
    N = length(dims)
    if N == 0
        return 1
    end
    @assert(N == length(I))

    ex = I[N] - 1
    for i in N-1:-1:1
        if i > N
            ex = (I[i] - 1 + ex)
        else
            ex = (I[i] - 1 + dims[i]*ex)
        end
    end

    ex + 1
end

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

immutable LinearGaussianStats
    # μ = wᵀx where x[end] = 1
    w::Vector{Float64} # [nparents + 1]
    σ::Float64
end
n_sufficient_statistics(s::LinearGaussianStats) = length(s.w) + 1

type ConditionalLinearGaussianNode
    #=
    Defines a variable in a Conditional Linear Gaussian Bayesian Network

    p(x|pa(x)) ~ Normal(wᵀcont + b) ∀ disc

    There is one linear gaussian for each combination of discrete parents
    =#

    index::Int # index of this node in the assignment list
    stats::Matrix{LinearGaussianStats} # [n_nann_parent_instantiations × n_disc_parent_instantiations]
    parents_disc::Vector{Int} # index, in assignment_disc, of discrete parents
    parents_cont::Vector{Int} # index, in assignment_cont, of continuous parents, may include other nodes's vars
    parents_nann::Vector{Int} # index, in assignment_cont, of nannable parents

    parental_assignments_disc::Vector{Int} # discrete assignment, preallocated array for sub2ind calculation
    parental_assignments_nann::Vector{Int} # discrete assignment, preallocated array for sub2ind calculation
    parental_assignments_cont::Vector{Float64} # continuous assignment (last entry is 1.0) [cont | nann]
    parent_instantiation_counts_disc::Vector{Int} # number of possible instantiations for each discrete parent
    parent_instantiation_counts_nann::Vector{Int} # number of possible instantiations for each nannable parent (all 2's)
end

function n_sufficient_statistics(node::ConditionalLinearGaussianNode)
    retval = 0
    for s in node.stats
        retval += n_sufficient_statistics(s)
    end
    retval
end

function _get_normal(
    node::ConditionalLinearGaussianNode,
    assignment_disc::Vector{Int},    # note: full array into which we must index with node.parents_disc
    assignment_cont::Vector{Float64} # note: full array into which we must index with node.parents_cont
    )

    # println("assignment_disc: ", assignment_disc)
    # println("assignment_cont: ", assignment_cont)

    for (i,v) in enumerate(node.parents_nann)
        node.parental_assignments_nann[i] = isnan(assignment_cont[v]) ? 1 : 2
    end
    q = sub2ind_vec(node.parent_instantiation_counts_nann, node.parental_assignments_nann)

    # if q == 0
    #     println("PROBLEM!")
    #     println("index: ", node.index)
    #     println("parental_assignments_nann: ", node.parental_assignments_nann)
    #     println("parent_instantiation_counts_nann: ", node.parent_instantiation_counts_nann)
    #     sleep(0.25)
    # end

    for (i,v) in enumerate(node.parents_disc)
        node.parental_assignments_disc[i] = assignment_disc[v]
    end
    j = sub2ind_vec(node.parent_instantiation_counts_disc, node.parental_assignments_disc)

    # if j == 0
    #     println("PROBLEM!")
    #     println("index: ", node.index)
    #     println("parental_assignments_disc: ", node.parental_assignments_disc)
    #     println("parent_instantiation_counts_disc: ", node.parent_instantiation_counts_disc)
    #     sleep(0.25)
    # end

    # println("q: ", q)
    # println("j: ", j)

    stat = node.stats[q,j]

    # println("stat: ", stat)

    μ = stat.w[end]

    i = 0
    for p in node.parents_cont
        i += 1
        μ += stat.w[i] * assignment_cont[p]
    end
    for p in node.parents_nann
        if !isnan(assignment_cont[p])
            i += 1
            μ += stat.w[i] * assignment_cont[p]
        end
    end
    @assert(i == length(stat.w)-1)

    normal = Normal(μ, stat.σ)
end

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

type LinearGaussianBayesianNetBehavior <: AbstractVehicleBehavior

    targets::ModelTargets
    extractor_disc :: FeatureSubsetExtractor
    extractor_cont :: FeatureSubsetExtractor
    clamper_cont ::DataClamper
    clamper_act :: DataClamper

    sample_lat_first::Bool # used to ensure topological ordering is preserved
    node_lat::ConditionalLinearGaussianNode
    node_lon::ConditionalLinearGaussianNode
    assignment_disc::Vector{Int} # [d]
    assignment_cont::Vector{Float64} # [c + 2], preallocated memory, with cont | lat | lon

    function LinearGaussianBayesianNetBehavior(
        targets::ModelTargets,
        extractor_disc::FeatureSubsetExtractor,
        extractor_cont::FeatureSubsetExtractor,
        clamper_cont::DataClamper,
        clamper_act::DataClamper,
        sample_lat_first::Bool,
        node_lat::ConditionalLinearGaussianNode,
        node_lon::ConditionalLinearGaussianNode,
        )

        retval = new()

        retval.targets = deepcopy(targets)
        retval.extractor_disc = extractor_disc
        retval.extractor_cont = extractor_cont
        retval.clamper_cont = clamper_cont
        retval.clamper_act = clamper_act
        retval.sample_lat_first = sample_lat_first
        retval.node_lat = node_lat
        retval.node_lon = node_lon
        retval.assignment_disc = Array(Int, length(retval.extractor_disc.indicators))
        retval.assignment_cont = Array(Float64, length(retval.extractor_cont.indicators)+2)

        retval
    end
end
function Base.print(io::IO, LB::LinearGaussianBayesianNetBehavior)

    symbols_disc = convert(Vector{Symbol}, map(f->symbol(f), LB.extractor_disc.indicators))
    symbols_cont = convert(Vector{Symbol},[map(f->symbol(f), LB.extractor_cont.indicators);
                                           symbol(LB.targets.lat); symbol(LB.targets.lon)])

    if !isempty(LB.extractor_cont.indicators)
        num_cont = sum(f->!couldna(f), LB.extractor_cont.indicators)
    else
        num_cont = 0
    end

    println(io, "Linear Gaussian Bayesian Network Behavior")
    println(io, "\ttargets: ", LB.targets)
    println(io, "\tfeatures: ")
    println(io, "\t\tdiscrete:   ", symbols_disc)
    println(io, "\t\tcontinuous: ", symbols_cont)
    println(io, "\tnumber of suffstats for lat: ", n_sufficient_statistics(LB.node_lat))
    println(io, "\tnumber of suffstats for lon: ", n_sufficient_statistics(LB.node_lon))
    println(io, "\tparents lat: ", symbols_disc[LB.node_lat.parents_disc])
    println(io, "\t             ", symbols_cont[LB.node_lat.parents_cont])
    println(io, "\t             ", symbols_cont[LB.node_lat.parents_nann])
    println(io, "\tparents lon: ", symbols_disc[LB.node_lon.parents_disc])
    println(io, "\t             ", symbols_cont[LB.node_lon.parents_cont])
    println(io, "\t             ", symbols_cont[LB.node_lon.parents_nann])
    println(io, "\tsample lat first: ", LB.sample_lat_first)
end
AutomotiveDrivingModels.get_targets(behavior::LinearGaussianBayesianNetBehavior) = behavior.targets
function AutomotiveDrivingModels.get_indicators(behavior::LinearGaussianBayesianNetBehavior)
    retval = unique([behavior.extractor_disc.indicators; behavior.extractor_cont.indicators])
end

type LB_TrainParams <: AbstractVehicleBehaviorTrainParams

    targets::ModelTargets
    indicators::Vector{AbstractFeature} # list of all potential indicators

    ridge_regression_constant::Float64
    min_σ_lat::Float64 # minimum standard deviation for lateral target
    min_σ_lon::Float64 # minimum standard deviation for longitudinal target

    max_parents::Int # maximum number of parents per node
    verbosity::Int

    function LB_TrainParams(;
        targets::ModelTargets = ModelTargets(Features.FUTUREDESIREDANGLE, Features.FUTUREACCELERATION),
        indicators::Vector{AbstractFeature} = AbstractFeature[],
        ridge_regression_constant::Float64=0.001,
        min_σ_lat::Float64=1e-4,
        min_σ_lon::Float64=1e-5,
        max_parents::Int=4,
        verbosity::Int=0,
        )

        @assert(!in(targets.lat, indicators))
        @assert(!in(targets.lon, indicators))

        retval = new()

        retval.ridge_regression_constant = ridge_regression_constant
        retval.targets = targets
        retval.indicators = indicators
        retval.min_σ_lat = min_σ_lat
        retval.min_σ_lon = min_σ_lon

        retval.max_parents = max_parents
        retval.verbosity = verbosity

        retval
    end
end
function Base.print(io::IO, p::LB_TrainParams)
    println(io, "LB Train Params")
    println(io, "\ttargets: ", p.targets)
    println(io, "\tindicators: ", map(f->symbol(f), p.indicators))
    println(io, "\tridge_regression_constant: ", p.ridge_regression_constant)
    println(io, "\tmin_σ_lat:      ", p.min_σ_lat)
    println(io, "\tmin_σ_lon:      ", p.min_σ_lon)
    println(io, "\tmax_parents:    ", p.max_parents)
end

type LB_PreallocatedData <: AbstractVehicleBehaviorPreallocatedData

    Y::Matrix{Float64} # [2 × m] # two targets (lat, lon)
    X_disc::Matrix{Int} # [d × m]
    X_cont::Matrix{Float64} # [c × m]
    X_nann::Matrix{Float64} # [n × m]

    features_disc  :: Vector{AbstractFeature}
    features_cont  :: Vector{AbstractFeature}
    features_nann  :: Vector{AbstractFeature}
    clamper_cont   :: DataClamper
    clamper_nann   :: DataClamper
    clamper_act    :: DataClamper

    function LB_PreallocatedData(dset::ModelTrainingData2, params::LB_TrainParams)

        retval = new()

        targets = params.targets
        indicators = params.indicators
        trainingframes = dset.dataframe
        nframes = nrow(trainingframes)
        nindicators = length(indicators)

        X = Array(Float64, nindicators, nframes)
        Y = Array(Float64, 2, nframes)
        pull_design_and_target_matrices!(X, Y, trainingframes, targets, indicators)

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

        features_disc_indeces = find(f->isint(f), indicators)
        features_cont_indeces = find(f->!isint(f) && !couldna(f), indicators)
        features_nann_indeces = find(f->!isint(f) &&  couldna(f), indicators)

        retval.Y = Y
        retval.X_disc = round(Int, X[features_disc_indeces, :]) .+ 1
        retval.X_cont = X[features_cont_indeces, :]
        retval.X_nann = X[features_nann_indeces, :]
        retval.features_disc = indicators[features_disc_indeces]
        retval.features_cont = indicators[features_cont_indeces]
        retval.features_nann = indicators[features_nann_indeces]

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

        retval.clamper_cont = DataClamper(retval.X_cont)
        retval.clamper_nann = DataClamper(retval.X_nann)
        retval.clamper_act = DataClamper(retval.Y)

        retval
    end
end
function preallocate_learning_data(
    dset::ModelTrainingData2,
    params::LB_TrainParams)

    LB_PreallocatedData(dset, params)
end

function _cast_discrete_to_int!(behavior::LinearGaussianBayesianNetBehavior)
    for (i,v) in enumerate(behavior.extractor_disc.x)
        behavior.assignment_disc[i] = round(Int, v)
    end
    behavior
end
function _copy_to_assignment!(behavior::LinearGaussianBayesianNetBehavior)

    copy!(behavior.assignment_cont, behavior.clamper_cont.x)

    # # temporarily set NaN values to 0.0
    # for i in 1:length(behavior.assignment_cont)
    #     if isnan(behavior.assignment_cont[i])
    #         behavior.assignment_cont[i] = 0.0
    #     end
    # end

    behavior
end
function _sample_from_node!(behavior::LinearGaussianBayesianNetBehavior, node::ConditionalLinearGaussianNode)
    normal = _get_normal(node,
                         behavior.assignment_disc,
                         behavior.assignment_cont)
    action_lat = rand(normal)
    behavior.assignment_cont[node.index] = action_lat

    action_lat
end
function _process_obs(behavior::LinearGaussianBayesianNetBehavior)
    process!(behavior.clamper_cont) # TODO - ensure that clamper_cont is tied to extractor_cont
    _cast_discrete_to_int!(behavior)
    _copy_to_assignment!(behavior)
    behavior
end
function _observe_on_runlog(
    behavior::LinearGaussianBayesianNetBehavior,
    runlog::RunLog,
    sn::StreetNetwork,
    colset::UInt,
    frame::Int
    )

    observe!(behavior.extractor_cont, runlog, sn, colset, frame)
    observe!(behavior.extractor_disc, runlog, sn, colset, frame)
    behavior.extractor_disc.x += 1 # compensate for fact that observations start at 0

    _process_obs(behavior)
end
function _observe_on_dataframe(
    behavior::LinearGaussianBayesianNetBehavior,
    dataframe::DataFrame,
    frameind::Integer,
    )

    observe!(behavior.extractor_cont, dataframe, frameind)
    observe!(behavior.extractor_disc, dataframe, frameind)
    behavior.extractor_disc.x += 1 # compensate for fact that observations start at 0

    _process_obs(behavior)
end
function _set_and_process_action!(behavior::LinearGaussianBayesianNetBehavior, action_lat::Float64, action_lon::Float64)
    behavior.clamper_act.x[1] = action_lat
    behavior.clamper_act.x[2] = action_lon
    process!(behavior.clamper_act)
    behavior.assignment_cont[behavior.node_lat.index] = behavior.clamper_act.x[1]
    behavior.assignment_cont[behavior.node_lon.index] = behavior.clamper_act.x[2]
    behavior
end

function select_action(
    behavior::LinearGaussianBayesianNetBehavior,
    runlog::RunLog,
    sn::StreetNetwork,
    colset::UInt,
    frame::Int
    )

    _observe_on_runlog(behavior, runlog, sn, colset, frame)

    if behavior.sample_lat_first
        behavior.clamper_act.x[1] = _sample_from_node!(behavior, behavior.node_lat)
        behavior.clamper_act.x[2] = _sample_from_node!(behavior, behavior.node_lon)
    else
        behavior.clamper_act.x[2] = _sample_from_node!(behavior, behavior.node_lon)
        behavior.clamper_act.x[1] = _sample_from_node!(behavior, behavior.node_lat)
    end

    process!(behavior.clamper_act) # clamp
    action_lat = behavior.clamper_act.x[1]
    action_lon = behavior.clamper_act.x[2]

    (action_lat, action_lon)
end

function _calc_action_loglikelihood(behavior::LinearGaussianBayesianNetBehavior)

    # NOTE: observation and setting action must have already occured

    logl = 0.0
    normal = _get_normal(behavior.node_lat, behavior.assignment_disc, behavior.assignment_cont)
    logl += logpdf(normal, behavior.clamper_act.x[1])

    normal = _get_normal(behavior.node_lon, behavior.assignment_disc, behavior.assignment_cont)
    logl += logpdf(normal, behavior.clamper_act.x[2])

    logl
end
function calc_action_loglikelihood(
    behavior::LinearGaussianBayesianNetBehavior,
    runlog::RunLog,
    sn::StreetNetwork,
    colset::UInt,
    frame::Int,
    action_lat::Float64,
    action_lon::Float64,
    )

    _observe_on_runlog(behavior, runlog, sn, colset, frame)
    _set_and_process_action!(behavior, action_lat, action_lon)
    _calc_action_loglikelihood(behavior)
end
function calc_action_loglikelihood(
    behavior::LinearGaussianBayesianNetBehavior,
    features::DataFrame,
    frameind::Integer,
    )

    action_lat = features[frameind, symbol(behavior.targets.lat)]::Float64
    action_lon = features[frameind, symbol(behavior.targets.lon)]::Float64

    _observe_on_dataframe(behavior, features, frameind)
    _set_and_process_action!(behavior, action_lat, action_lon)
    _calc_action_loglikelihood(behavior)
end

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

immutable NodeInTraining
    index::Int # index of this node in Y
    target_as_parent::Int # if != 0, index of other target in Y
    parents_disc::Vector{Int} # index, in X_disc, of discrete parents (always in sorted order)
    parents_cont::Vector{Int} # index, in X_cont, of continuous parents (always in sorted order)
    parents_nann::Vector{Int} # index, in X_nann, of NaNable parents

    NodeInTraining(index::Int) = new(index, 0, Int[], Int[], Int[])
    NodeInTraining(i::Int, t_as_p::Int, pdisc::Vector{Int}, pcont::Vector{Int}, pnann::Vector{Int}) = new(i, t_as_p, pdisc, pcont, pnann)
end
function Base.copy(node::NodeInTraining)
    NodeInTraining(node.index, node.target_as_parent,
                   copy(node.parents_disc), copy(node.parents_cont),
                   copy(node.parents_nann))
end
function set_new_parent_disc(node::NodeInTraining, parent::Int)
    parents_disc = sort!(unique(push!(copy(node.parents_disc), parent)))
    NodeInTraining(node.index, node.target_as_parent, parents_disc, copy(node.parents_cont), copy(node.parents_nann))
end
function set_new_parent_cont(node::NodeInTraining, parent::Int)
    parents_cont = sort!(unique(push!(copy(node.parents_cont), parent)))
    NodeInTraining(node.index, node.target_as_parent, copy(node.parents_disc), parents_cont, copy(node.parents_nann))
end
function set_new_parent_nann(node::NodeInTraining, parent::Int)
    parents_nann = sort!(unique(push!(copy(node.parents_nann), parent)))
    NodeInTraining(node.index, node.target_as_parent, copy(node.parents_disc), copy(node.parents_cont), parents_nann)
end
function set_new_parent_targ(node::NodeInTraining, parent::Int)
    @assert(node.target_as_parent == 0)
    NodeInTraining(node.index, parent, copy(node.parents_disc), copy(node.parents_cont), copy(node.parents_nann))
end
function set_new_parent(node::NodeInTraining, parent::Int, sym::Symbol)
    if sym == :disc
        set_new_parent_disc(node, parent)
    elseif sym == :cont
        set_new_parent_cont(node, parent)
    elseif sym == :nann
        set_new_parent_nann(node, parent)
    elseif sym == :targ
        set_new_parent_targ(node, parent)
    else
        error("unknown parent set type $sym")
    end
end
function Base.hash(node::NodeInTraining, h::UInt=one(UInt))
    hash(node.index, hash(node.parents_disc, hash(node.parents_cont, hash(node.target_as_parent, hash(node.parents_nann, h)))))
end
function Base.(:(==))(A::NodeInTraining, B::NodeInTraining)
    A.index            == B.index &&
    A.target_as_parent == B.target_as_parent &&
    A.parents_disc     == B.parents_disc &&
    A.parents_cont     == B.parents_cont &&
    A.parents_nann     == B.parents_nann
end
function Base.show(io::IO, node::NodeInTraining)
    println("Node In Training:")
    println("\tindex:            ", node.index)
    println("\ttarget_as_parent: ", node.target_as_parent)
    println("\tparents_disc:     ", node.parents_disc)
    println("\tparents_cont:     ", node.parents_cont)
    println("\tparents_nann:     ", node.parents_nann)
end
nparents_total(node::NodeInTraining) = length(node.parents_disc) + length(node.parents_cont) + length(node.parents_nann)

type TrainingData

    Y :: Matrix{Float64}
    X_disc :: Matrix{Float64}
    X_cont :: Matrix{Float64}
    X_nann :: Matrix{Float64}

    disc_parent_instantiations::Vector{Int}
    score_cache::Dict{NodeInTraining, Float64} # node -> (component_score)

    function TrainingData(
        training_data::ModelTrainingData2,
        preallocated_data::LB_PreallocatedData,
        foldset::FoldSet,
        )

        retval = new()

        retval.Y = copy_matrix_fold(preallocated_data.Y, foldset)
        retval.X_disc = copy_matrix_fold(preallocated_data.X_disc, foldset)::Matrix{Int}
        retval.X_cont = copy_matrix_fold(preallocated_data.X_cont, foldset)::Matrix{Float64}
        retval.X_nann = copy_matrix_fold(preallocated_data.X_nann, foldset)::Matrix{Float64}

        retval.disc_parent_instantiations = convert(Vector{Int}, map(f->round(Int, upperbound(f))+1, preallocated_data.features_disc))
        retval.score_cache = Dict{NodeInTraining, Float64}()

        retval
    end
end

nsamples(data) = size(data.Y, 1)

function _init_parent_stuff_nann(node::NodeInTraining)
    nparents_nann = length(node.parents_nann) + (node.target_as_parent != 0)

    # 1 -> all being nan, 2 -> first being cont, rest all nan, 3 -> 2nd cont, 4->1st 2 cont, etc.
    # last one is the optional target as parent
    parent_instantiation_counts_nann = fill(2, nparents_nann)

    # returns 1 if nparents_nann=0, which is what we want
    n_nann_parent_instantiations = prod(parent_instantiation_counts_nann)

    # NaN -> 1, cont -> 2
    parental_assignments_nann = Array(Int, nparents_nann)

    (parent_instantiation_counts_nann, n_nann_parent_instantiations, parental_assignments_nann)
end
function _init_parent_stuff_disc(node::NodeInTraining, data::TrainingData)

    nparents_disc = length(node.parents_disc)
    parent_instantiation_counts_disc = Array(Int, nparents_disc)
    n_disc_parent_instantiations = 1
    for (i,p) in enumerate(node.parents_disc)
        parent_instantiation_counts_disc[i] = data.disc_parent_instantiations[p]
        n_disc_parent_instantiations *= data.disc_parent_instantiations[p]
    end

    parental_assignments_disc = Array(Int, nparents_disc)

    (parent_instantiation_counts_disc, n_disc_parent_instantiations, parental_assignments_disc)
end

function _get_nann_assignment_index!(
    i::Int,
    node::NodeInTraining, # the indeces of the discrete parents
    data::TrainingData,
    parental_assignments_nann::Vector{Int}, # preallocated vector to hold the parental assignments
    parent_instantiation_counts_nann::Vector{Int} # all twos
    )

    for (j,p) in enumerate(node.parents_nann)
        parental_assignments_nann[j] = isnan(data.X_nann[i,p]) ? 1 : 2 # if NAN, then 1, if continuous, then 2
    end
    if node.target_as_parent != 0
        parental_assignments_nann[end] = isnan(data.Y[i,node.target_as_parent]) ? 1 : 2
    end

    sub2ind_vec(parent_instantiation_counts_nann, parental_assignments_nann)
end
function _get_disc_assignment_index!(
    i::Int,
    node::NodeInTraining, # the indeces of the discrete parents
    data::TrainingData,
    parental_assignments_disc::Vector{Int}, # preallocated vector to hold the parental assignments
    parent_instantiation_counts_disc::Vector{Int} # size of each parent's domain
    )

    for (j,p) in enumerate(node.parents_disc)
        parental_assignments_disc[j] = data.X_disc[i,p]
    end
    sub2ind_vec(parent_instantiation_counts_disc, parental_assignments_disc)
end
function _pull_cont_data!(i::Int, x::Vector{Float64}, node::NodeInTraining, data::TrainingData)

    y = data.Y[i, node.index]

    j = 0
    for p in node.parents_cont
        j += 1
        x[j] = data.X_cont[i,p]
    end
    for p in node.parents_nann
        v = data.X_nann[i,p]
        if !isnan(v)
            j += 1
            x[j] = v
        end
    end
    if node.target_as_parent != 0
        v = data.Y[i,node.target_as_parent]
        if !isnan(v)
            j += 1
            x[j] = v
        end
    end

    if j != length(x)-1
        println("node")
        println(node)
        println("x: ", x)
        println("j: ", j)
        println("expected: ", length(x)-1)
        println("target as parent: ", node.target_as_parent)
        sleep(0.5)
    end

    @assert(j == length(x)-1)

    y
end
function _update_ridge_regression!(
    lhs::Vector{Float64},
    rhs::Matrix{Float64},
    x::Vector{Float64},
    y::Float64
    )

    lenx = length(x)

    for i in 1 : lenx
        lhs[i] += x[i]*y

        for j in 1 : lenx
            rhs[i,j] += x[i]*x[j]
        end
    end

    nothing
end
function _calc_logl_contribution(w::Vector{Float64}, σ::Float64, x::Vector{Float64}, y::Float64)
    μ = dot(w, x)
    logpdf(Normal(μ, σ), y)
end

function _perform_ridge_regression(
    node::NodeInTraining,
    data::TrainingData,
    params::LB_TrainParams,
    parent_instantiation_counts_nann::Vector{Int},
    n_nann_parent_instantiations::Int,
    parental_assignments_nann::Vector{Int},
    parent_instantiation_counts_disc::Vector{Int},
    n_disc_parent_instantiations::Int,
    parental_assignments_disc::Vector{Int},
    )

    # perform ridge-regression and stdev calc online
    #   Xᵀy = (λI + XᵀX)w

    # LHS is the Xᵀy vector [nparents+1]
    # RHS is the λI + XᵀX matrix [nparents+1 × nparents+1]

    λ = params.ridge_regression_constant
    σ_min = node.index == 1 ? params.min_σ_lat : params.min_σ_lon

    nparents_disc = length(node.parents_disc)
    nparents_cont = length(node.parents_cont)
    nparents_nann = length(node.parents_nann)

    x_arr     = Array(Vector{Float64}, n_nann_parent_instantiations)
    LHS_arr   = Array(Vector{Float64}, n_nann_parent_instantiations, n_disc_parent_instantiations)
    RHS_arr   = Array(Matrix{Float64}, n_nann_parent_instantiations, n_disc_parent_instantiations)
    σ_var_arr = Array(Variance,        n_nann_parent_instantiations, n_disc_parent_instantiations)

    # dim_node = 0
    for n in 1 : n_nann_parent_instantiations
        n_nann_cont_parents = count_ones(n-1)

        x_len = nparents_cont + n_nann_cont_parents + 1 # one for the bias value
                                                        # note that target_as_parent is already in n_nann_cont_parents

        # x: [cont | nann | target_as_parent | bias]
        x_arr[n] = Array(Float64, x_len)
        x_arr[n][end] = 1.0 # set bias to 1.0
        # dim_node += n_disc_parent_instantiations * (x_len + 1) # one for the standard deviation

        for i in 1 : n_disc_parent_instantiations
            LHS_arr[n, i] = zeros(Float64, x_len)
            RHS_arr[n, i] = diagm(fill(λ, x_len))
            σ_var_arr[n, i] = Variance()
        end
    end

    # perform learning pass through data
    for i in 1 : nsamples(data)
        q = _get_nann_assignment_index!(i, node, data, parental_assignments_nann, parent_instantiation_counts_nann)
        x = x_arr[q]
        y = _pull_cont_data!(i, x, node, data)
        k = _get_disc_assignment_index!(i, node, data, parental_assignments_disc, parent_instantiation_counts_disc)
        _update_ridge_regression!(LHS_arr[q,k], RHS_arr[q,k], x, y)
    end

    # solve ridge regressions
    #  w = (λI + XᵀX)⁻¹ Xᵀy

    w_arr = Array(Vector{Float64}, n_nann_parent_instantiations, n_disc_parent_instantiations)
    for q in 1 : n_nann_parent_instantiations
        n_nann_cont_parents = count_ones(q-1)
        w_len = n_nann_cont_parents + nparents_cont + 1 # one for the 1.0 value

        for k in 1 : n_disc_parent_instantiations
            try
                w_arr[q,k] = RHS_arr[q,k] \ LHS_arr[q,k]
                @assert(findfirst(v->isnan(v), w_arr[q,k]) == 0)
            catch
                w_arr[q,k] = zeros(w_len)
            end

            @assert(findfirst(v->isnan(v), w_arr[q,k]) == 0)
        end
    end

    # perform σ fit pass
    for i in 1 : nsamples(data)
        q = _get_nann_assignment_index!(i, node, data, parental_assignments_nann, parent_instantiation_counts_nann)
        x = x_arr[q]
        y = _pull_cont_data!(i, x, node, data)
        k = _get_disc_assignment_index!(i, node, data, parental_assignments_disc, parent_instantiation_counts_disc)

        w = w_arr[q, k]
        y_pred = dot(w, x)
        fit!(σ_var_arr[q,k], y - y_pred)
    end

    # back out σ
    σ_arr = Array(Float64, n_nann_parent_instantiations, n_disc_parent_instantiations)
    for q in 1 : n_nann_parent_instantiations
        for k in 1 : n_disc_parent_instantiations
            σ_var = σ_var_arr[q,k]
            if StatsBase.nobs(σ_var) > 1
                σ_arr[q, k] = max(std(σ_var), σ_min)
            else
                σ_arr[q, k] = σ_min
            end
        end
    end

    (w_arr, σ_arr, x_arr)
end

function _get_component_score(node::NodeInTraining, data::TrainingData, params::LB_TrainParams)

    # Using BIC score
    # see: “Ideal Parent” Structure Learning for Continuous Variable Bayesian Networks
    #
    # BIC = max_θ  l(D|G, θ) - log(m)/2 * dim(G)
    #
    #   where max_θ  l(D|G, θ) is found using ridge regression
    #         and dim(G) is the number of parameters in G
    #
    #  Calling this function causes it to return:
    #    max_θ  l(D|G, θ) - log(m)/2*dim(node)
    #
    #   which can be added to other components to get the full BIC score

    m = nsamples(data)
    @assert(size(data.X_disc,1) == m)
    @assert(size(data.X_cont,1) == m)
    @assert(size(data.X_nann,1) == m)

    parent_instantiation_counts_nann, n_nann_parent_instantiations, parental_assignments_nann = _init_parent_stuff_nann(node)
    parent_instantiation_counts_disc, n_disc_parent_instantiations, parental_assignments_disc = _init_parent_stuff_disc(node, data)

    w_arr, σ_arr, x_arr = _perform_ridge_regression(node, data, params,
                            parent_instantiation_counts_nann, n_nann_parent_instantiations, parental_assignments_nann,
                            parent_instantiation_counts_disc, n_disc_parent_instantiations, parental_assignments_disc)

    # compute the log likelihood
    logl = 0.0
    for i in 1 : m
        q = _get_nann_assignment_index!(i, node, data, parental_assignments_nann, parent_instantiation_counts_nann)
        x = x_arr[q]
        y = _pull_cont_data!(i, x, node, data)
        k = _get_disc_assignment_index!(i, node, data, parental_assignments_disc, parent_instantiation_counts_disc)
        logl += _calc_logl_contribution(w_arr[q,k], σ_arr[q,k], x, y)
    end

    n_suff_stats = 0
    for w in w_arr
        n_suff_stats += length(w) + 1
    end

    logl - log(m)*n_suff_stats/2
end
function _get_component_score!(node::NodeInTraining, data::TrainingData, params::LB_TrainParams)

    if !haskey(data.score_cache, node)
        data.score_cache[node] = _get_component_score(node, data, params)
    end
    data.score_cache[node]
end

function _greedy_hillclimb_iter_on_node(node::NodeInTraining, data::TrainingData, params::LB_TrainParams, other_node_target_as_parent::Int)

    if params.verbosity > 1
        println("_greedy_hillclimb_iter_on_node: ")
        println(node)
    end

    start_score = data.score_cache[node]
    best_score = start_score
    parent_to_add = -1
    add_sym = :none

    n_indicators_cont = size(data.X_cont, 2)

    if nparents_total(node) ≥ params.max_parents
        if params.verbosity > 1
            println("skipping due to max_parents = $(params.max_parents) reached")
        end
        return (-1.0, parent_to_add, add_sym)
    end

    # try adding discrete edges
    for p in 1 : size(data.X_disc, 2)
        new_score = _get_component_score!(set_new_parent_disc(node, p), data, params)

        if params.verbosity > 1
            @printf("trying disc %5d  -> Δscore = %12.6f\n", p, new_score - start_score)
        end

        if new_score > best_score
            best_score = new_score
            parent_to_add = p
            add_sym = :disc
        end
    end

    # try adding continuous edges
    for p in 1 : n_indicators_cont
        new_score = _get_component_score!(set_new_parent_cont(node, p), data, params)

        if params.verbosity > 1
            @printf("trying cont %5d  -> Δscore = %12.6f\n", p, new_score - start_score)
        end

        if new_score > best_score
            best_score = new_score
            parent_to_add = p
            add_sym = :cont
        end
    end

    # try adding nannable edges
    for p in 1 : size(data.X_nann, 2)
        new_score = _get_component_score!(set_new_parent_nann(node, p), data, params)

        if params.verbosity > 1
            @printf("trying nann %5d  -> Δscore = %12.6f\n", p, new_score - start_score)
        end

        if new_score > best_score
            best_score = new_score
            parent_to_add = p
            add_sym = :nann
        end
    end

    # try adding other targets as parents
    if node.target_as_parent == 0 && other_node_target_as_parent == 0
        p = node.index == 1 ? 2 : 1

        new_score = _get_component_score!(set_new_parent_targ(node, p), data, params)

        if params.verbosity > 1
            @printf("trying targ %5d  -> Δscore = %12.6f\n", p, new_score - start_score)
        end

        if new_score > best_score
            best_score = new_score
            parent_to_add = p
            add_sym = :targ
        end
    end

    Δscore = best_score - start_score

    (Δscore, parent_to_add, add_sym)
end
function _build_linear_gaussian_node(
    node::NodeInTraining,
    data::TrainingData,
    params::LB_TrainParams,
    ind_old_to_new_disc::Dict{Int, Int},
    ind_old_to_new_cont::Dict{Int, Int},
    ind_old_to_new_nann::Dict{Int, Int},
    )

    index = length(ind_old_to_new_cont) + length(ind_old_to_new_nann) + node.index

    nparents_disc = length(node.parents_disc)
    parents_disc = Array(Int, nparents_disc)
    for (i,p) in enumerate(node.parents_disc)
        parents_disc[i] = ind_old_to_new_disc[p]
    end

    nparents_cont = length(node.parents_cont)
    parents_cont = Array(Int, nparents_cont)
    for (i,p) in enumerate(node.parents_cont)
        parents_cont[i] = ind_old_to_new_cont[p]
    end

    nparents_nann = length(node.parents_nann)
    parents_nann = Array(Int, nparents_nann)
    for (i,p) in enumerate(node.parents_nann)
        parents_nann[i] = ind_old_to_new_nann[p] + length(ind_old_to_new_cont)
    end
    if node.target_as_parent != 0
        nparents_nann += 1
        push!(parents_nann, node.target_as_parent+length(ind_old_to_new_cont)+length(ind_old_to_new_nann))
    end

    m = nsamples(data)
    @assert(size(data.X_disc,1) == m)
    @assert(size(data.X_cont,1) == m)
    @assert(size(data.X_nann,1) == m)

    parental_assignments_disc = Array(Int, nparents_disc)
    parental_assignments_cont = Array(Float64, nparents_cont+nparents_nann+1)

    parent_instantiation_counts_nann, n_nann_parent_instantiations, parental_assignments_nann = _init_parent_stuff_nann(node)
    parent_instantiation_counts_disc, n_disc_parent_instantiations, parental_assignments_disc = _init_parent_stuff_disc(node, data)

    w_arr, σ_arr, x_arr = _perform_ridge_regression(node, data, params,
                            parent_instantiation_counts_nann, n_nann_parent_instantiations, parental_assignments_nann,
                            parent_instantiation_counts_disc, n_disc_parent_instantiations, parental_assignments_disc)

    stats = Array(LinearGaussianStats, n_nann_parent_instantiations, n_disc_parent_instantiations)

    for n in 1 : n_nann_parent_instantiations
        for i in 1 : n_disc_parent_instantiations

            w = w_arr[n,i]
            σ = σ_arr[n,i]
            @assert(!isnan(σ))
            stats[n,i] = LinearGaussianStats(w, σ)
        end
    end

    ConditionalLinearGaussianNode(index, stats,
                       parents_disc, parents_cont, parents_nann,
                       parental_assignments_disc,
                       parental_assignments_nann,
                       parental_assignments_cont,
                       parent_instantiation_counts_disc,
                       parent_instantiation_counts_nann,
                    )
end

function train(
    training_data::ModelTrainingData2,
    preallocated_data::LB_PreallocatedData,
    params::LB_TrainParams,
    foldset::FoldSet,
    )

    @assert(findfirst(v->isnan(v), preallocated_data.Y) == 0)
    @assert(findfirst(v->isinf(v), preallocated_data.Y) == 0)

    # -----------------------------------------
    # run structure learning
    #  - always add the next best edge

    node_lat = NodeInTraining(1)
    node_lon = NodeInTraining(2)

    data = TrainingData(training_data, preallocated_data, foldset)

    _get_component_score!(node_lat, data, params)
    _get_component_score!(node_lon, data, params)

    finished = false
    while !finished

        Δscore_lat, parent_to_add_lat, add_sym_lat = _greedy_hillclimb_iter_on_node(node_lat, data, params, node_lon.target_as_parent)
        Δscore_lon, parent_to_add_lon, add_sym_lon = _greedy_hillclimb_iter_on_node(node_lon, data, params, node_lat.target_as_parent)

        if max(Δscore_lat, Δscore_lon) < 0.001
            # no improvement
            finished = true
        else
            if Δscore_lat > Δscore_lon

                node_lat = set_new_parent(node_lat, parent_to_add_lat, add_sym_lat)
                @assert(haskey(data.score_cache, node_lat))
                if params.verbosity > 0
                    println("\n===========================================")
                    @printf("adding to lat: %4d   %4s\n", parent_to_add_lat, string(add_sym_lat))
                    println(node_lat)
                    println("===========================================\n")
                end
            else
                node_lon = set_new_parent(node_lon, parent_to_add_lon, add_sym_lon)
                @assert(haskey(data.score_cache, node_lon))
                if params.verbosity > 0
                    println("\n===========================================")
                    @printf("adding to lon: %4d   %4s\n", parent_to_add_lon, string(add_sym_lon))
                    println(node_lon)
                    println("===========================================\n")
                end
            end
        end
    end

    # -----------------------------------------
    # now build the model

    ind_old_disc = sort!(unique([node_lat.parents_disc; node_lon.parents_disc]))
    extractor_disc = FeatureSubsetExtractor(preallocated_data.features_disc[ind_old_disc])
    ind_old_to_new_disc = Dict{Int,Int}()
    for (o,n) in enumerate(ind_old_disc)
        ind_old_to_new_disc[n] = o
    end

    ind_old_cont = sort!(unique([node_lat.parents_cont; node_lon.parents_cont]))
    ind_old_to_new_cont = Dict{Int,Int}()
    for (o,n) in enumerate(ind_old_cont)
        ind_old_to_new_cont[n] = o
    end

    # println("ind_old_cont: ", ind_old_cont)

    ind_old_nann = sort!(unique([node_lat.parents_nann; node_lon.parents_nann]))
    ind_old_to_new_nann = Dict{Int,Int}()
    for (n,o) in enumerate(ind_old_nann)
        ind_old_to_new_nann[o] = n
    end

    # println("ind_old_nann: ", ind_old_nann)

    extractor_cont = FeatureSubsetExtractor([preallocated_data.features_cont[ind_old_cont];
                                             preallocated_data.features_nann[ind_old_nann]])
    clamper_cont = DataClamper(extractor_cont.x,
                               [preallocated_data.clamper_cont.f_lo[ind_old_cont];
                                preallocated_data.clamper_nann.f_lo[ind_old_nann]],
                               [preallocated_data.clamper_cont.f_hi[ind_old_cont];
                                preallocated_data.clamper_nann.f_hi[ind_old_nann]])

    # println("extractor indicators: ", map(f->symbol(f), extractor_cont.indicators))
    # println("clamper lo: ", clamper_cont.f_lo)

    model_node_lat = _build_linear_gaussian_node(node_lat, data, params,
                                                 ind_old_to_new_disc, ind_old_to_new_cont, ind_old_to_new_nann)
    model_node_lon = _build_linear_gaussian_node(node_lon, data, params,
                                                 ind_old_to_new_disc, ind_old_to_new_cont, ind_old_to_new_nann)

    if params.verbosity > 0
        println("final node lat: ")
        println(model_node_lat)
        println("\n")
        println("final node lon: ")
        println(model_node_lon)
    end

    sample_lat_first = (node_lon.target_as_parent == 0)

    LinearGaussianBayesianNetBehavior(
        params.targets, extractor_disc, extractor_cont,
        clamper_cont, preallocated_data.clamper_act,
        sample_lat_first, model_node_lat, model_node_lon
        )
end


text:
### Calculate the line-of-sight (nautical) miles for each flight and the total travelled by each passenger

NMR.clean_split_dir()
jobid = 4

# this parses records in AComp_Passenger_data.cs into JSON, handling errors
NMR.runjob(NMR.Nmr(jobid, ["AComp_Passenger_data.csv"], NMR.mapper_parserecordacomp, "acomplineofsight.json"))
psgr = NMR.runcombiner("acomplineofsight.json", NMR.combiner_parsejson)

# this gives us list of flights with [flightid, [srcairport, dstairport]]
NMR.runjob(NMR.Nmr(jobid, ["acomplineofsight.json"], NMR.reducer_lineofsightpassenger, "lineofsightpassenger.json"))
# combine results
flights = NMR.runcombiner("lineofsightpassenger.json", NMR.combiner_parsejson, "lineofsightpassenger.json")

# This gives us Airport Lat Long
NMR.runjob(NMR.Nmr(jobid, ["Top30_airports_LatLong.csv"], NMR.mapper_airportlatlon, "airportlatlon.json"))
# combine results
latlon = NMR.runcombiner("airportlatlon.json", NMR.combiner_parsejson, "airportlatlon.json")

function airport_latlon(airportid)
 for a in latlon
  a = JSON.parse(a)
  a["airportcode"] == airportid || continue
  return a["lat"], a["lon"]
 end
end

function getnauticalmiles()
 flight_miles = Dict{String, Float64}()
 for j in flights
  p = JSON.parse(j)
  srcairport = p[2][1]
  dstairport = p[2][2]
  srclatlon = airport_latlon(srcairport)
  dstlatlon = airport_latlon(dstairport)
  if srclatlon === nothing
   warn("Unable to find location of airport $srcairport")
   continue
  end
  if dstlatlon === nothing
   warn("Unable to find location of airport $dstairport")
   continue
  end
  miles = NMR.haversine(airport_latlon(srcairport), airport_latlon(dstairport))
  flight_miles[p[1]] = miles
 end
 flight_miles
end
flight_miles = getnauticalmiles()

passenger_miles = Dict{String, Float64}()
for p in psgr
 j = JSON.parse(p)
 passenger = j["passengerid"]
 flight = j["flightid"]
 if ! (flight in keys(flight_miles))
  continue
 end
 if passenger in keys(passenger_miles)
  passenger_miles[passenger] += flight_miles[flight]
 else
  passenger_miles[passenger] = flight_miles[flight]
 end
end

header = "passenger,miles"
output = ["$(k),$(passenger_miles[k])" for k in keys(passenger_miles)]
NMR.save_job_output(jobid, header, output)


text:
b = SpinBasis(1//2)
up = spinup(b)
down = spindown(b)
# id = identityoperator(b)

@test qubit("010") ≈ qubit("010")
@test qubit("010")' ≈ qubit("010")'
@test (qubit("010")')' ≈ qubit("010")
@test sigmax() ≈ sigmax()
@test sigmax() ≈ full(sigmax())
@test full(sigmax()) ≈ sigmax()
@test full(sigmax()) ≈ full(sigmax())

@test sigmax()' == dagger(sigmax())
@test full(sigmax())' == dagger(full(sigmax()))

@test outer(spinup(b), spinup(b)) == spinup(b) ⊗ dagger(spinup(b))
@test inner(spinup(b), spinup(b)) == dagger(spinup(b)) * spinup(b)
n = 5; @test ghz(n) == (basisstate(b^n, 1) + basisstate(b^n, 2^n)) / sqrt(2)

srand(2018)
for i in 1:10
    tmpr = rand(1:5)
    @test sparse_spinallup_dm(tmpr) |> full == dm(basisstate(b^tmpr, 1))
end
for i in 1:10
    tmpr = rand(1:5)
    @test sparse_spinalldown_dm(tmpr) |> full == dm(basisstate(b^tmpr, 2^tmpr))
end

srand(2018)
@test begin
    outcome, state = measure(hadamard() * up, 1)
    outcome == 1 && state ≈ down
end
begin
    @show outcome, state = measure(hadamard() * up, 1)
    outcome == 0 && state ≈ up
end

@test qubit("010") == up ⊗ down ⊗ up
@test qubit("010")' == dagger(qubit("010"))
@test qudit("012", 3) == basisstate(NLevelBasis(3), 1) ⊗ basisstate(NLevelBasis(3), 2) ⊗ basisstate(NLevelBasis(3), 3)
@test quxit("012", [2,2,3]) == basisstate(NLevelBasis(2), 1) ⊗ basisstate(NLevelBasis(2), 2) ⊗ basisstate(NLevelBasis(3), 3)


text:
#####
##### A simple minmax player to be used as a baseline
#####

"""
A simple implementation of the minmax tree search algorithm, to be used as
a baseline against AlphaZero. Heuristic board values are provided by the
[`GameInterface.heuristic_value`](@ref) function.
"""
module MinMax

import ..GI, ..GameInterface, ..AbstractPlayer, ..think

function current_player_value(white_reward, white_playing) :: Float64
  if iszero(white_reward)
    return 0.
  else
    v = Inf * sign(white_reward)
    return white_playing ? v : - v
  end
end

# Return the value of a state for the player playing
function value(game, depth)
  wr = GI.white_reward(game)
  wp = GI.white_playing(game)
  if isnothing(wr)
    if depth == 0
      return GI.heuristic_value(game)
    else
      return maximum(qvalue(game, a, depth)
        for a in GI.available_actions(game))
    end
  else
    return current_player_value(wr, wp)
  end
end

function qvalue(game, action, depth)
  @assert isnothing(GI.white_reward(game))
  next = copy(game)
  GI.play!(next, action)
  pswitch = GI.white_playing(game) != GI.white_playing(next)
  nextv = value(next, depth - 1)
  return pswitch ? - nextv : nextv
end

minmax(game, actions, depth) = argmax([qvalue(game, a, depth) for a in actions])

"""
    MinMax.Player{Game} <: AbstractPlayer{Game}

A stochastic minmax player, to be used as a baseline.

    MinMax.Player{Game}(;depth, τ=0.)

The minmax player explores the game tree exhaustively at depth `depth`
to build an estimate of the Q-value of each available action. Then, it
chooses an action as follows:

- If there are winning moves (with value `Inf`), one of them is picked
  uniformly at random.
- If all moves are losing (with value `-Inf`), one of them is picked
  uniformly at random.

Otherwise,

- If the temperature `τ` is zero, a move is picked uniformly among those
  with maximal Q-value (there is usually only one choice).
- If the temperature `τ` is nonzero, the probability of choosing
  action ``a`` is proportional to ``e^{\\frac{q_a}{Cτ}}`` where ``q_a`` is the
  Q value of action ``a`` and ``C`` is the maximum absolute value of all
  finite Q values, making the decision invariant to rescaling of
  [`GameInterface.heuristic_value`](@ref).
"""
struct Player{G} <: AbstractPlayer{G}
  depth :: Int
  τ :: Float64
  Player{G}(;depth, τ=0.) where G = new{G}(depth, τ)
end

function think(p::Player, game, turn=nothing)
  actions = GI.available_actions(game)
  n = length(actions)
  qs = [qvalue(game, a, p.depth) for a in actions]
  winning = findall(==(Inf), qs)
  if isempty(winning)
    notlosing = findall(>(-Inf), qs)
    best = argmax(qs)
    if isempty(notlosing)
      π = ones(n)
    elseif iszero(p.τ)
      π = zeros(n)
      all_best = findall(==(qs[best]), qs)
      π[all_best] .= 1.
    else
      qmax = qs[best]
      @assert qmax > -Inf
      C = maximum(abs(qs[a]) for a in notlosing) + eps()
      π = exp.((qs .- qmax) ./ C)
      π .^= (1 / p.τ)
    end
  else
    π = zeros(n)
    π[winning] .= 1.
  end
  π ./= sum(π)
  return actions, π
end

end


text:
using PyCall
@pyimport matplotlib.pyplot as plt
@pyimport matplotlib.lines as lines
@pyimport matplotlib.animation as anim

export SeisPlot,
       waveAnim

include("SeisPlot.jl")
include("waveAnim.jl")

max_stars_repo_name: WneleiGao/ElasticWave


text:
export crossentropy

doc"""
    crossentropy(p, q)

Cross-entropy function between p and q.

```math
f(x) = -\sum_{x} p(x) \log q(x)
```

* p::Var: `Var` of Vector{Int} or Matrix{Float}. If p is `Vector{Int}` and p[i] == 0, returns 0.
* q::Var: `Var` of Matrix{Float}

```julia
p = Var(rand(0:10,5))
q = softmax(Var(rand(Float32,10,5)))
y = crossentropy(p, q)
```
"""
function crossentropy(p::Var, q::Var)
    configure!(p, q)
    Var(crossentropy(p.data,q.data), (crossentropy,p,q))
end

function crossentropy(p::Vector{Int}, q::Matrix{T}) where T
    y = Array{T}(length(p))
    @inbounds for i = 1:length(p)
        y[i] = p[i] > 0 ? -log(q[p[i],i]) : T(0)
    end
    y
end

function crossentropy(p::Vector{T}, q::Vector{T}) where T
    @assert length(p) == length(q)
    y = T[0]
    @inbounds for i = 1:length(p)
        y[1] -= p[i] * log(q[i])
    end
    y
end

function addgrad!(y::Var, ::typeof(crossentropy), p::Var, q::Var)
    @assert isvoid(p.grad)
    isvoid(q.grad) && return
    ∇crossentropy!(y.grad, p.data, q.data, q.grad)
end

function ∇crossentropy!(gy::Vector{T}, p::Vector{I}, q::Matrix{T}, gq::Matrix{T}) where {I<:Integer,T}
    @inbounds for i = 1:length(p)
        if p[i] > 0
            gq[p[i],i] -= gy[i] / q[p[i],i]
        end
    end
end

function ∇crossentropy!(gy::Vector{T}, p::Vector{T}, q::Vector{T}, gq::Vector{T}) where T
    @assert length(gy) == 1
    @inbounds for i = 1:length(p)
        gq[i] -= gy[1] * p[i] / q[i]
    end
end


text:
using IntegratedHU
using Documenter

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

makedocs(;
    modules=[IntegratedHU],
    authors="Dale <djblack@uci.edu> and contributors",
    repo="https://github.com/Dale-Black/IntegratedHU.jl/blob/{commit}{path}#{line}",
    sitename="IntegratedHU.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Dale-Black.github.io/IntegratedHU.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Dale-Black/IntegratedHU.jl",
)


text:
@testset "Element" begin
    # Need some Nodes
    node1 = Node{2}(1,   Vec{2}((2.3, 4.5)), [1,7,9], 3)
    node3 = Node{2}(3,   Vec{2}((2.3, 4.5)), [1,7,9], 3)
    node13 = Node{2}(13, Vec{2}((2.3, 4.5)), [1,7,9], 3)
    node18 = Node{2}(18, Vec{2}((2.3, 4.5)), [1,7,9], 3)

    # Test Full Construction
    element = Element(11, [node1, node3, node13, node18], [12,17,81], 3)

    @test element.id == 11
    @test element.nodes[2].id == 3
    @test element.dofs == [12,17,81]

    # Test DofObject Interface
    @test dofs(element) == [12,17,81]
    @test processor_id(element) == 3
end


text:
function integrator(dudt::Function, u0, tspan, p, solver, args...; kwargs...)
    prob = ODEProblem(dudt, u0, tspan, p)
    return init(prob, alg, args...; kwargs...)
end


text:
using Test, FinancialMonteCarlo;
@show "Black Scholes Model"
S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
d_ = FinancialMonteCarlo.Curve([0.009999, 0.01], T);
D = 90.0;

Nsim = 10000;
Nstep = 30;
sigma = 0.2;
mc = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.SobolMode());
toll = 0.8

spotData1 = ZeroRate(r);

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
EUBin = BinaryEuropeanOption(T, K)
AMData = AmericanOption(T, K)
AmBin = BinaryEuropeanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T)
AsianFixedStrikeData = AsianFixedStrikeOption(T, K)
Model = BlackScholesProcess(sigma, Underlying(S0, d));

display(Model)

@show FwdPrice = pricer(Model, spotData1, mc, FwdData);
@show EuPrice = pricer(Model, spotData1, mc, EUData);
@show EuBinPrice = pricer(Model, spotData1, mc, EUBin);
@show AmPrice = pricer(Model, spotData1, mc, AMData);
@show AmBinPrice = pricer(Model, spotData1, mc, AmBin);
@show BarrierPrice = pricer(Model, spotData1, mc, BarrierData);
@show AsianPrice1 = pricer(Model, spotData1, mc, AsianFloatingStrikeData);

@test abs(FwdPrice - 99.1078451563562) < toll
@test abs(EuPrice - 8.43005524824866) < toll
@test abs(AmPrice - 8.450489415187354) < toll
@test abs(BarrierPrice - 7.5008664470880735) < toll
@test abs(AsianPrice1 - 4.774451704549382) < toll


text:
### Transition model  ############################################################################
function POMDPs.transition(pomdp::SingleOCFPOMDP, s::SingleOCFState, a::SingleOCFAction, dt::Float64 = pomdp.ΔT)

    ## ego transition
    # longitudinal motion
    v_ego = s.ego_v + a.acc*pomdp.ΔT
    sp_ego_v = clamp(v_ego, 0, v_ego)
    x_delta_ego = s.ego_v*pomdp.ΔT + 0.5*a.acc*pomdp.ΔT^2

    # lateral motion, 0.5m in 0.5s should be fine
    y_delta_ego = a.lateral_movement * pomdp.ΔT  # simplified movement, not exactly correct -> a=1, t = 1 --> 1m movement to the left
    sp_ego_y = s.ego_y + y_delta_ego
    sp_ego_y = clamp(sp_ego_y, pomdp.EGO_Y_MIN, pomdp.EGO_Y_MAX )

    ## pedestrian transition
    # absent or not
    if ( !is_state_absent(pomdp, s) )   # pedestrian is not absent

        states = SingleOCFState[]
        sizehint!(states, 50);
        probs = Float64[]
        sizehint!(probs, 50);

        # ped transition --> variation in acceleration and theta
        for a_ped in pomdp.PED_A_RANGE
            sp_ped_v = s.ped_v + a_ped * pomdp.ΔT
            sp_ped_v = clamp(sp_ped_v, 0, pomdp.PED_V_MAX)      # speed limitation
            if a_ped == 0
                a_ped_prob = 3.0
            else
                a_ped_prob = 1.0
            end

            for theta_n in pomdp.PED_THETA_NOISE
                sp_ped_theta = s.ped_theta + theta_n
                @fastmath begin
                sp_ped_s = s.ped_s + sp_ped_v*cos(sp_ped_theta) * pomdp.ΔT - x_delta_ego
                sp_ped_T = s.ped_T + sp_ped_v*sin(sp_ped_theta) * pomdp.ΔT - y_delta_ego  # a_lat corresponds to lateral velocity --> a_lat == v_lat
                end
                sp_ped_s = clamp(sp_ped_s, pomdp.S_MIN, pomdp.S_MAX)
                sp_ped_T = clamp(sp_ped_T, pomdp.T_MIN, pomdp.T_MAX)

           #=
                (ego_y_state_space,ego_v_state_space) = getEgoDataInStateSpace(pomdp, sp_ego_y, sp_ego_v)
                state_vector = SVector(sp_ped_s, sp_ped_T, sp_ped_theta, sp_ped_v) # looks faster
                @inbounds ind, weight = interpolants(pomdp.state_space_grid_ped, state_vector)
                for i=1:length(ind)
                    if weight[i] > 0.1
                        state = pomdp.state_space_ped[ind[i]]
                        push!(states, SingleOCFState(ego_y_state_space, ego_v_state_space, state.ped_s, state.ped_T, state.ped_theta, state.ped_v))
                        push!(probs, weight[i]*a_ped_prob)
                    end
                end
=#
                # find the continious state in the state space
                state_vector = SVector(sp_ego_y, sp_ego_v, sp_ped_s, sp_ped_T, sp_ped_theta, sp_ped_v) # looks faster
                @inbounds ind, weight = interpolants(pomdp.state_space_grid, state_vector)
                for i=1:length(ind)
                    if weight[i] > 0.1
                        state = pomdp.state_space[ind[i]]
                        push!(states, state)
                        push!(probs, weight[i]*a_ped_prob)
                    end
                end

            end
        end

        # add roughening
        if length(probs) > 1
            normalize!(probs, 1)
            probs .+= maximum(probs)
            normalize!(probs,1)
        end
        return SingleOCFBelief(states,probs)

    else
        # pedestrian is absent
        return initBeliefAbsentPedestrianBorder(pomdp, sp_ego_y, sp_ego_v)
    end

end


text:
"""
Pads so that the length of the result is a multiple of block size in bytes.
A block size of 10 represents 10 bytes that make up a single block.

Padding for block ciphers is specified in the PKCS#7 standard and in RFC 5652.
Method can be
- :zero
- :pkcs7
- :iso

"""
function pad(bytes::Vector{UInt8}, block_size::Int64, method=:pkcs7)
    length_to_add = block_size - (length(bytes) % block_size)

    padding = repeat([0x00], length_to_add)
    if method == :zero
        padding = repeat([0x00], length_to_add)
    elseif method == :pkcs7
        padding = repeat([UInt8(length_to_add)], length_to_add)
    elseif method == :iso
        padding = vcat(0x80, repeat([0x00], length_to_add-1))
    elseif method == :ansi
        padding = vcat(repeat([0x00], length_to_add-1), [UInt8(length_to_add)])
    end
    return vcat(bytes, padding)
end

"""
This is not vulnerable to padded oracle attacks because it does not fail
if the padding is invalid.
"""
function unpad(bytes::Vector{UInt8}, method=:pkcs7)
    if method == :zero
        throw(PaddingError(method, "Does not support unpad."))
    elseif method == :pkcs7
        return bytes[1:(end-bytes[end])]
    elseif method == :iso
        while padding[end] != 0x80
            padding = padding[1:(end-1)]
        end
        return padding[1:(end-1)]
    elseif method == :ansi
        return bytes[1:(end-bytes[end])]
    end
end

"""
Galois Field (256) Multiplication of two Bytes
"""
function gmul(a::UInt8, b::UInt8)
    p = 0x00
    for _ in 1:8
        if b & 1 != 0x00
            p = xor(p, a)
        end

        hi_bit_set = (a & 0x80) != 0x00
        a <<= 0x01
        if hi_bit_set
            a = xor(a, 0x1b)
        end
        b >>= 0x01
    end
    return p
end

function eval_poly(polynomial, x, modulus)
    accum = zero(eltype(polynomial))
    for coeff in reverse(polynomial)
        accum *= x
        accum += coeff
        accum %= modulus
    end
    return accum
end

function lagrange_interpolate(x, xs, ys, modulus)
    k = length(xs)
    nums = []
    dens = []
    for i in 1:k
        curr = xs[i]
        others = vcat(xs[1:(i-1)], xs[(i+1):end])
        push!(nums, prod(x .- others))
        push!(dens, prod(curr .- others))
    end

    den = prod(dens)
    num = sum(nums[i] * den * ys[i] * invmod(dens[i], modulus) % modulus for i in 1:k)
    return (num * invmod(den, modulus) % modulus + modulus) % modulus
end


text:
# mathematical operators between two Timestamp objects
for op in [:+, :-, :*, :/]
    @eval begin
        function ($op)(ts1::Timestamp, ts2::Timestamp)
            tstamp = abs(ts1.timestamp - ts2.timestamp)
            vals   = ($op)(ts1.value, ts2.value)
            Timestamp(tstamp, vals)
        end # function
    end # eval
end # loop

# mathematical operators between a Timestamp object and Int,Float64
for op in [:+, :-, :*, :/]
    @eval begin
        function ($op)(ts1::Timestamp, num::Union(Int,Float64))
            vals = ($op)(ts1.value, num)
            Timestamp(ts1.timestamp, vals)
        end # function
    end # eval
end # loop

# comparison operators between two Timestamp objects
#for op in [:>, :<, :==, :>=, :<=]
for op in [:>, :<, :>=, :<=]
    @eval begin
        function ($op)(ts1::Timestamp, ts2::Timestamp)
            tstamp = abs(ts1.timestamp - ts2.timestamp)
            vals   = ($op)(ts1.value, ts2.value)
            Timestamp(tstamp, vals)
        end # function
    end # eval
end # loop

# comparison operators between a Timestamp object and Int,Float64
#for op in [:>, :<, :==, :>=, :<=]
for op in [:>, :<, :>=, :<=]
    @eval begin
        function ($op)(ts1::Timestamp, num::Union(Int,Float64))
            vals = ($op)(ts1.value, num)
            Timestamp(ts1.timestamp, vals)
        end # function
    end # eval
end # loop


text:
######################################
"""
    plt = plotsamr2d(tiles::AbstractVector{VisClaw.AMRGrid}; AMRlevel=[], wind::Bool=false, region="", kwargs...)

    plotsamr2d!(plt::Plots.Plot, tiles::AbstractVector{VisClaw.AMRGrid}; AMRlevel=[], wind::Bool=false, region="", kwargs...)

Function: plot values of AMR grids in two-dimension
"""
function plotsamr2d!(plt, tiles::AbstractVector{VisClaw.AMRGrid}; AMRlevel=[], wind::Bool=false, region="", kwargs...)

    # check arg
    if isa(tiles[1], VisClaw.SurfaceHeight)
        var = :eta
    elseif isa(tiles[1], VisClaw.Velocity)
        var = :vel
    elseif isa(tiles[1], VisClaw.Storm)
        if wind
            var = :u
        else
            var = :slp
        end
    else
        error("Invalid argument")
    end

    # parse keyword args
    kwdict = KWARG(kwargs)
    # -----------------------------
    # linetype
    seriestype, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_seriestype, :heatmap)
    # -----------------------------
    # color
    seriescolor, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_seriescolor, :auto)
    # -----------------------------
    # color axis
    clims, kwdict = VisClaw.parse_clims(kwdict)
    if clims == nothing
        vals = getfield.(tiles, var)
        clims = (
                 minimum(minimum.(v->isnan(v) ?  Inf : v, vals)),
                 maximum(maximum.(v->isnan(v) ? -Inf : v, vals))
                 )
        #clims = (minimum(minimum.(filter.(!isnan, vals))),
        #         maximum(maximum.(filter.(!isnan, vals))))
    end
    # -----------------------------
    # background_color_inside
    bginside, kwdict =  VisClaw.kwarg_default(kwdict, VisClaw.parse_bgcolor_inside, Plots.RGB(.7,.7,.7))
    if seriestype == :surface; bginside = Plots.RGBA(0.0,0.0,0.0,0.0); end
    # -----------------------------
    # colorbar
    cb, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_colorbar, :none)
    # -----------------------------
    # colorbar_title
    cbtitle, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_colorbar_title, "")
    # -----------------------------
    # xlims, ylims
    xlims, kwdict = VisClaw.parse_xlims(kwdict)
    ylims, kwdict = VisClaw.parse_ylims(kwdict)
    # -----------------------------
    ## set range
    if isa(region, VisClaw.AbstractTopo); xlims=extrema(region.x); ylims=extrema(region.y); end
    if isa(region, VisClaw.Region); xlims=region.xlims; ylims=region.ylims; end

    # Too fine grids are not plotted
    if isempty(AMRlevel) && xlims==nothing && ylims==nothing
        AMRlevel = 1:3
    end

    ## the number of tiles
    ntile = length(tiles)
    ##
    nplot_org = plt.n

    ## display each tile
    for i = 1:ntile

        ## skip when the AMR level of this tile doesn't match any designated level
        if !isempty(AMRlevel)
            if isempty(findall(tiles[i].AMRlevel .== AMRlevel)); continue; end
        end

        ## set the boundary
        x = [tiles[i].xlow, tiles[i].xlow+tiles[i].dx*tiles[i].mx]
        y = [tiles[i].ylow, tiles[i].ylow+tiles[i].dy*tiles[i].my]

        ## grid info
        xvec = collect(Float64, x[1]-0.5tiles[i].dx:tiles[i].dx:x[2]+0.5tiles[i].dx+1e-4)
        yvec = collect(Float64, y[1]-0.5tiles[i].dy:tiles[i].dy:y[2]+0.5tiles[i].dy+1e-4)

        ## check whether the tile is on the domain
        if xlims!=nothing
            if (xvec[end] < xlims[1]) | (xlims[2] < xvec[1]); continue; end
        end
        if ylims!=nothing
            if (yvec[end] < ylims[1]) | (ylims[2] < yvec[1]); continue; end
        end

        ## adjust data
        val = zeros(tiles[i].my+2,tiles[i].mx+2)
        if !wind
            val[2:end-1,2:end-1] = getfield(tiles[i], var)
        else
            val[2:end-1,2:end-1] = sqrt.(getfield(tiles[i], :u).^2 .+ getfield(tiles[i], :v).^2)
        end
        val[2:end-1,1] = val[2:end-1,2]
        val[2:end-1,end] = val[2:end-1,end-1]
        val[1,:] = val[2,:]
        val[end,:] = val[end-1,:]

        ## plot
        plt = Plots.plot!(plt, xvec, yvec, val, seriestype=seriestype, c=seriescolor, clims=clims, colorbar=false)

    end

    # check the number of added series
    if plt.n == nplot_org
        println("Nothing to be plotted. Check xlims and ylims you specified.")
        plt = Plots.plot!(plt, 0:5, 0:5, repeat([NaN], inner=(6,6)), label="", c=seriescolor, clims=clims, colorbar=cb)
        return plt
    end

    ## xlims, ylims
    x1, x2, y1, y2 = VisClaw.getlims(tiles)
    xrange = (x1, x2)
    yrange = (y1, y2)
    xlims = xlims==nothing ? xrange : xlims
    ylims = ylims==nothing ? yrange : ylims
    plt = Plots.plot!(plt, xlims=xlims, ylims=ylims)

    if !isempty(kwdict); plt = Plots.plot!(plt; kwdict...); end

    ## colorbar
    for i = nplot_org+1:plt.n; plt.series_list[i].plotattributes[:colorbar_entry] = false; end
    if (cb !== :none) && (cb !== false)
        plt.series_list[nplot_org+1].plotattributes[:colorbar_entry] = true
    end

    ## Appearance
    plt = Plots.plot!(plt, axis_ratio=:equal, grid=false, bginside=bginside, colorbar=cb, colorbar_title=cbtitle)

    ## return value
    return plt
end
######################################
"""
$(@doc plotsamr2d!)
"""
plotsamr2d(tiles; kwargs...) = plotsamr2d!(Plots.plot(), tiles; kwargs...)
######################################

#######################################
"""
    gridnumber!(plt::Plots.Plot, tiles; AMRlevel::AbstractVector=[], font::Plots.Font=Plots.font(12, :hcenter, :black), xlims=(), ylims=())

Function: add the grid numbers of tiles
"""
function gridnumber!(plt, tiles; AMRlevel::AbstractVector=[],
                     font::Plots.Font=Plots.font(12, :hcenter, :black),
                     xlims=(), ylims=())

    # Too fine grids are not plotted
    if isempty(AMRlevel) && isempty(xlims) && isempty(ylims)
        AMRlevel = 1:3
    end

    # get current lims
    if plt.n != 0
        if isempty(xlims)
            xlims = Plots.xlims(plt)
        end
        if isempty(ylims)
            ylims = Plots.ylims(plt)
        end
    end

    ## the number of tiles
    ntile = length(tiles)
    for i = 1:ntile

        ## skip when the AMR level of this tile doesn't match any designated level
        if !isempty(AMRlevel)
            if isempty(findall(tiles[i].AMRlevel .== AMRlevel)); continue; end
        end

        ## set the boundary
        x = [tiles[i].xlow, tiles[i].xlow+tiles[i].dx*tiles[i].mx]
        y = [tiles[i].ylow, tiles[i].ylow+tiles[i].dy*tiles[i].my]
        ann = @sprintf("%02d", tiles[i].gridnumber)

        ## check whether the tile is on the domain
        if !isempty(xlims)
            if (mean(x) < xlims[1]) | (xlims[2] < mean(x)); continue; end
        end
        if !isempty(ylims)
            if (mean(y) < ylims[1]) | (ylims[2] < mean(y)); continue; end
        end

        ## annotations
        plt = Plots.plot!(plt; annotations=(mean(x),mean(y), Plots.text(ann, font)))
    end
    return plt
end
#######################################

#######################################
"""
    tilebound!(plt, tiles; AMRlevel=[], kwargs...)

Function: draw boundaries of tiles
"""
function tilebound!(plt, tiles; AMRlevel=[], kwargs...)

    # parse keyword args
    kwdict = KWARG(kwargs)

    # -----------------------------
    # linestyle
    linestyle, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_linestyle, :solid)
    # -----------------------------
    # linecolor
    linecolor, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_linecolor, :black)
    # -----------------------------
    # xlims, ylims
    xlims, kwdict = VisClaw.parse_xlims(kwdict)
    ylims, kwdict = VisClaw.parse_ylims(kwdict)
    # get current lims
    if plt.n != 0
        if xlims == nothing
            xlims = Plots.xlims(plt)
        end
        if ylims == nothing
            ylims = Plots.ylims(plt)
        end
    end
    # -----------------------------

    # Too fine grids are not plotted
    if isempty(AMRlevel) && xlims==nothing && ylims==nothing; AMRlevel = 1:3; end

    ## the number of tiles
    ntile = length(tiles)
    for i = 1:ntile

        ## skip when the AMR level of this tile doesn't match any designated level
        if !isempty(AMRlevel); if isempty(findall(tiles[i].AMRlevel .== AMRlevel)); continue; end; end

        ## set the boundary
        x = [tiles[i].xlow, tiles[i].xlow+tiles[i].dx*tiles[i].mx]
        y = [tiles[i].ylow, tiles[i].ylow+tiles[i].dy*tiles[i].my]

        ## check whether the tile is on the domain
        if xlims!=nothing
            if (x[2] < xlims[1]) | (xlims[2] < x[1]); continue; end
        end
        if ylims!=nothing
            if (y[2] < ylims[1]) | (ylims[2] < y[1]); continue; end
        end

        plt = Plots.plot!(plt,
              [x[1], x[1], x[2], x[2], x[1]],
              [y[1], y[2], y[2], y[1], y[1]];
              label="", linestyle=linestyle, linecolor=linecolor, kwdict...)
    end
    return plt
end
#######################################

#######################################
"""
    plts = plotsamr(amrs::VisClaw.AMR, timeindex=1:amrs.nstep; kwargs...)

Function: plot VisClaw.AMR data. The keyword arguments follow [`plotsamr2d`](@ref)
"""
function plotsamr(amrs::VisClaw.AMR, timeindex=1:amrs.nstep; kwargs...)
    ## plot time-series
    plt = Array{Plots.Plot}(undef, length(timeindex))
    for i = timeindex
        if isempty(amrs.amr[i]); plt[i-timeindex[1]+1] = Plots.plot(); continue; end
        plt[i-timeindex[1]+1] = VisClaw.plotsamr2d(amrs.amr[i]; kwargs...)
    end
    ## return plots
    return plt
end
#############################################


text:
module FMIImport

using FMICore

using FMICore: fmi2Component, fmi3Instance

# functions that have (currently) no better place

# Receives one or an array of values and converts it into an Array{typeof(value)} (if not already).
function prepareValue(value)
    if isa(value, Array) && length(size(value)) == 1
        return value
    else
        return [value]
    end

    @assert false "prepareValue(...): Unknown dimension of structure `$dim`."
end
export prepareValue, prepareValueReference

# wildcards for how a user can pass a fmi[X]ValueReference
fmi2ValueReferenceFormat = Union{Nothing, String, Array{String,1}, fmi2ValueReference, Array{fmi2ValueReference,1}, Int64, Array{Int64,1}, Symbol}
fmi3ValueReferenceFormat = Union{Nothing, String, Array{String,1}, fmi3ValueReference, Array{fmi3ValueReference,1}, Int64, Array{Int64,1}}
export fmi2ValueReferenceFormat, fmi3ValueReferenceFormat

### FMI2 ###

include("FMI2_convert.jl")
include("FMI2_c.jl")
include("FMI2_int.jl")
include("FMI2_ext.jl")
include("FMI2_md.jl")
include("FMI2_fmu_to_md.jl")

# FMI2_c.jl
export fmi2CallbackLogger, fmi2CallbackAllocateMemory, fmi2CallbackFreeMemory, fmi2CallbackStepFinished
export fmi2ComponentState, fmi2ComponentStateModelSetableFMUstate, fmi2ComponentStateModelUnderEvaluation, fmi2ComponentStateModelInitialized
export fmi2Instantiate, fmi2FreeInstance!, fmi2GetTypesPlatform, fmi2GetVersion
export fmi2SetDebugLogging, fmi2SetupExperiment, fmi2EnterInitializationMode, fmi2ExitInitializationMode, fmi2Terminate, fmi2Reset
export fmi2GetReal!, fmi2SetReal, fmi2GetInteger!, fmi2SetInteger, fmi2GetBoolean!, fmi2SetBoolean, fmi2GetString!, fmi2SetString
export fmi2GetFMUstate!, fmi2SetFMUstate, fmi2FreeFMUstate!, fmi2SerializedFMUstateSize!, fmi2SerializeFMUstate!, fmi2DeSerializeFMUstate!
export fmi2GetDirectionalDerivative!, fmi2SetRealInputDerivatives, fmi2GetRealOutputDerivatives
export fmi2DoStep, fmi2CancelStep, fmi2GetStatus!, fmi2GetRealStatus!, fmi2GetIntegerStatus!, fmi2GetBooleanStatus!, fmi2GetStringStatus!
export fmi2SetTime, fmi2SetContinuousStates, fmi2EnterEventMode, fmi2NewDiscreteStates!, fmi2EnterContinuousTimeMode, fmi2CompletedIntegratorStep!
export fmi2GetDerivatives!, fmi2GetEventIndicators!, fmi2GetContinuousStates!, fmi2GetNominalsOfContinuousStates!

# FMI2_convert.jl
export fmi2StringToValueReference, fmi2ValueReferenceToString, fmi2ModelVariablesForValueReference
export fmi2GetSolutionState, fmi2GetSolutionTime, fmi2GetSolutionValue

# FMI2_int.jl
# almost everything exported in `FMI2_c.jl`
export fmi2GetReal, fmi2GetInteger, fmi2GetString, fmi2GetBoolean
export fmi2GetFMUstate, fmi2SerializedFMUstateSize, fmi2SerializeFMUstate, fmi2DeSerializeFMUstate
export fmi2GetDirectionalDerivative
export fmi2GetStartValue, fmi2SampleDirectionalDerivative, fmi2CompletedIntegratorStep

# FMI2_ext.jl
export fmi2Unzip, fmi2Load, loadBinary, fmi2Reload, fmi2Unload, fmi2Instantiate!
export fmi2SampleDirectionalDerivative!
export fmi2GetJacobian, fmi2GetJacobian!, fmi2GetFullJacobian, fmi2GetFullJacobian!
export fmi2Get, fmi2Get!, fmi2Set

# FMI2_md.jl
export fmi2LoadModelDescription
export fmi2GetDefaultStartTime, fmi2GetDefaultStopTime, fmi2GetDefaultTolerance, fmi2GetDefaultStepSize
export fmi2GetModelName, fmi2GetGUID, fmi2GetGenerationTool, fmi2GetGenerationDateAndTime, fmi2GetVariableNamingConvention, fmi2GetNumberOfEventIndicators, fmi2GetNumberOfStates, fmi2IsCoSimulation, fmi2IsModelExchange
export fmi2DependenciesSupported, fmi2GetModelIdentifier, fmi2CanGetSetState, fmi2CanSerializeFMUstate, fmi2ProvidesDirectionalDerivative

# FMI2_fmu_to_md.jl
# everything exported in `FMI2_md.jl`

### FMI3 ###

include("FMI3_c.jl")
include("FMI3_convert.jl")
include("FMI3_int.jl")
include("FMI3_ext.jl")
include("FMI3_md.jl")
include("FMI3_fmu_to_md.jl")

# FMI3_c.jl
export fmi3CallbackLogger, fmi3CallbackIntermediateUpdate, fmi3CallbackClockUpdate
export fmi3InstanceState
export fmi3InstantiateModelExchange, fmi3InstantiateCoSimulation, fmi3InstantiateScheduledExecution, fmi3FreeInstance!
export fmi3GetVersion, fmi3SetDebugLogging
export fmi3EnterInitializationMode, fmi3ExitInitializationMode, fmi3Terminate, fmi3Reset
export fmi3GetFloat32!, fmi3SetFloat32, fmi3GetFloat64!, fmi3SetFloat64
export fmi3GetInt8!, fmi3SetInt8, fmi3GetUInt8!, fmi3SetUInt8, fmi3GetInt16!, fmi3SetInt16, fmi3GetUInt16!, fmi3SetUInt16, fmi3GetInt32!, fmi3SetInt32, fmi3GetUInt32!, fmi3SetUInt32, fmi3GetInt64!, fmi3SetInt64, fmi3GetUInt64!, fmi3SetUInt64
export fmi3GetBoolean!, fmi3SetBoolean, fmi3GetString!, fmi3SetString, fmi3GetBinary!, fmi3SetBinary, fmi3GetClock!, fmi3SetClock
export fmi3GetFMUState!, fmi3SetFMUState, fmi3FreeFMUState!, fmi3SerializedFMUStateSize!, fmi3SerializeFMUState!, fmi3DeSerializeFMUState!
export fmi3SetIntervalDecimal, fmi3SetIntervalFraction, fmi3GetIntervalDecimal!, fmi3GetIntervalFraction!, fmi3GetShiftDecimal!, fmi3GetShiftFraction!
export fmi3ActivateModelPartition
export fmi3GetNumberOfVariableDependencies!, fmi3GetVariableDependencies!
export fmi3GetDirectionalDerivative!, fmi3GetAdjointDerivative!, fmi3GetOutputDerivatives!
export fmi3EnterConfigurationMode, fmi3ExitConfigurationMode, fmi3GetNumberOfContinuousStates!, fmi3GetNumberOfEventIndicators!, fmi3GetContinuousStates!, fmi3GetNominalsOfContinuousStates!
export fmi3EvaluateDiscreteStates, fmi3UpdateDiscreteStates, fmi3EnterContinuousTimeMode, fmi3EnterStepMode
export fmi3SetTime, fmi3SetContinuousStates, fmi3GetContinuousStateDerivatives, fmi3GetEventIndicators!, fmi3CompletedIntegratorStep!, fmi3EnterEventMode, fmi3DoStep!

# FMI3_convert.jl
export fmi3StringToValueReference, fmi3ValueReferenceToString, fmi3ModelVariablesForValueReference

# FMI3_int.jl
export fmi3GetFloat32, fmi3GetFloat64
export fmi3GetInt8, fmi3GetUInt8, fmi3GetInt16, fmi3GetUInt16, fmi3GetInt32, fmi3GetUInt32, fmi3GetInt64, fmi3GetUInt64
export fmi3GetBoolean, fmi3GetString, fmi3GetBinary, fmi3GetClock
export fmi3GetFMUState, fmi3SerializedFMUStateSize, fmi3SerializeFMUState, fmi3DeSerializeFMUState
export fmi3GetDirectionalDerivative

# FMI3_ext.jl
export fmi3Unzip, fmi3Load, fmi3Unload, fmi3InstantiateModelExchange!, fmi3InstantiateCoSimulation!, fmi3InstantiateScheduledExecution!

# FMI3_md.jl
export fmi3LoadModelDescription
export fmi3GetModelName, fmi3GetInstantiationToken, fmi3GetGenerationTool, fmi3GetGenerationDateAndTime, fmi3GetVariableNamingConvention
export fmi3IsCoSimulation, fmi3IsModelExchange, fmi3IsScheduledExecution
export fmi3GetModelIdentifier, fmi3CanGetSetState, fmi3CanSerializeFMUstate, fmi3ProvidesDirectionalDerivatives, fmi3ProvidesAdjointDerivatives

# FMI3_fmu_to_md.jl
# everything exported in `FMI2_md.jl`

###

fmi2Struct = Union{FMU2, FMU2Component}
fmi3Struct = Union{FMU3, FMU3Component}
export fmi2Struct, fmi3Struct

fmi2StructMD = Union{FMU2, FMU2Component, fmi2ModelDescription}
fmi3StructMD = Union{FMU3, FMU3Component, fmi3ModelDescription}
export fmi2StructMD, fmi3StructMD

end # module


text:
struct PlusMinusBridge{T, F <: MOI.AbstractVectorFunction,
                       G <: MOI.AbstractVectorFunction,
                       ST <: MOI.AbstractVectorSet} <: MOIB.Constraint.AbstractBridge
    plus::MOI.ConstraintIndex{F, ST}
    minus::MOI.ConstraintIndex{G, ST}
end

function MOIB.Constraint.bridge_constraint(
    ::Type{PlusMinusBridge{T, F, G, ST}},
    model::MOI.ModelLike,
    f::MOI.AbstractVectorFunction,
    set::PlusMinusSet{ST}) where {T, F, G, ST}

    plus = MOI.add_constraint(model, f, set.set)
    minus = MOI.add_constraint(model, MOIU.operate(-, T, f), set.set)
    return PlusMinusBridge{T, F, G, ST}(plus, minus)
end

function MOI.supports_constraint(::Type{PlusMinusBridge{T}},
                                 ::Type{<:MOI.AbstractVectorFunction},
                                 ::Type{<:PlusMinusSet}) where T
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:PlusMinusBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{<:PlusMinusBridge{T, F, F, ST}}) where {T, F, ST}
    return [(F, ST)]
end
function MOIB.added_constraint_types(::Type{<:PlusMinusBridge{T, F, G, ST}}) where {T, F, G, ST}
    return [(F, ST), (G, ST)]
end
function MOIB.Constraint.concrete_bridge_type(
    ::Type{<:PlusMinusBridge{T}},
    F::Type{<:MOI.AbstractVectorFunction},
    ::Type{<:PlusMinusSet{ST}}) where {T, ST}

    G = MOI.Utilities.promote_operation(-, T, F)
    return PlusMinusBridge{T, F, G, ST}
end

# Attributes, Bridge acting as an model
function MOI.get(::PlusMinusBridge{T, F, F, ST},
                 ::MOI.NumberOfConstraints{F, ST}) where {T, F, ST}
    return 2
end
function MOI.get(::PlusMinusBridge{T, F, G, ST},
                 ::MOI.NumberOfConstraints{F, ST}) where {T, F, G, ST}
    return 1
end
function MOI.get(::PlusMinusBridge{T, F, G, ST},
                 ::MOI.NumberOfConstraints{G, ST}) where {T, F, G, ST}
    return 1
end
function MOI.get(bridge::PlusMinusBridge{T, F, F, ST},
                 ::MOI.ListOfConstraintIndices{F, ST}) where {T, F, ST}
    return [bridge.plus, bridge.minus]
end
function MOI.get(bridge::PlusMinusBridge{T, F, G, ST},
                 ::MOI.ListOfConstraintIndices{F, ST}) where {T, F, G, ST}
    return [bridge.plus]
end
function MOI.get(bridge::PlusMinusBridge{T, F, G, ST},
                 ::MOI.ListOfConstraintIndices{G, ST}) where {T, F, G, ST}
    return [bridge.minus]
end

# Indices
function MOI.delete(model::MOI.ModelLike, bridge::PlusMinusBridge)
    MOI.delete(model, bridge.plus)
    MOI.delete(model, bridge.minus)
end


text:
#Atom factors for specific atom types
type atomfactor
    R::Float64
    B::Float64
end

#Constant directory with atom information
atomfactorlist = Dict( "C" => atomfactor(0.7, 6.0),
"N" => atomfactor(0.65, 7.0),
"S" => atomfactor(1.0, 16.0),
"O" => atomfactor(0.6, 8.0),
"default" => atomfactor(0.7, 6.0),
"P" => atomfactor(0.65, 15.0),
"H" => atomfactor(0.3, 1.0),
"D" => atomfactor(0.3, 1.0)
)

"""Atom type with position and type"""
type atomType
    index::Int64
    special_type::String
    pos::Array{Float64}
    t::String
    specific_type::String
    bfactor::Float64
    group::String
end

"""Loads the atom coordinates from a pdb file and returns centered list with atomic coordinates"""
function loadPDB(filename::String; center::Bool=true)
    list = atomType[]

    open(filename, "r") do file
        #print(readstring(file))
        for line in readlines(file)
            if length(line) > 4 && line[1:4] == "ATOM"
                idx = parse(Int64,line[9:11])
                special_type = line[13:16]
                x = float(line[31:38])
                y = float(line[39:46])
                z = float(line[47:54])
                bfac = float(line[62:66])
                t = line[78:78]
                specific_type = line[13:16]
                group = line[18:20]
                push!(list,atomType(idx,special_type,Float64[x,y,z],t,specific_type,bfac,group))
            end
        end
    end
    if center
        return centerStructure(list)
    else
        return list
    end
end

"""Centers list of atoms to the center of mass"""
function centerStructure(atomlist)
    center = zeros(3)
    for atom in atomlist
        center = center + atom.pos
    end

    center = center / length(atomlist)

    # shift by center
    for atom in atomlist
        atom.pos = atom.pos - center
    end
    return atomlist
end

"""Calculates density for a given point in spherical coordinates"""
function getDensity_cartesian(atomlist, pos::Array{Float64,1})
    result = 0.0
    for atom in atomlist
        rsq = norm(pos-atom.pos)^2
        if rsq < 5 && atom.t != "1"
            atomtype = atomfactorlist[atom.t]
            B = atomtype.B
            R = atomtype.R
            result += B/((R*sqrt(2.0*pi))^3) *exp(-rsq/(2.0*R^2))
        end
    end
    return result
end

"""Calculates intensity for a given point in spherical coordinates"""
function getIntensity_cartesian(atomlist, pos::Array{Float64,1})

    k2 = norm(pos)^2
    sumc = 0
    sums = 0
    for atom in atomlist
        R = atomfactorlist[atom.t].R
        B = atomfactorlist[atom.t].B
        expk = exp( -0.5 * R*R*k2)
        sumc += B* cos( dot(pos, atom.pos)) * expk
        sums += B* sin( dot(pos, atom.pos)) * expk
    end
    return sumc^2 + sums^2
end

"""Calculates intensity for a given point in spherical coordinates"""
function getFourier_cartesian(atomlist, pos::Array{Float64,1})

    k2 = norm(pos)^2
    sumc = 0
    sums = 0
    for atom in atomlist
        R = atomfactorlist[atom.t].R
        B = atomfactorlist[atom.t].B
        expk = exp( -0.5 * R*R*k2)
        sumc += B* cos( dot(pos, atom.pos)) * expk
        sums += B* sin( dot(pos, atom.pos)) * expk
    end
    return sumc + 1i*sums
end

##################################################################
# Structure handling
##################################################################

"""Loads a pdb file and creates a spherical harmonics description of the electron denisty, the Fourier transform and the Intensity (absolute square FT)"""
function createSphericalHarmonicsStructure(moleculepdb, LMAX::Integer,KMAX::Integer,rmax::Float64)
    println("Initializing SH structure from $moleculepdb with LMAX=$LMAX, KMAX=$KMAX, rmax=$rmax.")
    atomlist = loadPDB(moleculepdb)

    density = getDensity(atomlist, LMAX, KMAX, rmax)
    fourierTransform = forward(density)
    intensity = absoluteSquare(fourierTransform)

    return density, fourierTransform, intensity
end

"""Loads a pdb file and creates a cubic description of the electron denisty, the Fourier transform and the Intensity (absolute square FT)"""
function createCubicStructure(moleculepdb, cubesize::Integer, rmax::Float64)
    println("Initializing cubic structure from $moleculepdb with cubesize=$cubesize, rmax=$rmax.")
    atomlist = loadPDB(moleculepdb)

    #Calculate optional cubic representation
    densCube = getDensityCube(atomlist, cubesize, rmax)
    fourierCube = forward(densCube)
    intensCube = absoluteSquare(fourierCube)

    return densCube, fourierCube, intensCube
end

"""Expands a cubic volume in a spherical harmonics expansion"""
function cubeToSphericalHarmonics(volume::CubeVolume, KMAX::Int64, LMAX::Int64)
    rmax = volume.rmax
    surflist = Array{Complex{Float64},1}[]
    dr = rmax/KMAX
    for k = 1:KMAX
        r = k*dr
        push!(surflist, getSurface((phi,theta)->getVolumeInterpolated(volume, sphericalToCartesian(r,phi,theta)), LMAX))
    end
    surfVolume = SurfaceVolume(surflist, LMAX, KMAX, rmax)

    return getSphericalHarmonicsVolume(surfVolume)

end


text:
using VisualGrowthDynamics
using Test

@testset "VisualGrowthDynamics.jl" begin
    # Write your own tests here.
end


text:
"""
    module IrosBase

Essential tools used in the rest of Iros.
"""
module IrosBase

include("Shapes.jl")
using .Shapes
export Circle,
Ellipse,
DeformedDisk,
Annulus,
Square,
Rectangle,
AbstractShape

include("DielectricFunctions.jl")
using .DielectricFunctions
export DielectricFunction

include("Lattices.jl")
using .Lattices
export Lattices

include("Boundaries.jl")
using .Boundaries
export Boundary,
noBC,
DirichletBC,
NeumannBC,
FloquetBC,
MatchedBC,
PML,
cPML,
noBL

include("Tessellations.jl")
using .Tessellations

include("Domains.jl")
using .Domains
export Domain,
Cavity,
Resonator,
Void,
Dielectric,
Waveguide,
Lead

include("Simulations.jl")
using .Simulations
export Simulation,
smooth!,
update_dielectric!,
update_pump!

include("Laplacians.jl")
using .Laplacians
export Laplacian

end # module


text:
###################################
## Example code
## for Bayesian Poisson Mixture Model

push!(LOAD_PATH,".")
using PyPlot, PyCall
import PoissonMixtureModel

"""
Visualize data & estimation in 2D space.
"""
function visualize_2D(X::Matrix{Float64}, S::Matrix{Float64}, S_est::Matrix{Float64}, text)
    cmp = get_cmap("jet")

    K1 = size(S, 1)
    K2 = size(S_est, 1)
    col1 = [pycall(cmp.o, PyAny, Int(round(val)))[1:3] for val in linspace(0,255,K1)]
    col2 = [pycall(cmp.o, PyAny, Int(round(val)))[1:3] for val in linspace(0,255,K2)]

    f, (ax1, ax2) = subplots(1,2,num=text)
    f[:clf]()
    f, (ax1, ax2) = subplots(1,2,num=text)

    for k in 1 : K1
        ax1[:scatter](X[1, S[k,:].==1], X[2, S[k,:].==1], color=col1[k])
    end
    ax1[:set_title]("truth")

    for k in 1 : K2
        ax2[:scatter](X[1, S_est[k,:].==1], X[2, S_est[k,:].==1], color=col2[k])
    end

    ax2[:set_title]("estimation")
end

function draw_hist(ax, X, S, label)
    counts, bins, patches = ax[:hist](X', 20)
    for i in 1 : length(patches)
        if counts[i] > 0
            S_tmp = S[:,bins[i] .<= X[1,:] .<= bins[i+1]]
            S_sum = sum(S_tmp, 2) / sum(S_tmp)
            patches[i][:set_facecolor]((S_sum[1], 0, S_sum[2]))
        end
    end
    ax[:set_title](label)
end

"""
Visualize data & estimation using 1D histogram.
"""
function visualize_1D(X::Matrix{Float64}, S::Matrix{Float64}, S_est::Matrix{Float64})
    # separated figures
    f1, ax1 = subplots(1,1,num="observation")
    f2, ax2 = subplots(1,1,num="estimation")
    f1[:clf]()
    f2[:clf]()
    _, ax1 = subplots(1,1,num="observation")
    _, ax2 = subplots(1,1,num="estimation")
    ax1[:hist](X', 20)
    ax1[:set_title]("observation")
    draw_hist(ax2, X, S_est, "estimation")
end

"""
Run a test script for 1D data clustering.
"""
function test_1D()
    ## set model
    D = 1 # data dimension, must be 1.
    K = 2 #  number of mixture components, must be 2.
    alpha = 100.0 * ones(K)
    cmp = [PoissonMixtureModel.Gam(1.0*ones(D), 0.01) for i in 1 : K]
    bpmm = PoissonMixtureModel.BPMM(D, K, alpha, cmp)

    ## generate data
    N = 1000
    pmm = PoissonMixtureModel.sample_PMM(bpmm)
    X, S = PoissonMixtureModel.sample_data(pmm, N)

    ## inference
    max_iter = 100
    tic()
    S_est, post_bpmm, VB = PoissonMixtureModel.learn_VI(X, bpmm, max_iter)
    #S_est, post_bpmm, VB = PoissonMixtureModel.learn_GS(X, bpmm, max_iter)
    #S_est, post_bpmm, VB = PoissonMixtureModel.learn_CGS(X, bpmm, max_iter)
    toc()

    ## plot
    visualize_1D(X, S, S_est)

    figure("ELBO")
    clf()
    plot(VB)
    ylabel("ELBO")
    xlabel("iterations")
    show()
end

"""
Run a test script for 2D data clustering.
"""
function test_2D()
    ## set model
    D = 2 # data dimension, must be 2.
    K = 8 # number of mixture components
    #K = 5

    alpha = 100.0 * ones(K)
    cmp = [PoissonMixtureModel.Gam(1.0*ones(D), 0.01) for i in 1 : K]
    bpmm = PoissonMixtureModel.BPMM(D, K, alpha, cmp)

    ## generate data
    N = 300
    pmm = PoissonMixtureModel.sample_PMM(bpmm)
    X, S = PoissonMixtureModel.sample_data(pmm, N)

    ## inference
    max_iter = 100
    tic()
    S_est, post_bpmm, VB = PoissonMixtureModel.learn_VI(X, bpmm, max_iter)
    #S_est, post_bpmm, VB = PoissonMixtureModel.learn_GS(X, bpmm, max_iter)
    #S_est, post_bpmm, VB = PoissonMixtureModel.learn_CGS(X, bpmm, max_iter)
    toc()

    ## plot
    visualize_2D(X, S, PoissonMixtureModel.winner_takes_all(S_est), "2D plot")

    # VB check
    figure("ELBO")
    clf()
    plot(VB)
    ylabel("ELBO")
    xlabel("iterations")
    show()
end

test_1D()
#test_2D()


text:
struct ModelSetup{T<:AbstractFloat,Tprog<:AbstractFloat}
    parameters::Parameter
    grid::Grid{T,Tprog}
    constants::Constants{T,Tprog}
    forcing::Forcing{T}
end