text:
"""
$(TYPEDEF)

Structure containing plot information for Plotly.js
"""
mutable struct PlutoPlotlyPlot  <: AbstractPlutoVistaBackend

    # command list passed to javascript
    jsdict::Dict{String,Any}
    args

    # size in canvas coordinates
    w::Float64
    h::Float64

    update::Bool
    # uuid for identifying html element
    uuid::UUID
    PlutoPlotlyPlot(::Nothing)=new()
end

"""
$(TYPEDSIGNATURES)

Create a plotly plot.

"""
function PlutoPlotlyPlot(;resolution=(300,300), kwargs...)
    p=PlutoPlotlyPlot(nothing)
    p.uuid=uuid1()
    p.jsdict=Dict{String,Any}("cmdcount" => 0)
    p.w=resolution[1]
    p.h=resolution[2]
    default_args=(limits=(1,-1),
                  xlimits=(1,-1),
                  xlabel="",
                  ylabel="",
                  title="",
                  xscale=:linear,
                  yscale=:linear,
                  legend=:none,
                  titlefontsize=12,
                  axisfontsize=10,
                  tickfontsize=10,
                  legendfontsize=10,
                  colorbarticks=:default,
                  clear=false)
    p.args=merge(default_args,kwargs)
    p.update=false
    p
end

const plutoplotlyplot = read(joinpath(@__DIR__, "..", "assets", "plutoplotlyplot.js"), String)

"""
$(TYPEDSIGNATURES)

Show plotly plot.
"""
function Base.show(io::IO, ::MIME"text/html", p::PlutoPlotlyPlot)
    result="""

        """
    # updating only works when the div remains as output from another cell
    # so we can't create a plot and update it in the cell with the draing commands
    if !p.update
        result=result*"""

        """
    end
    p.update=true
    write(io,result)
end

#translate Julia attribute symbols to pyplot-speak
const mshapes=Dict(
    :dtriangle => "triangle-down",
    :utriangle => "triangle-up",
    :rtriangle => "triangle-right",
    :ltriangle => "triangle-left",
    :circle => "circle",
    :square => "square",
    :cross => "cross",
    :+ => "cross",
    :xcross => "x",
    :x => "x",
    :diamond => "diamond",
    :star5 => "star",
    :pentagon => "pentagon",
    :hexagon => "hexagon",
    :none => "none"
)

function axisargs!(p,args)
    xscale= args[:xscale]==:identity ? :linear :  args[:xscale]
    yscale= args[:yscale]==:identity ? :linear :  args[:yscale]

    parameter!(p,"ylimits",collect(Float32,args[:limits]))
    parameter!(p,"xlimits",collect(Float32,args[:xlimits]))
    parameter!(p,"xlabel",args[:xlabel])
    parameter!(p,"ylabel",args[:ylabel])
    parameter!(p,"xaxis",String(xscale))
    parameter!(p,"yaxis",String(yscale))
    parameter!(p,"title",args[:title])
    parameter!(p,"legendfontsize",args[:legendfontsize])
    parameter!(p,"axisfontsize"   ,args[:axisfontsize])
    parameter!(p,"titlefontsize"  ,args[:titlefontsize])
    parameter!(p,"tickfontsize"   ,args[:tickfontsize])
end

"""
$(SIGNATURES)

1D plotly.js plot
"""
function plot!(p::PlutoPlotlyPlot,x,y; kwargs...)

    default_args=(label="",
                  color=:auto,
                  linewidth=2,
                  linestyle=:solid,
                  markersize=6,
                  markercount=10,
                  markertype=:none)
    args=merge(p.args,default_args)
    args=merge(args,kwargs)

    if args[:clear]
        p.jsdict=Dict{String,Any}("cmdcount" => 0)
    end

    command!(p,"plot")
    parameter!(p,"x",collect(x))
    parameter!(p,"y",collect(y))
    parameter!(p,"label",args[:label])
    parameter!(p,"linewidth",args[:linewidth])
    parameter!(p,"markercount",args[:markercount])
    parameter!(p,"markertype",mshapes[args[:markertype]])
    parameter!(p,"markersize",args[:markersize])
    parameter!(p,"linestyle",String(args[:linestyle]))

    axisargs!(p,args)

    if args[:legend]==:none || args[:legend]==""
        parameter!(p,"showlegend",0)
    else
        parameter!(p,"showlegend",1)
        slegend=String(args[:legend])
        parameter!(p,"legendxpos",slegend[1:1])
        parameter!(p,"legendypos",slegend[2:2])
    end

    parameter!(p,"clear",args[:clear] ? 1 : 0)

    if args[:color] == :auto
        parameter!(p,"color","auto")
    else
        rgb=RGB(args[:color])
        rgb=[rgb.r,rgb.g,rgb.b]
        rgb=UInt8.(floor.(rgb*255))
        parameter!(p,"color",rgb)
    end
    p
end

"""
$(SIGNATURES)

Experimental. Plot piecewise linear function on  triangular grid given as "heatmap" and
with isolines using Plotly's mesh3d.
"""
function tricontour!(p::PlutoPlotlyPlot,pts, tris,f;kwargs...)

    default_args=(colormap=:viridis, levels=0, aspect=1)
    args=merge(p.args,default_args)
    args=merge(args,kwargs)

    levels,crange,colorbarticks=GridVisualize.isolevels(args,f)

    zval=0.0
    p.jsdict=Dict{String,Any}("cmdcount" => 0)

    command!(p,"tricontour")
    (fmin,fmax)=extrema(f)

    parameter!(p,"x",pts[1,:])
    parameter!(p,"y",pts[2,:])
    parameter!(p,"z",fill(zval,length(f)))
    parameter!(p,"f",f)
    parameter!(p,"aspect",[1.0,args[:aspect],1.0])

    axisargs!(p,args)

    stops=collect(0:0.01:1)
    rgb=reinterpret(Float64,get(colorschemes[args[:colormap]],stops,(0,1)))
    rgb=UInt8.(floor.(rgb*255))

    parameter!(p,"cstops",stops)
    parameter!(p,"colors",rgb)
    parameter!(p,"colorbarticks",colorbarticks)

    parameter!(p,"i",tris[1,:].-1)
    parameter!(p,"j",tris[2,:].-1)
    parameter!(p,"k",tris[3,:].-1)

    (fmin,fmax)=extrema(f)

    if levels==0
        parameter!(p,"iso_x","none")
        parameter!(p,"iso_y","none")
        parameter!(p,"iso_z","none")
    else
        iso_pts0=GridVisualize.marching_triangles(pts,tris,f,levels)
        niso_pts=length(iso_pts0)
        iso_pts=vcat(reshape(reinterpret(Float32,iso_pts0),(2,niso_pts)))

        iso_x=Float32[]
        iso_y=Float32[]
        iso_z=Float32[]
        ipt=0
        for i=1:niso_pts//2
            push!(iso_x,iso_pts[1,ipt+1])
            push!(iso_x,iso_pts[1,ipt+2])
            push!(iso_x,NaN32)
            push!(iso_y,iso_pts[2,ipt+1])
            push!(iso_y,iso_pts[2,ipt+2])
            push!(iso_y,NaN32)
            push!(iso_z,zval)
            push!(iso_z,zval)
            push!(iso_z,NaN32)
            ipt+=2
        end
        parameter!(p,"iso_x",iso_x)
        parameter!(p,"iso_y",iso_y)
        parameter!(p,"iso_z",iso_z)
    end

end

"""
$(SIGNATURES)

Experimental. Plot heatmap and isolines on rectangular grid defined by X and Y
using Plotly's native contour plot.
"""
function contour!(p::PlutoPlotlyPlot,X,Y,f; kwargs...)
    default_args=(colormap=:viridis, isolines=11,  aspect=1)
    args=merge(p.args,default_args)
    args=merge(args,kwargs)

    isolines=args[:isolines]
    colormap=args[:colormap]

    limits=args.limits

    p.jsdict=Dict{String,Any}("cmdcount" => 0)
    command!(p,"contour")
    parameter!(p,"x",collect(X))
    parameter!(p,"y",collect(Y))
    parameter!(p,"z",vec(f))
    stops=collect(0:0.01:1)
    rgb=reinterpret(Float64,get(colorschemes[colormap],stops,(0,1)))
    rgb=UInt8.(floor.(rgb*255))
    parameter!(p,"cstops",stops)
    parameter!(p,"colors",rgb)
    parameter!(p,"aspect",[1.0,args[:aspect],1.0])

    axisargs!(p,args)

    if limits[1]>limits[2]
        if isa(isolines,Number)
            (fmin,fmax)=extrema(f)
        else
            (fmin,fmax)=extrema(isolines)
        end
    else
        (fmin,fmax)=limits
    end
    if isa(isolines,Number)
        niso=max(2,isolines)
    else
        niso=length(isolines)
    end

    parameter!(p,"costart",fmin)
    parameter!(p,"coend",fmax)
    parameter!(p,"cosize",(fmax-fmin)/(niso-1))
    p
end

###############################################################
# Experimental part
# It turns out that plotly is significantly slower
# than VTK for 3D data.

"""
$(SIGNATURES)

Experimental. Plot piecewise linear function on  triangular grid given by points and triangles
as matrices
"""
function triplot!(p::PlutoPlotlyPlot,pts, tris,f; kwargs...)
    default_args=(colormap=:viridis, isolines=11,  aspect=1)
    args=merge(p.args,default_args)
    args=merge(args,kwargs)
    axisargs!(p,args)

    p.jsdict=Dict{String,Any}("cmdcount" => 0)

    command!(p,"triplot")

    parameter!(p,"x",pts[1,:])
    parameter!(p,"y",pts[2,:])
    parameter!(p,"z",f)

    parameter!(p,"i",tris[1,:].-1)
    parameter!(p,"j",tris[2,:].-1)
    parameter!(p,"k",tris[3,:].-1)
end

function triupdate!(p::PlutoPlotlyPlot,pts,tris,f; kwargs...)
    default_args=(colormap=:viridis, isolines=11,  aspect=1)
    args=merge(p.args,default_args)
    args=merge(args,kwargs)
    axisargs!(p,args)

    p.jsdict=Dict{String,Any}("cmdcount" => 0)
    command!(p,"triupdate")

    # make 3D points from 2D points by adding function value as
    # z coordinate
    parameter!(p,"x",pts[1,:])
    parameter!(p,"y",pts[2,:])
    parameter!(p,"z",f)

    parameter!(p,"i",tris[1,:].-1)
    parameter!(p,"j",tris[2,:].-1)
    parameter!(p,"k",tris[3,:].-1)
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    EstimationSolution

A solution to a spatial estimation problem.
"""
struct EstimationSolution{D<:AbstractDomain}
  domain::D
  mean::Dict{Symbol,<:AbstractVector}
  variance::Dict{Symbol,<:AbstractVector}
end

EstimationSolution(domain, mean, variance) =
  EstimationSolution{typeof(domain)}(domain, mean, variance)

"""
    getindex(solution, var)

Return estimation solution for specific variable `var`
as a named tuple (mean=m, variance=v).
"""
function Base.getindex(solution::EstimationSolution, var::Symbol)
  (mean=solution.mean[var], variance=solution.variance[var])
end

function Base.getindex(solution::EstimationSolution{<:RegularGrid}, var::Symbol)
  sz = size(solution.domain)
  M = reshape(solution.mean[var], sz)
  V = reshape(solution.variance[var], sz)
  (mean=M, variance=V)
end

# ------------
# IO methods
# ------------
function Base.show(io::IO, solution::EstimationSolution)
  dim = ndims(solution.domain)
  print(io, "$(dim)D EstimationSolution")
end

function Base.show(io::IO, ::MIME"text/plain", solution::EstimationSolution)
  println(io, solution)
  println(io, "  domain: ", solution.domain)
  print(  io, "  variables: ", join(keys(solution.mean), ", ", " and "))
end

ext: jl


text:
"""
    NeuralStatistician(IE, [clength::Int, pc], enc_c, cond_z, enc_z, dec)

Neural Statistician model from https://arxiv.org/abs/1606.02185
paper Towards a Neural Statistician.

Acts on bag data (sets of instances).

# Arguments
* `IE`: instance encoder (trainable neural network)
* `pc`: MvNormal prior p(c)
* `clength`: dimension of context prior MvNormal distribution
* `enc_c`: context encoder q(c|D)
* `cond_z`: conditional p(z|c)
* `enc_z`: instance encoder q(z|x,c)
* `dec`: decoder p(x|z)

# Example
Create a Neural Statistician model:
```julia-repl
julia> idim, vdim, cdim, zdim = 5, 3, 2, 4
julia> instance_enc = Chain(Dense(idim,15,swish),Dense(15,vdim))
julia> enc_c = SplitLayer(vdim,[cdim,1])
julia> enc_c_dist = ConditionalMvNormal(enc_c)
julia> cond_z = SplitLayer(cdim,[zdim,1])
julia> cond_z_dist = ConditionalMvNormal(cond_z)
julia> enc_z = SplitLayer(cdim+vdim,[zdim,1])
julia> enc_z_dist = ConditionalMvNormal(enc_z)
julia> dec = SplitLayer(zdim,[idim,1])
julia> dec_dist = ConditionalMvNormal(dec)

julia> model = NeuralStatistician(instance_enc, cdim, enc_c_dist, cond_z_dist, enc_z_dist, dec_dist)

julia> bag = randn(idim,12)
julia> loss(x) = -elbo(model,x)
julia> loss(bag)
10430.707315113537
```
"""
struct NeuralStatistician{IE,pc <: ContinuousMultivariateDistribution,qc <: ConditionalMvNormal,pz <: ConditionalMvNormal,qz <: ConditionalMvNormal,D <: ConditionalMvNormal} # <: AbstractNS
    instance_encoder::IE
    prior_c::pc
    encoder_c::qc
    conditional_z::pz
    encoder_z::qz
    decoder::D
end

Flux.@functor NeuralStatistician

function Flux.trainable(m::NeuralStatistician)
    (instance_encoder = m.instance_encoder, encoder_c = m.encoder_c, conditional_z = m.conditional_z, encoder_z = m.encoder_z, decoder = m.decoder)
end

function NeuralStatistician(IE, clength::Int, enc_c::ConditionalMvNormal, cond_z::ConditionalMvNormal, enc_z::ConditionalMvNormal, dec::ConditionalMvNormal)
    W = first(Flux.params(enc_c))
    μ_c = fill!(similar(W, clength), 0)
    σ_c = fill!(similar(W, clength), 1)
    prior_c = DistributionsAD.TuringMvNormal(μ_c, σ_c)
    NeuralStatistician(IE, prior_c, enc_c, cond_z, enc_z, dec)
end

"""
    elbo(m::NeuralStatistician,x::AbstractArray; β1=1.0, β2=1.0)

Neural Statistician log-likelihood lower bound.

For a Neural Statistician model, simply create a loss
function as

    `loss(x) = -elbo(model,x)`

where `model` is a NeuralStatistician type.

The β terms scale the KLDs:
* β1: KL[q(c|D) || p(c)]
* β2: KL[q(z|c,x) || p(z|c)]
"""
function elbo(m::NeuralStatistician, x::AbstractArray;β1=1.0,β2=1.0)
    # instance network
    v = m.instance_encoder(x)
    p = mean(v, dims=2)

    # sample latent for context
    c = rand(m.encoder_c, p)

    # sample latent for instances
    h = hcat([vcat(v[1:end,i], c) for i in 1:size(v, 2)]...)
    z = rand(m.encoder_z, h)

    # 3 terms - likelihood, kl1, kl2
    llh = mean(logpdf(m.decoder, x, z))
    kld1 = mean(kl_divergence(condition(m.encoder_c, v), m.prior_c))
    kld2 = mean(kl_divergence(condition(m.encoder_z, h), condition(m.conditional_z, c)))
    llh - β1 * kld1 - β2 * kld2
end

function Base.show(io::IO, m::NeuralStatistician)
    IE = repr(m.instance_encoder)
    IE = sizeof(IE) > 70 ? "($(IE[1:70 - 3])...)" : IE
    pc = repr(m.prior_c)
    pc = sizeof(pc) > 70 ? "($(pc[1:70 - 3])...)" : pc
    E_c = repr(m.encoder_c)
    E_c = sizeof(E_c) > 70 ? "($(E_c[1:70 - 3])...)" : E_c
    cond_z = repr(m.conditional_z)
    cond_z = sizeof(cond_z) > 70 ? "($(cond_z[1:70 - 3])...)" : E_z
    E_z = repr(m.encoder_z)
    E_z = sizeof(E_z) > 70 ? "($(E_z[1:70 - 3])...)" : E_z
    D = repr(m.decoder)
    D = sizeof(D) > 70 ? "($(D[1:70 - 3])...)" : D

    msg = """$(nameof(typeof(m))):
     instance_encoder = $(IE)
     prior_c          = $(pc)
     encoder_c        = $(E_c)
     conditional_z    = $(cond_z)
     encoder_z        = $(E_z)
     decoder          = $(D)
    """
    print(io, msg)
end

# needs to extend kl_divergence from IPMeasures
# (m::KLDivergence)(p::ConditionalDists.BMN, q::ConditionalDists.BMN) = IPMeasures._kld_gaussian(p,q)


text:
o1 = OrderedPair(0,1)
o2 = OrderedPair(0.0,1.0)
o3 = OrderedPair(1,2)
o4 = OrderedPair(0.5,0.8)
o5 = OrderedPair(-2.333,1.66)

@test o1.left == 0
@test o1.right == 1
@test o2.left == 0
@test o2.right == 1
@test o4.left == 0.5
@test o4.right == 0.8
@test o5.left == -2.333
@test o5.right == 1.66

@test string(o1) == "(0,1)"
@test string(o4) == "(0.5,0.8)"
@test string(o5) == "(-2.333,1.66)"

@test o1 == o2
@test o1 == o1
@test o1 != o3
@test o2 != o4

@test o1 < o3
@test o3 > o2
@test o5 < o1
@test o5 < o4


text:
module NamespacedEnums

export @nsenum

"""
    @nsenum [T=T] EnumName[::BaseType] value1[=x] value2[=y]

Create a namespaced enum.
This enum works similarly to a regular `Enum`, but its values are hidden behind a module.
To access the values, use `EnumName.value1`.
To access the enum itself, use `EnumName.T`.

```julia-repl
julia> using NamespacedEnums

julia> @nsenum Fruit apple=1 orange=2 kiwi=3;

julia> f(x::Fruit.T) = "I'm a Fruit with value: \$(Int(x))";

julia> f(apple)
ERROR: UndefVarError: apple not defined

julia> f(Fruit.apple)
"I'm a Fruit with value: 1"
```
"""
macro nsenum(exs...)
    return esc(nsenum(exs...))
end

function nsenum(exs...)
    T, name_idx = exs[1] isa Expr && exs[1].head === :(=) ? (exs[1].args[2], 2) : (:T, 1)
    name = exs[name_idx]
    if exs[name_idx] isa Expr
        T = Expr(:(::), T, exs[name_idx].args[2])
        name = name.args[1]
    end
    vals = exs[name_idx+1:end]
    M = gensym(name)
    return Expr(
        :block,
        Expr(:toplevel, :(module $M @enum $T $(vals...) end)),
        :(const $name = $M),
    )
end

end


text:
using Parameters
import Base: length

abstract type AbstractSpecies end
export AbstractSpecies

struct EmptySpecies <: AbstractSpecies end

@with_kw struct Species{T<:AbstractThermo,N<:AbstractDiffusivity} <: AbstractSpecies
    name::String
    index::Integer
    inchi::String = ""
    smiles::String = ""
    adjlist::String = ""
    thermo::T
    atomnums::Dict{String,Int64} = Dict()
    bondnum::Int64=-1
    diffusion::N = EmptyDiffusivity()
    radius::Float64 = 0.0
    radicalelectrons::Int64 = -100
    molecularweight::Float64 = 0.0
end
export Species

length(r::T) where {T<:AbstractSpecies}= 1
export length


text:
module Exponential

using ..Auxiliary, ..TensorTrain

export trigevalmask, trigeval, trigrefmask, trigrefmask2, trigdiffmask, trigdiff, trigdec, trigdeceval!
export cosfactor, cosdec

function trigevalmask(t::Vector{T}, ν::Vector{T}) where {T<:AbstractFloat}
	n = length(t)
	r = length(ν)
	if r == 0
		throw(ArgumentError("the number of frequencies should be positive"))
	end
	V = zeros(T, n, 2*r)
	for α ∈ 1:r
		for i ∈ 1:n
			V[i,2*α-1] = cos(π*ν[α]*t[i])
			V[i,2*α  ] = sin(π*ν[α]*t[i])
		end
	end
	return V
end

function trigeval(t::Vector{T}, ν::Vector{T}, c::Vector{T}) where {T<:AbstractFloat}
	n = length(t)
	r = length(ν)
	if r == 0
		throw(ArgumentError("the number of frequencies should be positive"))
	end
	if length(c) ≠ 2*r
		throw(ArgumentError("ν and c are incompatible in size"))
	end
	u = zeros(T, length(t))
	for i ∈ 1:n
		for α ∈ 1:r
			u[i] += c[2*α-1]*cos(π*ν[α]*t[i])+c[2*α]*sin(π*ν[α]*t[i])
		end
	end
	return u
end

function trigrefmask(η::T, ν::Vector{T}) where {T<:AbstractFloat}
	r = length(ν)
	if r == 0
		throw(ArgumentError("the number of frequencies should be positive"))
	end
	W = zeros(T, 2*r, 2*r)
	for α ∈ 1:r
		W[2*α-1,2*α-1] =  cos(π*ν[α]*η)
		W[2*α,  2*α-1] = -sin(π*ν[α]*η)
		W[2*α-1,2*α  ] =  sin(π*ν[α]*η)
		W[2*α,  2*α  ] =  cos(π*ν[α]*η)
	end
	return W
end

function trigrefmask2(ν::Vector{T}) where {T<:AbstractFloat}
	r = length(ν)
	if r == 0
		throw(ArgumentError("the number of frequencies should be positive"))
	end
	W1 = trigrefmask(convert(T, -1/2), ν)
	W2 = trigrefmask(convert(T, 1/2), ν)
	W = [W1[:] W2[:]]
	W = reshape(W, 2*r, 2*r, 2)
	return W
end

function trigdiffmask(ν::Vector{T}) where {T<:AbstractFloat}
	r = length(ν)
	if r == 0
		throw(ArgumentError("the number of frequencies should be positive"))
	end
	W = zeros(T, 2*r, 2*r)
	for α ∈ 1:r
		W[2*α-1,2*α  ] =  π*ν[α]
		W[2*α,  2*α-1] = -π*ν[α]
	end
	return W
end

function trigdiff(ν::Vector{T}, c::Vector{T}) where {T<:AbstractFloat}
	r = length(ν)
	if length(c) ≠ 2*r
		throw(ArgumentError("ν and c are incompatible in size"))
	end
	d = zeros(T, 2*r)
	for α ∈ 1:r
		d[2*α-1] =  π*ν[α]*c[2*α]
		d[2*α]   = -π*ν[α]*c[2*α-1]
	end
	return d
end

function cosfactor(τ::Vector{T}, c::T) where {T<:AbstractFloat}
	n = length(τ)
	if n == 0
		throw(ArgumentError("τ should be nonempty"))
	end
	U = zeros(T, 2, n, 2)
	for i ∈ 1:n
		α = c*cos(τ[i])
		β = c*sin(τ[i])
		U[1,i,1] =  α
		U[2,i,1] =  β
		U[1,i,2] = -β
		U[2,i,2] =  α
	end
	return U
end

function cosdec(τ::Vector{Vector{T}}, c::Vector{T}) where {T<:AbstractFloat}
	L = length(τ)
	if length(c) ≠ L
		throw(ArgumentError("τ and c should be of the same length"))
	end
	n = [ size(τ[ℓ]) for ℓ ∈ 1:L ]
	if any(τ .== 0)
		throw(ArgumentError("all elements of τ should be nonempty"))
	end
	U = [ cosfactor(τ[ℓ], c[ℓ]) for ℓ ∈ 1:L ]
	U[1] = U[1][1:1,:,:]
	U[L] = U[L][:,:,1:1]
	U = dec(U)
	return U
end

function trigdec(ν::Vector{T}, c::Vector{T}, L::Int; major::String="last") where {T<:AbstractFloat}
	if major ∉ ("first", "last")
		throw(ArgumentError("major should be either \"last\" (default) or \"first\""))
	end
	r = length(ν)
	if length(c) ≠ 2*r
		throw(ArgumentError("ν and c are incompatible in size"))
	end
	U = dec(T, 1, L+1)
	U[1] = factor(Matrix{T}(c'), 1, [])
	for ℓ ∈ 1:L
		U[ℓ+1] = permutedims(trigrefmask2(ν/2^(ℓ-1)), (2,3,1))
	end
	(major == "last") && decreverse!(U)
	return U
end

function trigdeceval!(U::Dec{T,N}, t::Vector{T}, ν::Vector{T}; major::String="last") where {T<:AbstractFloat,N}
	if major ∉ ("first", "last")
		throw(ArgumentError("major should be either \"last\" (default) or \"first\""))
	end
	if N ≠ 3
		throw(ArgumentError("only one mode index per factor is currently supported"))
	end
	n = length(t)
	r = length(ν)
	L = declength(U)
	p = decrank(U)
	p = (major == "last") ? p[1] : p[L+1]
	if p ≠ 2*r
		throw(ArgumentError("the ", (major == "last") ? "first" : "last", " rank of U and the size of ν are incompatible"))
	end
	W = trigevalmask(t, ν)
	if major == "last"
		W = reshape(W, 1, n, 2*r)
		W = Array{T,N}(W)
		decpushfirst!(U, W)
	else
		W = reshape(transpose(W), 2*r, n, 1)
		W = Array{T,N}(W)
		decpush!(U, W)
	end
	return U
end

end


text:
using Test
using GlobalSensitivityAnalysis

ATOL_sample = 1e-8
ATOL_sobol = 1e-9 # sobol method indices
ATOL_delta = 1e-3 # delta method indices
ATOL_CI = 1e-2

# see: https://www.nature.com/articles/sdata2018187#Sec5
include(joinpath(@__DIR__, "../src/quantile_matching/WRS.jl"))

@testset "Test Utils" begin
    include("test_utils.jl")
end

@testset "Test Sobol Method" begin

    @testset "Compare to SALib.py" begin
        include("sobol/test_sobol_sobolsequence.jl")
        include("sobol/test_sobol_ishigami.jl")
        include("sobol/test_sobol_borehole.jl")
        include("sobol/test_sobol_ci.jl")
    end

    @testset "Unit Testing" begin
        include("sobol/test_sobol_unit.jl")
    end
end

@testset "Test Delta Method" begin

    @testset "Compare to SALib.py" begin
        include("delta/test_delta_ishigami.jl")
        include("delta/test_delta_borehole.jl")
    end

    @testset "Unit Testing" begin
        include("delta/test_delta_unit.jl")
    end
end

max_stars_repo_name: lrennels/SAJulia

max_issues_repo_name: lrennels/SAJulia


text:
module HTSLib

export eachrecord

using Printf: @printf

include("htslib/htslib.jl")

"""
The version of htslib.
"""
const HTSLIB_VERSION = VersionNumber(unsafe_string(htslib.hts_version()))

# Record view
# -----------

"""
A light-weight wrapper of records.
"""
struct RecordView{T}
    ptr::Ptr{T}  # pointer to a htslib-allocated record struct
end

function Base.show(io::IO, view::RecordView{htslib.bam1_t})
    ptr = view.ptr
    bam = unsafe_load(ptr)
    @printf(io, "%s(<%d:%d@%p>)", summary(view), bam.core.tid, bam.core.pos, ptr)
end

@inline function Base.getproperty(view::RecordView{htslib.bam1_t}, name::Symbol)
    bam = unsafe_load(getfield(view, :ptr))
    if name == :pos
        return bam.core.pos
    elseif name == :tid
        return bam.core.tid
    elseif name == :bin
        return bam.core.bin
    elseif name == :qual
        return bam.core.qual
    elseif name == :l_extranul
        return bam.core.l_extranul
    elseif name == :flag
        return bam.core.flag
    elseif name == :l_qname
        return bam.core.l_qname
    elseif name == :n_cigar
        return bam.core.n_cigar
    elseif name == :l_qseq
        return bam.core.l_qseq
    elseif name == :mtid
        return bam.core.mtid
    elseif name == :mpos
        return bam.core.mpos
    elseif name == :isize
        return bam.core.isize
    else
        return getfield(view, name)
    end
end

# Iterator
# --------

mutable struct RecordIterator
    file::Ptr{htslib.htsFile}
    header::Ptr{htslib.sam_hdr_t}
    index::Ptr{htslib.hts_idx_t}
    iterator::Ptr{htslib.hts_itr_t}
    record::Ptr{htslib.bam1_t}

    function RecordIterator(file, header, index, iterator, record)
        iter = new(file, header, index, iterator, record)
        finalizer(iter) do iter
            if iter.record != C_NULL
                htslib.bam_destroy1(iter.record)
                iter.record = C_NULL
            end
            if iter.iterator != C_NULL
                htslib.hts_itr_destroy(iter.iterator)
                iter.iterator = C_NULL
            end
            if iter.index != C_NULL
                htslib.hts_idx_destroy(iter.index)
                iter.index = C_NULL
            end
            if iter.header != C_NULL
                htslib.sam_hdr_destroy(iter.header)
                iter.header = C_NULL
            end
            if iter.file != C_NULL
                res = htslib.hts_close(iter.file)
                @assert res == 0
                iter.file = C_NULL
            end
        end
        return iter
    end
end

function eachrecord(filepath::AbstractString)
    # TODO: more careful memory management
    file = htslib.hts_open(filepath, "r")
    file == C_NULL && error("failed to open $(filepath)")
    header = htslib.sam_hdr_read(file)
    header == C_NULL && error("failed to read the header of $(filepath)")
    record = htslib.bam_init1()
    record == C_NULL && error("failed to allocate a new record")
    return RecordIterator(file, header, C_NULL, C_NULL, record)
end

function eachrecord(filepath::AbstractString, region::AbstractString)
    # TODO: more careful memory management
    file = htslib.hts_open(filepath, "r")
    file == C_NULL && error("failed to open $(filepath)")
    header = htslib.sam_hdr_read(file)
    header == C_NULL && error("failed to read the header of $(filepath)")
    index_format = htslib.HTS_FMT_BAI
    index = htslib.hts_idx_load(filepath, index_format)
    index == C_NULL && error("failed to load the index file for $(filepath)")
    iterator = htslib.sam_itr_querys(index, header, region)
    iterator == C_NULL && error("failed to parse the region query string: $(repr(region))")
    record = htslib.bam_init1()
    record == C_NULL && error("failed to allocate a new record")
    return RecordIterator(file, header, index, iterator, record)
end

function Base.iterate(iter::RecordIterator, ::Nothing = nothing)
    @assert iter.index == C_NULL && iter.iterator == C_NULL ||
            iter.index != C_NULL && iter.iterator != C_NULL
    if iter.index == C_NULL
        res = htslib.sam_read1(iter.file, iter.header, iter.record)
    else
        res = htslib.sam_itr_next(iter.file, iter.iterator, iter.record)
    end
    @assert res ≥ -1
    if res ≥ 0
        # one or more records
        return RecordView(iter.record), nothing
    else
        # no records
        finalize(iter)  # early finalization
        return nothing
    end
end

end # module


text:
"""
Module that provides concrete implementations of the [`AbstractDVec`](@ref) interface.

- [`DVec`](@ref): basic [`AbstractDVec`](@ref)
- [`InitiatorDVec`](@ref): allows storing information about initiator status

See [`Interfaces`](@ref).
"""
module DictVectors

using Random
using LinearAlgebra
import SplittablesBase

using ..Interfaces
import ..Interfaces: deposit!, storage, StochasticStyle, default_style, freeze

export zero!, add!, deposit!, storage, walkernumber
export DVec, InitiatorDVec
export AbstractProjector, NormProjector, Norm2Projector, UniformProjector, Norm1ProjectorPPop

# The idea is to do linear algebra with data structures that are not
# subtyped to AbstractVector, much in the spirit of KrylovKit.jl.
# In particular we provide concrete data structures with the aim of being
# suitable for use with KrylovKit. From the manual:

# KrylovKit does not assume that the vectors involved in the problem are actual
# subtypes of AbstractVector. Any Julia object that behaves as a vector is
# supported, so in particular higher-dimensional arrays or any custom user type
# that supports the following functions (with v and w two instances of this type
# and α, β scalars (i.e. Number)):
#
# Base.eltype(v): the scalar type (i.e. <:Number) of the data in v
# Base.similar(v, [T::Type<:Number]): a way to construct additional similar vectors, possibly with a different scalar type T.
# Base.copyto!(w, v): copy the contents of v to a preallocated vector w
# LinearAlgebra.mul!(w, v, α): out of place scalar multiplication; multiply vector v with scalar α and store the result in w
# LinearAlgebra.rmul!(v, α): in-place scalar multiplication of v with α; in particular with α = false, v is initialized with all zeros
# LinearAlgebra.axpy!(α, v, w): store in w the result of α*v + w
# LinearAlgebra.axpby!(α, v, β, w): store in w the result of α*v + β*w
# LinearAlgebra.dot(v,w): compute the inner product of two vectors
# LinearAlgebra.norm(v): compute the 2-norm of a vector

# It turns out, KrylovKit also needs
# *(v, α::Number)
# fill!(v, α)

include("delegate.jl")
include("abstractdvec.jl")
include("dvec.jl")
include("initiators.jl")

end # module


text:
# Unit test functions
function unit_test_queue()
    q = Queue(DLinkedList{Float64}())
    @test length(q.data) == 0

    return nothing
end

function unit_test_queue_isempty()
    q = Queue(DLinkedList{Float64}())
    res = isempty(q)

    return nothing
end

function unit_test_queue_enqueue()
    q = Queue(DLinkedList{Float64}())

    val1 = 1.0
    val2 = 2.0
    enqueue!(q, val1)
    enqueue!(q, val2)
    @test q.data.head.value == val1
    @test q.data.head.next.value == val2
    @test length(q) == 2

    return nothing
end

function unit_test_queue_dequeue()
    q = Queue(DLinkedList{Float64}())

    val1 = 1.0
    val2 = 2.0
    enqueue!(q, val1)
    enqueue!(q, val2)
    @test length(q) == 2
    @test q.data.head.value == val1
    @test q.data.head.next.value == val2
    @test dequeue!(q) == val1
    @test dequeue!(q) == val2
    @test length(q) == 0

    return nothing
end

function unit_test_queue_enqueuedequeue()
    res = []
    q = Queue(DLinkedList{Float64}())

    val1 = 1.0
    val2 = 2.0
    val3 = 3.0
    val4 = 4.0
    enqueue!(q, val1)
    enqueue!(q, val2)
    @test length(q) == 2
    @test q.data.head.value == val1
    @test q.data.head.next.value == val2
    @test dequeue!(q) == val1
    enqueue!(q, val3)
    @test dequeue!(q) == val2
    enqueue!(q, val4)
    @test dequeue!(q) == val3
    @test dequeue!(q) == val4
    @test length(q) == 0

    return nothing
end

function unit_test_queue_peek()
    res = []
    q = Queue(DLinkedList{Float64}())

    val1 = 1.0
    val2 = 2.0
    enqueue!(q, val1)
    enqueue!(q, val2)
    @test peek(q) == val1
    @test peek(q) == val1
    @test length(q) == 2

    return nothing
end

function unit_test_queue_length()
    res = []

    q = Queue(DLinkedList{Float64}())
    val = 1.0
    enqueue!(q, val)
    @test length(q) == 1
    enqueue!(q, val)
    @test length(q) == 2

    return nothing
end

function unit_test_queue_size()
    res = []

    q = Queue(DLinkedList{Float64}())
    val = 1.0
    enqueue!(q, val)
    @test size(q) == 1
    enqueue!(q, val)
    @test size(q) == 2

    return nothing
end

# Run the unit tests
@testset "Queue" begin
    unit_test_queue()
    unit_test_queue_isempty()
    unit_test_queue_enqueue()
    unit_test_queue_dequeue()
    unit_test_queue_enqueuedequeue()
    unit_test_queue_peek()
    unit_test_queue_length()
    unit_test_queue_size()
end


text:
###########################################
## DE
###########################################
function reproduction(status, parameters::DE, problem)
    @assert !isempty(status.population)

    N = parameters.N
    D = length(get_position(status.population[1]))

    strategy = parameters.strategy
    xBest = get_position(status.best_sol)
    population = status.population
    F = parameters.F
    CR = parameters.CR

    X = zeros(N,D)

    for i in 1:N
        x = get_position(population[i])
        u = DE_mutation(population, F, strategy, 1)
        v = DE_crossover(x, u, CR)
        evo_boundary_repairer!(v, xBest, problem.bounds)
        X[i,:] = v
    end

    X
end

###########################################
## ECA
###########################################
function reproduction(status, parameters::ECA, problem)
    @assert !isempty(status.population)

    N = parameters.N
    D = length(get_position(status.population[1]))

    X = zeros(N,D)

    for i in 1:N
        X[i,:] = ECA_operator(status.population,
                              parameters.K,
                              parameters.η_max;
                              i=i,
                             bounds = problem.bounds)
    end

    X
end

###########################################
## generic GA reproduction
###########################################
function reproduction(status, parameters::AbstractNSGA, problem)
    @assert !isempty(status.population)

    I = randperm(parameters.N)
    J = randperm(parameters.N)
    Q = zeros(parameters.N, size(problem.bounds, 2))
    for i = 1:2:parameters.N

        pa = tournament_selection(status.population, I[i])
        pb = tournament_selection(status.population, J[i])

        c1, c2 = GA_reproduction(get_position(pa),
                                 get_position(pb),
                                 problem.bounds;
                                 η_cr = parameters.η_cr,
                                 p_cr = parameters.p_cr,
                                 η_m = parameters.η_m,
                                 p_m = parameters.p_m)
        Q[i,:] = c1
        Q[i+1,:] = c2
    end

    Q
end

###########################################
## NSGA3 reproduction
###########################################
function reproduction(status, parameters::NSGA3, problem)
    @assert !isempty(status.population)

    I = randperm(parameters.N)
    Q = zeros(parameters.N, size(problem.bounds, 2))
    for i = 1:parameters.N ÷ 2

        pa = status.population[I[2i-1]]
        pb = status.population[I[2i]]

        c1, c2 = GA_reproduction(get_position(pa),
                                 get_position(pb),
                                 problem.bounds;
                                 η_cr = parameters.η_cr,
                                 p_cr = parameters.p_cr,
                                 η_m = parameters.η_m,
                                 p_m = parameters.p_m)
        Q[i,:] = c1
        Q[i+1,:] = c2
    end

    Q
end


text:
#Function for calcualting Jmax from Nitrogen per leaf area (Nₐ)
function Calc_Jmax(Nₐ::T,a::T,b::T,b_Jmax::T) where {T<:Float64}
    return max(b_Jmax*(a*Nₐ+b),0.0)
end

Calc_LAI(model::CCPHStruct) = model.treesize.Wf/model.treepar.LMA*model.treesize.N

#Calculate the irradiance incident on a leaf at irradiance I
Calc_Iᵢ(I::Float64,model::CCPHStruct) = I*model.treepar.k/(1-model.treepar.m)

#calcualte total conductance
Calc_gₜ(gₛ::Float64,model::CCPHStruct) = gₛ*model.treepar.r_gₛ

#Calculate per sapwood mass nitrogen concentration
Calc_Nₘ_w(Nₘ_f::Float64,model::CCPHStruct) = model.treepar.rW*Nₘ_f
#Calcualte per fine roots mass nitrogen concentration
Calc_Nₘ_r(Nₘ_f::Float64,model::CCPHStruct) = model.treepar.rR*Nₘ_f
#Calcualte per leaf area nitrogen concentration
Calc_Nₐ(Nₘ_f::Float64,model::CCPHStruct) = model.treepar.LMA*Nₘ_f

#Calculate total maintenance respiration
function Calc_Rₘ(Nₘ_f::T,Nₘ_w::T,Nₘ_r::T,Wf::T,Ww::T,Wr::T,model::CCPHStruct) where {T<:Float64}
    Rₘ = model.treepar.rₘ*(Nₘ_f*Wf+Nₘ_w*Ww+Nₘ_r*Wr)
    return Rₘ
end

Calc_NPP(P::T,Rₘ::T,model::CCPHStruct) where {T<:Float64} = model.treepar.y*(P-Rₘ)

#Foliage and fine root time depednent senescence
Calc_S_fr(Wf::T,Wr::T,model::CCPHStruct) where {T<:Float64} = Wf/model.treepar.Tf+Wr/model.treepar.Tr

#Tree height time derivative
function Calc_dH(NPP::T,S::T,Wf::T,Ww::T,Wr::T,H::T,Hs::T,model::CCPHStruct) where {T<:Float64}
    dH = (NPP-S)/( model.treepar.β₁*Ww/(model.treepar.β₁*H+model.treepar.β₂*Hs)+
    model.treepar.z*Wf/(H-Hs)+
    model.treepar.z*Ww/(H-Hs)+
    model.treepar.z*Wr/(H-Hs))

    return dH
end

#Function for calcualting fine root to sapwood area ratio based on carbon and nitrogen constraints
function Find_αr(Nₘ_f::T,Nₘ_w::T,Nₘ_r::T,P::T,model::CCPHStruct) where {T<:Float64}
    γC0 = model.treepar.y*(P-model.treepar.rₘ*(Nₘ_f*model.treesize.Wf+Nₘ_w*model.treesize.Ww))-model.treesize.Wf/model.treepar.Tf
    γC1 = model.treepar.y*model.treepar.rₘ*Nₘ_r*model.treesize.As+model.treesize.As/model.treepar.Tr
    γC2 = model.treepar.β₁*model.treesize.Ww/(model.treepar.β₁*model.treesize.H+model.treepar.β₂*model.treesize.Hs)+
        model.treepar.z*model.treesize.Wf/(model.treesize.H-model.treesize.Hs)+model.treepar.z*model.treesize.Ww/(model.treesize.H-model.treesize.Hs)
    γC3 = model.treepar.z*model.treesize.As/(model.treesize.H-model.treesize.Hs)

    γN0 = Nₘ_f*model.treesize.Wf/model.treepar.Tf
    γN1 = Nₘ_r*model.treesize.As/model.treepar.Tr
    γN2 = Nₘ_w*model.treepar.β₁*model.treesize.Ww/(model.treepar.β₁*model.treesize.H+model.treepar.β₂*model.treesize.Hs)+
    Nₘ_f*model.treepar.z*model.treesize.Wf/(model.treesize.H-model.treesize.Hs)+Nₘ_w*model.treepar.z*model.treesize.Ww/(model.treesize.H-model.treesize.Hs)
    γN3 = Nₘ_r*model.treepar.z*model.treesize.As/(model.treesize.H-model.treesize.Hs)

    γU0 = model.treepar.Nₛ*model.treesize.As*model.treepar.Kr
    γU1 = model.treepar.Kr
    γU2 = model.treesize.As*model.treesize.N

    θ₀ = (γC0*γN2+γN0*γC2)*γU1
    θ₁ = (γC0*γN3+γN1*γC2+γN0*γC3-γC1*γN2)*γU1+(γC0*γN2+γN0*γC2)*γU2-γC2*γU0
    θ₂ = (γC0*γN3+γN1*γC2+γN0*γC3-γC1*γN2)*γU2+(γN1*γC3-γC1*γN3)*γU1-γC3*γU0
    θ₃ = (γN1*γC3-γC1*γN3)*γU2

    if isnan(θ₀)||isnan(θ₁)||isnan(θ₂)||isnan(θ₃)
        error("Polynomial coefficients contains NaN")
    end
    if isinf(θ₀)||isinf(θ₁)||isinf(θ₂)||isinf(θ₃)
        error("Polynomial coefficients contains ∞")
    end

    αr_vec = solvecubic(θ₃, θ₂, θ₁, θ₀)
    filter!(x->isnan(x)==false&&isinf(x)==false&&isreal(x)&&x>0.0,αr_vec)

    isempty(αr_vec)==false||error("No feasible αr")

    return αr_vec
end

#Calcualte gain (nett carbon gain minus fine root growth, performance measure)
function Gain_fun(αr::T,K_cost::T,Nₘ_f::T,Nₘ_w::T,Nₘ_r::T,P::T,model::CCPHStruct) where {T<:Float64}
    #Fine root mass
    Wr = model.treesize.As*αr

    #Calculate total maintenance respiration
    R_m = Calc_Rₘ(Nₘ_f,Nₘ_w,Nₘ_r,model.treesize.Wf,model.treesize.Ww,Wr,model)

    #Net primary production
    NPP = Calc_NPP(P,R_m,model)

    #Foliage and fine root time depednent senescence
    S = Calc_S_fr(model.treesize.Wf,Wr,model)

    #Tree height time derivative
    dH = Calc_dH(NPP,S,model.treesize.Wf,model.treesize.Ww,Wr,model.treesize.H,model.treesize.Hs,model)

    #Fine root growth
    Gr = model.treepar.z*Wr/(model.treesize.H-model.treesize.Hs)*dH+Wr/model.treepar.Tr

    return (NPP-Gr)*K_cost^model.hydPar.i
end

#Run the Coupled Canopy Photosynthesis and Hydraulics model
function CCPH_run(gₛ::T,Nₘ_f::T,growthlength::T,model::CCPHStruct) where {T<:Float64}
    #Calculate per sapwood mass nitrogen concentration
    Nₘ_w = Calc_Nₘ_w(Nₘ_f,model)
    #Calcualte per fine roots mass nitrogen concentration
    Nₘ_r = Calc_Nₘ_r(Nₘ_f,model)
    #Calcualte per leaf area nitrogen concentration
    Nₐ = Calc_Nₐ(Nₘ_f,model)
    #Calculate Jmax
    Jmax = Calc_Jmax(Nₐ,model.treepar.a_Jmax,model.treepar.b_Jmax,model.photopar.b_Jmax)
    #Irradiance incident on a leaf at canopy top
    Iᵢ = Calc_Iᵢ(model.env.I₀,model)
    #Calculate LAI
    LAI = Calc_LAI(model)
    #calcualte total conductance
    gₜ = Calc_gₜ(gₛ,model)
    #Calculate per tree carbon assimilation
    P = GPP(gₜ,Iᵢ,Jmax,LAI,growthlength,model)
    #---Calculate cost factor of hydraulic failure---
    K_cost, Kₓₗ, ψ_c = Calc_K_cost(gₛ,model)
    #---Calculate posible fine root to sapwood area ratio, αr---
    αr_vec = Find_αr(Nₘ_f,Nₘ_w,Nₘ_r,P,model)
    #Calcualte performance measure
    Gain_vec = [Gain_fun(αr,K_cost,Nₘ_f,Nₘ_w,Nₘ_r,P,model) for αr in αr_vec]

    #Find feasible αr which maximizes performance
    filter!(x->isnan(x)==false&&isnan(x)==false&&isreal(x),Gain_vec)
    isempty(Gain_vec)==false||error("No feasible Gain")
    αr_val = αr_vec[argmax(Gain_vec)]

    modeloutput = CCPHOutput(P,αr_val,ψ_c,Kₓₗ,K_cost,maximum(Gain_vec))

    return modeloutput
end

function Trait_objective_fun(x::Array{T,1},growthlength::T,model::CCPHStruct) where {T<:Float64}
    gₛ,Nₘ_f = x

    modeloutput = CCPH_run(gₛ,Nₘ_f,growthlength,model)

    return -modeloutput.Gain
end

#Find optimal triats
function CCPHTraitmodel(growthlength::T,model::CCPHStruct;
    gₛ_guess::T=0.02,Nₘ_f_guess::T=0.012,gₛ_lim_lo::T=0.01,gₛ_lim_hi::T=0.5,
    Nₘ_f_lim_lo::T=0.007,Nₘ_f_lim_hi::T=0.05) where {T<:Float64}

    x0 = [gₛ_guess, Nₘ_f_guess]

    lower = [gₛ_lim_lo, Nₘ_f_lim_lo]
    upper = [gₛ_lim_hi, Nₘ_f_lim_hi]

    inner_optimizer = Optim.BFGS(linesearch=Optim.LineSearches.BackTracking())

    opt_trait = Optim.optimize(x->Trait_objective_fun(x,growthlength,model),
    lower,upper,x0,Optim.Fminbox(inner_optimizer))

    Optim.converged(opt_trait)||error("No optimal traits could be found")

    gₛ_opt = opt_trait.minimizer[1]
    Nₘ_f_opt = opt_trait.minimizer[2]

    return gₛ_opt,Nₘ_f_opt
end

#Differential equations describing the growth of the forest
function TreeStandDyn!(dy::Array{T,1},y::Array{T,1},model::CCPHStruct,
    t::T,growthlength::T,gₛ::T,Nₘ_f::T,αr::T) where {T<:Float64}
    Wf,Ww,H,Hs,As,B,N,Wr = y

    #Calcualte per leaf area nitrogen concentration
    Nₐ = Calc_Nₐ(Nₘ_f,model)
    #Calculate Jmax
    Jmax = Calc_Jmax(Nₐ,model.treepar.a_Jmax,model.treepar.b_Jmax,model.photopar.b_Jmax)
    #Irradiance incident on a leaf at canopy top
    Iᵢ = Calc_Iᵢ(model.env.I₀,model)
    #Calculate LAI
    LAI = Calc_LAI(model)
    #calcualte total conductance
    gₜ = Calc_gₜ(gₛ,model)
    #Calculate GPP
    P = GPP(gₜ,Iᵢ,Jmax,LAI,growthlength,model)

    #Calculate per sapwood mass nitrogen concentration
    Nₘ_w = Calc_Nₘ_w(Nₘ_f,model)
    #Calcualte per fine roots mass nitrogen concentration
    Nₘ_r = Calc_Nₘ_r(Nₘ_f,model)

    #Calculate total maintenance respiration
    R_m = Calc_Rₘ(Nₘ_f,Nₘ_w,Nₘ_r,Wf,Ww,Wr,model)

    NPP = Calc_NPP(P,R_m,model)

    #Foliage and fine root time depednent senescence
    S = Calc_S_fr(Wf,Wr,model)

    #Calculate height growth
    dH = Calc_dH(NPP,S,Wf,Ww,Wr,H,Hs,model)

    dN = 0 #Population density is constant

    #---Calcualte leaf performance at crown base---
    Δ_leaf = Calc_Δ_leaf(gₜ,Iᵢ,LAI,growthlength,Nₘ_f,Jmax,model)

    dHs = dH*(Δ_leaf<0) #No crown rise when leaf performance >= 0

    dAs = model.treepar.z*As/(H-Hs)*(dH-dHs)
    dB = model.treepar.z*As/(H-Hs)*dH

    dWf = model.treepar.αf*dAs
    dWr = αr*dAs
    dWw = model.treepar.ρw*dAs*(model.treepar.β₁*H+model.treepar.β₂*Hs)+
    model.treepar.ρw*As*(model.treepar.β₁*dH+model.treepar.β₂*dHs)

    dy[1] = dWf
    dy[2] = dWw
    dy[3] = dH
    dy[4] = dHs
    dy[5] = dAs
    dy[6] = dB
    dy[7] = dN
    dy[8] = dWr
end

#initiate weather parameters from WeatherTS
function Init_weather_par!(i::Integer,model::CCPHStruct,weatherts::WeatherTS,photo_kinetic::PhotoKineticRates)
    growthlength = weatherts.tot_annual_daylight[i] #Length of day light time during the gorwth period of a specific year (s)
    step_length = weatherts.daylight[i]/weatherts.tot_annual_daylight[i] #length of simulation step (proportion of a year)
    model.env.I₀ = weatherts.PAR[i]
    model.env.VPD = weatherts.VPD[i]
    model.env.θₛ = weatherts.θₛ[i]
    model.env.Tₐ = weatherts.temp[i]
    PhotoPar!(model.photopar,photo_kinetic,weatherts.temp[i])
    model.treepar.Xₜ = weatherts.acclimation_fac[i] #Account for acclimation

    return growthlength,step_length
end

#Simulate growth with weather time series
function CCPHStandGrowth!(model::CCPHStruct,weatherts::WeatherTS,
    photo_kinetic::PhotoKineticRates,nstep::Integer;
    K_cost_crit::T=0.12,gₛ_guess::T=0.02,Nₘ_f_guess::T=0.012)::CCPHTS  where {T<:Float64}

    ccphts = CCPHTS() #Init Time series struct
    CCPHTS!(ccphts,model.treesize)

    for i in 1:nstep

        #initiate weather parameters
        growthlength,step_length = Init_weather_par!(i,model,weatherts,photo_kinetic)

        gₛ_crit = Calc_K_costᵢₙᵥ(K_cost_crit,model)

        gₛ_opt,Nₘ_f_opt = CCPHTraitmodel(growthlength,model;
        gₛ_guess=gₛ_guess,Nₘ_f_guess=Nₘ_f_guess,gₛ_lim_hi=gₛ_crit)

        modeloutput = CCPH_run(gₛ_opt,Nₘ_f_opt,growthlength,model)

        CCPHTS!(ccphts,modeloutput,gₛ_opt,Nₘ_f_opt,gₛ_crit)

        y0 = [model.treesize.Wf,model.treesize.Ww,model.treesize.H,
        model.treesize.Hs,model.treesize.As,model.treesize.B,model.treesize.N,modeloutput.αr*model.treesize.As]

        #Simulate one timestep
        prob = DifferentialEquations.ODEProblem((dy::Array{Float64,1},y::Array{Float64,1},P_in::CCPHStruct,t::Float64)->
        TreeStandDyn!(dy,y,P_in,t,growthlength,gₛ_opt,Nₘ_f_opt,modeloutput.αr),y0,(0.0,step_length),model)
        sol = DifferentialEquations.solve(prob)

        #Update current size
        model.treesize.Wf,model.treesize.Ww,model.treesize.H,
        model.treesize.Hs,model.treesize.As,model.treesize.B,model.treesize.N,Wr = sol(step_length)

        #Update time series
        CCPHTS!(ccphts,model.treesize)
    end

    #initiate weather parameters
    growthlength,step_length = Init_weather_par!(nstep+1,model,weatherts,photo_kinetic)

    gₛ_crit = Calc_K_costᵢₙᵥ(K_cost_crit,model)

    gₛ_opt,Nₘ_f_opt = CCPHTraitmodel(growthlength,model;
    gₛ_guess=gₛ_guess,Nₘ_f_guess=Nₘ_f_guess,gₛ_lim_hi=gₛ_crit)

    modeloutput = CCPH_run(gₛ_opt,Nₘ_f_opt,growthlength,model)

    CCPHTS!(ccphts,modeloutput,gₛ_opt,Nₘ_f_opt,gₛ_crit)

    return ccphts
end


text:
##########################################################################
# The IndoorsNetwork is an indoors microcell network inspired by the ITU-R
# InH indoors network. The channel model is simplified; it is just using
# uncorrelated Rayleigh fading. The large-scale fading parameters should
# be close to the InH model though. The number of BSs has also been
# generalized compared to InH. The LoS is independent between BSs, for
# each MS. The same holds for the shadow fading. This is motivated by
# Sec. 3.3.1 in Deliverable 1.1.2 of IST-WINNER II.
#
# The typical transmit power is 21 dBm over 20 MHz (1200 subcarriers).
# This leads to 10*log10(10^(21/10)/1200) = -9.8 dBm for one 15 kHz subcarrier.

##########################################################################
# Network definition
type IndoorsNetwork{MS_t <: PhysicalMS, BS_t <: PhysicalBS, System_t <: System, PropagationEnvironment_t <: PropagationEnvironment} <: PhysicalNetwork
    MSs::Vector{MS_t}
    BSs::Vector{BS_t}

    system::System_t
    num_MSs_per_cell::Int
    propagation_environments::Dict{Symbol,PropagationEnvironment_t} # two keys only, :LoS and :NLoS
    corridor_length::Float64
    corridor_width::Float64
    guard_distance::Float64
    aux_network_params::AuxNetworkParams

    assignment::Assignment
end

# Convenience constructor without network params and assignments
IndoorsNetwork(MSs, BSs, system, num_MSs_per_cell, propagation_environments, corridor_length, corridor_width, guard_distance) =
    IndoorsNetwork(MSs, BSs, system, num_MSs_per_cell, propagation_environments, corridor_length, corridor_width, guard_distance, AuxNetworkParams(), Assignment())

Base.show(io::IO, x::IndoorsNetwork) =
    print(io, "Indoors(I = $(length(x.BSs)), Kc = $(x.num_MSs_per_cell), corridor_length = $(x.corridor_length), corridor_width = $(x.corridor_width), GD = $(x.guard_distance))")
Base.showcompact(io::IO, x::IndoorsNetwork) =
    print(io, "Indoors($(length(x.BSs)), $(x.num_MSs_per_cell), $(x.corridor_length), $(x.corridor_width), $(x.guard_distance))")

# The default parameter values are taken from ITU-R M.2135-1.
function setup_indoors_network(
    num_BSs, num_MSs_per_cell, num_MS_antennas, num_BS_antennas;
    system = SinglecarrierSystem(3.4e9, 15e3),
    propagation_environments = (@compat Dict(:LoS => SimpleLargescaleFadingEnvironment(16.9, 32.8 + 20*log10(3.4), 0, 3),
                                            :NLoS => SimpleLargescaleFadingEnvironment(43.3, 11.5 + 20*log10(3.4), 0, 4))),
    corridor_length = 120.,
    corridor_width = 50.,
    guard_distance = 3.,
    transmit_power = 10^(-9.8/10), transmit_powers = transmit_power*ones(Float64, num_BSs),
    BS_antenna_gain_params = [ OmnidirectionalAntennaParams(0) for i = 1:num_BSs ],
    user_priority = 1., user_priorities = ones(Float64, num_BSs*num_MSs_per_cell),
    num_streams = 1, num_streamss = num_streams*ones(Int, num_BSs*num_MSs_per_cell),
    MS_antenna_gain_dB = 0., MS_antenna_gains_dB = MS_antenna_gain_dB*ones(Float64, num_BSs*num_MSs_per_cell),
    receiver_noise_figure = 7., receiver_noise_figures = receiver_noise_figure*ones(Float64, num_BSs*num_MSs_per_cell))

    isa(num_MS_antennas, Vector) || (num_MS_antennas = num_MS_antennas*ones(Int, num_BSs*num_MSs_per_cell))
    isa(num_BS_antennas, Vector) || (num_BS_antennas = num_BS_antennas*ones(Int, num_BSs))

    # Place BSs
    y = corridor_width/2
    Δ = corridor_length/num_BSs
    BSs = Array(PhysicalBS, 0)
    for i = 1:num_BSs
        # BSs uniformly placed along corridor
        push!(BSs, PhysicalBS(num_BS_antennas[i], Position((i-1)*Δ + Δ/2, y), transmit_powers[i], BS_antenna_gain_params[i]))
    end

    MSs = [ PhysicalMS(num_MS_antennas[k], Position(0, 0), Velocity(0, 0), user_priorities[k], num_streamss[k], MS_antenna_gains_dB[k], receiver_noise_figures[k], SimpleLargescaleFadingEnvironmentState(zeros(Float64, num_BSs), falses(num_BSs))) for k = 1:num_BSs*num_MSs_per_cell ]

    IndoorsNetwork(MSs, BSs, system, num_MSs_per_cell,
        propagation_environments, corridor_length, corridor_width, guard_distance)
end

##########################################################################
# Standard cell assignment functions
function IDCellAssignment!(channel, network::IndoorsNetwork)
    Kc = network.num_MSs_per_cell; I = get_num_BSs(network)
    cell_assignment = Array(Int, I*Kc)

    for i = 1:I
        cell_assignment[(i-1)*Kc+1:i*Kc] = i
    end

    network.assignment = Assignment(cell_assignment, I)

    return AssignmentResults()
end

# Do this for now, even though the shadowing might actually give a certain MS
# better connection to a BS in another cell.
LargeScaleFadingCellAssignment!(channel, network::IndoorsNetwork) =
    IDCellAssignment!(channel, network)

##########################################################################
# Simulation functions
function draw_user_drop!{MS_t <: PhysicalMS, BS_t <: PhysicalBS, System_t <: System}(network::IndoorsNetwork{MS_t, BS_t, System_t, SimpleLargescaleFadingEnvironment})
    I = get_num_BSs(network); K = get_num_MSs(network)
    Δ = network.corridor_length/I

    # ITU-R M.2135-1, p. 33
    P_LoS = d -> begin
        if d <= 18
            return 1.
        elseif d < 37
            return exp(-(d - 18)/27)
        else
            return 0.5
        end
    end

    for k = 1:K
        i = div(k - 1, network.num_MSs_per_cell) + 1 # serving BS id

        # MS position
        while true
            x = (i-1)*Δ + Δ*rand()
            y = network.corridor_width*rand()
            candidate_position = Position(x, y)

            if get_distance(network.BSs[i].position, candidate_position) > network.guard_distance
                network.MSs[k].position = candidate_position
                break
            end
        end

        # Shadow fading and LoS/NLoS. Uncorrelated between BS links.
        shadow_realizations_dB = zeros(Float64, I)
        LoSs = BitArray(I)
        for j = 1:I
            if rand() <= P_LoS(get_distance(network.BSs[j].position, network.MSs[k].position))
                LoSs[j] = true
                shadow_realizations_dB[j] = network.propagation_environments[:LoS].shadow_sigma_dB*randn()
            else
                LoSs[j] = false
                shadow_realizations_dB[j] = network.propagation_environments[:NLoS].shadow_sigma_dB*randn()
            end
        end
        network.MSs[k].propagation_environment_state =
            SimpleLargescaleFadingEnvironmentState(shadow_realizations_dB, LoSs)
    end
end

function draw_channel{MS_t <: PhysicalMS, BS_t <: PhysicalBS}(network::IndoorsNetwork{MS_t, BS_t, SinglecarrierSystem, SimpleLargescaleFadingEnvironment})
    K = get_num_MSs(network); I = get_num_BSs(network)
    Ns = get_num_MS_antennas(network); Ms = get_num_BS_antennas(network)

    coefs = Array(Matrix{Complex128}, K, I)
    large_scale_fading_factor = Array(Float64, K, I)

    distances = get_distances(network)

    for k = 1:K
        for i = 1:I
            # Small scale fading
            coefs[k,i] = (1/sqrt(2))*(randn(Ns[k], Ms[i])
                                 + im*randn(Ns[k], Ms[i]))

            # Pathloss
            if network.MSs[k].propagation_environment_state.LoSs[i]
                pathloss_factor =
                    sqrt(10^(-(network.propagation_environments[:LoS].pathloss_beta + network.propagation_environments[:LoS].pathloss_alpha*log10(distances[k,i]))/10))
            else
                pathloss_factor =
                    sqrt(10^(-(network.propagation_environments[:NLoS].pathloss_beta + network.propagation_environments[:NLoS].pathloss_alpha*log10(distances[k,i]))/10))
            end

            # Shadow fading
            shadow_factor = sqrt(10^(network.MSs[k].propagation_environment_state.shadow_realizations_dB[i]/10))

            # BS antenna gain
            bs_antenna_gain = sqrt(get_antenna_gain(network.BSs[i].antenna_params))

            # MS antenna gain
            ms_antenna_gain = sqrt(10^(network.MSs[k].antenna_gain_dB/10))

            # Apply scale factors
            large_scale_fading_factor[k,i] = pathloss_factor*shadow_factor*bs_antenna_gain*ms_antenna_gain
            coefs[k,i] *= large_scale_fading_factor[k,i]
        end
    end

    return SinglecarrierChannel(coefs, Ns, Ms, K, I, large_scale_fading_factor)
end

##########################################################################
# Visualization functions
function plot_network_layout(network::IndoorsNetwork)
    I = get_num_BSs(network); K = get_num_MSs(network)
    Δ = network.corridor_length/I
    δ = network.corridor_width

    fig = PyPlot.figure()
    ax = fig[:add_subplot](1, 1, 1)

    # Corridor outline
    pos1 = network.BSs[1].position; pos2 = network.BSs[end].position
    sw = (pos1.x - Δ/2, pos1.y - δ/2); se = (pos2.x + Δ/2, pos2.y - δ/2)
    nw = (pos1.x - Δ/2, pos1.y + δ/2); ne = (pos2.x + Δ/2, pos2.y + δ/2)
    ax[:plot]([sw[1], se[1]], [sw[2], se[2]]; color="b", linestyle="-", linewidth=4)
    ax[:plot]([se[1], ne[1]], [se[2], ne[2]]; color="b", linestyle="-", linewidth=4)
    ax[:plot]([ne[1], nw[1]], [ne[2], nw[2]]; color="b", linestyle="-", linewidth=4)
    ax[:plot]([nw[1], sw[1]], [nw[2], sw[2]]; color="b", linestyle="-", linewidth=4)

    # BSs
    for i = 1:I
        pos = network.BSs[i].position
        ax[:plot](pos.x, pos.y; marker="x", color="b", markersize=10)

        # Cell outline
        if i < I
            ax[:plot]([pos.x + Δ/2, pos.x + Δ/2], [pos.y - δ/2, pos.y + δ/2]; color="b", linestyle="--", linewidth=2)
        end
    end

    # MSs
    for k = 1:K
        pos = network.MSs[k].position
        ax[:plot](pos.x, pos.y; marker="o", color="r", markersize=10)
    end

    display(fig)
end


text:
mutable struct BarplotGraphics{R<:Number} <: GraphicsArea
    bars::Vector{R}
    color::ColorType
    char_width::Int
    visible::Bool
    max_freq::Number
    max_len::Int
    symbols::AbstractVector{String}
    transform

    function BarplotGraphics(
        bars::AbstractVector{R},
        char_width::Int,
        visible::Bool,
        color::UserColorType,
        symbols::AbstractVector{S},
        transform,
    ) where {R,S<:Union{Char,String}}
        for s in symbols
            length(s) == 1 || throw(
                ArgumentError("Symbol has to be a single character, got: \"" * s * "\""),
            )
        end
        transform_func = transform isa Symbol ? getfield(Main, transform) : transform
        char_width = max(char_width, 10)
        max_freq, i = findmax(transform_func.(bars))
        max_len = length(string(bars[i]))
        char_width = max(char_width, max_len + 7)
        new{R}(
            bars,
            crayon_256_color(color),
            char_width,
            visible,
            max_freq,
            max_len,
            map(string, symbols),
            transform_func,
        )
    end
end

nrows(g::BarplotGraphics) = length(g.bars)
ncols(g::BarplotGraphics) = g.char_width

BarplotGraphics(
    bars::AbstractVector{R},
    char_width::Int,
    transform = :identity;
    visible::Bool = true,
    color::UserColorType = :green,
    symbols = ['■'],
) where {R<:Number} =
    BarplotGraphics(bars, char_width, visible, crayon_256_color(color), symbols, transform)

function addrow!(g::BarplotGraphics{R}, bars::AbstractVector{R}) where {R<:Number}
    append!(g.bars, bars)
    g.max_freq, i = findmax(g.transform.(g.bars))
    g.max_len = length(string(g.bars[i]))
    g
end

function addrow!(g::BarplotGraphics{R}, bar::Number) where {R<:Number}
    push!(g.bars, R(bar))
    g.max_freq, i = findmax(g.transform.(g.bars))
    g.max_len = length(string(g.bars[i]))
    g
end

function printrow(io::IO, g::BarplotGraphics, row::Int)
    0 < row <= nrows(g) || throw(ArgumentError("Argument \"row\" out of bounds: $row"))
    bar = g.bars[row]
    max_freq = g.max_freq
    max_bar_width = max(g.char_width - 2 - g.max_len, 1)
    val = g.transform(bar)
    nsyms = length(g.symbols)
    frac = float(max_freq > 0 ? max(val, zero(val)) / max_freq : 0)
    bar_head = round(Int, frac * max_bar_width, nsyms > 1 ? RoundDown : RoundNearestTiesUp)
    print_color(g.color, io, max_freq > 0 ? repeat(g.symbols[nsyms], bar_head) : "")
    if nsyms > 1
        rem = (frac * max_bar_width - bar_head) * (nsyms - 2)
        print_color(g.color, io, rem > 0 ? g.symbols[1 + round(Int, rem)] : " ")
        bar_head += 1  # padding, we printed one more char
    end
    bar_lbl = string(bar)
    if bar >= 0
        print_color(:normal, io, " ", bar_lbl)
        len = length(bar_lbl)
    else
        len = -1
    end
    pad_len = max(max_bar_width + 1 + g.max_len - bar_head - len, 0)
    print(io, repeat(" ", round(Int, pad_len)))
    nothing
end


text:
# types used by the module

#################### Node #######################
abstract Node

#################### Factor #######################
"Defines a factor, a term in the variational objective."
abstract Factor{N}

###################################################
# RandomNode <: Node
###################################################
immutable RandomNode{D <: Distribution} <: Node
    name::Symbol
    innerinds::Vector{Symbol}
    outerinds::Vector{Symbol}
    data::Array{D}

    function RandomNode(name, indices, data)
        # handle the case where node is scalar in outer indices
        if length(data) == 1
            inds = vcat(indices, :scalar)
        else
            inds = copy(indices)
        end

        # inner indices are assumed listed first
        ninds = length(inds)
        nouter = ndims(data)
        @assert nouter <= ninds "Indices do not match data shape"
        outerinds = inds[ninds - nouter + 1:end]
        innerinds = inds[1:ninds - nouter]
        ninner = length(innerinds)
        @assert nouter == length(outerinds) "Indices do not match data shape."
        if ninner > 0
            # length(size(data[1])) gives number of dimensions in
            # the random distribution D
            @assert ninner == length(size(data[1])) "Inner indices, if provided, must match distribution shape."
        end
        new(name, innerinds, outerinds, data)
    end
end

RandomNode{D <: Distribution}(name::Symbol, indices::Vector{Symbol}, ::Type{D}, pars...) = begin
    dims = map(size, pars)
    #=
    if only one parameter passed, check if there is a matching constructor;
        if so, feed directly to D
        if not, assume this is a collection of parameters to be mapped over
    else if all pars have the same size, assume these are arrays of
        params to be mapped over;
    otherwise, assume they are tuples to be fed directly to D
    =#
    if length(pars) == 1
        if method_exists(call, (Type{D}, map(typeof, pars)...))
            RandomNode{D}(name, indices, [D(pars...)])
        else
            RandomNode{D}(name, indices, map(D, pars...))
        end
    elseif length(unique(dims)) == 1 && dims[1] != ()
        RandomNode{D}(name, indices, map(D, pars...))
    else
        RandomNode{D}(name, indices, [D(pars...)])
    end
end

###################################################
# ConstantNode <: Node
###################################################
immutable ConstantNode{T} <: Node
    name::Symbol
    innerinds::Vector{Symbol}
    outerinds::Vector{Symbol}
    data::Array{T}

    function ConstantNode(name, indices, data)
        if !isempty(indices)
            @assert ndims(data) == length(indices) "Indices do not match data shape."
        end
        outerinds = indices
        innerinds = Symbol[]
        new(name, innerinds, outerinds, data)
    end
end
ConstantNode{T}(name::Symbol, indices::Vector{Symbol}, data::Array{T}) = ConstantNode{T}(name, indices, data)
ConstantNode(name::Symbol, indices::Vector{Symbol}, data::Number) = ConstantNode(name, [:scalar], [data])
ConstantNode{T}(data::Array{T}, indices::Vector{Symbol}) =
    ConstantNode(gensym("const"), indices, data)
ConstantNode(x::Number) = ConstantNode(gensym("const"), [:scalar], [x])

function show(io::IO, n::Node)
    T = typeof(n)
    print(io, T.name, "{", T.parameters[1].name.name, "}", " ", n.name)
    inners = filter(x -> x != :scalar, n.innerinds)
    outers = filter(x -> x != :scalar, n.outerinds)
    if !isempty(inners)
        print(io, "(")
        _printinds(io, n.innerinds)
        print(io, ")")
    end
    if !isempty(outers)
        print(io, "[")
        _printinds(io, n.outerinds)
        print(io, "]")
    end
end

function _printinds(io::IO, inds::Vector{Symbol})
    for idx in inds
        lastind = last(inds)
        if idx == lastind
            print(io, idx)
        else
            print(io, idx, ",")
        end
    end
end

###################################################
# FactorInds
###################################################
immutable FactorInds
    indices::Vector{Symbol}  # all fully outer indices
    maxvals::Vector{Int}  # maximum value of each index
    #=
    The next two dicts take care of a pair of necessary mappings:
    - inds_in_factor maps a node symbol to the (integer) subset of factor
        indices involving this node
    - inds_in_node maps a node symbol to the (integer) subset of node indices
        involved in this factor
    =#
    inds_in_factor::Dict{Symbol, Vector{Int}}
    inds_in_node::Dict{Symbol, Vector{Int}}
end

function show(io::IO, f::Factor)
    println(io, typeof(f))
    for v in fieldnames(f)
        if v ∉ [:inds, :namemap]
            print(io, v, ": ")
            println(io, getfield(f, v))
        end
    end
end

###################################################
# ExprNode <: Node
###################################################
immutable ExprNode <: Node
    name::Symbol
    ex::Expr  # expression defining node
    nodedict::Dict{Symbol, Node}  # nodes in the expression
    inds::FactorInds  # so we can treat this like a factor
    innerinds::Vector{Symbol}
    outerinds::Vector{Symbol}
    dims::Vector{Int}
end

function ExprNode(name::Symbol, ex::Expr, nodelist::Vector{Node})
    explicit_inners = collect(get_all_inds(ex))
    stripped_ex = strip_inds(ex)
    fi = get_structure(explicit_inners, nodelist...)
    nodedict = Dict([n.name => n for n in nodelist])
    inners = union([n.innerinds for n in nodelist]...)
    outers = union([n.outerinds for n in nodelist]...)
    allinds = union(inners, outers)
    outerinds = fi.indices  # fully outer inds
    innerinds = setdiff(allinds, outerinds)
    dims = fi.maxvals  # size of each fully outer index

    # now get rid of :scalar index (if it's not the only one)
    scalar_inds = findin(outerinds, [:scalar])
    if !isempty(scalar_inds) && length(outerinds) > 1
        scalar_index = scalar_inds[1]  # should only have 1 or 0 entries
        splice!(outerinds, scalar_index)
        splice!(dims, scalar_index)
    end

    ExprNode(name, stripped_ex, nodedict, fi, innerinds, outerinds, dims)
end

###################################################
# ExprDist <: Distribution
###################################################
immutable ExprDist{V <: Val} <: Distribution
    nodedict::Dict{Symbol, Distribution}
end

#################### Model #######################
"Defines a Variational Bayes model."
type VBModel
    # nodes maps symbols to the nodes/groups of nodes associated with them
    nodes::Vector{Node}

    # all factors in the graph
    factors::Vector{Factor}

    # dictionary linking all random variables to a list of tuples;
    # each tuple gives a factor and the name of the random variable
    # in that factor
    graph::Dict{Node, Vector{Tuple{Factor, Symbol}}}

    # dictionary containing the update strategy for each node
    update_strategy::Dict{Node, Symbol}

    VBModel(nodes, factors) = begin
        m = new(nodes, factors)

        # build an empty graph
        m.graph = Dict{Node, Vector{Tuple{Factor, Symbol}}}()
        for n in nodes
            push!(m.graph, n => [])
        end

        # register factors
        for f in factors
            register(f, m)
        end

        # calculate default update strategy
        m.update_strategy = Dict{Node, Symbol}()
        for n in nodes
            if isa(n, ConstantNode)
                m.update_strategy[n] = :constant
            elseif isa(n, RandomNode) && check_conjugate(n, m)
                m.update_strategy[n] = :conjugate
            else
                m.update_strategy[n] = :undefined
            end
        end

        m

    end
end

# register a factor with its associated nodes in the graph
function register(f::Factor, m::VBModel)
    for var in fieldnames(f)
        n = getfield(f, var)
        if isa(n, Node)
            push!(m.graph[n], (f, var))
        end
    end
end

function update!(m::VBModel)
    for n in m.nodes
        update!(n, m, Val{m.update_strategy[n]})
    end
end

function update!(n::Node, m::VBModel, ::Type{Val{:constant}})
end


text:
#=
Julia function to create a .csv file with a concordance of matches of a regular expression in a directory of .txt files (other file types in the directory are ignored). The surrounding context includes up to the number of characters specified in the 'context_size' argument or all available context within the paragraph if the specified number exceeds the number of characters available. The functions returns an array of strings, each string with tabs that separate the filename of the file in which a match was found, the paragraph number, the preceding context, the match, and the following context.

(c) 2018 Earl K. Brown, ekbrown@byu.edu
=#

function concord(input_dir, search_regex, context_size = 50)

    file_names = filter(r"\.txt$"i, readdir(input_dir))

    # initializes collector array and pushes headers to it
    all_output = ["FILE\tPARA\tPRE\tMATCH\tPOST"]

    # loops over files
    for i in file_names

        # reads in current file
        cur_file = string(input_dir, "/", i)
        f = open(cur_file)
        cur_text = readstring(f)
        close(f)

        # deletes tabs in the file and splits up the file into paragraphs
        cur_text = replace(cur_text, r"\t", s" ")
        cur_text = split(cur_text, "\n")

        # loops over paragraphs in current file
        para_num = 0  # paragraph counter
        for j in cur_text
            para_num += 1  # increment paragraph counter

            # searches for regex in current paragraph and returns an iterator object
            hits = eachmatch(search_regex, j)

            # loops over hits in current paragraph
            for k in hits
                offset_start = k.offset
                offset_end = k.offset + length(k.match)

                # preceding context
                pre = try
                    j[offset_start - context_size:offset_start - 1]
                catch BoundsError
                    j[1:offset_start - 1]
                end

                # following context
                post = try
                    j[offset_end:offset_end + context_size]
                catch BoundsError
                    j[offset_end:end]
                end

                # creates output string and pushes it to the collector array
                cur_output = string(i, "\t", para_num, "\t", pre, "\t", k.match, "\t", post)
                push!(all_output, cur_output)

            end  # next hit in current paragraph; end k loop
        end  # next paragraph; end j loop
    end  # next file; end i loop

    return(all_output)

end  # function definition

### tests function
input_dir = "/pathway/to/directory"
search_regex = r"[REGEX HERE]"i  # case-insensitive search

results = concord(input_dir, search_regex)

writedlm("/pathway/to/output.csv", results, '\n', quotes=false, comments=false)
println("All done!")


text:
using SymbolicCRN
using Documenter

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

makedocs(;
    modules=[SymbolicCRN],
    authors="Laura Brustenga i Moncusí <brust@math.ku.dk> and contributors",
    repo="https://github.com/LauraBMo/SymbolicCRN.jl/blob/{commit}{path}#{line}",
    sitename="SymbolicCRN.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://LauraBMo.github.io/SymbolicCRN.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/LauraBMo/SymbolicCRN.jl",
)


text:
using Test, NodeCall
using Pkg

if dirname(Pkg.project().path) != @__DIR__
    Pkg.activate(@__DIR__)
end

if !NodeCall.initialized()
    NodeCall.initialize()
end


text:
# 1) given a Vector{T}, T<:Number, construct an object of type S using the same underlying data
# e.g.
#  array = [1.0,2.0,3.0], S = Tuple{Complex{Float64},Float64} => (1.0 + 2.0im, 3.0)

# 2) for an Array{T,N}, and some 1 <= D <= N,
#   reconstruct S from a slice of array#

# The main difference between StructArrays and what we want is that we want 1 underlying (rather than # fields)
# there is no heterogeneity in the struct (can store recursive representations) of field types

# get_offset(array, S, offset) => S(...), next_offset

"""
    is_valid_basetype(::Type{T}, ::Type{S})

Determines whether an object of type `S` can be stored in an array with elements
of type `T` by recursively checking whether the non-empty fields of `S` can be
stored in such an array. If `S` is empty, this is always true.
"""
is_valid_basetype(::Type{T}, ::Type{S}) where {T, S} =
    sizeof(S) == 0 ||
    (fieldcount(S) > 0 && is_valid_basetype(T, fieldtypes(S)...))
is_valid_basetype(::Type{T}, ::Type{<:T}) where {T} = true
is_valid_basetype(::Type{T}, ::Type{S}, Ss...) where {T, S} =
    is_valid_basetype(T, S) && is_valid_basetype(T, Ss...)

"""
    check_basetype(::Type{T}, ::Type{S})

Check whether the types `T` and `S` have well-defined sizes, and whether an
object of type `S` can be stored in an array with elements of type `T`.
"""
function check_basetype(::Type{T}, ::Type{S}) where {T, S}
    if @generated
        if !isbitstype(T)
            estr = "Base type $T has indeterminate size"
            :(error($estr))
        elseif !isbitstype(S)
            estr = "Struct type $S has indeterminate size"
            :(error($estr))
        elseif !is_valid_basetype(T, S)
            estr = "Struct type $S cannot be represented using base type $T"
            :(error($estr))
        else
            :(return nothing)
        end
    else
        isbitstype(T) || error("Base type $T has indeterminate size")
        isbitstype(S) || error("Struct type $S has indeterminate size")
        is_valid_basetype(T, S) ||
            error("Struct type $S cannot be represented using base type $T")
        return nothing
    end
end

"""
    replace_basetype(::Type{T}, ::Type{T′}, ::Type{S})

Changes the type parameters of `S` to produce a new type `S′` such that, if
`is_valid_basetype(T, S)` is true, then so is `is_valid_basetype(T′, S′)`.
"""
replace_basetype(::Type{T}, ::Type{T′}, ::Type{S}) where {T, T′, S} =
    length(S.parameters) == 0 ? S :
    S.name.wrapper{replace_basetypes(T, T′, S.parameters...)...}
replace_basetype(::Type{T}, ::Type{T′}, ::Type{<:T}) where {T, T′} = T′
replace_basetype(::Type{T}, ::Type{T′}, value) where {T, T′} = value
replace_basetypes(::Type{T}, ::Type{T′}, value) where {T, T′} =
    (replace_basetype(T, T′, value),)
replace_basetypes(::Type{T}, ::Type{T′}, value, values...) where {T, T′} =
    (replace_basetype(T, T′, value), replace_basetypes(T, T′, values...)...)
# TODO: This could potentially lead to some annoying bugs, since it replaces
# type parameters instead of field types. So, if `S` has `Float64` as a
# parameter, `replace_basetype(Float64, Float32, S)` will replace that parameter
# with `Float32`, regardless of whether the parameter corresponds to any values
# stored in an object of type `S`. In other words, if a user utilizes types as
# type parameters and specializes their code on those parameters, this may
# change the results of their code.
# Note that there is no way to write `replace_basetype` using field types
# instead of type parameters, since replacing the field types of an object will
# change that object's type parameters, but there is no general way to map field
# types to type parameters.

"""
    parent_array_type(::Type{<:AbstractArray})

Returns the parent array type underlying any wrapper types, with all
dimensionality information removed.
"""
parent_array_type(::Type{<:Array{T}}) where {T} = Array{T}
parent_array_type(::Type{<:SubArray{T, N, A}}) where {T, N, A} =
    parent_array_type(A)

"""
    promote_parent_array_type(::Type{<:AbstractArray}, ::Type{<:AbstractArray})

Given two parent array types (without any dimensionality information), promote
both the element types and the array types themselves.
"""
promote_parent_array_type(::Type{Array{T1}}, ::Type{Array{T2}}) where {T1, T2} =
    Array{promote_type(T1, T2)}

"""
    StructArrays.bypass_constructor(T, args)

Create an instance of type `T` from a tuple of field values `args`, bypassing
possible internal constructors. `T` should be a concrete type.
"""
@generated function bypass_constructor(::Type{T}, args) where {T}
    vars = ntuple(_ -> gensym(), fieldcount(T))
    assign = [
        :($var::$(fieldtype(T, i)) = getfield(args, $i)) for
        (i, var) in enumerate(vars)
    ]
    construct = Expr(:new, :T, vars...)
    Expr(:block, assign..., construct)
end

"""
    fieldtypeoffset(T,S,i)

Similar to `fieldoffset(S,i)`, but gives result in multiples of `sizeof(T)` instead of bytes.
"""
fieldtypeoffset(::Type{T}, ::Type{S}, i) where {T, S} =
    Int(div(fieldoffset(S, i), sizeof(T)))

"""
    typesize(T,S)

Similar to `sizeof(S)`, but gives the result in multiples of `sizeof(T)`.
"""
typesize(::Type{T}, ::Type{S}) where {T, S} = div(sizeof(S), sizeof(T))

# TODO: this assumes that the field struct zero type is the same as the backing
# zero'd out memory, which should be true in all "real world" cases
# but is something that should be revisited
@inline function _mzero!(out::MArray{S, T, N, L}, FT) where {S, T, N, L}
    TdivFT = DataLayouts.typesize(FT, T)
    Base.GC.@preserve out begin
        @inbounds for i in 1:(L * TdivFT)
            Base.unsafe_store!(
                Base.unsafe_convert(Ptr{FT}, Base.pointer_from_objref(out)),
                zero(FT),
                i,
            )
        end
    end
    return out
end

"""
    get_struct(array, S[, offset=0])

Construct an object of type `S` from the values of `array`, optionally offset by `offset` from the start of the array.
"""
function get_struct(array::AbstractArray{T}, ::Type{S}, offset) where {T, S}
    if @generated
        tup = :(())
        for i in 1:fieldcount(S)
            push!(
                tup.args,
                :(get_struct(
                    array,
                    fieldtype(S, $i),
                    offset + fieldtypeoffset(T, S, $i),
                )),
            )
        end
        return quote
            Base.@_propagate_inbounds_meta
            bypass_constructor(S, $tup)
        end
    else
        Base.@_propagate_inbounds_meta
        args = ntuple(fieldcount(S)) do i
            get_struct(array, fieldtype(S, i), offset + fieldtypeoffset(T, S, i))
        end
        return bypass_constructor(S, args)
    end
end

# recursion base case: hit array type is the same as the struct leaf type
@propagate_inbounds function get_struct(
    array::AbstractArray{S},
    ::Type{S},
    offset,
) where {S}
    return array[offset + 1]
end

@inline function get_struct(array::AbstractArray{T}, ::Type{S}) where {T, S}
    @inbounds get_struct(array, S, 0)
end

function set_struct!(array::AbstractArray{T}, val::S, offset) where {T, S}
    if @generated
        ex = quote
            Base.@_propagate_inbounds_meta
        end
        for i in 1:fieldcount(S)
            push!(
                ex.args,
                :(set_struct!(
                    array,
                    getfield(val, $i),
                    offset + fieldtypeoffset(T, S, $i),
                )),
            )
        end
        push!(ex.args, :(return val))
        return ex
    else
        Base.@_propagate_inbounds_meta
        for i in 1:fieldcount(S)
            set_struct!(
                array,
                getfield(val, i),
                offset + fieldtypeoffset(T, S, i),
            )
        end
        return val
    end
end

@propagate_inbounds function set_struct!(
    array::AbstractArray{S},
    val::S,
    offset,
) where {S}
    array[offset + 1] = val
    val
end

@inline function set_struct!(array, val)
    @inbounds set_struct!(array, val, 0)
end

if VERSION >= v"1.7.0-beta1"
    # For complex nested types (ex. wrapped SMatrix) we hit a recursion limit and de-optimize
    # We know the recursion will terminate due to the fact that bitstype fields
    # cannot be self referential so there are no cycles in get/set_struct (bounded tree)
    # TODO: enforce inference termination some other way
    if hasfield(Method, :recursion_relation)
        dont_limit = (args...) -> true
        for m in methods(get_struct)
            m.recursion_relation = dont_limit
        end
        for m in methods(set_struct!)
            m.recursion_relation = dont_limit
        end
    end
end


text:
using ValidatedNumerics

import Base.mod2pi

function standard_map(X::IntervalBox, k = 1.0)
    p, θ = X

    p′ = p + k*sin(θ)
    @show p′
    θ′ = mod2pi( θ + p′ )
    p′ = mod2pi(p′)

    @show p′, θ′
    @show typeof(p′), typeof(θ′)

    IntervalBox(p′, θ′)
end

function IntervalBox{T}(X::Vector{Interval{T}}, Y::Vector{Interval{T}})
    vec([IntervalBox(x, y) for x in X, y in Y])
end

function iterate(f, n, x)
    for i in 1:n
        x = f(x)
    end
    x
end

import Base.mod
function mod(X::Interval, width::Real)

    @show X, width

    X /= width

    if diam(X) >= 1.

        return [Interval(0, 1) * width]
    end

    a = X.lo - floor(X.lo)
    b = X.hi - floor(X.hi)

    if a < b
        return [Interval(a, b)*width]

    end

    return [Interval(0, b)*width, Interval(a, 1)*width]

end

function mod(X::IntervalBox, width::Real)
    x, y = X

    xx = mod(x, width)
    yy = mod(y, width)

    vec([IntervalBox(x, y) for x in xx, y in yy])
end

mod2pi{T}(x::Interval{T}) = mod(x, ValidatedNumerics.two_pi(T))

mod2pi{T}(X::Vector{Interval{T}}) = vcat(map(mod2pi, X)...)


text:
module DepParsing

importall ..JukaiNLP
using Merlin
using Compat
using ProgressMeter
using TransitionParser: beamsearch, max_violation!, to_seq

include("my_adagrad.jl")
include("../iddict.jl")
include("depparser.jl")
include("token.jl")
include("arcstd.jl")
include("util.jl")
include("feedforward.jl")
include("structuredfeedforward.jl")
include("unlabeled_feedforward.jl")
include("struct_extras.jl")
include("accuracy.jl")
include("io.jl")
include("saver.jl")
include("extra/utils.jl")
include("extra/perceptron.jl")

export DepParser,
       Unlabeled,
       Labeled,
       FeedForward,
       UnlabeledFeedForward,
       LabelerFF,
       assign!,
       StructuredFeedForward,
       Perceptron,
       readconll,
       train!,
       decode,
       evaluate,
       toconll,
       initmodel!,
       MyAdaGrad

end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
#!/usr/bin/env julia
using Printf

function ζ(N::Int64,a::Int64)
  s = 0.0
  @simd for i ∈ 1:N
    s+=(1.0/i)^a
  end
  s
end

N = 16000000000
@time ζ_2 = ζ(N,2)
π_squared_over_6 = π^2 / 6.0;

@printf "      ζ(2) = %18.16f\n" ζ_2
@printf "π^2 / 6.0  = %18.16f\n" π_squared_over_6
@printf "     error = %18.16f\n" π_squared_over_6 - ζ_2

max_stars_repo_name: joelandman/rzf

max_issues_repo_name: joelandman/rzf

max_forks_repo_name: joelandman/rzf


text:
using MPI, Test
include("../testhelpers.jl")

@testset "DGmethods" begin
  tests = [
    (1,"integral_test.jl")
    (1,"integral_test_sphere.jl")
    (1, "Euler/isentropicvortex.jl")
    (1, "Euler/isentropicvortex-imex.jl")
    (1, "Euler/isentropicvortex-multirate.jl")
    (1, "advection_diffusion/pseudo1D_advection_diffusion.jl")
    (1, "compressible_Navier_Stokes/ref_state.jl")
    (1, "horizontal_integral_test.jl")
    (1, "courant.jl")
   ]

  runmpi(tests, @__FILE__)
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
### Julia OpenStreetMap Package ###
### MIT License                 ###
### Copyright 2014              ###

### Functions for plotting using the Winston package ###

typealias Styles @compat(Union{Style, Dict{Int,Style}})

### Generic Map Plot ###
function plotMap(nodes::@compat(Union{Dict{Int,LLA},Dict{Int,ENU}}) ;
                 highways::@compat(Union{@compat(Void),Dict{Int,Highway}}) = nothing,
                 buildings::@compat(Union{@compat(Void),Dict{Int,Building}}) = nothing,
                 features::@compat(Union{@compat(Void),Dict{Int,Feature}}) = nothing,
                 bounds::@compat(Union{@compat(Void),Bounds}) = nothing,
                 intersections::@compat(Union{@compat(Void),Dict{Int,Intersection}}) = nothing,
                 roadways=nothing,
                 cycleways=nothing,
                 walkways=nothing,
                 feature_classes::@compat(Union{@compat(Void),Dict{Int,Int}}) = nothing,
                 building_classes::@compat(Union{@compat(Void),Dict{Int,Int}}) = nothing,
                 route::@compat(Union{@compat(Void),Vector{Int},Vector{Vector{Int}}}) = nothing,
                 highway_style::Styles=Style(0x007CFF, 1.5, "-"),
                 building_style::Styles=Style(0x000000, 1, "-"),
                 feature_style::Styles=Style(0xCC0000, 2.5, "."),
                 route_style::@compat(Union{Style,Vector{Style}}) = Style(0xFF0000, 3, "-"),
                 intersection_style::Style=Style(0x000000, 3, "."),
                 width::Integer=500,
                 fontsize::Integer=0,
                 km::Bool=false,
                 realtime::Bool=false)

    # Chose labels according to point type and scale
    xlab, ylab = if isa(nodes, Dict{Int,LLA})
        "Longitude (deg)", "Latitude (deg)"
    elseif km
        "East (km)", "North (km)"
    else
        "East (m)", "North (m)"
    end

    # Waiting for Winston to add capability to force equal scales. For now:
    if VERSION.minor < 4
        height = isa(bounds, @compat(Void)) ? width : int(width / aspectRatio(bounds))
    else
        height = isa(bounds, @compat(Void)) ? width : round(Int, width / aspectRatio(bounds))
    end

    # Create the figure
    fignum = Winston.figure(name="OpenStreetMap Plot", width=width, height=height)

    if isa(bounds, @compat(Void))
        p = Winston.FramedPlot("xlabel", xlab, "ylabel", ylab)
    else # Limit plot to specified bounds
        Winston.xlim(bounds.min_x, bounds.max_x)
        Winston.ylim(bounds.min_y, bounds.max_y)

        if km && isa(nodes, Dict{Int,ENU})
            xrange = (bounds.min_x/1000, bounds.max_x/1000)
            yrange = (bounds.min_y/1000, bounds.max_y/1000)
        else
            xrange = (bounds.min_x, bounds.max_x)
            yrange = (bounds.min_y, bounds.max_y)
        end

        p = Winston.FramedPlot("xlabel", xlab, "ylabel", ylab, xrange=xrange, yrange=yrange)
    end

    # Iterate over all buildings and draw
    if !isa(buildings, @compat(Void))
        if !isa(building_classes, @compat(Void))
            if isa(building_style, Dict{Int,Style})
                drawWayLayer(p, nodes, buildings, building_classes, building_style, km, realtime)
            else
                drawWayLayer(p, nodes, buildings, building_classes, LAYER_BUILDINGS, km, realtime)
            end
        else
            for (key, building) in buildings
                # Get coordinates of all nodes for object
                coords = getNodeCoords(nodes, building.nodes, km)

                # Add line(s) to plot
                drawNodes(p, coords, building_style, realtime)
            end
        end
    end

    # Iterate over all highways and draw
    if !isa(highways, @compat(Void))
        if !(nothing == roadways == cycleways == walkways)
            if !isa(roadways, @compat(Void))
                if isa(highway_style, Dict{Int,Style})
                    drawWayLayer(p, nodes, highways, roadways, highway_style, km, realtime)
                else
                    drawWayLayer(p, nodes, highways, roadways, LAYER_STANDARD, km, realtime)
                end
            end
            if !isa(cycleways, @compat(Void))
                if isa(highway_style, Dict{Int,Style})
                    drawWayLayer(p, nodes, highways, cycleways, highway_style, km, realtime)
                else
                    drawWayLayer(p, nodes, highways, cycleways, LAYER_CYCLE, km, realtime)
                end
            end
            if !isa(walkways, @compat(Void))
                if isa(highway_style, Dict{Int,Style})
                    drawWayLayer(p, nodes, highways, walkways, highway_style, km, realtime)
                else
                    drawWayLayer(p, nodes, highways, walkways, LAYER_PED, km, realtime)
                end
            end
        else
            for (key, highway) in highways
                # Get coordinates of all nodes for object
                coords = getNodeCoords(nodes, highway.nodes, km)

                # Add line(s) to plot
                drawNodes(p, coords, highway_style, realtime)
            end
        end
    end

    # Iterate over all features and draw
    if !isa(features, @compat(Void))
        if !isa(feature_classes, @compat(Void))
            if isa(feature_style, Dict{Int,Style})
                drawFeatureLayer(p, nodes, features, feature_classes, feature_style, km, realtime)
            else
                drawFeatureLayer(p, nodes, features, feature_classes, LAYER_FEATURES, km, realtime)
            end
        else
            coords = getNodeCoords(nodes, collect(keys(features)), km)

            # Add feature point(s) to plot
            drawNodes(p, coords, feature_style, realtime)
        end
    end

    # Draw route
    if isa(route, Vector{Int})
        # Get coordinates of all nodes for route
        coords = getNodeCoords(nodes, route, km)

        # Add line(s) to plot
        drawNodes(p, coords, route_style, realtime)
    elseif isa(route, Vector{Vector{Int}})
        for k = 1:length(route)
            coords = getNodeCoords(nodes, route[k], km)
            if isa(route_style, Vector{Style})
                drawNodes(p, coords, route_style[k], realtime)
            else
                drawNodes(p, coords, route_style, realtime)
            end
        end
    end

    # Iterate over all intersections and draw
    if !isa(intersections, @compat(Void))
        coords = Array(Float64, length(intersections), 2)
        k = 1
        for key in keys(intersections)
            coords[k, :] = getNodeCoords(nodes, key, km)
            k += 1
        end

        # Add intersection(s) to plot
        drawNodes(p, coords, intersection_style, realtime)
    end

    if fontsize > 0
        attr = Compat.@Dict(:fontsize => fontsize)
        Winston.setattr(p.x1, "label_style", attr)
        Winston.setattr(p.y1, "label_style", attr)
        Winston.setattr(p.x1, "ticklabels_style", attr)
        Winston.setattr(p.y1, "ticklabels_style", attr)
    end

    display(p)

    # Return figure object (enables further manipulation)
    return p
end

### Draw layered Map ###
function drawWayLayer(p::Winston.FramedPlot, nodes::Dict, ways, classes, layer, km=false, realtime=false)
    for (key, class) in classes
        # Get coordinates of all nodes for object
        if haskey(ways,key)
            coords = getNodeCoords(nodes, ways[key].nodes, km)

            # Add line(s) to plot
            drawNodes(p, coords, layer[class], realtime)
        end
    end
end

### Draw layered features ###
function drawFeatureLayer(p::Winston.FramedPlot, nodes::Dict, features, classes, layer, km=false, realtime=false)

    for id in unique(values(classes))
        ids = Int[]

        for (key, class) in classes
            if class == id
                push!(ids, key)
            end
        end

        # Get coordinates of node for object
        coords = getNodeCoords(nodes, ids, km)

        # Add point to plot
        drawNodes(p, coords, layer[id], realtime)
    end
end

### Get coordinates of lists of nodes ###
# Nodes in LLA coordinates
function getNodeCoords(nodes::Dict{Int,LLA}, id_list, km=false)
    coords = Array(Float64, length(id_list), 2)

    for k = 1:length(id_list)
        loc = nodes[id_list[k]]
        coords[k, 1] = loc.lon
        coords[k, 2] = loc.lat
    end

    return coords
end

# Nodes in ENU coordinates
function getNodeCoords(nodes::Dict{Int,ENU}, id_list, km=false)
    coords = Array(Float64, length(id_list), 2)

    for k = 1:length(id_list)
        loc = nodes[id_list[k]]
        coords[k, 1] = loc.east
        coords[k, 2] = loc.north
    end

    if km
        coords /= 1000
    end

    return coords
end

### Draw a line between all points in a coordinate list ###
function drawNodes(p::Winston.FramedPlot, coords, style="k-", width=1, realtime=false)
    x = coords[:, 1]
    y = coords[:, 2]
    if length(x) > 1
        if realtime
            display(Winston.plot(p, x, y, style, linewidth=width))
        else
            Winston.plot(p, x, y, style, linewidth=width)
        end
    end
    nothing
end

### Draw a line between all points in a coordinate list given Style object ###
function drawNodes(p::Winston.FramedPlot, coords, line_style::Style, realtime=false)
    x = coords[:, 1]
    y = coords[:, 2]
    if length(x) > 1
        if realtime
            display(Winston.plot(p, x, y, line_style.spec, color=line_style.color, linewidth=line_style.width))
        else
            Winston.plot(p, x, y, line_style.spec, color=line_style.color, linewidth=line_style.width)
        end
    end
    nothing
end

### Compute approximate "aspect ratio" at mean latitude ###
function aspectRatio(bounds::Bounds{LLA})
    c_adj = cosd((bounds.min_y + bounds.max_y) / 2)
    range_y = bounds.max_y - bounds.min_y
    range_x = bounds.max_x - bounds.min_x

    return range_x * c_adj / range_y
end

### Compute exact "aspect ratio" ###
function aspectRatio(bounds::Bounds{ENU})
    range_y = bounds.max_y - bounds.min_y
    range_x = bounds.max_x - bounds.min_x

    return range_x / range_y
end


text:
module AtomicDFTTests
using AtomicDFT
using AxisArrays
using Unitful, UnitfulHartree
using Base.Test
using LibXC

@testset "Orbital Indexing" begin
  include("OrbitalIndexing.jl")
end

@testset "Object creation" begin
    include("ObjectCreation.jl")
end

@testset "Radial" begin
    include("Radial.jl")
end

# @testset "Potentials" begin
#     include("Potentials.jl")
# end

end


text:
#-----------------------------------------------------------------------
# JuMPeR  --  JuMP Extension for Robust Optimization
# http://github.com/IainNZ/JuMPeR.jl
#-----------------------------------------------------------------------
# Copyright (c) 2016: Iain Dunning
# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at http://mozilla.org/MPL/2.0/.
#-----------------------------------------------------------------------
# example/inventory.jl
# Implements the inventory management example from
#   A. Ben-Tal, A. Goryashko, E. Guslitzer, A. Nemirovski
#   "Adjustable Robust Solutions of Uncertain Linear Programs"
# Requires a linear optimization solver.
#-----------------------------------------------------------------------

using JuMP, JuMPeR, Gurobi

# Define problem parameters
I = 3           # Number of factories
T = 24          # Number of time periods
d_nom = 1000*[1 + 0.5*sin(π*(t-1)/12)  for t in 1:T]  # Nominal demand
θ = 0.20        # Uncertainty level
α = [1.0, 1.5, 2.0]  # Production costs
c = [α[i] * (1 + 0.5*sin(π*(t-1)/12))  for i in 1:I, t in 1:T]
P = 567         # Maximimum production per period
Q = 13600       # Maximumum production over all
Vmin = 500      # Minimum inventory at warehouse
Vmax = 2000     # Maximum inventory at warehouse
v1 = Vmin       # Initial inventory (not provided in paper)

# Setup robust model
invmgmt = RobustModel()

# Uncertain parameter: demand at each time stage lies in a interval
@uncertain(invmgmt, d_nom[t]*(1-θ) <= d[t=1:T] <= d_nom[t]*(1+θ))

# Decision: how much to produce at each factory at each time
# As this decision can be updated as demand is realized, we will use adaptive
# policy - in particular, an affine policy where production at time t is an
# affine function of the demand realized previously.
@adaptive(invmgmt, p[i=1:I,t=1:T], policy=Affine, depends_on=d[1:t-1])

# Objective: minimize total cost of production
@variable(invmgmt, F)  # Overall cost
@objective(invmgmt, Min, F)
@constraint(invmgmt, F >= sum{c[i,t] * p[i,t], i=1:I, t=1:T})

# Constraint: cannot exceed production limits
for i in 1:I, t in 1:T
    @constraint(invmgmt, p[i,t] >= 0)
    @constraint(invmgmt, p[i,t] <= P)
end
for i in 1:I
    @constraint(invmgmt, sum{p[i,t], t=1:T} <= Q)
end

# Constraint: cannot exceed inventory limits
for t in 1:T
    @constraint(invmgmt,
        v1 + sum{p[i,s], i=1:I, s=1:t} - sum{d[s],s=1:t} >= Vmin)
    @constraint(invmgmt,
        v1 + sum{p[i,s], i=1:I, s=1:t} - sum{d[s],s=1:t} <= Vmax)
end

# Solve
setsolver(invmgmt, GurobiSolver())
status = solve(invmgmt)

println(getobjectivevalue(invmgmt))


text:
using Revise #load this first
module MinVariance2
using Revise

#########################Packages#########################
using Juno, UnPack, Atom
using LinearAlgebra, Distributions, CSV, DataFrames, Statistics
using StatsBase, Random, CodecZlib, CodecZstd, XLSX, Dates

using Finometrics

#########################Globals#########################
#should only use these if signficant advantage to a parameter

#########################Constants#########################
const PARAMETER_FILE = "MVParameters"
const PARAM = Dict{Symbol,Any}()

#########################Files#########################
include("Parameters.jl")
include("Utilities.jl")

end # module

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
include("Misc/acb_root_ctx.jl")
include("Misc/Infinity.jl")
include("Misc/MSet.jl")
include("Misc/fmpr.jl")
include("Misc/FactoredElem.jl")
include("Misc/FactoredElemIntegers.jl")
include("Misc/Integer.jl")
include("Misc/Primes.jl")
include("Misc/coprime.jl")
include("Misc/Magma.jl")
include("Misc/Matrix.jl")
include("Misc/CRT.jl")  # for crt_env and modular_env
include("Misc/Poly.jl")
include("Misc/nmod_poly.jl")
include("Misc/MPoly.jl")
include("Misc/MPolyFactor.jl")
include("Misc/Residue.jl")
include("Misc/ResidueRingPoly.jl")
include("Misc/FiniteField.jl")
include("Misc/Arb.jl")
include("Misc/psibound.jl")
include("Misc/UnitsModM.jl")
include("Misc/toMagma.jl")
include("Misc/toNemo.jl")
include("Misc/Places.jl")
include("Misc/RatRecon.jl")
include("Misc/HNFstorjohann.jl")
include("Misc/RealField.jl")
include("Misc/serialize.jl")
include("Misc/Series.jl")
include("Misc/Screen.jl")
include("Misc/IJulia.jl")
include("Misc/Localization.jl")
include("Misc/OrdLocalization.jl")
include("Misc/jordan.jl")
include("Misc/PseudoPolynomial.jl")
include("Misc/UnsafeRational.jl")
include("Misc/FactoredElemWrapped.jl")
include("Misc/Fields.jl")
include("Misc/Cartesian.jl")
#include("Misc/RelFiniteField.jl")
#include("Misc/Plesken.jl")
#include("Misc/cohomology.jl")


text:
abstract type AbstractTween end

# ----------------------------------------------------------------------------
# The tween engine cannot directly change your objects attributes, since it doesn't
# about know them. Therefore, you need to let it know how to get and set the different
# attributes of your objects:
# **you need to implement the [AbstractTweenAccessor] interface for each object you will animate**.
abstract type AbstractTweenAccessor end
# Gets one or many values from the target object associated to the
# given tween type. It is used by the Tween Engine to determine starting
# values.
#
# [target] The target object of the tween.
# [tween] The tween interpolating on the target
# [returnValues] An array which should be modified by the tween.
# Returns a count of how many values are to be modified.
# getValues!(target::AbstractNode, tween::Tween, buffer::Array{Float64});

# This method is called by the Tween Engine each time a running tween
# associated with the current target object has been updated.
#
# [target] The target object of the tween.
# [tween] The tween interpolating on the target
# [tweenType] An integer representing the property index.
# [newValues] The new values determined by the Tween Engine.

# setValues!(target::AbstractNode, tween::Tween, Array{Float64} newValues)

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL