text:
# LS-RTM of the 2D Marmousi model using stochastic gradient descent
# Author: pwitte.slim@gmail.com
# Date: December 2018
#
# Warning: The examples requires ~40 GB of memory per shot if used without optimal checkpointing.
#

using ClusterManagers, Distributed, ArgParse; addprocs_slurm(2)

### Process command line args
function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--time_order", "-t"
            help = "Time order"
            arg_type = Int
            default = 2
        "--velmod_dir"
            help = "Velocity model directory"
            default = ""
        "--dobs_dir"
            help = "Observed data directory"
            default = ""
        "--output_dir"
            help = "Outpur data directory"
            default = ""
    end
    return parse_args(s)
end

@everywhere using JUDI, HDF5, PyPlot, JLD, SegyIO, Random, Statistics, LinearAlgebra, Logging
parsed_args = parse_commandline()

time_order = parsed_args["time_order"]
velmod_dir = parsed_args["velmod_dir"]
dobs_dir = parsed_args["dobs_dir"]
output_dir = parsed_args["output_dir"]

# Load migration velocity model
n1, d1, o1, m, m0, rho, rho0, qp, qp0 = read(h5open(velmod_dir, "r"), "n", "d", "o", "m", "m0", "rho", "rho0", "qp", "qp0")

# Set up model structure
n=(n1[1], n1[2]);
d=(d1[1], d1[2]);
o=(o1[1], o1[2]);

# rho0 = ones(Float32,n)

model0 = Model(n, d, o, m0, rho0, qp0);

# Load data
block = segy_read(dobs_dir)
d_lin = judiVector(block)   # linearized observed data

# Set up wavelet
f0 = 0.015
src_geometry = Geometry(block; key = "source", segy_depth_key = "SourceDepth")
wavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], f0)    # 30 Hz wavelet
q = judiVector(src_geometry, wavelet)

# Set up info structure
ntComp = get_computational_nt(q.geometry, d_lin.geometry, model0)  # no. of computational time steps
info = Info(prod(model0.n), d_lin.nsrc, ntComp)

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

# Setup operators
opt = Options(optimal_checkpointing=false, f0=f0)  # ~40 GB of memory per source w/o checkpointing
M = judiModeling(model0, q.geometry, d_lin.geometry; options=opt)
J = judiJacobian(M, q)

# Stochastic gradient
x = zeros(Float32, prod(model0.n))
# batchsize = 50
niter = 20
fval = zeros(Float32, niter)

# Main loop
for j = 1: niter
    println("Iteration: ", j)

    # Compute residual and gradient
    r = J*x - d_lin
    g = adjoint(J)*r

    g = reshape(g, model0.n)
    g[:, 1:50] .= 0f0;

    # Step size and update variable
    fval[j] = .5f0*norm(r)^2

    println("fval: ", fval[j])
    t = norm(r)^2/norm(g)^2

    global x = x .- t .* g
    x = reshape(x, model0.n)
    x[:, 1:50] .= 0f0;
end

fhistory_dir = joinpath(output_dir, "fhistory_viscoacoustic_sgd_sls_2nd.bin")
write(fhistory_dir, htol.(transpose(fval)));

x_m = x
x_m_reshape = reshape(x_m, model0.n)
imagem_dir = joinpath(output_dir, "xm_final_sls_2nd.bin")
write(imagem_dir, htol.(transpose(x_m_reshape)));


text:
module CategoricalAlgebra

using Reexport

include("FreeDiagrams.jl")
include("Limits.jl")
include("Subobjects.jl")
include("Sets.jl")
include("FinSets.jl")
include("Matrices.jl")
include("FinRelations.jl")
include("CSets.jl")
include("GraphCategories.jl")
include("StructuredCospans.jl")
include("CommutativeDiagrams.jl")
include("CatElements.jl")
include("DataMigration.jl")
include("DPO.jl")

@reexport using .FreeDiagrams
@reexport using .CommutativeDiagrams
@reexport using .Limits
@reexport using .Subobjects
@reexport using .CSets
@reexport using .StructuredCospans
@reexport using .CatElements
@reexport using .DataMigration
@reexport using .DPO

end

ext: jl


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct TGFixed <: TurbineGov
        efficiency::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of a fixed Turbine Governor that returns a fixed mechanical torque
 given by the product of P_ref*efficiency

# Arguments
- `efficiency::Float64`:  Efficiency factor that multiplies P_ref, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct TGFixed <: TurbineGov
    " Efficiency factor that multiplies P_ref"
    efficiency::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function TGFixed(efficiency, ext=Dict{String, Any}(), )
    TGFixed(efficiency, ext, Vector{Symbol}(), 0, InfrastructureSystemsInternal(), )
end

function TGFixed(; efficiency, ext=Dict{String, Any}(), )
    TGFixed(efficiency, ext, )
end

# Constructor for demo purposes; non-functional.
function TGFixed(::Nothing)
    TGFixed(;
        efficiency=0,
        ext=Dict{String, Any}(),
    )
end

"""Get TGFixed efficiency."""
get_efficiency(value::TGFixed) = value.efficiency
"""Get TGFixed ext."""
get_ext(value::TGFixed) = value.ext
"""Get TGFixed states."""
get_states(value::TGFixed) = value.states
"""Get TGFixed n_states."""
get_n_states(value::TGFixed) = value.n_states
"""Get TGFixed internal."""
get_internal(value::TGFixed) = value.internal


text:
using SpineInterface
using JuMP
using Cbc
using Plots

path = "sqlite:///$(@__DIR__)/input_data/input_data.sqlite"

using_spinedb(path)
# now we have a direct handle access entities of the sqlite db from the path

# set up a JuMP model
model = JuMP.Model()
set_optimizer(model, Cbc.Optimizer)
set_optimizer_attributes(model, "LogLevel" => 1, "PrimalTolerance" => 1e-7)

dates = map(x -> x[1], collect(price(commodity=commodity(:elec))))
n_dates = length(dates)

# Unit states, integer
@variable(model, unit_states[dates, node(), unit()], Bin, container=DenseAxisArray)

# units power/production
@variable(model, unit_powers[dates, node(), unit()], container=DenseAxisArray)

# Set unit powers to either 0, or between min and max
for n in node()
    for date in dates, uname in unit()
        if uname in unit_output_node(node=n)
            @constraint(model, unit_powers[date, n, uname] .<= unit_states[date, n, uname] .* capacity(unit=uname) .* max_load(unit=uname))
            @constraint(model, unit_powers[date, n, uname] .>= unit_states[date, n, uname] .* capacity(unit=uname) .* min_load(unit=uname))
        else
            @constraint(model, unit_powers[date, n, uname] .== 0)
        end
    end
end

# Dummy power in case demand cannot be met by production units
@variable(model, dummy_power[dates, node()])
for date in dates, n in node()
    @constraint(model, dummy_power[date, n] >= 0)
end

# General constraints
# Production should meet demand for nodes with balance..
for n in node()
    flow_val = flow(node=n)
    balance_val = balance(node=n)
    if typeof(flow_val) != Nothing && balance_val != 0
        flow_val = map(x -> x[2], collect(flow_val))
        for d in 1:n_dates
            @constraint(model, sum(unit_powers[dates[d], n, :]) + dummy_power[dates[d], n] .== -1 .* flow_val[d])
        end
    end
end

# Resource limitations
for n in node()
   for d in 1:n_dates, uname in unit()
       resource_type = source(unit=uname)
       if typeof(resource_type) != Nothing
           resource_eff = eff(unit=uname)
           local resource_flow = flow(node=node(resource_type))
           resource_flow = map(x -> x[2], collect(resource_flow))
           @constraint(model, unit_powers[dates[d], n, uname] .<= resource_flow[d].* resource_eff)
       end
   end
end

#VOM_cost = VOM_cost .* unit states?
@expression(model, vom_costs[d in 1:n_dates, n in node(), uname in unit()], unit_states[dates[d], n, uname] .* VOM_cost(unit=uname))

# fuel_cost = power ./ eff .* fuel price
@expression(model, fuel_costs[d in 1:n_dates, n in node(), uname in unit()], unit_powers[dates[d], n, uname] ./ eff(unit = uname) .* map(x -> x[2], collect(price(commodity=commodity_node(node=n)[1])))[d])

# Dummy power cost
@expression(model, dummy_cost[d in 1:n_dates, n in node()], dummy_power[dates[d], n] .* 1000000)

# sell profit
@expression(model, sell_profit[d in 1:n_dates, n in node(), uname in unit()], unit_powers[dates[d], n, uname] .* map(x -> x[2], collect(price(commodity=commodity_node(node=unit_output_node(unit=uname))[1])))[d])

#Objective function
@objective(model, Min, sum(vom_costs) + sum(fuel_costs) + sum(dummy_cost) - sum(sell_profit))

# optimize the model
optimize!(model)

max_stars_repo_name: vttresearch/Predicer

max_issues_repo_name: vttresearch/Predicer

max_forks_repo_name: vttresearch/Predicer


text:
module Reproducible

using Markdown

#-----------------------------------------------------------------------# build
"""
    build(path::String [, builddir::String]; frontmatter::String = "")
    build(path::Vector{String} [, builddir::String]; frontmatter::String = "")

Evaluate the markdown document(s) in `path` and put the output in `builddir`, beginning
with `frontmatter` (since Julia's markdown parser does not support it).
"""
function build(path::String, builddir::String = joinpath(dirname(path), "build"); frontmatter::String = "", toc=false)
    TempModule = Main.eval(:(module __Temp__; const __builddir__ = $builddir end))
    !isdir(builddir) && mkdir(builddir)
    file = touch(joinpath(builddir, basename(path)))
    open(file, "w") do io
        toc && write(io, maketoc(path))
        !isempty(frontmatter) && write(io, "---\n$(strip(frontmatter))\n---\n\n")
        content = Markdown.parse(read(path, String)).content
        for (i, x) in enumerate(content)
            @info "Reproducible Build: $i / $(length(content))"
            write(io, markdown2string(x, TempModule, builddir) * "\n")
        end
    end
    return file
end

function build(paths::Vector{String}, builddir::String = joinpath(dirname(paths[1]), "build"))
    files = []
    for path in paths
        push!(files, build(path, builddir))
    end
    files
end

function markdown2string(x #=Markdown=#, mod::Module, builddir::String)
    if isa(x, Markdown.Code)
        isa(Meta.parse(x.language), Expr) ? code2string(x, mod, builddir) : Markdown.plain(x)
    else
        Markdown.plain(x)
    end
end

function code2string(x::Markdown.Code, mod::Module, builddir::String)
    ex = Meta.parse(x.language)
    lang = ex.args[1]
    if lang == :julia
        v = Val(Meta.parse(x.language).args[2])  # renderer
        render(CodeBlock(x.code, mod), v; builddir=builddir, _kws(ex.args[3:end])...)
    else
        @warn "Reproducible doesn't know how to eval a code block with language $lang...skipping"
        Markdown.plain(x)
    end
end

# create named tuple from array of expressions e.g. :(hide = true)
function _kws(args)
    vals = [eval(a.args[2]) for a in args]
    nt = NamedTuple{tuple([a.args[1] for a in args]...), Tuple{eltype.(vals)...}}(tuple(vals...))
end

#-----------------------------------------------------------------------# CodeRow
"""
Object that represents a block of code.  Fields are:

- `input (String)`: the input block of code
- `output (Any)`: the result of `eval(parse(input))`
- `repldisplay (String)`: string of what would get sent to REPL
"""
struct CodeRow
    input::String
    output::Any
    repldisplay::String
end

#-----------------------------------------------------------------------# CodeBlock
"""
    CodeBlock(code::String, mod::Module)

Object to represent code for a **Reproducible** renderer.  The constructor parses/evaluates
the code inside the provided module.
"""
struct CodeBlock
    rows::Vector{CodeRow}
end
function CodeBlock(code::String, mod::Module, args...; display_size = displaysize(), limit = true)
    n = 1
    rows = CodeRow[]
    while n < length(code)
        nold = n
        ex, n = Meta.parse(code, n)
        input = code[nold:n-1]
        output = @eval(mod, $ex)
        io = IOContext(IOBuffer(), :display_size => display_size, :limit => limit, args...)
        show(io, MIME"text/plain"(), output)
        out = (endswith(strip(input), ';') || output == nothing) ? "" : String(take!(io.io))
        push!(rows, CodeRow(input, output, out))
    end
    CodeBlock(rows)
end
Base.getindex(o::CodeBlock, i) = o.rows[i]
codestring(o::CodeBlock) = join([r.input for r in o.rows])

#-----------------------------------------------------------------------# utils
block(s::String, lang="") = "```$lang\n$(strip(s))\n```\n"

#-----------------------------------------------------------------------# includes
include("toc.jl")
include("render.jl")

end # module


text:
module DirichletTest

using Test
using ForneyLab
using ForneyLab: outboundType, isApplicable, prod!, unsafeMean, unsafeLogMean, unsafeVar, vague, dims, naturalParams, standardDistribution
using ForneyLab: SPDirichletOutNP, VBDirichletOut, VBDirichletIn1
using SpecialFunctions: digamma

@testset "Dirichlet ProbabilityDistribution and Message construction" begin
    @test ProbabilityDistribution(Multivariate, Dirichlet, a=[2.0, 3.0]) == ProbabilityDistribution{Multivariate, Dirichlet}(Dict(:a=>[2.0, 3.0]))
    @test ProbabilityDistribution(MatrixVariate, Dirichlet, a=[2.0 3.0; 4.0 5.0]) == ProbabilityDistribution{MatrixVariate, Dirichlet}(Dict(:a=>[2.0 3.0; 4.0 5.0]))
    @test_throws Exception ProbabilityDistribution(Univariate, Dirichlet)
    @test ProbabilityDistribution(Dirichlet, a=[2.0, 3.0]) == ProbabilityDistribution{Multivariate, Dirichlet}(Dict(:a=>[2.0, 3.0]))
    @test ProbabilityDistribution(Dirichlet) == ProbabilityDistribution{Multivariate, Dirichlet}(Dict(:a=>[1.0, 1.0, 1.0]))
    @test Message(Dirichlet) == Message{Dirichlet, Multivariate}(ProbabilityDistribution{Multivariate, Dirichlet}(Dict(:a=>[1.0, 1.0, 1.0])))
    @test Message(Multivariate, Dirichlet) == Message{Dirichlet, Multivariate}(ProbabilityDistribution{Multivariate, Dirichlet}(Dict(:a=>[1.0, 1.0, 1.0])))
    @test Message(MatrixVariate, Dirichlet) == Message{Dirichlet, MatrixVariate}(ProbabilityDistribution{MatrixVariate, Dirichlet}(Dict(:a=>ones(3,3))))
    @test_throws Exception Message(Univariate, Dirichlet)
end

@testset "dims" begin
    @test dims(ProbabilityDistribution(Multivariate, Dirichlet, a=[2.0, 2.0, 2.0])) == (3,)
    @test dims(ProbabilityDistribution(MatrixVariate, Dirichlet, a=[2.0 2.0 2.0; 2.0 2.0 2.0])) == (2,3)
end

@testset "vague" begin
    @test vague(Dirichlet, (3,)) == ProbabilityDistribution(Dirichlet, a=ones(3))
    @test vague(Dirichlet, (2,3)) == ProbabilityDistribution(MatrixVariate, Dirichlet, a=ones(2,3))
end

@testset "unsafe mean and variance" begin
    @test unsafeMean(ProbabilityDistribution(Multivariate, Dirichlet, a=[2.0, 2.0])) == [0.5, 0.5]
    @test unsafeMean(ProbabilityDistribution(MatrixVariate, Dirichlet, a=[2.0 1.0; 2.0 3.0])) == [0.5 0.25; 0.5 0.75]
    @test unsafeLogMean(ProbabilityDistribution(Multivariate, Dirichlet, a=[2.0, 3.0])) == [digamma(2.0), digamma(3.0)] .- digamma(5.0)
    @test unsafeLogMean(ProbabilityDistribution(MatrixVariate, Dirichlet, a=[2.0 4.0; 3.0 5.0])) == [digamma(2.0) digamma(4.0); digamma(3.0) digamma(5.0)] - [digamma(5.0) digamma(9.0); digamma(5.0) digamma(9.0)]
    @test unsafeVar(ProbabilityDistribution(Multivariate, Dirichlet, a=[2.0, 2.0])) == [0.05, 0.05]
end

@testset "prod!" begin
    # Multivariate
    @test ProbabilityDistribution(Multivariate, Dirichlet, a=[2.0, 2.0]) * ProbabilityDistribution(Multivariate, Dirichlet, a=[2.0, 3.0]) == ProbabilityDistribution(Multivariate, Dirichlet, a=[3.0, 4.0])
    @test ProbabilityDistribution(Multivariate, Dirichlet, a=[1.0, 2.0, 3.0]) * ProbabilityDistribution(Multivariate, PointMass, m=[0.1, 0.8, 0.1]) == ProbabilityDistribution(Multivariate, PointMass, m=[0.1, 0.8, 0.1])
    @test ProbabilityDistribution(Multivariate, PointMass, m=[0.1, 0.8, 0.1]) * ProbabilityDistribution(Multivariate, Dirichlet, a=[1.0, 2.0, 3.0]) == ProbabilityDistribution(Multivariate, PointMass, m=[0.1, 0.8, 0.1])

    # MatrixVariate
    @test ProbabilityDistribution(MatrixVariate, Dirichlet, a=[2.0 2.0; 3.0 3.0]) * ProbabilityDistribution(MatrixVariate, Dirichlet, a=[2.0 3.0; 4.0 5.0]) == ProbabilityDistribution(MatrixVariate, Dirichlet, a=[3.0 4.0; 6.0 7.0])
    @test ProbabilityDistribution(MatrixVariate, Dirichlet, a=[1.0 3.0; 2.0 4.0]) * ProbabilityDistribution(MatrixVariate, PointMass, m=[0.1 0.3; 0.9 0.7]) == ProbabilityDistribution(MatrixVariate, PointMass, m=[0.1 0.3; 0.9 0.7])
    @test ProbabilityDistribution(MatrixVariate, PointMass, m=[0.1 0.3; 0.9 0.7]) * ProbabilityDistribution(MatrixVariate, Dirichlet, a=[1.0 3.0; 2.0 4.0]) == ProbabilityDistribution(MatrixVariate, PointMass, m=[0.1 0.3; 0.9 0.7])
end

@testset "log pdf" begin
    @test isapprox(logPdf(ProbabilityDistribution(Multivariate, Dirichlet, a=[0.2,3.0,1.5]),[2,3,7]), 3.2556382883760024)
    @test isapprox(logPdf(ProbabilityDistribution(MatrixVariate, Dirichlet, a=[0.2 1.4; 3.0 1.8]),[2 7; 3 3]), 3.0442561618507087)
end

@testset "natural parameters" begin
    # Multivariate
    d = ProbabilityDistribution(Multivariate, Dirichlet, a=[1.5, 3.0, 5.0])
    η = naturalParams(d)
    s = standardDistribution(Multivariate, Dirichlet, η=η)
    @test d.params[:a] == s.params[:a] # Test conversion consistency

    x = [[0.2, 0.6, 0.2], [0.1, 0.4, 0.5]]
    d_x = logPdf.([d], x)
    η_x = logPdf.(Multivariate, Dirichlet, x; η=η)
    @test isapprox(d_x, η_x) # Test pdf consistency

    # MatrixVariate
    d = ProbabilityDistribution(MatrixVariate, Dirichlet, a=[1.5 3.0; 5.0 4.0; 4.0 6.0])
    η = naturalParams(d)
    s = standardDistribution(MatrixVariate, Dirichlet, η=η, dims=(3,2))
    @test d.params[:a] == s.params[:a] # Test conversion consistency

    x = [[0.2 0.1; 0.6 0.4; 0.2 0.5], [0.1 0.4; 0.5 0.1; 0.4 0.5]]
    d_x = logPdf.([d], x)
    η_x = logPdf.(MatrixVariate, Dirichlet, x; η=η)
    @test isapprox(d_x, η_x) # Test pdf consistency
end

#-------------
# Update rules
#-------------

@testset "SPDirichletOutNP" begin
    @test SPDirichletOutNP <: SumProductRule{Dirichlet}
    @test outboundType(SPDirichletOutNP) == Message{Dirichlet}
    @test isApplicable(SPDirichletOutNP, [Nothing, Message{PointMass}])

    @test ruleSPDirichletOutNP(nothing, Message(Multivariate, PointMass, m=[2.0, 3.0])) == Message(Multivariate, Dirichlet, a=[2.0, 3.0])
    @test ruleSPDirichletOutNP(nothing, Message(MatrixVariate, PointMass, m=[2.0 3.0; 4.0 5.0])) == Message(MatrixVariate, Dirichlet, a=[2.0 3.0; 4.0 5.0])
end

@testset "VBDirichletOut" begin
    @test VBDirichletOut <: NaiveVariationalRule{Dirichlet}
    @test outboundType(VBDirichletOut) == Message{Dirichlet}
    @test isApplicable(VBDirichletOut, [Nothing, ProbabilityDistribution])

    @test ruleVBDirichletOut(nothing, ProbabilityDistribution(Multivariate, PointMass, m=[2.0, 3.0])) == Message(Multivariate, Dirichlet, a=[2.0, 3.0])
    @test ruleVBDirichletOut(nothing, ProbabilityDistribution(MatrixVariate, PointMass, m=[2.0 3.0; 4.0 5.0])) == Message(MatrixVariate, Dirichlet, a=[2.0 3.0; 4.0 5.0])
end

@testset "VBDirichletIn1" begin
    @test VBDirichletIn1 <: NaiveVariationalRule{Dirichlet}
    @test outboundType(VBDirichletIn1) == Message{Function}
    @test isApplicable(VBDirichletIn1, [ProbabilityDistribution, Nothing])
end

@testset "averageEnergy and differentialEntropy" begin
    # Multivariate
    @test isapprox(differentialEntropy(ProbabilityDistribution(Multivariate, Dirichlet, a=[2.0, 3.0])), averageEnergy(Dirichlet, ProbabilityDistribution(Multivariate, Dirichlet, a=[2.0, 3.0]), ProbabilityDistribution(Multivariate, PointMass, m=[2.0, 3.0])))
    @test isapprox(averageEnergy(Dirichlet, ProbabilityDistribution(Multivariate, Dirichlet, a=[4.0, 5.0]), ProbabilityDistribution(Multivariate, PointMass, m=[2.0, 3.0])), averageEnergy(Beta, ProbabilityDistribution(Univariate, Beta, a=4.0, b=5.0), ProbabilityDistribution(Univariate, PointMass, m=2.0), ProbabilityDistribution(Univariate, PointMass, m=3.0)))
    @test isapprox(averageEnergy(Dirichlet,ProbabilityDistribution(Multivariate,Dirichlet,a=[2.0,3.0,4.0]),ProbabilityDistribution(Multivariate,SampleList,s=[[1.0,2.0,1.0],[3.,3.,1.],[2.,2.,2]],w=[0.1,0.7,0.3])),0.5093876870003795)

    # MatrixVariate
    @test differentialEntropy(ProbabilityDistribution(MatrixVariate, Dirichlet, a=[2.0 3.0; 4.0 5.0])) == differentialEntropy(ProbabilityDistribution(Multivariate, Dirichlet, a=[2.0, 4.0])) + differentialEntropy(ProbabilityDistribution(Multivariate, Dirichlet, a=[3.0, 5.0]))
    @test averageEnergy(Dirichlet, ProbabilityDistribution(MatrixVariate, Dirichlet, a=[2.0 3.0; 4.0 5.0]), ProbabilityDistribution(MatrixVariate, PointMass, m=[6.0 7.0; 8.0 9.0])) == averageEnergy(Dirichlet, ProbabilityDistribution(Multivariate, Dirichlet, a=[2.0, 4.0]), ProbabilityDistribution(Multivariate, PointMass, m=[6.0, 8.0])) + averageEnergy(Dirichlet, ProbabilityDistribution(Multivariate, Dirichlet, a=[3.0, 5.0]), ProbabilityDistribution(Multivariate, PointMass, m=[7.0, 9.0]))
end

end # module


text:
include("types.jl")
include("volume.jl")
include("singlecomponent/singlecomponent.jl")
include("multicomponent/multicomponent.jl")


text:
module Gateway

using Sockets

export getgateway

function parse_linux_proc_net_route(file::IOStream)
	# /proc/net/route file:
	# Iface	Destination	Gateway 	Flags	RefCnt	Use	Metric	Mask		MTU	Window	IRTT
	# wlo1	00000000	0100A8C0	0003	0	0	600	00000000	0	0	0
	# wlo1	0000A8C0	00000000	0001	0	0	600	00FFFFFF	0	0	0

	sep = "\t"
	field = 14

	scanner = read(file, String)
	tokens = split(scanner, sep)
	if length(tokens) <= field
		throw(error("No gateway found"))
	end

	gateway = parse(UInt32, tokens[field], base=16)
	gateway = hton(gateway) # Convert hex address to big endian
	Sockets.IPv4(gateway)
end

function parse_windows_route(output::IOStream)
	# Windows "route print 0.0.0.0" output:
	# ===========================================================================
	# Interface List
	# 3 ...02 16 4a a8 72 ca ...... Intel(R) PRO/100 VE Network Connection
	# 1 ........................... Software Loopback Interface 1
	# ===========================================================================
	# IPv4 Route Table
	# ===========================================================================
	# Active Routes:
	# Network Destination        Netmask          Gateway       Interface  Metric
	#           0.0.0.0          0.0.0.0      192.168.0.1    192.168.0.100     20
	# ===========================================================================
	#
	# Get to "Active Routes" section and jump 2 lines below and pick address from 3rd column

	sep = "\n"
	column = 3

	scanner = read(output, String)
	tokens = split(output, sep)
	sep = 0
	for (idx, line) in enumerate(tokens)
		if sep == 3
			if length(tokens) <= idx  + 2
				throw(error("No gateway found"))
			end

			fields = split(tokens[idx+2])
			if length(fields) < 3
				throw(error("No gateway found"))
			end

			return Sockets.IPv4(fields[column])
		end

		if startswith(line, "=======")
			sep += 1
			continue
		end
	end
	throw(error("No gateway found"))
end

function parse_osx_route(output::IOStream)
	# Darwin route frame:
	#    route to: default
	# destination: default
	#        mask: default
	#     gateway: 192.168.1.1

	sep = "\n"
	column = 2

	scanner = read(output, String)
	tokens = split(scanner, sep)

	for line in tokens
		fields = split(line)
		if length(fields) >=2 && fields[column - 1] == "gateway:"
			return Sockets.IPv4(fields[column])
		end
	end

	throw(error("No gateway found"))
end

function parse_unix_netstat(output::IOStream)
	# For unix based OS such as *BSD, solaris etc
	# netstat -rn output:
	# Routing tables
	#
	# Internet:
	# Destination        Gateway            Flags      Netif Expire
	# default            10.88.88.2         UGS         em0
	# 10.88.88.0/24      link#1             U           em0
	# 10.88.88.148       link#1             UHS         lo0
	# 127.0.0.1          link#2             UH          lo0
	#
	# Internet6:
	# Destination                       Gateway                       Flags      Netif Expire
	# ::/96                             ::1                           UGRS        lo0
	# ::1                               link#2                        UH          lo0
	# ::ffff:0.0.0.0/96                 ::1                           UGRS        lo0
	# fe80::/10                         ::1                           UGRS        lo0
	# ...
	sep = "\n"
	column = 2
	scanner = read(output, String)
	tokens = split(scanner, sep)

	for line in tokens
		fields = split(line)
		if length(fields) >=2 && fields[column - 1] == "default"
			return Sockets.IPv4(fields[column])
		end
	end

	throw(error("No gateway found"))
end

function getgateway()
	if Sys.islinux()
		open("/proc/net/route") do file
			return parse_linux_proc_net_route(file)
		end
	elseif Sys.iswindows()
		output = read(`cmd /c route print 0.0.0.0`)
		return parse_windows_route(output)
	elseif Sys.isapple()
		output = read(`/sbin/route -n get 0.0.0.0`)
		return parse_osx_route(output)
	elseif Sys.isbsd()
		output = read(`netstat -rn`)
		return parse_unix_netstat(output)
	else
		println("Operating system not supported please file a issue on github")
	end
end

end #module


text:
@testset "Statistics" begin
  @testset "Solution" begin
    grid2D = RegularGrid{Float64}(3,3)
    realizations = Dict(:value => [i*ones(npoints(grid2D)) for i in 1:3])
    solution2D = SimulationSolution(grid2D, realizations)

    # mean
    mean2D = mean(solution2D)
    @test mean2D[1,:value] == 2.
    @test domain(mean2D) == solution2D.domain

    # variance
    variance2D = var(solution2D)
    @test variance2D[1,:value] == 1.
    @test domain(variance2D) == solution2D.domain

    # quantile (scalar)
    p = 0.5
    quantile2D = quantile(solution2D, p)
    @test quantile2D[1,:value] == 2.
    @test domain(quantile2D) == solution2D.domain

    # quantile (vector)
    ps = [0.0, 0.5, 1.0]
    quantiles2D = quantile(solution2D, ps)
    @test eltype(quantiles2D) <: SpatialStatistic

    # show methods
    @test sprint(show, mean2D) == "2D SpatialStatistic"
    @test sprint(show, MIME"text/plain"(), mean2D) == "2D SpatialStatistic\n  domain: 3×3 RegularGrid{Float64,2}\n  variables: value"
  end

  @testset "Data" begin
    # load data with bias towards large values (gold mine)
    sdata = readgeotable(joinpath(datadir,"clustered.csv"), coordnames=[:x,:y])

    # spatial mean
    μn = mean(sdata[:Au])
    μs = mean(sdata, :Au)
    @test abs(μn - 0.5) > abs(μs - 0.5)
    @test mean(sdata)[:Au] ≈ μs

    # spatial variance
    σn = var(sdata[:Au])
    σs = var(sdata, :Au)
    @test σn ≤ σs
    @test var(sdata)[:Au] ≈ σs

    # spatial quantile
    qn = quantile(sdata[:Au], 0.5)
    qs = quantile(sdata, :Au, 0.5)
    @test qn ≥ qs
    @test quantile(sdata, 0.5)[:Au] ≈ qs
  end
end


text:
using DataFrames
using CSV

cd("/Users/kiristern/Documents/GitHub/SDM/")

predators = CSV.read("data/preds/originals/Nasua_nasua2.csv")
predators = select(predators, [:species, :decimalLatitude, :decimalLongitude, :year])

# select only ursus arctos from imported all predators csv file
#df = predators[predators.species .== "Ursus arctos", :]

p1 = CSV.read("data/prey_class/Diplopoda2.csv")
p1 = select(p1,[:class, :decimalLatitude, :decimalLongitude, :year])
rename!(p1, :class => :species)

#create one dataframe with predators/prey dataframes
pred_prey = vcat(predators, p1)

pred_prey = dropmissing!(pred_prey, [:decimalLatitude, :decimalLongitude, :year])

CSV.write("data/pred_prey/Nasua_nasua.csv", pred_prey)

show(by(pred_prey, :species, :species => length), allcols=true)

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
using StructuralCausalModels, Test

ProjDir = @__DIR__
cd(ProjDir) #do

d_string = "dag {I -> Y; X -> I; X -> Y; Z -> I; Z -> X}"
dag = DAG("test_open_paths_03", d_string);

bs = basis_set(dag)

fname = joinpath(ProjDir, "test_open_paths_03.dot")
to_graphviz(dag, fname)
#Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

ap  = all_paths(dag, :X, :Y)
bp = backdoor_paths(dag, ap, :X)
adjs = adjustment_sets(dag, :X, :Y)

@testset "Open_path_03" begin

  @test bs[1][1] == :Z
  @test bs[1][2] == :Y
  @test bs[1][3:end] == [:X, :I]
  @test bp[1] == [:X, :Z, :I, :Y]
  @test adjs == [[:Z], [:I]]

end


text:
function finalinit!(params::Parameters)

  for i in 1:params.noofnodes
    params.nodes[i].clearance = max(0.0, params.nodes[i].holeradius - params.nodes[i].stringradius)
  end

  if params.remove_bit_constraints
    params.nodes[1].clearance = 999.0
  end

  if params.has_rig
    params.nodes[params.noofnodes].clearance = 0.0
  end

  params.bf[1, 1, 1] = cos(params.nodes[1].dip) * params.nodes[1].weight
  params.f[1, 1, 1] = params.bf[1, 1, 1]

  # Still figuring out if below adjustment is needed for bit
  params.fn[1] = -params.wob + sin(params.nodes[1].dip) * params.nodes[1].weight
  params.friction[1] =
    params.nodes[1].fc * sqrt(params.bf[1, 1, 1]^2 + params.bf[2, 1, 1]^2)
  params.f[3, 1, 1] = params.fn[1]

  if params.use_geometric_correction
    params.sc[:, :, 1] += params.fn[1] * params.gm[1:6, 7:12, 1]
    params.sa[:, :, 1] += params.fn[1] * params.gm[7:12, 1:6, 1]
    params.sb1[:, :, 1] += params.fn[1] * params.gm[7:12, 7:12, 1]
    params.sb2[:, :, 1] += params.fn[1] * params.gm[1:6, 1:6, 1]
    params.si[:, :, 1] = BHAtp.stiffnessi!(1, params.sb1[:, :, 1], params.sb2[:, :, 1])
  end

  for i in 2:params.noofnodes
    params.bf[1, 1, i] = cos(params.nodes[i].dip) * params.nodes[i].weight
    params.f[1, 1, i] = params.bf[1, 1, i]
    params.del[1, 1, i] = params.nodes[i].x - params.nodes[i].clearance
    params.del[2, 1, i] = params.nodes[i].y + 0.001 * params.nodes[i].clearance

    params.fn[i] = params.fn[i-1] + sin(params.nodes[i].dip) * params.nodes[i].weight
    params.friction[i] =
      params.nodes[i].fc * sqrt(params.bf[1, 1, i]^2 + params.bf[2, 1, i]^2)

    params.f[3, 1, i] = params.fn[i]

    local l = params.nodes[i].z - params.nodes[i-1].z
    params.del[3, 1, i] = params.del[3, 1, i-1] + params.fn[i] * l / params.elements[i-1].ea

    if params.use_geometric_correction
      params.sc[:, :, i] += params.fn[i] * params.gm[1:6, 7:12, i]
      params.sa[:, :, i] += params.fn[i] * params.gm[7:12, 1:6, i]
      params.sb1[:, :, i] += params.fn[i] * params.gm[7:12, 7:12, i]
      params.sb2[:, :, i] += params.fn[i] * params.gm[1:6, 1:6, i]
      params.si[:, :, i] = BHAtp.stiffnessi!(1, params.sb1[:, :, i], params.sb2[:, :, i])
    end

    params.del[4, 1, i] = (params.nodes[i].y - params.nodes[i-1].y) / l
    params.del[5, 1, i] = (params.nodes[i].x - params.nodes[i-1].x) / l
  end

  # Below vectors are initialized from the top
  params.nodes[params.noofnodes].tvd = 0.0
  params.nodes[params.noofnodes].xvert = 0.0
  params.nodes[params.noofnodes].deltainclination = 0.0
  params.nodes[params.noofnodes].deltaheading = 0.0
  params.nodes[params.noofnodes].dispew = 0.0
  params.nodes[params.noofnodes].dispns = 0.0
  for i in reverse(1:params.noofnodes-1)

    local dz = params.nodes[i].z - params.nodes[i+1].z
    params.nodes[i].tvd = params.nodes[i+1].tvd + dz * cos(params.nodes[i].inclination)
    params.nodes[i].xvert = params.nodes[i+1].xvert + dz * sin(params.nodes[i].inclination)

    # deltas expressed per 100 ft
    params.nodes[i].deltainclination =
      1200.0 * (params.nodes[i].inclination - params.nodes[i+1].inclination)/dz
    params.nodes[i].deltaheading =
      1200.0 * (params.nodes[i].heading - params.nodes[i+1].heading)/dz

    local dx = params.nodes[i].xvert - params.nodes[i+1].xvert
    params.nodes[i].dispns = params.nodes[i+1].dispns + dx * cos(params.nodes[i].heading)
    params.nodes[i].dispew = params.nodes[i+1].dispew + dx * sin(params.nodes[i].heading)
  end

  nothing
end


text:
"""
    ValLength(::Type{NTuple{N}}) = Val{N}
"""
ValLength(::Type{NTuple{N,T}}) where {N,T} = Val{N}

"""
    ValLength(::NTuple{N}) = Val(N)
"""
ValLength(::NTuple{N}) where {N} = Val(N)

ValLength(::CartesianIndex{N}) where {N} = Val(N)
ValLength(::Type{CartesianIndex{N}}) where {N} = Val{N}

push(s::Tuple, val) = (s..., val)

pushfirst(s::Tuple, val) = (val, s...)

pop(s::NTuple{N}) where {N} = ntuple(i -> s[i],
                                     Val(N-1))

popfirst(s::NTuple{N}) where {N} = ntuple(i -> s[i+1],
                                          Val(N-1))

# Permute some other type by perm
# (for example, tuple, MVector, etc.)
# as long as the constructor accepts a tuple
function permute(s::T,perm) where {T}
  return T(ntuple(i->s[perm[i]], ValLength(T)()))
end

sim(s::NTuple) = s

"""
    getperm(col1,col2)

Get the permutation that takes collection 2 to collection 1,
such that col2[p].==col1
"""
function getperm(s1, s2)
  return ntuple(i->findfirst(==(s1[i]),s2),Val(length(s1)))
end

"""
    getperm(col1,col2,col3)

Get the permutations that takes collections 2 and 3 to collection 1.
"""
function getperms(s,s1,s2)
  N = length(s)
  N1 = length(s1)
  N2 = length(s2)
  N1+N2≠N && error("Size of partial sets don't match with total set")
  perm1 = ntuple(i->findfirst(==(s1[i]),s),Val(N1))
  perm2 = ntuple(i->findfirst(==(s2[i]),s),Val(N2))
  isperm((perm1...,perm2...)) || error("Combined permutations are $((perm1...,perm2...)), not a valid permutation")
  return perm1,perm2
end

function invperm!(permres,perm)
  for i in 1:length(perm)
    permres[perm[i]] = i
  end
  return permres
end

function invperm(perm::NTuple{N,Int}) where {N}
  mpermres = MVector{N,Int}(undef)
  invperm!(mpermres,perm)
  return Tuple(mpermres)
end

function invperm(perm)
  permres = similar(perm)
  invperm!(permres,perm)
  return permres
end

"""
    is_trivial_permutation(P)

Determine if P is a trivial permutation.
"""
function is_trivial_permutation(P)
  #isperm(P) || error("Input is not a permutation")
  # TODO: use `all(n->P[n]==n,1:length(P))`?
  N = length(P)
  for n in 1:N
    P[n]!=n && return false
  end
  return true
end

# Combine a bunch of tuples
@inline flatten(x) = x
@inline flatten(x, y) = (x..., y...)
@inline flatten(x, y, z...) = (x..., flatten(y, z...)...)

function _deleteat(t,pos,i)
  i < pos && return t[i]
  return t[i+1]
end

function deleteat(t::NTuple{N},pos::Integer) where {N}
  return ntuple(i -> _deleteat(t,pos,i),Val(N-1))
end

deleteat(t::Tuple, I::Tuple{Int}) = deleteat(t, I[1])
function deleteat(t::Tuple, I::Tuple{Int, Int, Vararg{Int}})
    return deleteat_sorted(t, sort(I, rev = true))
end

deleteat_sorted(t::Tuple,pos::Int64) = deleteat(t,pos[1])
deleteat_sorted(t::Tuple,pos::Tuple{Int}) = deleteat(t,pos[1])
deleteat_sorted(t::Tuple,pos::NTuple{N,Int}) where {N} = deleteat_sorted(deleteat_sorted(t,pos[1]),Base.tail(pos))

# Make a slice of the block on the specified dimensions
# Make this a generic tupletools function (TupleTools.jl calls it getindices)
function getindices(t::Tuple,
                    I::NTuple{N,Int}) where {N}
  return ntuple(i->t[I[i]],Val(N))
end

# Taken from TupleTools.jl
"""
    sort(t::Tuple; lt=isless, by=identity, rev::Bool=false) -> ::Tuple
Sorts the tuple `t`.
"""
Base.sort(t::Tuple; lt=isless, by=identity, rev::Bool=false) = _sort(t, lt, by, rev)
@inline function _sort(t::Tuple, lt=isless, by=identity, rev::Bool=false)
    t1, t2 = _split(t)
    t1s = _sort(t1, lt, by, rev)
    t2s = _sort(t2, lt, by, rev)
    return _merge(t1s, t2s, lt, by, rev)
end
_sort(t::Tuple{Any}, lt=isless, by=identity, rev::Bool=false) = t
_sort(t::Tuple{}, lt=isless, by=identity, rev::Bool=false) = t

function _split(t::NTuple{N}) where N
    M = N>>1
    ntuple(i->t[i], M), ntuple(i->t[i+M], N-M)
end

function _merge(t1::Tuple, t2::Tuple, lt, by, rev)
    if lt(by(first(t1)), by(first(t2))) != rev
        return (first(t1), _merge(Base.tail(t1), t2, lt, by, rev)...)
    else
        return (first(t2), _merge(t1, Base.tail(t2), lt, by, rev)...)
    end
end
_merge(t1::Tuple{}, t2::Tuple, lt, by, rev) = t2
_merge(t1::Tuple, t2::Tuple{}, lt, by, rev) = t1

#function tail(t::NTuple{N})
#  return ntuple(i -> t[i+1],Val(N-1))
#end

function _insertat(t,pos,n_insert,val,i)
  if i < pos
    return t[i]
  elseif i > pos+n_insert-1
    return t[i-n_insert+1]
  end
  return val[i-pos+1]
end

"""
    insertat

Remove the value at pos and insert the elements in val
"""
function insertat(t::NTuple{N},
                  val::NTuple{M},
                  pos::Integer) where {N,M}
  return ntuple(i -> _insertat(t,pos,M,val,i),Val(N+M-1))
end

function insertat(t::NTuple{N}, val, pos::Integer) where {N}
  return insertat(t,tuple(val),pos)
end

function _insertafter(t,pos,n_insert,val,i)
  if i <= pos
    return t[i]
  elseif i > pos+n_insert
    return t[i-n_insert]
  end
  return val[i-pos]
end

"""
    insertafter(t, val, pos)

Insert the elements in val after the position pos
"""
function insertafter(t::NTuple{N},
                     val::NTuple{M},
                     pos::Integer) where {N,M}
  return ntuple(i -> _insertafter(t,pos,M,val,i),Val(N+M))
end

function insertafter(t::NTuple{N}, val, pos::Integer) where {N}
  return insertafter(t,tuple(val),pos)
end

"""
    isdisjoint(s1, s2)

Determine if s1 and s2 have no overlapping elements.
"""
function isdisjoint(s1,s2)
  for i1 ∈ 1:length(s1)
    for i2 ∈ 1:length(s2)
      s1[i1] == s2[i2] && return false
    end
  end
  return true
end

"""
    diff(t::Tuple)

For a tuple of length N, return a tuple of length N-1
where element i is t[i+1] - t[i].
"""
diff(t::NTuple{N}) where {N} =
  ntuple(i -> t[i+1] - t[i], Val(N-1))

function count_unique(labelsT1,labelsT2)
  count = 0
  for l1 ∈ labelsT1
    l1 ∉ labelsT2 && (count += 1)
  end
  return count
end

function count_common(labelsT1,labelsT2)
  count = 0
  for l1 ∈ labelsT1
    l1 ∈ labelsT2 && (count += 1)
  end
  return count
end

function intersect_positions(labelsT1,labelsT2)
  for i1 = 1:length(labelsT1)
    for i2 = 1:length(labelsT2)
      if labelsT1[i1] == labelsT2[i2]
        return i1,i2
      end
    end
  end
  return nothing
end

function is_replacement(labelsT1,labelsT2)
  return count_unique(labelsT1,labelsT2) == 1 &&
         count_common(labelsT1,labelsT2) == 1
end

function is_combiner(labelsT1,labelsT2)
  return count_unique(labelsT1,labelsT2) == 1 &&
         count_common(labelsT1,labelsT2) > 1
end

function is_uncombiner(labelsT1,labelsT2)
  return count_unique(labelsT1,labelsT2) > 1 &&
         count_common(labelsT1,labelsT2) == 1
end


text:
module JUGRNModelGenerator

    # include -
    include("Include.jl")

    # export functions -
    export check
    export build_julia_model_object
    export read_model_document
    export make_julia_model

    # export types -
    export VLJuliaModelObject
    export VLResult

end # module


text:
################################################################################
#
#   (8,3)c
#
################################################################################

# Implementation
# - implementation 1
# - labels <: Any
function getUnitcell_8_3_c(
			unitcell_type  :: Type{U},
			implementation :: Val{1}
		) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

	# return a new Unitcell
	return newUnitcell(
		# Type of the unitcell
		U,
		# lattice vectors
		Vector{Float64}[
			Float64[1., 0., 0.],
			Float64[-1/2., sqrt(3)/2., 0.],
			Float64[0., 0., 2/5.],
		],
		# sites
		S[
			newSite(S, Float64[-1/5., 4/(5*sqrt(3)), 1/10.], getDefaultLabelN(LS,1)),
			newSite(S, Float64[ 0., 7/(5*sqrt(3)), 1/10.], getDefaultLabelN(LS,2)),
			newSite(S, Float64[ 1/5., 4/(5*sqrt(3)), 1/10.], getDefaultLabelN(LS,3)),
			newSite(S, Float64[ 1/2., 1/(2*sqrt(3)), 3/10.], getDefaultLabelN(LS,4)),
			newSite(S, Float64[ 0., 1/sqrt(3), 1/10.], getDefaultLabelN(LS,5)),
			newSite(S, Float64[ 3/10., 7/(10*sqrt(3)), 3/10.], getDefaultLabelN(LS,6)),
			newSite(S, Float64[ 1/2., 1/(10*sqrt(3)), 3/10.], getDefaultLabelN(LS,7)),
			newSite(S, Float64[ 7/10., 7/(10*sqrt(3)), 3/10.], getDefaultLabelN(LS,8)),
		],
		# bonds
		B[
			newBond(B, 1, 5, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 2, 5, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 5, 3, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 3, 6, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 6, 4, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 4, 7, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 4, 8, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 5, 1, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 5, 2, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 3, 5, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 6, 3, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 4, 6, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 7, 4, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 8, 4, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 8, 1, getDefaultLabel(LB), (1, 0, 0)),
			newBond(B, 1, 8, getDefaultLabel(LB), (-1, 0, 0)),
			newBond(B, 8, 1, getDefaultLabel(LB), (1, 0, 1)),
			newBond(B, 1, 8, getDefaultLabel(LB), (-1, 0, -1)),
			newBond(B, 2, 7, getDefaultLabel(LB), (0, 1, 0)),
			newBond(B, 7, 2, getDefaultLabel(LB), (0, -1, 0)),
			newBond(B, 2, 7, getDefaultLabel(LB), (0, 1, -1)),
			newBond(B, 7, 2, getDefaultLabel(LB), (0, -1, 1)),
			newBond(B, 3, 6, getDefaultLabel(LB), (0, 0, -1)),
			newBond(B, 6, 3, getDefaultLabel(LB), (0, 0, 1)),
		]
	)
end

# Implementation
# - implementation 4 (Kitaev)
# - labels <: Any
function getUnitcell_8_3_c(
			unitcell_type  :: Type{U},
			implementation :: Val{4}
		) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

	# return a new Unitcell
	return newUnitcell(
		# Type of the unitcell
		U,
		# lattice vectors
		Vector{Float64}[
			Float64[1., 0., 0.],
			Float64[-1/2., sqrt(3)/2., 0.],
			Float64[0., 0., 2/5.],
		],
		# sites
		S[
			newSite(S, Float64[-1/5., 4/(5*sqrt(3)), 1/10.], getDefaultLabelN(LS,1)),
			newSite(S, Float64[ 0., 7/(5*sqrt(3)), 1/10.], getDefaultLabelN(LS,2)),
			newSite(S, Float64[ 1/5., 4/(5*sqrt(3)), 1/10.], getDefaultLabelN(LS,3)),
			newSite(S, Float64[ 1/2., 1/(2*sqrt(3)), 3/10.], getDefaultLabelN(LS,4)),
			newSite(S, Float64[ 0., 1/sqrt(3), 1/10.], getDefaultLabelN(LS,5)),
			newSite(S, Float64[ 3/10., 7/(10*sqrt(3)), 3/10.], getDefaultLabelN(LS,6)),
			newSite(S, Float64[ 1/2., 1/(10*sqrt(3)), 3/10.], getDefaultLabelN(LS,7)),
			newSite(S, Float64[ 7/10., 7/(10*sqrt(3)), 3/10.], getDefaultLabelN(LS,8)),
		],
		# bonds
		B[
			newBond(B, 1, 5, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 2, 5, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 5, 3, getDefaultLabelY(LB), (0, 0, 0)),
			newBond(B, 3, 6, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 6, 4, getDefaultLabelY(LB), (0, 0, 0)),
			newBond(B, 4, 7, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 4, 8, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 5, 1, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 5, 2, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 3, 5, getDefaultLabelY(LB), (0, 0, 0)),
			newBond(B, 6, 3, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 4, 6, getDefaultLabelY(LB), (0, 0, 0)),
			newBond(B, 7, 4, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 8, 4, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 8, 1, getDefaultLabelY(LB), (1, 0, 0)),
			newBond(B, 1, 8, getDefaultLabelY(LB), (-1, 0, 0)),
			newBond(B, 8, 1, getDefaultLabelZ(LB), (1, 0, 1)),
			newBond(B, 1, 8, getDefaultLabelZ(LB), (-1, 0, -1)),
			newBond(B, 2, 7, getDefaultLabelX(LB), (0, 1, 0)),
			newBond(B, 7, 2, getDefaultLabelX(LB), (0, -1, 0)),
			newBond(B, 2, 7, getDefaultLabelY(LB), (0, 1, -1)),
			newBond(B, 7, 2, getDefaultLabelY(LB), (0, -1, 1)),
			newBond(B, 3, 6, getDefaultLabelX(LB), (0, 0, -1)),
			newBond(B, 6, 3, getDefaultLabelX(LB), (0, 0, 1)),
		]
	)
end


text:
##################################################
#   Override default handling (default = inline) #
##################################################

get_latex_function(r::DataFrames.DataFrame) = mdtable
get_md_function(args::DataFrames.DataFrame) = mdtable

###############################################
#         Overload environment functions      #
###############################################

function mdtable(d::DataFrames.DataFrame; kwargs...)
    body = convert(Matrix, d)
    head = propertynames(d)
    mdtable(body; head=head, kwargs...)
end

function latextabular(d::DataFrames.DataFrame; kwargs...)
    body = convert(Matrix, d)
    head = propertynames(d)
    latextabular(body; head=head, kwargs...)
end

function latexarray(d::DataFrames.DataFrame; kwargs...)
    body = convert(Matrix, d)
    head = permutedims(propertynames(d))
    result = vcat(head, body)
    latexarray(result; kwargs...)
end


text:
using IterativeSolvers
using Base.Test

@testset "SVD Lanczos" begin

srand(1234567)

#Thick restart methods
@testset "Thick restart with method=$method" for method in (:ritz, :harmonic)
    for T in (Float32, Float64)
        @testset "Diagonal Matrix{$T}" begin
            n = 30
            ns = 5
            tol = 1e-5

            A = full(Diagonal(T[1.0 : n;]))
            q = ones(T, n) / √n
            σ, L, history = svdl(A, nsv=ns, v0=q, tol=tol, reltol=tol, maxiter=n, method=method, vecs=:none, log=true)
            @test isa(history, ConvergenceHistory)
            @test norm(σ - [n : -1.0 : n - 4;]) < 5^2 * 1e-5
            @test_throws ArgumentError svdl(A, nsv=ns, v0=q, tol=tol, reltol=tol, maxiter=n, method=:fakemethod, vecs=:none)

            #Check the singular vectors also
            Σ, L = svdl(A, nsv=ns, v0=q, tol=tol, reltol=tol, maxiter=n, method=method, vecs=:both)

            #The vectors should have the structure
            # [ 0  0 ...  0 ]
            # [    ...      ]
            # [ 0  0 ... ±1 ]
            # [ 0   ...     ]
            # [ 0 ±1 ...  0 ]
            # [±1  0 ...  0 ]
            # and furthermore the signs should be aligned across Σ[:U] and Σ[:V]
            for i = 1 : 5
                Σ[:U][end + 1 - i, i] -= sign(Σ[:U][end + 1 - i, i])
            end
            @test vecnorm(Σ[:U]) < σ[1] * √tol
            for i = 1 : 5
                Σ[:Vt][i, end + 1 - i] -= sign(Σ[:Vt][i, end + 1 - i])
            end
            @test vecnorm(Σ[:U]) < σ[1] * √tol
            @test norm(σ - Σ[:S]) < 2max(tol * ns * σ[1], tol)

            #Issue #55
            let
                σ1, _ = svdl(A, nsv=1, tol=tol, reltol=tol)
                @test abs(σ[1] - σ1[1]) < 2max(tol * σ[1], tol)
            end
        end

        @testset "Rectangular Matrix{$T}" begin
            srand(1)
            m = 300
            n = 200
            k = 5
            l = 10

            A = randn(T, m, n)
            q = randn(T, n) |> x -> x / norm(x)
            σ, L = svdl(A, nsv=k, k=l, v0=q, tol=1e-5, maxiter=30, method=method)
            @test norm(σ - svdvals(A)[1 : k]) < k^2 * 1e-5
        end
    end
end
end #svdl

@testset "BrokenArrowBidiagonal" begin
    B = IterativeSolvers.BrokenArrowBidiagonal([1, 2, 3], [1, 2], Int[])
    @test full(B) == [1 0 1; 0 2 2; 0 0 3]
    @test B[3,3] == 3
    @test B[2,3] == 2
    @test B[3,2] == 0
    @test B[1,3] == 1
    @test size(B) == (3,3)
    @test_throws ArgumentError size(B,3)
    @test_throws BoundsError B[1,5]
end


text:
# Algorithm 3.6
# Bezier Curve

function bezier(n, x, xp, xn, y, yp, yn)
	a = zeros(1, 3)
	b = zeros(1, 3)
	for i ∈ 1:(n - 1)
		a[1] = x[i]
		b[1] = y[i]
		a[2] = 3(xp[i] - x[i])
		b[2] = 3(yp[i] - y[i])
		a[3] = 3(x[i] + xn[i + 1] -2xp[i])
		b[3] = 3(y[i] + yn[i + 1] -2yp[i])
		a[4] = x[i + 1] - x[i] + 3xp[i] -3xn[i + 1]
		a[4] = y[i + 1] - y[i] + 3yp[i] -3yn[i + 1]
	end
end


text:
module TestGraphAlgorithms
using Test

using Catlab.Graphs

# Connectivity
##############

g = Graph(6)
add_edges!(g, [2,3], [4,5])
@test connected_components(g) == [[1], [2,4], [3,5], [6]]

# DAGs
######

# Topological sorting
#--------------------

# Discrete graph.
g = Graph(5)
@test topological_sort(g) == 1:5

# Path graph.
g = Graph(5)
add_edges!(g, 5:-1:2, 4:-1:1)
@test topological_sort(g) == 5:-1:1

# Diamond.
g = Graph(4)
add_edges!(g, [1,1,2,3], [2,3,4,4])
@test topological_sort(g) == 1:4

# Multiple edges.
g = Graph(3)
add_edges!(g, [1,1,1,2], [2,2,3,3])
@test topological_sort(g) == 1:3

# Error handling: cyclic graph.
g = Graph(3)
add_edges!(g, [1,2,3], [2,3,1])
@test_throws Exception topological_sort(g)

# Transitivity
#-------------

g = Graph(3)
add_edges!(g, [1,2], [2,3])
g₀ = copy(g)
@test transitive_reduction!(g) == g₀
add_edge!(g, 1, 3)
@test transitive_reduction!(g) == g₀

g = Graph(5)
add_edges!(g, [1,1,1,1,2,3,3,4], [2,3,4,5,4,4,5,5])
transitive_reduction!(g)
@test sort!(collect(zip(src(g), tgt(g)))) == [(1,2),(1,3),(2,4),(3,4),(4,5)]

# Enumerating paths
g = Graph(3)
add_edges!(g, [1,2,2], [2,3,3])

ep = enumerate_paths(g)
@test sort(collect(zip(ep[:src],ep[:tgt], ep[:eprops]))) == [
    (1, 1, [])
    (1, 2, [1])
    (1, 3, [1, 2])
    (1, 3, [1, 3])
    (2, 2, [])
    (2, 3, [2])
    (2, 3, [3])
    (3, 3, [])]
end


text:
"""
Auto-generated dynamics for solving phase and tangent dynamics together.
"""
struct PhaseTangentDynamics{F, J}
    "A callable in `(du, u, p, t)` format (inplace version for `ODEProblem`)."
    phase_dynamics::F

    "Temporary variable to store the Jacobian of the `phase_dynamics`."
    jacobian::J

    # TODO: store ForwardDiff config here
end

function fdiff_tangent_dynamics(phase_dynamics, x0)
    dimension = size(x0, 1)
    return PhaseTangentDynamics(
        phase_dynamics,
        similar(x0, (dimension, dimension)),
    )
end

"""Co-evolve phase- and tangent-sapce dynamics."""
@inline function (dyn::PhaseTangentDynamics)(du, u, args...)
    ForwardDiff.jacobian!(
        dyn.jacobian,
        (dx, x) -> dyn.phase_dynamics(dx, x, args...),
        (@view du[:, 1]),       # phase space derivative `dx` goes here
        (@view u[:, 1]),        # current phase space state `x`
    )
    A_mul_B!((@view du[:, 2:end]), dyn.jacobian, (@view u[:, 2:end]))
end