text:
module OperationsTest

import Random: randstring
import LibGit2
using Test
using UUIDs

using Pkg3
using Pkg3.Types

import Random: randstring
import LibGit2

include("utils.jl")

const TEST_PKG = (name = "Example", uuid = UUID("7876af07-990d-54b4-ab0e-23690620f79a"))

temp_pkg_dir() do project_path
    @testset "simple add and remove with preview" begin
        Pkg3.init(project_path)
        Pkg3.add(TEST_PKG.name; preview = true)
        @test !isinstalled(TEST_PKG)
        Pkg3.add(TEST_PKG.name)
        @test isinstalled(TEST_PKG)
        @eval import $(Symbol(TEST_PKG.name))
        Pkg3.rm(TEST_PKG.name; preview = true)
        @test isinstalled(TEST_PKG)
        Pkg3.rm(TEST_PKG.name)
        @test !isinstalled(TEST_PKG)
    end

    @testset "package with wrong UUID" begin
        @test_throws ResolverError Pkg3.add(PackageSpec(TEST_PKG.name, UUID(UInt128(1))))
    end

    @testset "adding and upgrading different versions" begin
        # VersionNumber
        Pkg3.add(PackageSpec(TEST_PKG.name, v"0.3"))
        @test Pkg3.installed()[TEST_PKG.name] == v"0.3"
        Pkg3.add(PackageSpec(TEST_PKG.name, v"0.3.1"))
        @test Pkg3.installed()[TEST_PKG.name] == v"0.3.1"
        Pkg3.rm(TEST_PKG.name)

        # VersionRange
        Pkg3.add(PackageSpec(TEST_PKG.name, VersionSpec(VersionRange("0.3.0-0.3.2"))))
        @test Pkg3.installed()[TEST_PKG.name] == v"0.3.2"
        Pkg3.up(; level = UpgradeLevel(:patch))
        @test Pkg3.installed()[TEST_PKG.name] == v"0.3.3"
        Pkg3.up(; level = UpgradeLevel(:minor))
        @test Pkg3.installed()[TEST_PKG.name].minor != 3
        Pkg3.rm(TEST_PKG.name)
    end

    @testset "testing" begin
        # TODO: Check that preview = true doesn't actually execute the test
        # TODO: Test-only dependencies
        Pkg3.add(TEST_PKG.name)
        Pkg3.test(TEST_PKG.name; coverage=true)
        pkgdir = Base.locate_package(Base.PkgId(TEST_PKG.uuid, TEST_PKG.name))
        # No coverage files being generated?
        @test_broken TEST_PKG.name * ".cov" in readdir(pkgdir)
        Pkg3.rm(TEST_PKG.name)
    end

    @testset "pinning / freeing" begin
        Pkg3.add(TEST_PKG.name)
        old_v = Pkg3.installed()[TEST_PKG.name]
        Pkg3.pin(PackageSpec(TEST_PKG.name, v"0.2"))
        @test Pkg3.installed()[TEST_PKG.name].minor == 2
        Pkg3.up(TEST_PKG.name)
        @test Pkg3.installed()[TEST_PKG.name].minor == 2
        Pkg3.free(TEST_PKG.name)
        Pkg3.up()
        @test Pkg3.installed()[TEST_PKG.name] == old_v
        Pkg3.rm(TEST_PKG.name)
    end

    @testset "develop / freeing" begin
        Pkg3.add(TEST_PKG.name)
        old_v = Pkg3.installed()[TEST_PKG.name]
        Pkg3.rm(TEST_PKG.name)
        mktempdir() do devdir
            withenv("JULIA_PKG_DEVDIR" => devdir) do
                Pkg3.REPLMode.pkgstr("develop $(TEST_PKG.name)")
                @test isinstalled(TEST_PKG)
                @test Pkg3.installed()[TEST_PKG.name] > old_v
                test_pkg_main_file = joinpath(devdir, TEST_PKG.name, "src", TEST_PKG.name * ".jl")
                @test isfile(test_pkg_main_file)
                # Pkg3 #152
                write(test_pkg_main_file,
                    """
                    module Example
                        export hello, domath
                        const example2path = joinpath(@__DIR__, "..", "deps", "deps.jl")
                        if !isfile(example2path)
                            error("Example is not installed correctly")
                        end
                        hello(who::String) = "Hello, \$who"
                        domath(x::Number) = x + 5
                    end
                    """)
                mkpath(joinpath(devdir, TEST_PKG.name, "deps"))
                write(joinpath(devdir, TEST_PKG.name, "deps", "build.jl"),
                    """
                    touch("deps.jl")
                    """
                )
                Pkg3.build(TEST_PKG.name)
                @test isfile(joinpath(devdir, TEST_PKG.name, "deps", "deps.jl"))
                Pkg3.test(TEST_PKG.name)
                Pkg3.free(TEST_PKG.name)
                @test Pkg3.installed()[TEST_PKG.name] == old_v
            end
        end
    end

    @testset "stdlibs as direct dependency" begin
        uuid_pkg = (name = "CRC32c", uuid = UUID("8bf52ea8-c179-5cab-976a-9e18b702a9bc"))
        Pkg3.add("CRC32c")
        @test haskey(Pkg3.installed(), uuid_pkg.name)
        Pkg3.up()
        Pkg3.test("CRC32c")
        Pkg3.rm("CRC32c")
    end

    @testset "package name in resolver errors" begin
        try
            Pkg3.add([PackageSpec(TEST_PKG.name, VersionSpec(v"55"))])
        catch e
            @test occursin(TEST_PKG.name, sprint(showerror, e))
        end
    end

    @testset "check logging" begin
        usage = Pkg3.TOML.parse(String(read(joinpath(Pkg3.logdir(), "manifest_usage.toml"))))
        @test any(x -> startswith(x, joinpath(project_path, "Manifest.toml")), keys(usage))
    end

    @testset "adding nonexisting packages" begin
        nonexisting_pkg = randstring(14)
        @test_throws CommandError Pkg3.add(nonexisting_pkg)
        @test_throws CommandError Pkg3.up(nonexisting_pkg)
    end

    Pkg3.rm(TEST_PKG.name)

    @testset "legacy CI script" begin
        mktempdir() do dir
            LibGit2.with(LibGit2.clone("https://github.com/JuliaLang/Example.jl", joinpath(dir, "Example.jl"))) do r
                cd(joinpath(dir, "Example.jl")) do
                    let Pkg = Pkg3
                        Pkg.clone(pwd())
                        Pkg.build("Example")
                        Pkg.test("Example"; coverage=true)
                        @test isfile(Pkg.dir("Example", "src", "Example.jl"))
                    end
                end
            end
        end
    end

    @testset "up in Project without manifest" begin
        mktempdir() do dir
            cp(joinpath(@__DIR__, "test_packages", "UnregisteredWithProject"), joinpath(dir, "UnregisteredWithProject"))
            cd(joinpath(dir, "UnregisteredWithProject")) do
                try
                    pushfirst!(LOAD_PATH, Base.parse_load_path("@"))
                    Pkg3.up()
                    @test haskey(Pkg3.installed(), "Example")
                finally
                    popfirst!(LOAD_PATH)
                end
            end
        end
    end
end

temp_pkg_dir() do project_path
    @testset "libgit2 downloads" begin
        Pkg3.add(TEST_PKG.name; use_libgit2_for_all_downloads=true)
        @test haskey(Pkg3.installed(), TEST_PKG.name)
        Pkg3.rm(TEST_PKG.name)
    end
end

include("repl.jl")

end # module


text:
#julia -L rollingHarness.jl test_Rolling.jl outPathStub d datastub
#passed arguments
	#ARGS[1] is partial path for output.
	#ARGS[2] is d must be one of 20, 50, 1000
	#ARGS[3] specifies the data stub:  AdjSales_NoWeekends_ShuffleWithinCol

using Distributed

const spath = ARGS[1]
const param_path = ARGS[3] # e.g. "AdjSales_NoWeekends"
const d = parse(Int, ARGS[2])
const s = .95

@assert d in [20, 50, 1000] "Only bins of size 20, 50, 1000 currently supported"

K_grid = vcat(round.(Int, 2 .^(3:.25:10)), 1115)
N_grid = [10, 20, 40]

outPath = "$(spath)_RossRolling_$(s)__$(d)_$(param_path)"

#First read in the data and parse it appropriately
ps_full = readdlm("../RossmanKaggleData/CleanedData/ps_full$(d).csv", ',')
supp_full = readdlm("../RossmanKaggleData/CleanedData/support$(d).csv", ',')

@assert minimum(ps_full) >=0 "ps_full has negative entries"

#Load shuffled data
tdata = readdlm("../RossmanKaggleData/CleanedData/$(param_path)_Binned$(d).csv", ',')

binned_data = tdata[2:end, 2:end]  #drop column header = stores, drop row header = dates
dates= tdata[2:end, 1] #keep track of dates for fun

#currently run single-threaded for ease
@time rollingTest(K_grid, supp_full, ps_full, binned_data, dates, outPath, N_grid, s, onlySAA=false, numTestDays=10)

ext: jl


text:
# todo: ^, fma, muladd, div, fld, cld, rem, mod, mod1, fld1

Base.typemin(::Type{Infinite}) = NegInf
Base.typemax(::Type{Infinite}) = PosInf

Base.:+(x::Infinite) = x
Base.:-(x::Infinite) = Infinite(!x.signbit)

Base.:+(x::Infinite, y::Infinite) = signbit(x)==signbit(y) ? x : throw(InfMinusInfError())
Base.:-(x::Infinite, y::Infinite) = x + -y

Base.:*(x::Infinite, y::Infinite) = Infinite(x.signbit ⊻ y.signbit)
Base.:/(x::Infinite, y::Infinite) = throw(DivideError())

Base.abs(x::Infinite) = PosInf

Base.:(//)(x::Infinite, y::Infinite) = throw(DivideError())
Base.:(//)(x::Infinite, y::Real) = //(promote(x,y)...)
Base.:(//)(x::Real, y::Infinite) = //(promote(x,y)...)


text:
# matmul.jl: Everything to do with dense matrix multiplication

# multiply by diagonal matrix as vector
function scale!(C::Matrix, A::Matrix, b::Vector)
    m, n = size(A)
    n==length(b) || throw(DimensionMismatch(""))
    for j = 1:n
        bj = b[j]
        for i = 1:m
            C[i,j] = A[i,j]*bj
        end
    end
    C
end

function scale!(C::Matrix, b::Vector, A::Matrix)
    m, n = size(A)
    m==length(b) || throw(DimensionMismatch(""))
    for j=1:n, i=1:m
        C[i,j] = A[i,j]*b[i]
    end
    C
end
scale(A::Matrix, b::Vector) = scale!(similar(A, promote_type(eltype(A),eltype(b))), A, b)
scale(b::Vector, A::Matrix) = scale!(similar(b, promote_type(eltype(A),eltype(b)), size(A)), b, A)

# Dot products

dot{T<:BlasReal}(x::Vector{T}, y::Vector{T}) = BLAS.dot(x, y)
dot{T<:BlasComplex}(x::Vector{T}, y::Vector{T}) = BLAS.dotc(x, y)
function dot{T<:BlasReal, TI<:Integer}(x::Vector{T}, rx::Union(UnitRange{TI},Range{TI}), y::Vector{T}, ry::Union(UnitRange{TI},Range{TI}))
    length(rx)==length(ry) || throw(DimensionMismatch(""))
    if minimum(rx) < 1 || maximum(rx) > length(x) || minimum(ry) < 1 || maximum(ry) > length(y)
        throw(BoundsError())
    end
    BLAS.dot(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry))
end
function dot{T<:BlasComplex, TI<:Integer}(x::Vector{T}, rx::Union(UnitRange{TI},Range{TI}), y::Vector{T}, ry::Union(UnitRange{TI},Range{TI}))
    length(rx)==length(ry) || throw(DimensionMismatch(""))
    if minimum(rx) < 1 || maximum(rx) > length(x) || minimum(ry) < 1 || maximum(ry) > length(y)
        throw(BoundsError())
    end
    BLAS.dotc(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry))
end
function dot(x::AbstractVector, y::AbstractVector)
    length(x)==length(y) || throw(DimensionMismatch(""))
    s = zero(eltype(x))*zero(eltype(y))
    for i=1:length(x)
        s += conj(x[i])*y[i]
    end
    s
end
dot(x::Number, y::Number) = conj(x) * y
Ac_mul_B(x::Vector, y::Vector) = [dot(x, y)]
At_mul_B{T<:Real}(x::Vector{T}, y::Vector{T}) = [dot(x, y)]
At_mul_B{T<:BlasComplex}(x::Vector{T}, y::Vector{T}) = [BLAS.dotu(x, y)]

# Matrix-vector multiplication

*{T<:Union(Float32,Integer,Rational)}(A::StridedMatrix{Float64}, X::StridedVector{T}) = A*convert(Vector{eltype(A)},X)
*{T<:Union(Float32,Complex64,Integer,Rational)}(A::StridedMatrix{Complex128}, X::StridedVector{T}) = A*convert(Vector{eltype(A)},X)
function *{T<:BlasFloat}(A::StridedMatrix{T}, X::StridedVector{T})
    gemv(similar(A, size(A,1)), 'N', A, X)
end

A_mul_B!{T<:BlasFloat}(y::StridedVector{T}, A::StridedMatrix{T}, x::StridedVector{T}) = gemv(y, 'N', A, x)
A_mul_B!(y::StridedVector, A::StridedMatrix, x::StridedVector) = generic_matvecmul(y, 'N', A, x)

At_mul_B{T<:Union(Float32,Integer,Rational)}(A::StridedMatrix{Float64}, X::StridedVector{T}) = At_mul_B(A,convert(Vector{eltype(A)},X))
At_mul_B{T<:Union(Float32,Complex64,Integer,Rational)}(A::StridedMatrix{Complex128}, X::StridedVector{T}) = At_mul_B(A,convert(Vector{eltype(A)},X))
function At_mul_B{T<:BlasFloat}(A::StridedMatrix{T}, x::StridedVector{T})
    gemv(similar(A, size(A, 2)), 'T', A, x)
end

At_mul_B!{T<:BlasFloat}(y::StridedVector{T}, A::StridedMatrix{T}, x::StridedVector{T}) = gemv(y, 'T', A, x)
At_mul_B!(y::StridedVector, A::StridedMatrix, x::StridedVector) = generic_matvecmul(y, 'T', A, x)

Ac_mul_B{T<:Union(Float32,Integer,Rational)}(A::StridedMatrix{Float64}, X::StridedVector{T}) = Ac_mul_B(A,convert(Vector{eltype(A)},X))
Ac_mul_B{T<:Union(Float32,Complex64,Integer,Rational)}(A::StridedMatrix{Complex128}, X::StridedVector{T}) = Ac_mul_B(A,convert(Vector{eltype(A)},X))
Ac_mul_B{T<:BlasReal}(A::StridedMatrix{T}, x::StridedVector{T}) = At_mul_B(A, x)
function Ac_mul_B{T<:BlasComplex}(A::StridedMatrix{T}, x::StridedVector{T})
    gemv(similar(A, size(A, 2)), 'C', A, x)
end

Ac_mul_B!{T<:BlasReal}(y::StridedVector{T}, A::StridedMatrix{T}, x::StridedVector{T}) = At_mul_B!(y, A, x)
Ac_mul_B!{T<:BlasComplex}(y::StridedVector{T}, A::StridedMatrix{T}, x::StridedVector{T}) = gemv(y, 'C', A, x)
Ac_mul_B!(y::StridedVector, A::StridedMatrix, x::StridedVector) = generic_matvecmul(y, 'C', A, x)

# Matrix-matrix multiplication

(*){T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}) = gemm_wrapper('N', 'N', A, B)
A_mul_B!{T<:BlasFloat}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) = gemm_wrapper(C, 'N', 'N', A, B)
A_mul_B!{T,S,R}(C::StridedMatrix{R}, A::StridedMatrix{T}, B::StridedMatrix{S}) = generic_matmatmul(C, 'N', 'N', A, B)

function At_mul_B{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T})
    is(A, B) ? syrk_wrapper('T', A) : gemm_wrapper('T', 'N', A, B)
end

At_mul_B!{T<:BlasFloat}(C::StridedMatrix{T}, A::StridedVecOrMat{T}, B::StridedMatrix{T}) = gemm_wrapper(C, 'T', 'N', A, B)
At_mul_B{T,S}(A::StridedMatrix{T}, B::StridedMatrix{S}) = generic_matmatmul('T', 'N', A, B)
At_mul_B!{T,S,R}(C::StridedMatrix{R}, A::StridedVecOrMat{T}, B::StridedMatrix{S}) = generic_matmatmul(C, 'T', 'N', A, B)

function A_mul_Bt{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T})
    is(A, B) ? syrk_wrapper('N', A) : gemm_wrapper('N', 'T', A, B)
end

A_mul_Bt!{T<:BlasFloat}(C::StridedVecOrMat{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) = gemm_wrapper(C, 'N', 'T', A, B)
A_mul_Bt{T,S}(A::StridedMatrix{T}, B::StridedMatrix{S}) = generic_matmatmul('N', 'T', A, B)
A_mul_Bt!{T,S,R}(C::StridedVecOrMat{R}, A::StridedMatrix{T}, B::StridedMatrix{S}) = generic_matmatmul(C, 'N', 'T', A, B)

At_mul_Bt{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}) = gemm_wrapper('T', 'T', A, B)
At_mul_Bt!{T<:BlasFloat}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) = gemm_wrapper(C, 'T', 'T', A, B)
At_mul_Bt{T,S}(A::StridedMatrix{T}, B::StridedMatrix{S}) = generic_matmatmul('T', 'T', A, B)
At_mul_Bt!{T,S,R}(C::StridedMatrix{R}, A::StridedMatrix{T}, B::StridedMatrix{S}) = generic_matmatmul(C, 'T', 'T', A, B)

Ac_mul_B{T<:BlasReal}(A::StridedMatrix{T}, B::StridedMatrix{T}) = At_mul_B(A, B)
Ac_mul_B!{T<:BlasReal}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) = At_mul_B!(C, A, B)

function Ac_mul_B{T<:BlasComplex}(A::StridedMatrix{T}, B::StridedMatrix{T})
    is(A, B) ? herk_wrapper('C', A) : gemm_wrapper('C', 'N', A, B)
end

Ac_mul_B!{T<:BlasComplex}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) = gemm_wrapper('C', 'N', A, B)
Ac_mul_B{T,S}(A::StridedMatrix{T}, B::StridedMatrix{S}) = generic_matmatmul('C', 'N', A, B)
Ac_mul_B!{T,S,R}(C::StridedMatrix{R}, A::StridedMatrix{T}, B::StridedMatrix{S}) = generic_matmatmul(C, 'C', 'N', A, B)

A_mul_Bc{T<:BlasReal}(A::StridedMatrix{T}, B::StridedMatrix{T}) = A_mul_Bt(A, B)
A_mul_Bc!{T<:BlasReal}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) = A_mul_Bt!(C, A, B)
function A_mul_Bc{T<:BlasComplex}(A::StridedMatrix{T}, B::StridedMatrix{T})
    is(A, B) ? herk_wrapper('N', A) : gemm_wrapper('N', 'C', A, B)
end
A_mul_Bc!{T<:BlasComplex}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) = gemm_wrapper(C, 'N', 'C', A, B)
A_mul_Bc{T,S}(A::StridedMatrix{T}, B::StridedMatrix{S}) = generic_matmatmul('N', 'C', A, B)
A_mul_Bc!{T,S,R}(C::StridedMatrix{R}, A::StridedMatrix{T}, B::StridedMatrix{S}) = generic_matmatmul(C, 'N', 'C', A, B)

Ac_mul_Bc{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}) = gemm_wrapper('C', 'C', A, B)
Ac_mul_Bc!{T<:BlasFloat}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) = gemm_wrapper(C, 'C', 'C', A, B)
Ac_mul_Bt{T,S}(A::StridedMatrix{T}, B::StridedMatrix{S}) = generic_matmatmul('C', 'C', A, B)
Ac_mul_Bt!{T,S,R}(C::StridedMatrix{R}, A::StridedMatrix{T}, B::StridedMatrix{S}) = generic_matmatmul(C, 'C', 'C', A, B)

# Supporting functions for matrix multiplication

function copytri!(A::StridedMatrix, uplo::Char, conjugate::Bool=false)
    n = chksquare(A)
    @chkuplo
    if uplo == 'U'
        for i = 1:(n-1), j = (i+1):n
            A[j,i] = conjugate ? conj(A[i,j]) : A[i,j]
        end
    elseif uplo == 'L'
        for i = 1:(n-1), j = (i+1):n
            A[i,j] = conjugate ? conj(A[j,i]) : A[j,i]
        end
    else
        throw(ArgumentError("second argument must be 'U' or 'L'"))
    end
    A
end

function gemv{T<:BlasFloat}(y::StridedVector{T}, tA::Char, A::StridedMatrix{T}, x::StridedVector{T})
    stride(A, 1)==1 || return generic_matvecmul(y, tA, A, x)

    if tA != 'N'
        (nA, mA) = size(A)
    else
        (mA, nA) = size(A)
    end

    nA==length(x) || throw(DimensionMismatch(""))
    mA==length(y) || throw(DimensionMismatch(""))
    mA == 0 && return zeros(T, 0)
    nA == 0 && return zeros(T, mA)
    return BLAS.gemv!(tA, one(T), A, x, zero(T), y)
end

function syrk_wrapper{T<:BlasFloat}(tA::Char, A::StridedMatrix{T})
    if tA == 'T'
        (nA, mA) = size(A)
        tAt = 'N'
    else
        (mA, nA) = size(A)
        tAt = 'T'
    end

    if mA == 0 || nA == 0; return zeros(T, mA, mA); end
    if mA == 2 && nA == 2; return matmul2x2(tA,tAt,A,A); end
    if mA == 3 && nA == 3; return matmul3x3(tA,tAt,A,A); end

    stride(A, 1) == 1 || (return generic_matmatmul(tA, tAt, A, A))
    copytri!(BLAS.syrk('U', tA, one(T), A), 'U')
end

function herk_wrapper{T<:BlasFloat}(tA::Char, A::StridedMatrix{T})
    if tA == 'C'
        (nA, mA) = size(A)
        tAt = 'N'
    else
        (mA, nA) = size(A)
        tAt = 'C'
    end

    if mA == 2 && nA == 2; return matmul2x2(tA,tAt,A,A); end
    if mA == 3 && nA == 3; return matmul3x3(tA,tAt,A,A); end

    stride(A, 1) == 1 || (return generic_matmatmul(tA, tAt, A, A))

    # Result array does not need to be initialized as long as beta==0
    #    C = Array(T, mA, mA)

    copytri!(BLAS.herk('U', tA, one(T), A), 'U', true)
end

function gemm_wrapper{T<:BlasFloat}(tA::Char, tB::Char,
                             A::StridedVecOrMat{T},
                             B::StridedMatrix{T})
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)
    C = similar(B, T, mA, nB)
    gemm_wrapper(C, tA, tB, A, B)
end

function gemm_wrapper{T<:BlasFloat}(C::StridedVecOrMat{T}, tA::Char, tB::Char,
                             A::StridedVecOrMat{T},
                             B::StridedMatrix{T})
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)

    nA==mB || throw(DimensionMismatch("*"))

    if mA == 0 || nA == 0 || nB == 0; return zeros(T, mA, nB); end
    if mA == 2 && nA == 2 && nB == 2; return matmul2x2(C,tA,tB,A,B); end
    if mA == 3 && nA == 3 && nB == 3; return matmul3x3(C,tA,tB,A,B); end

    stride(A, 1)==stride(B, 1)==1 || (return generic_matmatmul(C, tA, tB, A, B))
    BLAS.gemm!(tA, tB, one(T), A, B, zero(T), C)
end

# blas.jl defines matmul for floats; other integer and mixed precision
# cases are handled here

lapack_size(t::Char, M::AbstractVecOrMat) = (size(M, t=='N' ? 1:2), size(M, t=='N' ? 2:1))

function copy!{R,S}(B::AbstractMatrix{R}, ir_dest::UnitRange{Int}, jr_dest::UnitRange{Int}, tM::Char, M::AbstractMatrix{S}, ir_src::UnitRange{Int}, jr_src::UnitRange{Int})
    if tM == 'N'
        copy!(B, ir_dest, jr_dest, M, ir_src, jr_src)
    else
        Base.copy_transpose!(B, ir_dest, jr_dest, M, jr_src, ir_src)
        tM == 'C' && conj!(B)
    end
end

function copy_transpose!{R,S}(B::AbstractMatrix{R}, ir_dest::UnitRange{Int}, jr_dest::UnitRange{Int}, tM::Char, M::AbstractVecOrMat{S}, ir_src::UnitRange{Int}, jr_src::UnitRange{Int})
    if tM == 'N'
        Base.copy_transpose!(B, ir_dest, jr_dest, M, ir_src, jr_src)
    else
        copy!(B, ir_dest, jr_dest, M, jr_src, ir_src)
        tM == 'C' && conj!(B)
    end
end

# TODO: It will be faster for large matrices to convert to float,
# call BLAS, and convert back to required type.

# NOTE: the generic version is also called as fallback for
#       strides != 1 cases in libalg_blas.jl
(*){T,S}(A::AbstractMatrix{T}, B::AbstractVector{S}) = generic_matvecmul('N', A, B)

arithtype(T) = T
arithtype(::Type{Bool}) = Int

function generic_matvecmul{T,S}(tA::Char, A::AbstractMatrix{T}, B::AbstractVector{S})
    C = similar(B, promote_type(arithtype(T),arithtype(S)), size(A, tA=='N' ? 1 : 2))
    generic_matvecmul(C, tA, A, B)
end

function generic_matvecmul{T,S,R}(C::AbstractVector{R}, tA, A::AbstractMatrix{T}, B::AbstractVector{S})
    mB = length(B)
    mA, nA = lapack_size(tA, A)
    mB==nA || throw(DimensionMismatch("*"))
    mA==length(C) || throw(DimensionMismatch("*"))
    z = zero(R)

    Astride = size(A, 1)

    if tA == 'T'  # fastest case
        for k = 1:mA
            aoffs = (k-1)*Astride
            s = z
            for i = 1:nA
                s += A[aoffs+i] * B[i]
            end
            C[k] = s
        end
    elseif tA == 'C'
        for k = 1:mA
            aoffs = (k-1)*Astride
            s = z
            for i = 1:nA
                s += conj(A[aoffs+i]) * B[i]
            end
            C[k] = s
        end
        else # tA == 'N'
        fill!(C, z)
        for k = 1:mB
            aoffs = (k-1)*Astride
            b = B[k]
            for i = 1:mA
                C[i] += A[aoffs+i] * b
            end
        end
    end
    C
end

(*){T,S}(A::AbstractVector{S}, B::AbstractMatrix{T}) = reshape(A,length(A),1)*B

# NOTE: the generic version is also called as fallback for strides != 1 cases
#       in libalg_blas.jl
(*){T,S}(A::AbstractVecOrMat{T}, B::AbstractMatrix{S}) = generic_matmatmul('N', 'N', A, B)

function generic_matmatmul{T,S}(tA, tB, A::AbstractVecOrMat{T}, B::AbstractMatrix{S})
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)
    C = similar(B, promote_type(arithtype(T),arithtype(S)), mA, nB)
    generic_matmatmul(C, tA, tB, A, B)
end

const tilebufsize = 10800  # Approximately 32k/3
const Abuf = Array(Uint8, tilebufsize)
const Bbuf = Array(Uint8, tilebufsize)
const Cbuf = Array(Uint8, tilebufsize)

function generic_matmatmul{T,S,R}(C::AbstractVecOrMat{R}, tA, tB, A::AbstractVecOrMat{T}, B::AbstractMatrix{S})
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)
    mB==nA || throw(DimensionMismatch("*"))
    if size(C,1) != mA || size(C,2) != nB; throw(DimensionMismatch("*")); end

    if mA == nA == nB == 2; return matmul2x2(C, tA, tB, A, B); end
    if mA == nA == nB == 3; return matmul3x3(C, tA, tB, A, B); end

    @inbounds begin
    if isbits(R)
        tile_size = int(ifloor(sqrt(tilebufsize/sizeof(R))))
        sz = (tile_size, tile_size)
        Atile = pointer_to_array(convert(Ptr{R}, pointer(Abuf)), sz)
        Btile = pointer_to_array(convert(Ptr{R}, pointer(Bbuf)), sz)

        z = zero(R)

        if mA < tile_size && nA < tile_size && nB < tile_size
            Base.copy_transpose!(Atile, 1:nA, 1:mA, tA, A, 1:mA, 1:nA)
            copy!(Btile, 1:mB, 1:nB, tB, B, 1:mB, 1:nB)
            for j = 1:nB
                boff = (j-1)*tile_size
                for i = 1:mA
                    aoff = (i-1)*tile_size
                    s = z
                    for k = 1:nA
                        s += Atile[aoff+k] * Btile[boff+k]
                    end
                    C[i,j] = s
                end
            end
        else
            Ctile = pointer_to_array(convert(Ptr{R}, pointer(Cbuf)), sz)
            for jb = 1:tile_size:nB
                jlim = min(jb+tile_size-1,nB)
                jlen = jlim-jb+1
                for ib = 1:tile_size:mA
                    ilim = min(ib+tile_size-1,mA)
                    ilen = ilim-ib+1
                    fill!(Ctile, z)
                    for kb = 1:tile_size:nA
                        klim = min(kb+tile_size-1,mB)
                        klen = klim-kb+1
                        Base.copy_transpose!(Atile, 1:klen, 1:ilen, tA, A, ib:ilim, kb:klim)
                        copy!(Btile, 1:klen, 1:jlen, tB, B, kb:klim, jb:jlim)
                        for j=1:jlen
                            bcoff = (j-1)*tile_size
                            for i = 1:ilen
                                aoff = (i-1)*tile_size
                                s = z
                                for k = 1:klen
                                    s += Atile[aoff+k] * Btile[bcoff+k]
                                end
                                Ctile[bcoff+i] += s
                            end
                        end
                    end
                    copy!(C, ib:ilim, jb:jlim, Ctile, 1:ilen, 1:jlen)
                end
            end
        end
    else
        # Multiplication for non-plain-data uses the naive algorithm
        if tA == 'N'
            if tB == 'N'
                for i = 1:mA, j = 1:nB
                    Ctmp = A[i, 1]*B[1, j]
                    for k = 2:nA
                        Ctmp += A[i, k]*B[k, j]
                    end
                    C[i,j] = Ctmp
                end
            elseif tB == 'T'
                for i = 1:mA, j = 1:nB
                    Ctmp = A[i, 1]*B[j, 1]
                    for k = 2:nA
                        Ctmp += A[i, k]*B[j, k]
                    end
                    C[i,j] = Ctmp
                end
            else
                for i = 1:mA, j = 1:nB
                    Ctmp = A[i, 1]*conj(B[j, 1])
                    for k = 2:nA
                        Ctmp += A[i, k]*conj(B[j, k])
                    end
                    C[i,j] = Ctmp
                end
            end
        elseif tA == 'T'
            if tB == 'N'
                for i = 1:mA, j = 1:nB
                    Ctmp = A[1, i]*B[1, j]
                    for k = 2:nA
                        Ctmp += A[k, i]*B[k, j]
                    end
                    C[i,j] = Ctmp
                end
            elseif tB == 'T'
                for i = 1:mA, j = 1:nB
                    Ctmp = A[1, i]*B[j, 1]
                    for k = 2:nA
                        Ctmp += A[k, i]*B[j, k]
                    end
                    C[i,j] = Ctmp
                end
            else
                for i = 1:mA, j = 1:nB
                    Ctmp = A[1, i]*conj(B[j, 1])
                    for k = 2:nA
                        Ctmp += A[k, i]*conj(B[j, k])
                    end
                    C[i,j] = Ctmp
                end
            end
        else
            if tB == 'N'
                for i = 1:mA, j = 1:nB
                    Ctmp = conj(A[1, i])*B[1, j]
                    for k = 2:nA
                        Ctmp += conj(A[k, i])*B[k, j]
                    end
                    C[i,j] = Ctmp
                end
            elseif tB == 'T'
                for i = 1:mA, j = 1:nB
                    Ctmp = conj(A[1, i])*B[j, 1]
                    for k = 2:nA
                        Ctmp += conj(A[k, i])*B[j, k]
                    end
                    C[i,j] = Ctmp
                end
            else
                for i = 1:mA, j = 1:nB
                    Ctmp = conj(A[1, i])*conj(B[j, 1])
                    for k = 2:nA
                        Ctmp += conj(A[k, i])*conj(B[j, k])
                    end
                    C[i,j] = Ctmp
                end
            end
        end
    end
    end # @inbounds
    C
end

# multiply 2x2 matrices
function matmul2x2{T,S}(tA, tB, A::AbstractMatrix{T}, B::AbstractMatrix{S})
    matmul2x2(similar(B, promote_type(T,S), 2, 2), tA, tB, A, B)
end

function matmul2x2{T,S,R}(C::AbstractMatrix{R}, tA, tB, A::AbstractMatrix{T}, B::AbstractMatrix{S})
    if tA == 'T'
        A11 = A[1,1]; A12 = A[2,1]; A21 = A[1,2]; A22 = A[2,2]
    elseif tA == 'C'
        A11 = conj(A[1,1]); A12 = conj(A[2,1]); A21 = conj(A[1,2]); A22 = conj(A[2,2])
    else
        A11 = A[1,1]; A12 = A[1,2]; A21 = A[2,1]; A22 = A[2,2]
    end
    if tB == 'T'
        B11 = B[1,1]; B12 = B[2,1]; B21 = B[1,2]; B22 = B[2,2]
    elseif tB == 'C'
        B11 = conj(B[1,1]); B12 = conj(B[2,1]); B21 = conj(B[1,2]); B22 = conj(B[2,2])
    else
        B11 = B[1,1]; B12 = B[1,2]; B21 = B[2,1]; B22 = B[2,2]
    end
    C[1,1] = A11*B11 + A12*B21
    C[1,2] = A11*B12 + A12*B22
    C[2,1] = A21*B11 + A22*B21
    C[2,2] = A21*B12 + A22*B22
    C
end

# Multiply 3x3 matrices
function matmul3x3{T,S}(tA, tB, A::AbstractMatrix{T}, B::AbstractMatrix{S})
    matmul3x3(similar(B, promote_type(T,S), 3, 3), tA, tB, A, B)
end

function matmul3x3{T,S,R}(C::AbstractMatrix{R}, tA, tB, A::AbstractMatrix{T}, B::AbstractMatrix{S})
    if tA == 'T'
        A11 = A[1,1]; A12 = A[2,1]; A13 = A[3,1];
        A21 = A[1,2]; A22 = A[2,2]; A23 = A[3,2];
        A31 = A[1,3]; A32 = A[2,3]; A33 = A[3,3];
    elseif tA == 'C'
        A11 = conj(A[1,1]); A12 = conj(A[2,1]); A13 = conj(A[3,1]);
        A21 = conj(A[1,2]); A22 = conj(A[2,2]); A23 = conj(A[3,2]);
        A31 = conj(A[1,3]); A32 = conj(A[2,3]); A33 = conj(A[3,3]);
    else
        A11 = A[1,1]; A12 = A[1,2]; A13 = A[1,3];
        A21 = A[2,1]; A22 = A[2,2]; A23 = A[2,3];
        A31 = A[3,1]; A32 = A[3,2]; A33 = A[3,3];
    end

    if tB == 'T'
        B11 = B[1,1]; B12 = B[2,1]; B13 = B[3,1];
        B21 = B[1,2]; B22 = B[2,2]; B23 = B[3,2];
        B31 = B[1,3]; B32 = B[2,3]; B33 = B[3,3];
    elseif tB == 'C'
        B11 = conj(B[1,1]); B12 = conj(B[2,1]); B13 = conj(B[3,1]);
        B21 = conj(B[1,2]); B22 = conj(B[2,2]); B23 = conj(B[3,2]);
        B31 = conj(B[1,3]); B32 = conj(B[2,3]); B33 = conj(B[3,3]);
    else
        B11 = B[1,1]; B12 = B[1,2]; B13 = B[1,3];
        B21 = B[2,1]; B22 = B[2,2]; B23 = B[2,3];
        B31 = B[3,1]; B32 = B[3,2]; B33 = B[3,3];
    end

    C[1,1] = A11*B11 + A12*B21 + A13*B31
    C[1,2] = A11*B12 + A12*B22 + A13*B32
    C[1,3] = A11*B13 + A12*B23 + A13*B33

    C[2,1] = A21*B11 + A22*B21 + A23*B31
    C[2,2] = A21*B12 + A22*B22 + A23*B32
    C[2,3] = A21*B13 + A22*B23 + A23*B33

    C[3,1] = A31*B11 + A32*B21 + A33*B31
    C[3,2] = A31*B12 + A32*B22 + A33*B32
    C[3,3] = A31*B13 + A32*B23 + A33*B33

    C
end

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
module M

using Emporos
using Emporos.IBKR
using MacroTools

function main()
    ib = IBKR.IB()
    as, ps = open(ib, clientId=3) do ib
        as = IBKR.accounts(ib)
        ps = map(a -> IBKR.positions(ib, a), as)
        as, ps
    end
    (;ib, as, ps)
end

r = main()

end # module


text:
# Here goes all the stuff that connects this code to ExtendableGrids like
# - definition of additional adjacency types and their instantiation
# - connections between CellGeometries and FaceGeometries of their faces
# - formulas to compute Volumes and Normal Vectors etc.
# (some of this might become native in the ExtendableGrids module itself at some point)

# additional ElementGeometryTypes with parent information
#abstract type Vertex0DWithParent{Parent <: AbstractElementGeometry} <: Vertex0D end
#abstract type Vertex0DWithParents{Parent1 <: AbstractElementGeometry, Parent2 <: AbstractElementGeometry} <: Vertex0D end
#export Vertex0DWithParent, Vertex0DWithParents

#function AddParent(FEG::Type{<:Vertex0D}, CEG::Type{<:AbstractElementGeometry})
#    return Vertex0DWithParent{CEG}
#end

#abstract type Edge1DWithParent{Parent <: AbstractElementGeometry} <: Edge1D end
#abstract type Edge1DWithParents{Parent1 <: AbstractElementGeometry, Parent2 <: AbstractElementGeometry} <: Edge1D end
#export Edge1DWithParent, Edge1DWithParents

#function AddParent(FEG::Type{<:Edge1D}, CEG::Type{<:AbstractElementGeometry})
#    return Edge1DWithParent{CEG}
#end

const GridEGTypes = Vector{DataType}
const GridRegionTypes{Ti} = Union{VectorOfConstants{Ti}, Array{Ti,1}}

# additional ExtendableGrids adjacency types
abstract type CellEdges <: AbstractGridAdjacency end
abstract type CellFaces <: AbstractGridAdjacency end
abstract type CellFaceSigns <: AbstractGridAdjacency end
abstract type CellFaceOrientations <: AbstractGridAdjacency end
abstract type CellEdgeSigns <: AbstractGridAdjacency end
abstract type CellVolumes <: AbstractGridFloatArray1D end
abstract type UniqueCellGeometries <: AbstractElementGeometries end

abstract type FaceNodes <: AbstractGridAdjacency end
abstract type FaceVolumes <: AbstractGridFloatArray1D end
abstract type FaceCells <: AbstractGridAdjacency end
abstract type FaceEdges <: AbstractGridAdjacency end
abstract type FaceNormals <: AbstractGridFloatArray2D end
abstract type FaceGeometries <: AbstractElementGeometries end
abstract type FaceRegions <: AbstractElementRegions end
abstract type UniqueFaceGeometries <: AbstractElementGeometries end

abstract type BFaceFaces <: AbstractGridIntegerArray1D end
abstract type BFaceCellPos <: AbstractGridIntegerArray1D end # position of bface in adjacent cell
abstract type BFaceVolumes <: AbstractGridFloatArray1D end
abstract type UniqueBFaceGeometries <: AbstractElementGeometries end

abstract type EdgeNodes <: AbstractGridAdjacency end
abstract type EdgeVolumes <: AbstractGridFloatArray1D end
abstract type EdgeCells <: AbstractGridAdjacency end
abstract type EdgeTangents <: AbstractGridFloatArray2D end
abstract type EdgeRegions <: AbstractElementRegions end
abstract type EdgeGeometries <: AbstractElementGeometries end
abstract type UniqueEdgeGeometries <: AbstractElementGeometries end

abstract type BEdgeNodes <: AbstractGridAdjacency end
abstract type BEdgeEdges <: AbstractGridIntegerArray1D end
#abstract type BEdgeCellPos <: AbstractGridIntegerArray1D end # position of bface in adjacent cell
abstract type BEdgeVolumes <: AbstractGridFloatArray1D end
abstract type BEdgeRegions <: AbstractElementRegions end
abstract type BEdgeGeometries <: AbstractElementGeometries end
abstract type UniqueBEdgeGeometries <: AbstractElementGeometries end

abstract type NodePatchGroups <: AbstractGridIntegerArray1D end

## grid item types to dispatch certain things to the correct GridComponents
abstract type ITEMTYPE_NODE end
abstract type ITEMTYPE_CELL end
abstract type ITEMTYPE_FACE end
abstract type ITEMTYPE_BFACE end
abstract type ITEMTYPE_EDGE end
abstract type ITEMTYPE_BEDGE end

abstract type PROPERTY_NODES end
abstract type PROPERTY_VOLUME end
abstract type PROPERTY_REGION end
abstract type PROPERTY_GEOMETRY end
abstract type PROPERTY_UNIQUEGEOMETRY end

GridComponent4TypeProperty(::Type{ITEMTYPE_CELL},::Type{PROPERTY_NODES}) = CellNodes
GridComponent4TypeProperty(::Type{ITEMTYPE_CELL},::Type{PROPERTY_VOLUME}) = CellVolumes
GridComponent4TypeProperty(::Type{ITEMTYPE_CELL},::Type{PROPERTY_REGION}) = CellRegions
GridComponent4TypeProperty(::Type{ITEMTYPE_CELL},::Type{PROPERTY_GEOMETRY}) = CellGeometries
GridComponent4TypeProperty(::Type{ITEMTYPE_CELL},::Type{PROPERTY_UNIQUEGEOMETRY}) = UniqueCellGeometries

GridComponent4TypeProperty(::Type{ITEMTYPE_FACE},::Type{PROPERTY_NODES}) = FaceNodes
GridComponent4TypeProperty(::Type{ITEMTYPE_FACE},::Type{PROPERTY_VOLUME}) = FaceVolumes
GridComponent4TypeProperty(::Type{ITEMTYPE_FACE},::Type{PROPERTY_REGION}) = FaceRegions
GridComponent4TypeProperty(::Type{ITEMTYPE_FACE},::Type{PROPERTY_GEOMETRY}) = FaceGeometries
GridComponent4TypeProperty(::Type{ITEMTYPE_FACE},::Type{PROPERTY_UNIQUEGEOMETRY}) = UniqueFaceGeometries

GridComponent4TypeProperty(::Type{ITEMTYPE_BFACE},::Type{PROPERTY_NODES}) = BFaceNodes
GridComponent4TypeProperty(::Type{ITEMTYPE_BFACE},::Type{PROPERTY_VOLUME}) = BFaceVolumes
GridComponent4TypeProperty(::Type{ITEMTYPE_BFACE},::Type{PROPERTY_REGION}) = BFaceRegions
GridComponent4TypeProperty(::Type{ITEMTYPE_BFACE},::Type{PROPERTY_GEOMETRY}) = BFaceGeometries
GridComponent4TypeProperty(::Type{ITEMTYPE_BFACE},::Type{PROPERTY_UNIQUEGEOMETRY}) = UniqueBFaceGeometries

GridComponent4TypeProperty(::Type{ITEMTYPE_EDGE},::Type{PROPERTY_NODES}) = EdgeNodes
GridComponent4TypeProperty(::Type{ITEMTYPE_EDGE},::Type{PROPERTY_VOLUME}) = EdgeVolumes
GridComponent4TypeProperty(::Type{ITEMTYPE_EDGE},::Type{PROPERTY_REGION}) = EdgeRegions
GridComponent4TypeProperty(::Type{ITEMTYPE_EDGE},::Type{PROPERTY_GEOMETRY}) = EdgeGeometries
GridComponent4TypeProperty(::Type{ITEMTYPE_EDGE},::Type{PROPERTY_UNIQUEGEOMETRY}) = UniqueEdgeGeometries

GridComponent4TypeProperty(::Type{ITEMTYPE_BEDGE},::Type{PROPERTY_NODES}) = BEdgeNodes
GridComponent4TypeProperty(::Type{ITEMTYPE_BEDGE},::Type{PROPERTY_VOLUME}) = BEdgeVolumes
GridComponent4TypeProperty(::Type{ITEMTYPE_BEDGE},::Type{PROPERTY_REGION}) = BEdgeRegions
GridComponent4TypeProperty(::Type{ITEMTYPE_BEDGE},::Type{PROPERTY_GEOMETRY}) = BEdgeGeometries
GridComponent4TypeProperty(::Type{ITEMTYPE_BEDGE},::Type{PROPERTY_UNIQUEGEOMETRY}) = UniqueBEdgeGeometries

function get_facegrid(source_grid::ExtendableGrid{Tc,Ti}) where {Tc,Ti}
    facegrid=ExtendableGrid{Tc,Ti}()
    facegrid[Coordinates]=source_grid[Coordinates]
    facegrid[CellNodes]=source_grid[FaceNodes]
    facegrid[CoordinateSystem]=source_grid[CoordinateSystem]
    facegrid[CellGeometries]=source_grid[FaceGeometries]
    facegrid[UniqueCellGeometries]=source_grid[UniqueFaceGeometries]
    facegrid[CellRegions] = source_grid[FaceRegions]
    # todo: facegrid[CellFaces] = source_grid[FaceEdges]
    return facegrid
end

function get_bfacegrid(source_grid::ExtendableGrid{Tc,Ti})  where {Tc,Ti}
    bfacegrid=ExtendableGrid{Tc,Ti}()
    bfacegrid[Coordinates]=source_grid[Coordinates]
    bfacegrid[CellNodes]=source_grid[BFaceNodes]
    bfacegrid[CoordinateSystem]=source_grid[CoordinateSystem]
    bfacegrid[CellGeometries]=source_grid[BFaceGeometries]
    bfacegrid[UniqueCellGeometries]=source_grid[UniqueBFaceGeometries]
    bfacegrid[CellRegions] = source_grid[BFaceRegions]
    # todo: bfacegrid[CellFaces] = source_grid[BFaceEdges] (or sub-view of FaceEdges ?)
    return bfacegrid
end

function get_edgegrid(source_grid::ExtendableGrid{Tc,Ti}) where {Tc,Ti}
    edgegrid=ExtendableGrid{Tc,Ti}()
    edgegrid[Coordinates]=source_grid[Coordinates]
    edgegrid[CellNodes]=source_grid[EdgeNodes]
    edgegrid[CoordinateSystem]=source_grid[CoordinateSystem]
    edgegrid[CellGeometries]=source_grid[EdgeGeometries]
    edgegrid[UniqueCellGeometries]=source_grid[UniqueedgeGeometries]
    edgegrid[CellRegions] = source_grid[EdgeRegions]
    return edgegrid
end

# show function for ExtendableGrids and defined Components in its Dict
function showmore(io::IO, xgrid::ExtendableGrid{Tc,Ti}) where {Tc,Ti}

    dim = size(xgrid[Coordinates],1)
    nnodes = num_sources(xgrid[Coordinates])
    ncells = num_sources(xgrid[CellNodes])

	println("ExtendableGrid information");
    println("==========================");
	println("dim: $(dim)")
	println("nnodes: $(nnodes)")
    println("ncells: $(ncells)")
    if haskey(xgrid.components,FaceNodes)
        nfaces = num_sources(xgrid[FaceNodes])
        println("nfaces: $(nfaces)")
    else
        println("nfaces: (FaceNodes not instantiated)")
    end
    if haskey(xgrid.components,EdgeNodes)
        nedges = num_sources(xgrid[EdgeNodes])
        println("nedges: $(nedges)")
    else
        println("nedges: (EdgeNodes not instantiated)")
    end
    println("")
    println("Components");
    println("==========");
    for tuple in xgrid.components
        println("> $(tuple[1])")
    end
end

# FaceNodes = nodes for each face (implicitly defines the enumerations of faces)
function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{FaceNodes}) where {Tc,Ti}

    xCellNodes = xgrid[CellNodes]
    ncells = num_sources(xCellNodes)
    nnodes = num_sources(xgrid[Coordinates])
    xCellGeometries = xgrid[CellGeometries]

    # transpose CellNodes to get NodeCells
    xNodeCells = atranspose(xCellNodes)
    max_ncell4node::Ti = max_num_targets_per_source(xNodeCells)

    # instantiate new empty adjacency fields
    xFaceCells = zeros(Ti,0) # cells are appended and at the end rewritten into 2,nfaces array

    # find unique face enumeration rules
    EG = unique(xCellGeometries)
    dim::Ti = dim_element(EG[1])
    face_rules = Array{Array{Ti,2},1}(undef,length(EG))
    maxfacenodes::Ti = 0
    FEG = []
    for j = 1 : length(EG)
        face_rules[j] = local_cellfacenodes(EG[j])
        maxfacenodes = max(size(face_rules[j],2),maxfacenodes)
        for k = 1 : num_faces(EG[j])
            append!(FEG,[facetype_of_cellface(EG[j],j)])
        end
    end
    unique!(FEG)

    # check if only one type of cell geometry is present in grid
    if length(EG) == 1
        singleEG = true
    else
        singleEG = false
    end
    # check if only one type of face geometry is present in grid
    if length(FEG) == 1
        singleFEG = true
    else
        singleFEG = false
    end

    if singleFEG
        xFaceNodes = zeros(Ti,0)
    else
        xFaceNodes = VariableTargetAdjacency(Ti)
        xFaceGeometries::Array{DataType,1} = []
    end
    if singleEG == true && singleFEG == true
        # only one geometry type allows for much faster code
        xCellFaces = zeros(Ti,num_faces(EG[1]),ncells)
        xCellFaceSigns = zeros(Ti,num_faces(EG[1]),ncells)
    else
        xCellFaces = VariableTargetAdjacency(Ti)
        xCellFaceSigns = VariableTargetAdjacency(Ti)
        # pre-allocate xCellFaces
        cellEG = xCellGeometries[1]
        for cell = 1 : ncells
            cellEG = xCellGeometries[cell]
            append!(xCellFaces,zeros(Ti,num_faces(cellEG)))
            append!(xCellFaceSigns,zeros(Ti,num_faces(cellEG)))
        end
    end

    # temporary variables
    # pre-initialised ready to work for singleEG and singleFEG
    node::Ti = 0
    face::Ti = 0
    cell::Ti = 0
    cell2::Ti = 0
    cellEG = EG[1]
    cell2EG = EG[1]
    faceEG = FEG[1]
    faceEG2 = FEG[1]
    face_rule::Array{Ti,2} = face_rules[1]
    face_rule2::Array{Ti,2} = face_rules[1]
    iEG::Ti = 1
    nneighbours::Ti = 0
    faces_per_cell::Ti = num_faces(cellEG)
    faces_per_cell2::Ti = num_faces(cellEG)
    nodes_per_cellface::Ti = num_nodes(faceEG)
    current_item::Array{Ti,1} = zeros(Ti,maxfacenodes) # should be large enough to store largest nnodes per cellface
    flag4item::Array{Bool,1} = zeros(Bool,nnodes)
    no_neighbours_found::Bool = true
    same_face::Bool = false

    # loop over cells
    for cell = 1 : ncells
        # find EG index for geometry
        if !singleEG
            cellEG = xCellGeometries[cell]
            faces_per_cell = num_faces(cellEG)
            for j=1:length(EG)
                if cellEG == EG[j]
                    iEG = j
                    break;
                end
            end
            face_rule = face_rules[iEG]
        end

        # loop over cell faces
        for k = 1 : faces_per_cell

            # check if face is already known to cell
            if xCellFaces[k,cell] > 0
                continue;
            end

            # get face geometry
            if !singleFEG
                faceEG = facetype_of_cellface(cellEG, k)
                nodes_per_cellface = num_nodes(faceEG)
            end

            # flag face nodes and commons4cells
            for j = nodes_per_cellface:-1:1
                node = xCellNodes[face_rule[j,k],cell]
                current_item[j] = node
                flag4item[node] = true;
            end

            # get neighbours for first node
            nneighbours = num_targets(xNodeCells,node)

            # loop over neighbours
            no_neighbours_found = true
            for n = 1 : nneighbours
                cell2 = xNodeCells[n,node]

                # skip if cell2 is the same as cell
                if (cell == cell2)
                    continue;
                end

                # find face enumeration rule
                if !singleEG
                    cell2EG = xCellGeometries[cell2]
                    faces_per_cell2 = num_faces(cell2EG)
                    for j=1:length(EG)
                        if cell2EG == EG[j]
                            iEG = j
                            break;
                        end
                    end
                    face_rule2 = face_rules[iEG]
                end

                # loop over faces face2 of adjacent cell2
                for f2 = 1 : faces_per_cell2

                    # check if face f2 is already known to cell2
                    if xCellFaces[f2,cell2] != 0
                        continue;
                    end

                    # check if face f2 has same geometry
                    if !singleFEG
                        faceEG2 = facetype_of_cellface(cell2EG, f2)
                        if faceEG != faceEG2
                            continue;
                        end
                    end

                    # otherwise compare nodes of face and face2
                    same_face = true
                    for j = 1 : nodes_per_cellface
                        if flag4item[xCellNodes[face_rule2[j,f2],cell2]] == false
                            same_face = false
                            break;
                        end
                    end

                    # if all nodes are the same, register face
                    if (same_face)
                        no_neighbours_found = false
                        face += 1
                        push!(xFaceCells,cell)
                        push!(xFaceCells,cell2)
                        if singleEG == false
                            xCellFaces.colentries[xCellFaces.colstart[cell]+k-1] = face
                            xCellFaces.colentries[xCellFaces.colstart[cell2]+f2-1] = face
                            xCellFaceSigns.colentries[xCellFaceSigns.colstart[cell]+k-1] = 1
                            xCellFaceSigns.colentries[xCellFaceSigns.colstart[cell2]+f2-1] = -1
                            if singleFEG == false
                                push!(xFaceGeometries,faceEG)
                            end
                        else
                            xCellFaces[k,cell] = face
                            xCellFaces[f2,cell2] = face
                            xCellFaceSigns[k,cell] = 1
                            xCellFaceSigns[f2,cell2] = -1
                        end
                        append!(xFaceNodes,view(current_item,1:nodes_per_cellface))
                        break;
                    end
                end
            end

            # if no common neighbour cell is found, register face (boundary faces)
            if no_neighbours_found == true
                face += 1
                push!(xFaceCells,cell)
                push!(xFaceCells,0)
                if singleEG == false
                    xCellFaces.colentries[xCellFaces.colstart[cell]+k-1] = face
                    xCellFaceSigns.colentries[xCellFaceSigns.colstart[cell]+k-1] = 1
                    if singleFEG == false
                        push!(xFaceGeometries,faceEG)
                    end
                else
                    xCellFaces[k,cell] = face
                    xCellFaceSigns[k,cell] = 1
                end
                append!(xFaceNodes,view(current_item,1:nodes_per_cellface))
            end

            # reset flag4item
            for j = 1:nodes_per_cellface
                flag4item[current_item[j]] = false
            end
        end
    end

    if singleFEG
        xFaceNodes = reshape(xFaceNodes,(nodes_per_cellface,Ti(length(xFaceNodes)/nodes_per_cellface)))
        xgrid[FaceGeometries] = VectorOfConstants(facetype_of_cellface(EG[1], 1), face)
    else
        xgrid[FaceGeometries] = xFaceGeometries
    end
    xgrid[CellFaces] = xCellFaces
    xgrid[CellFaceSigns] = xCellFaceSigns
    xgrid[FaceCells] = reshape(xFaceCells,(2,face))
    xFaceNodes
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{NodePatchGroups}) where {Tc,Ti}
    xCellNodes::Adjacency = xgrid[CellNodes]
    xNodeCells = atranspose(xCellNodes)
    nnodes = size(xgrid[Coordinates],2)
    ncells = num_sources(xCellNodes)
    ncells4node::Ti = 0
    group4node = zeros(Ti,nnodes)
    cgroup::Ti = 0
    cell_in_group = zeros(Bool,ncells)
    take_into = false
    cgroup = 0
    while minimum(group4node) == 0
        cell_in_group .= false
        cgroup += 1
        for node = 1 : nnodes
            if group4node[node] == 0
                ncells4node = num_targets(xNodeCells,node)
                take_into = true
                for c = 1 : ncells4node
                    if cell_in_group[xNodeCells[c,node]] == true
                        take_into = false
                        break
                    end
                end
                if take_into
                    group4node[node] = cgroup
                    for c = 1 : ncells4node
                        cell_in_group[xNodeCells[c,node]] = true
                    end
                end
            end
        end
    end
    group4node
end

# FaceNodes = nodes for each face (implicitly defines the enumerations of faces)
function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{EdgeNodes}) where {Tc,Ti}

    xCellNodes = xgrid[CellNodes]
    ncells::Ti = num_sources(xCellNodes)
    nnodes::Ti = num_sources(xgrid[Coordinates])
    xCellGeometries = xgrid[CellGeometries]
    dim::Ti = dim_element(xCellGeometries[1])

    # in 1D edges = cells
    if dim == 1
        xgrid[EdgeNodes] = xgrid[CellNodes]
        xgrid[EdgeGeometries] = xgrid[CellGeometries]
        xgrid[CellEdges] = reshape(collect(Ti,1:size(xCellNodes,2)), (1,size(xCellNodes,2)))
        xgrid[EdgeCells] = xgrid[CellEdges]
        return xgrid[EdgeNodes]
    end

    # in 2D edges = faces
    if dim == 2
        xgrid[CellEdges] = xgrid[CellFaces]
        xgrid[EdgeCells] = xgrid[FaceCells]
        xgrid[EdgeNodes] = xgrid[FaceNodes]
        xgrid[EdgeGeometries] = xgrid[FaceGeometries]
        return xgrid[EdgeNodes]
#        return  prepare_edges!(xgrid) # keep compatible to VoronoiFVM
    end

    # transpose CellNodes to get NodeCells
    xNodeCells = atranspose(xCellNodes)
    max_ncell4node::Ti = max_num_targets_per_source(xNodeCells)

    # init EdgeCells
    xEdgeCells = VariableTargetAdjacency(Ti)

    # find unique edge enumeration rules for each cell geometry
    EG::Array{DataType,1} = unique(xCellGeometries)
    edge_rules::Array{Array{Ti,2},1} = Array{Array{Ti,2},1}(undef,length(EG))
    maxedgenodes::Ti = 0
    for j = 1 : length(EG)
        edge_rules[j] = local_celledgenodes(EG[j])
        maxedgenodes = max(size(edge_rules[j],1),maxedgenodes)
    end

    # pre-allocate xCellEdges and xCellEdgeSigns
    if length(EG) == 1
        singleEG = true
        xCellEdges = zeros(Ti,num_edges(EG[1]),ncells)
        xCellEdgeSigns = zeros(Ti,num_edges(EG[1]),ncells)
    else
        singleEG = false
        xCellEdges = VariableTargetAdjacency(Ti)
        cellEG = xCellGeometries[1]
        for cell = 1 : ncells
            cellEG = xCellGeometries[cell]
            append!(xCellEdges,zeros(Ti,num_edges(cellEG)))
        end
        xCellEdgeSigns = deepcopy(xCellEdges)
    end

    # init xEdgeNodes
    xEdgeNodes::Array{Ti,1} = zeros(Ti,0)

    edge_rule::Array{Ti,2} = edge_rules[1]
    edge_rule2::Array{Ti,2} = edge_rules[1]
    current_item::Array{Ti,1} = zeros(Ti,2)
    flag4item::Array{Bool,1} = zeros(Bool,nnodes)
    cellEG = EG[1]
    node::Ti = 0
    cell2::Ti = 0
    nneighbours::Ti = 0
    edges_per_cell::Ti = num_edges(EG[1])
    edges_per_cell2::Ti = num_edges(EG[1])
    common_nodes::Ti = 0
    cells_with_common_edge::Array{Ti,1} = zeros(Ti,max_ncell4node)
    pos_in_cells_with_common_edge::Array{Ti,1} = zeros(Ti,max_ncell4node)
    sign_in_cells_with_common_edge::Array{Ti,1} = zeros(Ti,max_ncell4node)
    ncells_with_common_edge::Ti = 0
    edge::Ti = 0
    iEG::Ti = 0

    # loop over cells
    for cell = 1 : ncells

        # find EG index for geometry
        if singleEG == false
            cellEG = xCellGeometries[cell]
            edges_per_cell = num_edges(cellEG)
            for j=1:length(EG)
                if cellEG == EG[j]
                    iEG = j
                    break;
                end
            end
            edge_rule = edge_rules[iEG]
        end

        # loop over cell edges
        for k = 1 : edges_per_cell

            # check if edge is already known to cell
            if xCellEdges[k,cell] > 0
                continue;
            end
            ncells_with_common_edge = 1
            cells_with_common_edge[1] = cell
            pos_in_cells_with_common_edge[1] = k
            sign_in_cells_with_common_edge[1] = 1

            # flag edge nodes and commons4cells
            for j = 2 : - 1 : 1
                node = xCellNodes[edge_rule[j,k],cell]
                current_item[j] = node
                flag4item[node] = true;
            end

            # get first node and its neighbours
            nneighbours = num_targets(xNodeCells,node)

            # loop over neighbours
            for n = 1 : nneighbours
                cell2 = xNodeCells[n,node]

                # skip if cell2 is the same as cell
                if (cell == cell2)
                    continue;
                end

                # loop over edges of cell2
                if singleEG == false
                    cellEG = xCellGeometries[cell2]
                    edges_per_cell2 = num_edges(cellEG)

                    # find edge enumeration rule
                    for j=1:length(EG)
                        if cellEG == EG[j]
                            iEG = j
                            break;
                        end
                    end
                    edge_rule2 = edge_rules[iEG]
                end

                for f2 = 1 : edges_per_cell2
                    # compare nodes of edge and edge2
                    common_nodes = 0
                    for j = 1 : 2
                        if flag4item[xCellNodes[edge_rule2[j,f2],cell2]]
                            common_nodes += 1
                        else
                            continue;
                        end
                    end

                    # if all nodes are the same, register edge
                    if (common_nodes == 2)
                        ncells_with_common_edge += 1
                        cells_with_common_edge[ncells_with_common_edge] = cell2
                        pos_in_cells_with_common_edge[ncells_with_common_edge] = f2
                        if xCellNodes[edge_rule2[1,f2],cell2] == current_item[1]
                            sign_in_cells_with_common_edge[ncells_with_common_edge] = 1
                        else
                            sign_in_cells_with_common_edge[ncells_with_common_edge] = -1
                        end
                    end
                end
            end

            # register edge
            edge += 1
            for c = 1 : ncells_with_common_edge
                xCellEdges[pos_in_cells_with_common_edge[c],cells_with_common_edge[c]] = edge
                xCellEdgeSigns[pos_in_cells_with_common_edge[c],cells_with_common_edge[c]] = sign_in_cells_with_common_edge[c]
            end
            append!(xEdgeCells,view(cells_with_common_edge,1:ncells_with_common_edge))
            append!(xEdgeNodes,current_item)

            #reset flag4item
            for j = 1 : 2
                flag4item[current_item[j]] = false
            end
        end
    end
    xgrid[CellEdges] = xCellEdges
    xgrid[EdgeCells] = xEdgeCells
    xgrid[CellEdgeSigns] = xCellEdgeSigns
    xgrid[EdgeGeometries] = VectorOfConstants(Edge1D,edge)
    reshape(xEdgeNodes,(2,edge))
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{CellFaceOrientations}) where {Tc,Ti}
    xCellFaceSigns = xgrid[CellFaceSigns]
    xCellGeometries = xgrid[CellGeometries]
    xCellFaces = xgrid[CellFaces]
    xFaceNodes = xgrid[FaceNodes]
    xCellNodes = xgrid[CellNodes]
    ncells = num_sources(xCellNodes)

    EG = unique(xCellGeometries)
    face_rules = Array{Array{Ti,2},1}(undef,length(EG))
    maxfacenodes = 0
    for j = 1 : length(EG)
        face_rules[j] = local_cellfacenodes(EG[j])
        maxfacenodes = max(size(face_rules[j],2),maxfacenodes)
    end
    singleEG = false
    if typeof(xCellFaceSigns) <: VariableTargetAdjacency
        xCellFaceOrientations = deepcopy(xCellFaceSigns)
    else
        singleEG = true
        xCellFaceOrientations = zeros(Ti, size(xCellFaceSigns,1), size(xCellFaceSigns,2))
    end

    cellEG = EG[1]
    face_rule = face_rules[1]
    ncellfaces::Ti = 0
    nfacenodes::Ti = 0
    face::Ti = 0
    facenodes = zeros(Ti,maxfacenodes)
    found_configuration::Bool = false
    n::Ti = 0
    iEG::Ti = 0
    for cell = 1 : ncells
        cellEG = xCellGeometries[cell]

        # find EG index for geometry
        for j=1:length(EG)
            if cellEG == EG[j]
                iEG = j
                break;
            end
        end
        face_rule = face_rules[iEG] # determines local enumeration of faces

        # determine orientation
        ncellfaces = num_targets(xCellFaces,cell)
        for j = 1 : ncellfaces
            face = xCellFaces[j,cell]
            nfacenodes = num_targets(xFaceNodes,face)
            if xCellFaceSigns[j,cell] == 1
                if singleEG == false
                    xCellFaceOrientations.colentries[xCellFaceOrientations.colstart[cell]+j-1] = 1
                else
                    xCellFaceOrientations[j, cell] = 1
                end
            else
                for k = 1 : nfacenodes
                    facenodes[nfacenodes + 1 - k] = xFaceNodes[k,face]
                end
                found_configuration = false
                n = 0
                while !found_configuration
                    n += 1
                    if facenodes[n] == xCellNodes[face_rule[1,j],cell]
                        found_configuration = true
                    end
                end
                n = mod(n-1,3) + 1
                if singleEG == false
                    xCellFaceOrientations.colentries[xCellFaceOrientations.colstart[cell]+j-1] = 1+n
                else
                    xCellFaceOrientations[j, cell] = 1+n
                end
            end
        end

    end
    return xCellFaceOrientations
end

# FaceEdges = Edges for each face (in 3D)
function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{FaceEdges}) where {Tc,Ti}
    xFaceNodes = xgrid[FaceNodes]
    xFaceCells = xgrid[FaceCells]
    xEdgeNodes = xgrid[EdgeNodes]
    xCellEdges = xgrid[CellEdges]
    xFaceEdges = VariableTargetAdjacency(Ti)
    xFaceGeometries = xgrid[FaceGeometries]
    nfaces::Ti = num_sources(xFaceNodes)
    nnodes::Ti = size(xgrid[Coordinates],2)

    # find unique edge enumeration rules
    EG = unique(xFaceGeometries)
    edge_rules = Array{Array{Ti,2},1}(undef,length(EG))
    maxedgenodes = 0
    for j = 1 : length(EG)
        edge_rules[j] = local_cellfacenodes(EG[j]) # face edges are faces of the face
        maxedgenodes = max(size(edge_rules[j],2),maxedgenodes)
    end
    edge_rule::Array{Ti,2} = edge_rules[1]
    if length(EG) == 1
        singleEG = true
        xFaceEdges = zeros(Ti,num_faces(EG[1]),nfaces)
    else
        singleEG = false
        faceEG = xCellGeometries[1]
        for face = 1 : nfaces
            faceEG = xCellGeometries[cell]
            append!(xCellEdges,zeros(Ti,num_faces(cellEG)))
        end
    end

    nfacenodes::Ti = num_targets(xFaceNodes,1)
    ncelledges::Ti = 0
    nedgenodes::Ti = 0
    nfaceedges::Ti = num_edges(EG[1])
    node::Ti = 0
    cell::Ti = 0
    edge::Ti = 0
    faceedge::Ti = 0
    edge_is_in_face::Bool = false
    found_pos::Bool = false
    flag4face = zeros(Bool,nnodes)
    flag4edge = zeros(Bool,nnodes)
    faceEG = EG[1]
    iEG::Ti = 1
    pos::Ti = 0
    for face = 1 : nfaces

        if singleEG == false
            faceEG = xFaceGeometries[face]
            nfaceedges = num_edges(faceEG)
            # find EG index for geometry
            for j=1:length(EG)
                if faceEG == EG[j]
                    iEG = j
                    break;
                end
            end
            edge_rule = edge_rules[iEG] # determines local enumeration of face edges
            nfacenodes = num_targets(xFaceNodes,face)
        end

        # mark nodes of face
        for j = 1 : nfacenodes
            node = xFaceNodes[j, face]
            flag4face[node] = true;
        end

        # find edges in first adjacent cell
        cell = xFaceCells[1,face]
        ncelledges = num_targets(xCellEdges,cell)
        faceedge = 0
        for cedge = 1 : ncelledges
            edge = xCellEdges[cedge,cell]
            nedgenodes = num_targets(xEdgeNodes,edge)
            edge_is_in_face = true
            for j = 1 : nedgenodes
                if flag4face[xEdgeNodes[j, edge]] == false
                    edge_is_in_face = false
                    break;
                end
            end

            if edge_is_in_face

                # register edge
                # it is important obey same local ordering as in edge_rule
                # to ensure correct dof handling (e.g. for P2-FEM on boundary)

                # mark nodes of edge
                for j = 1 : nedgenodes
                    node = xEdgeNodes[j, edge]
                    flag4edge[node] = true;
                end

                pos = 0
                found_pos = false
                while found_pos == false
                    pos += 1
                    found_pos = true
                    for k = 1 : nedgenodes
                        if flag4edge[xFaceNodes[edge_rule[k,pos], face]] == false
                            found_pos = false
                            break;
                        end
                    end
                end
                if singleEG
                    xFaceEdges[pos,face] = edge
                else
                    xFaceEdges.colentries[xFaceEdges.colstart[face]+pos-1] = edge
                end

                # reset flag4edge
                for j = 1 : nedgenodes
                    node = xEdgeNodes[j, edge]
                    flag4edge[node] = false;
                end

                faceedge += 1
                if faceedge == nfaceedges
                    break;
                end
            end
        end

        #reset flag4face
        for j = 1 : nfacenodes
            node = xFaceNodes[j, face]
            flag4face[node] = false;
        end

    end

    xFaceEdges
end

# CellFaces = faces for each cell
function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{CellFaces}) where {Tc,Ti}
    ExtendableGrids.instantiate(xgrid, FaceNodes)
    xgrid[CellFaces]
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{FaceGeometries}) where {Tc,Ti}
    ExtendableGrids.instantiate(xgrid, FaceNodes)
    xgrid[FaceGeometries]
end

# CellEdges = edges for each cell
function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{CellEdges}) where {Tc,Ti}
    ExtendableGrids.instantiate(xgrid, EdgeNodes)
    xgrid[CellEdges]
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{EdgeCells}) where {Tc,Ti}
    ExtendableGrids.instantiate(xgrid, EdgeNodes)
    xgrid[EdgeCells]
end

# CellEdgeSigns = edges for each cell
function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{CellEdgeSigns}) where {Tc,Ti}
    ExtendableGrids.instantiate(xgrid, EdgeNodes)
    xgrid[CellEdgeSigns]
end

# CellFaceSigns = orientation signs for each face on each cell
function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{CellFaceSigns}) where {Tc,Ti}
    ExtendableGrids.instantiate(xgrid, FaceNodes)
    xgrid[CellFaceSigns]
end

function collectVolumes4Geometries(T::Type{<:Real}, xgrid::ExtendableGrid{Tc,Ti}, ItemType) where {Tc,Ti}
    # get links to other stuff
    xCoordinates = xgrid[Coordinates]
    xCoordinateSystem = xgrid[CoordinateSystem]
    xItemNodes = xgrid[GridComponent4TypeProperty(ItemType,PROPERTY_NODES) ]
    xGeometries = xgrid[GridComponent4TypeProperty(ItemType,PROPERTY_GEOMETRY) ]
    EG = xgrid[GridComponent4TypeProperty(ItemType,PROPERTY_UNIQUEGEOMETRY) ]
    nitems::Ti = num_sources(xItemNodes)

    # get Volume4ElemType handlers
    handlers = Array{Function,1}(undef, length(EG))
    for j = 1: length(EG)
        handlers[j] = Volume4ElemType(xCoordinates, xItemNodes, EG[j], xCoordinateSystem)
    end

    # init Volumes
    xVolumes::Array{T,1} = zeros(T,nitems)

    # loop over items and call handlers
    iEG::Ti = 1
    itemEG = EG[1]
    for item = 1 : nitems
        if length(EG) > 1
            itemEG = xGeometries[item]
            for j=1:length(EG)
                if itemEG == EG[j]
                    iEG = j
                    break;
                end
            end
        end
        xVolumes[item] = handlers[iEG](item)
    end

    xVolumes
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{CellVolumes}) where {Tc,Ti}
    collectVolumes4Geometries(Tc, xgrid, ITEMTYPE_CELL)
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{FaceVolumes}) where {Tc,Ti}
    collectVolumes4Geometries(Tc, xgrid, ITEMTYPE_FACE)
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{BFaceVolumes}) where {Tc,Ti}
    collectVolumes4Geometries(Tc, xgrid, ITEMTYPE_BFACE)
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{EdgeVolumes}) where {Tc,Ti}
    collectVolumes4Geometries(Tc, xgrid, ITEMTYPE_EDGE)
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{BEdgeVolumes}) where {Tc,Ti}
    collectVolumes4Geometries(Tc, xgrid, ITEMTYPE_BEDGE)
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{BFaceFaces}) where {Tc,Ti}
    # get links to other stuff
    xCoordinates = xgrid[Coordinates]
    xFaceNodes = xgrid[FaceNodes]
    xBFaceNodes = xgrid[BFaceNodes]
    nnodes::Ti = num_sources(xCoordinates)
    nbfaces::Ti = num_sources(xBFaceNodes)

    # init BFaceFaces
    xBFaceFaces::Array{Ti,1} = zeros(Ti,nbfaces)
    #xBFaceGeometries = xgrid[BFaceGeometries]
    #if typeof(xBFaceGeometries) == VectorOfConstants{DataType}
    #    EG = xBFaceGeometries[1]
    #    xBFaceGeometries = Array{DataType,1}(undef,nbfaces)
    #    for j = 1 : nbfaces
    #        xBFaceGeometries[j] = EG
    #    end
    #end

    # transpose FaceNodes to get NodeFaces
    xNodeFaces = atranspose(xFaceNodes)

    flag4item::Array{Bool,1} = zeros(Bool,nnodes)
    nodes_per_bface::Ti = 0
    nodes_per_face::Ti = 0
    common_nodes::Ti = 0
    node::Ti = 0
    nneighbours::Ti = 0
    for bface = 1 : nbfaces
        nodes_per_bface = num_targets(xBFaceNodes,bface)
        for j = 1 : nodes_per_bface
            node = xBFaceNodes[j,bface]
            flag4item[node] = true
        end

        # get faces for last node of bface
        nneighbours = num_targets(xNodeFaces,node)

        # loop over faces and find the one that matches the bface
        for n = 1 : nneighbours
            face = xNodeFaces[n,node]
            nodes_per_face = num_targets(xFaceNodes,face)
            common_nodes = 0
            for k = 1 : nodes_per_face
                if flag4item[xFaceNodes[k,face]] == true
                    common_nodes += 1
                else
                    break
                end
            end
            if common_nodes == nodes_per_face
                xBFaceFaces[bface] = face
                break
            end
        end

        if xBFaceFaces[bface] == 0
            println("WARNING(BFaceFaces): found no matching face for bface $bface with nodes $(xBFaceNodes[:,bface])")
        end

        for j = 1 : nodes_per_bface
            flag4item[xBFaceNodes[j,bface]] = false
        end
    end

    # enforce that BFaceNodes have same ordering as FaceNodes
    newBFaceNodes = deepcopy(xBFaceNodes)
    for bface = 1 : nbfaces
        nodes_per_face = num_targets(xBFaceNodes,bface)
        for j = 1 : nodes_per_bface
            if typeof(xBFaceNodes) <: VariableTargetAdjacency
                newBFaceNodes.colentries[newBFaceNodes.colstart[bface]+j-1] = xFaceNodes[j,xBFaceFaces[bface]]
            else
                newBFaceNodes[j,bface] = xFaceNodes[j,xBFaceFaces[bface]]
            end
        end
    end
    xgrid[BFaceNodes] = newBFaceNodes

   # xgrid[BFaceGeometries] = xBFaceGeometries
    xBFaceFaces
end

function instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{BEdgeNodes}) where {Tc,Ti}

    dim = size(xgrid[Coordinates],1)
    if dim==1
        xgrid[BEdgeEdges]=zeros(Ti,0)
        xgrid[BEdgeNodes]=zeros(Ti,0,0)
        return xgrid[BEdgeNodes]
    end

    if dim==2
        xgrid[BEdgeEdges] = xgrid[BFaceFaces]
        xgrid[BEdgeNodes] = xgrid[BFaceNodes]
        return xgrid[BEdgeNodes]
    end

    xBFaceFaces = xgrid[BFaceFaces]
    xEdgeNodes = xgrid[EdgeNodes]
    xFaceEdges = xgrid[FaceEdges]
    xBFaceGeometries = xgrid[BFaceGeometries]
    nbfaces = length(xBFaceFaces)
    xBEdgeEdges = zeros(Ti,0)

    EG = Triangle2D
    edge::Ti = 0
    face::Ti = 0
    nfaceedges::Ti = 0

    for bface = 1 : nbfaces
        EG = xBFaceGeometries[bface]
        nfaceedges = num_edges(EG)
        face = xBFaceFaces[bface]
        for k = 1 : nfaceedges
            edge = xFaceEdges[k,face]
            if !(edge in xBEdgeEdges)
                push!(xBEdgeEdges, edge)
            end

        end
    end

    nbedges = length(xBEdgeEdges)
    xBEdgeNodes = zeros(Ti,2,nbedges)
    for bedge = 1 : nbedges, k = 1 : 2
        xBEdgeNodes[k,bedge] = xEdgeNodes[k,xBEdgeEdges[bedge]]
    end

    xgrid[BEdgeEdges] = xBEdgeEdges
    xBEdgeNodes
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{BEdgeEdges}) where {Tc,Ti}
    ExtendableGrids.instantiate(xgrid, BEdgeNodes)
    xgrid[BEdgeEdges]
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{FaceCells}) where {Tc,Ti}
    ExtendableGrids.instantiate(xgrid, FaceNodes)
    xgrid[FaceCells]
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{FaceRegions}) where {Tc,Ti}
    # interior faces get region number 0, boundary faces get their boundary region
    xBFaceFaces = xgrid[BFaceFaces]
    xBFaceRegions = xgrid[BFaceRegions]
    xFaceRegions = zeros(Ti,num_sources(xgrid[FaceNodes]))
    for j = 1 : length(xBFaceFaces)
        xFaceRegions[xBFaceFaces[j]] = xBFaceRegions[j]
    end
    xFaceRegions
end

### Not a good idea to have them as vectors of constants
# function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{EdgeRegions}) where {Tc,Ti}
#     return VectorOfConstants(Ti(0),num_sources(xgrid[EdgeNodes]))
# end

# function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{BEdgeRegions}) where {Tc,Ti}
#     return VectorOfConstants(Ti(0),num_sources(xgrid[BEdgeNodes]))
# end

# function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{BEdgeGeometries}) where {Tc,Ti}
#     return VectorOfConstants(Edge1D,num_sources(xgrid[BEdgeNodes]))
# end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{BFaceCellPos}) where {Tc,Ti}
    # gets the local position of a bface in the CellFaces row of its adjacent cell

    # get links to other stuff
    xCellFaces = xgrid[CellFaces]
    xFaceCells = xgrid[FaceCells]
    xBFaceFaces = xgrid[BFaceFaces]
    nbfaces = length(xBFaceFaces)

    # init BFaceFaces
    xBFaceCellPos = zeros(Ti,nbfaces)

    cface = 0
    cell = 0
    nfaces4cell = 0
    for bface = 1 : nbfaces
        cface = xBFaceFaces[bface]
        cell = xFaceCells[1,cface]
        nfaces4cell = num_targets(xCellFaces,cell)
        for face = 1 : nfaces4cell
            if cface == xCellFaces[face,cell]
                xBFaceCellPos[bface] = face
                break
            end
        end
    end

    xBFaceCellPos
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{FaceNormals}) where {Tc,Ti}
    dim::Int = size(xgrid[Coordinates],1)
    xCoordinates::Array{Tc,2} = xgrid[Coordinates]
    xFaceNodes::Adjacency{Ti} = xgrid[FaceNodes]
    nfaces::Int = num_sources(xFaceNodes)
    xFaceGeometries::GridEGTypes = xgrid[FaceGeometries]
    xCoordinateSystem = xgrid[CoordinateSystem]
    xFaceNormals::Array{Tc,2} = zeros(Tc,dim,nfaces)
    normal::Array{Tc,1} = zeros(Tc,dim)
    EG = xFaceGeometries[1]
    for face = 1 : nfaces
        EG = xFaceGeometries[face]
        Normal4ElemType!(normal,xCoordinates,xFaceNodes,face,EG,xCoordinateSystem)
        for k = 1 : dim
            xFaceNormals[k, face] = normal[k]
        end
    end
    xFaceNormals
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{EdgeTangents}) where {Tc,Ti}
    dim::Int = size(xgrid[Coordinates],1)
    xCoordinates::Array{Tc,2} = xgrid[Coordinates]
    xEdgeNodes::Adjacency{Ti} = xgrid[EdgeNodes]
    nedges::Int = num_sources(xEdgeNodes)
    xEdgeGeometries::GridEGTypes = xgrid[EdgeGeometries]
    xCoordinateSystem = xgrid[CoordinateSystem]
    xEdgeTangents::Array{Tc,2} = zeros(Tc,dim,nedges)
    EG = xEdgeGeometries[1]
    tangent::Array{Tc,1} = zeros(Tc,dim)
    for edge = 1 : nedges
        EG = xEdgeGeometries[edge]
        Tangent4ElemType!(tangent,xCoordinates,xEdgeNodes,edge,EG,xCoordinateSystem)
        for k = 1 : dim
            xEdgeTangents[k, edge] = tangent[k]
        end
    end
    xEdgeTangents
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{UniqueCellGeometries}) where {Tc,Ti}
    xUniqueCellGeometries = unique(xgrid[CellGeometries])
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{UniqueFaceGeometries}) where {Tc,Ti}
    xUniqueFaceGeometries = unique(xgrid[FaceGeometries])
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{UniqueBFaceGeometries}) where {Tc,Ti}
    xUniqueBFaceGeometries = unique(xgrid[BFaceGeometries])
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{UniqueEdgeGeometries}) where {Tc,Ti}
    xUniqueEdgeGeometries = unique(xgrid[EdgeGeometries])
end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid{Tc,Ti}, ::Type{UniqueBEdgeGeometries}) where {Tc,Ti}
    xUniqueBEdgeGeometries = unique(xgrid[BEdgeGeometries])
end


text:
module Architectures
# Shout out to https://github.com/CliMA/Oceananigans.jl/blob/master/src/Architectures.jl

export
    @hascuda,
    AbstractArchitecture, CPU, GPU,
    devi, array_type,
    architecture,
    default_architecture,
    synchronize_if_gpu

using CUDA

using KernelAbstractions
using CUDAKernels

"""
    AbstractArchitecture
Abstract supertype for architectures supported by Oceananigans.
"""
abstract type AbstractArchitecture end

"""
    CPU <: AbstractArchitecture
Run on a single-core of a CPU.
"""
struct CPU <: AbstractArchitecture end

"""
    GPU <: AbstractArchitecture
Run on a single NVIDIA CUDA GPU.
"""
struct GPU <: AbstractArchitecture end

"""
    @hascuda expr
A macro to compile and execute `expr` only if CUDA is installed and available. Generally used to
wrap expressions that can only be compiled if `CuArrays` and `CUDAnative` can be loaded.
"""
macro hascuda(expr)
    return has_cuda() ? :($(esc(expr))) : :(nothing)
end

devi(::CPU) = KernelAbstractions.CPU()
devi(::GPU) = CUDAKernels.CUDADevice()

         architecture(::Array)   = CPU()
@hascuda architecture(::CuArray) = GPU()

         array_type(::CPU) = Array
@hascuda array_type(::GPU) = CuArray

default_architecture = has_cuda() ? GPU() : CPU()

synchronize_if_gpu() = has_cuda() ? synchronize() : nothing

end


text:
import YAML

doc = YAML.load_file("tests/tests.yaml")

summary = []

for ex in doc["expressions"]
    result = try
        Meta.parse(ex)
        "OK"
    catch
        "Failed"
    end
    push!(summary, (ex, result))
    println(ex, " : ", result)
end

# for ex in doc["symbols"]
#     result = try
#         parse(ex)
#         "OK"
#     catch
#         "Failed"
#     end
#     push!(summary, (ex, result))
#     println(ex, " : ", result)
# end

# println(summary)

max_issues_repo_name: EconForge/Dolang


text:
include("recipes.jl")
# reduce sizes a bit
PyPlot.rc("font", size = 20) # set default fontsize
PyPlot.rc("axes", labelsize = 22)
PyPlot.rc("legend", fontsize = 18)

using Statistics, StatsBase
# using KernelDensity

function field_comparison_timeplot(f1, f2, l1, l2;
    c = "C0", cmap = "viridis", savename = nothing,
    titlestr = "$l1 vs $l2")

    fig, axs = subplots(2, 3, figsize = (23.39, 16.53)) # A2 page size
    fig.suptitle("temporal analysis | "*titlestr)

    for (i, (v, l)) in enumerate(zip((f1, f2), (l1, l2)))
        isnothing(v) && continue
        co = i == 1 ? c : "C1"
        ls = i == 1 ? "solid" : "dashed"
        gv = globalaverage(v)
        # global timeseries
        ax = axs[1, 1]
        μ, σ = timeagg(mean, gv), timeagg(std, gv)
        plot!(ax, gv, label = l*" μ±σ = $(rdspl(μ,4))±$(rdspl(σ,4))", c = co, ls = ls)
        ax.legend()
        ax.set_ylabel("spatial mean")
        ax.set_xlabel("")
        ax.set_xticklabels([])

        ax = axs[2, 1]
        σv = globalagg(std, v)
        μ, σ = timeagg(mean, σv), timeagg(std, σv)
        plot!(ax, σv, label = l*" μ±σ = $(rdspl(μ,4))±$(rdspl(σ,4))", c = co, ls = ls)
        ax.legend()
        ax.set_ylabel("spatial std")
        i == 1 && ax.set_xlabel("")

        # Yearly variability
        ax = axs[2, 3]
        yearly_variability_plot!(ax, globalaverage(v), i; color = co, ls=ls,
        label = l)
        i == 1 && ax.set_xlabel("")

        # Autocorrelation
        ax = axs[1, 2]
        lags = 0:15
        a = autocor(gv, lags)
        ax.plot(lags, a, label = l, c = co, ls = ls)
        ax.set_ylabel("autocorrelation")
        ax.legend()
        ax.grid("on")
    end

    # cross correlation / mutual information
    ax = axs[2, 2]
    gf1, gf2 = globalaverage.((f1, f2))
    lags = -15:15
    cc = crosscor(gf1, gf2, lags)
    ax.plot(lags, cc, color = "C1", label = "Cor($l1, $l2+lag)", ls = "dashed")
    ax.set_ylabel("crosscorrelation")
    mi, center, σ = mutualinformation(Array(gf1), Array(gf2), lags)
    ax.axhspan(center-3σ, center+3σ, label = "95% null", color = "C1", alpha = 0.25, zorder = -9)
    ax.axhline(center, lw = 1.0, ls = "dashed", color = "C1")
    ax.plot(lags, mi, label = "MI($l1, $l2+lag)", color = c)
    ax.set_xlabel("delay [months]")
    ax.grid("on")
    ax.legend()

    # Time trajectory
    ax = axs[1, 3]
    ax.scatter(gf2, gf1, c = 1:length(t), cmap = cmap, edgecolors = "k")
    ax.plot(gf2, gf1, color = c, alpha = 0.75, lw = 0.5, zorder = -9)
    ax.set_ylabel(l1)
    ax.set_xlabel(l2)
    # correlation textboxes:
    rMI = relativeMI(gf2, gf1)
    ρ = corspearman(gf2, gf1)
    ax.text(
        0.95, 0.95, "time trajectory\nρ=$(rdspl(ρ))\nrMI=$(rdspl(rMI))", va = "top",
        bbox = bbox2, transform = ax.transAxes, ha = "right"
    )

    fig.tight_layout()
    fig.subplots_adjust(hspace = 0.13, top = 0.93)

    cax, cb, = add_colorbar!(axs[1,3], 1:length(t);
    name = "months", location = "top", c = cmap)
    !isnothing(savename) && wsave(savename*"_time.png", fig)
    return fig
end

function yearly_variability_plot!(ax, v, i = 1, t = Array(dims(v, Time));
    doscatter=true, shift=true, ploterror = true, kwargs...)
    ut, fluct = yearly_variability(t, v)
    means = mean.(fluct)
    stds = std.(fluct)
    meanσ = mean(stds)
    max = (1:12)
    shift && (max .*= (i-1)*0.1)
    M = maximum(length(l) for l in fluct)
    if ploterror == true
        ax.errorbar(max, means, yerr = stds; marker = "s", ms = 5,
        label = get(kwargs, :label, "")*" \$\\langle\\sigma_\\mathrm{month}\\rangle\$ = $(rdspl(meanσ))",
        color = get(kwargs, :color, "C$(i-1)"),
        kwargs...)
    elseif ploterror == :fill
        ax.plot(max, means;
        label = get(kwargs, :label, "")*" \$\\langle\\sigma_\\mathrm{month}\\rangle\$ = $(rdspl(meanσ))",
        color = get(kwargs, :color, "C$(i-1)"),
        kwargs...)
        ax.fill_between(max, means .- stds, means .+ stds;
        color = get(kwargs, :color, "C$(i-1)"), alpha = 0.5,
        kwargs...)
    else
        ax.plot(max, means;
        label = get(kwargs, :label, "")*" \$\\langle\\sigma_\\mathrm{month}\\rangle\$ = $(rdspl(meanσ))",
        color = get(kwargs, :color, "C$(i-1)"),
        kwargs...)
    end
    if doscatter
        for (j, m) in enumerate(max)
            f = fluct[j]
            ax.plot(fill(m, length(f)), f, marker = "o", ls = "none",
            ms = 6, color = get(kwargs, :color, "C$(i-1)"),
            alpha = 5/M, mew = 0, #fillstyle = "none",
            )
        end
    end
    ax.legend()
    ax.set_ylabel("yearly variability")
    ax.set_xticks(2:3:12)
    ax.set_xticklabels([monthname(τ)[1:3] for τ in t[2:3:12]])
end

function field_comparison_spaceplot(f1, f2, l1, l2;
    c = "C0", cmap = "viridis", savename = nothing, reverselat = false,
    titlestr = "$l1 vs $l2")

    fig, axs = subplots(2, 3, figsize = (23.39, 16.53)) # A2 page size
    fig.suptitle("spatial analysis | "*titlestr)

    for (i, (v, l)) in enumerate(zip((f1, f2), (l1, l2)))
        co = i == 1 ? c : "C1"
        ls = i == 1 ? "solid" : "dashed"
        tv = timeaverage(v)
        if hasdim(v, Lat)
            lv = lataverage(tv)
            # global timeseries
            ax = axs[1, 1]
            μ, σ = mean(lv), std(lv)
            plot!(ax, lv, label = l*", μ±σ = $(rdspl(μ, 4))±$(rdspl(σ))", c = co, ls = ls)
            ax.legend()
            ax.set_ylabel("time & lat mean")
            i == 1 && ax.set_xlabel("")

            # Pdfs
            ax = axs[2, 1]
            val = lv .- mean(lv)
            kde = KernelDensity.kde(val; npoints = 2048*4)
            N = 50
            fs = range(minimum(val), maximum(val); length = 100)
            ax.hist(val.data, N, density = true, color = co, alpha = 0.5)
            ax.plot(fs, KernelDensity.pdf(kde, fs), color = co, label = l, ls = ls)
            ax.legend()
            ax.set_yticks([])
            ax.set_ylabel("pdf of time & lat mean")
        end

        # latitudes
        ax = axs[1, 2]
        zv = zonalaverage(tv)
        ax.plot(sind.(Array(dims(zv, Lat))), zv, ls = ls, color=co, label = l)
        ax.set_xlabel("")
        ax.set_ylabel("time & lon mean")
        ax.legend()

        ax = axs[i, 3]
        yv = timeaverage(v)
        if spacestructure(yv) == Grid()
            lats = sind.(Array(dims(yv, Lat)))
            i == 2 && reverselat && (lats .*= -1)
            p = ax.pcolormesh(Array(dims(yv, Lon)), lats, Array(yv)', cmap = cmap)
            ax.set_aspect("auto")
            fig.colorbar(p, ax = ax)
        elseif spacestructure(yv) == EqArea()
            lons = [x[1] for x in dims(yv, Coord)]
            lats = [x[2] for x in dims(yv, Coord)]
            i == 2 && reverselat && (lats .*= -1)
            ax.scatter(lons, lats, c = Array(yv), cmap = cmap)
            cax, cb = add_colorbar!(ax, range(minimum(yv), maximum(yv); length = 10))
        end
        μ = globalagg(mean, yv)
        σ = globalagg(std, yv)
        ax.text(-0.1, 1.0, "time mean "*l*", μ±σ = $(rdspl(μ))±$(rdspl(σ))",
        transform = ax.transAxes, va = "bottom")
        ax.set_ylabel("\$\\sin(\\phi)\$")
        i == 2 && ax.set_xlabel("longitude λ [ᵒ]")
    end

    # latitude difference
    f1lat = zonalaverage(timeaverage(f1)) |> Array
    f2lat = zonalaverage(timeaverage(f2)) |> Array
    lat = dims(zonalaverage(f1), Lat) |> Array
    ax = axs[2, 2]
    ax.plot(sind.(lat), f1lat .- f2lat, color = c, label = "\$\\psi\$ = NH - SH")
    ax.plot(sind.(lat), (f1lat .- f2lat) .* cosd.(lat), color = "C1",
            alpha = 0.75, label = "\$\\psi \\cos\\phi\$")
    ax.set_xlabel("\$\\sin(\\phi)\$")
    latavg = lataverage(DimensionalArray(f1lat .- f2lat, (Lat(lat),)))
    hemavg = sum(globalaverage.(timeaverage.((f1, f2))))/2
    ax.axhline(latavg, ls = "dashed", color = "C1",
    label = "\$\\langle\\psi\\rangle_\\phi\$ = $(rdspl(latavg)) ($(rdspl(100*abs(latavg)/hemavg))%)")
    ax.legend()
    ax.set_ylabel("time & lon mean")

    fig.tight_layout()
    fig.subplots_adjust(hspace = 0.22, top = 0.91, bottom = 0.08)
    !isnothing(savename) && wsave(savename*"_space.png", fig)
    return fig
end

using FFTW
function surrogate_plot(args...; normalize = false)
    fig = figure()
    ax1 = subplot(2, 1, 1)
    ax1.set_ylabel("\$x\$")
    ax2 = subplot(2, 3, 4)
    ax2.set_ylabel("AC(\$x\$)")
    ax2.set_xlabel("lag")
    ax3 = subplot(2, 3, 5)
    ax3.set_ylabel("\$|\\mathcal{F}(x) |^2\$")
    ax3.set_xlabel("period (in points)")
    ax4 = subplot(2, 3, 6)
    ax4.set_ylabel("pdf(\$x\$)")
    ax4.set_xlabel("\$x\$")
    for (i, x) in enumerate(args)
        normalize && (x = (x .- mean(x))./std(x))
        c = "C$(i-1)"
        ax1.plot(x, color = c, alpha = 0.75, ls = LINES[i])
        ax1.axhline(mean(x), color = c, alpha = 0.75, ls = LINES[i])

        a = autocor(x)
        ax2.plot(0:length(a)-1, a, c=c, alpha = 0.75, ls = LINES[i])

        P = abs.(rfft(x .- mean(x)))
        f = rfftfreq(length(x))
        ps = 1 ./ Float64.(f)
        r = 10:(length(ps))
        ax3.plot(ps[r], P[r] ./ maximum(P[r]), marker = "o",
        color = c, zorder = -99, ms = 4, ls = LINES[i], alpha = 0.75)
        ax3.set_yticks([])

        # kde = KernelDensity.kde(x; npoints = 2048*2)
        fs = range(minimum(x), maximum(x); length = 100)
        ax4.hist(x, 50, density = true, color = c, alpha = 0.5)
        # ax4.plot(fs, KernelDensity.pdf(kde, fs), color = c, alpha = 0.75, ls = LINES[i])
        ax4.set_yticks([])
    end
    fig.tight_layout()
    return fig, ax1
end

using Statistics, StatsBase
function diff_corr_plot(resnhA, resshA, t = 1:length(resnhA))
    lags, lagsa = -15:15, 0:15
    fig = figure()
    ax1 = subplot(2,1,1)
    ax1.plot(t, resnhA .- resshA, label = "difference")
    ax1.axhline(0, label="0", color="C1")
    ax1.set_ylabel("difference"); ax1.legend()

    ax2 = subplot(2,2,3)
    ax2.plot(lagsa, autocor(resnhA, lagsa), label = "AC, NH")
    ax2.plot(lagsa, autocor(resshA, lagsa), label = "AC, SH", ls = "dashed")
    ax2.set_ylabel("correlation"); ax2.legend()
    ax3 = subplot(2,2,4)
    ax3.plot(lags, crosscor(resnhA, resshA, lags), label = "CC, SH")
    mi, center, σ = mutualinformation(Array(resnhA), Array(resshA), lags)
    ax3.axhspan(center-2σ, center+2σ, label = "95% null", color = "C1", alpha = 0.25, zorder = -9)
    ax3.axhline(center, lw = 1.0, ls = "dashed", color = "C1")
    ax3.plot(lags, mi, label = "MI")
    ax3.grid("on")
    ax3.legend()
end

function corr_plot(resnhA, resshA, t = 1:length(resnhA); lx = "x", ly = "y")
    lags, lagsa = -15:15, 0:15
    fig = figure()
    ax1 = subplot(2,1,1)
    x = Array((resnhA .- mean(resnhA)) ./ std(resnhA))
    y = Array((resshA .- mean(resshA)) ./ std(resshA))
    ax1.plot(t, x, label = lx)
    ax1.plot(t, y, label = ly, ls = "--")
    ax1.set_ylabel("standarized signal"); ax1.legend()

    ax2 = subplot(2,2,3)
    ax2.plot(lagsa, autocor(x, lagsa), label = "AC, $(lx)")
    ax2.plot(lagsa, autocor(y, lagsa), label = "AC, $(ly)", ls = "dashed")
    ax2.set_ylabel("correlation"); ax2.legend()
    ax3 = subplot(2,2,4)
    ax3.plot(lags, crosscor(x, y, lags), label = "CC")
    mi, center, σ = mutualinformation(Array(x), Array(y), lags)
    ax3.axhspan(center-2σ, center+2σ, label = "95% null", color = "C1", alpha = 0.25, zorder = -9)
    ax3.axhline(center, lw = 1.0, ls = "dashed", color = "C1")
    ax3.plot(lags, mi, label = "MI")
    ax3.grid("on")
    ax3.legend()
end

max_stars_repo_name: Datseris/EarthAlbedoSymmetry

max_issues_repo_name: Datseris/EarthAlbedoSymmetry

max_forks_repo_name: Datseris/EarthAlbedoSymmetry


text:
using LinearAlgebra, ProgressMeter
import KitBase

function KitBase.write_vtk(ps, ctr)
    cdata = zeros(length(ctr))
    for i in eachindex(cdata)
        cdata[i] = ctr[i].w
    end
    KitBase.write_vtk(ps.points, ps.cellid, cdata)
end

begin
    cd(@__DIR__)

    # time
    tspan = (0.0, 0.2)
    cfl = 0.7

    # quadrature
    quadratureorder = 6
    points, triangulation = KitBase.octa_quadrature(quadratureorder)
    weights = KitBase.quadrature_weights(points, triangulation)
    nq = size(points, 1)
    vs = KitBase.UnstructVSpace(-1.0, 1.0, nq, points, weights)

    # geometry
    cells, points = KitBase.read_mesh("../assets/mesh/pointsource.su2")
    cellid = KitBase.extract_cell(cells)
    edgePoints, edgeCells, cellNeighbors = KitBase.mesh_connectivity_2D(cellid)
    cellType = KitBase.mesh_cell_type(cellNeighbors)
    cellArea = KitBase.mesh_area_2D(points, cellid)
    cellCenter = KitBase.mesh_center_2D(points, cellid)
    edgeCenter = KitBase.mesh_edge_center(points, edgePoints)
    cellEdges = KitBase.mesh_cell_edge(cellid, edgeCells)
    ps = KitBase.UnstructPSpace(
        cells,
        points,
        cellid,
        cellType,
        cellNeighbors,
        cellEdges,
        cellCenter,
        cellArea,
        edgePoints,
        edgeCells,
        edgeCenter,
    )

    # particle
    SigmaS = ones(size(cellid, 1))
    SigmaA = zeros(size(cellid, 1))
    SigmaT = SigmaS + SigmaA
end

ctr = Array{KitBase.ControlVolumeUS1F}(undef, size(ps.cellid, 1))
for i in eachindex(ctr)
    n = Vector{Float64}[]
    for j = 1:3
        push!(
            n,
            KitBase.unit_normal(
                ps.points[edgePoints[cellEdges[i, j], 1], :],
                ps.points[edgePoints[cellEdges[i, j], 2], :],
            ),
        )

        if dot(ps.edgeCenter[ps.cellEdges[i, j], :] .- ps.cellCenter[i, :], n[j]) < 0
            n[j] .= -n[j]
        end
    end

    phi = ones(nq) .* 1e-4
    if 0.49 <= ps.cellCenter[i, 1] <= 0.51 && 0.49 <= ps.cellCenter[i, 2] <= 0.51
        phi .= 10.0
    end
    #s2 = 0.03^2
    #phi .= max(1e-4, 1.0 / (4.0 * π * s2) * exp(-((ps.cellCenter[i, 1]-0.5)^2 + (ps.cellCenter[i, 2]-0.5)^2) / 4.0 / s2))

    w = sum(weights .* phi)
    dx = [
        KitBase.point_distance(
            cellCenter[i, :],
            ps.points[ps.cellid[i, 1], :],
            ps.points[ps.cellid[i, 2], :],
        ),
        KitBase.point_distance(
            cellCenter[i, :],
            ps.points[ps.cellid[i, 2], :],
            ps.points[ps.cellid[i, 3], :],
        ),
        KitBase.point_distance(
            cellCenter[i, :],
            ps.points[ps.cellid[i, 3], :],
            ps.points[ps.cellid[i, 1], :],
        ),
    ]

    ctr[i] = KitBase.ControlVolumeUS1F(n, cellCenter[i, :], dx, w, w, phi)
end

face = Array{KitBase.Interface2D1F}(undef, size(ps.edgePoints, 1))
for i in eachindex(face)
    len = norm(ps.points[edgePoints[i, 1], :] .- ps.points[edgePoints[i, 2], :])
    n = KitBase.unit_normal(ps.points[edgePoints[i, 1], :], ps.points[edgePoints[i, 2], :])

    if !(-1 in ps.edgeCells[i, :])
        n0 = ps.cellCenter[ps.edgeCells[i, 2], :] .- ps.cellCenter[ps.edgeCells[i, 1], :]
    else
        n0 = zero(n)
    end
    if dot(n, n0) < 0
        n .= -n
    end

    fw = 0.0
    ff = zeros(nq)

    face[i] = KitBase.Interface2D1F(len, n[1], n[2], fw, ff)
end

dt = 1.0 / 200 * cfl
nt = tspan[2] ÷ dt |> Int
@showprogress for iter = 1:20#nt
    @inbounds Threads.@threads for i in eachindex(face)
        velo = vs.u[:, 1] .* face[i].n[1] + vs.u[:, 2] .* face[i].n[2]
        if !(-1 in ps.edgeCells[i, :])
            KitBase.flux_kfvs!(
                face[i].ff,
                ctr[ps.edgeCells[i, 1]].f,
                ctr[ps.edgeCells[i, 2]].f,
                velo,
                dt,
            )
        end
    end

    @inbounds Threads.@threads for i in eachindex(ctr)
        if ps.cellType[i] == 0
            for j = 1:3
                dirc = sign(dot(ctr[i].n[j], face[ps.cellEdges[i, j]].n))
                @. ctr[i].f -=
                    dirc * face[ps.cellEdges[i, j]].ff * face[ps.cellEdges[i, j]].len /
                    ps.cellArea[i]
            end

            #integral = KitBase.discrete_moments(ctr[i].f, vs.weights)
            #integral /= 4.0 * π
            #@. ctr[i].f += (integral - ctr[i].f) * dt

            ctr[i].w = sum(ctr[i].f .* vs.weights)
        end
    end
end

KitBase.write_vtk(ps, ctr)

@inbounds Threads.@threads for i in eachindex(face)
    velo = vs.u[:, 1] .* face[i].n[1] + vs.u[:, 2] .* face[i].n[2]
    if !(-1 in ps.edgeCells[i, :])
        KitBase.flux_kfvs!(
            face[i].ff,
            ctr[ps.edgeCells[i, 1]].f,
            ctr[ps.edgeCells[i, 2]].f,
            velo,
            dt,
        )
    end
end

@inbounds Threads.@threads for i in eachindex(ctr)
    if ps.cellType[i] == 0
        for j = 1:3
            dirc = sign(dot(ctr[i].n[j], face[ps.cellEdges[i, j]].n))
            @. ctr[i].f -=
                dirc * face[ps.cellEdges[i, j]].ff * face[ps.cellEdges[i, j]].len /
                ps.cellArea[i]
        end

        #integral = KitBase.discrete_moments(ctr[i].f, vs.weights)
        #integral /= 4.0 * π
        #@. ctr[i].f += (integral - ctr[i].f) * dt

        ctr[i].w = sum(ctr[i].f .* vs.weights)
    end
end


text:
module FunctionsCollection

    function function1(x)
        return 2 * exp.(-2 * (x .- 1).^2) .- exp.(-(x .- 1).^2)
    end

    function function2(x)
        return x.^(1/2)
    end

    function function3(x)
        return exp.(-x) .* sin.(2 * x)
    end

    function function4(x)
        return log.(log.(x))
    end

    function function5(x)
        return 6 * exp.(-2 * x) + 2 * sin.(x) - cos.(x)
    end

end # FunctionsCollection


text: open(IO ->write(IO, read("file1", String)), "file2", "w")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
__precompile__()

include("plain.jl")


text:
#!/usr/bin/env julia
m = map(x -> map(x -> x[1], split(x, "")), readlines(stdin))
let c = 1, d = Dict(), e = [], f = []
    for (i, v) in enumerate(m)
        for (j, x) in enumerate(m[i])
            if m[i][j] in "<>^v"
                if haskey(d, (i, j))
                    println("$(i - 1),$(j - 1)")
                    exit()
                end
                d[(i, j)] = c
                c += 1
                push!(e, 0)
                if m[i][j] in "<>"
                    push!(f, m[i][j] == '<' ? 0 : 2)
                    m[i][j] = '-'
                else
                    push!(f, m[i][j] == '^' ? 1 : 3)
                    m[i][j] = '|'
                end
            end
        end
    end
    while true
        for k in sort(collect(keys(d)))
            u = d[k]
            delete!(d, k)
            i, j = k
            y, x = i, j
            if f[u] % 4 == 0
                x -= 1
            elseif f[u] % 4 == 1
                y -= 1
            elseif f[u] % 4 == 2
                x += 1
            else
                y += 1
            end
            if haskey(d, (y, x))
                println("$(x - 1),$(y - 1)")
                exit()
            end
            d[(y, x)] = u
            if m[y][x] == '+'
                if e[u] % 3 == 0
                    f[u] -= 1
                elseif e[u] %3 == 2
                    f[u] += 1
                end
                e[u] += 1
            elseif m[y][x] == '\\'
                if f[u] % 4 == 0
                    f[u] = 1
                elseif f[u] % 4 == 1
                    f[u] = 0
                elseif f[u] % 4 == 2
                    f[u] = 3
                else
                    f[u] = 2
                end
            elseif m[y][x] == '/'
                f[u] = 3 - (f[u] % 4)
            end
        end
    end
end


text:
"Index mapping variables to integers"
struct IndexDict
    index::Dict{Sym, Int}
end

IndexDict(vars::Vector{<:Sym}) = IndexDict(Dict(v => i for (i, v) in enumerate(vars)))

struct Model
    variables
    constraints
end

Model() = Model([], [])

push!(m::Model, var) = push!(m.variables, var)
push!(m::Model, c::Constraint) = push!(m.constraints, c)

function clauses(m::Model)
    variable_clauses = reduce(∪, clauses.(m.variables))
    constraint_clauses = reduce(∪, clauses.(m.constraints))

    return variable_clauses ∪ constraint_clauses
end

variables(m::Model) = m.variables
boolean_variables(m::Model) = reduce(vcat, (v.vars for v in variables(m)))

function solve(m::Model)
    p = SymbolicSATProblem(m)

    status, result_dict = solve(p)

    (status == :unsat) && return status, missing

    return status, Dict(v => decode(result_dict, v) for v in m.variables)
end


text:
module test_04_ConstantSignalsWithUnit

using ModiaPlot

include("test_base/test_04_ConstantSignalsWithUnit.jl")

end


text:
@testset "962.maximum-width-ramp.jl" begin
    @test max_width_ramp([41, 19, 31, 50, 7, 13, 43, 3, 4, 32, 16, 11, 38, 31, 27, 32, 29, 21, 2, 45, 42, 34, 17, 1, 16, 45, 48, 33, 19]) == 27
    @test max_width_ramp([6,0,8,2,1,5]) == 4
    @test max_width_ramp([9,8,1,0,1,9,4,0,4,1]) == 7
end


text: # some code here

max_stars_repo_name: tomtrogdon/demo

max_issues_repo_name: tomtrogdon/demo

max_forks_repo_name: tomtrogdon/demo


text:
@require Blink begin

function Blink.body!(p::Blink.Page, x::AbstractWidget)
    Blink.body!(p, WebIO.render(x))
end

function Blink.body!(p::Blink.Window, x::AbstractWidget)
    Blink.body!(p, WebIO.render(x))
end

end