text:
function o_piece(cur_container::AbstractContainer)
  cur_piece = Piece(
    cur_container,
    Block[],
    'o',
    "yellow",
    2,
    0,
    -1,
    -1
  )

  cur_block_coords = [
    (-1, +1),
    (-1, -1),
    (+1, -1),
    (+1, +1)
  ]

  make_blocks!(cur_piece, cur_block_coords)

  cur_piece
end


text:
"""
    AIC(sem_fit::SemFit)

Return the akaike information criterion.
"""
AIC(sem_fit::SemFit) = minus2ll(sem_fit) + 2n_par(sem_fit)

ext: jl


text:
module GenerateTests

using Compat.Test
using Documenter

@testset "Generate" begin
    mktempdir() do root
        let path = joinpath(root, "docs")
            Documenter.generate("DocumenterTestPackage", dir = path)
            @test isdir(path)
            @test isfile(joinpath(path, "mkdocs.yml"))
            @test isfile(joinpath(path, ".gitignore"))
            @test isfile(joinpath(path, "make.jl"))
            @test isdir(joinpath(path, "src"))
            @test isfile(joinpath(path, "src", "index.md"))
        end
    end

    # TODO: these tests should be reviewed. Documenter.generate() does not really
    # support Pkg3 / Julia 0.7 at the moment.
    @test_throws ErrorException Documenter.generate("Documenter")
    if VERSION < v"0.7.0-"
        @test_throws ErrorException Documenter.generate(randstring())
    else
        @test_throws MethodError Documenter.generate(randstring())
    end
end

end


text:
using DiffEqDocs, DifferentialEquations
using Base.Test

# Tests are just for docs generation


text:
abstract type EoSModel end
abstract type LatticeFluidModel <: EoSModel end
abstract type SAFTModel <: LatticeFluidModel end
abstract type GCSAFTModel <: SAFTModel end
abstract type CubicModel <: EoSModel end
abstract type ABCubicModel <: CubicModel end #cubics that have an exact polynomial form to solve Z roots, this excludes CPA
abstract type IdealModel <: EoSModel end
abstract type EmpiricHelmholtzModel <: EoSModel end
abstract type EoSParam end


text:
using Glob: glob

const NAME_REGEX = r"\n?\s*--\s*name:\s*(?<query_name>\w+[\^\*\!\<\#\$]{0,2}+)\n?"
const DOCSTRING_REGEX = r"^--\s*(?<docstring>.*)\n+"
const COMMENT_REGEX = r"\n?+\s*?--.*\n?"
const ARG_REGEX = r"(?<leading_char>[^:]):(?<symbol_name>\w+)"

"""
    $(TYPEDSIGNATURES)

Parse `queries_text` and insert extracted `SQLQuery`s into `queries`

"""
function _extract!(queries::AbstractDict, queries_text::AbstractString)

    # Find all query names in query_text
    name_matches = collect(eachmatch(NAME_REGEX, queries_text))

    for (i,nm) in enumerate(name_matches)
        operation, query_name = _interpret(nm[:query_name])

        if query_name in keys(queries)
            @warn "$(query_name) already defined and will be overwritten"
        end

        # Find start/end indices of query and extract
        istart = nm.offset + length(nm.match)
        iend = i < length(name_matches) ? name_matches[i+1].offset - 1 : length(queries_text)
        query = queries_text[istart:iend]

        # Extract the docstring if present
        docstring = ""
        doc_match = match(DOCSTRING_REGEX, query)
        while doc_match !== nothing
            docstring *= " $(string(strip(doc_match[:docstring])))"
            # docstring = doc_match === nothing ? "" : string(strip(doc_match[:docstring]))

            query = replace(strip(query), DOCSTRING_REGEX => "")
            doc_match = match(DOCSTRING_REGEX, query)
        end
        docstring = string(lstrip(docstring))

        # Strip out any comments
        query = replace(query, COMMENT_REGEX => " ")

        # Find any named arguments in the query and replace them with $1, $2 etc
        args = Symbol[]
        num_args = 0
        for m in eachmatch(ARG_REGEX, query)
            arg = Symbol(m[:symbol_name])
            occursin(":$(arg)", query) || continue
            if !(arg in args)      # argument already captured
                num_args += 1
                push!(args, arg)
            end
            query = replace(query, m.match => "$(m[:leading_char])\$$(num_args)")
        end

        # Replace newlines and extra spaces with a single space
        query = replace(query, r"\s+" => " ")

        # Insert the SQLQuery into the Dict
        queries[query_name] = SQLQuery(query_name, query, docstring, args, operation)
    end
    return nothing
end

"""
    $(TYPEDSIGNATURES)

Read files(s) at `path` and insert extracted `SQLQuery`s into
`queries`. If `path` is a directory then include files matching
`glob_pattern`.

"""
function load_path!(queries::AbstractDict, path; glob_pattern="*.sql")
    ispath(path) || throw(ArgumentError("$path does not exist"))
    files = isfile(path) ? [path] : glob(glob_pattern, path)
    for file in files
       _extract!(queries, read(file, String))
    end
    return nothing
end

"""
    $(TYPEDSIGNATURES)

Read files(s) at `path` and return extracted `SQLQuery`s. If `path` is
a directory then include files matching `glob_pattern`.

"""
function load_path(path; glob_pattern="*.sql")
    queries = Dict{Symbol, SQLQuery}()
    load_path!(queries, path; glob_pattern=glob_pattern)
    return queries
end

"""
    $(TYPEDSIGNATURES)

Parse `query_text` and insert extracted `SQLQuery`s into `queries`

"""
load_string!(queries::AbstractDict, query_text::AbstractString) = _extract!(queries, query_text)

"""
    $(TYPEDSIGNATURES)

Parse `query_text` and return extracted `SQLQuery`s

"""
function load_string(query_text::AbstractString)
    queries=Dict{Symbol, SQLQuery}()
    load_string!(queries, query_text)
    return queries
end


text:
include("src/TALib.jl")

using TALib
using TALib: TA_MAType_SMA

using DataFrames

print("Initialize: ")
retCode = Initialize()
println(retCode)

for f in [GetVersionString, GetVersionMajor, GetVersionMinor, GetVersionBuild, GetVersionDate, GetVersionTime]
    print(string(f) * ": ")
    println(f())
end

s_xml = FunctionDescriptionXML()
#println(s_xml)

println("Input")
inReal = [0, pi/2, pi, 3pi/2, 0, pi/2, pi, 3pi/2]
println(inReal)

println("COS")
outReal = COS(inReal)
println(outReal)

println("ACOS")
outReal = ACOS(outReal)
println(outReal)

println("SIN")
outReal = SIN(inReal)
println(outReal)

println("ASIN")
outReal = ASIN(outReal)
println(outReal)

println("TAN")
outReal = TAN(inReal)
println(outReal)

println("ATAN")
outReal = ATAN(outReal)
println(outReal)

filename = "test/ford_2012.csv"
println("Read '$filename'")
dfOHLCV = readtable(filename)
dfOHLCV[:Date] = Date(dfOHLCV[:Date])
println(dfOHLCV)
opn = Array(dfOHLCV[:Open])
hig = Array(dfOHLCV[:High])
low = Array(dfOHLCV[:Low])
cls = Array(dfOHLCV[:Close])
price = Array(dfOHLCV[:Close])
vol = Array(dfOHLCV[:Volume])
println(price)

println("MA")
#indic = MA(price)
indic = MA(price, time_period=30, ma_type=TA_MAType_SMA)
println(indic)

println(MA(dfOHLCV, price=:Close))

using PyPlot
#plot(dfOHLCV[:Date], dfOHLCV[:Close])
plot(dfOHLCV[:Date], dfOHLCV[:Close], dfOHLCV[:Date], indic)
#=
plot(dfOHLCV[:Date], dfOHLCV[:Open],
     dfOHLCV[:Date], dfOHLCV[:High],
     dfOHLCV[:Date], dfOHLCV[:Low],
     dfOHLCV[:Date], dfOHLCV[:Close],
)
=#

println("BBANDS")
println(price)
#indic = BBANDS(price)
#indic = BBANDS(price, time_period=30, deviations_up=2.0, deviations_down=2.0, ma_type=TA_MAType_SMA)
#println(indic)
indic = BBANDS(dfOHLCV, price=:Close)
println(indic)
#=
plot(dfOHLCV[:Date], dfOHLCV[:Close],
    dfOHLCV[:Date], indic[:UpperBand],
    dfOHLCV[:Date], indic[:MiddleBand],
    dfOHLCV[:Date], indic[:LowerBand],
)
=#

pause(10)

#angles = readdlm("test/angles.csv")
#angles = reshape(angles, length(angles))  # convert Array{Float64,2} to Array{Float64,1}

print("Shutdown: ")
retCode = Shutdown()
println(retCode)


text:
module DataGen

using DataFrames, Distributions, StatsBase

#methods
export gen_ts

# includes
include("tsr.jl")

end


text:
module DLX

export DancingLink, is_solved

"""
    DancingLink
"""
struct DancingLink
    prev::Union{Missing, DancingLink}
    board::Array{Integer, 2}
    members::Array
    step::Int
    slots::Array{Tuple, 1}
    selection::Int
end
function DancingLink(world::Tuple, members)
    step = 0
    board = zeros(Int16, world)
    slots = Tuple[]

    if !isempty(members)
        if sum(el -> el[1]*el[2], members) <= world[1] * world[2]
            # TODO: 못 넣을 때 이거 순서 바꾸면 가능해짐...
            sort!(members; by = el -> el[1]*el[2], rev = true)

            slots = search_slot(board, members[1])
            if !isempty(slots)
                step = 1
                board = fillboard(board, slots[1], members[step])
            end
        end
    end
    # return setp=0 if it's not feasible
    me = DancingLink(missing, board, members, step, slots, 1)
    solve(me)
end
"""
    DancingLink(worlds::Vector, members)

제일 큰 world에 전부 구겨넣기, 남는거만 다음 world로 전달
"""
function DancingLink(worlds::Vector, members)
    sort!(worlds; by = el -> el[1] * el[2], rev = true)

    solutions = []
    for w in worlds
        i = length(members)

        while true
            if i == 0
                push!(solutions, DancingLink(w, members))
                break
            end
            sol = DancingLink(w, members[1:i])
            if is_solved(sol)
                push!(solutions, sol)
                if i < length(members)
                    members = members[i+1:end]
                else
                    members = []
                end
                break
            else
                i -= 1
            end

        end
    end
    if length(members) > 0
        @warn "DancingLink - $(members) cannot fit into any world"
    end
    return solutions
end

function solve(me::DancingLink)
    if is_solved(me)
        me
    else
        if is_windable(me)
            solve(wind(me))
        elseif is_forwardable(me)
            solve(forward(me))
        elseif is_unwindable(me)
            solve(unwind(me))
        else
            me
        end
    end
end

function wind(me::DancingLink)
    board = me.board
    members = me.members

    step = me.step+1
    slots = search_slot(board, members[step])

    board = fillboard(board, slots[1], members[step])
    DancingLink(me, board, members, step, slots, 1)
end
function forward(me::DancingLink)
    selection = me.selection +1

    board = replace(me.board, maximum(me.board) => 0)
    board = fillboard(board, me.slots[selection], me.members[1])

    return DancingLink(me.prev, board, me.members, me.step, me.slots, selection)
end
function unwind(me::DancingLink)
    forward(me.prev)
end

function is_solved(x::DancingLink)
    x.step == length(x.members)
end
function is_windable(me::DancingLink)
    if me.step < length(me.members)
        slots = search_slot(me.board, me.members[me.step+1])
        if !isempty(slots)
            return true
        end
    end
    return false
end
is_forwardable(x::Missing) = false
function is_forwardable(x::DancingLink)
    x.selection < length(x.slots)
end
function is_unwindable(x::DancingLink)
    is_forwardable(x.prev)

end

function Base.show(io::IO, x::DancingLink)
    summary(io, x)
    print(io, " - ")
    display(x.board)
    print(io, "Members:[")
    for (i, el) in enumerate(x.members)
        cor = i <= x.step ? :green : :normal
        printstyled(io, el; color = cor)

        i == length(x.members) ? print(io, "]") : print(io, ",")
    end
end

#
function fillboard(board, cod, target)
    b = copy(board)
    b[cod[1], cod[2]] .= maximum(board) +1

    return b
end

"""
    search_slot(board, member_size)

if value iszero and top-left value is nonzero
"""
function search_slot(board, memeger::Tuple)
    a = _search(board, memeger)
    if memeger[1] != memeger[2]
        a = append!(a, _search(board, reverse(memeger)))
    end
    return a
end

function _search(board, member::Tuple)
    # Row=1이면 왼쪽이 notzero Col은 무관
    # Col=1이면 위가 notzero Row는 무관
    # 그외에는 왼쪽과 위가 모두 notzero
    cods = Tuple[]
    board_size = size(board)

    for row in 1:board_size[1]
        for col in 1:board_size[2]
            if iszero(board[row, col])
                rg_row = row:(row + member[1] - 1)
                rg_col = col:(col + member[2] - 1)
                if last(rg_row) <= board_size[1] && last(rg_col) <= board_size[2]
                    # 또한 이 영역이 전부 0이어야 한다
                    if iszero(sum(board[rg_row, rg_col]))
                        push!(cods, (rg_row, rg_col))
                    end
                end
            end
        end
    end
    return cods
end

end # module


text:
using FromFile
@from "../types.jl" import MyInt
@from "ops.jl" import Square

int_unwrap(x::MyInt)::Int = x.value
int_square_unwrap(x::MyInt)::Int = int_unwrap(Square.int_square(x))


text:
# image_geom.jl

using MIRT: ImageGeom, image_geom, cbct
using MIRT: image_geom_ellipse
using MIRTjim: jim
using FillArrays: Trues
using Unitful: m

using Test: @test, @testset, @test_throws, @inferred

@testset "construct" begin
    ig = @inferred ImageGeom()
    @test ig isa ImageGeom

    arg1 = ((2,), (3,), (0,))
    arg2 = ((2,3), (3m,4.), (0,1))
    arg3 = ((2,3,4), (3m,4.,1//2), (0,2//3,2.))

    ig1 = @inferred ImageGeom(arg1..., trues(arg1[1]...))
    ig2 = @inferred ImageGeom(arg2..., trues(arg2[1]...))
    ig3 = @inferred ImageGeom(arg3..., trues(arg3[1]...))
    @test ig3 isa ImageGeom

    ig1 = @inferred ImageGeom(arg1...)
    ig2 = @inferred ImageGeom(arg2...)
    ig3 = @inferred ImageGeom(arg3...)
    @test ig3 isa ImageGeom

    # _zero tests
    @test ig1.dz === zero(Int)
    @test ig2.dz === zero(Int32) # alert!
    @test ImageGeom((2,), (3.0,), (0,)).dz === zero(Float64)
    @test ImageGeom((2,3), (3.0,4.0f0), (0,0)).dz === zero(Float32)
    @test ImageGeom((2,3), (3.0m,4.0m), (0,0)).dz === zero(0.0m)
end

function image_geom_test2(ig::ImageGeom)
	# test 2D functions provided by the constructor
	ig.dim
	ig.dims
	ig.deltas
	ig.offsets
	ig.x
	ig.y
	ig.wx
	ig.wy
	ig.xg
	ig.yg
	ig.fovs
	ig.np
	ig.mask_outline
	ig.ones
	ig.zeros
	ig.u
	ig.v
	ig.ug
	ig.vg
	ig.fg
	@test ig.shape(vec(ig.ones)) == ig.ones
    # todo: collect due to https://github.com/JuliaArrays/FillArrays.jl/pull/149
	@test ig.embed(collect(ig.ones)[ig.mask]) == ig.mask
	@test ig.maskit(ig.ones) == ones(ig.np)
#	@inferred # todo
	ig.unitv()
	ig.unitv(j=4)
	ig.unitv(i=ntuple(i->1, length(ig.dim)))
	ig.unitv(c=ntuple(i->0, length(ig.dim)))
#= todo-i: why do these fail?
	@inferred image_geom_ellipse(8, 10, 1, 2)
	@inferred ig.circ()
	@inferred ig.plot()
	@inferred ig.unitv()
=#
	ig.circ()
	ig.plot(jim)
#	isinteractive() && gui()
	ig.unitv()
#= todo-i:
    @inferred ig.down(2)
    @inferred ig.over(2)
=#
	@test ig.down(2) isa ImageGeom
	@test ig.over(2) isa ImageGeom
	image_geom_ellipse(8, 10, 1, 2)
	true
end

function image_geom_test2()
	#@inferred
	image_geom(nx=16, dx=2, offsets=:dsp, mask=:all_but_edge_xy)
	@test_throws String image_geom(nx=16, dx=1, offsets=:bad)
	@test_throws String image_geom(nx=16, dx=1, mask=:bad) # mask type
	@test_throws String image_geom(nx=16, dx=1, mask=trues(2,2)) # mask size
	ig = image_geom(nx=16, dx=2)
	@test ig.mask == Trues(ig.dims)
	show(isinteractive() ? stdout : devnull, ig)
	show(isinteractive() ? stdout : devnull, MIME("text/plain"), ig)
	image_geom_test2(ig)
	ig = image_geom(nx=16, dx=2, mask=:circ)
	ig.over(2)
	ig.down(3) # test both even and non-even factors
	ig.help
	true
end

function image_geom_test3(ig::ImageGeom)
	@test image_geom_test2(ig)
	ig.wz
	ig.zg
	ig.mask_or
	@inferred ig.expand_nz(2)
	@inferred cbct(ig)
	@test_throws String image_geom(nx=1, dx=1, nz=2, offset_z=-1, offsets=:dsp)
	true
end

function image_geom_test3()
	#@inferred
	image_geom(nx=16, nz=4, dx=2, zfov=1)
	ig = image_geom(nx=16, nz=4, dx=2, dz=3)
	image_geom_test3(ig)
	true
end

@testset "2d" begin
	@test all(ImageGeom((2,), (3,), (0,)).mask)
	@test image_geom_test2()
end

@testset "3d" begin
	@test image_geom_test3()
end


text:
# Julia wrapper for header: /usr/include/scip/scip_heur.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeHeur(scip, name, desc, dispchar, priority, freq, freqofs, maxdepth, timingmask, usessubscip, heurcopy, heurfree, heurinit, heurexit, heurinitsol, heurexitsol, heurexec, heurdata)
    ccall((:SCIPincludeHeur, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, UInt8, Cint, Cint, Cint, Cint, SCIP_HEURTIMING, UInt32, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_HEURDATA}), scip, name, desc, dispchar, priority, freq, freqofs, maxdepth, timingmask, usessubscip, heurcopy, heurfree, heurinit, heurexit, heurinitsol, heurexitsol, heurexec, heurdata)
end

function SCIPincludeHeurBasic(scip, heur, name, desc, dispchar, priority, freq, freqofs, maxdepth, timingmask, usessubscip, heurexec, heurdata)
    ccall((:SCIPincludeHeurBasic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_HEUR}}, Cstring, Cstring, UInt8, Cint, Cint, Cint, Cint, SCIP_HEURTIMING, UInt32, Ptr{Cvoid}, Ptr{SCIP_HEURDATA}), scip, heur, name, desc, dispchar, priority, freq, freqofs, maxdepth, timingmask, usessubscip, heurexec, heurdata)
end

function SCIPsetHeurCopy(scip, heur, heurcopy)
    ccall((:SCIPsetHeurCopy, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_HEUR}, Ptr{Cvoid}), scip, heur, heurcopy)
end

function SCIPsetHeurFree(scip, heur, heurfree)
    ccall((:SCIPsetHeurFree, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_HEUR}, Ptr{Cvoid}), scip, heur, heurfree)
end

function SCIPsetHeurInit(scip, heur, heurinit)
    ccall((:SCIPsetHeurInit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_HEUR}, Ptr{Cvoid}), scip, heur, heurinit)
end

function SCIPsetHeurExit(scip, heur, heurexit)
    ccall((:SCIPsetHeurExit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_HEUR}, Ptr{Cvoid}), scip, heur, heurexit)
end

function SCIPsetHeurInitsol(scip, heur, heurinitsol)
    ccall((:SCIPsetHeurInitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_HEUR}, Ptr{Cvoid}), scip, heur, heurinitsol)
end

function SCIPsetHeurExitsol(scip, heur, heurexitsol)
    ccall((:SCIPsetHeurExitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_HEUR}, Ptr{Cvoid}), scip, heur, heurexitsol)
end

function SCIPfindHeur(scip, name)
    ccall((:SCIPfindHeur, libscip), Ptr{SCIP_HEUR}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPgetHeurs(scip)
    ccall((:SCIPgetHeurs, libscip), Ptr{Ptr{SCIP_HEUR}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNHeurs(scip)
    ccall((:SCIPgetNHeurs, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPsetHeurPriority(scip, heur, priority)
    ccall((:SCIPsetHeurPriority, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_HEUR}, Cint), scip, heur, priority)
end

function SCIPcreateDiveset(scip, diveset, heur, name, minreldepth, maxreldepth, maxlpiterquot, maxdiveubquot, maxdiveavgquot, maxdiveubquotnosol, maxdiveavgquotnosol, lpresolvedomchgquot, lpsolvefreq, maxlpiterofs, initialseed, backtrack, onlylpbranchcands, specificsos1score, divesetgetscore)
    ccall((:SCIPcreateDiveset, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_DIVESET}}, Ptr{SCIP_HEUR}, Cstring, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cint, Cint, UInt32, UInt32, UInt32, UInt32, Ptr{Cvoid}), scip, diveset, heur, name, minreldepth, maxreldepth, maxlpiterquot, maxdiveubquot, maxdiveavgquot, maxdiveubquotnosol, maxdiveavgquotnosol, lpresolvedomchgquot, lpsolvefreq, maxlpiterofs, initialseed, backtrack, onlylpbranchcands, specificsos1score, divesetgetscore)
end


text:
sts = SegmentedTimeSpan(
    start_time = 0.0,
    sampling_time = 1.0,
    segment_duration = 3600,
    num_of_segments = 24,
)

@test length(sts) == 24
@test minimum(sts[1]) == 0.0
@test maximum(sts[1]) == 3600.0 - 1.0

for idx in eachindex(sts)
    @test length(sts[idx]) == 3600
end

for cur_time_span in sts
    @test length(cur_time_span) == 3600
end

@test maximum(sts[end]) == 86400.0 - 1.0

sts = SegmentedTimeSpan(
    start_time = 0.0,
    sampling_time = 1,
    segment_duration = 3,
    num_of_segments = 4,
)

sequence = []
for cur_segment in sts
    append!(sequence, collect(cur_segment))
end

@test sequence == collect(0.0:1.0:11.0)


text:
using Test

@testset "Cameras" begin
    include("acquisition.jl")
    include("iteration.jl")
    include("acquired_image.jl")
end


text:
__precompile__(true)
module PGFMapPlots

using WorldMapData
using PGFPlots

"""
Extract lon/lat points from polygon shapes.
"""
function extract_points(polygon)
    n_points = length(polygon.points)
    lon = zeros(Float64, n_points)
    lat = zeros(Float64, n_points)

    for (idx, pnt) in enumerate(polygon.points)
        lon[idx] = polygon.points[idx].x
        lat[idx] = polygon.points[idx].y
    end

    return lon, lat
end

"""
Check if any point in the given polygon is within the lat/lon bounding area. If
it is return **true** return **false** otherwise
"""
function any_in_bounds(polygon; lblon::Real=-180.0, ublon::Real=180.0, lblat::Real=-90.0, ublat::Real=90.0)
    # Check if all points of polygon are within the given area
    if any(p -> (lblon <= p.x <= ublon) && (lblat <= p.y <= ublat), polygon.points)
        return true
    else
        return false
    end
end

export plotmap
"""
Generate a map plot.

Arguments:
- `llcrnrlon::Real` Lower left corner longitude [deg]
- `llcrnrlat::Real` Lower left corner latitude [deg]
- `urcrnrlon::Real` Upper right corner longitude [deg]
- `urcrnrlat::Real` Upper right corner latitude [deg]

Returns:
- `ax::Axis` Axis
"""
function plotmap(;llcrnrlon::Real=-180.0, llcrnrlat::Real=-90.0,
                urcrnrlon::Real=180.0, urcrnrlat::Real=90.0,
                landresolution::Union{Symbol, Nothing}=nothing,
                landlevel::Union{Integer, Nothing}=1,
                landstyle::String="no marks, black, solid",
                borderresolution::Union{Symbol, Nothing}=nothing,
                borderlevel::Union{Integer, Nothing}=1,
                borderstyle::String="no marks, black, solid",
                riverresolution::Union{Symbol, Nothing}=nothing,
                riverlevel::Union{Integer, Nothing}=1,
                riverstyle::String="no marks, black, solid",
                oceanstyle="",
                style::String="")

    lines = PGFPlots.Plots.Linear[]

    # Plot ocean points
    if oceanstyle != ""
        # Hard code land points
        x = Float64[llcrnrlon, llcrnrlon, urcrnrlon, urcrnrlon, llcrnrlon]
        y = Float64[llcrnrlat, urcrnrlat, urcrnrlat, llcrnrlat, llcrnrlat]
        push!(lines, PGFPlots.Plots.Linear(x, y, style=oceanstyle))
    end

    # Plot land shapes
    if landresolution != nothing

        # Load land data
        for lvl in 1:landlevel
            if (landresolution == :CRUDE) && (lvl == 4)
                # Skip plotting this combination since the data
                # doesn't exist in the GSHHS data
                continue
            end

            land_data = load_GSHHS_data(resolution=landresolution, level=lvl)

            for ply in land_data.shapes
                if any_in_bounds(ply, lblon=llcrnrlon, ublon=urcrnrlon, lblat=llcrnrlat, ublat=urcrnrlat)
                    push!(lines, PGFPlots.Plots.Linear(extract_points(ply)..., style=landstyle))
                end
            end
        end

        # Load antartica data
        ant_data = load_GSHHS_data(resolution=landresolution, level=6)

        for ply in ant_data.shapes
            if any_in_bounds(ply, lblon=llcrnrlon, ublon=urcrnrlon, lblat=llcrnrlat, ublat=urcrnrlat)
                push!(lines, PGFPlots.Plots.Linear(extract_points(ply)..., style=landstyle))
            end
        end
    end

    # Plot borders
    if borderresolution != nothing

        # Load land data
        for lvl in 1:borderlevel
            border_data = load_WDBII_border_data(resolution=borderresolution, level=lvl)

            for ply in border_data.shapes
                if any_in_bounds(ply, lblon=llcrnrlon, ublon=urcrnrlon, lblat=llcrnrlat, ublat=urcrnrlat)
                    push!(lines, PGFPlots.Plots.Linear(extract_points(ply)..., style=borderstyle))
                end
            end
        end
    end

    # Plot rivers
    if riverresolution != nothing

        # Load land data
        for lvl in 1:riverlevel
            river_data = load_WDBII_river_data(resolution=riverresolution, level=lvl)

            for ply in river_data.shapes
                if any_in_bounds(ply, lblon=llcrnrlon, ublon=urcrnrlon, lblat=llcrnrlat, ublat=urcrnrlat)
                    push!(lines, PGFPlots.Plots.Linear(extract_points(ply)..., style=riverstyle))
                end
            end
        end
    end

    # Set style if none is provided
    if style == ""
        style = "width=25cm,height=12.5cm"
    end

    return PGFPlots.Axis(lines, xmin=llcrnrlon, xmax=urcrnrlon, ymin=llcrnrlat, ymax=urcrnrlat, style=style)
end

end # module


text:
#!/usr/bin/julia

# Trizen
# Date: 19 August 2016
# https://github.com/trizen

# Count the number of partitions of n, using a recursive relation.

# See also: https://oeis.org/A000041
#           https://en.wikipedia.org/wiki/Partition_(number_theory)

function partitions_count(n::Int64, cache::Dict{Int,Int})

    n <= 1 && return n

    if haskey(cache, n)
        return cache[n]
    end

    sum_1 = 0
    for i in 1:Int64(floor((sqrt(24*n + 1) + 1)/6))
        sum_1 += ((-1)^(i-1) * partitions_count(n - div(i*(3*i - 1), 2), cache))
    end

    sum_2 = 0
    for i in 1:Int64(ceil((sqrt(24*n + 1) - 7)/6))
        sum_2 += ((-1)^(i-1) * partitions_count(n - div((-i) * (-3*i - 1), 2), cache))
    end

    x = (sum_1 + sum_2)
    cache[n] = x
    x
end

println("P(200) = ", partitions_count(200+1, Dict{Int64, Int64}()))        # 3972999029388

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts


text:
@testset "slfm" begin
    rng = MersenneTwister(123)
    FDM = FiniteDifferences.central_fdm(5, 1)
    N = 6
    in_dim = 3
    out_dim = 4
    x1IO = KernelFunctions.MOInputIsotopicByOutputs(
        [rand(rng, in_dim) for _ in 1:N], out_dim
    )
    x2IO = KernelFunctions.MOInputIsotopicByOutputs(
        [rand(rng, in_dim) for _ in 1:N], out_dim
    )
    x3IO = KernelFunctions.MOInputIsotopicByOutputs(
        [rand(rng, in_dim) for _ in 1:div(N, 2)], out_dim
    )
    x1IO = MOInput([rand(rng, in_dim) for _ in 1:N], out_dim)
    x2IO = MOInput([rand(rng, in_dim) for _ in 1:N], out_dim)

    k = LatentFactorMOKernel(
        [Matern32Kernel(), SqExponentialKernel(), FBMKernel()],
        IndependentMOKernel(GaussianKernel()),
        rand(rng, out_dim, 3),
    )
    @test k isa LatentFactorMOKernel
    @test k isa MOKernel
    @test k isa Kernel
    @test k(x1IO[1], x2IO[1]) isa Real

    @test kernelmatrix(k, x1IO, x2IO) ≈ kernelmatrix(k, collect(x1IO), collect(x2IO))

    TestUtils.test_interface(k, x1IO, x2IO, x3IO)

    x1IF = KernelFunctions.MOInputIsotopicByFeatures(x1IO.x, out_dim)
    x2IF = KernelFunctions.MOInputIsotopicByFeatures(x2IO.x, out_dim)
    x3IF = KernelFunctions.MOInputIsotopicByFeatures(x3IO.x, out_dim)

    TestUtils.test_interface(k, x1IF, x2IF, x3IF)

    @test string(k) == "Semi-parametric Latent Factor Multi-Output Kernel"
    @test repr("text/plain", k) == (
        "Semi-parametric Latent Factor Multi-Output Kernel\n" *
        "\tgᵢ: Matern 3/2 Kernel (metric = Euclidean(0.0))\n" *
        "\t\tSquared Exponential Kernel (metric = Euclidean(0.0))\n" *
        "\t\tFractional Brownian Motion Kernel (h = 0.5)\n" *
        "\teᵢ: Independent Multi-Output Kernel\n" *
        "\tSquared Exponential Kernel (metric = Euclidean(0.0))"
    )

    # AD test
    function test_slfm(A::AbstractMatrix, x1IO, x2IO)
        k = LatentFactorMOKernel(
            [Matern32Kernel(), SqExponentialKernel(), FBMKernel()],
            IndependentMOKernel(GaussianKernel()),
            A,
        )
        return k((x1IO, 1), (x2IO, 1))
    end

    k = LatentFactorMOKernel(
        [SqExponentialKernel(), SqExponentialKernel(), SqExponentialKernel()],
        IndependentMOKernel(GaussianKernel()),
        rand(rng, out_dim, 3),
    )

    test_ADs(k)
end


text:
using WeightedNetwork

using Test
using SparseArrays
using Graphs
using SimpleWeightedGraphs
using DataStructures
using LinearAlgebra
using Plots
using GraphPlot
using DataStructures
using ShowGraphviz
using Cairo
using Fontconfig
using Random
using Tokenize

const testdir = dirname(@__FILE__)

tests = [
    "generators/test_generators",
    "generators/test_graph_utils"
]

@testset "WeightedNetwork.jl" begin
    for t in tests
        tp = joinpath(testdir, "$(t).jl")
        include(tp)
    end
end

max_stars_repo_name: tragisch/NetworkOptimization

max_issues_repo_name: tragisch/NetworkOptimization

max_forks_repo_name: tragisch/NetworkOptimization


text:
using GeoData, Test, Statistics, Dates, Plots, NCDatasets, ArchGDAL
using GeoData: name, mode, window, DiskStack
testpath = joinpath(dirname(pathof(GeoData)), "../test/")
include(joinpath(testpath, "test_utils.jl"))

maybedownload("https://raw.githubusercontent.com/rspatial/raster/master/inst/external/rlogo.grd", "rlogo.grd")
maybedownload("https://github.com/rspatial/raster/raw/master/inst/external/rlogo.gri", "rlogo.gri")
path = joinpath(testpath, "data/rlogo")
@test isfile(path * ".grd")
@test isfile(path * ".gri")

@testset "Grd array" begin
    grdarray = GRDarray(path);

    @testset "array properties" begin
        @test grdarray isa GRDarray{Float32,3}
    end

    @testset "dimensions" begin
        @test length(val(dims(dims(grdarray), Lon))) == 101
        @test ndims(grdarray) == 3
        @test dims(grdarray) isa Tuple{<:Lon,<:Lat,<:Band}
        @test refdims(grdarray) == ()
        @test bounds(grdarray) == ((0.0, 101.0), (0.0, 77.0), (1, 3))
    end

    @testset "other fields" begin
        @test missingval(grdarray) == -3.4f38
        @test metadata(grdarray) isa GRDarrayMetadata
        @test name(grdarray) == Symbol("red:green:blue")
        @test label(grdarray) == "red:green:blue"
        @test units(grdarray) == nothing
        customgrdarray = GRDarray(path; name=:test, usercrs=EPSG(4326));
        @test name(customgrdarray) == :test
        @test label(customgrdarray) == "test"
        @test usercrs(dims(customgrdarray, Lat)) == EPSG(4326)
        @test usercrs(dims(customgrdarray, Lon)) == EPSG(4326)
        @test usercrs(customgrdarray) == EPSG(4326)
        proj = ProjString("+proj=merc +datum=WGS84")
        @test crs(dims(customgrdarray, Lat)) == proj
        @test crs(dims(customgrdarray, Lon)) == proj
        @test crs(customgrdarray) == proj
    end

    @testset "getindex" begin
        @test grdarray[Band(1)] isa GeoArray{Float32,2}
        @test grdarray[Lat(1), Band(1)] isa GeoArray{Float32,1}
        @test grdarray[Lon(1), Band(1)] isa GeoArray{Float32,1}
        @test grdarray[Lon(50), Lat(30), Band(1)] == 115.0f0
        @test grdarray[1, 1, 1] == 255.0f0
        @test grdarray[Lat(At(20)), Lon(At(20)), Band(3)] == 255.0f0
        @test grdarray[Lat(Contains(60)), Lon(Contains(20)), Band(1)] == 255.0f0
    end

    # @testset "setindex" begin
    #     A = grdarray[:, :, :]
    #     temp = grdarray[1, 1, 1]
    #     println(temp)
    #     @test temp != 100.0f0
    #     grdarray[1, 1, 1] = 100.0f0
    #     grdarray[:, :, :] = 100.0f0
    #     @test grdarray[1, 1, 1] == 100.0f0
    #     grdarray[1, 1, 1] = temp
    #     @test grdarray[1, 1, 1] == temp
    #     println("sum: ", sum(A .- grdarray[:, :, :]))
    #     temp = grdarray[Lon(20), Lat(10), Band(3)]
    #     println(temp)
    #     @test temp != 200.0f0
    #     grdarray[Lon(20), Lat(10), Band(3)] = 200.0f0
    #     @test grdarray[20, 10, 3] == 200.0f0
    #     grdarray[Lon(20), Lat(10), Band(3)] = temp
    # end

    @testset "selectors" begin
        geoarray = grdarray[Lat(Contains(3)), Lon(:), Band(1)]
        @test geoarray isa GeoArray{Float32,1}
        @test grdarray[Lon(Contains(20)), Lat(Contains(10)), Band(1)] isa Float32
    end

    @testset "conversion to GeoArray" begin
        geoarray = grdarray[Lon(1:50), Lat(1:1), Band(1)]
        @test size(geoarray) == (50, 1)
        @test eltype(geoarray) <: Float32
        @time geoarray isa GeoArray{Float32,1}
        @test dims(geoarray) isa Tuple{<:Lon,Lat}
        @test refdims(geoarray) isa Tuple{<:Band}
        @test metadata(geoarray) == metadata(grdarray)
        @test missingval(geoarray) == -3.4f38
        @test name(geoarray) == Symbol("red:green:blue")
    end

    @testset "save" begin
        @testset "2d" begin
            filename2 = tempname()
            write(filename2, GRDarray, grdarray[Band(1)])
            saved = GeoArray(GRDarray(filename2))
            # 1 band is added again on save
            @test size(saved) == size(grdarray[Band(1:1)])
            @test data(saved) == data(grdarray[Band(1:1)])
        end

        @testset "3d with subset" begin
            geoarray = GeoArray(grdarray)[1:100, 1:50, 1:2]
            filename = tempname()
            write(filename, GRDarray, geoarray)
            saved = GeoArray(GRDarray(filename))
            @test size(saved) == size(geoarray)
            @test refdims(saved) == ()
            @test bounds(saved) == bounds(geoarray)
            @test size(saved) == size(geoarray)
            @test missingval(saved) === missingval(geoarray)
            @test metadata(saved) != metadata(geoarray)
            @test metadata(saved)["creator"] == "GeoData.jl"
            @test all(metadata.(dims(saved)) .== metadata.(dims(geoarray)))
            @test name(saved) == name(geoarray)
            @test all(mode.(dims(saved)) .== mode.(dims(geoarray)))
            @test dims(saved) isa typeof(dims(geoarray))
            @test all(val.(dims(saved)) .== val.(dims(geoarray)))
            @test all(mode.(dims(saved)) .== mode.(dims(geoarray)))
            @test all(metadata.(dims(saved)) .== metadata.(dims(geoarray)))
            @test dims(saved) == dims(geoarray)
            @test all(data(saved) .=== data(geoarray))
            @test saved isa typeof(geoarray)
            @test data(saved) == data(geoarray)
        end

        @testset "to netcdf" begin
            filename2 = tempname()
            write(filename2, NCDarray, grdarray[Band(1)])
            saved = GeoArray(NCDarray(filename2))
            @test size(saved) == size(grdarray[Band(1)])
            @test replace_missing(saved, missingval(grdarray)) ≈ reverse(grdarray[Band(1)]; dims=Lat)
            @test replace_missing(saved, missingval(grdarray)) ≈ reverse(grdarray[Band(1)]; dims=Lat)
            @test_broken val(dims(saved, Lon)) ≈ val(dims(grdarray, Lon)) .+ 0.5
            @test_broken val(dims(saved, Lat)) ≈ val(dims(grdarray, Lat)) .+ 0.5
            @test_broken bounds(saved, Lat) == bounds(grdarray, Lat)
            @test_broken bounds(saved, Lon) == bounds(grdarray, Lon)
        end

        @testset "to gdal" begin
            # No Band
            gdalfilename = tempname() * ".tif"
            write(gdalfilename, GDALarray, grdarray[Band(1)])
            size(grdarray)
            gdalarray = GDALarray(gdalfilename)
            # @test convert(ProjString, crs(gdalarray)) == convert(ProjString, EPSG(4326))
            @test val(dims(gdalarray, Lon)) ≈ val(dims(grdarray, Lon))
            @test reverse(val(dims(gdalarray, Lat))) ≈ val(dims(grdarray, Lat))
            @test GeoArray(gdalarray) ≈ permutedims(grdarray[Band(1)], [Lon(), Lat()])
            # 3 Bands
            gdalfilename2 = tempname() * ".tif"
            write(gdalfilename2, GDALarray, grdarray)
            gdalarray2 = GDALarray(gdalfilename2)
            @test all(GeoArray(gdalarray2) .== GeoArray(grdarray))
            @test val(dims(gdalarray2, Band)) == 1:3
        end

    end

    @testset "show" begin
        sh = sprint(show, grdarray)
        # Test but don't lock this down too much
        @test contains(sh, "GRDarray")
        @test contains(sh, "Latitude")
        @test contains(sh, "Longitude")
        @test contains(sh, "Band")
    end

    @testset "plot" begin
        p = grdarray |> plot
    end

end

@testset "Grd stack" begin
    grdstack = GRDstack((a=path, b=path))

    @testset "indexing" begin
        @test grdstack[:a][Lat(20), Lon(20), Band(3)] == 70.0f0
        @test grdstack[:a][Lat([2,3]), Lon(40), Band(2)] == [240.0f0, 246.0f0]
    end

    @testset "child array properties" begin
        @test size(grdstack[:a]) == size(GeoArray(grdstack[:a])) == (101, 77, 3)
        @test grdstack[:a] isa GeoArray{Float32,3}
    end

    @testset "window" begin
        windowedstack = GRDstack((a=path, b=path); window=(Lat(1:5), Lon(1:5), Band(1)))
        @test window(windowedstack) == (Lat(1:5), Lon(1:5), Band(1))
        windowedarray = windowedstack[:a]
        @test windowedarray isa GeoArray{Float32,2}
        @test length.(dims(windowedarray)) == (5, 5)
        @test size(windowedarray) == (5, 5)
        @test windowedarray[1:3, 2:2] == reshape([255.0f0, 255.0f0, 255.0f0], 3, 1)
        @test windowedarray[1:3, 2] == [255.0f0, 255.0f0, 255.0f0]
        @test windowedarray[1, 2] == 255.0f0
        windowedstack = GRDstack((a=path, b=path); window=(Lat(1:5), Lon(1:5), Band(1:1)))
        windowedarray = windowedstack[:b]
        @test windowedarray[1:3, 2:2, 1:1] == reshape([255.0f0, 255.0f0, 255.0f0], 3, 1, 1)
        @test windowedarray[1:3, 2:2, 1] == reshape([255.0f0, 255.0f0, 255.0f0], 3, 1)
        @test windowedarray[1:3, 2, 1] == [255.0f0, 255.0f0, 255.0f0]
        @test windowedarray[1, 2, 1] == 255.0f0
        windowedstack = GRDstack((a=path, b=path); window=(Band(1),));
        windowedarray = windowedstack[:b]
        @test windowedarray[1:3, 2:2] == reshape([255.0f0, 255.0f0, 255.0f0], 3, 1)
        @test windowedarray[1:3, 2] == [255.0f0, 255.0f0, 255.0f0]
        @test windowedarray[30, 30] == 185.0f0
        windowedarray |> plot
    end

    # Stack Constructors
    @testset "conversion to GeoStack" begin
        stack = GeoStack(grdstack)
        @test Symbol.(Tuple(keys(grdstack))) == keys(stack)
        smallstack = GeoStack(grdstack; keys=(:a,))
        @test keys(smallstack) == (:a,)
    end

    if VERSION > v"1.1-"
        @testset "copy" begin
            geoarray = zero(GeoArray(grdstack[:a]))
            copy!(geoarray, grdstack, :a)
            # First wrap with GeoArray() here or == loads from disk for each cell.
            # we need a general way of avoiding this in all disk-based sources
            @test geoarray == GeoArray(grdstack[:a])
        end
    end

    @testset "save" begin
        geoarray = GeoArray(grdstack[:b])
        filename = tempname()
        write(filename, GRDarray, grdstack)
        base, ext = splitext(filename)
        filename_b = string(base, "_b", ext)
        saved = GeoArray(GRDarray(filename_b))
        @test typeof(saved) == typeof(geoarray)
        @test data(saved) == data(geoarray)
    end

end

@testset "Grd series" begin
    series = GeoSeries([path, path], (Ti,); childtype=GRDarray, childkwargs=(usercrs=EPSG(4326), name=:test))
    @test GeoArray(series[Ti(1)]) ==
        GeoArray(GRDarray(path; usercrs=EPSG(4326), name=:test))
    stacks = [DiskStack((a=path, b=path); childtype=GRDarray, childkwargs=(usercrs=EPSG(4326), name=:test))]
    series = GeoSeries(stacks, (Ti,))
    @test series[Ti(1)][:a] ==
        GeoArray(GRDarray(path; usercrs=EPSG(4326), name=:test))
    modified_series = modify(Array, series)
    @test typeof(modified_series) <: GeoSeries{<:GeoStack{<:NamedTuple{(:a,:b),<:Tuple{<:GeoArray{Float32,3,<:Tuple,<:Tuple,<:Array{Float32,3}},Vararg}}}}
end

nothing


text:
import SparseGrids:nXtoU,cc_dM,cc_M,cc_dg,cc_bf_q,cc_bf_l
import HADSGE:State,Stochastic
function updateTsg(M)
    D         = length(M.G.L)
    mL        = maximum(M.G.L)
    J         = ones(Int, mL + 1, D)
    B         = zeros(mL + 1, D)
    T = spzeros(length(M), length(M))
    nstoc = size(M.ProbWeights, 2)
    SP = clamp!(reshape(nXtoU(M.SP, M.G.bounds), length(M), size(M.ProbWeights, 2), size(M.SP, 2)), 0, 1)
    isstoc = [isa(x, Stochastic) for x ∈ State(M.variables)]
    for i ∈ 1:length(M)
        X = SP[i,:, :]
        for k ∈ 1:nstoc
            x = X[k,:]
            for d = (1:D)[isstoc]
                for l = 1:mL + 1
                    j     = clamp(round(Int, x[d] * (cc_dM(l)) + 1 / 2), 1, cc_dM(l))
                    B[l,d]     = cc_bf_l(x[d], cc_dg(l, j), Int16(cc_M(15)))
                    J[l,d]  = j
                end
            end
            mLstoc = Int[findfirst(B[:, isstoc][:, ii]) for ii = 1:length(sum(isstoc))]

            for d = (1:D)[!isstoc]myd
                for l = 1:mL + 1
                    j     = clamp(round(Int, x[d] * (cc_dM(l)) + 1 / 2), 1, cc_dM(l))
                    # B[l,d]     = cc_bf_q(x[d],cc_dg(l,j),Int16(cc_M(l)))
                    B[l,d]     = cc_bf_l(x[d], cc_dg(l, j), Int16(cc_M(mL + D - sum(mLstoc))))
                    J[l,d]  = j
                end
            end

            res = []
            for ii = 1:size(M.G.coverings, 1)
                bv, jv = [B[M.G.coverings[ii,d],d] for d = 1:D], [J[M.G.coverings[ii,d],d] for d = 1:D]
                b  = B[M.G.coverings[ii,D],D] * B[M.G.coverings[ii,1],1]
                id1 = J[M.G.coverings[ii,D],D] - 1
                for d = D - 1:-1:2
                    b *= B[M.G.coverings[ii,d],d]
                    id1 = id1 * M.G.coverings_dM[ii,d] + (J[M.G.coverings[ii,d],d] - 1)
                end
                id1 = (J[M.G.coverings[ii,1],1] - 1) + M.G.coverings_dM[ii,1] * id1 + 1 + M.G.coveringsloc[1][ii] - 1

                M.G.coverings[ii,isstoc] == mLstoc && b > 0 && push!(res, (id1, b))
            end
            for r in res
                 T[r[1],i] += r[2] * M.ProbWeights[i,k]
            end
        end
    end
    T
end

max_stars_repo_name: ZacLN/HADSGE

max_issues_repo_name: ZacLN/HADSGE

max_forks_repo_name: ZacLN/HADSGE


text:
# Options structure
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: May 2017
#

export Options

# Object for velocity/slowness models
type Options
	space_order::Integer
	retry_n::Integer
	limit_m::Bool
	buffer_size::Real
	save_rate::Union{Real,Array{Any,1}}
	save_data_to_disk::Bool
	file_path::String
	file_name::String
	free_surface::Bool
end

"""
    Options
		space_order::Integer
        retry_n::Integer
        limit_m::Bool
        buffer_size::Real
        save_rate::Real
        save_data_to_disk::Bool
        file_path::String
        file_name::String
        free_surface::Bool

Options structure for seismic modeling.

`space_order`: order of finite difference stencils of discretized wave equations (needs to be a multiple of 4)

`retry_n`: retry modeling operations in case of worker failure up to `retry_n` times

`limit_m`: for 3D modeling, limit modeling domain to area with receivers (saves memory)

`buffer_size`: if `limit_m=true`, define buffer area on each side of modeling domain (in meters)

`save_rate`: for gradients, save background wavefield at this rate (in milliseconds)

`save_data_to_disk`: if `true`, saves shot records as separate SEG-Y files

`file_path`: path to directory where data is saved

`file_name`: shot records will be saved as specified file name plus its source coordinates

`free_surface`: use free surface if `true`

Constructor
==========

All arguments are optional keyword arguments with the following default values:

    Options(;space_order=8, retry_n=0, limit_m=false, buffer_size=1e3, save_rate=4f0, save_data_to_disk=false,
	         file_path=pwd(), file_name="shot", free_surface=false)

"""
Options(;space_order=8,retry_n=0,limit_m=false,buffer_size=1e3, save_rate=[], save_data_to_disk=false, file_path="", file_name="shot", free_surface=false) =
	Options(space_order,retry_n,limit_m,buffer_size,save_rate,save_data_to_disk,file_path,file_name, free_surface)

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
abstract type AbstractUnconstrainedSolution <: AbstractSolution end
abstract type AbstractConstrainedSolution   <: AbstractSolution end
# multi-objective solution are constrained by default
abstract type AbstractMultiObjectiveSolution <: AbstractConstrainedSolution end

mutable struct xf_indiv <: AbstractUnconstrainedSolution # Single Objective
    x::Vector{Float64}
    f::Float64
end

mutable struct xfgh_indiv <: AbstractConstrainedSolution # Single Objective Constraied
    x::Vector{Float64}
    f::Float64
    g::Vector{Float64}
    h::Vector{Float64}
    sum_violations::Float64 # ∑ max(0,g) + ∑|h|
    is_feasible::Bool

end

function xfgh_indiv(
    x::Vector{Float64},
    f::Float64,
    g::Vector{Float64},
    h::Vector{Float64};
    sum_violations = 0.0,
    ε = 0.0
)
    if sum_violations <= 0.0
        sum_violations = violationsSum(g, h; ε=ε)
    end

    xfgh_indiv(x, f, g, h, sum_violations, sum_violations == 0.0)
end

mutable struct xFgh_indiv <: AbstractMultiObjectiveSolution# Single Objective Constraied
    x::Vector{Float64}
    f::Vector{Float64}
    g::Vector{Float64}
    h::Vector{Float64}
    rank::Int
    crowding::Float64
    sum_violations::Float64 # ∑ max(0,g) + ∑|h|
    is_feasible::Bool
end

function xFgh_indiv(
    x::Vector{Float64},
    f::Vector{Float64},
    g::Vector{Float64},
    h::Vector{Float64};
    rank = 0,
    crowding = 0.0,
    sum_violations = 0.0,
    ε = 0.0
)

    if sum_violations <= 0
        sum_violations = violationsSum(g, h;ε=ε)
    end
    xFgh_indiv(x, f, g, h, Int(rank), crowding, sum_violations, sum_violations == 0.0)
end

const Solution = Union{xf_indiv, xfgh_indiv, xFgh_indiv}
const Population = Array{Solution, 1}

############################################################
# Generate solutions depending on the objective function
# output
#############################################################

function generateChild(x::Vector{Float64}, fResult::Float64;ε=0.0)
    return xf_indiv(x, fResult)
end

function generateChild(x::Vector{Float64},
        fResult::Tuple{Float64,Array{Float64,1},Array{Float64,1}};
        ε=0.0
    )
    f, g, h = fResult
    return xfgh_indiv(x, f, g, h; ε=ε)
end

function generateChild(x::Vector{Float64},
        fResult::Tuple{Array{Float64,1},Array{Float64,1},Array{Float64,1}};
        ε = 0.0
    )
    f, g, h = fResult
    return xFgh_indiv(x, f, g, h;ε=ε)
end

function generate_population(N::Int, problem;ε=0.0)
    a = view(problem.bounds, 1, :)'
    b = view(problem.bounds, 2, :)'
    D = length(a)

    X = a .+ (b - a) .* rand(N, D)

    population = [ create_solution(X[i,:], problem; ε=ε) for i in 1:N]

    return population
end

"""
    Metaheuristics.create_child(x, fx)

Constructor for a solution depending on the result of `fx`.

## Example

```
julia> import Metaheuristics

julia> Metaheuristics.create_child(rand(3), 1.0)
| f(x) = 1
| solution.x = [0.2700437125780806, 0.5233263210622989, 0.12871108215859772]

julia> Metaheuristics.create_child(rand(3), (1.0, [2.0, 0.2], [3.0, 0.3]))
| f(x) = 1
| g(x) = [2.0, 0.2]
| h(x) = [3.0, 0.3]
| x = [0.9881102595664819, 0.4816273348099591, 0.7742585077942159]

julia> Metaheuristics.create_child(rand(3), ([-1, -2.0], [2.0, 0.2], [3.0, 0.3]))
| f(x) = [-1.0, -2.0]
| g(x) = [2.0, 0.2]
| h(x) = [3.0, 0.3]
| x = [0.23983577719146854, 0.3611544510766811, 0.7998754930109109]

julia> population = [ Metaheuristics.create_child(rand(2), (randn(2),  randn(2), rand(2))) for i = 1:100  ]
                           F space
          ┌────────────────────────────────────────┐
        2 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⠀⠂⠀⠀⠀⠀⠀⡇⠈⡀⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠘⠀⡇⠀⠀⠘⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠂⠀⠀⢀⠠⠐⠀⡇⠄⠁⠀⠀⠀⡀⠀⢁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠂⢈⠀⠈⡇⠀⡐⠃⠀⠄⠄⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠄⢐⠠⠀⡄⠀⠀⡇⠀⠂⠈⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠉⠉⠉⠉⠋⠉⠉⠉⠉⠉⠉⠙⢉⠉⠙⠉⠉⡏⠉⠉⠩⠋⠉⠩⠉⠉⠉⡉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉│
   f_2    │⠀⠀⠀⠀⠀⡀⠀⠀⠀⠄⠀⠀⡀⠀⠀⠂⠀⡇⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⠀⠐⡇⠠⠀⠀⠀⠈⢀⠄⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠄⠀⡀⠀⠂⡇⠐⠘⠈⠂⠀⠈⡀⠀⠀⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠄⠀⠀⠀⠀⠀⠂⠀⠂⠀⠀⡇⠀⠈⢀⠐⠀⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⠀⠀⢁⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       -3 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          └────────────────────────────────────────┘
          -3                                       4
                             f_1
```

"""
create_child(x, fx) = generateChild(x, fx)

function create_solution(x, problem::Problem; ε=0.0)
    problem.f_calls += 1
    return generateChild(x, problem.f(x), ε=ε)
end

# getters for the above structures

"""
    get_position(solution)

Get the position vector.
"""
get_position(solution::Solution) = solution.x

function positions(population::AbstractArray)
    if isempty(population)
        zeros(0,0)
    end

    D = length(get_position(population[1]))
    [ get_position(sol)[i] for sol in population, i in 1:D]
end
"""
    fval(solution)

Get the objective function value (fitness) of a solution.
"""
fval(solution::Solution) = solution.f
sum_violations(solution::T) where T <: AbstractConstrainedSolution = solution.sum_violations

"""
    gval(solution)

Get the inequality constraints of a solution.
"""
gval(solution::T) where T <: AbstractConstrainedSolution = solution.g

"""
    hval(solution)

Get the equality constraints of a solution.
"""
hval(solution::T) where T <: AbstractConstrainedSolution = solution.h

function fvals(population::AbstractArray{T}) where T <: AbstractMultiObjectiveSolution
    if isempty(population)
        return zeros(0, 0)
    end

    M = length(fval(population[1]))
    [ fval(sol)[i] for sol in population, i in 1:M]
end

fvals(population::AbstractArray{T}) where T <: AbstractSolution = fval.(population)

function gvals(population::AbstractArray{T}) where T <: AbstractConstrainedSolution
    if isempty(population)
        return zeros(0, 0)
    end

    M = length(gval(population[1]))
    [ gval(sol)[i] for sol in population, i in 1:M]
end

function hvals(population::AbstractArray{T}) where T <: AbstractConstrainedSolution
    if isempty(population)
        return zeros(0, 0)
    end

    M = length(hval(population[1]))
    [ hval(sol)[i] for sol in population, i in 1:M]
end

"""
    is_feasible(solution)
Returns `true` if solution is feasible, otherwise returns `false`.
"""
is_feasible(solution::T) where T <: AbstractConstrainedSolution = solution.is_feasible
is_feasible(solution::T) where T <: AbstractSolution = true

"""
    pareto_front(state::State)
Returns the non-dominated solutions in state.population.
"""
pareto_front(st::State) = pareto_front(st.population)

"""
    pareto_front(population::Array)
Returns non-dominated solutions.
"""
function pareto_front(population::AbstractArray{T}) where T <: AbstractMultiObjectiveSolution
    return fvals(get_non_dominated_solutions(population))
end


text:
####################################
### Types used in OpenStreetMapXSim package ###
####################################

"""
The `SimData` type is the type used to store all the crucial data  used in the OSM simulation module:

**Fields**
* `map_data` : data that have been processed from OpenStreetMap osm file
* `features` : dictionary of all features (shops, schools, etc.)
* `feature_classes` : dictionary mapping all features with the proper classes
* `feature_to_intersections` : dictionary mapping all features with the nearest intersection
* `DAs_to_intersection` : dictionary mapping all DA's with the nearest intersection
* `demographic_data` : dictionary containing a demographic data about each DA's
* `business_data` : array of dictionaries with a description of each business in the area
* `DAs_flow_dictionary` : dictionary mapping all DA's ID's with columns and rows of journey matrix
* `DAs_flow_matrix` : Journey matrix
* `googleapi_key` : non obligatory key for using a *Google Distances API*
"""
mutable struct SimData
	map_data::OpenStreetMapX.MapData
    features::Dict{Int,Tuple{String,String}}
    feature_classes::Dict{String,Int}
    feature_to_intersections::Dict{Int,Int}
    DAs_to_intersection::Dict{Int,Int}
	demographic_data::Dict{Int,Union{Dict{Symbol,Int}, Nothing}}
	business_data::Union{Array{Dict{Symbol,Union{String, Int,UnitRange{Int}}},1}, Nothing}
	DAs_flow_dictionary::Dict{Int,Int}
    DAs_flow_matrix::SparseArrays.SparseMatrixCSC{Int,Int}
	googleapi_key::Union{Nothing,String}
end

"""
The `Road` type is the type used to store the informations about the roads chosen by agents during the simulation

**Fields**
* `start_node` :  starting point of a road
* `fin_node` :  finish point of a road
* `waypoint` :  type of waypoint (school, recreation, etc.) chosen by agent during his trip (if exist)
* `mode` : field describing the way how agent selects the road (he can choose either fastest, shortest or based on Google Distances API road)
* `route` : list of intersections crossed during agent's trip
* `count` : how many agents choose this particular way
"""
mutable struct Road
    start_node::Int
    fin_node::Int
    waypoint::Union{String,Nothing}
    mode::String
    route::Array{Int,1}
    count::Int
end

"""
The `NodeStat` type is the type used to store the informations about each intersection (graph vertex) in simulation and  demographic profiles of agents crossing this particular intersection.

**Fields**
* `count` : indicates how many agents had crossed this intersection
* `latitude` : latitude of the vertex
* `longitude` : longitude of the vertex
* `agents_data` : a data frame with demographic profiles of agents
"""
mutable struct NodeStat
    count::Int
    latitude::Float64
    longitude::Float64
    agents_data::Union{DataFrames.DataFrame, Nothing}
end


text:
# generic fallback
function formatinfo(s::Spec, x::Any)
    str = string(x)
    size, meta = formatinfo(s, str)
    return size, (str, meta)
end

function formatfield(data::Vector{UInt8}, p::Int, s::Spec, ::Any, (str, meta))
    return formatfield(data, p, s, str, meta)
end


text:
# include plot3d of Winston
filePath = joinpath(Pkg.dir("Winston"), "src", "canvas3d.jl")
include(filePath)

function wstPlot(cop::PairCop, n=50::Int64;
                 f=Copulas.pdf::Function,
                 zMax=Inf::Float64)

    cutoff = zMax
    grid = linspace(0, 1, n)
    grid[1] = 0.001
    grid[end] = 0.999

    Plot3d.plot3d(
                  Plot3d.surf((u,v)->u,
                              (u,v)->min(f(cop, u, v)[1], cutoff),
                              (u,v)->v,
                              grid, grid))
end


text:
"""
    mutable struct SingleTimeSeries <: StaticTimeSeries
        name::String
        data::TimeSeries.TimeArray
        scaling_factor_multiplier::Union{Nothing, Function}
        internal::InfrastructureSystemsInternal
    end

A deterministic forecast for a particular data field in a Component.

# Arguments
- `name::String`: user-defined name
- `data::TimeSeries.TimeArray`: timestamp - scalingfactor
- `scaling_factor_multiplier::Union{Nothing, Function}`: Applicable when the time series
  data are scaling factors. Called on the associated component to convert the values.
- `internal::InfrastructureSystemsInternal`
"""
mutable struct SingleTimeSeries <: StaticTimeSeries
    "user-defined name"
    name::String
    "timestamp - scalingfactor"
    data::TimeSeries.TimeArray
    resolution::Dates.Period
    "Applicable when the time series data are scaling factors. Called on the associated component to convert the values."
    scaling_factor_multiplier::Union{Nothing, Function}
    internal::InfrastructureSystemsInternal
end

function SingleTimeSeries(;
    name,
    data,
    scaling_factor_multiplier = nothing,
    normalization_factor = 1.0,
    internal = InfrastructureSystemsInternal(),
)
    data = handle_normalization_factor(data, normalization_factor)
    return SingleTimeSeries(
        name,
        data,
        _get_resolution(data),
        scaling_factor_multiplier,
        internal,
    )
end

function _get_resolution(data::TimeSeries.TimeArray)
    if length(data) < 2
        throw(
            ConflictingInputsError(
                "Resolution can't be inferred from the data. Please select an appropiate constructor.",
            ),
        )
    end
    return TimeSeries.timestamp(data)[2] - TimeSeries.timestamp(data)[1]
end

"""
Construct SingleTimeSeries from a TimeArray or DataFrame.

# Arguments
- `name::AbstractString`: user-defined name
- `data::Union{TimeSeries.TimeArray, DataFrames.DataFrame}`: time series data
- `normalization_factor::NormalizationFactor = 1.0`: optional normalization factor to apply
  to each data entry
- `scaling_factor_multiplier::Union{Nothing, Function} = nothing`: If the data are scaling
  factors then this function will be called on the component and applied to the data when
  [`get_time_series_array`](@ref) is called.
- `timestamp = :timestamp`: If a DataFrame is passed then this must be the column name that
  contains timestamps.
"""
function SingleTimeSeries(
    name::AbstractString,
    data::Union{TimeSeries.TimeArray, DataFrames.DataFrame};
    normalization_factor::NormalizationFactor = 1.0,
    scaling_factor_multiplier::Union{Nothing, Function} = nothing,
    timestamp = :timestamp,
)
    if data isa DataFrames.DataFrame
        ta = TimeSeries.TimeArray(data; timestamp = timestamp)
    elseif data isa TimeSeries.TimeArray
        ta = data
    else
        error("fatal: $(typeof(data))")
    end

    return SingleTimeSeries(
        name = name,
        data = ta,
        scaling_factor_multiplier = scaling_factor_multiplier,
        normalization_factor = normalization_factor,
        internal = InfrastructureSystemsInternal(),
    )
end

"""
Construct SingleTimeSeries from a CSV file. The file must have a column that is the name of the
component.

# Arguments
- `name::AbstractString`: user-defined name
- `filename::AbstractString`: name of CSV file containing data
- `component::InfrastructureSystemsComponent`: component associated with the data
- `resolution::Dates.Period`: resolution of the time series
- `normalization_factor::NormalizationFactor = 1.0`: optional normalization factor to apply
  to each data entry
- `scaling_factor_multiplier::Union{Nothing, Function} = nothing`: If the data are scaling
  factors then this function will be called on the component and applied to the data when
  [`get_time_series_array`](@ref) is called.
"""
function SingleTimeSeries(
    name::AbstractString,
    filename::AbstractString,
    component::InfrastructureSystemsComponent,
    resolution::Dates.Period;
    normalization_factor::NormalizationFactor = 1.0,
    scaling_factor_multiplier::Union{Nothing, Function} = nothing,
)
    component_name = get_name(component)
    raw = read_time_series(SingleTimeSeries, filename, component_name)
    ta = make_time_array(raw, component_name, resolution)
    return SingleTimeSeries(
        name = name,
        data = ta,
        normalization_factor = normalization_factor,
        scaling_factor_multiplier = scaling_factor_multiplier,
    )
end

"""
Construct SingleTimeSeries after constructing a TimeArray from `initial_time` and
`time_steps`.
"""
function SingleTimeSeries(
    name::String,
    resolution::Dates.Period,
    initial_time::Dates.DateTime,
    time_steps::Int,
)
    data = TimeSeries.TimeArray(
        initial_time:resolution:(initial_time + resolution * (time_steps - 1)),
        ones(time_steps),
    )
    return SingleTimeSeries(; name = name, data = data)
end

function SingleTimeSeries(time_series::Vector{SingleTimeSeries})
    @assert !isempty(time_series)
    timestamps =
        collect(Iterators.flatten((TimeSeries.timestamp(get_data(x)) for x in time_series)))
    data = collect(Iterators.flatten((TimeSeries.values(get_data(x)) for x in time_series)))
    ta = TimeSeries.TimeArray(timestamps, data)

    time_series = SingleTimeSeries(
        name = get_name(time_series[1]),
        data = ta,
        scaling_factor_multiplier = time_series[1].scaling_factor_multiplier,
    )
    @debug "concatenated time_series" LOG_GROUP_TIME_SERIES time_series
    return time_series
end

function SingleTimeSeries(ts_metadata::SingleTimeSeriesMetadata, data::TimeSeries.TimeArray)
    return SingleTimeSeries(
        get_name(ts_metadata),
        data,
        get_resolution(ts_metadata),
        get_scaling_factor_multiplier(ts_metadata),
        InfrastructureSystemsInternal(get_time_series_uuid(ts_metadata)),
    )
end

function SingleTimeSeriesMetadata(ts::SingleTimeSeries)
    return SingleTimeSeriesMetadata(
        get_name(ts),
        get_resolution(ts),
        get_initial_timestamp(ts),
        get_uuid(ts),
        length(ts),
        get_scaling_factor_multiplier(ts),
    )
end

function SingleTimeSeries(info::TimeSeriesParsedInfo)
    data = make_time_array(info)
    return SingleTimeSeries(
        name = info.name,
        data = data,
        normalization_factor = info.normalization_factor,
        scaling_factor_multiplier = info.scaling_factor_multiplier,
    )
end

"""Get [`SingleTimeSeries`](@ref) `name`."""
get_name(value::SingleTimeSeries) = value.name
"""Get [`SingleTimeSeries`](@ref) `data`."""
get_data(value::SingleTimeSeries) = value.data
"""Get [`SingleTimeSeries`](@ref) `resolution`."""
get_resolution(value::SingleTimeSeries) = value.resolution
"""Get [`SingleTimeSeries`](@ref) `scaling_factor_multiplier`."""
get_scaling_factor_multiplier(value::SingleTimeSeries) = value.scaling_factor_multiplier
"""Get [`SingleTimeSeries`](@ref) `internal`."""
get_internal(value::SingleTimeSeries) = value.internal

"""Set [`SingleTimeSeries`](@ref) `name`."""
set_name!(value::SingleTimeSeries, val) = value.name = val
"""Set [`SingleTimeSeries`](@ref) `data`."""
set_data!(value::SingleTimeSeries, val) = value.data = val
"""Set [`SingleTimeSeries`](@ref) `scaling_factor_multiplier`."""
set_scaling_factor_multiplier!(value::SingleTimeSeries, val) =
    value.scaling_factor_multiplier = val
"""Set [`SingleTimeSeries`](@ref) `internal`."""
set_internal!(value::SingleTimeSeries, val) = value.internal = val

eltype_data(ts::SingleTimeSeries) = eltype(TimeSeries.values(ts.data))

get_initial_timestamp(time_series::SingleTimeSeries) =
    TimeSeries.timestamp(get_data(time_series))[1]

function get_array_for_hdf(ts::SingleTimeSeries)
    return transform_array_for_hdf(TimeSeries.values(ts.data))
end

function Base.getindex(time_series::SingleTimeSeries, args...)
    return SingleTimeSeries(time_series, getindex(get_data(time_series), args...))
end

Base.first(time_series::SingleTimeSeries) = head(time_series, 1)

Base.last(time_series::SingleTimeSeries) = tail(time_series, 1)

Base.firstindex(time_series::SingleTimeSeries) = firstindex(get_data(time_series))

Base.lastindex(time_series::SingleTimeSeries) = lastindex(get_data(time_series))

Base.lastindex(time_series::SingleTimeSeries, d) = lastindex(get_data(time_series), d)

Base.eachindex(time_series::SingleTimeSeries) = eachindex(get_data(time_series))

Base.iterate(time_series::SingleTimeSeries, n = 1) = iterate(get_data(time_series), n)

"""
Refer to TimeSeries.when(). Underlying data is copied.
"""
function when(time_series::SingleTimeSeries, period::Function, t::Integer)
    return SingleTimeSeries(time_series, TimeSeries.when(get_data(time_series), period, t))
end

"""
Return a time_series truncated starting with timestamp.
"""
function from(time_series::SingleTimeSeries, timestamp)
    return SingleTimeSeries(
        name = get_name(time_series),
        data = TimeSeries.from(get_data(time_series), timestamp),
    )
end

"""
Return a time_series truncated after timestamp.
"""
function to(time_series::SingleTimeSeries, timestamp)
    return SingleTimeSeries(
        name = get_name(time_series),
        data = TimeSeries.to(get_data(time_series), timestamp),
    )
end

"""
Return a time_series with only the first num values.
"""
function head(time_series::SingleTimeSeries)
    return SingleTimeSeries(time_series, TimeSeries.head(get_data(time_series)))
end

function head(time_series::SingleTimeSeries, num)
    return SingleTimeSeries(time_series, TimeSeries.head(get_data(time_series), num))
end

"""
Return a time_series with only the ending num values.
"""
function tail(time_series::SingleTimeSeries)
    return SingleTimeSeries(time_series, TimeSeries.tail(get_data(time_series)))
end

function tail(time_series::SingleTimeSeries, num)
    return SingleTimeSeries(time_series, TimeSeries.tail(get_data(time_series), num))
end

"""
Creates a new SingleTimeSeries from an existing instance and a subset of data.
"""
function SingleTimeSeries(time_series::SingleTimeSeries, data::TimeSeries.TimeArray)
    vals = []
    for (fname, ftype) in zip(fieldnames(SingleTimeSeries), fieldtypes(SingleTimeSeries))
        if ftype <: TimeSeries.TimeArray
            val = data
        elseif ftype <: InfrastructureSystemsInternal
            # Need to create a new UUID.
            val = InfrastructureSystemsInternal()
        else
            val = getfield(time_series, fname)
        end

        push!(vals, val)
    end

    return SingleTimeSeries(vals...)
end

get_columns(::Type{<:TimeSeriesMetadata}, ta::TimeSeries.TimeArray) = nothing

function make_time_array(
    time_series::SingleTimeSeries,
    start_time::Dates.DateTime;
    len::Union{Nothing, Int} = nothing,
)
    ta = get_data(time_series)
    first_time = first(TimeSeries.timestamp(ta))
    if start_time == first_time && (len === nothing || len == length(ta))
        return ta
    end

    resolution = Dates.Millisecond(get_resolution(ta))
    start_index = Int((start_time - first_time) / resolution) + 1
    end_index = start_index + len - 1
    return ta[start_index:end_index]
end

function SingleTimeSeriesMetadata(ts_metadata::DeterministicMetadata)
    return SingleTimeSeriesMetadata(
        name = get_name(ts_metadata),
        resolution = get_resolution(ts_metadata),
        initial_timestamp = get_initial_timestamp(ts_metadata),
        time_series_uuid = get_time_series_uuid(ts_metadata),
        length = get_count(ts_metadata) * get_horizon(ts_metadata),
        scaling_factor_multiplier = get_scaling_factor_multiplier(ts_metadata),
        internal = get_internal(ts_metadata),
    )
end


text:
module PumpProbePlotting

using PyPlot
using PumpProbeModels
using PumpProbeFitting

export plot

"""
Plots the fit result
"""
function PyPlot.plot(r::PumpProbeFitting.PumpProbeFit, m::PumpProbeModel, dls, wns)
    x = r.m.wavenumbers
    y = r.m.delaytimes
    zfit = model2data(r.m)
    zdata = fit2data(r)
    resid = reshape(r.f.resid, size(zdata))

    pump = r.m.pumpfunction.(r.m.pumptimes)
    t0 = r.m.parameters.t0

    plotmodel = false
    local zmodel
    if m !== r.m
        zmodel = model2data(m)
        plotmodel = true
    end

    # Get extrema of the data
    vmin = min(zdata..., zfit...)
    vmax = max(zdata..., zfit...)
    rmax = maximum(abs.(resid))     # Limits for residual plot

    # Selected delay times and wavenumbers of slices
    dlt = y[dls]
    wn  = x[wns]

    fh = figure(figsize=(12,12))

    subplot(331)
    pcolormesh(x,y,zdata, vmin=vmin, vmax=vmax)
    colorbar()
    title("Input Data")

    subplot(334)
    pcolormesh(x,y,zfit, vmin=vmin, vmax=vmax)
    colorbar()
    title("Fit")
    hlines(dlt, xlim()..., alpha=0.4)
    vlines(wn , ylim()..., alpha=0.4)

    subplot(337)
    title("Residuals")
    pcolormesh(x,y,resid, cmap=ColorMap("RdBu"), vmin=-rmax, vmax=rmax)
    colorbar()

    ND = length(dls)
    for (i,k) in dls |> enumerate
        subplot(ND,3,3i-1)
        plot(x, zdata[k,:], ".", color="k", alpha=0.3)
        plot(x, zfit[k,:], "-", color="k", label="$(round(dlt[i], digits=1)) ps")
        plot(x, zmodel[k,:], "-", color="k", alpha=0.3)
        legend()
        # Also plot
        # Remove Ticks
        i != length(dls) && gca().set_xticklabels([])
    end

    NW = length(wns)
    for (i,k) in wns |> enumerate
        subplot(NW+1,3,3i)
        plot(y, zdata[:,k], ".", color="k", alpha=0.3)
        plot(y, zfit[:,k], "-", color="k", label="$(round(wn[i], digits=1)) 1/cm")
        plot(y, zmodel[:,k], "-", color="k", alpha=0.3)
        vlines(t0, ylim()..., linestyle="--", alpha=0.3)
        legend()
        # Remove Ticks
        i != length(wns) && gca().set_xticklabels([])
    end

    # Also plot pumpfunction
    xlims = xlim()
    subplot(NW+1,3,3*(NW+1))
    plot(r.m.pumptimes .+ t0, pump, ".-", color="k", alpha=0.2, label="Pump Pulse")
    xlim(xlims...)
    legend()

    plt.tight_layout()
    fh
end

end


text:
module Solvent

using Adapt
using CUDA
using LinearAlgebra
using LazyArrays
using StaticArrays
using KernelAbstractions

array_device(::Union{Array, SArray, MArray}) = CPU()
array_device(::CuArray) = CUDADevice()
realview(x::Union{Array, SArray, MArray}) = x
realview(x::CuArray) = x

# Just for testing Solvent
LinearAlgebra.norm(A::MVector, p::Real, weighted::Bool) = norm(A, p)
LinearAlgebra.norm(A::MVector, weighted::Bool) = norm(A, 2, weighted)
LinearAlgebra.dot(A::MVector, B::MVector, weighted) = dot(A, B)
LinearAlgebra.norm(A::AbstractVector, p::Real, weighted::Bool) = norm(A, p)
LinearAlgebra.norm(A::AbstractVector, weighted::Bool) = norm(A, 2, weighted)
LinearAlgebra.dot(A::AbstractVector, B::AbstractVector, weighted) = dot(A, B)

export linearsolve!, settolerance!, prefactorize
export AbstractSystemSolver, AbstractIterativeSystemSolver

"""
    AbstractSystemSolver

This is an abstract type representing a generic linear solver.
"""
abstract type AbstractSystemSolver end

include(joinpath("Preconditioners", "Preconditioners.jl"))
include(joinpath("LinearSolvers", "LinearSolvers.jl"))

end  # End of module

max_stars_repo_name: CliMA/SystemSolvers

max_issues_repo_name: CliMA/SystemSolvers


text:
"""

`leftjoin(base, src, join_axis...)`

Left join a `LabeledArray` into another `LabeledArray`.
Note that the left array (base) can be multidimensional. The function creates a dictionary from the right array (`src`).

##### Arguments

* `base` : the left `LabeledArray`.
* `src` : the right `LabeledArray`.
* `join_axes...` can be one of the following forms:
    * integers for the directions in `src` along which to join. In this case, the keys is the `base` side are found by matching the field names in the join directions in `src` with those in `base`.
    * a list of integer=>integer or integer=>vector of arrays, each of the same shape as `base`.

Ultimately, `join_axes...` produces pairs of direction in `src` => vector of arrays, each of the shape of `base`. If the value in key=>value is an integer, the axis along that direction in `base` is taken, after broadcast. The field values are combined into a vector of arrays. If the right hand side is missing (i.e. just an integer), the field names in the axis along the integer direction are used to create an array for `base`.

##### Return

A left joined `LabeledArray`. The join is performed as follows: Given an `i=>arr` form as an element in `join_axes`, the keys in `i`th direction in `src` are used as keys and `arr` are used the keys in the `base` side to left join. The values will be the sliced subarrays for each value in the `join_axes`. Note that `join_axis...` chooses multiple axes for keys.
The output number of dimensions is `ndims(base) + ndims(src) - length(join_axes)`.
Note that when `join_axis` is empty, the result is the tensor product of `base` and `src` from `tensorprod`.

##### Examples

```julia
julia> b = larr(k=[:x :x :y;:z :u :v], axis1=[:x,:y], axis2=darr(r=[:x, :y, :z]))
2 x 3 LabeledArray

r |x |y |z
--+--+--+--
  |k |k |k
--+--+--+--
x |x |x |y
y |z |u |v

julia> s = larr(axis1=darr(k=[:x,:y,:z,:m,:n,:p]), b=[1,2,3,4,5,6])
6 LabeledArray

k |b
--+--
x |1
y |2
z |3
m |4
n |5
p |6

julia> leftjoin(b, s, 1)
2 x 3 LabeledArray

r |x   |y   |z
--+----+----+----
  |k b |k b |k b
--+----+----+----
x |x 1 |x 1 |y 2
y |z 3 |u   |v

julia> leftjoin(b, s, 1=>1)
2 x 3 LabeledArray

r |x   |y   |z
--+----+----+----
  |k b |k b |k b
--+----+----+----
x |x 1 |x 1 |y 1
y |z 2 |u 2 |v 2

julia> leftjoin(b, s, 1=>Any[nalift([:x :z :n;:y :m :p])])
2 x 3 LabeledArray

r |x   |y   |z
--+----+----+----
  |k b |k b |k b
--+----+----+----
x |x 1 |x 3 |y 5
y |z 2 |u 4 |v 6
```

"""
function leftjoin end

leftjoin{N}(base::LabeledArray{TypeVar(:V),N}, src::LabeledArray, join_axes::Int...) = begin
  # needs to find corresponding fields in base, so src axes for joining should be DictArray.
  for axis_index in join_axes
    @assert isa(src.axes[axis_index], DictArray)
  end
  base_key_arrays = AbstractArray{TypeVar(:V),N}[]
  for axis_index in join_axes
    for field in src.axes[axis_index].data.keys
      push!(base_key_arrays, selectfield(base, field))
    end
  end
  perform_leftjoin(base, src, (base_key_arrays...), join_axes)
end

leftjoin(base::LabeledArray, src::LabeledArray) = tensorprod(base, src)

leftjoin(base::DictArray, src::LabeledArray, join_axes::Int...) = leftjoin(LabeledArray(base), src, join_axes...)
leftjoin(base::DictArray, src::LabeledArray) = leftjoin(LabeledArray(base), src)
leftjoin(base::DictArray, src::LabeledArray, join_pairs::Pair{Int}...) = leftjoin(LabeledArray(base), src, join_pairs...)

leftjoin{N}(base::LabeledArray{TypeVar(:V),N}, src::LabeledArray, join_pairs::Pair{Int}...) = begin
  join_axes = Int[join_pair[1] for join_pair in join_pairs]
  axis_labeled_src = LabeledArray(src.data, ntuple(length(src.axes)) do d
    srcaxesd = src.axes[d]
    if in(d, join_axes) && !isa(srcaxesd, DictArray)
      create_dictarray_nocheck(create_ldict_nocheck(:dummy => srcaxesd))
    else
      srcaxesd
    end
  end)
  base_key_arrays = AbstractArray{TypeVar(:V),N}[]
  for (axis_index, axis_values0) in join_pairs
    axis_values = lift_axis_values0_helper1(axis_values0, base, src.axes[axis_index])
    push!(base_key_arrays, axis_values...)
  end
  perform_leftjoin(base, axis_labeled_src, (base_key_arrays...), (join_axes...))
end

lift_axis_values0_helper1(axis_values0::Integer, base, srcaxis) = lift_axis_values0_helper2(BroadcastAxis(base.axes[axis_values0], base, axis_values0), srcaxis)
lift_axis_values0_helper1(axis_values0::DictArray, base, srcaxis) = values(axis_values0)
lift_axis_values0_helper1(axis_values0, base, srcaxis) = axis_values0

lift_axis_values0_helper2(arr::DictArray, srcaxis::DictArray) = begin
  if keys(srcaxis) != keys(arr)
    throw(ArgumentError("field names do not match: $(keys(srcaxis)) vs $(keys(arr))"))
  end
  values(arr)
end
lift_axis_values0_helper2(arr::DictArray, srcaxis) = values(arr)
lift_axis_values0_helper2(arr, srcaxis) = Any[arr]

@generated perform_leftjoin{N,M,L,K}(base::LabeledArray{TypeVar(:T),N},
               src::LabeledArray{TypeVar(:U),M},
               base_key_arrays::NTuple{L,AbstractArray{TypeVar(:V),N}},
               src_key_axis_indices::NTuple{K,Int},
               concat_function::Function=default_concat_array_function) = begin
  rest_dims = M - K
  quote
    @assert $L == sum(i->length(src.axes[i].data), src_key_axis_indices)
    rest_src_axis_indices = filter(i->!(i in src_key_axis_indices), 1:$M) #length(src.axes))
    new_src_indices_ordering = [rest_src_axis_indices...;src_key_axis_indices...]
    permuted_src = permutedims_if_necessary(src, (new_src_indices_ordering...))
    src_axis_index_offset = length(rest_src_axis_indices)
    bymaps = ntuple($K) do i
      axis = permuted_src.axes[i+src_axis_index_offset] #src_key_axis_indices[i]]
      create_join_bymap(axis, Tuple{[eltype(v) for v in axis.data.values]...})
    end
    offset = 1
    bymaps_tuple_positions = ntuple($K) do i
      r = length(src.axes[i+src_axis_index_offset].data)
      offset0 = offset
      offset += r
      offset0
    end
    base_to_src_loc_map = similar(base_key_arrays[1], NTuple{K,Int})
    create_base_to_src_loc_map_all!(base_to_src_loc_map, base_key_arrays, ntuple(d->0,length(base_key_arrays)), bymaps_tuple_positions, bymaps)
    extra_dims = ntuple(d->length(permuted_src.axes[d]), $rest_dims)
    extra_axes = ntuple(d->permuted_src.axes[d], $rest_dims)
    src_mapped_to_base = LabeledArray(similar(src.data, (size(base_to_src_loc_map)...,extra_dims...)), (base.axes...,extra_axes...))
    srctobasedata = src_mapped_to_base.data
    permutedsrcdata = permuted_src.data

    fill_srctobasedata_join_helper!(srctobasedata, base_to_src_loc_map, permutedsrcdata)

    lifted_base = expand_dims(base, (), extra_axes)
    concat_function(lifted_base, src_mapped_to_base)
  end
end

"""

`innerjoin(base, src, join_axis...)`

Inner join an LabeledArray into another LabeledArray. `innerjoin` is different from `leftjoin` in that only elements in the left array that have the corresponding elements in the right array will be kept. Otherwise, the elements will be set to null. If the entire elements along some direction are null, they will be all removed in the output.
Note that the left array (base) can be multidimensional. The function creates a dictionary from the right array (`src`).

##### Arguments

* `base` : the left `LabeledArray`.
* `src` : the right `LabeledArray`.
* `join_axes...` can be one of the following forms:
    * integers for the directions in `src` along which to join.
    * a list of integer=>integer or integer=>vector of arrays, each array of the same shape as `base`.

Ultimately, `join_axes...` produces pairs of direction in `src` => vector of arrays, each of the shape of `base`. If the value in key=>value is an integer, the axis along that direction in `base` is taken, after broadcast. The field values are combined into a vector of arrays. If the right hand side is missing (i.e. just an integer), the field names in the axis along the integer direction are used to create an array for `base`.

##### Return

An inner joined `LabeledArray`. The join is performed as follows: Given an `i=>arr` form as an element in `join_axes`, the keys in `i`th direction in `src` are used as keys and `arr` are used the keys in the `base` side to inner join. The values will be the sliced subarrays for each value in the `join_axes`. Note that `join_axis...` chooses multiple axes for keys.
The output number of dimensions is `ndims(base) + ndims(src) - length(join_axes)`.
Note that when `join_axis` is empty, the result is the tensor product of `base` and `src` from `tensorprod`.

##### Examples

```julia
julia> b = larr(k=[:x :x :y;:z :u :v], axis1=[:x,:u], axis2=darr(r=[:x, :y, :z]))
2 x 3 LabeledArray

r |x |y |z
--+--+--+--
  |k |k |k
--+--+--+--
x |x |x |y
u |z |u |v

julia> s = larr(axis1=darr(k=[:x,:y,:z,:m,:n,:p]), b=[1,2,3,4,5,6])
6 LabeledArray

k |b
--+--
x |1
y |2
z |3
m |4
n |5
p |6

julia> innerjoin(b, s, 1)
2 x 3 LabeledArray

r |x   |y   |z
--+----+----+----
  |k b |k b |k b
--+----+----+----
x |x 1 |x 1 |y 2
u |z 3 |u   |v

julia> innerjoin(b, s, 1=>1)
1 x 3 LabeledArray

r |x   |y   |z
--+----+----+----
  |k b |k b |k b
--+----+----+----
x |x 1 |x 1 |y 1

julia> innerjoin(b, s, 1=>Any[nalift([:o :x :x;:q :r :y])])
2 x 2 LabeledArray

r |y   |z
--+----+----
  |k b |k b
--+----+----
x |x 1 |y 1
u |u   |v 2
```

"""
function innerjoin end

innerjoin{N}(base::LabeledArray{TypeVar(:V),N}, src::LabeledArray, join_axes::Int...) = begin
  # needs to find corresponding fields in base, so src axes for joining should be DictArray.
  for axis_index in join_axes
    @assert isa(src.axes[axis_index], DictArray)
  end
  base_key_arrays = AbstractArray{TypeVar(:V),N}[]
  for axis_index in join_axes
    for field in src.axes[axis_index].data.keys
      push!(base_key_arrays, selectfield(base, field))
    end
  end
  perform_innerjoin(base, src, (base_key_arrays...), join_axes)
end
innerjoin(base::LabeledArray, src::LabeledArray) = tensorprod(base, src)

innerjoin(base::DictArray, src::LabeledArray, join_axes::Int...) = innerjoin(LabeledArray(base), src, join_axes...)
innerjoin(base::DictArray, src::LabeledArray) = innerjoin(LabeledArray(base), src)
innerjoin(base::DictArray, src::LabeledArray, join_pairs::Pair{Int}...) = innerjoin(LabeledArray(base), src, join_pairs...)

innerjoin{N}(base::LabeledArray{TypeVar(:V),N}, src::LabeledArray, join_pairs::Pair{Int}...) = begin
  join_axes = Int[join_pair[1] for join_pair in join_pairs]
  axis_labeled_src = LabeledArray(src.data, ntuple(length(src.axes)) do d
    srcaxesd = src.axes[d]
    if in(d, join_axes) && !isa(srcaxesd, DictArray)
      create_dictarray_nocheck(create_ldict_nocheck(:dummy => srcaxesd))
    else
      srcaxesd
    end
  end)
  base_key_arrays = AbstractArray{TypeVar(:V),N}[]
  for (axis_index, axis_values0) in join_pairs
    axis_values = lift_axis_values0_helper1(axis_values0, base, src.axes[axis_index])
    push!(base_key_arrays, axis_values...)
  end
  perform_innerjoin(base, axis_labeled_src, (base_key_arrays...), (join_axes...))
end

@generated perform_innerjoin{N,M,L,K}(base::LabeledArray{TypeVar(:T),N},
               src::LabeledArray{TypeVar(:U),M},
               base_key_arrays::NTuple{L,AbstractArray{TypeVar(:V),N}},
               src_key_axis_indices::NTuple{K,Int},
               concat_function::Function=default_concat_array_function) = begin
  rest_dims = M - K
  quote
    @assert $L == sum(i->length(src.axes[i].data), src_key_axis_indices)
    rest_src_axis_indices = filter(i->!(i in src_key_axis_indices), 1:$M) #length(src.axes))
    new_src_indices_ordering = [rest_src_axis_indices...;src_key_axis_indices...]
    permuted_src = permutedims(src, (new_src_indices_ordering...))
    src_axis_index_offset = length(rest_src_axis_indices)
    bymaps = ntuple($K) do i
      axis = permuted_src.axes[i+src_axis_index_offset] #src_key_axis_indices[i]]
      create_join_bymap(axis, Tuple{[eltype(v) for v in axis.data.values]...})
    end
    offset = 1
    bymaps_tuple_positions = ntuple($K) do i
      r = length(src.axes[i+src_axis_index_offset].data)
      offset0 = offset
      offset += r
      offset0
    end
    base_to_src_loc_map_all = similar(base_key_arrays[1], NTuple{K,Int})
    create_base_to_src_loc_map_all!(base_to_src_loc_map_all, base_key_arrays, ntuple(d->0,length(base_key_arrays)), bymaps_tuple_positions, bymaps)
    base_axes_to_display = ntuple(d->falses(size(base,d)), $N)
    join_helper_populate_base_axes_to_display!(base_axes_to_display, base_to_src_loc_map_all)
    base_to_src_loc_map = base_to_src_loc_map_all[base_axes_to_display...]
    extra_dims = ntuple(d->length(permuted_src.axes[d]), $rest_dims)
    extra_axes = ntuple(d->permuted_src.axes[d], $rest_dims)
    src_mapped_to_base = LabeledArray(similar(src.data, (size(base_to_src_loc_map)...,extra_dims...)),
                                      ([axis[base_axes_to_display[i]] for (i,axis) in enumerate(base.axes)]...,extra_axes...))
    srctobasedata = src_mapped_to_base.data
    permutedsrcdata = permuted_src.data
    fill_srctobasedata_join_helper!(srctobasedata, base_to_src_loc_map, permutedsrcdata)

    lifted_base = expand_dims(base[base_axes_to_display...], (), extra_axes)
    concat_function(lifted_base, src_mapped_to_base)
  end
end

@generated join_helper_populate_base_axes_to_display!{T,K,N}(base_axes_to_display::T,
                                                             base_to_src_loc_map_all::AbstractArray{NTuple{K,Int},N}) = quote
  @nloops $N i base_to_src_loc_map_all begin
    if @nref($N,base_to_src_loc_map_all,i)[1] > 0
      @nexprs $N j->(base_axes_to_display[j][i_j] = true)
    end
  end
end

@generated fill_srctobasedata_join_helper!{N,K,MKN,M}(srctobasedata::AbstractArray{TypeVar(:KK),MKN},
                                                      base_to_src_loc_map::AbstractArray{NTuple{K,Int},N},
                                                      permutedsrcdata::AbstractArray{TypeVar(:TT),M}) = begin
  rest_dims = M - K
  @assert(MKN == M-K+N)
  quote
    @nloops $N i base_to_src_loc_map begin
      src_coords = @nref($N,base_to_src_loc_map,i)
      if src_coords[1] == 0
        setna!(srctobasedata, @ntuple($N,i)..., @ntuple($rest_dims, d->Colon())...)
      else
        setindex_nocheck!(srctobasedata, @nref($M,permutedsrcdata,d->d<=$rest_dims ? Colon() : src_coords[d-$rest_dims]), @ntuple($N,i)...,@ntuple($rest_dims,d->Colon())...)
      end
    end
  end
end

fill_srctobasedata_join_helper!{N,K,MKN,M}(srctobasedata::DictArray{TypeVar(:KK),MKN},
                                           base_to_src_loc_map::AbstractArray{NTuple{K,Int},N},
                                           permutedsrcdata::DictArray{TypeVar(:TT),M}) = begin
  for i in eachindex(srctobasedata.data.values, permutedsrcdata.data.values)
    fill_srctobasedata_join_helper!(srctobasedata.data.values[i], base_to_src_loc_map, permutedsrcdata.data.values[i])
  end
end

create_join_bymap{T}(axis::DictArray, ::Type{T}) = begin
  result = Dict{T,Int}()
  for i in 1:length(axis)
    result[getindexvalue(axis,i)] = i
  end
  result
end

@generated create_base_to_src_loc_map_all!{N,L,V1,A}(result::AbstractArray{NTuple{1,Int},N},
                                      base_key_arrays::A, #NTuple{L}, #::NTuple{L,AbstractArray{TypeVar(:V),N}},
                                      ::NTuple{L,Int},
                                      bymaps_tuple_positions::NTuple{1,Int},
                                      bymaps::Tuple{Dict{V1,Int}}) = quote
  zerocoords = (0,) # K=1
  bymaps1::Dict{V1,Int} = bymaps[1]
  for i in eachindex(result, base_key_arrays...)
    k::V1 = @ntuple $L d->base_key_arrays[d][i]
    result[i] = if haskey(bymaps1, k)
      (bymaps1[k],)
    else
      zerocoords
    end
  end
end

@generated create_base_to_src_loc_map_all!{K,N,L,A}(result::AbstractArray{NTuple{K,Int},N},
                                    base_key_arrays::A, #::NTuple{L}, #::NTuple{L,AbstractArray{TypeVar(:V),N}},
                                    ::NTuple{L,Int},
                                    bymaps_tuple_positions::NTuple{K,Int},
                                    bymaps) = quote
  zerocoords = ntuple(k->0, $K)
  for i in eachindex(result, base_key_arrays...)
    base_indices = @ntuple $L d->base_key_arrays[d][i]
    na_flag = false
    src_coords = ntuple($K) do i
      k = base_indices[bymaps_tuple_positions[i]:(i==$K ? L:bymaps_tuple_positions[i+1]-1)]
      bymapsi = bymaps[i]
      if haskey(bymapsi, k)
        bymapsi[k]
      else
        na_flag = true
        0
      end
    end
    result[i] = if na_flag
      zerocoords
    else
      src_coords
    end
  end
end

default_concat_array_function{N}(x::LabeledArray{TypeVar(:T),N}, y::LabeledArray{TypeVar(:T),N}) = begin
  tracker = Any[allfieldnames(x)...;allfieldnames(y)...] # need Array, not AbstractArrayWrapper, hence ...
  xdata = if isa(x.data, DictArray)
    x.data
  else
    create_dictarray_nocheck(create_ldict_nocheck(create_additional_fieldname(x, tracker) => x.data))
  end
  ydata = if isa(y.data, DictArray)
    y.data
  else
    create_dictarray_nocheck(create_ldict_nocheck(create_additional_fieldname(y, tracker) => y.data))
  end
  data = DictArray(LDict([xdata.data.keys;ydata.data.keys], [xdata.data.values;ydata.data.values]))
  axes = x.axes
  LabeledArray(data, axes)
end


text:
# Routines for interaction with spglib
# Note: spglib/C uses the row-major convention, thus we need to perform transposes
#       between julia and spglib (https://spglib.github.io/spglib/variable.html)
import Spglib
const SPGLIB = spglib_jll.libsymspg

function spglib_get_error_message()
    error_code = ccall((:spg_get_error_code, SPGLIB), Cint, ())
    return unsafe_string(ccall((:spg_get_error_message, SPGLIB), Cstring, (Cint,), error_code))
end

"""
Convert the DFTK atoms datastructure into a tuple of datastructures for use with spglib.
`positions` contains positions per atom, `numbers` contains the mapping atom
to a unique number for each indistinguishable element, `spins` contains
the ``z``-component of the initial magnetic moment on each atom, `mapping` contains the
mapping of the `numbers` to the element objects in DFTK and `collinear` whether
the atoms mark a case of collinear spin or not. Notice that if `collinear` is false
then `spins` is garbage.
"""
function spglib_atoms(atoms, magnetic_moments=[])
    n_attypes = isempty(atoms) ? 0 : sum(length(positions) for (typ, positions) in atoms)
    spg_numbers = Vector{Cint}(undef, n_attypes)
    spg_spins = Vector{Cdouble}(undef, n_attypes)
    spg_positions = Matrix{Cdouble}(undef, 3, n_attypes)

    arbitrary_spin = false
    offset = 0
    nextnumber = 1
    mapping = Dict{Int, Any}()
    for (iatom, (el, positions)) in enumerate(atoms)
        mapping[nextnumber] = el

        # Default to zero magnetic moment unless this is a case of collinear magnetism
        for (ipos, pos) in enumerate(positions)
            # assign the same number to all elements with this position
            spg_numbers[offset + ipos] = nextnumber
            spg_positions[:, offset + ipos] .= pos

            if !isempty(magnetic_moments)
                magmom = magnetic_moments[iatom][2][ipos]
                spg_spins[offset + ipos] = magmom[3]
                !iszero(magmom[1:2]) && (arbitrary_spin = true)
            end
        end
        offset += length(positions)
        nextnumber += 1
    end

    collinear = !isempty(magnetic_moments) && !arbitrary_spin && !all(iszero, spg_spins)
    (positions=spg_positions, numbers=spg_numbers, spins=spg_spins,
     mapping=mapping, collinear=collinear)
end

@timing function spglib_get_symmetry(lattice, atoms, magnetic_moments=[]; tol_symmetry=1e-5)
    lattice = Matrix{Float64}(lattice)  # spglib operates in double precision

    if isempty(atoms)
        # spglib doesn't like no atoms, so we default to
        # no symmetries (even though there are lots)
        return [Mat3{Int}(I)], [Vec3(zeros(3))]
    end

    # Ask spglib for symmetry operations and for irreducible mesh
    spg_positions, spg_numbers, spg_spins, _, collinear = spglib_atoms(atoms, magnetic_moments)

    # Maximal number of symmetry operations spglib searches for
    max_ops = max(384, 50 * length(spg_numbers))
    spg_rotations    = Array{Cint}(undef, 3, 3, max_ops)
    spg_translations = Array{Cdouble}(undef, 3, max_ops)
    if collinear
        spg_equivalent_atoms = Array{Cint}(undef, max_ops)
        spg_n_ops = ccall((:spg_get_symmetry_with_collinear_spin, SPGLIB), Cint,
                          (Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint}, Cint, Ptr{Cdouble},
                           Ptr{Cdouble}, Ptr{Cint}, Ptr{Cdouble}, Cint, Cdouble),
                          spg_rotations, spg_translations, spg_equivalent_atoms, max_ops, copy(lattice'),
                          spg_positions, spg_numbers, spg_spins, Cint(length(spg_numbers)), tol_symmetry)
    else
        spg_n_ops = ccall((:spg_get_symmetry, SPGLIB), Cint,
            (Ptr{Cint}, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Cint, Cdouble),
            spg_rotations, spg_translations, max_ops, copy(lattice'), spg_positions, spg_numbers,
            Cint(length(spg_numbers)), tol_symmetry)
    end

    # If spglib does not find symmetries give an error
    if spg_n_ops == 0
        err_message = spglib_get_error_message()
        error("spglib failed to get the symmetries. Check your lattice, use a " *
              "uniform BZ mesh or disable symmetries. Spglib reported : " * err_message)
    end

    # Note: Transposes are performed to convert between spglib row-major to julia column-major
    Stildes = [Mat3{Int}(spg_rotations[:, :, i]') for i in 1:spg_n_ops]
    τtildes = [rationalize.(Vec3{Float64}(spg_translations[:, i]), tol=tol_symmetry)
               for i in 1:spg_n_ops]

    # Checks: (A Stilde A^{-1}) is unitary
    for Stilde in Stildes
        Scart = lattice * Stilde * inv(lattice)  # Form S in cartesian coords
        if maximum(abs, Scart'Scart - I) > tol_symmetry
            error("spglib returned bad symmetries: Non-unitary rotation matrix.")
        end
    end

    # Check (Stilde, τtilde) maps atoms to equivalent atoms in the lattice
    for (Stilde, τtilde) in zip(Stildes, τtildes)
        for (elem, positions) in atoms
            for coord in positions
                diffs = [rationalize.(Stilde * coord + τtilde - pos, tol=5*tol_symmetry)
                         for pos in positions]

                # If all elements of a difference in diffs is integer, then
                # Stilde * coord + τtilde and pos are equivalent lattice positions
                if !any(all(isinteger, d) for d in diffs)
                    error("spglib returned bad symmetries: Cannot map the atom at position " *
                          "$coord to another atom of the same element under the symmetry " *
                          "operation (Stilde, τtilde):\n" *
                          "($Stilde, $τtilde)")
                end
            end
        end
    end

    Stildes, τtildes
end

function spglib_get_stabilized_reciprocal_mesh(kgrid_size, rotations::Vector;
                                               is_shift=Vec3(0, 0, 0),
                                               is_time_reversal=false,
                                               qpoints=[Vec3(0.0, 0.0, 0.0)])
    spg_rotations = cat([copy(Cint.(S')) for S in rotations]..., dims=3)
    nkpt = prod(kgrid_size)
    mapping = Vector{Cint}(undef, nkpt)
    grid_address = Matrix{Cint}(undef, 3, nkpt)

    nrot = length(rotations)
    n_kpts = ccall((:spg_get_stabilized_reciprocal_mesh, SPGLIB), Cint,
      (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint, Cint, Ptr{Cint}, Cint, Ptr{Cdouble}),
      grid_address, mapping, [Cint.(kgrid_size)...], [Cint.(is_shift)...], Cint(is_time_reversal),
      Cint(nrot), spg_rotations, Cint(length(qpoints)), Vec3{Float64}.(qpoints))

    return n_kpts, Int.(mapping), [Vec3{Int}(grid_address[:, i]) for i in 1:nkpt]
end

# Returns crystallographic conventional cell if primitive is false, else the primitive
# cell in the convention by spglib
function spglib_standardize_cell(lattice::AbstractArray{T}, atoms; correct_symmetry=true,
                                 primitive=false, tol_symmetry=1e-5) where {T}
    # Convert lattice and atoms to spglib and keep the mapping between our atoms
    spg_lattice = copy(Matrix{Float64}(lattice)')
    # and spglibs atoms
    spg_positions, spg_numbers, spg_spins, atommapping = spglib_atoms(atoms)

    # TODO This drops magnetic moments!
    # TODO What about time-reversal symmetry?

    # Ask spglib to standardize the cell (i.e. find a cell, which fits the spglib conventions)
    num_atoms = ccall((:spg_standardize_cell, SPGLIB), Cint,
      (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Cint, Cint, Cint, Cdouble),
      spg_lattice, spg_positions, spg_numbers, length(spg_numbers), Cint(primitive),
      Cint(!correct_symmetry), tol_symmetry)
    spg_lattice = copy(spg_lattice')

    newatoms = [(atommapping[iatom]
                 => T.(spg_positions[findall(isequal(iatom), spg_numbers), :]))
                for iatom in unique(spg_numbers)]
    Matrix{T}(spg_lattice), newatoms
end

# TODO merge this function into spglib_standardize_cell
"""
Returns crystallographic conventional cell according to the International Table of
Crystallography Vol A (ITA) in case `to_primitive=false`. If `to_primitive=true`
the primitive lattice is returned in the convention of the reference work of
Cracknell, Davies, Miller, and Love (CDML). Of note this has minor differences to
the primitive setting choice made in the ITA.
"""
function get_spglib_lattice(model; to_primitive=false)
    # TODO This drops magnetic moments!
    # TODO For time-reversal symmetry see the discussion in PR 496.
    #      https://github.com/JuliaMolSim/DFTK.jl/pull/496/files#r725203554
    #      Essentially this does not influence the standardisation,
    #      but it only influences the kpath.
    spg_positions, spg_numbers, _ = spglib_atoms(model.atoms)
    structure = Spglib.Cell(transpose(model.lattice), spg_positions, spg_numbers)
    spglib_lattice = Spglib.standardize_cell(structure, to_primitive=to_primitive).lattice
    Matrix(copy(spglib_lattice'))
end

function spglib_spacegroup_number(model)
    # Get spacegroup number according to International Tables for Crystallography (ITA)
    # TODO Time-reversal symmetry disabled? (not yet available in DFTK)
    # TODO Are magnetic moments passed?
    spg_positions, spg_numbers, _ = spglib_atoms(model.atoms)
    structure = Spglib.Cell(transpose(model.lattice), spg_positions, spg_numbers)
    spacegroup_number = Spglib.get_spacegroup_number(structure)
    spacegroup_number
end


text:
with_repo_access() do test_repo, path
    context("test packfile object exists") do
        @test exists(test_repo, Oid("41bc8c69075bbdb46c5c6f0566cc8cc5b46e8bd9"))
        @test exists(test_repo, Oid("f82a8eb4cb20e88d1030fd10d89286215a715396"))
    end

    context("test read packed object") do
        rawobj = read(test_repo, Oid("41bc8c69075bbdb46c5c6f0566cc8cc5b46e8bd9"))
        @test match(r"tree f82a8eb4cb20e88d1030fd10d89286215a715396", data(rawobj)) != nothing
        @test length(rawobj) == 230
        @test isa(rawobj, OdbObject{GitCommit})
    end

    context("test read packed header") do
         h = read_header(test_repo, Oid("41bc8c69075bbdb46c5c6f0566cc8cc5b46e8bd9"))
         @test h[:nbytes] == 230
         @test h[:type] == GitCommit
    end
end