text:
const ignorefirst = 10 # cm
const bigturn = π/3 # 60°
# const smallturn = π/93 # 60°
const s = 500

const Point = SVector{2, Float64}
point(::Missing) = missing
point(x::Instantaneous)= Point(x.data[1], x.data[2])
point(x::Point) = x

_getv(spl, k) = SVector{2, Float64}(derivative(spl, k))

function gettpindex(spl, ks)
    tp = ks[1]
    vlast = _getv(spl, ks[1])
    for k in Iterators.drop(ks, 1)
        v = _getv(spl, k)
        Δ = angle(vlast, v)
        tp = k
        Δ > bigturn && break
        vlast = v
    end
    return tp
end

function gettpknot(spl)
    ks = Dierckx.get_knots(spl)
    filter!(k -> norm(spl(k) - spl(0)) > ignorefirst, ks)
    isempty(ks) && return 0.0
    tp2 = gettpindex(spl, ks)
    # return tp2
    tp1 = copy(tp2)
    for k in ks
        k == tp2 && break
        tp1 = k
    end
    tp1 += 0.1
    if tp1 < tp2
        main = _getv(spl, tp1)
        for t in tp2:-0.3:tp1
            v = _getv(spl, t)
            Δ = angle(main, v)
            Δ < bigturn && return t
        end
    end
    return tp2
end

mutable struct TimedPoint
    xy::Point
    t::Float64
end

const PointCollection = StructVector{TimedPoint}
pointcollection(x::Missing, t₀) = StructVector{TimedPoint}(undef, 0)
pointcollection(x, t₀) = StructVector(TimedPoint(Point(i[1], i[2]), i[3] - t₀) for i in eachrow(x.data))

mutable struct Track
    coords::Vector{Point}
    t::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
    tp::Int
    rawcoords::StructArray{TimedPoint}
end

function filterdance(xy, Δt)
    xy2 = [xy[1,:]]
    t = [0.0]
    for p in eachrow(xy)
        if norm(p - xy2[end]) > 4
            push!(xy2, p)
            push!(t, t[end] + Δt)
        else
            t[end] += Δt
        end
    end
    t .-= t[1]
    return t, hcat(xy2...)
end

function Track(x::Prolonged)
    xyt = !issorted(x.data[:, 3]) ? sortslices(x.data, dims = 1, lt = (x, y) -> isless(x[3], y[3])) : x.data
    Δt = mean(trim(diff(xyt[:, 3]), prop = 0.1))
    t, xy = filterdance(xyt[:,1:2], Δt)
    spl = ParametricSpline(t, xy; s = s, k = 2)
    tl = range(0.0, step = Δt, stop = t[end])
    xyl = Point.(spl.(tl))
    tp = gettpknot(spl)
    i = findfirst(≥(tp), tl)
    if isnothing(i)
        i = length(tl)
    end
    raw = pointcollection((data = xyt, ), xyt[1,3])
    Track(xyl, tl, i, raw)
end

homing(t::Track) = t.coords[1:t.tp]
searching(t::Track) = t.coords[t.tp:end]
searchcenter(t::Track) = mean(searching(t))
turningpoint(t::Track) = t.coords[t.tp]

mutable struct Common
    feeder::Union{Missing, Point}
    nest::Union{Missing, Point}
    track::Track
    pellet::PointCollection
    fictive_nest::Point
    pickup::Union{Missing, Point}
    dropoff::Point
end


text:
#######################
# Convergence Metrics #
#######################

# Utilities

abschange(curr, prev) = Float64(abs(curr - prev))
relchange(curr, prev) = Float64(abs(curr - prev)/abs(curr))

maxdiff(x::AbstractArray, y::AbstractArray) = mapreduce((a, b) -> abs(a - b), max, x, y)
abschange(curr::T, prev) where {T<:AbstractArray} = maxdiff(curr, prev)
relchange(curr::T, prev) where {T<:AbstractArray} = maxdiff(curr, prev)/maximum(abs, curr)

# Single value

"""
Absolute difference convergence for single objective optimization.

This convergence metric allows to estimate an absolute difference between consecutive
states of the optimization algorithm, and triggers convergence when,

- `|f(x) - f(x')| < ε`

where `ε` is a tolerance value, `x` and `x'` are previous and current minimizers
found by the optimization algorithm.
"""
mutable struct AbsDiff{T} <: ConvergenceMetric
    tol::T
    Δ::Float64
    value::T
end
AbsDiff(tol::T) where {T<:AbstractFloat} = AbsDiff(tol, Inf, zero(T))
AbsDiff() = AbsDiff(1e-12)
description(m::AbsDiff) = "|f(x) - f(x')|"
function assess!(m::AbsDiff, state::AbstractOptimizerState)
    val = value(state)
    m.Δ = abschange(val, m.value)
    m.value = val
    converged(m)
end

"""
Relative difference convergence metric for single objective optimization.

This convergence metric allows to estimate a relative difference between consecutive
states of the optimization algorithm, and triggers convergence when,

- `|f(x) - f(x')|/|f(x')| < ε`

where `ε` is a tolerance value, `x` and `x'` are previous and current minimizers
found by the optimization algorithm.
"""
mutable struct RelDiff{T} <: ConvergenceMetric
    tol::T
    Δ::Float64
    value::T
end
RelDiff(tol::T) where {T<:AbstractFloat} = RelDiff(tol, Inf, zero(T))
RelDiff() = RelDiff(1e-12)
description(m::RelDiff) = "|f(x) - f(x')|/|f(x')|"
function assess!(m::RelDiff, state::AbstractOptimizerState)
    val = value(state)
    m.Δ = relchange(val, m.value)
    m.value = val
    converged(m)
end

"""
    gd(A,R)

Calculate a generational distance between set `A` and the reference set `R`.
This metric measures the convergence, i.e. closeness of the non-dominated solutions
to the Pareto front, of a population.

*Note:* Parameters are column-major matrices.
"""
function gd(A::AbstractMatrix, R::AbstractMatrix)
    da, na = size(A)
    dr, nr = size(R)
    (na == 0 || nr == 0) && return Inf
    sum = 0
    for a in eachcol(A)
        sum += minimum(norm(a-r) for r in eachcol(R))
    end
    sum/na
end

"""
    igd(S,R)

Calculate an inverted generational distance, [`gd`](@ref), between set `S` and the reference set `R`.
Parameters are column-major matrices.
"""
igd(S,R) = gd(R,S)

"""
(Inverse) generational distance convergence metric for multi-objective optimization.

This convergence metric allows to estimate a (inverse) generational distance between
consecutive states of the multi-objective optimization algorithm, and triggers convergence when,

- ``GD(S,R) = \\frac{1}{|S|}\\left(\\sum_{s \\in S} \\min_{r \\in R}||f(s) - f(r)||^2\\right)^{1/2} < \\varepsilon``

where `ε` is a tolerance value, `S` and `R` are discrete Pareto front approximations
for previous and current state of the optimization algorithm.

The inverse generational distance calculated as `IGD(S,R) = GR(R,S)`.
To use, `IGD` metric pass `true` value to the metric constructor argument, e.g.

- `GD()` create a generational distance metric
- `GD(true)` create an inverse generational metric
"""
mutable struct GD{T} <: ConvergenceMetric
    tol::T
    Δ::Float64
    value::Matrix{T}
    inverted::Bool
end
GD(tol::T, inv=false) where {T<:AbstractFloat} = GD(tol, Inf, zeros(T,0,0), inv)
GD(inv=false) = GD(1e-5, inv)
function description(m::GD)
    prefix = m.inverted ? "I" : ""
    "|$(prefix)GD(P) - $(prefix)GD(P')|"
end
function assess!(m::GD, state::AbstractOptimizerState)
    val = value(state)
    m.Δ = m.inverted ? igd(val, m.value) : gd(val, m.value)
    m.value = val
    converged(m)
end

"""
    spread(S,R)

Returns a diversity metric of a population of set `S` to the reference set `R`.
"""
function spread(S::AbstractMatrix, R::AbstractMatrix)
    n = size(S,2)
    m = size(R,2)
    Δₖ = [ minimum(norm(view(S,:,i)-view(R,:,j)) for i in 1:n if view(S,:,i) != view(R,:,j)) for j in 1:m ]
    Δ = mean(Δₖ)
    sum(abs.(Δₖ.-Δ))/m*Δ
end

function spread(S::AbstractMatrix)
    n = size(S,2)
    n == 1 && return NaN
    Δₖ = [ minimum(norm(view(S,:,i)-view(S,:,j)) for j in 1:n if view(S,:,i) != view(S,:,j)) for i in 1:n ]
    Δ = mean(Δₖ)
    sum(abs.(Δₖ.-Δ))/n*Δ
end

##########################
# Convergence Assessment #
##########################

assess_convergence(state::AbstractOptimizerState, method) =
    any(assess!(cm, state) for cm in metrics(method))

ext: jl


text:
# # Interactive simple plot example
# this exmaple uses Interact to build a GUI to plot data parsed from
# CSV file to DataFrame variable.
# the documentation was writen using the Literate and Weave Pkg

# ### load needed Pkg:
# `CSV` - to read the CSV file and parse the data to DataFrame variable [`for more info see docs`](https://juliadata.github.io/CSV.jl/stable/)
#-
# `DataFrames` - Pkg to work with DataFrames [`for more info see docs`](https://juliadata.github.io/DataFrames.jl/stable/)
#-
# `Interact` - Pkg tio create the GUI [`for more info see docs`](https://juliagizmos.github.io/Interact.jl/latest/)
#-
# `Plots` - Pkg to plot the data [`for more info see docs`](http://docs.juliaplots.org/latest/)

using CSV, DataFrames, Interact, Plots

# ### define the Observables for the Interact GUI
loadbutton = filepicker() # file explorer GUI, the output is the file path
columnbuttons = Observable{Any}(dom"div"()) #
data = Observable{Any}(DataFrame) #
plt = Observable{Any}(plot()) #
dropdown_x = Observable{Any}(dom"div"()) #
dropdown_y = Observable{Any}(dom"div"()) #

# ### using the map!(function, destination, collection...) function
# in this case the map! takes the 'loadbutton' which holds the CSV file name
# and pipe it to the CSV.read function which parse the file into DataFrame
# variable data
map!(CSV.read, data, loadbutton)

"""
$(SIGNATURES)

this function build the buttons of all the DataFrame columns
#Inputs
    * df - parsed data from CSV file into DataFrame format
"""
function makebuttons(df)
    buttons = button.(string.(names(df)))
    for (btn, name) in zip(buttons, names(df))
        map!(t -> histogram(df[name]), plt, btn)
    end
    dom"div"(hbox(buttons))
end

# this map! each column in data to columnbuttons variable using the 'makebuttons'
# function defined earlier
map!(makebuttons, columnbuttons, data)

# create the ui interface witch has virtecly aranged:
# 1. loadbutton
# 2. columnbuttons
# 3. plt
ui = dom"div"(loadbutton, columnbuttons, plt)

# open the GUI window
w = Window()

# parse the ui into the window
body!(w, ui);

max_stars_repo_name: kobied/interactTests

max_issues_repo_name: kobied/interactTests

max_forks_repo_name: kobied/interactTests


text:
@testset "Utils" begin
    # Use a random seed
    p = Parameters(0.4, 1.4737, 7^3, 1e-5, Float64)
    @test p.seed == 0

    # And change the diameter of the particles
    s = SimulationSystem(p, 3, Float64)
    potential = PseudoHS(50.0, 1.5, p.kT)
    @test potential.σ == 1.5
    # And just check that everything works as expected
    initialize!(s)
    # As well as particles can be moved
    move!(1000, s, potential; interval = 500)
end

@testset "DataFrames" begin
    p = Parameters(0.4, 1.4737, 7^3, 1e-5, Float64)
    s = SimulationSystem(p, 3, Float64)
    potential = PseudoHS(50.0, p.kT)
    initialize!(s)

    move!(1000, s, potential; interval = 500)
    # Save positions and forces to DataFrames
    posdf, forcesdf = todataframe(s.system.positions, s.system.forces)

    # Check consistency for positions
    @test posdf.x == s.system.positions[:, 1]
    @test posdf.y == s.system.positions[:, 2]
    @test posdf.z == s.system.positions[:, 3]

    # Check consistency for forces
    @test forcesdf.x == s.system.forces[:, 1]
    @test forcesdf.y == s.system.forces[:, 2]
    @test forcesdf.z == s.system.forces[:, 3]
end


text:
using Random, Plots, DataFrames, ColorSchemes

import Plots: plot

"""
    plot(𝐼::Individuals)

Plot the initial and final positions as scatter plot in x,y plane.
"""
function plot(𝐼::Individuals)
    🔴_by_t = groupby(𝐼.🔴, :t)
    if (sum(names(🔴_by_t).=="lon")==0)
        fig=scatter(🔴_by_t[1].x,🔴_by_t[1].y,c=:red,label="t0",marker = (:circle, stroke(0)))
        scatter!(🔴_by_t[end].x,🔴_by_t[end].y,c=:blue,label="t1",marker = (:circle, stroke(0)))
    else
        fig=scatter(🔴_by_t[1].lon,🔴_by_t[1].lat,c=:red,label="t0",marker = (:circle, stroke(0)))
        scatter!(🔴_by_t[end].lon,🔴_by_t[end].lat,c=:blue,label="t1",marker = (:circle, stroke(0)))
    end
    return fig
end

"""
    plot_paths(df::DataFrame,nn::Integer,dMax::Float64=0.)

Plot random subset of size nn trajectories / paths.
"""
function plot_paths(df::DataFrame,nn::Integer,dMax::Float64=0.)
   IDs = randperm(maximum(df.ID))
   COs=["w" "y" "g" "k"]

   plt=plot(leg=false)
   df_by_ID = groupby(df, :ID)
   ID_list=[df_by_ID[i][1,:ID] for i in 1:length(df_by_ID)]
   for ii=1:nn
      jj=findall(ID_list.==IDs[ii])[1]
      tmp=df_by_ID[jj]
      if dMax > 0.
         d=abs.(diff(tmp[!,:lon]))
         jj=findall(d .> dMax)
         tmp[jj,:lon].=NaN; tmp[jj,:lat].=NaN
         d=abs.(diff(tmp[!,:lat]))
         jj=findall(d .> dMax)
         tmp[jj,:lon].=NaN; tmp[jj,:lat].=NaN
      end
      CO=COs[mod(ii,4)+1]

      hasproperty(df,:z) ? plot!(tmp[!,:lon],tmp[!,:lat],tmp[!,:z],linewidth=0.3) : nothing
      !hasproperty(df,:z) ? plot!(tmp[!,:lon],tmp[!,:lat],linewidth=0.3) : nothing
      #plot!(tmp[!,:lon],tmp[!,:lat],tmp[!,:z],linewidth=0.3)
   end
   return plt
end

"""
    scatter_zcolor(df,zc; cam=(0, 90))

```
df=groupby(𝐼.🔴, :t)
scatter_zcolor(df[end],df[end].z)
```
"""
function scatter_zcolor(df,zc; cam=(0, 90))
    lo=extrema(df.lon); lo=round.(lo).+(-5.0,5.0)
    la=extrema(df.lat); la=round.(la).+(-5.0,5.0)
    de=extrema(zc); de=round.(de).+(-5.0,5.0)

    scatter(df.lon,df.lat,zc,zcolor = zc,
    markersize=2,markerstrokewidth=0.1,camera = cam,
    xlims=lo,ylims=la,zlims=de,clims=de,leg=false)
end

"""
    scatter_movie(𝐼; cam=(0, 90))

Animation using `scatter_zcolor()
```
𝐼,Γ=example3("OCCA", lon_rng=(-165.0,-155.0),lat_rng=(25.0,35.0), z_init=5.5,)
scatter_movie(𝐼,cam=(70, 70))
```
"""
function scatter_movie(𝐼; cam=(0, 90))
    df=groupby(𝐼.🔴, :t)
    return @gif for t in 1:length(df)
        scatter_zcolor(df[t],df[t].z;cam=cam)
   end
end

"""
    phi_scatter(ϕ,df)

```
phi_scatter(ϕ,df)
```
"""
function phi_scatter(ϕ,df)
    nx,ny=size(ϕ)
    contourf(-0.5 .+ (1:nx),-0.5 .+ (1:ny),
             transpose(ϕ),c = :blues,linewidth = 0.1)
    scatter!(df.x,df.y,markersize=4.0,c=:red,marker = (:circle, stroke(0)),
             xlims=(0,nx),ylims=(0,ny),leg=:none)
end

"""
    OceanDepthLog()

Compute Ocean depth logarithm on regular grid.
"""
function OceanDepthLog(Γ)
    lon=[i for i=20.:2.0:380., j=-79.:2.0:89.]
    lat=[j for i=20.:2.0:380., j=-79.:2.0:89.]
    DL=interp_to_lonlat(Γ["Depth"],Γ,lon,lat)
    DL[findall(DL.<0)].=0
    DL=transpose(log10.(DL))
    DL[findall((!isfinite).(DL))].=NaN
    return (lon=lon[:,1],lat=lat[1,:],fld=DL,rng=(1.5,5))
end

"""
    map(𝐼::Individuals,background::NamedTuple)

Plot initial and final positions, superimposed on a map of ocean depth log.
"""
function map(𝐼::Individuals,𝐵::NamedTuple)
    xlims=extrema(𝐵.lon)
    ylims=extrema(𝐵.lat)
    plt=contourf(𝐵.lon,𝐵.lat,𝐵.fld,clims=𝐵.rng,c = :ice,
    colorbar=false, xlims=xlims,ylims=ylims)

    🔴_by_t = groupby(𝐼.🔴, :t)
    lo=deepcopy(🔴_by_t[1].lon); lo[findall(lo.<xlims[1])]=lo[findall(lo.<xlims[1])].+360
    scatter!(lo,🔴_by_t[1].lat,markersize=1.5,c=:red,leg=:none,marker = (:circle, stroke(0)))
    lo=deepcopy(🔴_by_t[end].lon); lo[findall(lo.<xlims[1])]=lo[findall(lo.<xlims[1])].+360
    scatter!(lo,🔴_by_t[end].lat,markersize=1.5,c=:yellow,leg=:none,marker = (:dot, stroke(0)))

    return plt
end


text:
# Test Benders Algorithm for Mixed Integer Linear Programming
# Edward J. Xu, edxu96@outlook.com
# March 16th, 2019
# ---------------------------------------------------------------------------
push!(LOAD_PATH, "$(homedir())/Desktop/OptiGas, DTU42136")
cd("$(homedir())/Desktop/OptiGas, DTU42136")
using BendersMilp_EDXU
# ---------------------------------------------------------------------------
# Test Example 1: mat_a and mat_b being column vector
n_x = 1
n_y = 1
vec_max_y = hcat([10])
vec_c = hcat([5])
vec_f = hcat([-3])
vec_b = hcat([4; 0; -13])
mat_a = hcat([1; 2; 1])
mat_b = hcat([2; -1; -3])
BendersMilp(n_x = n_x,
            n_y = n_y,
            vec_max_y = vec_max_y,
            vec_c = vec_c,
            vec_f = vec_f,
            vec_b = vec_b,
            mat_a = mat_a,
            mat_b = mat_b,
            epsilon = 0,
            timesIterationMax = 5)
# ---------------------------------------------------------------------------
# Test Example 2: mat_a and mat_b being matrix
n_x_2 = 2
n_y_2 = 2
vec_max_y_2 = hcat([10; 10])
vec_c_2 = hcat([5; 3])
vec_f_2 = hcat([-3; 1])
vec_b_2 = hcat([4; 0; -13])
mat_a_2 = [1 3; 2 1; 1 -5]
mat_b_2 = [2 -4; -1 2; -3 1]
BendersMilp(n_x = n_x_2,
            n_y = n_y_2,
            vec_max_y = vec_max_y_2,
            vec_c = vec_c_2,
            vec_f = vec_f_2,
            vec_b = vec_b_2,
            mat_a = mat_a_2,
            mat_b = mat_b_2,
            epsilon = 0.000001,
            timesIterationMax = 10)
# ---------------------------------------------------------------------------
# Test Example 3: sub and ray problem
n_x_3 = 2
n_y_3 = 2
vec_max_y_3 = hcat([2; 2])
vec_c_3 = hcat([2; 6])
vec_f_3 = hcat([2; 3])
vec_b_3 = hcat([5; 4])
mat_a_3 = [-1 2; 1 -3]
mat_b_3 = [3 -1; 2 2]
BendersMilp(n_x = n_x_3,
            n_y = n_y_3,
            vec_max_y = vec_max_y_3,
            vec_c = vec_c_3,
            vec_f = vec_f_3,
            vec_b = vec_b_3,
            mat_a = mat_a_3,
            mat_b = mat_b_3,
            epsilon = 0.000001,
            timesIterationMax = 10)
# ---------------------------------------------------------------------------
# Test Example 4: what if when infeasibility in x
n_x_4 = 2
n_y_4 = 2
vec_max_y_4 = hcat([2; 2])
vec_c_4 = hcat([2; 6])
vec_f_4 = hcat([2; 3])
vec_b_4 = hcat([5; 4; 1])
mat_a_4 = [-1 2; 1 -3; -1 -1]
mat_b_4 = [3 -1; 2 2; 0 0]
BendersMilp(n_x = n_x_4,
            n_y = n_y_4,
            vec_max_y = vec_max_y_4,
            vec_c = vec_c_4,
            vec_f = vec_f_4,
            vec_b = vec_b_4,
            mat_a = mat_a_4,
            mat_b = mat_b_4,
            epsilon = 0.000001,
            timesIterationMax = 10)
# ---------------------------------------------------------------------------
# Test Example 5: what if when infeasibility in y
n_x_5 = 2
n_y_5 = 2
vec_max_y_5 = hcat([2; 2])
vec_c_5 = hcat([2; 6])
vec_f_5 = hcat([2; 3])
vec_b_5 = hcat([5; 4; 5])
mat_a_5 = [-1 2; 1 -3; 0 0]
mat_b_5 = [3 -1; 2 2; 1 1]
BendersMilp(n_x = n_x_5,
            n_y = n_y_5,
            vec_max_y = vec_max_y_5,
            vec_c = vec_c_5,
            vec_f = vec_f_5,
            vec_b = vec_b_5,
            mat_a = mat_a_5,
            mat_b = mat_b_5,
            epsilon = 0.000001,
            timesIterationMax = 10)


text:
module CuDiff

sqr(x) = (x*x)
cube(x) = (x*x*x)
inverse(x) = one(x) / x

include("./dual.jl")
# include("./single.jl")

function derivative(f, x, p...)
    dual = f(Dual(x, one(x)), p...)
    return dual.u, dual.du
end

function evaluate(f, x, p...)
    res = f(Single(x), p...)
    return map(x -> isa(Dual) ? x.u : x, res)
end

export derivative, evaluate

end # module


text:
using LinearAlgebra
using CompScienceMeshes
using SauterSchwab3D
using StaticArrays

const pI   = point(1,0,0)
const pII  = point(0,1,0)
const pIII = point(0,0,0)
const pIV  = point(0,0,1)
const qIV  = point(0,0,-1)

const P = simplex(pI,pII,pIV,pIII)
const Q = simplex(pI,pIII,pII)

const sing = SauterSchwab3D.singularity_detection(P,Q)

Accuracy2 = 18
cf_ref = CommonFace5D(sing,SauterSchwab3D._legendre(Accuracy2,0.0,1.0))

function integrand(x,y)
      return ((x-pI)'*(y-pIV))*exp(-im*1*norm(x-y))/(4pi*norm(x-y))
end

function INTEGRAND(u,v)
   j1 = volume(P) * factorial(dimension(P))
   j2 = volume(Q) * factorial(dimension(Q))
   x = barytocart(P,u)
   y = barytocart(Q,v)
   output = integrand(x,y)*j1*j2

   return output
end

print("Ref: ")
ref = sauterschwab_parameterized(INTEGRAND, cf_ref)
println(ref)
println()

res_tp =[]
res_sp =[]
res_gm =[]
n1 = []
n2 = []
n3 = []

for i in 2:1:15
   Accuracy = i
   cf = CommonFace5D(sing,SauterSchwab3D._legendre(Accuracy,0.0,1.0))

   int_tp = sauterschwab_parameterized(INTEGRAND, cf)

   num_pts = 9*length(cf.qps)^5
   push!(n1,num_pts)

   push!(res_tp,int_tp)
end

for i in 2:1:7
   Accuracy = i
   cf_s = CommonFace5D_S(sing,(SauterSchwab3D._legendre(Accuracy,0.0,1.0),
                          SauterSchwab3D._shunnham2D(Accuracy),
                          SauterSchwab3D._shunnham3D(Accuracy)))

   int_sp = sauterschwab_parameterized(INTEGRAND, cf_s)

   num_pts = length(cf_s.qps[1])*length(cf_s.qps[2])^2+8*length(cf_s.qps[2])*length(cf_s.qps[3])
   push!(n2,num_pts)
   push!(res_sp,int_sp)
end

for i in 2:1:12
   Accuracy = i
   cf_gm = CommonFace5D_S(sing,(SauterSchwab3D._legendre(Accuracy,0.0,1.0),
                          SauterSchwab3D._grundmannMoeller2D(2*Accuracy-1),
                          SauterSchwab3D._grundmannMoeller3D(2*Accuracy-1)))

   int_gm = sauterschwab_parameterized(INTEGRAND, cf_gm)

   num_pts = length(cf_gm.qps[1])*length(cf_gm.qps[2])^2+8*length(cf_gm.qps[2])*length(cf_gm.qps[3])
   push!(n3,num_pts)
   push!(res_gm,int_gm)
end

err_tp = norm.(res_tp.-ref)/norm(ref)
err_sp = norm.(res_sp.-ref)/norm(ref)
err_gm = norm.(res_gm.-ref)/norm(ref)

println(err_tp)
println(err_sp)
println(err_gm)

using Plots
plot( yaxis=:log, xaxis=:log, fontfamily="Times")
plot!(n1,err_tp, label="Gauss Tensor-Product",markershape=:circle)
plot!(n2,err_sp, label="Simplex Tensor-Product",markershape=:rect)
#plot!(n3,err_gm, label="Simplex-Product GM",markershape=:x)
plot!(xlims=(1e2,1e5),ylims=(1e-7,1))
plot!(xlabel="#Quad. pts/Func. evals", ylabel="Rel. Error.", title="Common Face 5D",legend=:bottomleft)

using BenchmarkTools

ref = -0.00463274359881397 + 0.002581200521272345im

cf = CommonFace5D(sing,SauterSchwab3D._legendre(4,0.0,1.0))

int_tp = sauterschwab_parameterized(INTEGRAND, cf)

num_pts = 9*length(cf.qps)^5
println("#Pts: ",num_pts)
err_tp = norm.(int_tp.-ref)/norm(ref)
println("Rel. Err: ",err_tp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, cf) end

cf_s = CommonFace5D_S(sing,(SauterSchwab3D._legendre(4,0.0,1.0),
                       SauterSchwab3D._shunnham2D(4),
                       SauterSchwab3D._shunnham3D(4)))

int_sp = sauterschwab_parameterized(INTEGRAND, cf_s)

num_pts = length(cf_s.qps[1])*length(cf_s.qps[2])^2+8*length(cf_s.qps[2])*length(cf_s.qps[3])

println("#Pts: ",num_pts)
err_sp = norm.(int_sp.-ref)/norm(ref)
println("Rel Err: ",err_sp)
@time  for i in 1:100  sauterschwab_parameterized(INTEGRAND, cf_s) end

#=
cf_gm = CommonFace5D_S((SauterSchwab3D._legendre(4,0.0,1.0),
                        SauterSchwab3D._grundmannMoeller2D(2*4-1),
                        SauterSchwab3D._grundmannMoeller3D(2*4-1)))

int_gm = sauterschwab_parameterized(INTEGRAND, cf_gm)

num_pts = length(cf_gm.qps[1])*length(cf_gm.qps[2])^2+8*length(cf_gm.qps[2])*length(cf_gm.qps[3])

println("#Pts: ",num_pts)
err_gm = norm.(int_gm.-ref)/norm(ref)
println("Rel Err: ",err_gm)
=#


text:
# This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT).

# Note: Most (if not all) of these types are temporary, to be replaced later
# on by generic types from RadiationDetectorSignals.jl (when the latter is
# mature enough).

# Internal units should be SI units
const internal_length_unit  = ConstructiveSolidGeometry.internal_length_unit
const internal_angle_unit   = ConstructiveSolidGeometry.internal_angle_unit
const internal_time_unit    = u"s"
const internal_voltage_unit = u"V"
const internal_energy_unit  = u"eV"
const internal_efield_unit  = internal_voltage_unit / internal_length_unit
const internal_charge_unit  = u"C"
const external_charge_unit  = u"e_au" # elementary charge - from UnitfulAtomic.jl

to_internal_units(x::Quantity{<:Real, dimension(internal_time_unit)})    = ustrip(uconvert(internal_time_unit,    x))
to_internal_units(x::Quantity{<:Real, dimension(internal_voltage_unit)}) = ustrip(uconvert(internal_voltage_unit, x))
to_internal_units(x::Quantity{<:Real, dimension(internal_efield_unit)})   = ustrip(uconvert(internal_efield_unit,  x))
to_internal_units(x::Quantity{<:Real, dimension(internal_energy_unit)})  = ustrip(uconvert(internal_energy_unit,  x))
to_internal_units(x::Quantity{<:Real, dimension(internal_charge_unit)})  = ustrip(uconvert(internal_charge_unit,  x))

from_internal_units(x::Real, unit::Unitful.Units{<:Any, dimension(internal_time_unit)})    = uconvert(unit, x * internal_time_unit)
from_internal_units(x::Real, unit::Unitful.Units{<:Any, dimension(internal_voltage_unit)}) = uconvert(unit, x * internal_voltage_unit)
from_internal_units(x::Real, unit::Unitful.Units{<:Any, dimension(internal_efield_unit)})  = uconvert(unit, x * internal_efield_unit)
from_internal_units(x::Real, unit::Unitful.Units{<:Any, dimension(internal_energy_unit)})  = uconvert(unit, x * internal_energy_unit)
from_internal_units(x::Real, unit::Unitful.Units{<:Any, dimension(internal_charge_unit)})  = uconvert(unit, x * internal_charge_unit)

# Internal function for now: We should also fano noise here (optionally)
_convert_internal_energy_to_external_charge(material) = inv(to_internal_units(material.E_ionisation)) * external_charge_unit

unit_conversion = Dict{String, Unitful.Units}(
    "nm" => u"nm", "um" => u"μm", "mm" => u"mm", "cm" => u"cm", "m" => u"m", #length
    "deg" => u"°","rad" => u"rad", #angle
    "V" => u"V", "kV" => u"kV", #potential
    "K" => u"K", "Kelvin" => u"K", "C" => u"°C", "Celsius" => u"°C", #temperature
)

const UnitTuple = NamedTuple{(:length, :angle, :potential, :temperature), Tuple{
                    Unitful.FreeUnits{<:Any, Unitful.𝐋},
                    Unitful.FreeUnits{<:Any, NoDims},
                    Unitful.FreeUnits{<:Any, Unitful.𝐋^2 * Unitful.𝐌 * Unitful.𝐈^-1 * Unitful.𝐓^-3},
                    Unitful.FreeUnits{<:Any, Unitful.𝚯}}}

function default_unit_tuple()::UnitTuple
    return (
        length = u"m",
        angle = u"°",
        potential = u"V",
        temperature = u"K"
    )
end

function construct_unit(ustring::String)::Unitful.Units
    try
        parsed_unit = uparse(ustring)
        @assert parsed_unit isa Unitful.Units
        return parsed_unit
    catch
        @assert ustring in keys(unit_conversion) "Unit string $(ustring) cannot be interpreted..."
        return unit_conversion[ustring]
    end
end

function construct_units(config_file_dict::AbstractDict)::UnitTuple
    dunits::NamedTuple = default_unit_tuple()
    if haskey(config_file_dict, "units")
        d = config_file_dict["units"]
        dunits = (
            length = haskey(d, "length") ? construct_unit(d["length"]) : dunits.length,
            angle  = haskey(d, "angle") ? construct_unit(d["angle"]) : dunits.angle,
            potential = haskey(d, "potential") ? construct_unit(d["potential"]) : dunits.potential,
            temperature = haskey(d, "temperature") ? construct_unit(d["temperature"]) : dunits.temperature
        )
    end
    matching_units::Bool = dimension.(values(dunits)) == (Unitful.𝐋, Unitful.NoDims, Unitful.𝐋^2 * Unitful.𝐌 * Unitful.𝐈^-1 * Unitful.𝐓^-3, Unitful.𝚯)
    @assert matching_units "Some of the units are not given in the right format:"*
        (dimension(dunits.length)      == Unitful.𝐋       ? "" : "\n  - units.length ($(dunits.length)) does not describe a length.")*
        (dimension(dunits.angle)       == Unitful.NoDims  ? "" : "\n  - units.angle ($(dunits.angle)) does not describe an angle.")*
        (dimension(dunits.potential)   == Unitful.𝐋^2 * Unitful.𝐌 * Unitful.𝐈^-1 * Unitful.𝐓^-3 ? "" : "\n  - units.potential ($(dunits.potential)) does not describe a potential.")*
        (dimension(dunits.temperature) == Unitful.𝚯       ? "" : "\n  - units.temperature ($(dunits.temperature)) does not describe a temperature.")

    return UnitTuple(dunits)
end

const MaybeWithUnits{T} = Union{T, Quantity{<:T}}
const RealQuantity = MaybeWithUnits{<:Real}

DetectorHitEvents = TypedTables.Table{
    <:NamedTuple{
        (:evtno, :detno, :thit, :edep, :pos),
        <:Tuple{
            Union{Integer, AbstractVector{<:Integer}},
            Union{Integer, AbstractVector{<:Integer}},
            AbstractVector{<:RealQuantity},
            AbstractVector{<:RealQuantity},
            AbstractVector{<:AbstractVector{<:RealQuantity}}
        }
    }
}


text:
module BumpOnTailSimulation

using Plots, LaTeXStrings
using ReducedBasisMethods
using Particles
using Particles.BumpOnTail
using Random

# parameters
const dt = 1e-1             # timestep
const T = 25                # final time
const nt = Int(div(T, dt))  # nb. of timesteps
const np = Int(5e3)         # nb. of particles
const nh = 16               # nb. of elements
const p = 3                 # spline degree

const nₚ₁ = 10
const nₚ₂ = 1
const nₚ₃ = 1
const nₚ₄ = 1
const nₚ₅ = 1
const nparam = nₚ₁ * nₚ₂ * nₚ₃ * nₚ₄ * nₚ₅

const vmax = +10
const vmin = -10

# bump-on tail instability - reference values: κ, ε, a, v₀, σ
params = (κ = 0.3,      # spatial perturbation wave number
          ε = 0.03,     # amplitude of spatial perturbation
          a = 0.1,      # fast particle share
          v₀= 4.5,      # velocity
          σ = 0.5,      # temperature
          χ = 1.0)

# sampling parameters
κ = Parameter(:κ,  0.1,  0.5,  nₚ₁)
ε = Parameter(:ε,  0.03, 0.03, nₚ₂)
a = Parameter(:a,  0.1,  0.1,  nₚ₃)
v₀= Parameter(:v₀, 4.5,  4.5,  nₚ₄)
σ = Parameter(:σ,  0.5,  0.5,  nₚ₅)

const pspace = ParameterSpace(κ, ε, a, v₀, σ)

const L = 2π/params.κ         # domain length
# const h = L/nh                # element width

function run()

    # integrator parameters
    IP = VPIntegratorParameters(dt, nt, nt+1, nh, np)

    # integrator cache
    IC = VPIntegratorCache(IP)

    # B-spline Poisson solver
    poisson = PoissonSolverPBSplines(p, nh, L)

    # set random generator seed
    Random.seed!(1234)

    # initial data
    particles = BumpOnTail.draw_accept_reject(np, params)
    # particles = BumpOnTail.draw_importance_sampling(np, params)

    # training set
    TS = TrainingSet(nt+1, poisson, particles, pspace)

    # loop over parameter set
    for p in eachindex(pspace)

        # get parameter tuple
        lparams = merge(pspace(p), (χ = pspace[p].κ / params.κ,))

        # integrate particles for parameter
        integrate_vp!(particles, poisson, lparams, IP, IC; save=true, given_phi=false)

        # copy solution
        TS.X[1,:,:,p] .= IC.X
        TS.V[1,:,:,p] .= IC.V
        TS.A[1,:,:,p] .= IC.A
        TS.Φ[1,:,:,p] .= IC.Φ

        # copy diagnostics
        TS.W[:,p] .= IC.W
        TS.K[:,p] .= IC.K
        TS.M[:,p] .= IC.M
    end

    # save results to HDF5
    h5save("../runs/BoT_Np5e4_k_010_050_np_10_T25.h5", TS, IP, poisson, pspace, params)

    # plot
    plot(IP.t, TS.W, linewidth = 2, xlabel = L"$n_t$", yscale = :log10, legend = :none,
        grid = true, gridalpha = 0.5, minorgrid = true, minorgridalpha = 0.2)
    savefig("../runs/BoT_Np5e4_k_010_050_np_10_T25_plot1.pdf")
    # TODO: Change filename to something meaningful!

    #
    α, β = get_regression_αβ(IP.t, TS.W, 2)

    #
    Wₗᵢₙ = zero(TS.W)
    for i in axes(Wₗᵢₙ,2)
        Wₗᵢₙ[:,i] .= exp.(α[i] .+ β[i] .* IP.t)
    end

    # plot
    plot(xlabel = L"$n_t$", yscale = :log10, ylims = (1E-3,1E1), legend = :none,
        grid = true, gridalpha = 0.5, minorgrid = true, minorgridalpha = 0.2)
    plot!(IP.t, TS.W[:,1:5], linewidth = 2, alpha = 0.25)
    plot!(IP.t, Wₗᵢₙ[:,1:5], linewidth = 2, alpha = 0.5)
    savefig("../runs/BoT_Np5e4_k_010_050_np_10_T25_plot2.pdf")
    # TODO: Change filename to something meaningful!

    # plot
    plot(xlabel = L"$n_t$", yscale = :log10, ylims = (1E-3,1E1), legend = :none,
        grid = true, gridalpha = 0.5, minorgrid = true, minorgridalpha = 0.2)
    plot!(IP.t, TS.W[:,6:10], linewidth = 2, alpha = 0.25)
    plot!(IP.t, Wₗᵢₙ[:,6:10], linewidth = 2, alpha = 0.5)
    savefig("../runs/BoT_Np5e4_k_010_050_np_10_T25_plot3.pdf")
    # TODO: Change filename to something meaningful!

end

end

using .BumpOnTailSimulation

BumpOnTailSimulation.run()


text:
const IterationResamplingTypes =
    Union{Holdout,Nothing,MLJBase.TrainTestPairs}

## TYPES AND CONSTRUCTOR

mutable struct DeterministicIteratedModel{M<:Deterministic} <: MLJBase.Deterministic
    model::M
    controls
    resampling # resampling strategy
    measure
    weights::Union{Nothing,Vector{<:Real}}
    class_weights::Union{Nothing,Dict{Any,<:Real}}
    operation
    retrain::Bool
    check_measure::Bool
    iteration_parameter::Union{Nothing,Symbol,Expr}
    cache::Bool
end

mutable struct ProbabilisticIteratedModel{M<:Probabilistic} <: MLJBase.Probabilistic
    model::M
    controls
    resampling # resampling strategy
    measure
    weights::Union{Nothing,AbstractVector{<:Real}}
    class_weights::Union{Nothing,Dict{Any,<:Real}}
    operation
    retrain::Bool
    check_measure::Bool
    iteration_parameter::Union{Nothing,Symbol,Expr}
    cache::Bool
end

const ERR_MISSING_TRAINING_CONTROL =
    ArgumentError("At least one control must be a training control "*
                  "(have type `$TrainingControl`) or be a "*
                  "custom control that calls IterationControl.train!. ")

const ERR_TOO_MANY_ARGUMENTS =
    ArgumentError("At most one non-keyword argument allowed. ")
const EitherIteratedModel{M} =
    Union{DeterministicIteratedModel{M},ProbabilisticIteratedModel{M}}
const ERR_NOT_SUPERVISED =
    ArgumentError("Only `Deterministic` and `Probabilistic` "*
                  "model types supported.")
const ERR_NEED_MEASURE =
    ArgumentError("Unable to deduce a default measure for specified model. "*
                  "You must specify `measure=...`. ")
const ERR_NEED_PARAMETER =
    ArgumentError("Unable to determine the name of your model's iteration "*
                  "parameter. Please specify `iteration_parameter=...`. This "*
                  "must be a `Symbol` or, in the case of a nested parameter, "*
                  "an `Expr` (as in `booster.nrounds`). ")
const ERR_MODEL_UNSPECIFIED = ArgumentError(
"Expecting atomic model as argument, or as keyword argument `model=...`, "*
    "but neither detected. ")

const WARN_POOR_RESAMPLING_CHOICE =
    "Training could be very slow unless "*
    "`resampling` is `Holdout(...)`, `nothing`, or "*
    "a vector of the form `[(train, test),]`, where `train` and `test` "*
    "are valid row indices for the data, as in "*
    "`resampling = [(1:100, 101:150),]`. "

const WARN_POOR_CHOICE_OF_PAIRS =
    "Training could be very slow unless you limit the number of `(train, test)` pairs "*
    "to one, as in resampling = [(1:100, 101:150),]. Alternatively, "*
    "use a `Holdout` resampling strategy. "

err_bad_iteration_parameter(p) =
    ArgumentError("Model to be iterated does not have :($p) as an iteration parameter. ")

"""
    IteratedModel(model=nothing,
                  controls=$CONTROLS_DEFAULT,
                  retrain=false,
                  resampling=Holdout(),
                  measure=nothing,
                  weights=nothing,
                  class_weights=nothing,
                  operation=predict,
                  verbosity=1,
                  check_measure=true,
                  iteration_parameter=nothing,
                  cache=true)

Wrap the specified `model <: Supervised` in the specified iteration
`controls`. Training a machine bound to the wrapper iterates a
corresonding machine bound to `model`. Here `model` should support
iteration.

To list all controls, do `MLJIteration.CONTROLS`. Controls are
summarized at
[https://alan-turing-institute.github.io/MLJ.jl/dev/getting_started/](https://alan-turing-institute.github.io/MLJ.jl/dev/controlling_iterative_models/)
but query individual doc-strings for details and advanced options. For
creating your own controls, refer to the documentation just cited.

To make out-of-sample losses available to the controls, the machine
bound to `model` is only trained on part of the data, as iteration
proceeds.  See details on training below. Specify `retrain=true`
to ensure the model is retrained on *all* available data, using the
same number of iterations, once controlled iteration has stopped.

Specify `resampling=nothing` if all data is to be used for controlled
iteration, with each out-of-sample loss replaced by the most recent
training loss, assuming this is made available by the model
(`supports_training_losses(model) == true`). Otherwise, `resampling`
must have type `Holdout` (eg, `Holdout(fraction_train=0.8, rng=123)`).

Assuming `retrain=true` or `resampling=nothing`,
`iterated_model` behaves exactly like the original `model` but with
the iteration parameter automatically selected. If
`retrain=false` (default) and `resampling` is not `nothing`, then
`iterated_model` behaves like the original model trained on a subset
of the provided data.

Controlled iteration can be continued with new `fit!` calls (warm
restart) by mutating a control, or by mutating the iteration parameter
of `model`, which is otherwise ignored.

### Training

Given an instance `iterated_model` of `IteratedModel`, calling
`fit!(mach)` on a machine `mach = machine(iterated_model, data...)`
performs the following actions:

- Assuming `resampling !== nothing`, the `data` is split into *train* and
  *test* sets, according to the specified `resampling` strategy, which
  must have type `Holdout`.

- A clone of the wrapped model, `iterated_model.model`, is bound to
  the train data in an internal machine, `train_mach`. If `resampling
  === nothing`, all data is used instead. This machine is the object
  to which controls are applied. For example, `Callback(fitted_params
  |> print)` will print the value of `fitted_params(train_mach)`.

- The iteration parameter of the clone is set to `0`.

- The specified `controls` are repeatedly applied to `train_mach` in
  sequence, until one of the controls triggers a stop. Loss-based
  controls (eg, `Patience()`, `GL()`, `Threshold(0.001)`) use an
  out-of-sample loss, obtained by applying `measure` to predictions
  and the test target values. (Specifically, these predictions are
  those returned by `operation(train_mach)`.)  If `resampling ===
  nothing` then the most recent training loss is used instead. Some
  controls require *both* out-of-sample and training losses (eg,
  `PQ()`).

- Once a stop has been triggered, a clone of `model` is bound to all
  `data` in a machine called `mach_production` below, unless
  `retrain == false` or `resampling === nothing`, in which case
  `mach_production` coincides with `train_mach`.

### Prediction

Calling `predict(mach, Xnew)` returns `predict(mach_production,
Xnew)`. Similar similar statements hold for `predict_mean`,
`predict_mode`, `predict_median`.

### Controls

A control is permitted to mutate the fields (hyper-parameters) of
`train_mach.model` (the clone of `model`). For example, to mutate a
learning rate one might use the control

    Callback(mach -> mach.model.eta = 1.05*mach.model.eta)

However, unless `model` supports warm restarts with respect to changes
in that parameter, this will trigger retraining of `train_mach` from
scratch, with a different training outcome, which is not recommended.

### Warm restarts

If `iterated_model` is mutated and `fit!(mach)` is called again, then
a warm restart is attempted if the only parameters to change are
`model` or `controls` or both.

Specifically, `train_mach.model` is mutated to match the current value
of `iterated_model.model` and the iteration parameter of the latter is
updated to the last value used in the preceding `fit!(mach)` call. Then
repeated application of the (updated) controls begin anew.

"""
function IteratedModel(args...;
                       model=nothing,
                       control=CONTROLS_DEFAULT,
                       controls=control,
                       resampling=MLJBase.Holdout(),
                       measures=nothing,
                       measure=measures,
                       weights=nothing,
                       class_weights=nothing,
                       operation=predict,
                       retrain=false,
                       check_measure=true,
                       iteration_parameter=nothing,
                       cache=true)

    length(args) < 2 || throw(ArgumentError("At most one non-keyword argument allowed. "))
    if length(args) === 1
        atom = first(args)
        model === nothing ||
            @warn "Using `model=$atom`. Ignoring specification `model=$model`. "
    else
        model === nothing && throw(ERR_MODEL_UNSPECIFIED)
        atom = model
    end

    if atom isa Deterministic
        iterated_model = DeterministicIteratedModel(atom,
                                                    controls,
                                                    resampling,
                                                    measure,
                                                    weights,
                                                    class_weights,
                                                    operation,
                                                    retrain,
                                                    check_measure,
                                                    iteration_parameter,
                                                    cache)
    elseif atom isa Probabilistic
        iterated_model = ProbabilisticIteratedModel(atom,
                                                    controls,
                                                    resampling,
                                                    measure,
                                                    weights,
                                                    class_weights,
                                                    operation,
                                                    retrain,
                                                    check_measure,
                                                    iteration_parameter,
                                                    cache)
    else
        throw(ERR_NOT_SUPERVISED)
    end

    message = clean!(iterated_model)
    isempty(message) || @warn message

    return iterated_model

end

function MLJBase.clean!(iterated_model::EitherIteratedModel)
    message = ""
    measure = iterated_model.measure
    if measure === nothing &&
        iterated_model.resampling !== nothing
        measure = MLJBase.default_measure(iterated_model.model)
        measure === nothing && throw(ERR_NEED_MEASURE)
    end
    iter = deepcopy(iterated_model.iteration_parameter)
    if iter === nothing
        iter = iteration_parameter(iterated_model.model)
        iter === nothing && throw(ERR_NEED_PARAMETER)
    end
    try
        MLJBase.recursive_getproperty(iterated_model.model,
                                      iter)
    catch
        throw(err_bad_iteration_parameter(iter))
    end

    resampling = iterated_model.resampling

    if !(resampling isa IterationResamplingTypes)
        message *= WARN_POOR_RESAMPLING_CHOICE
    end

    if resampling isa MLJBase.TrainTestPairs && length(resampling) !== 1
        message *= WARN_POOR_CHOICE_OF_PAIRS
    end

    training_control_candidates = filter(iterated_model.controls) do c
        c isa TrainingControl || !(c isa Control)
    end

    if isempty(training_control_candidates)
        throw(ERR_MISSING_TRAINING_CONTROL)
    end

    return message
end


text:
"""
Simple container representing a task in a taskgraph.
"""
struct TaskgraphNode
    "The name of this task."
    name::String
    metadata::Dict{String,Any}

    # Constructor
    TaskgraphNode(name; metadata = emptymeta()) = new(name, metadata)
end

"""
Simple container representing an edge in a taskgraph.
"""
struct TaskgraphEdge
    "Source task names."
    sources::Vector{String}
    "Sink task names."
    sinks::Vector{String}
    metadata::Dict{String,Any}

    function TaskgraphEdge(source, sink; metadata = Dict{String,Any}())
        sources = wrap_vector(source)
        sinks = wrap_vector(sink)
        return new(sources, sinks, metadata)
    end
end

"""
$(SIGNATURES)

Return the names of sources of a `TaskgraphEdge`.
"""
getsources(taskgraph_edge::TaskgraphEdge) = taskgraph_edge.sources

"""
$(SIGNATURES)

Return the names of sinks of a `TaskgraphEdge`.
"""
getsinks(taskgraph_edge::TaskgraphEdge) = taskgraph_edge.sinks

"""
Data structure encoding tasks and their relationships.
"""
struct Taskgraph
    "The name of the taskgraph"
    name::String

    "Nodes in the taskgraph. Type: [`Dict{String, TaskgraphNode}`](@ref TaskgraphNode)"
    nodes::Dict{String,TaskgraphNode}

    "Edges in the taskgraph. Type: [`Vector{TaskgraphEdge}`](@ref TaskgraphEdge)"
    edges::Vector{TaskgraphEdge}

    """
    Outgoing adjacency list mapping node names to edge indices.
    Type: `Dict{String, Vector{Int64}}`
    """
    node_edges_out::Dict{String,Vector{Int}}

    """
    Incoming adjacency list mapping node names to edge indices.
    Type: `Dict{String, Vector{Int64}}`
    """
    node_edges_in::Dict{String,Vector{Int}}

    function Taskgraph(name = "noname")
        return new(
            name,
            Dict{String,TaskgraphNode}(),
            TaskgraphEdge[],
            Dict{String,Vector{Int}}(),
            Dict{String,Vector{Int}}(),
        )
    end

    function Taskgraph(name::String, nodes, edges)
        if eltype(nodes) != TaskgraphNode
            typer = TypeError(
                :Taskgraph, "Incorrect Node Element Type", TaskgraphNode, eltype(nodes)
            )

            throw(typer)
        end
        if eltype(edges) != TaskgraphEdge
            typer = TypeError(
                :Taskgraph, "Incorrect Edge Element Type", TaskgraphEdge, eltype(edges)
            )

            throw(typer)
        end
        # First - create the dictionary to store the nodes. Nodes can be
        # accessed via their name.
        nodes = Dict(n.name => n for n in nodes)
        # Initialize the adjacency lists with an entry for each node. Initialize
        # the values to empty arrays of edges so down-stream algortihms won't
        # have to check if an adjacency list exists for a node.
        edges = collect(edges)
        node_edges_out = Dict(name => Int[] for name in keys(nodes))
        node_edges_in = Dict(name => Int[] for name in keys(nodes))
        # Iterate through all edges - grow adjacency lists correctly.
        for (index, edge) in enumerate(edges)
            for source in edge.sources
                push!(node_edges_out[source], index)
            end
            for sink in edge.sinks
                push!(node_edges_in[sink], index)
            end
        end
        # Return the data structure
        return new(name, nodes, edges, node_edges_out, node_edges_in)
    end
end

# Allow construction without any name.
Taskgraph(nodes, edges) = Taskgraph("noname", nodes, edges)

################################################################################
# METHODS FOR THE TASKGRAPH
################################################################################
# -- Some accessor methods.

"""
$(SIGNATURES)

Return an iterator of [`TaskgraphNode`](@ref) yielding all nodes in `taskgraph`.
"""
getnodes(taskgraph::Taskgraph) = values(taskgraph.nodes)

"""
$(SIGNATURES)

Return an iterator of [`TaskgraphEdge`](@ref) yielding all edges in `taskgraph`.
"""
getedges(taskgraph::Taskgraph) = taskgraph.edges

"""
$(SIGNATURES)

Return the [`TaskgraphNode`](@ref) in `taskgraph` with `name`.
"""
getnode(taskgraph::Taskgraph, name::String) = taskgraph.nodes[name]

"""
$(SIGNATURES)

Return the [`TaskgraphEdge`](@ref) in `taskgraph` with `index`.
"""
getedge(taskgraph::Taskgraph, index::Integer) = taskgraph.edges[index]

# -- helpful query methods.
"""
$(SIGNATURES)

Return an iterator yielding all names of nodes in `taskgraph`.
"""
nodenames(taskgraph::Taskgraph) = keys(taskgraph.nodes)

"""
$(SIGNATURES)

Return the number of nodes in `taskgraph`.
"""
num_nodes(taskgraph::Taskgraph) = length(taskgraph.nodes)

"""
$(SIGNATURES)

Return the number of edges in `taskgraph`.
"""
num_edges(taskgraph::Taskgraph) = length(taskgraph.edges)

"""
$(SIGNATURES)

Return `Vector{TaskgraphNode}` of sources for `edge`.
"""
function getsources(taskgraph::Taskgraph, edge::TaskgraphEdge)
    return (taskgraph.nodes[n] for n in getsources(edge))
end

"""
$(SIGNATURES)

Return `Vector{TaskgraphNode}` of sinks for `edge`.
"""
function getsinks(taskgraph::Taskgraph, edge::TaskgraphEdge)
    return (taskgraph.nodes[n] for n in getsinks(edge))
end

"""
$(SIGNATURES)

Add a `node` to `taskgraph`. Error if node already exists.
"""
function add_node(taskgraph::Taskgraph, node::TaskgraphNode)
    if haskey(taskgraph.nodes, node.name)
        error("Task $(node.name) already exists in taskgraph.")
    end
    taskgraph.nodes[node.name] = node
    # Create adjacency list entries for the new nodes
    taskgraph.node_edges_out[node.name] = TaskgraphEdge[]
    taskgraph.node_edges_in[node.name] = TaskgraphEdge[]
    return nothing
end

"""
$(SIGNATURES)

Add a `edge` to `taskgraph`.
"""
function add_edge(taskgraph::Taskgraph, edge::TaskgraphEdge)
    # Update the edge array
    push!(taskgraph.edges, edge)
    index = length(taskgraph.edges)
    # Update the adjacency lists.
    for source in edge.sources
        push!(taskgraph.node_edges_out[source], index)
    end
    for sink in edge.sinks
        push!(taskgraph.node_edges_in[sink], index)
    end
    return nothing
end

# Methods for accessing the adjacency lists
out_edges(t::Taskgraph, task::String) = [t.edges[i] for i in t.node_edges_out[task]]
out_edges(t::Taskgraph, task::TaskgraphNode) = out_edges(t, task.name)

out_edge_indices(t::Taskgraph, task::String) = t.node_edges_out[task]
out_edge_indices(t::Taskgraph, task::TaskgraphNode) = out_edge_indices(t, task.name)

in_edges(t::Taskgraph, task::String) = [t.edges[i] for i in t.node_edges_in[task]]
in_edges(t::Taskgraph, task::TaskgraphNode) = in_edges(t, task.name)

in_edge_indices(t::Taskgraph, task::String) = t.node_edges_in[task]
in_edge_indices(t::Taskgraph, task::TaskgraphNode) = in_edge_indices(t, task.name)

"""
$(SIGNATURES)

Return `true` if `taskgraph` has a task named `node`.
"""
hasnode(taskgraph::Taskgraph, node::String) = haskey(taskgraph.nodes, node)

"""
$(SIGNATURES)

Return `Set{String}` of names of unique nodes that are the sink of an edges
starting at `node`.
"""
function outnode_names(taskgraph::Taskgraph, node)
    edges = out_edges(taskgraph, node)
    node_names = Set{String}()
    for edge in edges
        union!(node_names, getsinks(edge))
    end
    return node_names
end

"""
$(SIGNATURES)

Return [`Vector{TaskgraphNode}`](@ref TaskgraphNode) of unique nodes that are
the sink of an edge starting at `node`.
"""
function outnodes(taskgraph::Taskgraph, node)
    names = outnode_names(taskgraph, node)
    return [getnode(taskgraph, name) for name in names]
end

"""
$(SIGNATURES)

Return `Set{String}` of names of unique nodes that are the source of an edges
ending at `node`.
"""
function innode_names(taskgraph::Taskgraph, node)
    edges = in_edges(taskgraph, node)
    node_names = Set{String}()
    for edge in edges
        union!(node_names, getsources(edge))
    end
    return node_names
end

"""
$(SIGNATURES)

Return [`Vector{TaskgraphNode}`](@ref TaskgraphNode) of unique nodes that are
the source of an edge ending at `node`.
"""
function innodes(taskgraph::Taskgraph, node)
    names = innode_names(taskgraph, node)
    return [getnode(taskgraph, name) for name in names]
end


text:
using CartesianGP
using Base.Test

reload("Utilities.jl")

@test output_mask(1) == convert(BitString,0b11)
@test output_mask(2) == convert(BitString,0b1111)
@test output_mask(3) == convert(BitString,0b11111111)
@test output_mask(4) == convert(BitString,0b1111111111111111)

@test std_input_context(1) == BitString[0b10]
@test std_input_context(2) == BitString[0b1100,0b1010]
@test std_input_context(3) == BitString[0b11110000,0b11001100,0b10101010]
@test std_input_context(4) == BitString[0b1111111100000000,0b1111000011110000,0b1100110011001100,0b1010101010101010]


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/InterfaceMechanics.jl/blob/master/LICENSE
using InterfaceMechanics, Test, StaticArrays

parameters = CouloumbParameterState(mu=0.2,
                                    elastic_slip=0.001,
                                    c0 = 0.001,
                                    p0 = 100.0)

h0 = -0.005 # Initial overclosure
dh = 0.006 # Overclosure
du = 1.0
drivers = CouloumbDriverState(displacements = SVector{3,Float64}([h0, 0.0, 0.0]))
dtime = 1.0
ddrivers = CouloumbDriverState(time = dtime,
                               displacements = SVector{3,Float64}([dh, 0.0, 0.0]))

interface = Couloumb(parameters=parameters,
                     drivers=drivers,
                     ddrivers=ddrivers)

@info("********* Starting simulation *********")
# Step 1: Initialize contact
@info "time = $(interface.drivers.time), traction = $(interface.variables.traction)"
integrate_interface!(interface)
update!(interface)
@info "time = $(interface.drivers.time), traction = $(interface.variables.traction)"
@test interface.variables.traction[1]>0
@test interface.variables.traction[2]==0
@test interface.variables.traction[3]==0

# Step 2: Tangential movement 1
ddrivers = CouloumbDriverState(time = dtime,
                               displacements = SVector{3,Float64}([0.0, du, 0.0]))
interface.ddrivers = ddrivers
integrate_interface!(interface)
update!(interface)
@info "time = $(interface.drivers.time), traction = $(interface.variables.traction)"
@test interface.variables.traction[1]>0
@test interface.variables.traction[2]>0
@test interface.variables.traction[3]==0

# Step 3: Tangential movement 2
ddrivers = CouloumbDriverState(time = dtime,
                               displacements = SVector{3,Float64}([0.0, -du, 0.0]))
interface.ddrivers = ddrivers
integrate_interface!(interface)
update!(interface)
@info "time = $(interface.drivers.time), traction = $(interface.variables.traction)"
@test interface.variables.traction[1]>0
@test interface.variables.traction[2]<0
@test interface.variables.traction[3]==0

# Step 3: Tangential movement 3
ddrivers = CouloumbDriverState(time = dtime,
                               displacements = SVector{3,Float64}([0.0, du/sqrt(2), du/sqrt(2)]))
interface.ddrivers = ddrivers
integrate_interface!(interface)
update!(interface)
@info "time = $(interface.drivers.time), traction = $(interface.variables.traction)"
@test interface.variables.traction[1]>0
@test interface.variables.traction[2]>0
@test interface.variables.traction[3]>0


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
    DomesticHotWaterLoad

There are many ways in which a DomesticHotWaterLoad can be defined:
1. When using either `doe_reference_name` or `blended_doe_reference_names` in an `ElectricLoad` one only needs to provide the input key "DomesticHotWaterLoad" in the `Scenario` (JSON or Dict). In this case the values from DoE reference names from the `ElectricLoad` will be used to define the `DomesticHotWaterLoad`.
2. One can provide the `doe_reference_name` or `blended_doe_reference_names` directly in the `DomesticHotWaterLoad` key within the `Scenario`. These values can be combined with the `annual_mmbtu` or `monthly_mmbtu` inputs to scale the DoE reference profile(s).
3. One can provide the `fuel_loads_mmbtu_per_hour` value in the `DomesticHotWaterLoad` key within the `Scenario`.

```julia
function DomesticHotWaterLoad(;
    doe_reference_name::String = "",
    city::String = "",
    blended_doe_reference_names::Array{String, 1} = String[],
    blended_doe_reference_percents::Array{<:Real,1} = Real[],
    annual_mmbtu::Union{Real, Nothing} = nothing,
    monthly_mmbtu::Array{<:Real,1} = Real[],
    fuel_loads_mmbtu_per_hour::Array{<:Real,1} = Real[],
    time_steps_per_hour::Int = 1,
    latitude::Real=0.0,
    longitude::Real=0.0
)
```
"""
struct DomesticHotWaterLoad
    loads_kw::Array{Real, 1}

    function DomesticHotWaterLoad(;
        doe_reference_name::String = "",
        city::String = "",
        blended_doe_reference_names::Array{String, 1} = String[],
        blended_doe_reference_percents::Array{<:Real,1} = Real[],
        annual_mmbtu::Union{Real, Nothing} = nothing,
        monthly_mmbtu::Array{<:Real,1} = Real[],
        # addressable_load_fraction,  # TODO
        fuel_loads_mmbtu_per_hour::Array{<:Real,1} = Real[],  # TODO shouldn't this be per_time_step?
        time_steps_per_hour::Int = 1,
        latitude::Real=0.0,
        longitude::Real=0.0
    )
        if length(fuel_loads_mmbtu_per_hour) > 0

            if !(length(fuel_loads_mmbtu_per_hour) / time_steps_per_hour ≈ 8760)
                @error "Provided domestic hot water load does not match the time_steps_per_hour."
            end

            loads_kw = fuel_loads_mmbtu_per_hour .* (MMBTU_TO_KWH * EXISTING_BOILER_EFFICIENCY)

        elseif !isempty(doe_reference_name)
            loads_kw = BuiltInDomesticHotWaterLoad(city, doe_reference_name, latitude, longitude, 2017, annual_mmbtu, monthly_mmbtu)

        elseif length(blended_doe_reference_names) > 1 &&
            length(blended_doe_reference_names) == length(blended_doe_reference_percents)
            loads_kw = blend_and_scale_doe_profiles(BuiltInDomesticHotWaterLoad, latitude, longitude, 2017,
                                                    blended_doe_reference_names, blended_doe_reference_percents, city,
                                                    annual_mmbtu, monthly_mmbtu)
        else
            error("Cannot construct DomesticHotWaterLoad. You must provide either [fuel_loads_mmbtu_per_hour],
                [doe_reference_name, city],
                or [blended_doe_reference_names, blended_doe_reference_percents, city].")
        end

        if length(loads_kw) < 8760*time_steps_per_hour
            loads_kw = repeat(loads_kw, inner=Int(time_steps_per_hour / (length(loads_kw)/8760)))
            @info "Repeating domestic hot water loads in each hour to match the time_steps_per_hour."
        end

        new(
            loads_kw
        )
    end
end

"""
    SpaceHeatingLoad

There are many ways in which a SpaceHeatingLoad can be defined:
1. When using either `doe_reference_name` or `blended_doe_reference_names` in an `ElectricLoad` one only needs to provide the input key "SpaceHeatingLoad" in the `Scenario` (JSON or Dict). In this case the values from DoE reference names from the `ElectricLoad` will be used to define the `SpaceHeatingLoad`.
2. One can provide the `doe_reference_name` or `blended_doe_reference_names` directly in the `SpaceHeatingLoad` key within the `Scenario`. These values can be combined with the `annual_mmbtu` or `monthly_mmbtu` inputs to scale the DoE reference profile(s).
3. One can provide the `fuel_loads_mmbtu_per_hour` value in the `SpaceHeatingLoad` key within the `Scenario`.

```julia
function SpaceHeatingLoad(;
    doe_reference_name::String = "",
    city::String = "",
    blended_doe_reference_names::Array{String, 1} = String[],
    blended_doe_reference_percents::Array{<:Real,1} = Real[],
    annual_mmbtu::Union{Real, Nothing} = nothing,
    monthly_mmbtu::Array{<:Real,1} = Real[],
    fuel_loads_mmbtu_per_hour::Array{<:Real,1} = Real[],
    time_steps_per_hour::Int = 1,
    latitude::Real=0.0,
    longitude::Real=0.0
)
```
"""
struct SpaceHeatingLoad
    loads_kw::Array{Real, 1}

    function SpaceHeatingLoad(;
        doe_reference_name::String = "",
        city::String = "",
        blended_doe_reference_names::Array{String, 1} = String[],
        blended_doe_reference_percents::Array{<:Real,1} = Real[],
        annual_mmbtu::Union{Real, Nothing} = nothing,
        monthly_mmbtu::Array{<:Real,1} = Real[],
        # addressable_load_fraction,  # TODO
        fuel_loads_mmbtu_per_hour::Array{<:Real,1} = Real[],
        time_steps_per_hour::Int = 1,
        latitude::Real=0.0,
        longitude::Real=0.0
    )
        if length(fuel_loads_mmbtu_per_hour) > 0

            if !(length(fuel_loads_mmbtu_per_hour) / time_steps_per_hour ≈ 8760)
                @error "Provided space heating load does not match the time_steps_per_hour."
            end

            loads_kw = fuel_loads_mmbtu_per_hour .* (MMBTU_TO_KWH * EXISTING_BOILER_EFFICIENCY)

        elseif !isempty(doe_reference_name)
            loads_kw = BuiltInSpaceHeatingLoad(city, doe_reference_name, latitude, longitude, 2017, annual_mmbtu, monthly_mmbtu)

        elseif length(blended_doe_reference_names) > 1 &&
            length(blended_doe_reference_names) == length(blended_doe_reference_percents)
            loads_kw = blend_and_scale_doe_profiles(BuiltInSpaceHeatingLoad, latitude, longitude, 2017,
                                                    blended_doe_reference_names, blended_doe_reference_percents, city,
                                                    annual_mmbtu, monthly_mmbtu)
        else
            error("Cannot construct BuiltInSpaceHeatingLoad. You must provide either [fuel_loads_mmbtu_per_hour],
                [doe_reference_name, city],
                or [blended_doe_reference_names, blended_doe_reference_percents, city].")
        end

        if length(loads_kw) < 8760*time_steps_per_hour
            loads_kw = repeat(loads_kw, inner=Int(time_steps_per_hour / (length(loads_kw)/8760)))
            @info "Repeating space heating loads in each hour to match the time_steps_per_hour."
        end

        new(
            loads_kw
        )
    end
end

function BuiltInDomesticHotWaterLoad(
    city::String,
    buildingtype::String,
    latitude::Real,
    longitude::Real,
    year::Int,
    annual_mmbtu::Union{<:Real, Nothing}=nothing,
    monthly_mmbtu::Union{Vector{<:Real}, Nothing}=nothing
    )
    dhw_annual_mmbtu = Dict(
        "Miami" => Dict(
            "FastFoodRest" => 53.47209411,
            "FullServiceRest" => 158.0518043,
            "Hospital" => 442.7295435,
            "LargeHotel" => 3713.248373,
            "LargeOffice" => 127.9412792,
            "MediumOffice" => 22.09603477,
            "MidriseApartment" => 158.0580017,
            "Outpatient" => 27.60091429,
            "PrimarySchool" => 105.3179165,
            "RetailStore" => 0.0,
            "SecondarySchool" => 250.1299246,
            "SmallHotel" => 242.232695,
            "SmallOffice" => 9.891779415,
            "StripMall" => 0.0,
            "Supermarket" => 17.94985187,
            "warehouse" => 0.0,
            "FlatLoad" => 333.0450133
        ),
        "Houston" => Dict(
            "FastFoodRest" => 62.56835989,
            "FullServiceRest" => 188.292814,
            "Hospital" => 530.6352726,
            "LargeHotel" => 4685.666667,
            "LargeOffice" => 160.8917808,
            "MediumOffice" => 25.9266894,
            "MidriseApartment" => 199.2544784,
            "Outpatient" => 32.87691943,
            "PrimarySchool" => 128.0705362,
            "RetailStore" => 0.0,
            "SecondarySchool" => 314.7654465,
            "SmallHotel" => 290.6293673,
            "SmallOffice" => 10.27839347,
            "StripMall" => 0.0,
            "Supermarket" => 19.86608717,
            "warehouse" => 0.0,
            "FlatLoad" => 415.6076756
        ),
        "Phoenix" => Dict(
            "FastFoodRest" => 57.34025418,
            "FullServiceRest" => 170.9086319,
            "Hospital" => 480.098265,
            "LargeHotel" => 4127.191046,
            "LargeOffice" => 141.8507451,
            "MediumOffice" => 23.71397275,
            "MidriseApartment" => 175.5949563,
            "Outpatient" => 29.83372212,
            "PrimarySchool" => 116.811664,
            "RetailStore" => 0.0,
            "SecondarySchool" => 285.2339344,
            "SmallHotel" => 262.8487714,
            "SmallOffice" => 10.05557471,
            "StripMall" => 0.0,
            "Supermarket" => 18.7637822,
            "warehouse" => 0.0,
            "FlatLoad" => 368.7653325
        ),
        "Atlanta" => Dict(
            "FastFoodRest" => 71.33170579,
            "FullServiceRest" => 217.4332205,
            "Hospital" => 615.3498557,
            "LargeHotel" => 5622.340656,
            "LargeOffice" => 192.7164525,
            "MediumOffice" => 29.62182675,
            "MidriseApartment" => 238.9315749,
            "Outpatient" => 37.9759973,
            "PrimarySchool" => 148.8362119,
            "RetailStore" => 0.0,
            "SecondarySchool" => 372.2083434,
            "SmallHotel" => 337.2688069,
            "SmallOffice" => 10.65138846,
            "StripMall" => 0.0,
            "Supermarket" => 21.71038069,
            "warehouse" => 0.0,
            "FlatLoad" => 494.7735263
        ),
        "LasVegas" => Dict(
            "FastFoodRest" => 63.63848459,
            "FullServiceRest" => 191.8494897,
            "Hospital" => 540.9697668,
            "LargeHotel" => 4800.331564,
            "LargeOffice" => 164.7154124,
            "MediumOffice" => 26.36796732,
            "MidriseApartment" => 204.1120165,
            "Outpatient" => 33.48190098,
            "PrimarySchool" => 131.9651451,
            "RetailStore" => 0.0,
            "SecondarySchool" => 327.441087,
            "SmallHotel" => 296.3578765,
            "SmallOffice" => 10.32392915,
            "StripMall" => 0.0,
            "Supermarket" => 20.08676069,
            "warehouse" => 0.0,
            "FlatLoad" => 425.7275876
        ),
        "LosAngeles" => Dict(
            "FastFoodRest" => 69.63212501,
            "FullServiceRest" => 211.7827529,
            "Hospital" => 598.9350422,
            "LargeHotel" => 5440.174033,
            "LargeOffice" => 186.6199083,
            "MediumOffice" => 28.91483286,
            "MidriseApartment" => 231.215325,
            "Outpatient" => 37.00823296,
            "PrimarySchool" => 142.8059487,
            "RetailStore" => 0.0,
            "SecondarySchool" => 352.7467563,
            "SmallHotel" => 328.1935523,
            "SmallOffice" => 10.58011717,
            "StripMall" => 0.0,
            "Supermarket" => 21.35337379,
            "warehouse" => 0.0,
            "FlatLoad" => 478.7476251
        ),
        "SanFrancisco" => Dict(
            "FastFoodRest" => 77.13092952,
            "FullServiceRest" => 236.7180594,
            "Hospital" => 671.40531,
            "LargeHotel" => 6241.842643,
            "LargeOffice" => 213.8445094,
            "MediumOffice" => 32.07909301,
            "MidriseApartment" => 265.1697301,
            "Outpatient" => 41.35500136,
            "PrimarySchool" => 160.4507431,
            "RetailStore" => 0.0,
            "SecondarySchool" => 401.395655,
            "SmallHotel" => 368.0979112,
            "SmallOffice" => 10.90004379,
            "StripMall" => 0.0,
            "Supermarket" => 22.9292287,
            "warehouse" => 0.0,
            "FlatLoad" => 546.4574286
        ),
        "Baltimore" => Dict(
            "FastFoodRest" => 78.2191761,
            "FullServiceRest" => 240.338156,
            "Hospital" => 681.9322322,
            "LargeHotel" => 6358.710286,
            "LargeOffice" => 217.7306132,
            "MediumOffice" => 32.52815422,
            "MidriseApartment" => 270.1195541,
            "Outpatient" => 41.96148216,
            "PrimarySchool" => 165.3116185,
            "RetailStore" => 0.0,
            "SecondarySchool" => 417.9512972,
            "SmallHotel" => 373.906416,
            "SmallOffice" => 10.94554028,
            "StripMall" => 0.0,
            "Supermarket" => 23.15795696,
            "warehouse" => 0.0,
            "FlatLoad" => 557.0507802
        ),
        "Albuquerque" => Dict(
            "FastFoodRest" => 76.9149868,
            "FullServiceRest" => 235.9992545,
            "Hospital" => 669.3128607,
            "LargeHotel" => 6219.08303,
            "LargeOffice" => 212.9944774,
            "MediumOffice" => 31.97726287,
            "MidriseApartment" => 264.2063457,
            "Outpatient" => 41.20639013,
            "PrimarySchool" => 162.1556119,
            "RetailStore" => 0.0,
            "SecondarySchool" => 409.1649863,
            "SmallHotel" => 366.9712928,
            "SmallOffice" => 10.88949351,
            "StripMall" => 0.0,
            "Supermarket" => 22.88525618,
            "warehouse" => 0.0,
            "FlatLoad" => 545.235078
        ),
        "Seattle" => Dict(
            "FastFoodRest" => 81.80231236,
            "FullServiceRest" => 252.2609525,
            "Hospital" => 716.6111323,
            "LargeHotel" => 6741.736717,
            "LargeOffice" => 230.8057849,
            "MediumOffice" => 34.04746055,
            "MidriseApartment" => 286.3412104,
            "Outpatient" => 44.07342164,
            "PrimarySchool" => 172.0233322,
            "RetailStore" => 0.0,
            "SecondarySchool" => 434.0806311,
            "SmallHotel" => 392.968915,
            "SmallOffice" => 11.09863592,
            "StripMall" => 0.0,
            "Supermarket" => 23.91178737,
            "warehouse" => 0.0,
            "FlatLoad" => 588.8601433
        ),
        "Chicago" => Dict(
            "FastFoodRest" => 84.2645196,
            "FullServiceRest" => 260.4454844,
            "Hospital" => 740.4172516,
            "LargeHotel" => 7005.083356,
            "LargeOffice" => 239.7065959,
            "MediumOffice" => 35.08184587,
            "MidriseApartment" => 297.4938584,
            "Outpatient" => 45.49600079,
            "PrimarySchool" => 179.4639347,
            "RetailStore" => 0.0,
            "SecondarySchool" => 456.8817409,
            "SmallHotel" => 406.0751832,
            "SmallOffice" => 11.2033023,
            "StripMall" => 0.0,
            "Supermarket" => 24.4292392,
            "warehouse" => 0.0,
            "FlatLoad" => 611.6276445
        ),
        "Boulder" => Dict(
            "FastFoodRest" => 83.95201542,
            "FullServiceRest" => 259.3997752,
            "Hospital" => 737.372005,
            "LargeHotel" => 6971.32924,
            "LargeOffice" => 238.572519,
            "MediumOffice" => 34.9486709,
            "MidriseApartment" => 296.06471,
            "Outpatient" => 45.31437164,
            "PrimarySchool" => 178.3378526,
            "RetailStore" => 0.0,
            "SecondarySchool" => 453.228537,
            "SmallHotel" => 404.4154946,
            "SmallOffice" => 11.18970855,
            "StripMall" => 0.0,
            "Supermarket" => 24.36505320,
            "warehouse" => 0.0,
            "FlatLoad" => 608.6556221
        ),
        "Minneapolis" => Dict(
            "FastFoodRest" => 89.48929949,
            "FullServiceRest" => 277.8184269,
            "Hospital" => 790.9262388,
            "LargeHotel" => 7563.607619,
            "LargeOffice" => 258.6874644,
            "MediumOffice" => 37.28641454,
            "MidriseApartment" => 321.1473562,
            "Outpatient" => 48.51884975,
            "PrimarySchool" => 191.9480118,
            "RetailStore" => 0.0,
            "SecondarySchool" => 491.5554097,
            "SmallHotel" => 433.8738637,
            "SmallOffice" => 11.42620649,
            "StripMall" => 0.0,
            "Supermarket" => 25.53218144,
            "warehouse" => 0.0,
            "FlatLoad" => 658.8635839
        ),
        "Helena" => Dict(
            "FastFoodRest" => 90.44011877,
            "FullServiceRest" => 280.9757902,
            "Hospital" => 800.0940058,
            "LargeHotel" => 7665.023574,
            "LargeOffice" => 262.1461576,
            "MediumOffice" => 37.68905029,
            "MidriseApartment" => 325.4421541,
            "Outpatient" => 49.09222188,
            "PrimarySchool" => 193.4573283,
            "RetailStore" => 0.0,
            "SecondarySchool" => 494.7393735,
            "SmallHotel" => 438.9398731,
            "SmallOffice" => 11.46564268,
            "StripMall" => 0.0,
            "Supermarket" => 25.72866824,
            "warehouse" => 0.0,
            "FlatLoad" => 667.2021224
        ),
        "Duluth" => Dict(
            "FastFoodRest" => 98.10641517,
            "FullServiceRest" => 306.4772907,
            "Hospital" => 874.2611723,
            "LargeHotel" => 8484.906093,
            "LargeOffice" => 290.0193773,
            "MediumOffice" => 40.92475821,
            "MidriseApartment" => 360.161261,
            "Outpatient" => 53.53681127,
            "PrimarySchool" => 211.2386551,
            "RetailStore" => 0.0,
            "SecondarySchool" => 543.3733772,
            "SmallHotel" => 479.7414481,
            "SmallOffice" => 11.79316054,
            "StripMall" => 0.0,
            "Supermarket" => 27.3451629,
            "warehouse" => 0.0,
            "FlatLoad" => 736.3678114
        ),
        "Fairbanks" => Dict(
            "FastFoodRest" => 108.5335945,
            "FullServiceRest" => 341.1572799,
            "Hospital" => 975.1062178,
            "LargeHotel" => 9600.267161,
            "LargeOffice" => 327.8820873,
            "MediumOffice" => 45.32138512,
            "MidriseApartment" => 407.3910855,
            "Outpatient" => 59.6203514,
            "PrimarySchool" => 234.2595741,
            "RetailStore" => 0.0,
            "SecondarySchool" => 604.6838786,
            "SmallHotel" => 535.2525234,
            "SmallOffice" => 12.23744003,
            "StripMall" => 0.0,
            "Supermarket" => 29.53958045,
            "warehouse" => 0.0,
            "FlatLoad" => 830.07826
        )
    )
    if isempty(city)
        city = find_ashrae_zone_city(latitude, longitude)
    end
    if !(buildingtype in default_buildings)
        error("buildingtype $(buildingtype) not in $(default_buildings).")
    end
    if isnothing(annual_mmbtu)
        annual_mmbtu = dhw_annual_mmbtu[city][buildingtype]
    end
    built_in_load("domestic_hot_water", city, buildingtype, year, annual_mmbtu, monthly_mmbtu)
end

function BuiltInSpaceHeatingLoad(
    city::String,
    buildingtype::String,
    latitude::Real,
    longitude::Real,
    year::Int,
    annual_mmbtu::Union{<:Real, Nothing}=nothing,
    monthly_mmbtu::Union{Vector{<:Real}, Nothing}=nothing,
    )
    spaceheating_annual_mmbtu = Dict(
        "Miami" => Dict(
            "FastFoodRest" => 5.426780867,
            "FullServiceRest" => 12.03181471,
            "Hospital" => 6248.413294,
            "LargeHotel" => 198.0691407,
            "LargeOffice" => 168.9731637,
            "MediumOffice" => 0.036985655,
            "MidriseApartment" => 38.70606161,
            "Outpatient" => 2559.185872,
            "PrimarySchool" => 49.78021153,
            "RetailStore" => 12.12015432,
            "SecondarySchool" => 203.5185485,
            "SmallHotel" => 9.098564901,
            "SmallOffice" => 0.312524873,
            "StripMall" => 20.73216748,
            "Supermarket" => 101.2785324,
            "warehouse" => 56.0796017,
            "FlatLoad" => 605.2352137
        ),
        "Houston" => Dict(
            "FastFoodRest" => 85.49111065,
            "FullServiceRest" => 199.7942842,
            "Hospital" => 8732.10385,
            "LargeHotel" => 1307.035548,
            "LargeOffice" => 2229.971744,
            "MediumOffice" => 16.25994314,
            "MidriseApartment" => 386.0269973,
            "Outpatient" => 2829.324307,
            "PrimarySchool" => 469.2532935,
            "RetailStore" => 289.0470815,
            "SecondarySchool" => 2011.1678969999998,
            "SmallHotel" => 108.9825885,
            "SmallOffice" => 19.55157672,
            "StripMall" => 292.23235389999996,
            "Supermarket" => 984.7374347000001,
            "warehouse" => 475.9377273,
            "FlatLoad" => 1277.307359
        ),
        "Phoenix" => Dict(
            "FastFoodRest" => 57.89972381,
            "FullServiceRest" => 147.2569493,
            "Hospital" => 9382.021026,
            "LargeHotel" => 896.790817,
            "LargeOffice" => 1584.061452,
            "MediumOffice" => 1.922551528,
            "MidriseApartment" => 290.9887152,
            "Outpatient" => 3076.340876,
            "PrimarySchool" => 305.573525,
            "RetailStore" => 208.66330580000002,
            "SecondarySchool" => 1400.638544,
            "SmallHotel" => 83.98084516,
            "SmallOffice" => 9.988210938,
            "StripMall" => 230.16060699999997,
            "Supermarket" => 972.3008295,
            "warehouse" => 362.42249280000004,
            "FlatLoad" => 1188.188154
        ),
        "Atlanta" => Dict(
            "FastFoodRest" => 168.8402371,
            "FullServiceRest" => 379.5865464,
            "Hospital" => 10467.659959999999,
            "LargeHotel" => 2427.6589879999997,
            "LargeOffice" => 3624.593975,
            "MediumOffice" => 49.00635733,
            "MidriseApartment" => 718.9316697,
            "Outpatient" => 3186.250588,
            "PrimarySchool" => 931.7212450999999,
            "RetailStore" => 627.2489826000001,
            "SecondarySchool" => 3968.4936420000004,
            "SmallHotel" => 202.26124219999997,
            "SmallOffice" => 42.74797302,
            "StripMall" => 615.2240506,
            "Supermarket" => 1880.5304489999999,
            "warehouse" => 930.9449202,
            "FlatLoad" => 1888.856302
        ),
        "LasVegas" => Dict(
            "FastFoodRest" => 100.0877773,
            "FullServiceRest" => 247.21791319999997,
            "Hospital" => 9100.302056,
            "LargeHotel" => 1500.581408,
            "LargeOffice" => 2479.152321,
            "MediumOffice" => 5.220181581,
            "MidriseApartment" => 487.43122850000003,
            "Outpatient" => 2924.8220460000002,
            "PrimarySchool" => 499.6562223,
            "RetailStore" => 386.0185744,
            "SecondarySchool" => 2277.7410649999997,
            "SmallHotel" => 138.4427074,
            "SmallOffice" => 19.16330622,
            "StripMall" => 389.30494280000005,
            "Supermarket" => 1479.302604,
            "warehouse" => 579.7671637999999,
            "FlatLoad" => 1413.3882199999998
        ),
        "LosAngeles" => Dict(
            "FastFoodRest" => 40.90390152,
            "FullServiceRest" => 97.94277036,
            "Hospital" => 10346.1713,
            "LargeHotel" => 707.848762,
            "LargeOffice" => 1458.148818,
            "MediumOffice" => 0.12342009699999999,
            "MidriseApartment" => 265.2851759,
            "Outpatient" => 3417.120585,
            "PrimarySchool" => 318.73600980000003,
            "RetailStore" => 175.104083,
            "SecondarySchool" => 1198.276619,
            "SmallHotel" => 72.42852638,
            "SmallOffice" => 5.898878347,
            "StripMall" => 193.18730269999998,
            "Supermarket" => 1040.273464,
            "warehouse" => 323.96697819999997,
            "FlatLoad" => 1228.8385369999999
        ),
        "SanFrancisco" => Dict(
            "FastFoodRest" => 127.22328700000001,
            "FullServiceRest" => 362.48645889999995,
            "Hospital" => 11570.9155,
            "LargeHotel" => 1713.3629670000003,
            "LargeOffice" => 2690.1191,
            "MediumOffice" => 3.8159670660000002,
            "MidriseApartment" => 648.4472797999999,
            "Outpatient" => 3299.0539519999998,
            "PrimarySchool" => 818.2159102,
            "RetailStore" => 569.8081034,
            "SecondarySchool" => 3414.148347,
            "SmallHotel" => 189.0244446,
            "SmallOffice" => 27.53039453,
            "StripMall" => 526.2320428,
            "Supermarket" => 2301.616069,
            "warehouse" => 675.6758453,
            "FlatLoad" => 1808.604729
        ),
        "Baltimore" => Dict(
            "FastFoodRest" => 305.2671204,
            "FullServiceRest" => 657.1337578,
            "Hospital" => 11253.61694,
            "LargeHotel" => 3731.0254619999996,
            "LargeOffice" => 5109.311943,
            "MediumOffice" => 116.8101842,
            "MidriseApartment" => 1132.964052,
            "Outpatient" => 3285.227941,
            "PrimarySchool" => 1428.239177,
            "RetailStore" => 1068.034778,
            "SecondarySchool" => 6557.634924,
            "SmallHotel" => 346.3683857,
            "SmallOffice" => 63.29818348,
            "StripMall" => 1075.39546,
            "Supermarket" => 2929.182261,
            "warehouse" => 1568.722061,
            "FlatLoad" => 2539.2645399999997
        ),
        "Albuquerque" => Dict(
            "FastFoodRest" => 199.73581399999998,
            "FullServiceRest" => 398.5712205,
            "Hospital" => 8371.240776999999,
            "LargeHotel" => 2750.8382260000003,
            "LargeOffice" => 3562.0023950000004,
            "MediumOffice" => 47.49307973,
            "MidriseApartment" => 805.0965778,
            "Outpatient" => 2971.868562,
            "PrimarySchool" => 981.4176700999999,
            "RetailStore" => 755.4523907,
            "SecondarySchool" => 4338.227865999999,
            "SmallHotel" => 232.2194443,
            "SmallOffice" => 43.25360481,
            "StripMall" => 760.0982018,
            "Supermarket" => 2302.228741,
            "warehouse" => 1151.250885,
            "FlatLoad" => 1854.437216
        ),
        "Seattle" => Dict(
            "FastFoodRest" => 255.5992711,
            "FullServiceRest" => 627.5634984000001,
            "Hospital" => 11935.157290000001,
            "LargeHotel" => 3343.683348,
            "LargeOffice" => 5266.970348,
            "MediumOffice" => 28.97979768,
            "MidriseApartment" => 1117.5465470000001,
            "Outpatient" => 3468.128914,
            "PrimarySchool" => 1263.541878,
            "RetailStore" => 952.2758742000001,
            "SecondarySchool" => 6367.850187,
            "SmallHotel" => 310.8087307,
            "SmallOffice" => 49.34878545,
            "StripMall" => 969.1074739000001,
            "Supermarket" => 3004.1844929999997,
            "warehouse" => 1137.398514,
            "FlatLoad" => 2506.1340600000003
        ),
        "Chicago" => Dict(
            "FastFoodRest" => 441.93439000000006,
            "FullServiceRest" => 888.3312571,
            "Hospital" => 12329.57943,
            "LargeHotel" => 5104.848129,
            "LargeOffice" => 7706.028917000001,
            "MediumOffice" => 216.01411800000002,
            "MidriseApartment" => 1482.040156,
            "Outpatient" => 3506.5381090000005,
            "PrimarySchool" => 2006.0002120000001,
            "RetailStore" => 1472.8704380000001,
            "SecondarySchool" => 8962.172873,
            "SmallHotel" => 479.4653436000001,
            "SmallOffice" => 94.19308949,
            "StripMall" => 1497.556168,
            "Supermarket" => 3696.2112950000005,
            "warehouse" => 2256.477231,
            "FlatLoad" => 3258.766323
        ),
        "Boulder" => Dict(
            "FastFoodRest" => 306.8980525,
            "FullServiceRest" => 642.8843574,
            "Hospital" => 9169.381845,
            "LargeHotel" => 3975.1080020000004,
            "LargeOffice" => 5027.882454,
            "MediumOffice" => 124.26913059999998,
            "MidriseApartment" => 1098.944993,
            "Outpatient" => 3087.969786,
            "PrimarySchool" => 1356.396807,
            "RetailStore" => 1086.9187570000001,
            "SecondarySchool" => 6268.036872,
            "SmallHotel" => 342.77800099999996,
            "SmallOffice" => 65.95714912,
            "StripMall" => 1093.093638,
            "Supermarket" => 2966.790122,
            "warehouse" => 1704.8648210000001,
            "FlatLoad" => 2394.8859239999997
        ),
        "Minneapolis" => Dict(
            "FastFoodRest" => 588.8854722,
            "FullServiceRest" => 1121.229499,
            "Hospital" => 13031.2313,
            "LargeHotel" => 6359.946704,
            "LargeOffice" => 10199.279129999999,
            "MediumOffice" => 394.1525556,
            "MidriseApartment" => 1814.148381,
            "Outpatient" => 3661.1462229999997,
            "PrimarySchool" => 2600.964302,
            "RetailStore" => 1869.8106289999998,
            "SecondarySchool" => 11963.323859999999,
            "SmallHotel" => 618.0427338999999,
            "SmallOffice" => 128.12525349999999,
            "StripMall" => 1952.731917,
            "Supermarket" => 4529.776664,
            "warehouse" => 3231.223746,
            "FlatLoad" => 4004.001148
        ),
        "Helena" => Dict(
            "FastFoodRest" => 468.8276835,
            "FullServiceRest" => 934.8994934,
            "Hospital" => 10760.57411,
            "LargeHotel" => 5554.910785,
            "LargeOffice" => 7373.056709,
            "MediumOffice" => 239.8330306,
            "MidriseApartment" => 1531.102079,
            "Outpatient" => 3390.42972,
            "PrimarySchool" => 2097.777112,
            "RetailStore" => 1494.85988,
            "SecondarySchool" => 9535.484059,
            "SmallHotel" => 499.85992930000003,
            "SmallOffice" => 98.85818175,
            "StripMall" => 1604.0043970000002,
            "Supermarket" => 3948.5338049999996,
            "warehouse" => 2504.784991,
            "FlatLoad" => 3252.362248
        ),
        "Duluth" => Dict(
            "FastFoodRest" => 738.1353594999999,
            "FullServiceRest" => 1400.36692,
            "Hospital" => 14179.84149,
            "LargeHotel" => 7781.9012760000005,
            "LargeOffice" => 12504.64187,
            "MediumOffice" => 468.2112216,
            "MidriseApartment" => 2204.85149,
            "Outpatient" => 3774.3233130000003,
            "PrimarySchool" => 3160.1200719999997,
            "RetailStore" => 2298.8242920000002,
            "SecondarySchool" => 14468.64346,
            "SmallHotel" => 772.5386662000001,
            "SmallOffice" => 155.8350887,
            "StripMall" => 2411.847491,
            "Supermarket" => 5587.977185,
            "warehouse" => 3962.122014,
            "FlatLoad" => 4741.886326
        ),
        "Fairbanks" => Dict(
            "FastFoodRest" => 1245.3608279999999,
            "FullServiceRest" => 2209.293209,
            "Hospital" => 20759.042680000002,
            "LargeHotel" => 12298.7791,
            "LargeOffice" => 23214.51532,
            "MediumOffice" => 949.8812392000001,
            "MidriseApartment" => 3398.039504,
            "Outpatient" => 4824.076322999999,
            "PrimarySchool" => 6341.861225,
            "RetailStore" => 3869.670979,
            "SecondarySchool" => 25619.149269999998,
            "SmallHotel" => 1264.41064,
            "SmallOffice" => 297.08593010000004,
            "StripMall" => 3934.89178,
            "Supermarket" => 8515.422039000001,
            "warehouse" => 6882.6512680000005,
            "FlatLoad" => 7851.508208
        )
    )
    if isempty(city)
        city = find_ashrae_zone_city(latitude, longitude)
    end
    if !(buildingtype in default_buildings)
        error("buildingtype $(buildingtype) not in $(default_buildings).")
    end
    if isnothing(annual_mmbtu)
        annual_mmbtu = spaceheating_annual_mmbtu[city][buildingtype]
    end
    built_in_load("space_heating", city, buildingtype, year, annual_mmbtu, monthly_mmbtu)
end


text:
# Autogenerated wrapper script for Pixman_jll for aarch64-apple-darwin
export libpixman

JLLWrappers.@generate_wrapper_header("Pixman")
JLLWrappers.@declare_library_product(libpixman, "@rpath/libpixman-1.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libpixman,
        "lib/libpixman-1.0.40.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module GLSLBuildins

# not in base, so inferable stump must be defined
gl_FragCoord() = zero(Vec2f0)
# replace rule for transpiler
gl_FragCoord(backend::GL) = TGlobal("gl_FragCoord", Vec2f0)

gl_InstanceID() = Int32(0)
gl_InstanceID(backend::GL) = TGlobal("gl_InstanceID", Vec2f0)

#Base.sin(x) already defined
sin{T}(backend::GL, x::Type{T}) = TExpr("sin", T)

EmitVertex() = nothing
sin{T}(backend::GL) = TExpr("EmitVertex", T)

end


text:
# This file is a part of MGVI.jl, licensed under the MIT License (MIT).

using MGVI

using Distributions
using LinearAlgebra
using Random

Test.@testset "test_normal_mvnormal_jac" begin

    epsilon = 1E-5
    Random.seed!(145)

    true_params = randn(4)

    function _common_params(p)
        μ1 = p[1] - p[2]
        μ2 = p[2] * p[1]
        σ1 = p[3]^2 + 10.
        σ2 = (p[3] + p[4])^2 + 10.

        μ1, σ1, μ2, σ2
    end

    function normal_model(p)
        μ1, σ1, μ2, σ2 = _common_params(p)
        n1 = Normal(μ1, σ1)
        n2 = Normal(μ2, σ2)
        Product([n1, n2])
    end
    _flat_normal_model = MGVI.flat_params ∘ normal_model

    # NOTE: models are not the same. Σ[1, 1] = σ1^2, while here we assign σ1.
    # We do this because in this way it is easier to compare jacobians.
    # While in Normal, parameters are μ and σ, in MvNormal corresponding paramters
    # are μ and σ^2.
    function mvnormal_model(p)
        μ1, σ1, μ2, σ2 = _common_params(p)
        MvNormal([μ1, μ2], [σ1 0.;0. σ2])
    end
    _flat_mvnormal_model = MGVI.flat_params ∘ mvnormal_model

    normal_full_jac = FullJacobianFunc(_flat_normal_model)(true_params)
    normal_fwdder_jac = FwdDerJacobianFunc(_flat_normal_model)(true_params)
    normal_fwdrevad_jac = FwdRevADJacobianFunc(_flat_normal_model)(true_params)

    mvnormal_full_jac = FullJacobianFunc(_flat_mvnormal_model)(true_params)
    mvnormal_fwdder_jac = FwdDerJacobianFunc(_flat_mvnormal_model)(true_params)
    mvnormal_fwdrevad_jac = FwdRevADJacobianFunc(_flat_mvnormal_model)(true_params)

    for i in 1:min(size(normal_full_jac)...)

        vec = rand(size(normal_full_jac, 2))

        # We have to reorder jacobian for MvNormal and pick only relevant parts
        # While stacked Normals have parameters [μ1, σ1, μ2, σ2]
        # MvNormal is parametrized as [μ1, μ2, Σ11, Σ12, Σ22]
        # Moreover. In our example Σ12 = 0. Thus we skip this element
        mv_reorder = [1, 3, 2, 5]

        Test.@test norm(normal_full_jac*vec - (mvnormal_full_jac*vec)[mv_reorder]) < epsilon
        Test.@test norm(normal_fwdder_jac*vec - (mvnormal_fwdder_jac*vec)[mv_reorder]) < epsilon
        Test.@test norm(normal_fwdrevad_jac*vec - (mvnormal_fwdrevad_jac*vec)[mv_reorder]) < epsilon

    end

end

Test.@testset "test_normal_mvnormal_logpdf_der" begin

    epsilon = 1E-5
    Random.seed!(145)

    dim = 2
    data = randn(dim)
    pardim = dim + dim
    true_params = randn(pardim)

    function _common_params(p)
        μ1 = p[1] - p[2]
        μ2 = p[2] * p[1]
        σ1 = p[3]^2 + 10.
        σ2 = (p[3] + p[4])^2 + 10.

        μ1, σ1, μ2, σ2
    end

    function normal_model(p)
        μ1, σ1, μ2, σ2 = _common_params(p)
        n1 = Normal(μ1, σ1)
        n2 = Normal(μ2, σ2)
        Product([n1, n2])
    end

    # est_res_sampler = MGVI._create_residual_sampler(normal_model, true_params;
                                               # residual_sampler=ImplicitResidualSampler,
                                               # jacobian_func=FwdRevADJacobianFunc,
                                               # residual_sampler_options=(;))
    # residual_samples = rand(Random.GLOBAL_RNG, est_res_sampler, 5)

    residual_samples = [
        1.57327   -0.545016   -0.468532  -0.148169  -0.0476087
        0.889662  -0.32206    -0.776208  -1.3703    -0.262721
        0.535246  -0.0152493   0.847152   0.723876  -0.0277677
        2.13597    0.252785   -0.278254   1.11853   -0.189659
    ]

    normal_par_dim = dim + dim
    normal_kl(params::AbstractVector) = MGVI.mgvi_kl(normal_model, data, residual_samples, params)
    normal_grad = zeros(pardim)
    MGVI._gradient_for_optim(normal_kl)(normal_grad, true_params)

    function mvnormal_model(p)
        μ1, σ1, μ2, σ2 = _common_params(p)
        MvNormal([μ1, μ2], [σ1^2 0.;0. σ2^2])
    end
    mvnormal_kl(params::AbstractVector) = MGVI.mgvi_kl(mvnormal_model, data, residual_samples, params)
    mvnormal_grad = zeros(pardim)
    MGVI._gradient_for_optim(mvnormal_kl)(mvnormal_grad, true_params)

    Test.@test norm(normal_grad - mvnormal_grad) < epsilon

end


text:
function printUpdate(o...)
    if isinteractive()
        print("\u1b[1G")
        print(o...)
        print("\u1b[K")
    end
end

function printlnFlush(o...)
    println(o...)
    if !isinteractive()
        flush(STDOUT)
    end
end

function saveNetwork(weights,filename)
    w = map(Array,weights)
	save(pretrainedPath("$(filename).jld"),"w",w,compress=true)
end

function loadNetwork(filename)
    w = load(pretrainedPath(filename),"w")
    w = map(KnetArray,w)
end

function pixelwiseSoftloss(ypred,y)
    -sum(y.*logp(ypred,3)) / (size(y,1)*size(y,2)*size(y,4))
end

function pixelwiseAccuracy(ypred,y)
    ncorrect = sum(y.*(ypred.==maximum(ypred,3)))
    ncount = size(ypred,1)*size(ypred,2)*size(ypred,4)
    return ncorrect,ncount
end

max_stars_repo_name: eeyrdkl/MScThesis

max_issues_repo_name: eeyrdkl/MScThesis

max_forks_repo_name: eeyrdkl/MScThesis


text:
#=
   Provides calc_model_rv(theta, time)
   Computes the velocity of the star due to the perturbations of multiple planets, as the linear superposition of the circular orbit induced by each planet, i.e., neglecting eccentricity and mutual planet-planet interactions
=#

# Calculate Keplerian velocity of star due to one planet (with parameters displaced by offset)
function calc_rv_circ_one_planet{T1,T2}( theta::Array{T1,1}, time::T2; plid::Integer = 1 )
  (P,K,M0) = extract_PKM0(theta,plid=plid)
  #(P,Kc,Ks) = extract_PKcKs(theta,plid=plid)
  n = 2pi/P
  #M = mod2pi(time*n-M0)
  M = time*n-M0
  return K*cos(M)
  #return Kc*cos(M)+Ks*sin(M)
end

# Calculate Keplerian velocity of star due to num_pl planets
function calc_rv_circ_multi_planet{T1,T2}( theta::Array{T1,1}, time::T2)
  zdot = zero(T1)
  for plid in 1:num_planets(theta)
      zdot += calc_rv_circ_one_planet(theta,time,plid=plid)
  end
  return zdot
end

# Assumes model parameters = [ Period_1, K_1, M0_1,   Period_2, K_2, M0_2, ...  C ]
function calc_model_rv{T1,T2}( theta::Array{T1,1}, t::T2; obsid::Integer = 1, tol::Real = 1.e-8)
  Offset = extract_rvoffset(theta, obsid=obsid)
  calc_rv_circ_multi_planet(theta, t) + Offset
end

export calc_model_rv


text:
# Read in file of helper functions
# This file includes: cleantext, countwords
include("../word_count_helpers.jl")

# Load the file names
dataLoc = "../../../data/word_count/";
fnames = dataLoc.*readdir(dataLoc)

# Grab the argument that is passed in
# This is the index into fnames for this process
task_id = parse(Int,ARGS[1])
num_tasks = parse(Int,ARGS[2])

# Check to see if the index is valid (so the program exits cleanly if the wrong indices are passed)
for i in task_id+1:num_tasks:length(fnames)

    # Read in file and clean the text
    f = open(fnames[i])
    text=cleantext(readlines(f))

    # Count number of times each word appears
    wordCounts = countwords(text)

    # Sort and print the top 5 words with their counts
    rankedwords = sort(collect(wordCounts), by=x->x[2], rev=true)
    println(rankedwords[1:5,:])

end

max_stars_repo_name: llsc-supercloud/teaching-examples

max_issues_repo_name: llsc-supercloud/teaching-examples

max_forks_repo_name: llsc-supercloud/teaching-examples


text:
function Base.truncate(uv::TheoreticalFittedUncertainScalar,
    constraint::TruncateStd; n_draws::Int = 10000)
    m = mean(uv.distribution.distribution)
    s = std(uv.distribution.distribution)
    lower_bound = m - s
    upper_bound = m + s

    truncated(uv.distribution, lower_bound, upper_bound)
end


text:
filename = ARGS[1]
input = readlines(filename)
orbits = Dict()
for orbit in input
    planet, moon = split(orbit, ')')
    orbits[String(moon)] = String(planet)
end

function get_paths()
    paths = []
    for planet in keys(orbits)
        i = planet
        path = [i]
        while haskey(orbits, i)
            i = orbits[i]
            push!(path, i)
        end
        if "SAN" in path || "YOU" in path
            push!(paths, path)
        end
    end
    return paths
end

function find_orbital_transfer(paths)
    path_one = pop!(paths)
    path_two = pop!(paths)
    transfers = -2 # offset planet jumping
    meet_at = nothing # common root

    # find common root and distance from planet a
    for planet in path_one
        if planet in path_two
            meet_at = planet
            break
        end
        transfers += 1
    end

    # find distance from planet b
    for planet in path_two
        if planet == meet_at
            break
        end
        transfers += 1
    end

    return transfers
end

println(
    find_orbital_transfer(
        get_paths()
    )
)


text:
using PtFEM

nf_path = joinpath(dirname(@__FILE__), "Ex57.1.nf.dat")
loads_path = joinpath(dirname(@__FILE__), "Ex57.1.loads.dat")

data = Dict(
  # Solid(ndim, nst, nxe, nye, nze, nip, direction=:r, finite_element(nod, nodof))
  :struc_el => Solid(3, 6, 20, 60, 40, 8, Hexahedron(20, 3)),
  :properties => [
    100.0 0.3;
     ],
  :x_coords => [
    0.0000,  0.0250,  0.0500,  0.0750,  0.1000,  0.1250,  0.1500,  0.1750,  0.2000,
    0.2250,  0.2500,  0.2750,  0.3000,  0.3250,  0.3500,  0.3750,  0.4000,  0.4250,
    0.4500,  0.4750,  0.5000
  ],
  :y_coords => [
    0.0000, 0.0500, 0.1000, 0.1500, 0.2000, 0.2500, 0.3000, 0.3500, 0.4000,
    0.4500, 0.5000, 0.5500, 0.6000, 0.6500, 0.7000, 0.7500, 0.8000, 0.8500,
    0.9000, 0.9500, 1.0000, 1.0500, 1.1000, 1.1500, 1.2000, 1.2500, 1.3000,
    1.3500, 1.4000, 1.4500, 1.5000, 1.5500, 1.6000, 1.6500, 1.7000, 1.7500,
    1.8000, 1.8500, 1.9000, 1.9500, 2.0000, 2.0500, 2.1000, 2.1500, 2.2000,
    2.2500, 2.3000, 2.3500, 2.4000, 2.4500, 2.5000, 2.5500, 2.6000, 2.6500,
    2.7000, 2.7500, 2.8000, 2.8500, 2.9000, 2.9500, 3.0000
  ],
  :z_coords => [
    0.0000, -0.0500, -0.1000, -0.1500, -0.2000, -0.2500, -0.3000, -0.3500, -0.4000,
   -0.4500, -0.5000, -0.5500, -0.6000, -0.6500, -0.7000, -0.7500, -0.8000, -0.8500,
   -0.9000, -0.9500, -1.0000, -1.0500, -1.1000, -1.1500, -1.2000, -1.2500, -1.3000,
   -1.3500, -1.4000, -1.4500, -1.5000, -1.5500, -1.6000, -1.6500, -1.7000, -1.7500,
   -1.8000, -1.8500, -1.9000, -1.9500, -2.0000
  ],
  :support => PtFEM.read_nf_file(nf_path),
  :loaded_nodes => PtFEM.read_loads_file(loads_path),
  :cg_tol => 1.0e-5,
  :cg_limit => 2000
)

@time m = p56(data)
println()