text:
module NetCDFIO

using NCDatasets

float = Union{Float64, Float32}

mutable struct MapInfo{T <: float}
    nx :: Integer
    ny :: Integer
    lsize :: Integer

    xc :: AbstractArray{T, 2}
    yc :: AbstractArray{T, 2}

    mask :: AbstractArray{T, 2}
    area :: AbstractArray{T, 2}
    frac :: AbstractArray{T, 2}

    missing_value :: T

    function MapInfo{T}(
        filename::String;
        missing_value::T = 1e20
    ) where T <: float

        ds = Dataset(filename, "r")

        _mask = ds["mask"][:, :]
        _area = ds["area"][:, :]
        _frac = ds["frac"][:, :]
        _xc  = ds["xc"][:, :]
        _yc  = ds["yc"][:, :]
        _nx  = ds.dim["ni"]
        _ny  = ds.dim["nj"]

        close(ds)

        _mask[_mask .!= 0] .= 1.0

        return new{T}(
            _nx, _ny, _nx * _ny,
            _xc, _yc,
            _mask, _area, _frac,
            missing_value,
        )

    end

end

function createNCFile(
    mi::MapInfo{T},
    filename::String;
) where T <: float

    ds = Dataset(filename, "c")

    defDim(ds, "time", Inf)
    defDim(ds, "ni", mi.nx)
    defDim(ds, "nj", mi.ny)
    defDim(ds, "n", mi.lsize)

    close(ds)

    write2NCFile(mi, filename, "xc", mi.xc;     time_exists=false, missing_value=mi.missing_value)
    write2NCFile(mi, filename, "yc", mi.yc;     time_exists=false, missing_value=mi.missing_value)
    write2NCFile(mi, filename, "mask", mi.mask; time_exists=false, missing_value=mi.missing_value)
    write2NCFile(mi, filename, "frac", mi.frac; time_exists=false, missing_value=mi.missing_value)
    write2NCFile(mi, filename, "area", mi.area; time_exists=false, missing_value=mi.missing_value)

end

"""
    This function serves as multiple variable appending
"""
function write2NCFile(
    mi          :: MapInfo{T},
    filename    :: String,
    vars        :: Dict;
    time        :: Union{Nothing, UnitRange, Integer} = nothing,
    time_exists :: Bool = true,
    missing_value :: Union{T, Nothing} = nothing,
) where T <: float

    for (varname, var) in vars
        #println("varname: ", varname)
        write2NCFile(
            mi, filename, varname, var;
            time=time,
            time_exists=time_exists,
            missing_value=missing_value
        )
    end
end

function write2NCFile(
    mi          :: MapInfo{T},
    filename    :: String,
    varname     :: String,
    var         :: AbstractArray{T};
    time        :: Union{Nothing, UnitRange, Integer} = nothing,
    time_exists :: Bool = true,
    missing_value :: Union{T, Nothing} = nothing,
) where T <: float

    local ds_var

    ds = Dataset(filename, "a")

    # Create variable if it is not in the file yet
    if ! ( varname in keys(ds) )

        ds_var = defVar(ds, varname, T, (time_exists) ? ("ni", "nj", "time") : ("ni", "nj") )

        if missing_value != nothing
            ds_var.attrib["_FillValue"] = missing_value
        end
    else
        ds_var = ds[varname]
    end

    if time_exists
        # Append data if time axis exists

        # If this is a static 2d data, then
        # make the third dimension as time
        if length(size(var)) == 2
            var = view(var, :, :, 1)
        end

        append_time_len = size(var, 3)

        if time == nothing # append at last

            old_time_len = size(ds_var, 3)
            _beg = 1 + old_time_len
            _end = _beg + append_time_len - 1

            ds_var[:, :, _beg:_end] = var

        elseif typeof(time) <: Integer

            _beg = time
            _end = _beg + append_time_len - 1

            ds_var[:, :, _beg:_end] = var

        else  # Assumed UnitRange

            ds_var[:, :, time] = var

        end

    else
        # Simply dump data if time axis does not exist
        # such as coordinate information lat, lon... etc.
        ds_var[:] = var
    end

    close(ds)
end

    #=
    time = defVar(ds, "time", Float64, ("time",))
    time[:] = convert(Array{Float64, 1}, collect(1:time_len))

    v = defVar(ds, varname, eltype(var), ("time",))
    v.attrib["_FillValue"] = missing_value
    v[:] = var
    println("Missing_value:", missing_value)
    =#

function appendNCFile()
end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
## ExoplanetsSysSim/src/transit_observations_unused.jl
## (c) 2015 Eric B. Ford

abstract EphemerisAbstract                     # Abstract, so we can generalize to allow for TTVs.  # TODO OPT: Make sure not a big speed penalty
immutable EphemerisLinear <: EphemerisAbstract  # Currently unused
  P::Float64
  t0::Float64
end

immutable TransitShape                       # Currently Unused
  depth::Float64        # Fractional
  duration::Float64     # days               #  QUERY: What definition of duration here?  Since we're going to use expressions from Price & Rogers for measurement uncertainties, let's treat this as the Full-width, half-max-duration  (until further notice)
  b::Float64            # in stellar radii
end

immutable TransitParameter{EphemT}           # Currently Unusued
  ephem::EphemT
  shape::TransitShape
end

if false  # Warning only works if TransitPlanetObs is mutable.  For testing purposes only
function set!(obs::TransitPlanetObs, P::Float64, t0::Float64, d::Float64, D::Float64)
  obs.period = P
  obs.t0 = t0
  obs.depth = d
  obs.duration = D
end
end

ext: jl


text:
export trajectory

export DFFile

struct DFFile <: MRIFile
  datafilename::String
  trajfilename::String
end

# reading

function trajectory(f::DFFile; kargs...)
  numSamplingPerProfile, numProfiles, nodes = open(f.trajfilename,"r") do fd
    tmp1, numSamplingPerProfile, numProfiles, tmp2 = Int64.(read!(fd,Array{Int32}(undef,4)))
    nodes = read!(fd, Array{Float32}(undef, 2, numSamplingPerProfile,numProfiles))
    return numSamplingPerProfile, numProfiles, nodes
  end

  return Trajectory(reshape(nodes,2,:), numProfiles, numSamplingPerProfile; kargs...)
end

function rawdata(f::DFFile)
  data = open(f.datafilename,"r") do fd
    tmp1, numSamplingPerProfile, numProfiles, numSlices = Int64.(read!(fd,Array{Int32}(undef,4)))
    data = read!(fd, Array{ComplexF32}(undef,numSamplingPerProfile,numProfiles,numSlices))
    return data
  end

  return convert(Array{ComplexF64},data)
end

function rawdata(f::DFFile, slice::Int)
  data = open(f.datafilename,"r") do fd
    # tmp1,numSamplingPerProfile,numProfiles,numSlices= read(fd,Int32,4)
    tmp1,numSamplingPerProfile,numProfiles,numSlices= read!(fd,Array{Int32}(undef,4))
    # data = read(fd,ComplexF32,numSamplingPerProfile,numProfiles,numSlices)
    data = read!(fd, Array{ComplexF32}(undef,numSamplingPerProfile,numProfiles,numSlices))
    return data
  end

  return convert(Array{ComplexF64}, data[:,:,slice])
end


text:
@testset "Read features" begin
    g = path_graph(5)
    features_graph = get_features_graph(g)
    @test features_graph["nb_edge"] == 4
    @test features_graph["nb_vertex"] == 5
    @test features_graph["degree_max"] == 2
    @test features_graph["degree_min"] == 1
    @test features_graph["degree_mean"] == 8 / 5
    @test features_graph["global_clustering_coefficient"] == 0
    @test features_graph["density"] == (2 * 4) / (5  * (5 - 1))
    @test features_graph["radius"] == 2
    @test features_graph["diameter"] == 4

    path = "./test/data/case6ww.m"
    network = PowerFlowNetwork(path, "MATPOWERM")
    features_opf = get_features_opf(network)

    @test (features_opf["PD_max"], features_opf["PD_min"]) == (70, 0)
    @test (features_opf["QD_max"], features_opf["QD_min"]) == (70, 0)
    @test (features_opf["GS_max"], features_opf["GS_min"]) == (0, 0)
    @test (features_opf["BS_max"], features_opf["BS_min"]) == (0, 0)
    @test (features_opf["VM_max"], features_opf["VM_min"]) == (1.07, 1)
    @test (features_opf["VA_max"], features_opf["VA_min"]) == (0, 0)
    @test (features_opf["VMAX_max"], features_opf["VMAX_min"]) == (1.07, 1.05)
    @test (features_opf["VMIN_max"], features_opf["VMIN_min"]) == (1.07, 0.95)

    # gen
    @test (features_opf["PG_max"], features_opf["PG_min"]) == (60, 0)
    @test (features_opf["QG_max"], features_opf["QG_min"]) == (0, 0)
    @test (features_opf["QMAX_max"], features_opf["QMAX_min"]) == (100, 100)
    @test (features_opf["QMIN_max"], features_opf["QMIN_min"]) == (-100, -100)
    @test (features_opf["PMAX_max"], features_opf["PMAX_min"]) == (200, 150)
    @test (features_opf["PMIN_max"], features_opf["PMIN_min"]) == (50, 37.5)

    # branch
    @test (features_opf["BR_R_max"], features_opf["BR_R_min"]) == (0.2, 0.02)
    @test (features_opf["BR_X_max"], features_opf["BR_X_min"]) == (0.4, 0.1)
    @test (features_opf["BR_B_max"], features_opf["BR_B_min"]) == (0.08, 0.02)
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
# variable manipulation
export @zeros, @ones

"""
Create zeros of specific type.

```julia
julia> @i function f(x)
           @zeros Float64 a b c
           # do something
       end
```
"""
macro zeros(T, args...)
    esc(Expr(:block, map(x->:($x ← zero($T)), args)...))
end

macro ones(T, args...)
    esc(Expr(:block, map(x->:($x ← one($T)), args)...))
end

function NiLangCore.chfield(a::AbstractArray, ::typeof(vec), val)
    reshape(val, size(a)...)
end


text:
function generate_problem_dictionary(path_to_parameters_file::String)::Dict{String,Any}

    # initialize -
    problem_dictionary = Dict{String,Any}()

    try

        # load the TOML parameters file -
        toml_dictionary = TOML.parsefile(path_to_parameters_file)

        # --- PUT STUFF INTO problem_dictionary ----
        problem_dictionary["Fructose_6_Phosphate_concentration"] = toml_dictionary["biophysical_constants"]["Fructose_6_Phosphate_concentration"]
        problem_dictionary["Adenosine_Tri_Phosphate_concentration"] = toml_dictionary["biophysical_constants"]["Adenosine_Tri_Phosphate_concentration"]
        problem_dictionary["Phosphofructokinase_concentration"] = toml_dictionary["biophysical_constants"]["Phosphofructokinase_concentration"]
        problem_dictionary["F6P_binding_constant"] = toml_dictionary["biophysical_constants"]["F6P_binding_constant"] * 1000   # convert to uM
        problem_dictionary["ATP_binding_constant"] = toml_dictionary["biophysical_constants"]["ATP_binding_constant"] * 1000   # convert to uM
        problem_dictionary["Catalytic_rate"] = toml_dictionary["biophysical_constants"]["Catalytic_rate"]

        # return -
        return problem_dictionary
    catch error
        throw(error)
    end
end


text:
using Documenter, CompactBases, LinearAlgebra

isdefined(Main, :NOPLOTS) && NOPLOTS || include("plots.jl")

makedocs(;
    modules=[CompactBases],
    format=Documenter.HTML(assets = ["assets/latex.js"],
                           mathengine = Documenter.MathJax()),
    pages = [
        "Home" => "index.md",
        "Overview" => "overview.md",
        "Inner products & norms" => "inner_products.md",
        "Operators" => [
            "Diagonal operators" => "diagonal_operators.md",
        ],
        "Densities" => "densities.md",
        "Finite-differences" => [
            "Overview" => "fd_overview.md",
            "Non-uniform grids" => "fd_non_uniform.md",
        ],
        "B-splines" => [
            "Theory" => "bsplines/theory.md",
            "Usage" => [
                "Basis creation" => "bsplines/usage.md",
                "Knot sets" => "bsplines/knot_sets.md",
                "Splines" => [
                    "Spline creation & evaluation" => "bsplines/splines.md",
                    "Function approximation" => "bsplines/function_approximation.md",
                ],
                "Approximating operators" => "bsplines/operators.md",
                "Examples" => [
                    "Differentiating functions" => "bsplines/differentiation.md",
                    "Ordinary differential equations" => "bsplines/odes.md",
                    "Eigenproblems" => "bsplines/eigenproblems.md"
                ]
            ],
        ],
    ],
    repo="https://github.com/JuliaApproximation/CompactBases.jl/blob/{commit}{path}#L{line}",
    sitename="CompactBases.jl",
    authors="Stefanos Carlström <stefanos.carlstrom@gmail.com>",
    doctest=false
)

deploydocs(;
    repo="github.com/JuliaApproximation/CompactBases.jl",
)


text:
using Profile
using Cthulhu
using ChainPartitioners
using ChainPartitioners: StepHint, Step, Next, Same, Prev
using SparseArrays
using MatrixDepot

function stepstop(f::F, n) where {F}
    @inbounds begin
        f(1, 1, 1)
        for j′ = 2:n
            Step(f)(Same(1), Next(j′), Same(1))
            f(1, j′, 1)
            for j = 2:j′
                Step(f)(Next(j), Same(j′), Same(1))
            end
        end
    end
end

function main()
    work_model = (AffineWorkModel(0, 10, 1), "work_model")
    net_model = (AffineConnectivityModel(0, 10, 1, 100), "net_model")
    sym_model = (AffineMonotonizedSymmetricConnectivityModel(0, 0, 1, 100, 90), "sym_model")
    comm_model = (AffinePrimaryConnectivityModel(0, 10, 1, 0, 100), "comm_model")
    local_model = (AffineSecondaryConnectivityModel(0, 10, 1, 0, 100), "local_model")
    col_block_model = (ColumnBlockComponentCostModel{Int}(3, (w) -> 1 + w), "col_block_model")
    block_model = (BlockComponentCostModel{Int}(1, 3, (1, identity), (1, identity)), "block_model")

    A = SparseMatrixCSC(mdopen("Boeing/ct20stif").A)
    (m, n) = size(A)
    f = oracle_stripe(StepHint(), net_model[1], A, partition_stripe(A', 32, EquiSplitter()))

    stepstop(f, n)
    @profile for i = 1:10 stepstop(f, n) end
    Profile.print()
end

main()


text:
export AsyncScheduler

import Base: show, similar

"""
    AsyncScheduler

`AsyncScheduler` executes scheduled actions asynchronously and uses `Channel` object to order different actions on a single asynchronous task
"""
struct AsyncScheduler{N} <: AbstractScheduler end

Base.show(io::IO, ::AsyncScheduler) = print(io, "AsyncScheduler()")

function AsyncScheduler(size::Int = typemax(Int))
    return AsyncScheduler{size}()
end

Base.similar(::AsyncScheduler{N}) where N = AsyncScheduler{N}()

makeinstance(::Type{D}, ::AsyncScheduler{N}) where { D, N } = AsyncSchedulerInstance{D}(N)

instancetype(::Type{D}, ::Type{<:AsyncScheduler}) where D = AsyncSchedulerInstance{D}

struct AsyncSchedulerDataMessage{D}
    data :: D
end

struct AsyncSchedulerErrorMessage
    err
end

struct AsyncSchedulerCompleteMessage end

const AsyncSchedulerMessage{D} = Union{AsyncSchedulerDataMessage{D}, AsyncSchedulerErrorMessage, AsyncSchedulerCompleteMessage}

mutable struct AsyncSchedulerInstance{D}
    channel        :: Channel{AsyncSchedulerMessage{D}}
    isunsubscribed :: Bool
    subscription   :: Teardown

    AsyncSchedulerInstance{D}(size::Int = typemax(Int)) where D = begin
        return new(Channel{AsyncSchedulerMessage{D}}(size), false, voidTeardown)
    end
end

isunsubscribed(instance::AsyncSchedulerInstance) = instance.isunsubscribed
getchannel(instance::AsyncSchedulerInstance) = instance.channel

function dispose(instance::AsyncSchedulerInstance)
    if !isunsubscribed(instance)
        instance.isunsubscribed = true
        close(instance.channel)
        @async begin
            unsubscribe!(instance.subscription)
        end
    end
end

function __process_channeled_message(instance::AsyncSchedulerInstance{D}, message::AsyncSchedulerDataMessage{D}, actor) where D
    on_next!(actor, message.data)
end

function __process_channeled_message(instance::AsyncSchedulerInstance, message::AsyncSchedulerErrorMessage, actor)
    on_error!(actor, message.err)
    dispose(instance)
end

function __process_channeled_message(instance::AsyncSchedulerInstance, message::AsyncSchedulerCompleteMessage, actor)
    on_complete!(actor)
    dispose(instance)
end

struct AsyncSchedulerSubscription{ H <: AsyncSchedulerInstance } <: Teardown
    instance :: H
end

Base.show(io::IO, ::AsyncSchedulerSubscription) = print(io, "AsyncSchedulerSubscription()")

as_teardown(::Type{ <: AsyncSchedulerSubscription}) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::AsyncSchedulerSubscription)
    dispose(subscription.instance)
    return nothing
end

function scheduled_subscription!(source, actor, instance::AsyncSchedulerInstance)
    subscription = AsyncSchedulerSubscription(instance)

    channeling_task = @async begin
        while !isunsubscribed(instance)
            message = take!(getchannel(instance))
            if !isunsubscribed(instance)
                __process_channeled_message(instance, message, actor)
            end
        end
    end

    subscription_task = @async begin
        if !isunsubscribed(instance)
            tmp = on_subscribe!(source, actor, instance)
            if !isunsubscribed(instance)
                subscription.instance.subscription = tmp
            else
                unsubscribe!(tmp)
            end
        end
    end

    bind(getchannel(instance), channeling_task)

    return subscription
end

function scheduled_next!(actor, value::D, instance::AsyncSchedulerInstance{D}) where { D }
    put!(getchannel(instance), AsyncSchedulerDataMessage{D}(value))
end

function scheduled_error!(actor, err, instance::AsyncSchedulerInstance)
    put!(getchannel(instance), AsyncSchedulerErrorMessage(err))
end

function scheduled_complete!(actor, instance::AsyncSchedulerInstance)
    put!(getchannel(instance), AsyncSchedulerCompleteMessage())
end


text:
module TestProcesses
include("preamble.jl")

@testset "mapfoldl" begin
    # https://clojure.org/reference/transducers#_transduce
    xf = Filter(isodd) |> Map(inc)
    @testset for xs in iterator_variants(0:4)
        @test mapfoldl(xf, +, 0, xs) == 6
        @test mapfoldl(xf, +, 100, xs) == 106
    end

    # https://clojuredocs.org/clojure.core/transduce
    xf = Filter(isodd) |> Take(10)
    @testset "$(typeof(xs))" for xs in iterator_variants(0:1000)
        @test mapfoldl(xf, push!, Int[], xs) == 1:2:19
        @test mapfoldl(xf, +, 0, xs) == 100
        @test mapfoldl(xf, +, 17, xs) == 117
        @test mapfoldl(xf, string, "", xs) == "135791113151719"

        @test transduce(xf, push!, Int[], xs) == Reduced(1:2:19)
        @test transduce(xf, +, 0, xs) == Reduced(100)
        @test transduce(xf, +, 17, xs) == Reduced(117)
        @test transduce(xf, string, "", xs) == Reduced("135791113151719")
    end
end

@testset "eduction" begin
    # https://clojuredocs.org/clojure.core/eduction
    xf = Filter(isodd) |> Take(5)
    @testset "$(typeof(xs))" for xs in iterator_variants(0:1000)
        @test collect(xf, xs) == 1:2:9
        @test collect(eduction(xf, xs)) == 1:2:9
        @test collect(Map(x -> 2x), eduction(xf, xs)) == 2(1:2:9)
    end

    @testset for xs in iterator_variants(1:5)
        @test collect(xf, xs) == 1:2:5
        @test collect(eduction(xf, xs)) == 1:2:5
        @test collect(Map(x -> 2x), eduction(xf, xs)) == 2(1:2:5)
    end

    ed = eduction(xf, 1:5)
    @test finaltype(ed.rf) === Int
end

@testset "append!" begin
    # https://clojuredocs.org/clojure.core/into#example-57294b20e4b050526f331420
    xf = Map(x -> x + 2) |> Filter(isodd)
    @testset for xs in iterator_variants(0:9)
        @test append!(xf, [-1, -2], xs) == [-1, -2, 3, 5, 7, 9, 11]
    end
end

@testset "map!" begin
    xs0 = 1:10
    @testset for src in [xs0, collect(xs0)]
        dest = similar(xs0)
        @test map!(Scan(+), dest, src) === dest
        @test dest == cumsum(xs0)
    end

    xs0 = 1:5
    @testset for src in [xs0, collect(xs0)]
        dest = zero(src)
        @test map!(Filter(isodd) |> Scan(+), dest, src) == [1, 0, 4, 0, 9]
    end
end

@testset "copy!" begin
    xs0 = 1:10
    @testset for src in [xs0, collect(xs0)]
    # @testset for src in iterator_variants(xs0)  # TODO: fix
        dest = similar(xs0)
        @test copy!(Scan(+), dest, src) === dest
        @test dest == cumsum(src)
    end

    xs0 = 1:5
    @testset for src in [xs0, collect(xs0)]
    # @testset for src in iterator_variants(xs0)  # TODO: fix
        dest = zero(xs0)
        @test copy!(Filter(isodd) |> Scan(+), dest, src) == [1, 4, 9]
    end
end

@testset "simple_transduce" begin
    xf = PartitionBy(x -> x > 0) |>
        Filter(xs -> mean(abs, xs) < 1.0) |>
        Map(prod)

    for _ in 1:100
        xs = randn(100)
        @test transduce(xf, +, 0.0, xs) == simple_transduce(xf, +, 0.0, xs)
    end
end

end  # module


text:
using BitTwiddlingConvenienceFunctions
using Documenter

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

makedocs(;
    modules=[BitTwiddlingConvenienceFunctions],
    authors="Julia Computing",
    repo="https://github.com/JuliaSIMD/BitTwiddlingConvenienceFunctions.jl/blob/{commit}{path}#{line}",
    sitename="BitTwiddlingConvenienceFunctions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaSIMD.github.io/BitTwiddlingConvenienceFunctions.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaSIMD/BitTwiddlingConvenienceFunctions.jl",
    devbranch="main",
)


text:
abstract AbstractPolicyTrait

immutable MHRunnerType <: AbstractPolicyTrait
    trait::Symbol
    MHRunnerType(s::Symbol) = (@assert in(s,[:standard,:generalized]) ; new(s))
end

immutable InitializeFrom <: AbstractPolicyTrait
    trait::Symbol
    InitializeFrom(s::Symbol) = (@assert in(s,[:default,:prior]) ; new(s))
end

immutable ProposeFrom <: AbstractPolicyTrait
    trait::Symbol
    ProposeFrom(s::Symbol) = (@assert in(s,[:indicator,:auxiliary]) ; new(s))
end

immutable IndicatorType <: AbstractPolicyTrait
    trait::Symbol
    IndicatorType(s::Symbol) = (@assert in(s,[:stationary,:cyclical]) ; new(s))
end

immutable JobSegments <: AbstractPolicyTrait
    trait::Symbol
    JobSegments(s::Symbol) = (@assert in(s,[:procs,:workers,:test,:none]) ; new(s))
end

immutable ChainType <: AbstractPolicyTrait
    trait::Symbol
    ChainType(s::Symbol) = (@assert in(s,[:standard,:gradient]) ; new(s))
end

immutable StoreDuring <: AbstractPolicyTrait
    trait::Symbol
    StoreDuring(s::Symbol) = (@assert in(s,[:burnin,:main,:all]) ; new(s))
end

@inline _traitname(::Type{MHRunnerType}) = :mhrunner
@inline _traitname(::Type{InitializeFrom}) = :initialize
@inline _traitname(::Type{ProposeFrom}) = :propose
@inline _traitname(::Type{IndicatorType}) = :indicator
@inline _traitname(::Type{JobSegments}) = :jobsegments
@inline _traitname(::Type{ChainType}) = :chain
@inline _traitname(::Type{StoreDuring}) = :store
@inline _traitnametype(t) = Val{_traitname(t)}

trait(n::Symbol,t::Symbol) = _trait(Val{n},t)
_trait(::Type{_traitnametype(MHRunnerType)},t::Symbol)  = MHRunnerType(t)
_trait(::Type{_traitnametype(InitializeFrom)},t::Symbol) = InitializeFrom(t)
_trait(::Type{_traitnametype(ProposeFrom)},t::Symbol) = ProposeFrom(t)
_trait(::Type{_traitnametype(IndicatorType)},t::Symbol) = IndicatorType(t)
_trait(::Type{_traitnametype(JobSegments)},t::Symbol) = JobSegments(t)
_trait(::Type{_traitnametype(ChainType)},t::Symbol) = ChainType(t)
_trait(::Type{_traitnametype(StoreDuring)},t::Symbol) = StoreDuring(t)

traitvalue(a::AbstractPolicyTrait) = a.trait
traittype(a::AbstractPolicyTrait) = Val{a.trait}


text:
# Demultiplexer
# =============
#
# Sequence demultiplexer based on DNA barcodes.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# Trie-like data structure for DNA barcodes.
immutable BarcodeTrie
    nodes::Vector{Int64}
end

function BarcodeTrie(barcodes, ids)
    @assert length(barcodes) == length(ids)
    if !issorted(barcodes)
        error("DNA barcodes must be sorted")
    end
    for i in 1:endof(barcodes)-1
        if barcodes[i] == barcodes[i+1]
            error("duplicated sequences")
        end
    end
    nodes = Int64[1]
    build_trie!(nodes, 1, barcodes, ids, 1:endof(barcodes), 1)
    return BarcodeTrie(nodes)
end

# Build a trie recursively.
function build_trie!(nodes, root, barcodes, ids, range, j)
    n = length(range)
    if n == 0
        return
    elseif n == 1
        i = first(range)
        if j > length(barcodes[i]) + 1
            # leaf
            nodes[root] = ids[i] << 1
            return
        end
    end

    # allocate six consecutive nodes {'#', 'A', 'C', 'G', 'T', 'N' (or others)}
    resize!(nodes, length(nodes) + 6)
    for l in endof(nodes)-5:endof(nodes)
        nodes[l] = 0
    end
    base = endof(nodes) - 6

    # calculate ranges for the next level
    starts = Vector{Int}(6)
    i = first(range)
    starts[1] = i  # '#'
    while i ≤ last(range) && length(barcodes[i]) < j
        i += 1
    end
    for k in 1:4  # ACGT
        starts[k+1] = i
        while i ≤ last(range) && findfirst(ACGT, barcodes[i][j]) == k
            i += 1
        end
    end
    starts[6] = i  # N

    # set nodes
    nodes[root] = (base << 1) | 1
    for k in 1:6
        if k == 6
            r = starts[k]:last(range)
        else
            r = starts[k]:starts[k+1]-1
        end
        build_trie!(nodes, base + k, barcodes, ids, r, j + 1)
    end
end

# Find the longest match of `seq` prefixes in `trie`.
function findbarcode(trie::BarcodeTrie, seq)
    s = 1
    for nt in seq
        base = trie.nodes[s] >> 1
        if nt == DNA_A
            t = base + 2
        elseif nt == DNA_C
            t = base + 3
        elseif nt == DNA_G
            t = base + 4
        elseif nt == DNA_T
            t = base + 5
        else
            t = base + 6
        end
        if trie.nodes[t] & 1 == 0
            # leaf
            break
        end
        s = t
    end
    s = (trie.nodes[s] >> 1) + 1
    return trie.nodes[s] >> 1
end

type Demultiplexer{S}
    # barcode tries with ≥ 0 errors in ascending order
    # (i.e. tries[1]: no errors, tries[2]: one error, ...)
    tries::Vector{BarcodeTrie}

    # original barcodes
    barcodes::Vector{S}

    # distance (:hamming or :levenshtein)
    distance::Symbol
end

function Base.show(io::IO, demultiplexer::Demultiplexer)
    println(io, summary(demultiplexer), ":")
    println(io, "  distance: ", demultiplexer.distance)
    println(io, "  number of barcodes: ", length(demultiplexer.barcodes))
      print(io, "  number of correctable errors: ", length(demultiplexer.tries) - 1)
end

"""
    Demultiplexer(barcodes::Vector{DNASequence};
                  n_max_errors::Integer=1,
                  distance::Symbol=:hamming)

Create a demultiplexer object from `barcodes`.

# Arguments
* `barcodes`: the sorted list of DNA barcodes; integer IDs are assigned in this order.
* `n_max_errors=1`: the number of maximum correctable errors.
* `distance=:hamming`: the distance metric (`:hamming` or `:levenshtein`).
"""
function Demultiplexer(barcodes::Vector{DNASequence};
                       n_max_errors::Integer=1,
                       distance::Symbol=:hamming)
    if n_max_errors < 0
        error("n_max_errors must be non-negative")
    elseif distance ∉ (:hamming, :levenshtein)
        error("distance must be either :hamming or :levenshtein")
    end

    # check barcodes
    for barcode in barcodes
        if hasambiguity(barcode)
            error("barcode must be umambiguous")
        end
    end

    # build trie trees for each error level (0-n_max_errors)
    tries = BarcodeTrie[]
    for m in 0:n_max_errors
        # generate "erroneous" barcodes
        mutated_barcodes = DNASequence[]
        ids = Int[]
        for (i, barcode) in enumerate(barcodes)
            if distance == :hamming
                circle = hamming_circle(barcode, m)
            else
                @assert distance == :levenshtein
                circle = levenshtein_circle(barcode, m)
            end
            append!(mutated_barcodes, circle)
            append!(ids, collect(repeated(i, length(circle))))
        end
        ord = sortperm(mutated_barcodes)
        push!(tries, BarcodeTrie(mutated_barcodes[ord], ids[ord]))
    end

    return Demultiplexer(tries, barcodes, distance)
end

"""
    demultiplex(demultiplexer::Demultiplexer,
                seq::Sequence,
                linear_search_fallback::Bool=false) -> (index, distance)

Return a barcode index that matches `seq` with least errors and its distance.

The upper limit of the number of maximum errors is bounded by the `n_max_errors`
parameter of `demultiplexer`. When `linear_search_fallback` is `true`, this
function tries to find the best matching barcodes using linear search and
returns one of them at random.
"""
function demultiplex(demultiplexer::Demultiplexer, seq::Sequence, linear_search_fallback::Bool=false)
    if eltype(seq) != DNANucleotide
        error("sequence must be a DNA sequence")
    end
    for (k, trie) in enumerate(demultiplexer.tries)
        i = findbarcode(trie, seq)
        if i != 0
            return i, k - 1
        end
    end
    if !linear_search_fallback
        # not found
        return 0, -1
    end

    # find the best matching barcode using linear search
    i_min = Int[]
    dist_min = typemax(Int)
    for (i, barcode) in enumerate(demultiplexer.barcodes)
        if demultiplexer.distance == :hamming
            dist = hamming_distance(barcode, seq[1:length(barcode)])
        elseif demultiplexer.distance == :levenshtein
            dist = sequencelevenshtein_distance(barcode, seq)
        else
            assert(false)
        end
        if dist ≤ dist_min
            if dist < dist_min
                empty!(i_min)
            end
            push!(i_min, i)
            dist_min = dist
        end
    end
    @assert !isempty(i_min)

    return rand(i_min), dist_min
end

function Base.getindex(demultiplexer::Demultiplexer, i::Integer)
    return demultiplexer.barcodes[i]
end

# Generate a list of sequences s.t. `hamming_distance(seq, seq′) == m`.
function hamming_circle(seq, m)
    if m == 0
        return [seq]
    end
    ret = DNASequence[]
    for ps in combinations(1:endof(seq), m)
        for rs in product(repeated(1:4, m)...)
            seq′ = copy(seq)
            for (p, r) in zip(ps, rs)
                if findfirst(ACGT, seq[p]) ≤ r
                    r += 1
                end
                seq′[p] = ACGTN[r]
            end
            push!(ret, seq′)
        end
    end
    return ret
end

# Calculate the Hamming distance between `seq1` and `seq2`.
function hamming_distance(seq1, seq2)
    @assert length(seq1) == length(seq2)
    n = 0
    for (x, y) in zip(seq1, seq2)
        n += x != y
    end
    return n
end

# Generate a list of sequences s.t. `levenshtein_distance(seq, seq′) == m`.
function levenshtein_circle(seq, m)
    if m == 0
        return [seq]
    end
    @assert m > 0
    seqs = DNASequence[]
    # substitution
    append!(seqs, hamming_circle(seq, 1))
    # deletion
    for i in 1:endof(seq)
        push!(seqs, deleteat!(copy(seq), i))
    end
    # insertion
    for i in 1:endof(seq), nt in ACGTN
        if nt != seq[i]
            push!(seqs, insert!(copy(seq), i, nt))
        end
    end
    ball = copy(seqs)
    for seq′ in unique(seqs)
        append!(ball, levenshtein_circle(seq′, m - 1))
    end
    ball = sort!(unique(ball))
    ret = DNASequence[]
    for seq′ in ball
        if levenshtein_distance(seq, seq′) == m
            push!(ret, seq′)
        end
    end
    return ret
end

# Calculate the Levenshtein distance between `seq1` and `seq2`.
function levenshtein_distance(seq1, seq2)
    m = length(seq1)
    n = length(seq2)
    dist = Matrix{Int}(m + 1, n + 1)

    dist[1,1] = 0
    for i in 1:m
        dist[i+1,1] = i
    end

    for j in 1:n
        dist[1,j+1] = j
        for i in 1:m
            dist[i+1,j+1] = min(
                dist[i,j+1] + 1,
                dist[i+1,j] + 1,
                dist[i,j]   + ifelse(seq1[i] == seq2[j], 0, 1))
        end
    end

    return dist[m+1,n+1]
end

# Calculate the Sequence-Levenshtein distance between `seq1` and `seq2`.
# Buschmann, Tilo, and Leonid V. Bystrykh. "Levenshtein error-correcting
# barcodes for multiplexed DNA sequencing." BMC bioinformatics 14.1 (2013): 272.
function sequencelevenshtein_distance(seq1, seq2)
    m = length(seq1)
    n = length(seq2)
    dist = Matrix{Int}(m + 1, n + 1)

    dist[1,1] = 0
    for i in 1:m
        dist[i+1,1] = i
    end

    for j in 1:n
        dist[1,j+1] = j
        for i in 1:m
            dist[i+1,j+1] = min(
                dist[i,j+1] + 1,
                dist[i+1,j] + 1,
                dist[i,j]   + ifelse(seq1[i] == seq2[j], 0, 1))
        end
    end

    mindist = dist[end,end]
    for i in 0:m
        mindist = min(mindist, dist[i+1,end])
    end
    for j in 0:n
        mindist = min(mindist, dist[end,j+1])
    end

    return mindist
end


text:
struct Grid

    r0 :: Vector{Float64}
    r1 :: Vector{Float64}
    points :: Matrix{Float64}
    Proj :: SparseMatrixCSC
    dx :: Float64
    dy :: Float64
    Nk :: Int64
    Nz :: Int64
    _Nx :: Int64
    _Ny :: Int64
    _Nbuffer :: Int64
    _nan_outside_boundaries :: Matrix{Float64}

end

# TODO: parallelize this calculation
function Grid(is_inside::Function, r0::Vector{Float64}, r1::Vector{Float64}, Nx::Int64, Ny::Int64, Nz::Int64; Nbuffer=100)

    deltaX = r1[1] - r0[1]
    deltaY = r1[2] - r0[2]
    dx = deltaX / Nx
    dy = deltaY / Ny

    _Nx = Nx + 2*Nbuffer
    _Ny = Ny + 2*Nbuffer

    points = zeros(Float64, (2, Nx*Ny))
    proj_rows = Int32[i for i=1:Nx*Ny]
    proj_cols = zeros(Int32, Nx*Ny)
    proj_vals = ones(Int32, Nx*Ny)
    _nan_outside_boundaries = fill(NaN, (Nx, Ny))
    k = 0

    for j=1:Ny
        for i=1:Nx

            x = r0[1] + (i-1)*dx
            y = r0[2] + (j-1)*dy
            if is_inside(x, y)
                k += 1
                points[:,k] = [x, y]
                proj_cols[k] = i + _Nx*(j-1+Nbuffer) + Nbuffer
                _nan_outside_boundaries[i,j] = 1.0
            end

        end
    end

    Proj = sparse(proj_rows[1:k], proj_cols[1:k], proj_vals[1:k], k, _Nx*_Ny)
    return Grid(r0, r1, points[:,1:k], Proj, dx, dy, k, Nz, _Nx, _Ny, Nbuffer, _nan_outside_boundaries)
end

# constructor for cases where the "do-block" syntax might make the code more readable
Grid(f::Function, Nx, Ny, Nz; Nbuffer=100) = Grid(f()..., Nx, Ny, Nz; Nbuffer=Nbuffer)

@recipe function f(v::Vector{Float64}, grd::Grid, z::Int64)

    vp = reshape(v, grd.Nk, :)[:,z]
    vp_cart = transpose(grd.Proj) * vp
    V = reshape(vp_cart, grd._Nx, grd._Ny)
    Vplot = transpose(V[grd._Nbuffer+1:grd._Nx-grd._Nbuffer, grd._Nbuffer+1:grd._Ny-grd._Nbuffer] .* grd._nan_outside_boundaries)

    Nx = grd._Nx - 2*grd._Nbuffer
    Ny = grd._Ny - 2*grd._Nbuffer
    x = LinRange(grd.r0[1], grd.r1[1], Nx)
    y = LinRange(grd.r0[2], grd.r1[2], Ny)

    xlims --> (x[1], x[end])
    ylims --> (y[1], y[end])
    size --> (1000, 1000)
    aspect_ratio --> :equal

    x, y, Vplot
end

function f_to_grid(f::Function, grd::Grid)

    vec = Float64[f(grd.points[:,k]...) for k=1:grd.Nk]

    return vcat([vec for _=1:grd.Nz]...)
end


text:
module Google_53

using Printf
using XEB
using OMEinsum, OMEinsumContractionOrders
using CUDA
using Graphs
using Statistics

L = sort!([52, 37, 35, 32, 22, 11, 31, 21, 8, 24, 7, 1,
    29, 18, 5, 26, 15, 6, 40, 25, 16, 44, 42, 51, 53, 48, 46])
R = setdiff(1:53, L)

D = 20
layout = google_layout_53(53, D)
cuts = XEB.generate_cut(layout, L, 1, D)
XEB.simplify!(layout, cuts);
@show flux(g)
flux_g = flux(g)
connected_components(layout)
ec_L, ts_L, ids_size_L = to_ein_code(layout, L; enable_cuda = true)
ec_R, ts_R, ids_size_R = to_ein_code(layout, R; enable_cuda = true)

# ec_L_greedy = optimize_code(ec_L, ids_size_L, GreedyMethod())
ec_L_greedy = optimize_greedy(ec_L, ids_size_L; nrepeat = 100)
@show timespace_complexity(ec_L_greedy, ids_size_L)
ec_L_opt = ec_L_greedy
# ec_R_greedy = optimize_code(ec_R, ids_size_R, GreedyMethod())
ec_R_greedy = optimize_greedy(ec_R, ids_size_R; nrepeat = 100)
@show timespace_complexity(ec_R_greedy, ids_size_R)
ec_R_opt = ec_R_greedy

# ec_R_opt = optimize_code(ec_R_greedy, ids_size_R, TreeSA())
# timespace_complexity(ec_R_opt, ids_size_R)
CUDA.@time ec_L_opt(ts_L...);
CUDA.@time ec_R_opt(ts_R...);
xebs_L = Float64[]
xebs_R = Float64[]
maxs_L = Float64[]
maxs_R = Float64[]

for i = 1:10000
    g = google_layout_53(53, 20);
    cuts = XEB.generate_cut(g, L, 1, 20)
    XEB.simplify!(g, cuts);
    _, ts_L, _ = to_ein_code(g, L; haar = false, enable_cuda = true)
    _, ts_R, _ = to_ein_code(g, R; haar = false, enable_cuda = true)
    p_L = ec_L_opt(ts_L...);
    p_R = ec_R_opt(ts_R...);
    # @show sum(p_L), sum(p_R)
    # p_L = p_L / sum(p_L);
    # p_R = p_R / sum(p_R);
    xeb_L = 2^(length(ec_L.iy))*sum(p_L .* p_L) - 1
    push!(xebs_L, xeb_L)
    xeb_R = 2^(length(ec_R.iy))*sum(p_R .* p_R) - 1
    push!(xebs_R, xeb_R)
    push!(maxs_L, (2^(length(ec_L.iy))*maximum(p_L) - 1))
    push!(maxs_R, (2^(length(ec_R.iy))*maximum(p_R) - 1))
    @printf "XEB_max = %.3e = %.3e + %.3e (flux = %.3e)\n" (mean(maxs_L)+mean(maxs_R)) mean(maxs_L) mean(maxs_R) 2.0^(-flux_g)
    # println("$(i): max_L = $(mean(maxs_L)), max_R = $(mean(maxs_R)), max = $(mean(maxs_R)+mean(maxs_L))")
    # println("$(i): XEB_L = $(mean(xebs_L)), XEB_R = $(mean(xebs_R))")
end

std(xebs_L) / 100
std(xebs_R) / 100

end


text:
"""
    HyperSphere{N, T}

A `HyperSphere` is a generalization of a sphere into N-dimensions.
A `center` and radius, `r`, must be specified.
"""
struct HyperSphere{N,T} <: GeometryPrimitive{N,T}
    center::Point{N,T}
    r::T
end

"""
    Circle{T}

An alias for a HyperSphere of dimension 2. (i.e. `HyperSphere{2, T}`)
"""
const Circle{T} = HyperSphere{2,T}

"""
    Sphere{T}

An alias for a HyperSphere of dimension 3. (i.e. `HyperSphere{3, T}`)
"""
const Sphere{T} = HyperSphere{3,T}

HyperSphere{N}(p::Point{N,T}, number) where {N,T} = HyperSphere{N,T}(p, convert(T, number))

widths(c::HyperSphere{N,T}) where {N,T} = Vec{N,T}(radius(c) * 2)
radius(c::HyperSphere) = c.r
origin(c::HyperSphere) = c.center

Base.minimum(c::HyperSphere{N,T}) where {N,T} = Vec{N,T}(origin(c)) - Vec{N,T}(radius(c))
Base.maximum(c::HyperSphere{N,T}) where {N,T} = Vec{N,T}(origin(c)) + Vec{N,T}(radius(c))

function Base.in(x::AbstractPoint, c::HyperSphere)
    return norm(origin(c) - x) ≤ radius(c)
end

centered(S::Type{HyperSphere{N,T}}) where {N,T} = S(Vec{N,T}(0), T(0.5))
function centered(::Type{T}) where {T<:HyperSphere}
    return centered(HyperSphere{ndims_or(T, 3),eltype_or(T, Float32)})
end

function coordinates(s::Circle, nvertices=64)
    rad = radius(s)
    inner(fi) = Point(rad * sin(fi + pi), rad * cos(fi + pi)) .+ origin(s)
    return (inner(fi) for fi in LinRange(0, 2pi, nvertices))
end

function texturecoordinates(s::Circle, nvertices=64)
    return coordinates(Circle(Point2f(0.5), 0.5f0), nvertices)
end

function coordinates(s::Sphere, nvertices=24)
    θ = LinRange(0, pi, nvertices)
    φ = LinRange(0, 2pi, nvertices)
    inner(θ, φ) = Point(cos(φ) * sin(θ), sin(φ) * sin(θ), cos(θ)) .* s.r .+ s.center
    return ivec((inner(θ, φ) for θ in θ, φ in φ))
end

function texturecoordinates(s::Sphere, nvertices=24)
    ux = LinRange(0, 1, nvertices)
    return ivec(((φ, θ) for θ in reverse(ux), φ in ux))
end

function faces(sphere::Sphere, nvertices=24)
    return faces(Rect(0, 0, 1, 1), (nvertices, nvertices))
end

function normals(s::Sphere{T}, nvertices=24) where {T}
    return coordinates(Sphere(Point{3,T}(0), 1), nvertices)
end


text:
reshapeVec2Mat(vec::Vector, rows::Int) = reshape(vec, rows, round(Int,length(vec)/rows))
# function reshapeVec2Mat(vec::Vector, rows::Int)
#   M = reshape(vec, rows, round(Int,length(vec)/rows))
#   return ndims(M) < 2 ? (M')' : M
# end

# get vertex from factor graph according to label symbol "x1"
getVert(fgl::FactorGraph, lbl::Symbol; api::DataLayerAPI=dlapi, nt::Symbol=:var) = api.getvertex(fgl, lbl, nt=nt)
getVert(fgl::FactorGraph, id::Int; api::DataLayerAPI=dlapi) = api.getvertex(fgl, id)

# TODO -- upgrade to dedicated memory location in Graphs.jl
# see JuliaArchive/Graphs.jl#233
getData(v::Graphs.ExVertex) = v.attributes["data"]
# Convenience functions
getData(fgl::FactorGraph, lbl::Symbol; api::DataLayerAPI=dlapi) = getData(getVert(fgl, lbl, api=api))
getData(fgl::FactorGraph, id::Int; api::DataLayerAPI=dlapi) = getData(getVert(fgl, id, api=api))

function setData!(v::Graphs.ExVertex, data)
  v.attributes["data"] = data
  nothing
end

function getVal(v::Graphs.ExVertex)
  return getData(v).val
end
function getVal(v::Graphs.ExVertex, idx::Int)
  return getData(v).val[:,idx]
end

# Convenience function to get values for given variable label
function getVal(fgl::FactorGraph, lbl::Symbol; api::DataLayerAPI=dlapi)
  #getVal(dlapi.getvertex(fgl, lbl))
  getVal(getVert(fgl, lbl, api=api))
end
function getVal(fgl::FactorGraph, exvertid::Int; api::DataLayerAPI=dlapi)
  # getVal(dlapi.getvertex(fgl, exvertid))
  getVal(getVert(fgl, exvertid, api=api))
end

function getNumPts(v::Graphs.ExVertex)
  return size(getData(v).val,2)
end

function getfnctype(data)
  if typeof(data).name.name == :VariableNodeData
    return VariableNodeData
  end
  data.fnc.usrfnc!
end

function getfnctype(vertl::Graphs.ExVertex)
  data = getData(vertl)
  getfnctype(data)
end

function getfnctype(fgl::FactorGraph, exvertid::Int; api::DataLayerAPI=dlapi)
  #
  # data = getData(fgl, exvertid, api=api)
  # data.fnc.usrfnc!
  getfnctype(getVert(fgl, exvertid, api=api))
end

# setVal! assumes you will update values to database separate, this used for local graph mods only
function setVal!(v::Graphs.ExVertex, val::Array{Float64,2})
  getData(v).val = val
  # v.attributes["data"].val = val
  nothing
end
function getBWVal(v::Graphs.ExVertex)
  return getData(v).bw
end
function setBW!(v::Graphs.ExVertex, bw::Array{Float64,2})
  getData(v).bw = bw
  # v.attributes["data"].bw = bw
  nothing
end
function setVal!(v::Graphs.ExVertex, val::Array{Float64,2}, bw::Array{Float64,2})
  setVal!(v,val)
  setBW!(v,bw)
  nothing
end
function setVal!(v::Graphs.ExVertex, val::Array{Float64,2}, bw::Vector{Float64})
  setVal!(v,val,reshape(bw,length(bw),1)) #(bw')')
  nothing
end
function setVal!(fg::FactorGraph, sym::Symbol, val::Array{Float64,2}; api::DataLayerAPI=localapi)
  setVal!(api.getvertex(fg, sym), val)
end
function setValKDE!(v::Graphs.ExVertex, val::Array{Float64,2})
  p = kde!(val)
  setVal!(v,val,getBW(p)[:,1]) # TODO -- this can be little faster
  nothing
end
function setValKDE!(v::Graphs.ExVertex, em::EasyMessage)
  setVal!(v, em.pts, em.bws ) # getBW(p)[:,1]
  nothing
end
function setValKDE!(v::Graphs.ExVertex, p::BallTreeDensity)
  pts = getPoints(p)
  setVal!(v, pts, getBW(p)[:,1]) # BUG ...al!(., val, . ) ## TODO -- this can be little faster
  nothing
end
setVal!(v::Graphs.ExVertex, em::EasyMessage) = setValKDE!(v, em)
setVal!(v::Graphs.ExVertex, p::BallTreeDensity) = setValKDE!(v, p)

function kde!(em::EasyMessage)
  return kde!(em.pts,em.bws)
end

# TODO -- there should be a better way, without retrieving full vertex
getOutNeighbors(fgl::FactorGraph, v::ExVertex; api::DataLayerAPI=dlapi, needdata::Bool=false, ready::Int=1,backendset::Int=1 ) = api.outneighbors(fgl, v, needdata=needdata, ready=ready, backendset=backendset )
getOutNeighbors(fgl::FactorGraph, vertid::Int; api::DataLayerAPI=dlapi, needdata::Bool=false, ready::Int=1,backendset::Int=1 ) = api.outneighbors(fgl, api.getvertex(fgl,vertid), needdata=needdata, ready=ready, backendset=backendset )

function updateFullVert!(fgl::FactorGraph, exvert::ExVertex;
            api::DataLayerAPI=IncrementalInference.dlapi,
            updateMAPest::Bool=false  )
  #
  warn("use of updateFullVert! should be clarified for local or remote operations.")
  api.updatevertex!(fgl, exvert, updateMAPest=updateMAPest)
end

function setDefaultNodeData!(v::Graphs.ExVertex,
                             initval::Array{Float64,2},
                             stdev::Array{Float64,2},
                             dodims::Int,
                             N::Int,
                             dims::Int;
                             gt=nothing, initialized::Bool=true,
                             softtype=nothing)
  # TODO review and refactor this function, exists as legacy from pre-v0.3.0
  # this should be the only function allocating memory for the node points (unless number of points are changed)
  data = nothing
  if initialized
    if size(initval,2) < N && size(initval, 1) == dims
      warn("setDefaultNodeData! -- deprecated use of stdev.")
      p = kde!(initval,diag(stdev));
      pN = resample(p,N)
    elseif size(initval,2) < N && size(initval, 1) != dims
      println("Node value memory allocated but not initialized")
      pN = kde!(randn(dims, N));
    else
      pN = kde!(initval)
    end
    # dims = size(initval,1) # rows indicate dimensions
    sp = Int[0;] #round.(Int,linspace(dodims,dodims+dims-1,dims))
    gbw = getBW(pN)[:,1]
    gbw2 = Array{Float64}(length(gbw),1)
    gbw2[:,1] = gbw[:]
    pNpts = getPoints(pN)
    data = VariableNodeData(initval, stdev, pNpts,
                            gbw2, Int[], sp,
                            dims, false, 0, Int[], gt, softtype, true) #initialized
  else
      sp = round.(Int,linspace(dodims,dodims+dims-1,dims))
      data = VariableNodeData(initval, stdev, zeros(dims, N),
                              zeros(dims,1), Int[], sp,
                              dims, false, 0, Int[], gt, softtype, false) #initialized
  end
  #
  setData!(v, data)
  # v.attributes["data"] = data

  # what about dedicated user variable metadata
  # p.factormetadata

  nothing
end

"""
    $(SIGNATURES)

Initialize a new Graphs.ExVertex which will be added to some factor graph.
"""
function addNewVarVertInGraph!(fgl::FactorGraph,
            vert::Graphs.ExVertex,
            id::Int,
            lbl::Symbol,
            ready::Int,
            smalldata  )
  #
  vert.attributes = Graphs.AttributeDict() #fg.v[fg.id]
  vert.attributes["label"] = string(lbl) #fg.v[fg.id]
  fgl.IDs[lbl] = id

  # used for cloudgraph solving
  vert.attributes["ready"] = ready
  vert.attributes["backendset"] = 0

  # store user metadata
  vert.attributes["smalldata"] = smalldata
  nothing
end

"""
$(SIGNATURES)

Add a node (variable) to a graph. Use this over the other dispatches.
"""
function addNode!(fg::FactorGraph,
                  lbl::Symbol,
                  softtype::T;
                  N::Int=100,
                  autoinit::Bool=true,  # does init need to be separate from ready? TODO
                  ready::Int=1,
                  labels::Vector{<:AbstractString}=String[],
                  api::DataLayerAPI=dlapi,
                  uid::Int=-1,
                  smalldata=""  ) where {T <:InferenceVariable}
  #
  currid = fg.id+1
  if uid==-1
    fg.id=currid
  else
    currid = uid
  end
  # dims = dims != -1 ? dims : st.dims
  lblstr = string(lbl)
  vert = ExVertex(currid,lblstr)
  addNewVarVertInGraph!(fg, vert, currid, lbl, ready, smalldata)
  # dlapi.setupvertgraph!(fg, vert, currid, lbl) #fg.v[currid]
  # dodims = fg.dimID+1
  setDefaultNodeData!(vert, zeros(softtype.dims,N), zeros(0,0), 0, N, softtype.dims, initialized=!autoinit, softtype=softtype) # dodims

  vnlbls = union(string.(labels), softtype.labels, String["VARIABLE";])
  push!(vnlbls, fg.sessionname)

  api.addvertex!(fg, vert, labels=vnlbls)

  # fg.dimID+=dims # DONE -- drop this, rows indicate dimensions, move to last dimension
  push!(fg.nodeIDs, currid)

  vert
end

"""
$(SIGNATURES)

Add a node (variable) to a graph. Use this over the other dispatches.
"""
function addNode!(fg::FactorGraph,
                  lbl::Symbol,
                  softtype::Type{<:InferenceVariable};
                  N::Int=100,
                  autoinit::Bool=true,
                  ready::Int=1,
                  labels::Vector{<:AbstractString}=String[],
                  api::DataLayerAPI=dlapi,
                  uid::Int=-1,
                  # dims::Int=-1,
                  smalldata=""  )
  #
  addNode!(fg,
           lbl,
           softtype();
           N=N,
           autoinit=autoinit,
           ready=ready,
           labels=labels,
           api=api,
           uid=uid,
           smalldata=smalldata  )
end

function getVal(vA::Array{Graphs.ExVertex,1})
  warn("getVal(::Vector{ExVertex}) is obsolete, use getVal.(ExVertex) instead.")
  len = length(vA)
  vals = Array{Array{Float64,2},1}()
  cols = Array{Int,1}()
  push!(cols,0)
  rows = Array{Int,1}()
  for v in vA
      push!(vals, getVal(v))
      c = size(vals[end],2)
      r = size(vals[end],1)
      push!(cols, floor(Int,c))
      push!(rows, floor(Int,r))
  end
  cols = cumsum(cols)
  sc = cols[end]
  rw = floor(Int,rows[1])
  val = Array{Float64,2}(rw, sc)
  for i in 1:(len-1)
      val[:,(cols[i]+1):cols[i+1]] = vals[i]
  end
  val[:,(cols[len]+1):cols[len+1]] = vals[len] # and the last one
  return val
end

"""
    $(SIGNATURES)

Prepare the particle arrays `ARR` to be used for approximate convolution.  This function ensures that ARR has te same dimensions among all the parameters.  Function returns with ARR[sfidx] pointing at newly allocated deepcopy of the existing values in getVal(Xi[.index==solvefor]).  Return values `sfidx` is the element in ARR where `Xi.index==solvefor` and `maxlen` is length of all (possibly resampled) `ARR` contained particles.  Note `Xi` is order sensitive.
"""
function prepareparamsarray!(ARR::Array{Array{Float64,2},1},
            Xi::Vector{Graphs.ExVertex},
            N::Int,
            solvefor::Int  )
  #
  LEN = Int[]
  maxlen = N
  count = 0
  sfidx = 0

  for xi in Xi
    push!(ARR, getVal(xi))
    len = size(ARR[end], 2)
    push!(LEN, len)
    if len > maxlen
      maxlen = len
    end
    count += 1
    if xi.index == solvefor
      sfidx = count #xi.index
    end
  end
  SAMP=LEN.<maxlen
  for i in 1:count
    if SAMP[i]
      ARR[i] = KernelDensityEstimate.sample(getKDE(Xi[i]), maxlen)[1]
    end
  end

  # TODO --rather define reusable memory for the proposal
  # we are generating a proposal distribution, not direct replacement for existing memory and hence the deepcopy.
  if sfidx > 0 ARR[sfidx] = deepcopy(ARR[sfidx]) end
  return maxlen, sfidx
end

function parseusermultihypo(multihypo::Void)
  verts = Symbol[]
  mh = nothing
  return mh
end
function parseusermultihypo(multihypo::Union{Tuple,Vector{Float64}})
  mh = nothing
  if multihypo != nothing
    multihypo2 = Float64[multihypo...]
    # verts = Symbol.(multihypo[1,:])
    for i in 1:length(multihypo)
      if multihypo[i] > 0.999999
        multihypo2[i] = 0.0
      end
    end
    mh = Categorical(Float64[multihypo2...] )
  end
  return mh
end

function prepgenericconvolution(
            Xi::Vector{Graphs.ExVertex},
            usrfnc::T;
            multihypo::Union{Void, Distributions.Categorical}=nothing,
            threadmodel=MultiThreaded  ) where {T <: FunctorInferenceType}
      # multiverts::Vector{Symbol}=Symbol[]
  #
  ARR = Array{Array{Float64,2},1}()
  maxlen, sfidx = prepareparamsarray!(ARR, Xi, 0, 0)
  fldnms = fieldnames(usrfnc)
  zdim = typeof(usrfnc) != GenericMarginal ? size(getSample(usrfnc, 2)[1],1) : 0
  ccw = CommonConvWrapper(
          usrfnc,
          zeros(1,0),
          zdim,
          ARR,
          specialzDim = sum(fldnms .== :zDim) >= 1,
          partial = sum(fldnms .== :partial) >= 1,
          hypotheses=multihypo,
          threadmodel=threadmodel
        )
  #
  for i in 1:Threads.nthreads()
    ccw.cpt[i].factormetadata.variableuserdata = []
    ccw.cpt[i].factormetadata.solvefor = :null
    for xi in Xi
      push!(ccw.cpt[i].factormetadata.variableuserdata, getData(xi).softtype)
    end
  end
  return ccw
end

function setDefaultFactorNode!(
      fgl::FactorGraph,
      vert::Graphs.ExVertex,
      Xi::Vector{Graphs.ExVertex},
      usrfnc::T;
      multihypo::Union{Void,Tuple,Vector{Float64}}=nothing,
      threadmodel=MultiThreaded) where {T <: Union{FunctorInferenceType, InferenceType}}
  #
  ftyp = typeof(usrfnc) # maybe this can be T
  # @show "setDefaultFactorNode!", usrfnc, ftyp, T
  mhcat = parseusermultihypo(multihypo)
  # gwpf = prepgenericwrapper(Xi, usrfnc, getSample, multihypo=mhcat)
  ccw = prepgenericconvolution(Xi, usrfnc, multihypo=mhcat, threadmodel=threadmodel)

  m = Symbol(ftyp.name.module)

  # experimental wip
  data_ccw = FunctionNodeData{CommonConvWrapper{T}}(Int[], false, false, Int[], m, ccw)
  vert.attributes["data"] = data_ccw

  # existing interface
  # data = FunctionNodeData{GenericWrapParam{T}}(Int[], false, false, Int[], m, gwpf)
  # vert.attributes["data"] = data

  nothing
end

function addNewFncVertInGraph!(fgl::FactorGraph, vert::Graphs.ExVertex, id::Int, lbl::Symbol, ready::Int)
  vert.attributes = Graphs.AttributeDict() #fg.v[fg.id]
  vert.attributes["label"] = lbl #fg.v[fg.id]
  # fgl.f[id] = vert #  -- not sure if this is required, using fg.g.vertices
  fgl.fIDs[lbl] = id # fg.id,

  # used for cloudgraph solving
  vert.attributes["ready"] = ready
  vert.attributes["backendset"] = 0

  # for graphviz drawing
  vert.attributes["shape"] = "point"
  vert.attributes["width"] = 0.2
  nothing
end
addNewFncVertInGraph!{T <: AbstractString}(fgl::FactorGraph, vert::Graphs.ExVertex, id::Int, lbl::T, ready::Int) =
    addNewFncVertInGraph!(fgl,vert, id, Symbol(lbl), ready)

function isInitialized(vert::Graphs.ExVertex)::Bool
  return getData(vert).initialized
end
function isInitialized(fgl::FactorGraph, vsym::Symbol)::Bool
  # TODO, make cloudgraphs work and make faster by avoiding all the getVerts
  isInitialized(getVert(fgl, vsym))
end

"""
    $(SIGNATURES)

initialize destination variable nodes based on this factor in factor graph, fg, generally called
during addFactor!. Destination factor is first (singletons) or second (dim 2 pairwise) variable vertex in Xi.
"""
function doautoinit!(fgl::FactorGraph,
                     Xi::Vector{Graphs.ExVertex};
                     api::DataLayerAPI=dlapi,
                     singles::Bool=true,
                     N::Int=100)
  # Mighty inefficient function, since we only need very select fields nearby from a few neighboring nodes
  # do double depth search for variable nodes
  # TODO this should maybe stay localapi only...
  for xi in Xi
    if !isInitialized(xi)
      vsym = Symbol(xi.label)
      neinodes = ls(fgl, vsym)
      if (length(neinodes) > 1 || singles) # && !isInitialized(xi)
        # Which of the factors can be used for initialization
        useinitfct = Symbol[]
        # println("Consider all pairwise factors connected to $vsym...")
        for xifct in neinodes #potntlfcts
          xfneivarnodes = lsf(fgl, xifct)
          for vsym2 in xfneivarnodes
            # println("find all variables that are initialized for $vsym2")
            vert2 = getVert(fgl, vsym2)
            if (isInitialized(vert2) && sum(useinitfct .== xifct) == 0 ) || length(xfneivarnodes) == 1
              # OR singleton  TODO get faster version of isInitialized for database version
              # println("adding $xifct to init factors list")
              push!(useinitfct, xifct)
            end
          end
        end
        # println("Consider all singleton (unary) factors to $vsym...")

        # calculate the predicted belief over $vsym
        pts = predictbelief(fgl, vsym, useinitfct, api=api)
        setValKDE!(xi, pts)
        getData(xi).initialized = true
        api.updatevertex!(fgl, xi, updateMAPest=false)
      end
    end
  end

  nothing
end

"""
    $(SIGNATURES)

initialize destination variable nodes based on this factor in factor graph, fg, generally called
during addFactor!.  Destination factor is first (singletons) or second (dim 2 pairwise) variable vertex in Xi.
"""
function doautoinit!(fgl::FactorGraph,
                     xsyms::Vector{Symbol};
                     api::DataLayerAPI=dlapi,
                     singles::Bool=true,
                     N::Int=100)
  #
  verts = getVert.(fgl, xsyms, api=api)
  doautoinit!(fgl, verts, api=api, singles=singles, N=N)
end
function doautoinit!(fgl::FactorGraph,
                     xsym::Symbol;
                     api::DataLayerAPI=dlapi,
                     singles::Bool=true,
                     N::Int=100)
  #
  doautoinit!(fgl, [getVert(fgl, xsym, api=api);], api=api, singles=singles, N=N)
end

function ensureAllInitialized!(fgl::FactorGraph; api::DataLayerAPI=dlapi)
  xx, xl = ls(fgl)
  allvarnodes = union(xx, xl)
  for vsym in allvarnodes
    if !isInitialized(fgl, vsym)
      println("$vsym is not initialized, and will do so now...")
      doautoinit!(fgl, Graphs.ExVertex[getVert(fgl, vsym, api=api);], api=api, singles=true)
    end
  end
  nothing
end

function assembleFactorName(fgl::FactorGraph, Xi::Vector{Graphs.ExVertex}; maxparallel::Int=50)
  namestring = ""
  for vert in Xi #f.Xi
    namestring = string(namestring,vert.attributes["label"])
  end
  for i in 1:maxparallel
    tempnm = string(namestring, "f$i")
    if !haskey(fgl.fIDs, Symbol(tempnm))
      namestring = tempnm
      break
    end
    i != maxparallel ? nothing : error("Cannot currently add more than $(maxparallel) factors in parallel, please open an issue if this is too restrictive.")
  end
  return namestring
end

"""
    $(SIGNATURES)

Add factor with user defined type <: FunctorInferenceType to the factor graph object.  Define whether the automatic initialization of variables should be performed.  Use order sensitive `multihypo` keyword argument to define if any variables are related to data association uncertainty.
"""
function addFactor!(fgl::FactorGraph,
      Xi::Vector{Graphs.ExVertex},
      usrfnc::R;
      multihypo::Union{Void,Tuple,Vector{Float64}}=nothing,
      ready::Int=1,
      api::DataLayerAPI=dlapi,
      labels::Vector{T}=String[],
      uid::Int=-1,
      autoinit::Bool=true,
      threadmodel=MultiThreaded  ) where
        {R <: Union{FunctorInferenceType, InferenceType},
         T <: AbstractString}
  #
  currid = fgl.id+1
  if uid==-1
    fgl.id=currid
  else
    currid = uid
  end
  namestring = assembleFactorName(fgl, Xi)
  # fgl.id+=1
  newvert = ExVertex(currid,namestring)
  addNewFncVertInGraph!(fgl, newvert, currid, namestring, ready)
  setDefaultFactorNode!(fgl, newvert, Xi, deepcopy(usrfnc), multihypo=multihypo, threadmodel=threadmodel)
  push!(fgl.factorIDs,currid)

  for vert in Xi
    push!(getData(newvert).fncargvID, vert.index)
    # push!(newvert.attributes["data"].fncargvID, vert.index)
  end

  fnlbls = deepcopy(labels)
  fnlbls = union(fnlbls, String["FACTOR";])
  push!(fnlbls, fgl.sessionname)
  # TODO -- multiple accesses to DB with this method, must refactor!
  newvert = api.addvertex!(fgl, newvert, labels=fnlbls)  # used to be two be three lines up ##fgl.g
  for vert in Xi
    api.makeaddedge!(fgl, vert, newvert)
  end

  # TODO change this operation to update a conditioning variable
  if autoinit
    doautoinit!(fgl, Xi, api=api, singles=false)
  end

  return newvert
end

"""
    $(SIGNATURES)

Add factor with user defined type <: FunctorInferenceType to the factor graph object.  Define whether the automatic initialization of variables should be performed.  Use order sensitive `multihypo` keyword argument to define if any variables are related to data association uncertainty.
"""
function addFactor!(
      fgl::FactorGraph,
      xisyms::Vector{Symbol},
      usrfnc::R;
      multihypo::Union{Void,Tuple,Vector{Float64}}=nothing,
      ready::Int=1,
      api::DataLayerAPI=dlapi,
      labels::Vector{T}=String[],
      uid::Int=-1,
      autoinit::Bool=true,
      threadmodel=MultiThreaded  ) where
        {R <: Union{FunctorInferenceType, InferenceType},
         T <: AbstractString}
  #
  verts = Vector{Graphs.ExVertex}()
  for xi in xisyms
      push!( verts, api.getvertex(fgl,xi) )
  end
  addFactor!(fgl, verts, usrfnc, multihypo=multihypo, ready=ready, api=api, labels=labels, uid=uid, autoinit=autoinit, threadmodel=threadmodel )
end

function prtslperr(s)
  println(s)
  sleep(0.1)
  error(s)
end

# for computing the Bayes Net-----------------------------------------------------
function getEliminationOrder(fg::FactorGraph; ordering::Symbol=:qr)
    s = fg.nodeIDs
    lens = length(s)
    sf = fg.factorIDs
    lensf = length(sf)
    adjm, dictpermu = adjacency_matrix(fg.g,returnpermutation=true)
    permuteds = Vector{Int}(lens)
    permutedsf = Vector{Int}(lensf)
    for j in 1:length(dictpermu)
      semap = 0
      for i in 1:lens
        if dictpermu[j] == s[i]
          permuteds[i] = j#dictpermu[j]
          semap += 1
          if semap >= 2  break; end
        end
      end
      for i in 1:lensf
        if dictpermu[j] == sf[i]
          permutedsf[i] = j#dictpermu[j]
          semap += 1
          if semap >= 2  break; end
        end
      end
    end

      A=convert(Array{Int},adjm[permutedsf,permuteds]) # TODO -- order seems brittle
      p = Int[]
      if ordering==:chol
        p = cholfact(A'A,:U,Val{true})[:p] #,pivot=true
      elseif ordering==:qr
        q,r,p = qr(A,Val{true})
      else
        prtslperr("getEliminationOrder -- cannot do the requested ordering $(ordering)")
      end

    # we need the IDs associated with the Graphs.jl and our Type fg
    return  dictpermu[permuteds[p]] # fg.nodeIDs[p]
end

# lets create all the vertices first and then deal with the elimination variables thereafter
function addBayesNetVerts!(fg::FactorGraph, elimOrder::Array{Int,1})
  for p in elimOrder
    vert = getVert(fg, p, api=localapi)
    @show vert.label, getData(vert).BayesNetVertID
    if getData(vert).BayesNetVertID == 0
      fg.bnid+=1
      vert.attributes["data"].BayesNetVertID = p
      localapi.updatevertex!(fg, vert)
    else
      println("addBayesNetVerts -- something is very wrong, should not have a Bayes net vertex")
    end
  end
end

function addConditional!(fg::FactorGraph, vertID::Int, lbl, Si)
  bnv = getVert(fg, vertID, api=localapi) #fg.v[vertID]
  bnvd = getData(bnv) # bnv.attributes["data"]
  bnvd.separator = Si
  for s in Si
    push!(bnvd.BayesNetOutVertIDs, s)
  end
  localapi.updatevertex!(fg, bnv)
  nothing
end

function addChainRuleMarginal!(fg::FactorGraph, Si)
  lbls = String[]
  genmarg = GenericMarginal()
  Xi = Graphs.ExVertex[]
  for s in Si
    push!(Xi, getVert(fg, s, api=localapi))
  end
  println("adding marginal to")
  for x in Xi @show x.index end
  addFactor!(fg, Xi, genmarg, api=localapi, autoinit=false)
  nothing
end

# TODO -- Cannot have any CloudGraph calls at this level, must refactor
function rmVarFromMarg(fgl::FactorGraph, fromvert::Graphs.ExVertex, gm::Array{Graphs.ExVertex,1})
  for m in gm
    # get all out edges
    # get neighbors
    for n in localapi.outneighbors(fgl, m)
      if n.index == fromvert.index
        alleids = m.attributes["data"].edgeIDs
        i = 0
        for id in alleids
          i+=1
          edge = localapi.getedge(fgl, id)
          if edge != nothing # hack to avoid dictionary use case
            if edge.SourceVertex.exVertexId == m.index || edge.DestVertex.exVertexId == m.index
              warn("removing edge $(edge.neo4jEdgeId), between $(m.index) and $(n.index)")
              localapi.deleteedge!(fgl, edge)
              m.attributes["data"].edgeIDs = alleids[[collect(1:(i-1));collect((i+1):length(alleids))]]
              localapi.updatevertex!(fgl, m)
            end
          end
        end
      end
    end
    # if 0 edges, delete the marginal
    if length(localapi.outneighbors(fgl, m)) <= 1
      warn("removing vertex id=$(m.index)")
      localapi.deletevertex!(fgl,m)
    end
  end
  nothing
end

function buildBayesNet!(fg::FactorGraph, p::Array{Int,1})
    addBayesNetVerts!(fg, p)
    for v in p
      println()
      println("Eliminating $(v)")
      println("===============")
      println()
      # which variable are we eliminating

      # all factors adjacent to this variable
      fi = Int[]
      Si = Int[]
      gm = ExVertex[]
      # TODO -- optimize outneighbor calls like this
      vert = localapi.getvertex(fg,v)
      for fct in localapi.outneighbors(fg, vert)
        if (getData(fct).eliminated != true)
          push!(fi, fct.index)
          for sepNode in localapi.outneighbors(fg, fct)
            if sepNode.index != v && length(findin(sepNode.index,Si)) == 0
              push!(Si,sepNode.index)
            end
          end
          getData(fct).eliminated = true #fct.attributes["data"].eliminated = true
          localapi.updatevertex!(fg, fct) # TODO -- this might be a premature statement
        end

        if typeof(getData(fct).fnc) == GenericMarginal
          push!(gm, fct)
        end
      end

      if v != p[end]
        addConditional!(fg, v, "", Si)
        # not yet inserting the new prior p(Si) back into the factor graph
      end

      tuv = localapi.getvertex(fg, v) # TODO -- This may well through away valuable data
      tuv.attributes["data"].eliminated = true # fg.v[v].
      localapi.updatevertex!(fg, tuv)

      # TODO -- remove links from current vertex to any marginals
      rmVarFromMarg(fg, vert, gm)

      #add marginal on remaining variables... ? f(xyz) = f(x | yz) f(yz)
      # new function between all Si
      addChainRuleMarginal!(fg, Si)

    end
    nothing
end

# some plotting functions on the factor graph
function stackVertXY(fg::FactorGraph, lbl::String)
    v = dlapi.getvertex(fg,lbl)
    vals = getVal(v)
    X=vec(vals[1,:])
    Y=vec(vals[2,:])
    return X,Y
end

function getKDE(v::Graphs.ExVertex)
  return kde!(getVal(v), getBWVal(v)[:,1])
end

function getVertKDE(v::Graphs.ExVertex)
  return getKDE(v)
end

function getVertKDE(fgl::FactorGraph, id::Int; api::DataLayerAPI=dlapi)
  v = api.getvertex(fgl,id)
  return getKDE(v)
end
function getVertKDE(fgl::FactorGraph, lbl::Symbol; api::DataLayerAPI=dlapi)
  v = api.getvertex(fgl,lbl)
  return getKDE(v)
end

function drawCopyFG(fgl::FactorGraph)
  fgd = deepcopy(fgl)
  for (sym,vid) in fgd.IDs
    delete!(fgd.g.vertices[vid].attributes,"data")
    !haskey(fgd.g.vertices[vid].attributes,"frtend") ? nothing : delete!(fgd.g.vertices[vid].attributes,"frtend")
  end
  for (sym,vid) in fgd.fIDs
    delete!(fgd.g.vertices[vid].attributes,"data")
    !haskey(fgd.g.vertices[vid].attributes,"frtend") ? nothing : delete!(fgd.g.vertices[vid].attributes,"frtend")
  end
  return fgd
end

function writeGraphPdf(fgl::FactorGraph;
                       pdfreader::Union{Void, String}="evince",
                       filename::AS="/tmp/fg.pdf"  ) where {AS <: AbstractString}
  #
  fgd = drawCopyFG(fgl)
  println("Writing factor graph file")
  dotfile = split(filename, ".pdf")[1]*".dot"
  fid = open(dotfile,"w")
  write(fid,Graphs.to_dot(fgd.g))
  close(fid)
  run(`dot $(dotfile) -Tpdf -o $(filename)`)
  try
    pdfreader != nothing ? (@async run(`$(pdfreader) $(filename)`)) : nothing
  catch e
    warn("not able to show $(filename) with pdfreader=$(pdfreader). Exception e=$(e)")
  end
  nothing
end

function expandEdgeListNeigh!(fgl::FactorGraph,
                              vertdict::Dict{Int,Graphs.ExVertex},
                              edgedict::Dict{Int,Graphs.Edge{Graphs.ExVertex}})
  #asfd
  for vert in vertdict
    for newedge in out_edges(vert[2],fgl.g)
      if !haskey(edgedict, newedge.index)
        edgedict[newedge.index] = newedge
      end
    end
  end

  nothing
end

# dictionary of unique vertices from edgelist
function expandVertexList!(fgl::FactorGraph,
  edgedict::Dict{Int,Graphs.Edge{Graphs.ExVertex}},
  vertdict::Dict{Int,Graphs.ExVertex})

  # go through all source and target nodes
  for edge in edgedict
    if !haskey(vertdict, edge[2].source.index)
      vertdict[edge[2].source.index] = edge[2].source
    end
    if !haskey(vertdict, edge[2].target.index)
      vertdict[edge[2].target.index] = edge[2].target
    end
  end
  nothing
end

function edgelist2edgedict(edgelist::Array{Graphs.Edge{Graphs.ExVertex},1})
  edgedict = Dict{Int,Graphs.Edge{Graphs.ExVertex}}()
  for edge in edgelist
    edgedict[edge.index] = edge
  end
  return edgedict
end

# TODO -- convert to use add_vertex! instead, since edges type must be made also
function addVerticesSubgraph(fgl::FactorGraph,
    fgseg::FactorGraph,
    vertdict::Dict{Int,Graphs.ExVertex})

    for vert in vertdict
      fgseg.g.vertices[vert[1]] = vert[2]
      if haskey(fgl.v,vert[1])
        fgseg.g.vertices[vert[1]] = vert[2]
        fgseg.IDs[Symbol(vert[2].label)] = vert[1]

        # add edges going in opposite direction
        elr = Graphs.out_edges(vert[2], fgl.g)
        len = length(elr)
        keeprm = trues(len)
        j = 0
        for i in 1:len
          if !haskey(vertdict, elr[i].target.index) # a function node in set, so keep ref
            keeprm[i] = false
            j+=1
          end
        end
        if j < len
          elridx = elr[1].source.index
          fgseg.g.inclist[elridx] = elr[keeprm]
        end
      elseif haskey(fgl.f, vert[1])
        fgseg.f[vert[1]] = vert[2] # adding element to subgraph
        fgseg.fIDs[Symbol(vert[2].label)] = vert[1]
        # get edges associated with function nodes and push edges onto incidence list
        el = Graphs.out_edges(vert[2], fgl.g)
        elidx = el[1].source.index
        fgseg.g.inclist[elidx] = el # okay because treating function nodes only
        fgseg.g.nedges += length(el)
      else
        error("Unknown type factor graph vertex type, something is wrong")
      end
    end
    nothing
end

# NOTICE, nodeIDs and factorIDs are not brough over by this method yet
# must sort out for incremental updates
function genSubgraph(fgl::FactorGraph,
    vertdict::Dict{Int,Graphs.ExVertex})
    # edgedict::Dict{Int,Graphs.Edge{Graphs.ExVertex}},

  fgseg = FactorGraph() # new handle for just a segment of the graph
  fgseg.g = Graphs.inclist(Graphs.ExVertex,is_directed=false)

  fgseg.v = Dict{Int,Graphs.ExVertex}()
  fgseg.f = Dict{Int,Graphs.ExVertex}()
  fgseg.IDs = Dict{AbstractString,Int}()
  fgseg.fIDs = Dict{AbstractString,Int}()

  # TODO -- convert to use empty constructor since Graphs.incdict now works
  fgseg.g.vertices = Array{Graphs.ExVertex,1}(length(fgl.g.vertices))
  fgseg.g.inclist = Array{Array{Graphs.Edge{Graphs.ExVertex},1},1}(length(fgl.g.inclist))

  addVerticesSubgraph(fgl, fgseg, vertdict)

  fgseg.id = fgl.id
  fgseg.bnid = fgl.bnid
  fgseg.dimID = fgl.dimID

  return fgseg
end

function getShortestPathNeighbors(fgl::FactorGraph;
    from::Graphs.ExVertex=nothing,
    to::Graphs.ExVertex=nothing,
    neighbors::Int=0 )

  edgelist = shortest_path(fgl.g, ones(num_edges(fgl.g)), from, to)
  vertdict = Dict{Int,Graphs.ExVertex}()
  edgedict = edgelist2edgedict(edgelist)
  expandVertexList!(fgl, edgedict, vertdict) # grow verts
  for i in 1:neighbors
    expandEdgeListNeigh!(fgl, vertdict, edgedict) # grow edges
    expandVertexList!(fgl, edgedict, vertdict) # grow verts
  end
  return vertdict
end

function subgraphShortestPath(fgl::FactorGraph;
    from::Graphs.ExVertex=nothing,
    to::Graphs.ExVertex=nothing,
    neighbors::Int=0  )

  vertdict = getShortestPathNeighbors(fgl, from=from, to=to, neighbors=neighbors)
  return genSubgraph(fgl, vertdict)
end

# explore all shortest paths combinations in verts, add neighbors and reference subgraph
function subgraphFromVerts(fgl::FactorGraph,
    verts::Dict{Int,Graphs.ExVertex};
    neighbors::Int=0  )

    allverts = Dict{Int,Graphs.ExVertex}()
    allkeys = collect(keys(verts))
    len = length(allkeys)
    # union all shortest path combinations in a vertdict
    for i in 1:len, j in (i+1):len
      from = verts[allkeys[i]]
      to = verts[allkeys[j]]
      vertdict = getShortestPathNeighbors(fgl, from=from, to=to, neighbors=neighbors)
      for vert in vertdict
        if !haskey(allverts, vert[1])
          allverts[vert[1]] = vert[2]
        end
      end
    end

  return genSubgraph(fgl, allverts)
end

# explore all shortest paths combinations in verts, add neighbors and reference subgraph
# Using unique index into graph data structure
function subgraphFromVerts(fgl::FactorGraph,
    verts::Array{Int,1};
    neighbors::Int=0  )

  vertdict = Dict{Int,Graphs.ExVertex}()
  for vert in verts
    vertdict[vert] = fgl.g.vertices[vert]
  end

  return subgraphFromVerts(fgl,vertdict,neighbors=neighbors)
end

# explore all shortest paths combinations in verts, add neighbors and reference subgraph
# Using unique index into graph data structure
function subgraphFromVerts(fgl::FactorGraph,
    verts::Array{String,1};
    neighbors::Int=0  )

  vertdict = Dict{Int,Graphs.ExVertex}()
  for vert in verts
    id = -1
    if haskey(fgl.IDs, vert)
      id = fgl.IDs[Symbol(vert)]
    else
      error("FactorGraph01 only supports variable node subgraph search at this point")
    end
    vertdict[id] = fgl.g.vertices[id]
  end

  return subgraphFromVerts(fgl,vertdict,neighbors=neighbors)
end


text:
let
    @test 1 == 1
end


text:
# Preprocessing to narrow what we schedule

# Apply this to coverage
function get_window(data)
    start_index = length(data[1])
    end_index = 1

    # Scan 1: Find indices
    for datum in data
        first_nonzero = 1
        last_nonzero = length(datum)

        for i=1:length(datum)
            first_nonzero = i
            if datum[i] != 0
                break
            end
        end

        for j=length(datum):-1:1
            last_nonzero = j
            if datum[j] != 0
                break
            end
        end

        # Avoid case of all zeros
        if first_nonzero < last_nonzero
            if first_nonzero < start_index
                start_index = first_nonzero
            end
            if last_nonzero > end_index
                end_index = last_nonzero
            end
        end
    end

    # Check for valid window
    if start_index >= end_index
        # No window detected
        return data
    end

    # time_delta is used for adjusting time_between_shifts
    time_delta = (start_index - 1) + (length(data[1]) - end_index)
    # time_delta is what should be applied to env
    return start_index, end_index, time_delta
end

# Input can be availability, coverage, etc - by the week
function apply_window(input, start_index, end_index)
    input = deepcopy(input) # passed by reference
    for t in 1:length(input)
        input[t] = input[t][start_index:end_index]
    end
    return input
end

# Remove the window from shifts.
function remove_window(shifts, start_index)
    shifts = deepcopy(shifts)
    offset = start_index - 1 # one-indexed language!
    for e in keys(shifts)
        for shift in shifts[e]
            shift["start"] += offset
        end
    end
    return shifts
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
module Serializer

import SearchLight

function serialize end

function deserialize end

function serializables(m::Type{T})::Vector{Symbol} where {T<:SearchLight.AbstractModel}
  Symbol[]
end

end


text:
module TextGrid

include("features/prattFeatures.jl")
export feature

include("features/returnCMD.jl")
export command

include("TextGrid/extractor.jl")
export extractTextGrid

end


text:
# http://sunlightlabs.github.io/datacommons
global INFLUENCE_EXPLORER_API = URI("http://transparencydata.com/")

# -------

abstract Entity

type Org <: Entity
    id
end

type Ind <: Entity
    id
end

type Pol <: Entity
    id
end

print(io::IO, e::Entity) = print(io, e.id)

# -------

macro sunlight_method__id_limit_cycle(method_name, entity_type, route_fn)
    quote
        function $(esc(symbol(method_name)))(auth::String, entity::($entity_type); limit = nothing, cycle = nothing, options...)
            route = $route_fn(entity)

            args = Dict()
            limit != nothing && (args["limit"] = limit)
            cycle != nothing && (args["cycle"] = cycle)

            sunlight_get(auth, INFLUENCE_EXPLORER_API, route, args; options...)
        end

        $(esc(symbol(method_name)))(entity; auth = "", options...) = $(esc(symbol(method_name)))(auth, entity; options...)
    end
end

# -------

macro sunlight_method__id_cycle(method_name, entity_type, route_fn)
    quote
        function $(esc(symbol(method_name)))(auth::String, entity::($entity_type); cycle = nothing, options...)
            route = $route_fn(entity)

            args = Dict()
            cycle != nothing && (args["cycle"] = cycle)

            sunlight_get(auth, INFLUENCE_EXPLORER_API, route, args; options...)
        end

        $(esc(symbol(method_name)))(entity; auth = "", options...) = $(esc(symbol(method_name)))(auth, entity; options...)
    end
end

# -------

macro sunlight_method__id(method_name, entity_type, route_fn)
    quote
        function $(esc(symbol(method_name)))(auth::String, entity::($entity_type); options...)
            route = $route_fn(entity)
            sunlight_get(auth, INFLUENCE_EXPLORER_API, route, Dict(); options...)
        end

        $(esc(symbol(method_name)))(entity; auth = "", options...) = $(esc(symbol(method_name)))(auth, entity; options...)
    end
end

# -------

party_breakdown(::String, ::String; options...) = error("'party_breakdown(entity_id::String)' is ambiguous; the method is defined for Organizations and Individuals. Try passing in `Org(entity_id)` or `Ind(entity_id)`")
party_breakdown(::String; options...) = error("'party_breakdown(entity_id::String)' is ambiguous; the method is defined for Organizations and Individuals. Try passing in `Org(entity_id)` or `Ind(entity_id)`")

lobbying_issues(::String, ::String; options...) = error("'lobbying_issues(entity_id::String)' is ambiguous; the method is defined for Organizations and Individuals. Try passing in `Org(entity_id)` or `Ind(entity_id)`")
lobbying_issues(::String; options...) = error("'lobbying_issues(entity_id::String)' is ambiguous; the method is defined for Organizations and Individuals. Try passing in `Org(entity_id)` or `Ind(entity_id)`")

fec_summary(::String, ::String; options...) = error("'fec_summary(entity_id::String)' is ambiguous; the method is defined for Organizations and Politicians. Try passing in `Org(entity_id)` or `Pol(entity_id)`")
fec_summary(::String, options...) = error("'fec_summary(entity_id::String)' is ambiguous; the method is defined for Organizations and Politicians. Try passing in `Org(entity_id)` or `Pol(entity_id)`")

# -------

@sunlight_method__id_limit_cycle "top_contributors" Union(String, Pol) id->"/api/1.0/aggregates/pol/$id/contributors.json"

@sunlight_method__id_limit_cycle "top_industries" Union(String, Pol) id->"/api/1.0/aggregates/pol/$id/contributors/industries.json"

@sunlight_method__id_cycle "unknown_industries" Union(String, Pol) id->"/api/1.0/aggregates/pol/$id/contributors/industries_unknown.json"

@sunlight_method__id_limit_cycle "top_sectors" Union(String, Pol) id->"/api/1.0/aggregates/pol/$id/contributors/sectors.json"

@sunlight_method__id_cycle "local_breakdown" Union(String, Pol) id->"/api/1.0/aggregates/pol/$id/contributors/local_breakdown.json"

@sunlight_method__id_cycle "contributor_breakdown" Union(String, Pol) id->"/api/1.0/aggregates/pol/$id/contributors/type_breakdown.json"

@sunlight_method__id "fec_summary" Union(String, Pol) id->"/api/1.0/aggregates/pol/$id/fec_summary.json"

@sunlight_method__id "fec_independent_expenditures" Union(String, Pol) id->"/api/1.0/aggregates/pol/$id/fec_indexp.json"

@sunlight_method__id_limit_cycle "top_recipient_organizations" Union(String, Ind) id->"/api/1.0/aggregates/indiv/$id/recipient_orgs.json"

@sunlight_method__id_limit_cycle "top_recipient_politicians" Union(String, Ind) id->"/api/1.0/aggregates/indiv/$id/recipient_pols.json"

@sunlight_method__id_cycle "party_breakdown" Union(String, Ind) id->"/api/1.0/aggregates/indiv/$id/recipients/party_breakdown.json"

@sunlight_method__id_limit_cycle "lobbying_registrants" Union(String, Ind) id->"/api/1.0/aggregates/indiv/$id/registrants.json"

@sunlight_method__id_limit_cycle "lobbying_clients" Union(String, Ind) id->"/api/1.0/aggregates/indiv/$id/clients.json"

@sunlight_method__id_limit_cycle "lobbying_issues" Union(String, Ind) id->"/api/1.0/aggregates/indiv/$id/issues.json"

@sunlight_method__id_limit_cycle "top_recipients" Union(String, Org) id->"/api/1.0/aggregates/org/$id/recipients.json"

@sunlight_method__id_limit_cycle "pac_recipients" Union(String, Org) id->"/api/1.0/aggregates/org/$id/recipient_pacs.json"

@sunlight_method__id_cycle "party_breakdown" Union(String, Org) id->"/api/1.0/aggregates/org/$id/recipients/party_breakdown.json"

@sunlight_method__id_cycle "state_federal_breakdown" Union(String, Org) id->"/api/1.0/aggregates/org/$id/recipients/level_breakdown.json"

@sunlight_method__id_limit_cycle "lobbing_registrants" Union(String, Org) id->"/api/1.0/aggregates/org/$id/registrants.json"

@sunlight_method__id_limit_cycle "lobbying_issues" Union(String, Org) id->"/api/1.0/aggregates/org/$id/issues.json"

@sunlight_method__id_limit_cycle "bills" Union(String, Org) id->"/api/1.0/aggregates/org/$id/bills.json"

@sunlight_method__id_limit_cycle "lobbyists" Union(String, Org) id->"/api/1.0/aggregates/org/$id/lobbyists.json"

@sunlight_method__id_limit_cycle "registrant_clients" Union(String, Org) id->"/api/1.0/aggregates/org/$id/registrant/clients.json"

@sunlight_method__id_limit_cycle "registrant_issues" Union(String, Org) id->"/api/1.0/aggregates/org/$id/registrant/issues.json"

@sunlight_method__id_limit_cycle "registrant_bills" Union(String, Org) id->"/api/1.0/aggregates/org/$id/registrant/bills.json"

@sunlight_method__id_limit_cycle "registrant_lobbyists" Union(String, Org) id->"/api/1.0/aggregates/org/$id/registrant/lobbyists.json"

@sunlight_method__id_limit_cycle "mentions_in_regulations" Union(String, Org) id->"/api/1.0/aggregates/org/$id/regulations_text.json"

@sunlight_method__id_limit_cycle "regulatory_comment_submissions" Union(String, Org) id->"/api/1.0/aggregates/org/$id/regulations_submitter.json"

@sunlight_method__id_limit_cycle "faca_memberships" Union(String, Org) id->"/api/1.0/aggregates/org/$id/faca.json"

@sunlight_method__id "fec_summary" Union(String, Org) id->"/api/1.0/aggregates/org/$id/fec_summary.json"

@sunlight_method__id_limit_cycle "top_industry_organizations" String id->"/api/1.0/aggregates/industry/$id/orgs.json"

# -------

function entity_search(auth::String, search_str; entity_type = nothing, options...)
    args = {
        "search" => search_str
    }
    entity_type != nothing && (args["type"] = entity_type)

    sunlight_get(auth, INFLUENCE_EXPLORER_API, "/api/1.0/entities.json", args; options...)
end

entity_search(search_str; auth = "", options...) = entity_search(auth, search_str; options...)

# -------

function entity_info(auth::String, entity_id; cycle = nothing, options...)
    args = Dict()
    cycle != nothing && (args["cycle"] = cycle)

    sunlight_get(auth, INFLUENCE_EXPLORER_API, "/api/1.0/entities/$entity_id.json", args; options...)
end

entity_info(entity_id; auth = "", options...) = entity_info(auth, entity_id; options...)

# -------

function top_politicians(auth::String; limit = 16, cycle = nothing, options...)
    args = Dict()
    cycle != nothing && (args["cycle"] = cycle)

    sunlight_get(auth, INFLUENCE_EXPLORER_API, "/api/1.0/aggregates/pols/top_$limit.json", args; options...)
end

top_politicians(; auth = "", options...) = top_politicians(auth; options...)

# -------

function top_individuals(auth::String; limit = 16, cycle = nothing, options...)
    args = Dict()
    cycle != nothing && (args["cycle"] = cycle)

    sunlight_get(auth, INFLUENCE_EXPLORER_API, "/api/1.0/aggregates/indivs/top_$limit.json", args; options...)
end

top_individuals(; auth = "", options...) = top_individuals(auth; options...)

# -------

function top_organizations(auth::String; limit = 16, cycle = nothing, options...)
    args = Dict()
    cycle != nothing && (args["cycle"] = cycle)

    sunlight_get(auth, INFLUENCE_EXPLORER_API, "/api/1.0/aggregates/orgs/top_$limit.json", args; options...)
end

top_organizations(; auth = "", options...) = top_organizations(auth; options...)

# -------

function top_industries(auth::String; limit = 16, cycle = nothing, options...)
    args = Dict()
    cycle != nothing && (args["cycle"] = cycle)

    sunlight_get(auth, INFLUENCE_EXPLORER_API, "/api/1.0/aggregates/industries/top_$limit.json", args; options...)
end

top_industries(; auth = "", options...) = top_industries(auth; options...)


text:
#!/bin/bash
# -*- mode: julia -*-
#=
JULIA="${JULIA:-julia --color=yes --startup-file=no}"
export JULIA_PROJECT="$(dirname ${BASH_SOURCE[0]})"

set -ex
${JULIA} -e 'using Pkg; Pkg.instantiate();
             Pkg.develop(PackageSpec(path=pwd()))'
exec ${JULIA} "${BASH_SOURCE[0]}" "$@"
=#

using Documenter, JuliaPythonSnippets

makedocs(;
    modules=[JuliaPythonSnippets],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "arrays.md",
        "introspection.md",
        "miscellaneous.md",
    ],
    repo="https://github.com/tkf/julia-python-snippets/blob/{commit}{path}#L{line}",
    sitename="julia-python-snippets",
    authors="Takafumi Arakaki",
    assets=[],
    strict=true,
)

deploydocs(;
    repo="github.com/tkf/julia-python-snippets",
)

max_stars_repo_name: tkf/julia-python-snippets

max_issues_repo_name: tkf/julia-python-snippets

max_forks_repo_name: tkf/julia-python-snippets


text:
# ---
# notebook: false
# ---


text:
using Test
using TarIterators
using BoundedStreams
using Tar
using CodecZlib
using TranscodingStreams

@testset "TarIterators" begin include("tariterators.jl"); end


text:
# OpticalInterface.jl is included in Geometry.jl so is omitted here

include("OpticalRay.jl")
include("Emitter.jl")
include("Emitter2.jl")
include("Fresnel.jl")
include("Paraxial.jl")
include("Grating.jl")
include("LensAssembly.jl")
include("HierarchicalImage.jl")
include("OpticalSystem.jl")
include("Lenses.jl")
include("Eye.jl")
include("Fields.jl")


text:
using Dates

"""
    {{meta}}

Plug in specific meta information for a blog page. The `meta` local page
variable should be given as a list of tuples of pairs like so:
```
@def meta = [("property"=>"og:video", "content"=>"http://example.com/"),
             ("propery"=>"og:title", "content"=>"The Rock")]
```
"""
function hfun_meta()
    m = locvar(:meta)
    io = IOBuffer()
    for tuple in locvar(:meta)
        write(io, "<meta ")
        for (prop, val) in tuple
            write(io, "$prop=\"$val\" ")
        end
        write(io, ">")
    end
    return String(take!(io))
end

"""
    {{blogposts}}

Plug in the list of blog posts contained in the `/blog/` folder.
"""
function hfun_blogposts()
    curyear = year(Dates.today())
    io = IOBuffer()
    for year in curyear:-1:2012
        ys = "$year"
        year < curyear && write(io, "\n**$year**\n")
        for month in 12:-1:1
            ms = "0"^(month < 10) * "$month"
            base = joinpath("blog", ys, ms)
            isdir(base) || continue
            posts = filter!(p -> endswith(p, ".md"), readdir(base))
            days  = zeros(Int, length(posts))
            lines = Vector{String}(undef, length(posts))
            for (i, post) in enumerate(posts)
                ps  = splitext(post)[1]
                url = "/blog/$ys/$ms/$ps/"
                surl = strip(url, '/')
                title = pagevar(surl, :title)
                pubdate = pagevar(surl, :published)
                if isnothing(pubdate)
                    date    = "$ys-$ms-01"
                    days[i] = 1
                else
                    date    = Date(pubdate, dateformat"d U Y")
                    days[i] = day(date)
                end
                lines[i] = "\n[$title]($url) $date \n"
            end
            # sort by day
            foreach(line -> write(io, line), lines[sortperm(days, rev=true)])
        end
    end
    # markdown conversion adds `<p>` beginning and end but
    # we want to  avoid this to avoid an empty separator
    r = Franklin.fd2html(String(take!(io)), internal=true)
    return r
end

"""
    {{recentblogposts}}

Input the 3 latest blog posts.
"""
function hfun_recentblogposts()
    curyear = Dates.Year(Dates.today()).value
    ntofind = 3
    nfound  = 0
    recent  = Vector{Pair{String,Date}}(undef, ntofind)
    for year in curyear:-1:2019
        for month in 12:-1:1
            ms = "0"^(1-div(month, 10)) * "$month"
            base = joinpath("blog", "$year", "$ms")
            isdir(base) || continue
            posts = filter!(p -> endswith(p, ".md"), readdir(base))
            days  = zeros(Int, length(posts))
            surls = Vector{String}(undef, length(posts))
            for (i, post) in enumerate(posts)
                ps       = splitext(post)[1]
                surl     = "blog/$year/$ms/$ps"
                surls[i] = surl
                pubdate  = pagevar(surl, :published)
                days[i]  = isnothing(pubdate) ?
                                1 : day(Date(pubdate, dateformat"d U Y"))
            end
            # go over month post in antichronological orders
            sp = sortperm(days, rev=true)
            for (i, surl) in enumerate(surls[sp])
                recent[nfound + 1] = (surl => Date(year, month, days[sp[i]]))
                nfound += 1
                nfound == ntofind && break
            end
            nfound == ntofind && break
        end
        nfound == ntofind && break
    end
    #
    io = IOBuffer()
    for (surl, date) in recent
        url   = "/$surl/"
        title = pagevar(surl, :title)
        sdate = "$(day(date)) $(monthname(date)) $(year(date))"
        blurb = pagevar(surl, :rss)
        write(io, """

$title

$date

$blurb

            """)
    end
    return String(take!(io))
end

"""
    {{redirect url}}

Creates a HTML layout for a redirect to `url`.
"""
function hfun_redirect(url)
    s = """

    <!doctype html>

        <meta http-equiv="refresh" content="0; url=$(url[1])">

    """
    return s
end


text:
using IcanProj
using NetCDF
using ProgressMeter

# Settings

path = "/data04/jmg/fsm_simulations"

time_start = DateTime(2008,1,1)

# Load soil parameters

soil_param = read_soil_params(path)

# Dimensions

dim_time = 365*4*8
dim_space = size(soil_param, 1)

# Auxiliary variables

time_str = Dates.format.([time_start + i*Dates.Hour(3) for i in 0:dim_time-1], "yyyy-mm-dd HH:MM:SS")
lon = convert(Array{Float64}, soil_param[:lon])
lat = convert(Array{Float64}, soil_param[:lat])
id = convert(Array{Int64}, soil_param[:gridcel])

# Create netcdf files

id_desc = "ind_senorge"

var = ["rainf",
       "snowf",
       "tair",
       "iswr",
       "ilwr",
       "pres",
       "rhum",
       "wind"]

var_atts = [Dict("units" => "mm/tstep"),
            Dict("units" => "mm/tstep"),
            Dict("units" => "C"),
            Dict("units" => "W/m2"),
            Dict("units" => "W/m2"),
            Dict("units" => "kPa"),
            Dict("units" => "%"),
            Dict("units" => "m/s")]

fn = joinpath.(path, "netcdf", var .* "_1km.nc")

for i in 1:length(var)

    create_netcdf(fn[i], var[i], var_atts[i], dim_time, dim_space, time_str, id, id_desc)

end

# Loop over files

@showprogress "Computing for ..." for i in 1:size(soil_param,1)

    lat_str = @sprintf("%0.5f", lat[i])
    lon_str = @sprintf("%0.5f", lon[i])

    file_src = joinpath(path, "results/metdata_$(lat_str)_$(lon_str)")

    rainf, snowf, tair, iswr, ilwr, pres, rhum, wind = read_mtclim_fsm(file_src)

    ncwrite(rainf[1:dim_time], fn[1], var[1], start=[1,i], count=[-1,1])
    ncwrite(snowf[1:dim_time], fn[2], var[2], start=[1,i], count=[-1,1])
    ncwrite(tair[1:dim_time],  fn[3], var[3], start=[1,i], count=[-1,1])
    ncwrite(iswr[1:dim_time],  fn[4], var[4], start=[1,i], count=[-1,1])
    ncwrite(ilwr[1:dim_time],  fn[5], var[5], start=[1,i], count=[-1,1])
    ncwrite(pres[1:dim_time],  fn[6], var[6], start=[1,i], count=[-1,1])
    ncwrite(rhum[1:dim_time],  fn[7], var[7], start=[1,i], count=[-1,1])
    ncwrite(wind[1:dim_time],  fn[8], var[8], start=[1,i], count=[-1,1])

end

ncclose()