text:
struct VX
    id::Int
end

@testset "Basic Operations" begin

    g = OptimizingIR.Graphs.DirectedGraph(VX)

    @test OptimizingIR.Graphs.ne(g) == 0
    @test OptimizingIR.Graphs.nv(g) == 0

    v1 = VX(1)
    v2 = VX(2)

    OptimizingIR.Graphs.add_vertex!(g, v1)
    @test !OptimizingIR.Graphs.has_edge(g, v1, v1)

    @test OptimizingIR.Graphs.ne(g) == 0
    @test OptimizingIR.Graphs.nv(g) == 1

    OptimizingIR.Graphs.add_edge!(g, v1, v1)
    @test OptimizingIR.Graphs.has_edge(g, v1, v1)
    @test OptimizingIR.Graphs.ne(g) == 1
    OptimizingIR.Graphs.remove_edge!(g, v1, v1)
    @test !OptimizingIR.Graphs.has_edge(g, v1, v1)
    @test OptimizingIR.Graphs.ne(g) == 0

    OptimizingIR.Graphs.add_vertex!(g, v2)
    @test OptimizingIR.Graphs.ne(g) == 0
    @test OptimizingIR.Graphs.nv(g) == 2

    OptimizingIR.Graphs.add_edge!(g, v1, v2)
    @test OptimizingIR.Graphs.has_edge(g, v1, v2)
    @test !OptimizingIR.Graphs.has_edge(g, v2, v1)
    @test OptimizingIR.Graphs.ne(g) == 1
    @test OptimizingIR.Graphs.nv(g) == 2

    OptimizingIR.Graphs.add_edge!(g, v2, v2)
    @test OptimizingIR.Graphs.has_edge(g, v2, v2)
    @test !OptimizingIR.Graphs.has_edge(g, v2, v1)
    @test OptimizingIR.Graphs.ne(g) == 2
    @test OptimizingIR.Graphs.nv(g) == 2

    println(g)
end

@testset "CFG Example" begin
    cfg = OptimizingIR.Graphs.DirectedGraph(VX)

    for i in 0:8
        OptimizingIR.Graphs.add_vertex!(cfg, VX(i))
    end

    OptimizingIR.Graphs.add_edge!(cfg, VX(0), VX(1))
    OptimizingIR.Graphs.add_edge!(cfg, VX(1), VX(2))
    OptimizingIR.Graphs.add_edge!(cfg, VX(1), VX(5))
    OptimizingIR.Graphs.add_edge!(cfg, VX(2), VX(3))
    OptimizingIR.Graphs.add_edge!(cfg, VX(3), VX(4))
    OptimizingIR.Graphs.add_edge!(cfg, VX(5), VX(6))
    OptimizingIR.Graphs.add_edge!(cfg, VX(5), VX(8))
    OptimizingIR.Graphs.add_edge!(cfg, VX(6), VX(7))
    OptimizingIR.Graphs.add_edge!(cfg, VX(8), VX(7))
    OptimizingIR.Graphs.add_edge!(cfg, VX(7), VX(3))
    OptimizingIR.Graphs.add_edge!(cfg, VX(3), VX(1))

    @test OptimizingIR.Graphs.ne(cfg) == 11
    @test OptimizingIR.Graphs.nv(cfg) == 9

    # successors
    @test collect(OptimizingIR.Graphs.each_successor(cfg, VX(0))) == [ VX(1) ]
    @test OptimizingIR.Graphs.successors(cfg, VX(1)) == Set([ VX(2), VX(5) ])
    @test OptimizingIR.Graphs.successors(cfg, VX(1)) == Set([ VX(5), VX(2) ])
    @test OptimizingIR.Graphs.successors(cfg, VX(2)) == Set([ VX(3) ])
    @test OptimizingIR.Graphs.successors(cfg, VX(3)) == Set([ VX(1), VX(4) ])
    @test OptimizingIR.Graphs.successors(cfg, VX(4)) == Set()
    @test OptimizingIR.Graphs.successors_count(cfg, VX(4)) == 0
    @test eltype(OptimizingIR.Graphs.successors(cfg, VX(4))) == VX
    @test OptimizingIR.Graphs.successors(cfg, VX(5)) == Set([ VX(6), VX(8) ])
    @test OptimizingIR.Graphs.successors(cfg, VX(6)) == Set([ VX(7) ])
    @test OptimizingIR.Graphs.successors(cfg, VX(8)) == Set([ VX(7) ])
    @test OptimizingIR.Graphs.successors(cfg, VX(7)) == Set([ VX(3) ])
    @test OptimizingIR.Graphs.successors_count(cfg, VX(7)) == 1

    # predecessors
    @test OptimizingIR.Graphs.predecessors(cfg, VX(0)) == Set()
    @test OptimizingIR.Graphs.predecessors_count(cfg, VX(0)) == 0
    @test eltype(OptimizingIR.Graphs.predecessors(cfg, VX(0))) == VX
    @test OptimizingIR.Graphs.predecessors(cfg, VX(1)) == Set([ VX(0), VX(3) ])
    @test OptimizingIR.Graphs.predecessors(cfg, VX(2)) == Set([ VX(1) ])
    @test OptimizingIR.Graphs.predecessors(cfg, VX(3)) == Set([ VX(2), VX(7) ])
    @test OptimizingIR.Graphs.predecessors(cfg, VX(4)) == Set([ VX(3) ])
    @test OptimizingIR.Graphs.predecessors(cfg, VX(5)) == Set([ VX(1) ])
    @test OptimizingIR.Graphs.predecessors(cfg, VX(6)) == Set([ VX(5) ])
    @test OptimizingIR.Graphs.predecessors(cfg, VX(7)) == Set([ VX(6), VX(8) ])
    @test OptimizingIR.Graphs.predecessors_count(cfg, VX(7)) == 2

    # BFS
    @test OptimizingIR.Graphs.reachable_vertices(cfg, VX(4), true) == Set([ VX(4) ])
    @test OptimizingIR.Graphs.reachable_vertices(cfg, VX(3), true) == Set([ VX(1), VX(2), VX(3), VX(4), VX(5), VX(6), VX(7), VX(8) ])
    @test OptimizingIR.Graphs.reachable_vertices(cfg, VX(1), false) == Set([ VX(0), VX(1), VX(2), VX(3), VX(5), VX(6), VX(7), VX(8) ])

    # dominance
    @test OptimizingIR.Graphs.dominator_sets(cfg) == [Set([1]), Set([2, 1]), Set([2, 3, 1]), Set([4, 2, 1]), Set([4, 2, 5, 1]), Set([2, 6, 1]), Set([7, 2, 6, 1]), Set([2, 8, 6, 1]), Set([9, 2, 6, 1])]
    @test OptimizingIR.Graphs.strict_dominator_sets(cfg) == [Set(), Set([1]), Set([2, 1]), Set([2, 1]), Set([4, 2, 1]), Set([2, 1]), Set([2, 6, 1]), Set([2, 6, 1]), Set([2, 6, 1])]
    @test OptimizingIR.Graphs.idominators(cfg) == [nothing, 1, 2, 2, 4, 2, 6, 6, 6]
    @test OptimizingIR.Graphs.dominance_frontiers(cfg) == [Set([]), Set([2]), Set([4]), Set([2]), Set([]), Set([4]), Set([8]), Set([4]), Set([8])]
end


text:
using PyPlot, JLD2, Polynomials, PyCall, GLMakie, Images
include("../utils/Simulation.jl")
include("../utils/plotting.jl")

n = (64, 64, 64)
K = zeros(Float32, n)
K[:,:,1:36] .= 15f0
K[:,:,37:38] .= 1f-3
K[:,:,39:end] .= 5000f0

#K = Float32.(imfilter(K, Kernel.gaussian((2,2,2))))

phi = 0.36f0*ones(Float32, n)

for i = 1:n[1]
    for j = 1:n[2]
        for k = 1:n[3]
            p = Polynomial([-0.0314^2*K[i,j,k],2*0.0314^2*K[i,j,k],-0.0314^2*K[i,j,k],1.527^2])
            phi[i,j,k] = minimum(real(roots(p)[findall(real(roots(p)).== roots(p))]))
        end
    end
end

d = (100f0, 100f0, 24f0)
qinj = (d[1]*n[1]/2, d[2]*n[2]/2, d[3]*n[3]+(341-256)*6f0-500f0)
qrate = 7
time = 80
nt = 10

sat, p = TwoPhase(K, phi, qinj, qrate, d, time, nt; o=(0f0,0f0,(341-256)*6f0))

figure();imshow(K[:,32,:]')
figure();imshow(sat[end][:,32,:]')

#plot_3D(reverse(sat[end],dims=3))
#plot_3D(reverse(K,dims=3))
#plot_3D(reverse(phi,dims=3))

ext: jl

max_stars_repo_name: slimgroup/TwoPhaseSpeedTest

max_issues_repo_name: slimgroup/TwoPhaseSpeedTest

max_forks_repo_name: slimgroup/TwoPhaseSpeedTest


text:
### A Pluto.jl notebook ###
# v0.12.11

using Markdown
using InteractiveUtils

# ╔═╡ c9f9318a-1e89-11eb-2939-f3a5f1289a12
begin
	import Pkg; Pkg.add("ForwardDiff")
	Pkg.add("Plots")
	using ForwardDiff
	using Plots
end

# ╔═╡ 652e2558-1e7a-11eb-2483-cb1203507751
begin
	using LinearAlgebra

	using InstantiateFromURL

end

# ╔═╡ e30d5e5c-1e7b-11eb-0173-5dbf06d875bb
md"""
# Introductory Examples
### Solutions to the Exercise"""

# ╔═╡ 02f30026-1e7c-11eb-1bdf-713cc1950a22
md"""
##### `Exercise 1 - 4 are quite naive and will only help you to get familiarity with the platform` """

# ╔═╡ c0b1a1c0-1e8e-11eb-11e4-9187b5708dad
html"""

"""

# ╔═╡ 0297fd2c-1e7c-11eb-20de-ed44583cc810
md"""
##### Exercise 5
Simulate and plot the correlated time series

xt+1 = αxt + ϵt+1 wherex0 =0 and t= 0,…,n.
The sequence of shocks {ϵt} is assumed to be iid and standard normal.
Set n=200 and α=0.9."""

# ╔═╡ c5a9b326-1e7c-11eb-043e-6b700836e7df
begin
	path=accumulate(1:200;init=0) do old,i
		old*0.9 +randn()
	end
	plot(path,legend=false)
end

# ╔═╡ 79411808-1e7e-11eb-1073-d9d8e2b1a662
html"""

"""

# ╔═╡ 55e7aa80-1e7d-11eb-3447-810cb162aa9b
md"""
##### Exercise 6
Plot three simulated time series, one for each of the cases α=0, α=0.8 and α=0.98.

(The figure will illustrate how time series with the same one-step-ahead conditional volatilities, as these three processes have, can have very different unconditional volatilities)"""

# ╔═╡ 96c557f0-1e7d-11eb-0577-1bef484e817a
let p=plot()
	for α in [0,0.8,0.98]
		pa=accumulate(1:200;init=1) do old,i
			old*α +randn()
		end
		plot!(p,pa,label="α = $α")
	end
	p
end

# ╔═╡ 9438bd3c-1e7e-11eb-038d-c759ea6cf550
html"""

"""

# ╔═╡ b83ea3d6-1e7e-11eb-3c70-f5e2349bd425
md"""
##### Exercise 7
*This exercise is more challenging.*

Take a random walk, starting from x0=1
xt+1 = αxt + σϵt+1 where x0=1 and t= 0,…,tmax

* Furthermore, assume that the xtmax=0 (i.e. at tmax, the value drops to zero, regardless of its current state).
* The sequence of shocks {ϵt} is assumed to be iid and standard normal.

* For a given path {xt} define a first-passage time as Ta=min{t|xt≤a}, where by the assumption of the process Ta≤tmax.
Start with σ=0.2,α=1.0.
calculate the first-passage time, T0, for 100 simulated random walks – to a tmax=200 and plot a histogram
plot the sample mean of T0 from the simulation for α∈{0.8,1.0,1.2}"""

# ╔═╡ ae782114-1e7f-11eb-31bb-9f1b009f6853
let left= plot(), right=plot(), σ=0.2
	bigg=[]
	for j in [0.8,1.0,1.2]
		result_=Vector{}()
		for i ∈ 1:100
			pathing=accumulate(1:200;init=1) do old, i
				j*old + σ*randn()
			end
			push!(result_,findfirst(x->x<0.0,pathing))
		end
		push!(bigg,replace(result_,nothing=>200))
	end
	left=bar(bigg[2],bins=:scott,legend=false)
	right=bar(sum.(bigg)./100,bins=:scott,label=[0.8,1.0,1.2])
	plot(left,right)
end

# ╔═╡ a56d8f26-1e89-11eb-2c21-791307ccc9c2
html"""

"""

# ╔═╡ 6167139e-1e89-11eb-26f4-55cbdbc70de2
md"""
##### Exercise 8
*This exercise is more challenging.*

The root of a univariate function f(⋅) is an x such that f(x)=0.
One solution method to find local roots of smooth functions is called Newton’s method.
Starting with an x0 guess, a function f(⋅) and the first-derivative f′(⋅), the algorithm is to repeat

xn+1=xn−f(xn)f′(xn)

until |xn+1−xn| is below a tolerance

Use a variation of the fixedpointmap code to implement Newton’s method, where the function would accept arguments f, f_prime, x_0, tolerance, maxiter.
Test it with f(x)=(x−1)3 and another function of your choice where you can analytically find the derivative."""

# ╔═╡ 3928f3c2-1e8c-11eb-2244-3f9c7c2857d4
begin
	f(x) = (x-1)^3
	f_d(x) = ForwardDiff.derivative(f, x)

	f(1.7),f_d(1.7)
end

# ╔═╡ 4b431106-1e8b-11eb-00cd-814c8e98e0b7
let
	maxitr=200
	ini=rand(1:10)
	second=0
	for i in 1:maxitr
		second= ini - (f(ini)/f_d(ini))
		if abs(second-ini)< 1e-7
			break
		else
			ini=second
		end
	end
	second
end

# ╔═╡ Cell order:
# ╟─e30d5e5c-1e7b-11eb-0173-5dbf06d875bb
# ╟─02f30026-1e7c-11eb-1bdf-713cc1950a22
# ╟─c0b1a1c0-1e8e-11eb-11e4-9187b5708dad
# ╟─0297fd2c-1e7c-11eb-20de-ed44583cc810
# ╟─c5a9b326-1e7c-11eb-043e-6b700836e7df
# ╟─79411808-1e7e-11eb-1073-d9d8e2b1a662
# ╟─55e7aa80-1e7d-11eb-3447-810cb162aa9b
# ╠═96c557f0-1e7d-11eb-0577-1bef484e817a
# ╟─9438bd3c-1e7e-11eb-038d-c759ea6cf550
# ╟─b83ea3d6-1e7e-11eb-3c70-f5e2349bd425
# ╠═ae782114-1e7f-11eb-31bb-9f1b009f6853
# ╟─a56d8f26-1e89-11eb-2c21-791307ccc9c2
# ╟─6167139e-1e89-11eb-26f4-55cbdbc70de2
# ╠═3928f3c2-1e8c-11eb-2244-3f9c7c2857d4
# ╠═4b431106-1e8b-11eb-00cd-814c8e98e0b7
# ╠═c9f9318a-1e89-11eb-2939-f3a5f1289a12
# ╠═652e2558-1e7a-11eb-2483-cb1203507751

max_stars_repo_path: getting started with julia/Introductory Example.jl

max_stars_repo_name: nehamsoni/Quantitative-Economics

max_issues_repo_name: nehamsoni/Quantitative-Economics

max_forks_repo_name: nehamsoni/Quantitative-Economics


text:
# This file was generated, do not modify it. # hide
using Gnuplot, Random
Random.seed!(123)
goffset = "set offsets graph .05, graph .05, graph .05, graph .05" # hide
@gp rand(30) rand(30) "w p pt 6 ps 2 lc '#D55E00' lw 2 t 'marker'"
@gp :- title = "scatter plot" xlabel ="X" ylabel = "Y"
@gp :- "set key b l box opaque"
save(term="pngcairo size 600,400", output="./code/plt_ex11.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
import Random
using Distributions
using Dates: DateTime, Day
using UnPack: @unpack
import Bijectors

abstract type EpiRt end
## ===========================================================================
# Random Walk

struct RandomWalkDistribution{Tn, Ts} <: ContinuousMultivariateDistribution
    n::Tn
    s::Ts
end

Random.rand(rng::Random.AbstractRNG, d::RandomWalkDistribution{Tn, Ts}) where {Tn, Ts} = begin
    @unpack n, s = d
    step_prior = Normal(zero(Ts), s)
    x = Vector{Ts}(undef, n)

    ######## random walk
    x[1] = rand(step_prior)
    for i in 2:n
        x[i] = x[i-1] + rand(step_prior)
    end
    return x
end

Distributions.logpdf(d::RandomWalkDistribution{Tn, Ts}, x::AbstractVector{T}) where {Tn, Ts, T} =
    logpdf( MvNormal( d.n-1, d.s ), diff(x) )# + logpdf( Normal( zero(Ts), x[1] ) )
Bijectors.bijector(d::RandomWalkDistribution) = Bijectors.Identity{1}()

## ===========================================================================
"""
random walk with variable step interval
"""
struct RWModel{Tno, Tnt, Tns, Tin, Tli, Til} <: EpiRt
    nobs        ::Tno
    ntotal      ::Tnt
    nsteps      ::Tns
    index       ::Tin
    link        ::Tli
    invlink     ::Til
end

function stepindex(n, rwstep)
    index = Vector{Int64}(undef, n)
    nsteps = 1
    for i in eachindex(index)
        index[i] = nsteps
        i % rwstep == 0 && ( nsteps += 1 )
    end
    index
end

RWModel(
    dates     ::AbstractVector;
    endobs    = dates[end],
    rwstep    ::Int = 1,
    link      = exp,
    invlink   = log,
) = begin

	#@assert all( dates[1]:Day(1):dates[end] .== dates )
	ntotal  = length(dates)
	nobs    = findfirst(==(endobs), dates)
    index   = stepindex(nobs-1, rwstep)
    nsteps  = length( unique( index ) )

	RWModel(
        nobs,
		ntotal,
		nsteps,
		index,
		link,
		invlink
	)
end

(epirt::RWModel{Tno, Tnt, Tns, Tin, Tli, Til})(
	θ, ::Val{false}
) where {Tno, Tnt, Tns, Tin, Tli, Til} = begin
    @unpack rw, r0 = θ
    @unpack nobs, link = epirt

    rt         = Vector{typeof(r0)}(undef, nobs)
    rt[1]      = r0
    rt[2:nobs] = link.(rw)
    return rt
end

(epirt::RWModel{Tno, Tnt, Tns, Tin, Tli, Til})(
	θ, ::Val{true}
) where {Tno, Tnt, Tns, Tin, Tli, Til} = begin
    @unpack rw, r0 = θ
    @unpack ntotal, nobs, link = epirt

    rt              = Vector{typeof(r0)}(undef, ntotal)
    rt[1]           = r0
    rt[2:nobs]      = link.(rw)
    rt[nobs+1:end] .= rt[nobs]
    return rt
end

## =============================================================================

struct RWModelDistribution2{Tr0, Tσ, Tbm} <: ContinuousMultivariateDistribution
	r0_prior   ::Tr0
	σ_rw_prior ::Tσ
	rwmodel::Tbm
end

# RWModelDistribution2(
# 	r0_prior   ::Distribution,
# 	σ_rw_prior ::Distribution,
# 	args...; kwargs...
# ) = begin
# 	RWModelDistribution2(
# 		r0_prior, σ_rw_prior,
# 		RWModel(args...; kwargs...)
# 	)
# end

Base.rand(rng::Random.AbstractRNG, epirt::RWModelDistribution2) = begin
	@unpack r0_prior, σ_rw_prior, rwmodel = epirt
    @unpack nsteps, invlink = rwmodel

    ######## draw from univariate prior
    σ_rw = rand(rng, σ_rw_prior )
    r0   = rand(rng, r0_prior )

    ######## draw from random walk prior
    rw_prior = RandomWalkDistribution(nsteps, σ_rw)
    rw = rand(rng, rw_prior)
    return (; σ_rw, r0, rw)
end

Distributions.logpdf(epirt::RWModelDistribution2, θ) = begin
    @unpack rw, r0, σ_rw = θ
    @unpack r0_prior, σ_rw_prior, rwmodel = epirt
	@unpack nsteps, invlink = rwmodel

    rw_prior   = RandomWalkDistribution(nsteps, σ_rw)
    step_prior = Normal(zero(σ_rw), σ_rw) #NOTE same prior as in RandomWalkDistribution
    step = first(rw) - invlink(r0)

    target  = logpdf(σ_rw_prior, σ_rw)
    target += logpdf(r0_prior, r0)
    target += logpdf(step_prior, step)
    target += logpdf(rw_prior, rw)
    return target
end

Bijectors.bijector(d::RWModelDistribution2) = begin
	@unpack r0_prior, σ_rw_prior = d
	bj_r0 = Bijectors.bijector(r0_prior)
	bj_σ  = Bijectors.bijector(σ_rw_prior)
	bj_rw = Bijectors.Identity{1}()
	return Bijectors.Stacked( (bj_r0, bj_σ, bj_rw) )
end

## ===========================================================================
"""
random walk with constant early phase
"""
struct ConstantPhaseRWModel{Tr0, Tr1, Tσ, Tnc, Tno, Tnt, Tns, Tin, Tli, Til} <: EpiRt
    r0_prior    ::Tr0
    r1_prior    ::Tr1
    σ_rw_prior  ::Tσ
    nconst      ::Tnc
    nobs        ::Tno
    ntotal      ::Tnt
    nsteps      ::Tns
    index       ::Tin
    link        ::Tli
    invlink     ::Til
end

ConstantPhaseRWModel(
    dates     ::Vector{DateTime},
    startrw   ::DateTime;
    endobs    ::DateTime = dates[end],
    r0_prior  ::ContinuousUnivariateDistribution = truncated(Normal(3.5, 1.), 0, Inf),
    r1_prior  ::ContinuousUnivariateDistribution = truncated(Normal(1.,0.5), 0, Inf),
    σ_rw_prior::ContinuousUnivariateDistribution = truncated(Normal(0.02,0.01), 0, Inf),
    rwstep    ::Int = 1,
    link = exp,
    invlink = log,
) = begin

    @assert all( dates[1]:Day(1):dates[end] .== dates )
    ntotal  = length(dates)
    nconst  = findfirst(==(startrw), dates)-1
    nobs    = findfirst(==(endobs), dates)
    index   = stepindex(nobs-nconst-1, rwstep)
    nsteps  = length( unique( index ) )

	ConstantPhaseRWModel(
		r0_prior,
		r1_prior,
		σ_rw_prior,
		nconst,
		nobs,
		ntotal,
		nsteps,
		index,
		link,
		invlink
	)
end

(epirt::ConstantPhaseRWModel{Tr0, Tr1, Tσ, Tnc, Tno, Tnt, Tns, Tin, Tli, Til})(
	θ, ::Val{false}
) where {Tr0, Tr1, Tσ, Tnc, Tno, Tnt, Tns, Tin, Tli, Til} = begin
    @unpack rw, r0, r1 = θ
    @unpack nconst, nobs, link = epirt

    rt = Vector{typeof(r0)}(undef, nobs)
    rt[1:nconst]     .= r0
    rt[nconst+1]      = r1
    rt[nconst+2:end]  = link.(rw)
    return rt
end

(epirt::ConstantPhaseRWModel{Tr0, Tr1, Tσ, Tnc, Tno, Tnt, Tns, Tin, Tli, Til})(
	θ, ::Val{true}
) where {Tr0, Tr1, Tσ, Tnc, Tno, Tnt, Tns, Tin, Tli, Til} = begin
    @unpack rw, r0, r1 = θ
    @unpack nconst, ntotal, nobs, link, index = epirt

    rt = Vector{typeof(r0)}(undef, ntotal)
    rt[1:nconst]      .= r0
    rt[nconst+1]       = r1
    rt[nconst+2:nobs]  = link.(rw[index])
    rt[nobs+1:end]    .= rt[nobs]
    return rt
end

Base.rand(rng::Random.AbstractRNG, epirt::ConstantPhaseRWModel) = begin
    @unpack r0_prior, r1_prior, σ_rw_prior, nsteps, invlink = epirt

    ######## draw from univariate prior
    σ_rw = rand(rng, σ_rw_prior )
    r0   = rand(rng, r0_prior )
    r1   = rand(rng, r1_prior )

    ######## draw from random walk prior
    rw0 = invlink(r1)
    rw_prior = RandomWalkDistribution(nsteps, σ_rw)
    rw = rand(rng, rw_prior )
    return (; σ_rw, r0, r1, rw)
end

Distributions.logpdf(epirt::ConstantPhaseRWModel, θ) = begin
    @unpack rw, r0, r1, σ_rw = θ
    @unpack r0_prior, r1_prior, σ_rw_prior, nsteps, invlink = epirt

    rw_prior   = RandomWalkDistribution(nsteps, σ_rw)
    step_prior = Normal(zero(σ_rw), σ_rw) #NOTE same prior as in RadomWalk
    step = rw[1] - invlink(r1)

    target  = logpdf(σ_rw_prior, σ_rw)
    target += logpdf(r0_prior, r0)
    target += logpdf(r1_prior, r1)
    target += logpdf(step_prior, step)
    target += logpdf(rw_prior, rw)
    return target
end

## ===========================================================================

# """
# # Arguments
#
# formula: An object of class formula which determines the linear predictor for the reproduction rates. The left hand side must take the form R(group, date), where group and date variables. group must be a factor vector indicating group membership (i.e. country, state, age cohort), and date must be a vector of class Date. This is syntactic sugar for the reproduction number in the given group at the given date.
#
# link: The link function. This must be either "identity", "log", or a call to scaled_logit.
#
# center: If TRUE then the covariates for the regression are centered to have mean zero. All of the priors are then interpreted as prior on the centered covariates. Defaults to FALSE.
#
# prior: Same as in stan_glm. In addition to the rstanarm provided priors, a shifted_gamma can be used. Note: If autoscale=TRUE in the call to the prior distribution then automatic rescaling of the prior may take place.
#
# prior_intercept: Same as in stan_glm. Prior for the regression intercept (if it exists).
#
# prior_covariance: Same as in stan_glmer. Only used if the formula argument specifies random effects.
#
# # Details
#
# epirt has a formula argument which defines the linear predictor, an argument link defining the link function, and additional arguments to specify priors on parameters making up the linear predictor.
# A general R formula gives a symbolic description of a model. It takes the form y ~ model, where y is the response and model is a collection of terms separated by the + operator. model fully defines a linear predictor used to predict y. In this case, the “response” being modeled are reproduction numbers which are unobserved. epirt therefore requires that the left hand side of the formula takes the form R(group, date), where group and date refer to variables representing the region and date respectively. The right hand side can consist of fixed effects, random effects, and autocorrelation terms.
#
# """
# function epirt(
#   formula,
#   link = "log",
#   center = FALSE,
#   prior = Normal(1, 0.5),
#   prior_intercept = Normal(1., 0.5),
#   prior_covariance = rstanarm::decov(scale = 0.5)
# ) end
#
#
# """
# A call to rw can be used in the 'formula' argument of epim, allowing random walks for the reproduction number.
#
# time:	An optional name defining the random walk time periods for each date and group. This must be a column name found in the data argument to epim. If not specified, determined by the dates column implied by the formula argument to epim is used.
#
# gr: Same as for time, except this defines the grouping to use for the random walks. A separate walk is defined for each group. If not specified a common random walk is used for all groups.
#
# prior_scale: The steps of the walks are independent zero mean normal with an unknown scale hyperparameter. This scale is given a half-normal prior. prior_scale sets the scale parameter of this hyperprior.
# """
# rw(time = "date", gr = "country", prior_scale = 0.2) = RWModel(time, gr, prior_scale)

# epirt(
#   formula,
#   link = "log",
#   prior = Normal(1, 0.5),
#   prior_intercept = Normal(1., 0.5),
#   prior_covariance = rstanarm::decov(scale = 0.5),
#   ...
# )


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, BinaryBuilderBase, Pkg, Base.BinaryPlatforms

name = "xtrx"
version = v"0.1.2"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/xtrx-sdr/libxtrxdsp", "eec28640c0ebd5639b642f07b310a0a0d02d9834"),
    GitSource("https://github.com/xtrx-sdr/libxtrxll", "1b6eddfbedc700efb6f7e3c3594e43ac6ff29ea4"),
    GitSource("https://github.com/xtrx-sdr/libxtrx", "acb0b1cf7ab92744034767a04c1d4b4c281b840f"),
    GitSource("https://github.com/xtrx-sdr/liblms7002m", "b07761b7386181f0e6a35158456b75bce14f2aca"),
    DirectorySource("bundled"),
]

# Bash recipe for building across all platforms
script = raw"""

# Apply our patches
cd ${WORKSPACE}/srcdir
for f in ${WORKSPACE}/srcdir/patches/*.patch; do
    atomic_patch -p1 ${f}
done

# liblms7002m requires python3 and cheetah to perform build-time templating
apk add python3 py3-cheetah

if [[ ${bb_target} == x86_64-* ]] || [[ ${bb_target} == i686-* ]]; then
    XTRX_ARCH=x86
elif [[ ${bb_target} == aarch64-* ]] || [[ ${bb_target} == arm* ]]; then
    XTRX_ARCH=arm
fi

CMAKE_ARGS=(
    "-DCMAKE_INSTALL_PREFIX=${prefix}"
    "-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN}"
    "-DCMAKE_BUILD_TYPE=Release"

    # CC_ARCH isn't set properly; just override
    "-DFORCE_ARCH=${XTRX_ARCH}"

    # Install utilities to `bin`
    "-DXTRXDSP_UTILS_DIR=${bindir}"
    "-DXTRXLL_UTILS_DIR=${bindir}"
)

mkdir ${WORKSPACE}/srcdir/liblms7002m_build
cd ${WORKSPACE}/srcdir/liblms7002m_build
cmake ${CMAKE_ARGS[@]} ${WORKSPACE}/srcdir/liblms7002m
make -j${nproc}
make install

mkdir ${WORKSPACE}/srcdir/libxtrxdsp_build
cd ${WORKSPACE}/srcdir/libxtrxdsp_build
cmake ${CMAKE_ARGS[@]} ${WORKSPACE}/srcdir/libxtrxdsp
make -j${nproc}
make install

mkdir ${WORKSPACE}/srcdir/libxtrxll_build
cd ${WORKSPACE}/srcdir/libxtrxll_build
cmake ${CMAKE_ARGS[@]} -DENABLE_PCIE=ON ${WORKSPACE}/srcdir/libxtrxll
make -j${nproc}
make install

mkdir ${WORKSPACE}/srcdir/libxtrx_build
cd ${WORKSPACE}/srcdir/libxtrx_build
cmake ${CMAKE_ARGS[@]} ${WORKSPACE}/srcdir/libxtrx
make -j${nproc}
make install
install_license ${WORKSPACE}/srcdir/libxtrx/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# Disable everything except glibc Linux for now, and only intel/arm processors
filter!(p -> Sys.islinux(p) && libc(p) == "glibc" && proc_family(p) ∈ ("intel", "arm"), platforms)
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("test_xtrxll", :test_xtrxll),
    ExecutableProduct("test_xtrxflash", :test_xtrxflash),
    ExecutableProduct("test_xtrx", :test_xtrxflash),
    LibraryProduct("libxtrx", :libxtrx),
    LibraryProduct("libxtrxll", :libxtrxll),
    LibraryProduct("libxtrxdsp", :libxtrxdsp),
    LibraryProduct("liblms7compact", :liblms7compact),
    LibraryProduct("libXTRXSupport", :libXTRXSupport, ["lib/SoapySDR/modules0.8"]),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
    Dependency("soapysdr_jll"),
]

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

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
module IAE

export iae

"""
```julia
iae(error)
```
IAE (Integral Absolute Error) from error vector.
Inputs:
error -> error vector.
"""
function iae(error::Vector{T}) where T
    sum(abs.(error))
end

"""
```julia
iae(error)
```
IAE (Integral Absolute Error) from vector of error vector.
Inputs:
error -> vector of error vector.
"""
function iae(error::Vector{Vector{T}}) where T
    map(x->sum(abs.(x)), error)
end

"""
```julia
iae(x, t, desired)
```
IAE (Integral Absolute Error) from error betwwen desired value and controler output.
Inputs:
x -> controler output vector;
t -> time vector;
desired -> desired value function.
"""
function iae(x::Vector{T}, t::Vector{T}, desired::Function) where T
    x = x .- map(desired, t)
    sum(abs.(x))
end

end


text:
using Printf, BenchmarkTools, Trixi

function run_benchmarks(benchmark_run; levels=0:5, polydeg=3)
  runtimes = zeros(length(levels))
  for (idx,initial_refinement_level) in enumerate(levels)
    result = benchmark_run(; initial_refinement_level, polydeg)
    display(result)
    runtimes[idx] = result |> median |> time # in nanoseconds
  end
  return (; levels, runtimes, polydeg)
end

function tabulate_benchmarks(args...; kwargs...)
  result = run_benchmarks(args...; kwargs...)
  println("#Elements | Runtime in seconds")
  for (level,runtime) in zip(result.levels, result.runtimes)
    @printf("%9d | %.2e\n", 4^level, 1.0e-9 * runtime)
  end
  for (level,runtime) in zip(result.levels, result.runtimes)
    @printf("%.16e\n", 1.0e-9 * runtime)
  end
end

function benchmark_euler(; initial_refinement_level=1, polydeg=3)

  γ = 1.4
  equations = CompressibleEulerEquations2D(γ)

  surface_flux = flux_ranocha
  volume_flux  = flux_ranocha
  solver = DGSEM(polydeg, surface_flux, VolumeIntegralFluxDifferencing(volume_flux))

  coordinates_min = (-2.0, -2.0)
  coordinates_max = ( 2.0,  2.0)
  mesh = TreeMesh(coordinates_min, coordinates_max,
                  initial_refinement_level=initial_refinement_level,
                  n_cells_max=100_000)

  semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_weak_blast_wave, solver)

  t0 = 0.0
  u0 = compute_coefficients(t0, semi)
  du = similar(u0)

  @benchmark Trixi.rhs!($du, $u0, $semi, $t0)
end

# versioninfo(verbose=true)
@show Threads.nthreads()
tabulate_benchmarks(benchmark_euler, levels=0:8)


text:
function est_mmvmnfa(model::MMVMNFAModel, Mquantities::MMVMNFAMutableQuantities, Squantities::MMVMNFAStaticQuantities)

    #p = size(Y,2)
    #n = size(Y,1)
    model.w_esteps = zeros(Squantities.n,Squantities.n_esteps,model.g)
    model.Ψ = reshape(repeat(Squantities.psi_initial,model.g),(size(Squantities.psi_initial)[1],model.g))
    Squantities.logL_calc!(model, Mquantities, Squantities)
    if model.status == "failed"
        model.status = "Failed: got a NaN in a log-likelihod or E-step for Ψ after initialisation"
        model.niter = 0
        return
    elseif model.status == "failedE"
        model.status = "Failed: E-step for Ψ failed at initialisation" #probably caused by a bad numerical algorithm ... such as finite differences for GH
        model.niter = 0
        return
    end

    model.niter = Squantities.itmax
    Mquantities.logL_old = model.logL
    model.logL_history[1] = model.logL
    @inbounds for niter in 1:Squantities.itmax
        model.niter = niter
        Squantities.Mstep_calc!(model, Mquantities, Squantities) #Need to update this function
        if model.status == "failed"
            model.status = "Failed: got a complex eigenvalue in the M-steps after $niter iterations"
            model.niter = niter
            return
        end

        Squantities.logL_calc!(model, Mquantities, Squantities)
        if model.status == "failed"
            model.status = "Failed: got a NaN in a log-likelihod or E-step for Ψ after $niter iterations"
            model.niter = niter
            break
        elseif model.status == "failedE"
            model.status = "Failed: E-step for Ψ failed after $niter iterations" #probably caused by a bad numerical algorithm ... such as finite differences for GH
            model.niter = niter
            break
        end
        if ((Squantities.conv_measure == "diff") && (abs(model.logL - Mquantities.logL_old) < Squantities.tol))
            model.niter = niter
            break
        end

        if ((Squantities.conv_measure == "ratio") && (abs((model.logL - Mquantities.logL_old)/model.logL) < Squantities.tol))
            model.niter = niter
            break
        end

        Mquantities.logL_old = model.logL
        model.logL_history[niter + 1] = model.logL
    end #for
    #model.niter = Squantities.itmax
end

function est_msmnfa(model::MMVMNFAModel, Mquantities::MMVMNFAMutableQuantities, Squantities::MMVMNFAStaticQuantities)

    model.w_esteps = zeros(Squantities.n,Squantities.n_esteps,model.g)
    model.Ψ = reshape(repeat(Squantities.psi_initial,model.g),(size(Squantities.psi_initial)[1],model.g))
    model.beta = zeros(Squantities.p,model.g)

    Squantities.logL_calc!(model, Mquantities, Squantities)
    if model.status == "failed"
        model.status = "Failed: got a NaN in a log-likelihod or E-step for Ψ after initialisation"
        return
    end
    Mquantities.logL_old = model.logL
    model.logL_history[1] = model.logL
    @inbounds for niter in 1:Squantities.itmax
        Squantities.Mstep_calc!(model, Mquantities, Squantities)
        if model.status == "failed"
            model.status = "Failed: got a complex eigenvalue in the M-steps after $niter iterations"
            return
        end

        Squantities.logL_calc!(model, Mquantities, Squantities)
        if model.status == "failed"
            model.status = "Failed: got a NaN in a log-likelihod or E-step for Ψ after $niter iterations"
            break
        end
        if ((Squantities.conv_measure == "diff") && (abs(model.logL - Mquantities.logL_old) < Squantities.tol))
            model.niter = niter
            break
        end

        if ((Squantities.conv_measure == "ratio") && (abs((model.logL - Mquantities.logL_old)/model.logL) < Squantities.tol))
            model.niter = niter
            break
        end
        model.logL_history[niter+1] = model.logL
        Mquantities.logL_old = model.logL

    end #for
    #if model == "failedModel"
    #    return "failedModel"
    #else
    #    return model
    #end
end

function est_mfa(model::MMVMNFAModel, Mquantities::MMVMNFAMutableQuantities, Squantities::MMVMNFAStaticQuantities)

    #p = size(Y,2)
    #n = size(Y,1)
    model.beta =  zeros(Squantities.p,model.g)
    #init_para.beta = zeros(p,init_para.g)

    Squantities.logL_calc!(model, Mquantities, Squantities)
    if model.status == "failed"
        model.status = "Failed: got a NaN in a log-likelihod or E-step for Ψ after initialisation"
        return
    end
    Mquantities.logL_old = model.logL
    model.logL_history[1] = model.logL
    @inbounds for niter in 1:Squantities.itmax
        Squantities.Mstep_calc!(model, Mquantities, Squantities) #Need to update this function
        if model.status == "failed"
            model.status = "Failed: got a complex eigenvalue in the M-steps after $niter iterations"
            return
        end

        Squantities.logL_calc!(model, Mquantities, Squantities)
        if model.status == "failed"
            model.status = "Failed: got a NaN in a log-likelihod or E-step for Ψ after $niter iterations"
            return
        end
        if ((Squantities.conv_measure == "diff") && (abs(model.logL - Mquantities.logL_old) < Squantities.tol))
            model.niter = niter
            model.logL_history[niter + 1] = model.logL
            break
        end

        if ((Squantities.conv_measure == "ratio") && (abs((model.logL - Mquantities.logL_old)/model.logL) < Squantities.tol))
            model.niter = niter
            model.logL_history[niter + 1] = model.logL
            break
        end

        Mquantities.logL_old = model.logL
        model.logL_history[niter + 1] = model.logL
    end #for
end

max_stars_repo_name: john-c-davey/FactorMixtures

max_issues_repo_name: john-c-davey/FactorMixtures

max_forks_repo_name: john-c-davey/FactorMixtures


text:
# list of point position + velocity.  need to figure out what they
# "spell"...maybe can use min distance from origin as a heuristic here?  it
# seems like all the inputs are large but have opposite signs from their
# velocies (on a scale of like 10_000x or something).  or sum of pairwise
# distances?
#
# or just find the relative position/velocity of each pair of point and use that
# to get the overall (sum) pariwise distance as a function of time then find the
# minimum...
#
# ((a_0 - b_0) + (a_v - b_v)*t)^2 = Δ_0^2 + 2 Δ_0 Δ_v t + Δ_v^2 t^2
# derivative of that is 2Δ_0Δ_v + 2Δ_v^2 t => max/min at t = Δ0 / Δv
#
# that's just for one.  the quadratic coefficients sum but need to do the sum
# first though...

mutable struct Star
    x::Tuple{Int,Int}
    v::Tuple{Int,Int}
end

function parseline(line)
    regex = r"position=< *(-?\d*), *(-?\d*)> velocity=< *(-?\d*), *(-?\d*)>"
    x, y, vx, vy = parse.(Int, match(regex, line).captures)
    Star((x,y), (vx, vy))
end

test_input = parseline.(split(strip(
    """position=< 9,  1> velocity=< 0,  2>
position=< 7,  0> velocity=<-1,  0>
position=< 3, -2> velocity=<-1,  1>
position=< 6, 10> velocity=<-2, -1>
position=< 2, -4> velocity=< 2,  2>
position=<-6, 10> velocity=< 2, -2>
position=< 1,  8> velocity=< 1, -1>
position=< 1,  7> velocity=< 1,  0>
position=<-3, 11> velocity=< 1, -2>
position=< 7,  6> velocity=<-1, -1>
position=<-2,  3> velocity=< 1,  0>
position=<-4,  3> velocity=< 2,  0>
position=<10, -3> velocity=<-1,  1>
position=< 5, 11> velocity=< 1, -2>
position=< 4,  7> velocity=< 0, -1>
position=< 8, -2> velocity=< 0,  1>
position=<15,  0> velocity=<-2,  0>
position=< 1,  6> velocity=< 1,  0>
position=< 8,  9> velocity=< 0, -1>
position=< 3,  3> velocity=<-1,  1>
position=< 0,  5> velocity=< 0, -1>
position=<-2,  2> velocity=< 2,  0>
position=< 5, -2> velocity=< 1,  2>
position=< 1,  4> velocity=< 2,  1>
position=<-2,  7> velocity=< 2, -2>
position=< 3,  6> velocity=<-1, -1>
position=< 5,  0> velocity=< 1,  0>
position=<-6,  0> velocity=< 2,  0>
position=< 5,  9> velocity=< 1, -2>
position=<14,  7> velocity=<-2,  0>
position=<-3,  6> velocity=< 2, -1>
"""), '\n'))

input = parseline.(readlines("10.input"))

move(star::Star, t::Int) = Star(star.x .+ t .* star.v, star.v)

# compute the coefficients of the quadratic expression of squared distance for a pair
coefs(a::Star, b::Star) = (sum((a.x .- b.x).^2),
                           sum(2 .* (a.x .- b.x) .* (a.v .- b.v)),
                           sum((a.v .- b.v).^2)

function tmin(input)
    a,b,c = reduce((x,y)->x.+y, coefs(a,b) for a in input for b in input)
    tmin = round(Int, -b / (2c))
end

function draw(stars::Vector{Star})
    mins, maxs = foldl(((minx,maxx),x) -> (min.(minx,x), max.(maxx,x)),
                       s.x for s in stars,
                       init=(typemax(Int), typemin(Int)))
    img = fill(' ', maxs .- mins .+ 1)
    for s in stars
        img[(s.x .- mins .+ 1)...] = '█'
    end
    cols, rows = axes(img)
    println()
    for r in rows
        for c in cols
            print(img[c,r])
        end
        println()
    end
    img
end

star1(input) = draw(move.(input, tmin(input)))

star1(test_input);

star1(input);

# star 2:
tmin(input)


text:
## labelmap

function labelmap!(dict::Dict{T,Vector{Int}}, idx::Int, elem) where {T}
    labelmap!(dict, idx, T(elem))
end

function labelmap!(dict::Dict{T,Vector{Int}}, idx::Int, elem::T) where {T}
    if !haskey(dict, elem)
        push!(dict, elem => [idx])
    else
        push!(dict[elem], idx)
    end
    dict
end

function labelmap!(dict::Dict{T,Vector{Int}}, idx::AbstractVector{Int}, itr::S) where {T,S}
    for (i, elem) in zip(idx, itr)
        labelmap!(dict, i, elem)
    end
    dict
end

function labelmap(itr::T) where {T}
    dict = Dict{eltype(T),Vector{Int}}()
    for (idx, elem) in enumerate(itr)
        labelmap!(dict, idx, elem)
    end
    dict
end

## labelfreq

function labelfreq!(dict::Dict{T,Int}, elem) where {T}
    labelfreq!(dict, T(elem))
end

function labelfreq!(dict::Dict{T,Int}, elem::T) where {T}
    cnt = get(dict, elem, 0)
    dict[elem] = cnt + 1
    dict
end

function labelfreq!(dict::Dict{T,Int}, iter::AbstractVector) where {T}
    for elem in iter
        labelfreq!(dict, T(elem))
    end
    dict
end

function rarelabel(labels::Array,threshold::Integer,replace_label::Symbol)
    label_freq = labelfreq(labels)
    for i = 1:length(labels)
        if label_freq[labels[i]]<=threshold
            labels[i]=replace_label
        end
    end
    return labels
end

labelfreq(itr::T) where {T} = countmap(itr, alg = :dict)::Dict{eltype(T),Int}
labelfreq(dict::Dict{T,Vector{Int}})  where {T} = Dict(k => length(v) for (k,v) in dict)::Dict{T,Int}

## Errors

labelmap(A::AbstractMatrix)  = throw(ArgumentError("labelmap not supported for matrices"))
labelfreq(A::AbstractMatrix) = throw(ArgumentError("labelfreq not supported for matrices"))

## General Dict

labelenc(x::Dict) = labelenc(label(x))
nlabel(x::Dict) = length(keys(x))
label(x::Dict) = _arrange_label(collect(keys(x)))


text:
using ..Atmos
using ..ConfigTypes
using ..DiagnosticsMachine

@diagnostics_group(
    "AtmosLESDefault",      # name
    AtmosLESConfigType,     # configuration type
    Nothing,                # params type
    (_...) -> nothing,      # initialization function
    NoInterpolation,        # if/when to interpolate
    # various horizontal averages
    u,
    v,
    w,
    rho,
    temp,
    pres,
    thd,
    et,
    ei,
    ht,
    hi,
    w_ht_sgs,
    # moisture related
    qt,
    ql,
    qi,
    qv,
    thv,
    thl,
    w_qt_sgs,
    # for variances and co-variances
    uu,
    vv,
    ww,
    www,
    eiei,
    wu,
    wv,
    wrho,
    wthd,
    wei,
    qtqt,
    thlthl,
    wqt,
    wql,
    wqi,
    wqv,
    wthv,
    wthl,
    qtthl,
    qtei,
    #cld_top, TODO
    #cld_base, TODO
    #cld_cover, TODO
    #lwp, TODO
    #rwp, TODO
)


text:
function raw_to_counts(raw)
    x_values = convert(Array{Float64}, unique(sort(raw[:, 2])))

    shots = Array(Array{UInt32, 1}, length(x_values))
    for i in 1:length(x_values)
        shots[i] = UInt32[]
    end

    for i in 1:(size(raw)[1])
        x_index = searchsorted(x_values, raw[i, 2])[1]
        for val in raw[i, 4:end]
            if val == ""
                break
            end

            counts = convert(UInt32, val)
            if counts & (1 << 31) != 0
                continue
            end

            push!(shots[x_index], counts)
        end
    end

    x_values, shots
end

function last_complete_row_2d(raw, y_column)
    lines_per_y = 0
    for (i, y) in enumerate(slice(raw, :, y_column))
        if y != raw[1, y_column]
            lines_per_y = i - 1
            break
        end
    end

    row_count = size(raw)[1]
    if lines_per_y == 0
       lines_per_y = row_count
    end

    if row_count % lines_per_y == 0
        return row_count, lines_per_y
    end

    for (i, y) in enumerate(slice(raw, row_count:-1:1, y_column))
        if y != raw[row_count, y_column]
            current = row_count + 1 - i
            return current, lines_per_y
        end
    end

    error("Dataset has only one y value, cannot determine whether it is complete.")
end

function raw_to_counts_2d(raw, lines_per_y::Int64)
    @assert(size(raw)[1] % lines_per_y == 0, "Inconsistent number of measurements per iterator value")

    y_values = convert(Array{Float64}, sort(raw[1:lines_per_y:end, 3]))
    x_values = convert(Array{Float64}, unique(sortrows(slice(raw, 1:lines_per_y, 1:2))[:, 2]))

    shots = Array(Array{UInt32, 1}, length(x_values), length(y_values))
    for j in 1:length(y_values)
        for i in 1:length(x_values)
            shots[i, j] = UInt32[]
        end
    end

    for i in 1:(size(raw)[1])
        x_index = searchsorted(x_values, raw[i, 2])[1]
        y_index = searchsorted(y_values, raw[i, 3])[1]
        for val in raw[i, 5:end]
            if val == ""
                break
            end

            counts = convert(UInt32, val)
            if counts & (1 << 31) != 0
                continue
            end

            push!(shots[x_index, y_index], counts)
        end
    end

    x_values, y_values, shots
end

function walk_match(dir, prefix, suffix)
    if isempty(dir)
        dir = "."
    end
    entries = readdir(dir)
    for e in entries
        if !startswith(e, prefix)
            continue
        end

        path = joinpath(dir, e)
        if endswith(path, suffix)
            return path
        end

        # Recurse into child directories.
        if isdir(path)
            p = walk_match(path, "", suffix)
            if !isempty(p)
                return p
            end
        end
    end
    return ""
end

function find_file(datapath, suffix)
    if endswith(datapath, suffix)
        return datapath
    end

    dirname, start = splitdir(datapath)
    path = walk_match(dirname, start, suffix)
    if isempty(path)
        error("No file found in '$datapath' ending in '$suffix'.")
    end
    return path
end

function read_shots_log(datapath)
    raw_to_counts(readcsv(find_file(datapath, "_shots_log.txt")))
end

function read_shots_log_2d(datapath; y_reversed=false)
    data = readcsv(find_file(datapath, "_2D_shots_log.txt"))
    last_row, lines_per_y = last_complete_row_2d(data, 3)

    if y_reversed
        clipped_data = data[last_row:-1:1, :]
    else
        clipped_data = data[1:last_row, :]
    end

    raw_to_counts_2d(clipped_data, lines_per_y)
end

function read_thresholded_2d(datapath)
    # Read in file, skipping header line
    data = readcsv(find_file(datapath, "_plot_0.csv"))[2:end, :]

    last_row, lines_per_y = last_complete_row_2d(data, 2)

    y_values = convert(Array{Float64}, sort(data[1:lines_per_y:end, 2]))
    x_values = convert(Array{Float64}, unique(sort(data[1:lines_per_y, 1])))

    signal = zeros(Float64, length(x_values), length(y_values))
    for i in 1:(size(data)[1])
        x_index = searchsorted(x_values, data[i, 1])[1]
        y_index = searchsorted(y_values, data[i, 2])[1]
        signal[x_index, y_index] += data[i, 3]
    end

    signal ./= (lines_per_y / length(x_values))

    x_values, y_values, signal
end


text:
include("./imports.jl")

@testset ExtendedTestSet "DataAugmentation.jl" begin
    @testset ExtendedTestSet "base.jl" begin
        include("base.jl")
    end
    @testset ExtendedTestSet "testing.jl" begin
        include("testing.jl")
    end
    @testset ExtendedTestSet "buffered.jl" begin
        include("buffered.jl")
    end
    @testset ExtendedTestSet "projective/base.jl" begin
        include("projective/base.jl")
    end
    @testset ExtendedTestSet "projective/compose.jl" begin
        include("projective/compose.jl")
    end
    @testset ExtendedTestSet "projective/crop.jl" begin
        include("projective/crop.jl")
    end
    @testset ExtendedTestSet "projective/affine.jl" begin
        include("projective/affine.jl")
    end
    @testset ExtendedTestSet "projective/warp.jl" begin
        include("projective/warp.jl")
    end
    @testset ExtendedTestSet "items/image.jl" begin
        include("items/image.jl")
    end
    @testset ExtendedTestSet "items/keypoints.jl" begin
        include("items/keypoints.jl")
    end
    @testset ExtendedTestSet "items/mask.jl" begin
        include("items/mask.jl")
    end
    @testset ExtendedTestSet "preprocessing.jl" begin
        include("preprocessing.jl")
    end
    @testset ExtendedTestSet "oneof.jl" begin
        include("oneof.jl")
    end
    @testset ExtendedTestSet "visualization.jl" begin
        include("visualization.jl")
    end
    @testset ExtendedTestSet "rowtransforms.jl" begin
        include("rowtransforms.jl")
    end
end


text:
function propagation(settings::PyObject, data::PyObject, x_obs::Array{Float64, 1}, msap_source, r, x, z, c_bar, rho_0, I_0, beta)

    # Impedance at the observer
    I_0_obs = 409.74

    # Size of inputs
    n_t = size(r)[1]

    # Get type of input variables
    T = eltype(msap_source)
    msap_prop = zeros(T, (n_t, settings.N_f))

    # Apply spherical spreading and characteristic impedance effects to the MSAP
    # Source: Zorumski report 1982 part 1. Chapter 5.1 Equation 1
    z = max.(z, 0)
    if settings.direct_propagation
        msap_prop = msap_source .* ((settings.r_0^2 ./ r.^2) .* (I_0_obs./I_0))
    else
        msap_prop = msap_source
    end

    # Split mean-square acoustic pressure in frequency sub-bands
    if settings.N_b > 1
        msap_sb = split_subbands(settings, msap_prop)
    end

    # Apply atmospheric absorption on sub-bands
    if settings.absorption
        # Calculate average absorption factor between observer and source
        f_abs = LinearInterpolation((data.abs_alt, data.abs_freq), data.abs)
        alpha_f = f_abs.(z.*ones(T, (1, settings.N_b*settings.N_f)), reshape(data.f_sb, (1,settings.N_f*settings.N_b)).*ones(T, (n_t,1)))

        # Calculate absorption (convert dB to Np: 1dB is 0.115Np)
        # Source: Zorumski report 1982 part 1. Chapter 5.1 Equation 14
        msap_prop = msap_sb .* exp.(-2 * 0.115 * alpha_f .* (r .- settings.r_0))
    end

    # Apply ground effects on sub-bands
    if settings.groundeffects
        # Empirical lateral attenuation for microphone on sideline
        if (settings.lateral_attenuation == true) && (x_obs[2] != 0)
            # Lateral attenuation factor
            Lambda = lateral_attenuation(settings, beta, x_obs)

            # Compute elevation angle from "center-line" observer.
            r_cl = sqrt.((x.-x_obs[1]).^2 .+ 1. ^2 .+ z.^2)
            beta_cl = asin.(z./r_cl) * 180. / pi

            # Ground reflection factor for center-line
            G = ground_reflections(settings, data, x_obs, r_cl, beta_cl, c_bar, rho_0)

            # Apply ground effects
            msap_prop = msap_prop .* (G .* Lambda)
        else
            # Compute ground effects factor
            G = ground_reflections(settings, data, x_obs, r, beta, c_bar, rho_0)

            # Multiply subbands with G
            msap_prop = msap_prop .* G
        end
    end

    # Recombine the mean-square acoustic pressure in the frequency sub-bands
    if settings.N_b > 1
        msap_prop = transpose(reshape(sum(reshape(transpose(msap_prop), (settings.N_b, settings.N_f, n_t)), dims=1), (settings.N_f, n_t)))
    end

    return msap_prop
end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
using JGUI
using Base.Test

w = window(title="notebook")
nb = notebook(w); push!(w, nb)
connect(nb, "valueChanged", value -> println("On tab $value"))

b1 = button(nb, "one")
b2 = button(nb, "two")
b3 = button(nb, "three")

push!(nb, b1, "one")
push!(nb, b2, "two")
push!(nb, b3, "three")

setValue(nb, 3)
@assert nb[:value] == 3

pop!(nb, b2)
@assert length(nb) == 2
@assert nb[1][:value] == "one"

raise(w)


text:
# EXCEPTIONS

function squareroot(x)
    try
        return sqrt(x)
    catch
        return sqrt(x + 0im)
    end
end

squareroot(2)
squareroot(-2)

max_stars_repo_name: milankl/misc

max_issues_repo_name: milankl/misc

max_forks_repo_name: milankl/misc


text:
#=
district_cooling:
- Julia version:
- Author: qgis1
- Date: 2019-03-28
=#
module district_cooling
    using JuMP
    using Cbc

    using CSV
    using DelimitedFiles
    csv_sep = ','

    function district_cooling_simulation(input_folder, output_folder, tech_infos, network_id, print_function = print)
        m = Model(with_optimizer(Cbc.Optimizer, logLevel=1, seconds=600))

        ####### Load distributions #########

        DEM=readdlm(string(input_folder, "\\DEM_time.csv"), csv_sep)
        #DEM=DEM[4381:6570]

        Electricity_price=readdlm(string(input_folder, "\\Electricity_price_time.csv"), csv_sep)
        #Electricity_price=Electricity_price[4381:6570]

        Heat_exchanger_specific=readdlm(string(input_folder, "\\Heat_exchanger_specific_time.csv"), csv_sep)
        #Heat_exchanger_specific=Heat_exchanger_specific[4381:6570]

        Heat_exchanger_specific_cool=readdlm(string(input_folder, "\\Heat_exchanger_specific_time_cool.csv"), csv_sep)
        #Heat_exchanger_specific_cool=Heat_exchanger_specific_cool[4381:6570]

        eta_HP = readdlm(string(input_folder, "\\eta_HP_1.csv"), csv_sep)
        #eta_HP = eta_HP[4381:6570]
        eta_HP_2 = readdlm(string(input_folder, "\\eta_HP_2.csv"), csv_sep)
        #eta_HP_2 = eta_HP_2[4381:6570]
        eta_HP_3 = readdlm(string(input_folder, "\\eta_HP_3.csv"), csv_sep)
        #eta_HP_3 = eta_HP_3[4381:6570]

        eta_HP_waste_heat_1=readdlm(string(input_folder, "\\COP_heat_pump_temperature_group_1.csv"), csv_sep)
        #eta_HP_waste_heat_1=eta_HP_waste_heat_1[4381:6570]
        eta_HP_waste_heat_2=readdlm(string(input_folder, "\\COP_heat_pump_temperature_group_2.csv"), csv_sep)
        #eta_HP_waste_heat_2=eta_HP_waste_heat_2[4381:6570]
        eta_HP_waste_heat_3=readdlm(string(input_folder, "\\COP_heat_pump_temperature_group_3.csv"), csv_sep)
        #eta_HP_waste_heat_3=eta_HP_waste_heat_3[4381:6570]

        Available_waste_heat_heat_pump_1=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_temperature_group_1.csv"), csv_sep)
        #Available_waste_heat_heat_pump_1=Available_waste_heat_heat_pump_1[4381:6570]
        Available_waste_heat_heat_pump_2=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_temperature_group_2.csv"), csv_sep)
        #Available_waste_heat_heat_pump_2=Available_waste_heat_heat_pump_2[4381:6570]
        Available_waste_heat_heat_pump_3=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_temperature_group_3.csv"), csv_sep)
        #Available_waste_heat_heat_pump_3=Available_waste_heat_heat_pump_3[4381:6570]

        ####### Technology input data ########

        ## Heat only boiler ##

        # Installed capacity [MW]
        P_max_HOB = tech_infos["P_max_HOB"] #2
        # Technical minimum [MW]
        Tecnical_minimum_HOB= tech_infos["Tecnical_minimum_HOB"] #0.2
        P_min_HOB = Tecnical_minimum_HOB*P_max_HOB
        # Efficiency [-]
        eta_HOB = tech_infos["eta_HOB"] #0.8
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HOB= tech_infos["Percentage_ramp_up_down_HOB"] #0.5
        ramp_up_down_HOB = P_max_HOB*Percentage_ramp_up_down_HOB
        # Variable cost [EUR/MWh]
        cost_var_HOB = tech_infos["cost_var_HOB"] #5.4
        # Fuel cost [EUR/MWh]
        cost_fuel_HOB = tech_infos["cost_fuel_HOB"] #10

        ## Heat only boiler 2 ##

        # Installed capacity [MW]
        P_max_HOB_2 = tech_infos["P_max_HOB_2"] #2
        # Technical minimum [MW]
        Tecnical_minimum_HOB_2= tech_infos["Tecnical_minimum_HOB_2"] #0.3
        P_min_HOB_2 = Tecnical_minimum_HOB_2*P_max_HOB_2
        # Efficiency [-]
        eta_HOB_2 = tech_infos["eta_HOB_2"] #0.8
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HOB_2= tech_infos["Percentage_ramp_up_down_HOB_2"] #0.5
        ramp_up_down_HOB_2 = P_max_HOB_2*Percentage_ramp_up_down_HOB_2
        # Variable cost [EUR/MWh]
        cost_var_HOB_2 = tech_infos["cost_var_HOB_2"] #5.4
        # Fuel cost [EUR/MWh]
        cost_fuel_HOB_2 = tech_infos["cost_fuel_HOB_2"] #12.5

        ## Heat only boiler 3 ##

        # Installed capacity [MW]
        P_max_HOB_3 = tech_infos["P_max_HOB_3"] #2
        # Technical minimum [MW]
        Tecnical_minimum_HOB_3= tech_infos["Tecnical_minimum_HOB_3"] #0.2
        P_min_HOB_3 = Tecnical_minimum_HOB_3*P_max_HOB_3
        # Efficiency [-]
        eta_HOB_3 = tech_infos["eta_HOB_3"] #0.8
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HOB_3= tech_infos["Percentage_ramp_up_down_HOB_3"] #0.6
        ramp_up_down_HOB_3 = P_max_HOB_3*Percentage_ramp_up_down_HOB_3
        # Variable cost [EUR/MWh]
        cost_var_HOB_3 = tech_infos["cost_var_HOB_3"] #3.4
        # Fuel cost [EUR/MWh]
        cost_fuel_HOB_3 = tech_infos["cost_fuel_HOB_3"] #20

        ## Heat pump ##

        # Installed capacity [MW]
        P_max_HP = tech_infos["P_max_HP"] #2
        # Technical minimum [MW]
        Technical_minimum_HP= tech_infos["Technical_minimum_HP"] #0.2
        P_min_HP = Technical_minimum_HP*P_max_HP
        # Efficiency [-]
        eta_HP = eta_HP
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP= tech_infos["Percentage_ramp_up_down_HP"] #0.9
        ramp_up_down_HP = Percentage_ramp_up_down_HP*P_max_HP
        # Variable cost [EUR/MWh]
        cost_var_HP = tech_infos["cost_var_HP"] #0.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP = Electricity_price

        ## Heat pump 2 ##

        # Installed capacity [MW]
        P_max_HP_2 = tech_infos["P_max_HP_2"] #2
        # Technical minimum [MW]
        Technical_minimum_HP_2= tech_infos["Technical_minimum_HP_2"] #0.1
        P_min_HP_2 = Technical_minimum_HP_2*P_max_HP_2
        # Efficiency [-]
        eta_HP_2 = eta_HP_2
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_2= tech_infos["Percentage_ramp_up_down_HP_2"] #0.95
        ramp_up_down_HP_2 = Percentage_ramp_up_down_HP_2*P_max_HP_2
        # Variable cost [EUR/MWh]
        cost_var_HP_2 = tech_infos["cost_var_HP_2"] #0.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_2 = Electricity_price

        ## Heat pump 3 ##

        # Installed capacity [MW]
        P_max_HP_3 = tech_infos["P_max_HP_3"] #2
        # Technical minimum [MW]
        Technical_minimum_HP_3= tech_infos["Technical_minimum_HP_3"] #0.1
        P_min_HP_3 = Technical_minimum_HP_3*P_max_HP_3
        # Efficiency [-]
        eta_HP_3 = eta_HP_3
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_3= tech_infos["Percentage_ramp_up_down_HP_3"] #0.95
        ramp_up_down_HP_3 = Percentage_ramp_up_down_HP_3*P_max_HP_3
        # Variable cost [EUR/MWh]
        cost_var_HP_3 = tech_infos["cost_var_HP_3"] #0.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_3 = Electricity_price

        ## Cogeneration ##

        # Installed capacity [MW]
        P_max_CHP = tech_infos["P_max_CHP"] #0
        # Technical minimum [MW]
        Technical_minimum_CHP= tech_infos["Technical_minimum_CHP"] #0.0
        P_min_CHP = Technical_minimum_CHP*P_max_CHP
        # Efficiency [-]
        eta_CHP_th = tech_infos["eta_CHP_th"] #0.65
        # Power-to-heat ratio [-]
        cb_CHP = tech_infos["cb_CHP"] #0.3
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_CHP= tech_infos["Percentage_ramp_up_down_CHP"] #1
        ramp_up_down_CHP = Percentage_ramp_up_down_CHP*P_max_CHP
        # Variable cost [EUR/MWh]
        cost_var_CHP = tech_infos["cost_var_CHP"] #3.9
        # Fuel cost [EUR/MWh]
        cost_fuel_CHP = tech_infos["cost_fuel_CHP"] #20
        # Price paid to CHP for electrcity production
        CHP_electricity_price=Electricity_price

        ## Cogeneration 2 ##

        # Installed capacity [MW]
        P_max_CHP_2 = tech_infos["P_max_CHP_2"] #0
        # Technical minimum [MW]
        Technical_minimum_CHP_2= tech_infos["Technical_minimum_CHP_2"] #0.4
        P_min_CHP_2 = Technical_minimum_CHP_2*P_max_CHP_2
        # Efficiency [-]
        eta_CHP_th_2 = tech_infos["eta_CHP_th_2"] #0.65
        # Power-to-heat ratio [-]
        cb_CHP_2 = tech_infos["cb_CHP_2"] #0.3
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_CHP_2= tech_infos["Percentage_ramp_up_down_CHP_2"] #0.2
        ramp_up_down_CHP_2 = Percentage_ramp_up_down_CHP_2*P_max_CHP_2
        # Variable cost [EUR/MWh]
        cost_var_CHP_2 = tech_infos["cost_var_CHP_2"] #3.9
        # Fuel cost [EUR/MWh]
        cost_fuel_CHP_2 = tech_infos["cost_fuel_CHP_2"] #15
        # Price paid to CHP for electrcity production
        CHP_electricity_price_2=Electricity_price

        ## Absorption Heat pump ##

        # Installed capacity [MW]
        P_max_HP_abs = tech_infos["P_max_HP_abs"] #9
        # Technical minimum [MW]
        Technical_minimum_HP_abs= tech_infos["Technical_minimum_HP_abs"] #0.1
        P_min_HP_abs = Technical_minimum_HP_abs*P_max_HP_abs
        # Efficiency [-]
        eta_HP_abs = tech_infos["eta_HP_abs"] #0.9
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_abs= tech_infos["Percentage_ramp_up_down_HP_abs"] #1
        ramp_up_down_HP_abs = Percentage_ramp_up_down_HP_abs*P_max_HP_abs
        # Variable cost [EUR/MWh]
        cost_var_HP_abs = tech_infos["cost_var_HP_abs"] #0.1

        ## Waste heat - heat exchangers ##

        # Installed capacity [MW]
        P_HEX_capacity = tech_infos["P_HEX_capacity"] #3
        P_max_HEX=P_HEX_capacity*Heat_exchanger_specific
        # Technical minimum [MW]
        P_min_HEX= tech_infos["P_min_HEX"] #0
        # Variable cost [EUR/MWh]
        cost_var_HEX = tech_infos["cost_var_HEX"] #0.25

        ## Waste cooling - heat exchangers ##

        # Installed capacity [MW]
        P_HEX_capacity_cool = tech_infos["P_HEX_capacity_cool"] #2
        P_max_HEX_cool=P_HEX_capacity_cool*Heat_exchanger_specific_cool
        # Technical minimum [MW]
        P_min_HEX_cool= tech_infos["P_min_HEX_cool"] #0
        # Variable cost [EUR/MWh]
        cost_var_HEX_cool = tech_infos["cost_var_HEX_cool"] #0.25

        ## Waste heat - heat pumps - temperature group 1 ##

        # Installed capacity [MW]
        P_max_HP_waste_heat_1 = tech_infos["P_max_HP_waste_heat_1"] #1
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_1 = tech_infos["Technical_minimum_HP_waste_heat_1"] #0.0
        P_min_HP_waste_heat_1 = Technical_minimum_HP_waste_heat_1*P_max_HP_waste_heat_1
        # Efficiency [-]
        eta_HP_waste_heat_1 = eta_HP_waste_heat_1
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_1= tech_infos["Percentage_ramp_up_down_HP_waste_heat_1"] #0.95
        ramp_up_down_HP_waste_heat_1 = Percentage_ramp_up_down_HP_waste_heat_1*P_max_HP_waste_heat_1
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_1 = tech_infos["cost_var_HP_waste_heat_1"] #0.7
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_waste_heat_1 = Electricity_price

        ## Waste heat - heat pumps - temperature group 2 ##

        # Installed capacity [MW]
        P_max_HP_waste_heat_2 = tech_infos["P_max_HP_waste_heat_2"] #1
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_2 = tech_infos["Technical_minimum_HP_waste_heat_2"] #0.1
        P_min_HP_waste_heat_2 = Technical_minimum_HP_waste_heat_2*P_max_HP_waste_heat_2
        # Efficiency [-]
        eta_HP_waste_heat_2 = eta_HP_waste_heat_2
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_2= tech_infos["Percentage_ramp_up_down_HP_waste_heat_2"] #0.95
        ramp_up_down_HP_waste_heat_2 = Percentage_ramp_up_down_HP_waste_heat_2*P_max_HP_waste_heat_2
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_2 = tech_infos["cost_var_HP_waste_heat_2"] #1.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_waste_heat_2 = Electricity_price

        ## Waste heat - heat pumps - temperature group 3 ##

        # Installed capacity [MW]
        P_max_HP_waste_heat_3 = tech_infos["P_max_HP_waste_heat_3"] #1
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_3 = tech_infos["Technical_minimum_HP_waste_heat_3"] #0.1
        P_min_HP_waste_heat_3 = Technical_minimum_HP_waste_heat_3*P_max_HP_waste_heat_3
        # Efficiency [-]
        eta_HP_waste_heat_3 = eta_HP_waste_heat_3
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_3= tech_infos["Percentage_ramp_up_down_HP_waste_heat_3"] #0.95
        ramp_up_down_HP_waste_heat_3 = Percentage_ramp_up_down_HP_waste_heat_3*P_max_HP_waste_heat_3
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_3 = tech_infos["cost_var_HP_waste_heat_3"] #1.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_waste_heat_3 = Electricity_price

        ## Thermal energy storage ##

        # Thermal energy storage size [MWh]
        TES_size = tech_infos["TES_size"] #200
        # Minimum percentage of state of charge of thermal energy storage [-]
        SOC_min = tech_infos["SOC_min"] #0.05
        # Thermal energy storage state of charge at the start/end of the year [-]
        TES_start_end = tech_infos["TES_start_end"] #0.05
        # Thermal storage charge/discharge capacity [MW]
        TES_charge_discharge_time = tech_infos["TES_charge_discharge_time"] #5
        TES_in_out_max = TES_size/TES_charge_discharge_time
        TES_loss = tech_infos["TES_loss"] #0.1/(24*7)

        ######## Constraints and variables ########

        ### Variables ###

        # Technologies' operation
        @variable(m, P_max_HOB >= HOB[1:length(DEM)] >= 0)
        @variable(m, P_max_HOB_2 >= HOB_2[1:length(DEM)] >= 0)
        @variable(m, P_max_HOB_3 >= HOB_3[1:length(DEM)] >= 0)

        @variable(m, P_max_HP >= HP[1:length(DEM)] >= 0)
        @variable(m, P_max_HP_2 >= HP_2[1:length(DEM)] >= 0)
        @variable(m, P_max_HP_3 >= HP_3[1:length(DEM)] >= 0)

        @variable(m, P_max_CHP >= CHP[1:length(DEM)] >= 0)
        @variable(m, P_max_CHP_2 >= CHP_2[1:length(DEM)] >= 0)

        @variable(m, P_max_HP_abs >= HP_abs[1:length(DEM)] >= 0)

        @variable(m, HEX[1:length(DEM)] >= 0)

        @variable(m, HEX_cool[1:length(DEM)] >= 0)

        @variable(m, P_max_HP_waste_heat_1 >= HP_waste_heat_1[1:length(DEM)] >= 0)
        @variable(m, P_max_HP_waste_heat_2 >= HP_waste_heat_2[1:length(DEM)] >= 0)
        @variable(m, P_max_HP_waste_heat_3 >= HP_waste_heat_3[1:length(DEM)] >= 0)

        @variable(m, TES_in_out_max >=TES[1:length(DEM)]>= -TES_in_out_max)
        @variable(m, TES_size >= SOC[1:length(DEM)] >= TES_size*SOC_min)

        #@variable(m, maximum(DEM) >= Error[1:length(DEM)] >= 0)
        @variable(m, 0 >= Error[1:length(DEM)] >= 0)

        # Binary variables, on/off

        @variable(m, OnOffHOB[1:length(DEM)], Bin)
        @variable(m, OnOffHOB_2[1:length(DEM)], Bin)
        @variable(m, OnOffHOB_3[1:length(DEM)], Bin)

        @variable(m, OnOffHP[1:length(DEM)], Bin)
        @variable(m, OnOffHP_2[1:length(DEM)], Bin)
        @variable(m, OnOffHP_3[1:length(DEM)], Bin)

        @variable(m, OnOffHP_abs[1:length(DEM)], Bin)

        @variable(m, OnOffCHP[1:length(DEM)], Bin)
        @variable(m, OnOffCHP_2[1:length(DEM)], Bin)

        @variable(m, OnOffHP_waste_heat_1[1:length(DEM)], Bin)
        @variable(m, OnOffHP_waste_heat_2[1:length(DEM)], Bin)
        @variable(m, OnOffHP_waste_heat_3[1:length(DEM)], Bin)

        @variable(m, OnOffHP_waste_heat_seasonal[1:length(DEM)], Bin)

        ### Constraints ###

        # Supply + TES = DEM
        @constraint(m, HP+HP_2+HP_3+HEX_cool+HP_abs+HP_waste_heat_1+HP_waste_heat_2+HP_waste_heat_3 - TES + Error .== DEM)

        @constraint(m, HP_abs.==(HOB+HOB_2+HOB_3+CHP+CHP_2+HEX)*eta_HP_abs)

        # Thermal storage constraints
        t=2:(length(DEM)-1)
        @constraint(m, SOC[t]-SOC[t .- 1]-TES[t]+SOC[t]*TES_loss .== 0)
        @constraint(m, SOC[length(DEM)-1]+TES[length(DEM)] == TES_size*TES_start_end)
        @constraint(m, SOC[1]-TES[1] == TES_size*TES_start_end)
        @constraint(m, SOC[length(DEM)] == TES_size*TES_start_end)

        # Waste heat constraints - can't utilize more than available
        @constraint(m, HEX .<=P_max_HEX)
        @constraint(m, HEX_cool .<=P_max_HEX_cool)

        @constraint(m, HP_waste_heat_1 .<= Available_waste_heat_heat_pump_1)
        @constraint(m, HP_waste_heat_2 .<= Available_waste_heat_heat_pump_2)
        @constraint(m, HP_waste_heat_3 .<= Available_waste_heat_heat_pump_3)

        # OnOff constraints - technical minimum
        @constraint(m, OnOffHOB*P_min_HOB-HOB.<=0)
        @constraint(m, -OnOffHOB*P_max_HOB+HOB.<=0)
        @constraint(m, OnOffHOB_2*P_min_HOB_2-HOB_2.<=0)
        @constraint(m, -OnOffHOB_2*P_max_HOB_2+HOB_2.<=0)
        @constraint(m, OnOffHOB_3*P_min_HOB_3-HOB_3.<=0)
        @constraint(m, -OnOffHOB_3*P_max_HOB_3+HOB_3.<=0)

        @constraint(m, OnOffHP*P_min_HP-HP.<=0)
        @constraint(m, -OnOffHP*P_max_HP+HP.<=0)
        @constraint(m, OnOffHP_2*P_min_HP_2-HP_2.<=0)
        @constraint(m, -OnOffHP_2*P_max_HP_2+HP_2.<=0)
        @constraint(m, OnOffHP_3*P_min_HP_3-HP_3.<=0)
        @constraint(m, -OnOffHP_3*P_max_HP_3+HP_3.<=0)

        @constraint(m, OnOffHP_abs*P_min_HP_abs-HP_abs.<=0)
        @constraint(m, -OnOffHP_abs*P_max_HP_abs+HP_abs.<=0)

        @constraint(m, OnOffCHP*P_min_CHP-CHP.<=0)
        @constraint(m, -OnOffCHP*P_max_CHP+CHP.<=0)
        @constraint(m, OnOffCHP_2*P_min_CHP_2-CHP_2.<=0)
        @constraint(m, -OnOffCHP_2*P_max_CHP_2+CHP_2.<=0)

        @constraint(m, OnOffHP_waste_heat_1*P_min_HP_waste_heat_1-HP_waste_heat_1.<=0)
        @constraint(m, -OnOffHP_waste_heat_1*P_max_HP_waste_heat_1+HP_waste_heat_1.<=0)
        @constraint(m, OnOffHP_waste_heat_2*P_min_HP_waste_heat_2-HP_waste_heat_2.<=0)
        @constraint(m, -OnOffHP_waste_heat_2*P_max_HP_waste_heat_2+HP_waste_heat_2.<=0)
        @constraint(m, OnOffHP_waste_heat_3*P_min_HP_waste_heat_3-HP_waste_heat_3.<=0)
        @constraint(m, -OnOffHP_waste_heat_3*P_max_HP_waste_heat_3+HP_waste_heat_3.<=0)

        # Ramp-up-down constraint for each technology

        i=2:(length(DEM)-1)

        @constraint(m, HOB[i .- 1]-HOB[i].<=ramp_up_down_HOB*OnOffHOB[i]+(1 .- OnOffHOB[i])*P_min_HOB)
        @constraint(m, HOB[i]-HOB[i .- 1].<=ramp_up_down_HOB*OnOffHOB[i .- 1]+(1 .- OnOffHOB[i .- 1])*P_min_HOB)
        @constraint(m, HOB_2[i .- 1]-HOB_2[i].<=ramp_up_down_HOB_2*OnOffHOB_2[i]+(1 .- OnOffHOB_2[i])*P_min_HOB_2)
        @constraint(m, HOB_2[i]-HOB_2[i .- 1].<=ramp_up_down_HOB_2*OnOffHOB_2[i .- 1]+(1 .- OnOffHOB_2[i .- 1])*P_min_HOB_2)
        @constraint(m, HOB_3[i .- 1]-HOB_3[i].<=ramp_up_down_HOB_3*OnOffHOB_3[i]+(1 .- OnOffHOB_3[i])*P_min_HOB_3)
        @constraint(m, HOB_3[i]-HOB_3[i .- 1].<=ramp_up_down_HOB_3*OnOffHOB_3[i .- 1]+(1 .- OnOffHOB_3[i .- 1])*P_min_HOB_3)

        @constraint(m, HP[i .- 1]-HP[i].<=ramp_up_down_HP*OnOffHP[i]+(1 .- OnOffHP[i])*P_min_HP)
        @constraint(m, HP[i]-HP[i .- 1].<=ramp_up_down_HP*OnOffHP[i .- 1]+(1 .- OnOffHP[i .- 1])*P_min_HP)
        @constraint(m, HP_2[i .- 1]-HP_2[i].<=ramp_up_down_HP_2*OnOffHP_2[i]+(1 .- OnOffHP_2[i])*P_min_HP_2)
        @constraint(m, HP_2[i]-HP_2[i .- 1].<=ramp_up_down_HP_2*OnOffHP_2[i .- 1]+(1 .- OnOffHP_2[i .- 1])*P_min_HP_2)
        @constraint(m, HP_3[i .- 1]-HP_3[i].<=ramp_up_down_HP_3*OnOffHP_3[i]+(1 .- OnOffHP_3[i])*P_min_HP_3)
        @constraint(m, HP_3[i]-HP_3[i .- 1].<=ramp_up_down_HP_3*OnOffHP_3[i .- 1]+(1 .- OnOffHP_3[i .- 1])*P_min_HP_3)

        @constraint(m, HP_abs[i .- 1]-HP_abs[i].<=ramp_up_down_HP_abs*OnOffHP_abs[i]+(1 .- OnOffHP_abs[i])*P_min_HP_abs)
        @constraint(m, HP_abs[i]-HP_abs[i .- 1].<=ramp_up_down_HP_abs*OnOffHP_abs[i .- 1]+(1 .- OnOffHP_abs[i .- 1])*P_min_HP_abs)

        @constraint(m, CHP[i .- 1]-CHP[i].<=ramp_up_down_CHP*OnOffCHP[i]+(1 .- OnOffCHP[i])*P_min_CHP)
        @constraint(m, CHP[i]-CHP[i .- 1].<=ramp_up_down_CHP*OnOffCHP[i .- 1]+(1 .- OnOffCHP[i .- 1])*P_min_CHP)
        @constraint(m, CHP_2[i .- 1]-CHP_2[i].<=ramp_up_down_CHP_2*OnOffCHP_2[i]+(1 .- OnOffCHP_2[i])*P_min_CHP_2)
        @constraint(m, CHP_2[i]-CHP_2[i .- 1].<=ramp_up_down_CHP_2*OnOffCHP_2[i .- 1]+(1 .- OnOffCHP_2[i .- 1])*P_min_CHP_2)

        @constraint(m, HP_waste_heat_1[i .- 1]-HP_waste_heat_1[i].<=ramp_up_down_HP_waste_heat_1*OnOffHP_waste_heat_1[i]+(1 .- OnOffHP_waste_heat_1[i])*P_min_HP_waste_heat_1)
        @constraint(m, HP_waste_heat_1[i]-HP_waste_heat_1[i .- 1].<=ramp_up_down_HP_waste_heat_1*OnOffHP_waste_heat_1[i .- 1]+(1 .- OnOffHP_waste_heat_1[i .- 1])*P_min_HP_waste_heat_1)
        @constraint(m, HP_waste_heat_2[i .- 1]-HP_waste_heat_1[i].<=ramp_up_down_HP_waste_heat_2*OnOffHP_waste_heat_2[i]+(1 .- OnOffHP_waste_heat_2[i])*P_min_HP_waste_heat_2)
        @constraint(m, HP_waste_heat_2[i]-HP_waste_heat_2[i .- 1].<=ramp_up_down_HP_waste_heat_2*OnOffHP_waste_heat_2[i .- 1]+(1 .- OnOffHP_waste_heat_2[i .- 1])*P_min_HP_waste_heat_2)
        @constraint(m, HP_waste_heat_3[i .- 1]-HP_waste_heat_3[i].<=ramp_up_down_HP_waste_heat_3*OnOffHP_waste_heat_3[i]+(1 .- OnOffHP_waste_heat_3[i])*P_min_HP_waste_heat_3)
        @constraint(m, HP_waste_heat_3[i]-HP_waste_heat_3[i .- 1].<=ramp_up_down_HP_waste_heat_3*OnOffHP_waste_heat_3[i .- 1]+(1 .- OnOffHP_waste_heat_3[i .- 1])*P_min_HP_waste_heat_3)

        ###### Objective function ######

        # Minimize total running costs

        @objective(m,Min, sum(Error)*1e10
                          + sum(HOB)*(cost_var_HOB .+ cost_fuel_HOB/eta_HOB)
                          + sum(HOB_2)*(cost_var_HOB_2 .+ cost_fuel_HOB_2/eta_HOB_2)
                          + sum(HOB_3)*(cost_var_HOB_3 .+ cost_fuel_HOB_3/eta_HOB_3)
                          + sum(HP.*(cost_var_HP .+ cost_fuel_HP./eta_HP))
                          + sum(HP_2.*(cost_var_HP_2 .+ cost_fuel_HP_2./eta_HP_2))
                          + sum(HP_3.*(cost_var_HP_3 .+ cost_fuel_HP_3./eta_HP_3))
                          + sum(HP_abs)*cost_var_HP_abs
                          + sum(CHP)*(cost_var_CHP .+ cost_fuel_CHP/eta_CHP_th)
                          + sum(CHP_2)*(cost_var_CHP_2 .+ cost_fuel_CHP_2/eta_CHP_th_2)
                          + sum(HEX)*cost_var_HEX+sum(HEX_cool)*cost_var_HEX_cool
                          + sum(HP_waste_heat_1.*(cost_var_HP_waste_heat_1 .+ cost_fuel_HP_waste_heat_1./eta_HP_waste_heat_1))
                          + sum(HP_waste_heat_2.*(cost_var_HP_waste_heat_2 .+ cost_fuel_HP_waste_heat_2./eta_HP_waste_heat_2))
                          + sum(HP_waste_heat_3.*(cost_var_HP_waste_heat_3 .+ cost_fuel_HP_waste_heat_3./eta_HP_waste_heat_3))
                          - sum(CHP*cb_CHP.*CHP_electricity_price)
                          - sum(CHP_2*cb_CHP_2.*CHP_electricity_price_2))

        optimize!(m)
        status = termination_status(m)
        if primal_status(m) == MOI.NO_SOLUTION
            print_function("No solution found.")
            return
        end

        ##### Result visualization ####

        Result_HOB=value.(HOB)
        Result_HOB_2=value.(HOB_2)
        Result_HOB_3=value.(HOB_3)

        Result_HP=value.(HP)
        Result_HP_2=value.(HP_2)
        Result_HP_3=value.(HP_3)

        Result_HP_abs=value.(HP_abs)

        Result_CHP=value.(CHP)
        Result_CHP_2=value.(CHP_2)

        Result_HEX=value.(HEX)
        Result_HEX_cool=value.(HEX_cool)

        Result_HP_waste_heat_1=value.(HP_waste_heat_1)
        Result_HP_waste_heat_2=value.(HP_waste_heat_2)
        Result_HP_waste_heat_3=value.(HP_waste_heat_3)

        Result_SOC=value.(SOC)

        Result_Error=value.(Error)

        ## Export results ##

        writedlm(string(output_folder, "\\Result_HOB_", network_id, ".csv"), Result_HOB, csv_sep)
        writedlm(string(output_folder, "\\Result_HOB_2_", network_id, ".csv"), Result_HOB_2, csv_sep)
        writedlm(string(output_folder, "\\Result_HOB_3_", network_id, ".csv"), Result_HOB_3, csv_sep)

        writedlm(string(output_folder, "\\Result_HP_", network_id, ".csv"), Result_HP, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_2_", network_id, ".csv"), Result_HP_2, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_3_", network_id, ".csv"), Result_HP_3, csv_sep)

        writedlm(string(output_folder, "\\Result_HP_abs_", network_id, ".csv"), Result_HP_abs, csv_sep)

        writedlm(string(output_folder, "\\Result_CHP_", network_id, ".csv"), Result_CHP, csv_sep)
        writedlm(string(output_folder, "\\Result_CHP_2_", network_id, ".csv"), Result_CHP_2, csv_sep)

        writedlm(string(output_folder , "\\Result_HEX_", network_id, ".csv"), Result_HEX, csv_sep)

        writedlm(string(output_folder, "\\Result_HEX_cool_", network_id, ".csv"), Result_HEX_cool, csv_sep)

        writedlm(string(output_folder, "\\Result_HP_waste_heat_1_", network_id, ".csv"), Result_HP_waste_heat_1, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_waste_heat_2_", network_id, ".csv"), Result_HP_waste_heat_2, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_waste_heat_3_", network_id, ".csv"), Result_HP_waste_heat_3, csv_sep)

        writedlm(string(output_folder, "\\Result_DEM_", network_id, ".csv"), DEM, csv_sep)

        writedlm(string(output_folder, "\\Result_SOC_", network_id, ".csv"), Result_SOC, csv_sep)

        writedlm(string(output_folder, "\\Result_Error_", network_id, ".csv"), Result_Error, csv_sep)

    end
    export district_cooling_simulation
end

max_stars_repo_name: Planheat/Planheat-Tool

max_issues_repo_name: Planheat/Planheat-Tool

max_forks_repo_name: Planheat/Planheat-Tool


text:
"""
Model to approximate cross products of node voltages

```
wdcr[(i,j)] <= wdc[i]*wdc[j]
```
"""
function constraint_voltage_dc(pm::AbstractWRModel, n::Int, cnd::Int)
    wdc = PowerModels.var(pm, n, cnd, :wdc)
    wdcr = PowerModels.var(pm, n, cnd, :wdcr)

    for (i,j) in PowerModels.ids(pm, n, :buspairsdc)
        JuMP.@constraint(pm.model, wdcr[(i,j)]^2 <= wdc[i]*wdc[j])
        # InfrastructureModels.relaxation_product(pm.model, wdc[i], wdc[j], wdcr_, wdcr_)
    end
end

function constraint_voltage_dc(pm::AbstractWRConicModel, n::Int, cnd::Int)
    wdc = PowerModels.var(pm, n, cnd, :wdc)
    wdcr = PowerModels.var(pm, n, cnd, :wdcr)

    for (i,j) in PowerModels.ids(pm, n, :buspairsdc)
        relaxation_complex_product_conic(pm.model, wdc[i], wdc[j], wdcr[(i,j)])
    end
end

"""
Limits dc branch current

```
p[f_idx] <= wdc[f_bus] * Imax
```
"""
function constraint_dc_branch_current(pm::AbstractWRModel, n::Int, cnd::Int, f_bus, f_idx, ccm_max, p)
    p_dc_fr = PowerModels.var(pm, n, cnd, :p_dcgrid, f_idx)
    wdc_fr = PowerModels.var(pm, n, cnd, :wdc, f_bus)

    @constraint(pm.model, p_dc_fr <= wdc_fr * ccm_max * p^2)
end


text:
export build_wire_coordinates, get_minimum_intersection_distance

function build_wire_coordinates(wire_directions)
    path = []
    current_x = 0
    current_y = 0
    dx = Dict('L'=> -1, 'R'=> 1, 'U'=> 0, 'D'=> 0)
    dy = Dict('L'=> 0, 'R'=> 0, 'U'=> 1, 'D'=> -1)
    for direction_instruction in wire_directions
        direction = direction_instruction[1]
        count = parse(Int, direction_instruction[2:end])
        step_x = dx[direction]
        step_y = dy[direction]
        for _ = 1:count
            current_x += step_x
            current_y += step_y
            push!(path, (current_x, current_y))
        end
    end
    return path
end

function get_minimum_intersection_distance(intersections)
    minimum_distance = nothing
    for row in eachrow(intersections)
        distance = sum(abs.(row[1]))
        if isnothing(minimum_distance) || minimum_distance > distance
            minimum_distance = distance
        end
    end
    return minimum_distance
end

d03_inputpath = joinpath(@__DIR__, "input.txt")
d03_inputread(path) = readdlm(path, ',')

function d03_answer1(input)
    coords1 = build_wire_coordinates(input[1, :])
    coords2 = build_wire_coordinates(input[2, :])
    intersections = intersect(coords1, coords2)
    return get_minimum_intersection_distance(intersections)
end

function d03_answer2(input)
    return 0
end


text:
module ConjugateGradientTests

using Test, LinearAlgebra, LazyAlgebra, ConjugateGradient

VERBOSE = true
io = (VERBOSE ? stdout : ConjugateGradient.quiet)
dims = (29,11)
nrows, ncols = dims
@testset "least-square fit ($T)" for T in (Float32, Float64)
    # State the problem (least square fit).
    H = rand(T, dims) # direct model
    x = rand(T, ncols) # ground truth
    y = H*x + T(0.01)*rand(T, nrows) # noisy data

    # LHS matrix and RHS vector of the normal equations.
    A = H'*H
    b = H'*y

    # Least-squares solution using Julia linear algebra.
    x0 = A\b

    # Tolerances.
    rtol = 1e-3
    atol = rtol*vnorm2(x0)

    # Solve with conjugate gradient algorithm.
    n = length(b)
    ws = ConjugateGradient.Context(b;
                                   precond = true,
                                   maxiter = n+1,
                                   ftol = (0.0, 0.0),
                                   gtol = (0.0, 0.0),
                                   xtol = (0.0, 0.0))
    x1 = fill!(similar(b), 0)
    VERBOSE && println("\n# Run the conjugate gradient from 0:")
    status = ConjugateGradient.solve!(x1, A, b, ws, io)
    @test status ≥ 0
    VERBOSE && println("rel. err. for x1: ", vnorm2(x1 - x0)/vnorm2(x0))
    @test vnorm2(x1 - x0) ≤ atol

    # Simple diagonal preconditioner.
    q = convert(Vector{T}, diag(A))
    for i in eachindex(q)
        q[i] = 1/q[i]
    end
    M = Diag(q)

    x2 = fill!(similar(b), 0)
    VERBOSE && println("\n# Run the preconditioned conjugate gradient from 0:")
    status = ConjugateGradient.solve!(x2, A, b, M, ws, io)
    @test status ≥ 0
    VERBOSE && println("rel. err. for x1: ", vnorm2(x1 - x0)/vnorm2(x0))
    @test vnorm2(x1 - x0) ≤ atol
end

end # module


text:
using TemporalGPs: posterior_and_lml, predict, predict_marginals

@testset "linear_gaussian_conditionals" begin
    Dlats = [1, 3]
    Dobss = [1, 2]
    # Dlats = [3]
    # Dobss = [2]
    storages = [
        (name="dense storage Float64", val=ArrayStorage(Float64)),
    ]
    Q_types = [
        Val(:dense),
        Val(:diag),
    ]

    @testset "SmallOutputLGC (Dlat=$Dlat, Dobs=$Dobs, Q=$(Q_type), $(storage.name))" for
        Dlat in Dlats,
        Dobs in Dobss,
        Q_type in Q_types,
        storage in storages

        println("SmallOutputLGC (Dlat=$Dlat, Dobs=$Dobs, Q=$(Q_type), $(storage.name))")

        rng = MersenneTwister(123456)
        x = random_gaussian(rng, Dlat, storage.val)
        model = random_small_output_lgc(rng, Dlat, Dobs, Q_type, storage.val)

        test_interface(
            rng, model, x;
            check_adjoints=true,
            check_infers=true,
            check_allocs=storage.val isa SArrayStorage,
        )

        Q_type == Val(:diag) && @testset "missing data" begin

            # Generate observation with a missing.
            y = conditional_rand(rng, model, rand(rng, x))
            y_missing = Vector{Union{Missing, eltype(y)}}(undef, length(y))
            y_missing .= y
            y_missing[1] = missing

            # Compute logpdf and posterior under model given the missing.
            x_post, lml = posterior_and_lml(x, model, y_missing)

            # Modify the model and compute the logpdf and posterior.
            new_model = SmallOutputLGC(
                model.A[2:end, :], model.a[2:end], model.Q[2:end, 2:end],
            )
            y_new = y[2:end]
            x_post_new, lml_new = posterior_and_lml(x, new_model, y_new)

            # Verify that both things give the same answer.
            @test x_post ≈ x_post_new
            @test lml ≈ lml_new atol=1e-8 rtol=1e-8

            # Check that everything infers and AD gives the right answer.
            @inferred posterior_and_lml(x, model, y_missing)
            x̄ = adjoint_test(posterior_and_lml, (x, model, y_missing))
            @test x̄[2].Q isa NamedTuple{(:diag, )}
        end
    end

    @testset "LargeOutputLGC (Dlat=$Dlat, Dobs=$Dobs, Q=$(Q_type), $(storage.name))" for
        Dlat in Dlats,
        Dobs in Dobss,
        Q_type in Q_types,
        storage in storages

        println("LargeOutputLGC (Dlat=$Dlat, Dobs=$Dobs, Q=$(Q_type), $(storage.name))")

        rng = MersenneTwister(123456)
        x = random_gaussian(rng, Dlat, storage.val)
        model = random_large_output_lgc(rng, Dlat, Dobs, Q_type, storage.val)

        @testset "consistency with SmallOutputLGC" begin
            y = rand(rng, TemporalGPs.predict(x, model))
            vanilla_model = TemporalGPs.SmallOutputLGC(model.A, model.a, model.Q)
            x_vanilla, lml_vanilla = posterior_and_lml(x, vanilla_model, y)
            x_large, lml_large = posterior_and_lml(x, model, y)

            @test x_vanilla.m ≈ x_large.m
            @test x_vanilla.P ≈ x_large.P
            @test lml_vanilla ≈ lml_large

            @test predict(x, vanilla_model) ≈ predict(x, model)

            Q_type == Val(:diag) && @testset "missing data" begin

                # Create missing data.
                y_missing = Vector{Union{Missing, eltype(y)}}(undef, length(y))
                y_missing .= y
                y_missing[1] = missing

                # Compute posterior and lml under both SmallOutputLGC and LargeOutputLGC.
                x_post_vanilla, lml_vanilla = posterior_and_lml(x, vanilla_model, y_missing)
                x_post_large, lml_large = posterior_and_lml(x, model, y_missing)

                # Check that they give roughly the same answer.
                @test x_post_vanilla ≈ x_post_large
                @test lml_vanilla ≈ lml_large rtol=1e-8 atol=1e-8

                # Check that everything infers and AD gives the right answer.
                @inferred posterior_and_lml(x, model, y_missing)
                x̄ = adjoint_test(posterior_and_lml, (x, model, y_missing))
                @test x̄[2].Q isa NamedTuple{(:diag, )}
            end
        end

        test_interface(
            rng, model, x;
            check_adjoints=true,
            check_infers=true,
            check_allocs=storage.val isa SArrayStorage,
        )
    end

    @testset "ScalarOutputLGC (Dlat=$Dlat, ($storage.name))" for
        Dlat in Dlats,
        storage in [
            (name="dense storage Float64", val=ArrayStorage(Float64)),
            (name="static storage Float64", val=SArrayStorage(Float64)),
        ]

        println("ScalarOutputLGC (Dlat=$Dlat, ($storage.name))")

        rng = MersenneTwister(123456)
        x = random_gaussian(rng, Dlat, storage.val)
        model = random_scalar_output_lgc(rng, Dlat, storage.val)

        @testset "consistency with LGC" begin
            vanilla_model = lgc_from_scalar_output_lgc(model)
            y_vanilla = rand(rng, predict(x, vanilla_model))
            x_vanilla, lml_vanilla = posterior_and_lml(x, vanilla_model, y_vanilla)
            x_scalar, lml_scalar = posterior_and_lml(x, model, only(y_vanilla))

            @test x_vanilla.m ≈ x_scalar.m
            @test x_vanilla.P ≈ x_scalar.P
            @test lml_vanilla ≈ lml_scalar
        end

        test_interface(
            rng, model, x;
            check_adjoints=true,
            check_infers=true,
            check_allocs=storage.val isa SArrayStorage,
        )
    end

    Dmids = [1, 3]

    @testset "BottleneckLGC (Din=$Din, Dmid=$Dmid, Dout=$Dout, Q=$(Q_type))" for
        Din in Dlats,
        Dout in Dobss,
        Dmid in Dmids,
        Q_type in Q_types

        println("BottleneckLGC (Din=$Din, Dmid=$Dmid, Dout=$Dout, Q=$(Q_type))")

        storage = ArrayStorage(Float64)
        rng = MersenneTwister(123456)
        x = random_gaussian(rng, Din, storage)
        model = random_bottleneck_lgc(rng, Din, Dmid, Dout, Q_type, storage)

        @test TemporalGPs.dim_out(model) == Dout
        @test TemporalGPs.dim_in(model) == Din

        test_interface(
            rng, model, x;
            check_adjoints=true,
            check_infers=true,
            check_allocs=false,
        )

        @testset "consistency with SmallOutputLGC" begin
            vanilla_model = small_output_lgc_from_bottleneck(model)

            @test predict(x, vanilla_model) ≈ predict(x, model)
            @test predict_marginals(x, vanilla_model) ≈ predict_marginals(x, model)

            y = rand(rng, predict(x, model))
            x_vanilla, lml_vanilla = posterior_and_lml(x, vanilla_model, y)
            x_bottle, lml_bottle = posterior_and_lml(x, model, y)
            @test x_vanilla.P ≈ x_bottle.P rtol=1e-6
            @test x_vanilla.m ≈ x_bottle.m rtol=1e-6
            @test lml_vanilla ≈ lml_bottle rtol=1e-6

            Q_type == Val(:diag) && @testset "missing data" begin

                # Create missing data.
                y_missing = Vector{Union{Missing, eltype(y)}}(undef, length(y))
                y_missing .= y
                y_missing[1] = missing

                # Compute posterior and lml under both SmallOutputLGC and LargeOutputLGC.
                x_post_vanilla, lml_vanilla = posterior_and_lml(x, vanilla_model, y_missing)
                x_post_large, lml_large = posterior_and_lml(x, model, y_missing)

                # Check that they give roughly the same answer.
                @test x_post_vanilla ≈ x_post_large rtol=1e-8 atol=1e-8
                @test lml_vanilla ≈ lml_large rtol=1e-8 atol=1e-8

                # Check that everything infers and AD gives the right answer.
                @inferred posterior_and_lml(x, model, y_missing)
                x̄ = adjoint_test(posterior_and_lml, (x, model, y_missing))
                @test x̄[2].fan_out.Q isa NamedTuple{(:diag, )}
            end
        end
    end
end


text:
struct Arc
    u::Any
    v::Any
    key::Any
end

Arc(t::Tuple) = Arc(t[1], t[2], t[3])

source(a::Arc) = a.u
target(a::Arc) = a.v
key(a::Arc) = a.key

add_arc!(g::Graph, a::Arc) = add_arc!(g, a.u, a.v, a.key)
rem_arc!(g::Graph, a::Arc) = rem_arc!(g, a.u, a.v, a.key)

function add_arcs_from!(g::Graph, abunch::Vector{Arc})
    keylist = []
    for a in abunch
        res = add_arc!(g, a)
        push!(keylist, res)
    end
    return keylist
end

function rem_arcs_from!(g::Graph, abunch::Vector{Arc})
    removed = 0
    for a in abunch
        remove = try
            rem_arc!(g, a) # if a doesn't comply,
        catch end # just pass.
        if remove
            removed += 1
        end
    end
    return removed
end

function ==(a1::Arc, a2::Arc)
    return source(a1) == source(a2) && target(a1) == target(a2) && key(a1) == key(a2)
end

getindex(g::Graph, a::Arc) = g.fadj[a.u][a.v][a.key]
getindex(g::Graph, u, v, key) = g[Arc(u, v, key)]
getindex(g::Graph, a::Arc, attrib::Symbol) = g[a][attrib]
getindex(g::Graph, u, v, key, attrib::Symbol) = g[Arc(u, v, key)][attrib]

"""
    setindex!(g, val, a, attrib)
setindex! are exclusive for setting arc attributes only.
To delete an attribute, use delete!(g[a], attrib).
"""
function setindex!(g::Graph, val, a::Arc, attrib::Symbol)
    if has_arc(g, a)
        set!(g.fadj[a.u][a.v][a.key], attrib, val)
        set!(g.badj[a.v][a.u][a.key], attrib, val)
        return g[a,attrib]
    else
        return nothing
    end
end

setindex!(g::Graph, val, u, v, key, attrib::Symbol) = setindex!(
    g, val, Arc(u, v, key), attrib
)

"""
    has_arc(g, a)
has_arc checks wether elementary arc a (u,v,key) exists in graph.
"""
function has_arc(g::Graph, a::Arc)
    try
        g[a]
    catch KeyError
        return false
    end
    return true
end

"""
    ArcIter
An iterator of the arcs in the graph.
"""
struct ArcIter
    len::Integer
    ulist::Vector
    vlist::Vector
    keylist::Vector
end

function iterate(ait::ArcIter, state=1)
    if state <= ait.len
        res = Arc(ait.ulist[state], ait.vlist[state], ait.keylist[state])
        state += 1
        return res, state
    else
        return nothing
    end
end

eltype(::ArcIter) = Arc
length(ait::ArcIter) = ait.len
keys(ait::ArcIter) = 1:ait.len

function rand(ait::ArcIter)
    r = rand(1:ait.len)
    res = (ait.ulist[r], ait.vlist[r], ait.keylist[r])
    return res
end

function ==(ait1::ArcIter, ait2::ArcIter)
    return ait1.len == ait2.len && ait1.ulist == ait2.ulist && ait1.vlist ==
    ait2.vlist && ait1.keylist == ait2.keylist
end

"""
    ArcIter(g)
constructs an iterator of all arc in graph.
"""
function ArcIter(g::Graph, ::Colon, ::Colon)
    ulist = []; sizehint!(ulist, g.na)
    vlist = []; sizehint!(vlist, g.na)
    keylist = []; sizehint!(keylist, g.na)
    @inbounds for n in keys(g.fadj), nbrs in keys(g.fadj[n]), key in keys(g.fadj[n][nbrs])
        push!(ulist, n)
        push!(vlist, nbrs)
        push!(keylist, key)
    end
    return ArcIter(g.na, ulist, vlist, keylist)
end

"""
    ArcIter(g, u)
constructs an iterator of all arc starting from u (single or multiple node).
ArcIter is precise and will raise an error if an element in u is not in graph.
"""
function ArcIter(g::Graph, u, ::Colon)
    len = 0
    ulist = []
    vlist = []
    keylist = []
    @inbounds for n in u, nbrs in keys(g.fadj[n]), key in keys(g.fadj[n][nbrs])
        push!(ulist, n)
        push!(vlist, nbrs)
        push!(keylist, key)
        len += 1
    end
    return ArcIter(len, ulist, vlist, keylist)
end

"""
    ArcIter(g, :, v)
constructs an iterator of all arc ending at v.
ArcIter is precise and will raise an error if an element in v is not in graph.
"""
function ArcIter(g::Graph, ::Colon, v)
    len = 0
    ulist = []
    vlist = []
    keylist = []
    @inbounds for nbrs in v, n in keys(g.badj[nbrs]), key in keys(g.badj[nbrs][n])
        push!(ulist, n)
        push!(vlist, nbrs)
        push!(keylist, key)
        len += 1
    end
    return ArcIter(len, ulist, vlist, keylist)
end

"""
    ArcIter(g, u, v)
constructs an iterator of all arc inside u-v.
ArcIter is precise and will raise an error if an element in v is not in graph.
"""
function ArcIter(g::Graph, u, v)
    len = length(g.fadj[u][v])
    ulist = []; sizehint!(ulist, len)
    vlist = []; sizehint!(vlist, len)
    keylist = []; sizehint!(keylist, len)
    @inbounds for key in keys(g.fadj[u][v])
        push!(ulist, u)
        push!(vlist, v)
        push!(keylist, key)
    end
    return ArcIter(len, ulist, vlist, keylist)
end

"""
    arcs(g, u, v)
Returns the arcs of a graph g starting from u, ending at v.
The default value of u and v are Colon() to collect every arc.
If u is an iterable container, v must be Colon, vice versa.
"""
arcs(g::Graph, u=:, v=:) = ArcIter(g, u, v)

function getindex(ait::ArcIter, inds)
    if isa(inds, Integer)
        return Arc(ait.ulist[inds], ait.vlist[inds], ait.keylist[inds])
    else
        res = Arc[]
        for i in inds
            push!(res, getindex(ait, i))
        end
        return res
    end
end

firstindex(ait::ArcIter) = 1
lastindex(ait::ArcIter) = ait.len

function has_attrib(g::Graph, attrib::Symbol)
    for a in arcs(g)
        if !(attrib in keys(g[a]))
            return false
        end
    end
    return true
end


text:
facts("Volume Image") do
    fname = joinpath(dirname(@__FILE__), "../../data", "test-3d.dat")
    t = read_VolumeImage(fname)
    t2 = read_VolumeImage(fname)

    context("Reading") do
	t = read_VolumeImage(fname)
	@fact isa(t, EEG.VolumeImage) --> true
    end

    context("Create") do
        # Array
	n = VolumeImage(t.data, t.units, t.x, t.y, t.z, t.t, t.method, t.info, t.coord_system)
	@fact isa(n, EEG.VolumeImage) --> true
        @fact isequal(n.data, t.data) --> true

        # Vector
        d = zeros(length(t.data))
        x = zeros(length(t.data))
        y = zeros(length(t.data))
        z = zeros(length(t.data))
        s = zeros(length(t.data))  # t is already taken
        i = 1
        for xi in 1:length(t.x)
            for yi in 1:length(t.y)
                for zi in 1:length(t.z)
                    for ti in 1:length(t.t)
                        d[i] = float(t.data[xi, yi, zi, ti])
                        x[i] = float(t.x[xi])
                        y[i] = float(t.y[yi])
                        z[i] = float(t.z[zi])
                        s[i] = float(t.t[ti])
                        i += 1
                    end
                end
            end
        end
        n2 = VolumeImage(d, t.units, x, y, z, s, t.method, t.info, t.coord_system)
        @fact isa(n2, EEG.VolumeImage) --> true
        @fact isequal(n2.data, t.data) --> true
        @fact isequal(n2.x, t.x) --> true
        @fact isequal(n2.y, t.y) --> true
        @fact isequal(n2.z, t.z) --> true
        @fact isequal(n2.t, t.t) --> true

    end

    context("Maths") do

	@fact isequal(t, t2) --> true

        @fact t + t --> t * 2
	t4 = t * 4
	@fact t4 / 2 --> t + t
	@fact t4 / 2 --> t4 - t - t

	t2 = read_VolumeImage(fname)
	t2.units = "A/m^3"
	@fact_throws t + t2

    end

    context("Min/Max") do

	@fact maximum(t) --> 33.2692985535
	@fact minimum(t) --> -7.5189352036

	b = t * 2

	@fact b.data[1, 1, 1] --> 2 * b.data[1, 1, 1]
	@fact b.data[2, 2, 2] --> 2 * b.data[2, 2, 2]

	c = mean([b, t])

	@fact size(c.data) --> size(b.data)
	@fact c.data[1, 1, 1] --> (b.data[1, 1, 1] + t.data[1, 1, 1]) / 2

	@fact maximum([b, t]) --> maximum(t) * 2
	@fact minimum([b, t]) --> minimum(b)

	n = normalise([b, t])

	@fact maximum(n) --> less_than_or_equal(1.0)
	@fact minimum(n) --> greater_than_or_equal(-1.0)

    end

    context("Dimension checks") do
    	t2 = deepcopy(t); t2.x = t2.x[1:3]
    	@fact_throws KeyError EEG.dimensions_equal(t, t2)
    	t2 = deepcopy(t); t2.y = t2.y[1:3]
    	@fact_throws KeyError EEG.dimensions_equal(t, t2)
    	t2 = deepcopy(t); t2.z = t2.z[1:3]
    	@fact_throws KeyError EEG.dimensions_equal(t, t2)
    	t2 = read_VolumeImage(joinpath(dirname(@__FILE__), "../../data", "test-4d.dat"))
    	@fact_throws KeyError EEG.dimensions_equal(t, t2)
    end

    context("Average") do
    	s = mean(t)
    end

    context("Find dipoles") do

    	t = read_VolumeImage(fname)
    	dips = find_dipoles(mean(t))
    	@fact size(dips) --> (9,)

    	dips = EEG.new_dipole_method(mean(t))
    	@fact size(dips) --> (9,)

    	fname = joinpath(dirname(@__FILE__), "../../data", "test-4d.dat")
    	t2 = read_VolumeImage(fname)
    	dips = find_dipoles(mean(t2))
    	@fact size(dips) --> (3,)

    	# Test dipoles are returned in order of size
    	@fact issorted([dip.size for dip in dips], rev = true) --> true
    end

    context("Best Dipole") do

        dips = find_dipoles(mean(t))
        # Plots.pyplot(size=(1400, 400))
        # p = plot(t, c = :inferno)
        # p = plot(p, Talairach(-0.04, 0.01, 0.02))
        bd = best_dipole(Talairach(-0.05, 0, 0.01), dips)
        #p = plot(p, Talairach(-0.05, 0, 0.01), c = :red)

        @fact float(bd.x) --> roughly(-0.0525, atol =  0.001)
        @fact float(bd.y) --> roughly(-0.00378, atol =  0.001)
        @fact float(bd.z) --> roughly(0.0168099, atol =  0.001)

        # Take closest
        bd = best_dipole(Talairach(-0.05, 0, 0.01), dips, maxdist = 0.0015)
        @fact float(bd.x) --> roughly(-0.0525, atol =  0.001)
        @fact float(bd.y) --> roughly(-0.00378, atol =  0.001)
        @fact float(bd.z) --> roughly(0.0168099, atol =  0.001)

        # Take only valid dipole
        dists = [euclidean(Talairach(-0.05, 0, 0.01), dip) for dip=dips]
        bd = best_dipole(Talairach(-0.05, 0, 0.01), dips, maxdist = 0.015)
        @fact float(bd.x) --> roughly(-0.0525, atol =  0.001)
        @fact float(bd.y) --> roughly(-0.00378, atol =  0.001)
        @fact float(bd.z) --> roughly(0.0168099, atol =  0.001)

        bd = best_dipole(Talairach(-0.05, 0, 0.01), Dipole[])
        @fact isnan(bd) --> true

    end

    context("Show") do
    	show(t)
    	show(normalise(t))
    	t.info["Regularisation"] = 1.2
    	show(t)
    end

    context("Plotting") do
    	EEG.plot(mean(t))
    	EEG.plot(mean(t), min_val = 0, max_val = 50)
        EEG.plot(mean(t), elp = joinpath(dirname(@__FILE__), "../../data", "test.elp"))
    	p = EEG.plot(mean(t), threshold = 24 )

        context("Overlay dipole") do

            dips = find_dipoles(mean(t))
            p = plot(p, dips)
            p = plot(p, dips[1], c = :red)

        end
    end
end