text:
function testlut(
    # The SA structure that is being placed.
    sa_struct,
    # Instantied move_generator. It is assumed that this generator has
    # already been initialized.
    move_generator::SA.CachedMoveGenerator,
    # Class index currently under test.
    idx::Integer,
    # List of addresses of tiles that it is expected that this class
    # maps to.
    addrs,
    # Distance function between addresses. Useful for using this same
    # function for Cartesian architectures and hexagonal architectures.
    dist::Function,
)
    # Get the move dictionary form the Move Generagor for this index .
    move_dict = move_generator.moves[idx]

    # Make sure that the keys of this dictionary match the expected
    # addresses.
    @test sort(collect(keys(move_dict))) == sort(addrs)

    # Test all limits for this geometry.
    for limit in SA.distancelimit(move_generator, sa_struct):-1:1
        # Update the move generator.
        SA.update!(move_generator, sa_struct, limit)

        # Iterate over all valid addresses. We'll recompute what the
        # collection of moves should look like from scratch and verify
        # the match.
        for (address, lut) in move_dict
            # Compute the expected set of addresses from scratch.
            expected_addresses = sort(
                [a for a in addrs if dist(a, address) <= limit];
                # Sort addresses by distances from "address"
                # Then apply the normal CartesianIndex sort.
                lt = (x, y) ->
                    (dist(x, address) < dist(y, address)) ||
                        (dist(x, address) == dist(y, address) && x < y),
            )
            @test expected_addresses == lut.targets[1:(lut.idx)]

            # Test that all of the addresses in the targets are within the
            # distance limit and none of the addresses outside this target
            # are within the distance limit.
            norms_below = map(x -> dist(x, address), lut.targets[1:(lut.idx)])
            @test all(norms_below .<= limit)

            # Test the remainder of the array.
            norms_above = map(x -> dist(x, address), lut.targets[(lut.idx + 1):end])
            @test all(norms_above .> limit)

            # Test random move genertion.
            # Only test (and thus fail) if generation fails.
            for _ in 1:1000
                rand_address = rand(lut)
                if rand_address ∉ expected_addresses
                    @test in(rand_address, expected_addresses)
                end
            end
        end
    end
end

function get_colors(m::Map, sa_struct::SA.SAStruct)
    # Decode the classes.
    color_to_class_idx = Dict{SquareColor,Int}()

    tasktable_rev = Dict(v => k for (k, v) in sa_struct.tasktable)
    for (idx, node) in enumerate(sa_struct.nodes)
        class = SA.getclass(node)
        taskname = tasktable_rev[idx]
        color = getnode(m.taskgraph, taskname).metadata["color"]

        # Initialize and check the colors.
        if !haskey(color_to_class_idx, color)
            color_to_class_idx[color] = class
        else
            @test color_to_class_idx[color] == class
        end
    end

    return color_to_class_idx
end

################################################################################
@testset "Testing Move Generators" begin
    import Base.Iterators.product

    # Use the Chessboard module
    # Build a 4x4 rectangular architecture with chessboard coloring.
    board_dims = 4
    A = architecture(board_dims, Rectangle2D(), ChessboardColor())
    # Construct a taskgraph that will entierly fill the architecture.
    # Color it so 1/4 the tasks are Black, 1/4 are White, and 1/2 are Gray.
    T = taskgraph(board_dims^2, board_dims^2 + board_dims, Quarters())

    m = Map(Chess(), A, T)

    # Build the SA Struct
    sa_struct = SA.SAStruct(m)

    # SA here should have 3 classes for Black, White, and Gray.
    # We need to figure out which class index refers set of processors.
    @test length(sa_struct.maptable.mask) == 3

    color_to_class_idx = get_colors(m, sa_struct)

    # Should have found all of the colors.
    @test length(color_to_class_idx) == 3

    # Chess Board:
    #
    # B W B W
    # W B W B
    # B W B W
    # W B W B

    # Initialize a CachedMoveGenerator.
    move_generator = SA.CachedMoveGenerator{SA.location_type(sa_struct.maptable)}()

    # Test that the maximum distance from corner to corner is correct.
    @test SA.distancelimit(move_generator, sa_struct) == 2 * (board_dims - 1)
    # Initialize the move generator.
    SA.initialize!(move_generator, sa_struct)

    # Convenience Definition
    CI = CartesianIndex

    euclidean(x, y) = sum(abs.((x - y).I))

    #################
    # TESTING WHITE #
    #################
    idx = color_to_class_idx[White]

    addrs = [CI(i, j) for (i, j) in product(1:board_dims, 1:board_dims) if isodd(i + j)]
    testlut(sa_struct, move_generator, idx, addrs, euclidean)

    #################
    # TESTING BLACK #
    #################
    idx = color_to_class_idx[Black]

    addrs = [CI(i, j) for (i, j) in product(1:board_dims, 1:board_dims) if iseven(i + j)]
    testlut(sa_struct, move_generator, idx, addrs, euclidean)

    ################
    # TESTING GRAY #
    ################
    idx = color_to_class_idx[Gray]

    addrs = reshape([CI(i, j) for (i, j) in product(1:board_dims, 1:board_dims)], :)
    testlut(sa_struct, move_generator, idx, addrs, euclidean)

    SA.place!(sa_struct; move_attempts = 20)
    @test SA.verify_placement(m, sa_struct)
end

@testset "Big Move Generator Test" begin
    # In this test, we're going to try the 3D, 2D, hexagonal architectures with
    # the Chessboard and Hash color schemes. In all cases, the
    # CachedMoveGenerator should generate correct results ... hopefully.

    # Distance between two coordinates using rectilinear coordinates.
    euclidean(x, y) = sum(abs.((x - y).I))
    function hexagonal(a, b)
        # Break cartesian indices into their components.
        x0, y0 = a[1], a[2]
        x1, y1 = b[1], b[2]

        # Find the distance between y coordinages
        d = abs(y1 - y0)

        # Compute how many diagonal steps we can take given the y displacement.
        #
        # Because of the awkward offset in the hexagonal architecture, we have
        # to check if we're on an even column or not.
        if isodd(y0)
            x_hi = x0 + (d >> 1)
            x_lo = x0 - ((d + 1) >> 1)
        else
            x_hi = x0 + ((d + 1) >> 1)
            x_lo = x0 - (d >> 1)
        end

        # Check the location of the "x" coordinate relative to the number of
        # diagonal steps we could have taken. If the "x" coordinate is outside
        # of the bound, we have to take some more hops to get there. Otherwise,
        # it "x" is inside that bound, we don't have to take any extra steps.
        if x1 < x_lo
            dist = d + x_lo - x1
        elseif x1 > x_hi
            dist = d + x1 - x_hi
        else
            dist = d
        end

        return dist
    end

    # Build the maps.
    maps = [
        (
            Map(
                Chess(),
                architecture(4, Rectangle2D(), ChessboardColor()),
                taskgraph(16, 20, Quarters()),
            ),
            euclidean,
        ),
        (
            Map(
                Chess(),
                architecture(4, Rectangle2D(), HashColor()),
                taskgraph(16, 20, Quarters()),
            ),
            euclidean,
        ),
        (
            Map(
                Chess(),
                architecture(3, Rectangle3D(), ChessboardColor()),
                taskgraph(27, 40, Quarters()),
            ),
            euclidean,
        ),
        (
            Map(
                Chess(),
                architecture(3, Rectangle3D(), HashColor()),
                taskgraph(27, 40, Quarters()),
            ),
            euclidean,
        ),
        (
            Map(
                Chess(),
                architecture(4, Hexagonal2D(), ChessboardColor()),
                taskgraph(16, 24, Quarters()),
            ),
            hexagonal,
        ),
    ]

    for (m, dist) in maps
        sa_struct = SA.SAStruct(m)

        # Get the color to class index mapping
        color_to_class_idx = get_colors(m, sa_struct)
        # Should have found all of the colors.
        @test length(color_to_class_idx) == 3

        # Initialize a CachedMoveGenerator.
        move_generator = SA.CachedMoveGenerator{SA.location_type(sa_struct.maptable)}()
        # Initialize the move generator.
        SA.initialize!(move_generator, sa_struct)

        # Iterate through all colors
        for color in (White, Black, Gray)
            # Get the class index for this color
            idx = color_to_class_idx[color]
            # Search through all the addresses in the architecture. Get the ones
            # where the mappable component has the color in question.
            if color == Gray
                addrs = collect(addresses(m.toplevel))
            else
                addrs = [
                    a for a in addresses(m.toplevel) if
                    search_metadata!(m.toplevel[a], "color", color, ==)
                ]
            end

            # Test the move generator
            testlut(sa_struct, move_generator, idx, addrs, dist)
        end
    end
end


text:
# wrappers of low-level functionality

function curandGetProperty(property::libraryPropertyType)
  value_ref = Ref{Cint}()
  curandGetProperty(property, value_ref)
  value_ref[]
end

version() = VersionNumber(curandGetProperty(CUDA.MAJOR_VERSION),
                          curandGetProperty(CUDA.MINOR_VERSION),
                          curandGetProperty(CUDA.PATCH_LEVEL))

ext: jl


text:
module BlockBandedMatrices
using BlockArrays, BandedMatrices, ArrayLayouts, FillArrays, SparseArrays, MatrixFactorizations
using LinearAlgebra

import Base: getindex, setindex!, checkbounds, @propagate_inbounds, convert,
                        +, *, -, /, \, strides, zeros, size,
                        unsafe_convert, fill!, length, first, last,
                        eltype, getindex, to_indices, to_index,
                        reindex, _maybetail, tail, @_propagate_inbounds_meta,
                        ==, axes, copyto!, similar, OneTo, replace_in_print_matrix

import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, DefaultArrayStyle, Broadcasted, broadcasted,
                        materialize, materialize!

import LinearAlgebra: UniformScaling, isdiag, rmul!, lmul!, ldiv!, rdiv!,
                        AbstractTriangular, AdjOrTrans, HermOrSym, StructuredMatrixStyle,
                        qr, qr!, QRPackedQ
import LinearAlgebra.BLAS: BlasInt, BlasFloat, @blasfunc, libblas, BlasComplex, BlasReal
import LinearAlgebra.LAPACK: chktrans, chkdiag, liblapack, chklapackerror, checksquare, chkstride1,
                    chkuplo
import MatrixFactorizations: ql, ql!, QLPackedQ
import SparseArrays: sparse

import ArrayLayouts: @lazymul, MatMulMatAdd, MatMulVecAdd, BlasMatLmulVec,
                    triangularlayout, UpperTriangularLayout, TriangularLayout, MatLdivVec,
                    triangulardata, sublayout, @lazyldiv, @lazylmul,
                    AbstractColumnMajor, DenseColumnMajor, ColumnMajor,
                    DiagonalLayout, materialize!, MulAdd, mul

import BlockArrays: BlockSizes, nblocks, blocksize, blockcheckbounds, global2blockindex,
                        Block, BlockSlice, getblock, unblock, setblock!, globalrange,
                        _unblock, _find_block, BlockIndexRange, blocksizes, cumulsizes,
                        AbstractBlockSizes, sizes_from_blocks, DefaultBlockSizes

import BandedMatrices: isbanded, bandwidths, bandwidth, banded_getindex, colrange,
                        inbands_setindex!, inbands_getindex, banded_setindex!,
                        banded_generic_axpy!,
                        BlasFloat, banded_dense_axpy!, MemoryLayout,
                        BandedColumnMajor,
                        BandedSubBandedMatrix, bandeddata, tribandeddata,
                        _BandedMatrix, colstart, colstop, rowstart, rowstop,
                        BandedStyle, _fill_lmul!,
                        _banded_colval, _banded_rowval, _banded_nzval # for sparse

export BandedBlockBandedMatrix, BlockBandedMatrix, BlockSkylineMatrix, blockbandwidth, blockbandwidths,
        subblockbandwidth, subblockbandwidths, Ones, Zeros, Fill, Block, BlockTridiagonal, isblockbanded

include("AbstractBlockBandedMatrix.jl")
include("broadcast.jl")
include("BlockSkylineMatrix.jl")
include("BandedBlockBandedMatrix.jl")

include("linalg.jl")
include("blockskylineqr.jl")

include("interfaceimpl.jl")
include("triblockbanded.jl")
include("adjtransblockbanded.jl")

end # module


text:
function threadrandinit!(sptr::StackPointer, pcg_vector::Vector{PtrPCG{P}}) where {P}
    W = VectorizationBase.pick_vector_width(Float64)
#    myprocid = 0#myid()-1
    local_stack_size = LOCAL_STACK_SIZE[]
    stack_ptr = sptr
    for n ∈ 1:NTHREADS[]
        rng = VectorizedRNG.random_init_pcg!(PtrPCG{P}(stack_ptr))#, P*(n-1)*myprocid)
        if n > length(pcg_vector)
            push!(pcg_vector, rng)
        else
            pcg_vector[n] = rng
        end
        stack_ptr += local_stack_size
    end
    sptr + (2P+1)*8W
end


text:
using Base.Test

using DataFrames

import Celeste: GalsimBenchmark

# Ensure GalSim test images are available.
wd = pwd()
cd(joinpath(Pkg.dir("Celeste"), "benchmark", "galsim"))
run(`make fetch`)
cd(wd)

if test_long_running
    GALSIM_CASES_EXERCISED = [
        "simple_star",
        "star_with_noise",
        "angle_and_axis_ratio_1",
        "galaxy_with_all",
        "galaxy_with_noise",
    ]
else
    GALSIM_CASES_EXERCISED = [
        "star_with_noise",
    ]
end

function assert_estimates_are_close(benchmark_results)
    for row in eachrow(benchmark_results)
        if ismissing(row[:truth])
            continue
        end
        if row[:variable] == :is_star
            maximum_error = 0.1
        elseif row[:variable] == :gal_frac_dev
            maximum_error = 0.2
        elseif row[:variable] == :gal_angle_deg
            maximum_error = 5
        else
            maximum_error = 0.2 * abs(row[:truth])
        end
        @test isapprox(row[:truth], row[:estimate],
                       atol=maximum_error)
    end
end

@testset "galsim benchmarks" begin
    @testset "single-source inference" begin
        truth, results = GalsimBenchmark.run_benchmarks(
            test_case_names=GALSIM_CASES_EXERCISED,
            joint_inference=false
        )
        assert_estimates_are_close(
            GalsimBenchmark.truth_comparison_df(truth, results))
    end

    @testset "GalSim benchmark tests, joint inference" begin
        truth, results = GalsimBenchmark.run_benchmarks(
            test_case_names=GALSIM_CASES_EXERCISED,
            joint_inference=true
        )
        assert_estimates_are_close(
            GalsimBenchmark.truth_comparison_df(truth, results))
    end
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct VersionDef <: ProtoType
    producer::Int32
    min_consumer::Int32
    bad_consumers::Base.Vector{Int32}
    VersionDef(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct VersionDef
const __pack_VersionDef = Symbol[:bad_consumers]
meta(t::Type{VersionDef}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_VersionDef, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

export VersionDef


text:
module Manifolds

import ManifoldsBase:
    _access_nested,
    _read,
    _write,
    allocate,
    allocate_result,
    allocate_result_type,
    allocation_promotion_function,
    array_value,
    base_manifold,
    check_point,
    check_point__transparent,
    check_vector,
    decorated_manifold,
    decorator_transparent_dispatch,
    default_decorator_dispatch,
    distance,
    dual_basis,
    embed,
    embed!,
    exp!,
    exp!__intransparent,
    get_basis,
    get_component,
    get_coordinates,
    get_coordinates!,
    get_embedding,
    get_iterator,
    get_vector,
    get_vector!,
    get_vectors,
    gram_schmidt,
    injectivity_radius,
    inner,
    inner__intransparent,
    is_point,
    is_vector,
    inverse_retract,
    inverse_retract!,
    log,
    log!,
    manifold_dimension,
    mid_point,
    mid_point!,
    number_eltype,
    number_of_coordinates,
    power_dimensions,
    project,
    project!,
    representation_size,
    retract,
    retract!,
    set_component!,
    vector_space_dimension,
    vector_transport_direction,
    vector_transport_direction!,
    vector_transport_to,
    vector_transport_to!,
    zero_vector,
    zero_vector!,
    CotangentSpace,
    TangentSpace
import Base: copyto!, convert, in, isapprox, isempty, length, showerror

using Base.Iterators: repeated
using Distributions
using Einsum: @einsum
using FiniteDifferences
using HybridArrays
using LightGraphs
using LinearAlgebra
using ManifoldsBase
using ManifoldsBase:
    ℝ,
    ℂ,
    ℍ,
    AbstractBasis,
    AbstractIsometricEmbeddingType,
    AbstractNumbers,
    AbstractOrthogonalBasis,
    AbstractOrthonormalBasis,
    AbstractPowerManifold,
    AbstractPowerRepresentation,
    AbstractVectorTransportMethod,
    AbstractLinearVectorTransportMethod,
    ApproximateInverseRetraction,
    ApproximateRetraction,
    DifferentiatedRetractionVectorTransport,
    ComponentManifoldError,
    CompositeManifoldError,
    CotangentSpaceType,
    CoTFVector,
    DefaultManifold,
    DefaultOrDiagonalizingBasis,
    DiagonalizingBasisData,
    DifferentiatedRetractionVectorTransport,
    FVector,
    NLsolveInverseRetraction,
    InversePowerRetraction,
    PowerManifold,
    PowerManifoldNested,
    PowerRetraction,
    PowerVectorTransport,
    TangentSpaceType,
    TCoTSpaceType,
    TFVector,
    VectorSpaceType,
    VeeOrthogonalBasis,
    @decorator_transparent_fallback,
    @decorator_transparent_function,
    @decorator_transparent_signature,
    @invoke_maker,
    _euclidean_basis_vector,
    _extract_val,
    combine_allocation_promotion_functions,
    is_decorator_transparent,
    is_default_decorator,
    manifold_function_not_implemented_message,
    rep_size_to_colons,
    size_to_tuple
using Markdown: @doc_str
using Random
using RecipesBase
using RecipesBase: @recipe, @series
using Colors: RGBA
using Requires
using SimpleWeightedGraphs: AbstractSimpleWeightedGraph, get_weight
using StaticArrays
using Statistics
using StatsBase
using StatsBase: AbstractWeights

include("utils.jl")

include("atlases.jl")
include("cotangent_space.jl")
include("differentiation.jl")
include("riemannian_diff.jl")

include("statistics.jl")

include("product_representations.jl")

include("manifolds/VectorBundle.jl")

include("distributions.jl")
include("projected_distribution.jl")

# It's included early to ensure visibility of `Identity`
include("groups/group.jl")

include("manifolds/MetricManifold.jl")
include("manifolds/ProductManifold.jl")
include("manifolds/PowerManifold.jl")
include("manifolds/GraphManifold.jl")

include("manifolds/Euclidean.jl")
include("manifolds/Lorentz.jl")

include("manifolds/CenteredMatrices.jl")
include("manifolds/CholeskySpace.jl")
include("manifolds/Circle.jl")
include("manifolds/Elliptope.jl")
include("manifolds/FixedRankMatrices.jl")
include("manifolds/GeneralizedGrassmann.jl")
include("manifolds/GeneralizedStiefel.jl")
include("manifolds/Grassmann.jl")
include("manifolds/Hyperbolic.jl")
include("manifolds/MultinomialDoublyStochastic.jl")
include("manifolds/MultinomialSymmetric.jl")
include("manifolds/ProbabilitySimplex.jl")
include("manifolds/PositiveNumbers.jl")
include("manifolds/ProjectiveSpace.jl")
include("manifolds/Rotations.jl")
include("manifolds/SkewHermitian.jl")
include("manifolds/Spectrahedron.jl")
include("manifolds/Stiefel.jl")
include("manifolds/StiefelEuclideanMetric.jl")
include("manifolds/StiefelCanonicalMetric.jl")
include("manifolds/Sphere.jl")
include("manifolds/SphereSymmetricMatrices.jl")
include("manifolds/Symmetric.jl")
include("manifolds/SymmetricPositiveDefinite.jl")
include("manifolds/SymmetricPositiveDefiniteLinearAffine.jl")
include("manifolds/SymmetricPositiveDefiniteLogCholesky.jl")
include("manifolds/SymmetricPositiveDefiniteLogEuclidean.jl")
include("manifolds/SymmetricPositiveSemidefiniteFixedRank.jl")

# Product or power based manifolds
include("manifolds/Torus.jl")
include("manifolds/Multinomial.jl")
include("manifolds/Oblique.jl")
include("manifolds/EssentialManifold.jl")

include("groups/metric.jl")
include("groups/group_action.jl")
include("groups/group_operation_action.jl")
include("groups/array_manifold.jl")
include("groups/product_group.jl")
include("groups/semidirect_product_group.jl")

include("groups/general_linear.jl")
include("groups/special_linear.jl")
include("groups/translation_group.jl")
include("groups/special_orthogonal.jl")
include("groups/circle_group.jl")

include("groups/translation_action.jl")
include("groups/rotation_action.jl")

include("groups/special_euclidean.jl")

@doc raw"""
    Base.in(p, M::AbstractManifold; kwargs...)
    p ∈ M

Check, whether a point `p` is a valid point (i.e. in) a [`AbstractManifold`](@ref) `M`.
This method employs [`is_point`](@ref) deaticating the error throwing option.
"""
Base.in(p, M::AbstractManifold; kwargs...) = is_point(M, p, false; kwargs...)

@doc raw"""
    Base.in(p, TpM::TangentSpaceAtPoint; kwargs...)
    X ∈ TangentSpaceAtPoint(M,p)

Check whether `X` is a tangent vector from (in) the tangent space $T_p\mathcal M$, i.e.
the [`TangentSpaceAtPoint`](@ref) at `p` on the [`AbstractManifold`](@ref) `M`.
This method uses [`is_vector`](@ref) deactivating the error throw option.
"""
function Base.in(X, TpM::TangentSpaceAtPoint; kwargs...)
    return is_vector(base_manifold(TpM), TpM.point, X, false; kwargs...)
end

function __init__()
    @require FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" begin
        using .FiniteDiff
        include("finite_diff.jl")
    end

    @require ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" begin
        using .ForwardDiff
        include("forward_diff.jl")
    end

    @require OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" begin
        using .OrdinaryDiffEq: ODEProblem, AutoVern9, Rodas5, solve
        include("ode.jl")
    end

    @require NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" begin
        using .NLsolve: NLsolve
        include("nlsolve.jl")
    end

    @require Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" begin
        using .Test: Test
        include("tests/tests_general.jl")
        export test_manifold
        include("tests/tests_group.jl")
        export test_group, test_action
        @require ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" begin
            include("tests/tests_forwarddiff.jl")
            export test_forwarddiff
        end

        @require ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" begin
            include("tests/tests_reversediff.jl")
            export test_reversediff
        end
    end

    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
        using RecipesBase: @recipe, @series
        using Colors: RGBA
        include("recipes.jl")
    end

    @require RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" begin
        @require Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" begin
            using .RecipesBase: @recipe, @series
            using Colors: RGBA
            include("recipes.jl")
        end
    end
    return nothing
end

#
export CoTVector, AbstractManifold, AbstractManifoldPoint, TVector, AbstractManifold
export AbstractSphere, AbstractProjectiveSpace
export Euclidean,
    ArrayProjectiveSpace,
    ArraySphere,
    CenteredMatrices,
    CholeskySpace,
    Circle,
    Elliptope,
    EssentialManifold,
    FixedRankMatrices,
    GeneralizedGrassmann,
    GeneralizedStiefel,
    Grassmann,
    Hyperbolic,
    Lorentz,
    MultinomialDoubleStochastic,
    MultinomialMatrices,
    MultinomialSymmetric,
    Oblique,
    PositiveArrays,
    PositiveMatrices,
    PositiveNumbers,
    PositiveVectors,
    ProbabilitySimplex,
    ProjectiveSpace,
    Rotations,
    SkewHermitianMatrices,
    SkewSymmetricMatrices,
    Spectrahedron,
    Sphere,
    SphereSymmetricMatrices,
    Stiefel,
    SymmetricMatrices,
    SymmetricPositiveDefinite,
    SymmetricPositiveSemidefiniteFixedRank,
    Torus
export HyperboloidPoint, PoincareBallPoint, PoincareHalfSpacePoint, SVDMPoint
export HyperboloidTVector, PoincareBallTVector, PoincareHalfSpaceTVector, UMVTVector
export AbstractNumbers, ℝ, ℂ, ℍ

# decorator manifolds
export AbstractDecoratorManifold
export ValidationManifold, ValidationMPoint, ValidationTVector, ValidationCoTVector
export CotangentBundle,
    CotangentSpaceAtPoint, CotangentBundleFibers, CotangentSpace, FVector
export AbstractPowerManifold,
    AbstractPowerRepresentation,
    ArrayPowerRepresentation,
    NestedPowerRepresentation,
    NestedReplacingPowerRepresentation,
    PowerManifold
export ProductManifold
export GraphManifold, GraphManifoldType, VertexManifold, EdgeManifold
export ProjectedPointDistribution, ProductRepr, TangentBundle, TangentBundleFibers
export TangentSpace, TangentSpaceAtPoint, VectorSpaceAtPoint, VectorSpaceType, VectorBundle
export VectorBundleFibers
export AbstractVectorTransportMethod,
    DifferentiatedRetractionVectorTransport, ParallelTransport, ProjectedPointDistribution
export PoleLadderTransport, SchildsLadderTransport
export PowerVectorTransport, ProductVectorTransport
export AbstractEmbeddedManifold
export AbstractMetric,
    RiemannianMetric,
    LorentzMetric,
    EmbeddedManifold,
    EuclideanMetric,
    LinearAffineMetric,
    LogCholeskyMetric,
    LogEuclideanMetric,
    MinkowskiMetric,
    PowerMetric,
    ProductMetric,
    CanonicalMetric,
    MetricManifold
export AbstractAtlas, RetractionAtlas
export AbstractEmbeddingType, AbstractIsometricEmbeddingType
export DefaultEmbeddingType, DefaultIsometricEmbeddingType, TransparentIsometricEmbedding
export AbstractVectorTransportMethod, ParallelTransport, ProjectionTransport
export AbstractRetractionMethod,
    CayleyRetraction,
    ExponentialRetraction,
    QRRetraction,
    PolarRetraction,
    ProjectionRetraction,
    SoftmaxRetraction,
    PadeRetraction,
    ProductRetraction,
    PowerRetraction
export AbstractInverseRetractionMethod,
    ApproximateInverseRetraction,
    LogarithmicInverseRetraction,
    QRInverseRetraction,
    PolarInverseRetraction,
    ProjectionInverseRetraction,
    SoftmaxInverseRetraction
export AbstractEstimationMethod,
    GradientDescentEstimation,
    CyclicProximalPointEstimation,
    GeodesicInterpolation,
    GeodesicInterpolationWithinRadius,
    ExtrinsicEstimation
export CachedBasis,
    DefaultBasis,
    DefaultOrthogonalBasis,
    DefaultOrthonormalBasis,
    DiagonalizingOrthonormalBasis,
    InducedBasis,
    ProjectedOrthonormalBasis
export ComponentManifoldError, CompositeManifoldError
export ×,
    allocate,
    allocate_result,
    base_manifold,
    bundle_projection,
    check_point,
    check_vector,
    christoffel_symbols_first,
    christoffel_symbols_second,
    christoffel_symbols_second_jacobian,
    convert,
    complex_dot,
    decorated_manifold,
    det_local_metric,
    distance,
    dual_basis,
    einstein_tensor,
    embed,
    embed!,
    exp,
    exp!,
    flat,
    flat!,
    gaussian_curvature,
    geodesic,
    get_default_atlas,
    get_component,
    get_embedding,
    hat,
    hat!,
    induced_basis,
    incident_log,
    injectivity_radius,
    inner,
    inverse_local_metric,
    inverse_retract,
    inverse_retract!,
    isapprox,
    is_group_decorator,
    is_decorator_transparent,
    is_default_metric,
    is_default_decorator,
    is_point,
    is_vector,
    isapprox,
    kurtosis,
    local_metric,
    local_metric_jacobian,
    log,
    log!,
    log_local_metric_density,
    manifold_dimension,
    metric,
    mean,
    mean!,
    mean_and_var,
    mean_and_std,
    median,
    median!,
    mid_point,
    mid_point!,
    minkowski_metric,
    moment,
    norm,
    normal_tvector_distribution,
    number_eltype,
    one,
    power_dimensions,
    project,
    project!,
    projected_distribution,
    real_dimension,
    ricci_curvature,
    ricci_tensor,
    representation_size,
    retract,
    retract!,
    riemann_tensor,
    set_component!,
    sharp,
    sharp!,
    shortest_geodesic,
    skewness,
    std,
    sym_rem,
    submanifold,
    submanifold_component,
    submanifold_components,
    uniform_distribution,
    var,
    vector_space_dimension,
    vector_transport_along,
    vector_transport_along!,
    vector_transport_direction,
    vector_transport_direction!,
    vector_transport_to,
    vector_transport_to!,
    vee,
    vee!,
    zero_vector,
    zero_vector!
# Lie group types & functions
export AbstractGroupAction,
    AbstractGroupOperation,
    AbstractGroupManifold,
    ActionDirection,
    AdditionOperation,
    CircleGroup,
    GeneralLinear,
    GroupManifold,
    GroupOperationAction,
    Identity,
    InvariantMetric,
    LeftAction,
    LeftInvariantMetric,
    MultiplicationOperation,
    ProductGroup,
    ProductOperation,
    RightAction,
    RightInvariantMetric,
    RotationAction,
    SemidirectProductGroup,
    SpecialEuclidean,
    SpecialLinear,
    SpecialOrthogonal,
    TranslationGroup,
    TranslationAction
export affine_matrix,
    apply,
    apply!,
    apply_diff,
    apply_diff!,
    base_group,
    center_of_orbit,
    has_approx_invariant_metric,
    compose,
    compose!,
    direction,
    g_manifold,
    group_exp,
    group_exp!,
    group_log,
    group_log!,
    has_biinvariant_metric,
    has_invariant_metric,
    identity,
    identity!,
    inv,
    inv!,
    invariant_metric_dispatch,
    inverse_apply,
    inverse_apply!,
    inverse_apply_diff,
    inverse_apply_diff!,
    inverse_translate,
    inverse_translate!,
    inverse_translate_diff,
    inverse_translate_diff!,
    make_identity,
    optimal_alignment,
    optimal_alignment!,
    screw_matrix,
    switch_direction,
    translate,
    translate!,
    translate_diff,
    translate_diff!
# Orthonormal bases
export AbstractBasis,
    AbstractOrthonormalBasis,
    DefaultOrthonormalBasis,
    DiagonalizingOrthonormalBasis,
    ProjectedOrthonormalBasis,
    CachedBasis,
    DiagonalizingBasisData,
    ProductBasisData,
    PowerBasisData
export OutOfInjectivityRadiusError
export get_basis,
    get_coordinates, get_coordinates!, get_vector, get_vector!, get_vectors, number_system
# differentiation
export AbstractDiffBackend,
    AbstractRiemannianDiffBackend, FiniteDifferencesBackend, RiemannianONBDiffBackend
export diff_backend, diff_backend!, diff_backends
# atlases and charts
export get_point, get_point!, get_parameters, get_parameters!

end # module


text:
using CellularPotts
using GLMakie

M = ModelParameters(
    graphDimension = (100,100),
    isPeriodic = true,
    cellTypes = ["SimpleCell"],
    cellCounts = [1],
    cellVolumes = [500],
    penalties = [
        AdhesionPenalty([0 20; 20 50]),
        VolumePenalty([500],[1])
    ],
    temperature = 10.0
)

CPM = CellPotts(M)

CellGUI(CPM)


text:
# - Julia: High Performance Programming
# - Module 1 Chapter 2

using Allnix: add, power, fpower, fplus
using Base

@testset "add" begin
# int comparison
@test add(1,2) == 3

# - Floating point comparison - #

# type \approx and use tab complete to get ≈
@test add(1.0, 2.0) ≈ 3.0
# or use isapprox() function
@test isapprox(3.0, add(1.0,2.0))

# Exact floating point in some cases
@test add(1.0, 2.0) == 3.0
end

@testset "power" begin
@test power(2,2) == 4
@test power(2.0, 2.0) == 4.0
@test fpower(2.0, 2.0) == 4.0
end

@testset "fplus" begin
@test fplus(1,2,3) == 6
end

# ---------------------- #
# - Variable arguments - #
# ---------------------- #
@testset "varargs" begin
args = [1,2]
# Notice the ellipses
@test add(args...) == 3
end

# --------------------- #
# - Ranges and arrays - #
# --------------------- #
@testset "ones" begin
length = 100
v = ones(Float64, length)
@test Base.length(v) == length
for i in v
    @test i == 1.0
end
end

@testset "zeros" begin
length = 100
v = zeros(Float64, length)
@test Base.length(v) == length
for i in v
    @test i == 0.0
end
end

@testset "linspace" begin
start = 0
stop = 100
# IMPORTANT: n is the number of nodes not intervals
n = 201

@test linspace(start,stop,n) == 0.:0.5:100.
end

@testset "fill" begin
v = Vector{Float64}()
Base.resize!(v, 100)
### - fill! - ###
fill!(v, 13.0)
for i in v
    @test i == 13.0
end
end

# ------------------------------------------ #
# - Deep copy of array with primitive type - #
# ------------------------------------------ #
@testset "copy" begin
count = 100
v = Vector{Float64}(linspace(1,100,count))
w = copy(v)
# - Expect deep copy - #
@test pointer_from_objref(v) != pointer_from_objref(w)
@test length(w) == count
for i in 1:count
    @test w[i] == v[i]
end
end

@testset "Common Array Operators" begin
count = 100
# source = [1, 2, ... 100]
source = Vector{Float64}(linspace(1,100.,count))

####  append!  ####
v = copy(source)
append!(v, 13.0)
@test length(v) == 1 + count
@test v[1+count] == 13.0

v = copy(source)
b = copy(source)
append!(v, b)
@test length(v) == 2*count
@test v[2] == 2.0
@test v[101] == 1.0
@test v[200] == 100

### push! ###
v = copy(source)
w = push!(v, 101) # return v reference
@test length(v) == count + 1
# w and v point to the same array
@test pointer_from_objref(w) == pointer_from_objref(v)

#### pop! ####
# - Remove the last item and returns it
v = copy(source)
@test pop!(v) == 100.0
@test length(v) == count - 1

### shift! ###
# - Remove the first item and returns it
v = copy(source)
@test shift!(v) == 1.0
# v = 2,...10
@test length(v) == count - 1

### splice! ###
v = copy(source)
splice!(v, 55)
@test v[55] == 56

### in ###
v = copy(source)
@test in(55, v)

### sort ###
v = Vector{Float64}(linspace(100.,1.,count))
@test v[1] == 100.0
sort!(v) # mutable version
@test v[1] == 1.0

v = Vector{Float64}(linspace(100.,1.,count))
@test v[1] == 100.0
t = sort(v) # immutable version
@test v[1] == 100.0
@test t[1] == 1.0
end

@testset "dot operator" begin
v = Vector{Float64}()
Base.resize!(v, 100)

# - Assign value to every element - #
v .= 13.0

for i in v
    @test i == 13.0
end

@testset "pointer" begin
a::Vector{Float64} = 1:10.0::Float64
b = a
@test b == a
### pointer_from_objref ###
a_ptr = pointer_from_objref(a)
b_ptr = pointer_from_objref(b)
@test b_ptr == a_ptr

c::Vector{Float64} = 1:10.0::Float64
# compare values
@test b == c
c_ptr = pointer_from_objref(c)
@test b_ptr != c_ptr

@test b[10] == 10.0
a[10] = 13.0
@test b[10] == 13.0
@test b == a

end

end


text:
type FDEuropeanEngine{B <: AbstractBlackScholesProcess} <: AbstractFDVanillaEngine
  process::B
  timeSteps::Int
  gridPoints::Int
  timeDependent::Bool
  exerciseDate::Date
  finiteDifferenceOperator::TridiagonalOperator
  intrinsicValues::SampledCurve
  BCs::Vector{BoundaryCondition}
  sMin::Float64
  center::Float64
  sMax::Float64
  prices::SampledCurve
  fdEvolverFunc::Function
end

function FDEuropeanEngine(process::AbstractBlackScholesProcess, fdEvolverFunc::Function, timeSteps::Int = 100,
                          gridPoints::Int = 100, timeDependent::Bool = false)
  prices = SampledCurve(gridPoints)
  finiteDifferenceOperator, intrinsicValues, BCs = gen_fd_vanilla_engine_params(gridPoints)
  sMin = center = sMax = 0.0
  exerciseDate = Date()

  return FDEuropeanEngine(process, timeSteps, gridPoints, timeDependent, exerciseDate, finiteDifferenceOperator,
                          intrinsicValues, BCs, sMin, center, sMax, prices, fdEvolverFunc)
end

type FDBermudanEngine{B <: AbstractBlackScholesProcess} <: FDMultiPeriodEngine
  process::B
  timeSteps::Int
  gridPoints::Int
  timeDependent::Bool
  exerciseDate::Date
  finiteDifferenceOperator::TridiagonalOperator
  intrinsicValues::SampledCurve
  BCs::Vector{BoundaryCondition}
  sMin::Float64
  center::Float64
  sMax::Float64
  prices::SampledCurve
  stoppingTimes::Vector{Float64}
  timeStepPerPeriod::Int
  fdEvolverFunc::Function
end

function FDBermudanEngine(process::AbstractBlackScholesProcess, fdEvolverFunc::Function, timeSteps::Int = 100,
                          gridPoints::Int = 100, timeDependent::Bool = false)
  prices = SampledCurve(gridPoints)
  finiteDifferenceOperator, intrinsicValues, BCs = gen_fd_vanilla_engine_params(gridPoints)
  sMin = center = sMax = 0.0
  exerciseDate = Date()
  stoppingTimes = Vector{Float64}()
  timeStepPerPeriod = timeSteps

  return FDBermudanEngine(process, timeSteps, gridPoints, timeDependent, exerciseDate, finiteDifferenceOperator,
                          intrinsicValues, BCs, sMin, center, sMax, prices, stoppingTimes, timeStepPerPeriod, fdEvolverFunc)
end

type FDAmericanEngine{B <: AbstractBlackScholesProcess} <: FDStepConditionEngine
  process::B
  timeSteps::Int
  gridPoints::Int
  timeDependent::Bool
  exerciseDate::Date
  finiteDifferenceOperator::TridiagonalOperator
  intrinsicValues::SampledCurve
  BCs::Vector{BoundaryCondition}
  sMin::Float64
  center::Float64
  sMax::Float64
  prices::SampledCurve
  controlPrices::SampledCurve
  controlBCs::Vector{BoundaryCondition}
  controlOperator::TridiagonalOperator
  fdEvolverFunc::Function
end

function FDAmericanEngine(process::AbstractBlackScholesProcess, fdEvolverFunc::Function, timeSteps::Int = 100,
                          gridPoints::Int = 100, timeDependent::Bool = false)
  # build engine
  prices = SampledCurve()
  controlPrices = SampledCurve(gridPoints)
  finiteDifferenceOperator, intrinsicValues, BCs = gen_fd_vanilla_engine_params(gridPoints)
  sMin = center = sMax = 0.0
  exerciseDate = Date()
  controlBCs = Vector{BoundaryCondition}(2)
  controlOperator = TridiagonalOperator()

  return FDAmericanEngine(process, timeSteps, gridPoints, timeDependent, exerciseDate, finiteDifferenceOperator,
                          intrinsicValues, BCs, sMin, center, sMax, prices, controlPrices, controlBCs,
                          controlOperator, fdEvolverFunc)
end

function gen_fd_vanilla_engine_params(gridPoints::Int)
  finiteDifferenceOperator = TridiagonalOperator()
  intrinsicValues = SampledCurve(gridPoints)
  BCs = Vector{BoundaryCondition}(2)

  return finiteDifferenceOperator, intrinsicValues, BCs
end

get_residual_time(pe::AbstractFDVanillaEngine) = get_time(pe.process, pe.exerciseDate)

function safe_grid_points(gridPoints::Int, residualTime::Float64)
  minGridPoints = 10
  minGridPointsPerYear = 2
  return max(gridPoints, residualTime > 1.0 ? round(Int, floor(minGridPoints + (residualTime - 1.0) * minGridPointsPerYear)) : minGridPoints)
end

function set_grid_limits!(pe::AbstractFDVanillaEngine, center::Float64, t::Float64)
  center > 0.0 || error("negative or null underlying given")
  t > 0.0 || error("negative or zero residual time")

  pe.center = center
  newGridPoints = safe_grid_points(pe.gridPoints, t)
  if newGridPoints > QuantLib.Math.get_size(pe.intrinsicValues)
    pe.intrinsicValues = SampledCurve(newGridPoints)
  end

  volSqrtTime = sqrt(black_variance(pe.process.blackVolatility, t, pe.center))

  # the prefactor fine-tunes performance at small volatilities
  prefactor = 1.0 + 0.02 / volSqrtTime
  minMaxFactor = exp(4.0 * prefactor * volSqrtTime)

  pe.sMin = pe.center / minMaxFactor
  pe.sMax = pe.center * minMaxFactor

  return pe
end

function ensure_strike_in_grid!(pe::AbstractFDVanillaEngine, opt::VanillaOption)
  striked_payoff = opt.payoff
  requiredGridValue = striked_payoff.strike
  safetyZoneFactor = 1.1

  if pe.sMin > requiredGridValue / safetyZoneFactor
    pe.sMin = requiredGridValue / safetyZoneFactor
    # enforce central placement of underlying
    pe.sMax = pe.center / (pe.sMin / pe.center)
  end

  if pe.sMax < requiredGridValue * safetyZoneFactor
    pe.sMax = requiredGridValue * safetyZoneFactor
    # enforce central placement of underlying
    pe.sMin = pe.center / (pe.sMax / pe.center)
  end

  return pe
end

function set_grid_limits!(pe::AbstractFDVanillaEngine, opt::VanillaOption)
  set_grid_limits!(pe, state_variable(pe.process).value, get_residual_time(pe))
  ensure_strike_in_grid!(pe, opt)

  return pe
end

function initialize_initial_condition!(pe::AbstractFDVanillaEngine, opt::VanillaOption)
  set_log_grid!(pe.intrinsicValues, pe.sMin, pe.sMax)
  sample!(pe.intrinsicValues, opt.payoff)

  return pe
end

function initialize_boundary_conditions!(pe::AbstractFDVanillaEngine)
  pe.BCs[1] = build_NeumannBC(pe.intrinsicValues.values[2] - pe.intrinsicValues.values[1], LowerSide())
  pe.BCs[2] = build_NeumannBC(pe.intrinsicValues.values[end] - pe.intrinsicValues.values[end - 1], UpperSide())

  return pe
end

initialize_operator!(pe::AbstractFDVanillaEngine) = pe.finiteDifferenceOperator = get_operator(pe.process, pe.intrinsicValues.grid, get_residual_time(pe), pe.timeDependent)

function _calculate!(pe::FDEuropeanEngine, opt::EuropeanOption)
  pe.exerciseDate = opt.exercise.dates[end]
  payoff = opt.payoff

  set_grid_limits!(pe, opt)
  initialize_initial_condition!(pe, opt)
  initialize_operator!(pe)
  initialize_boundary_conditions!(pe)

  model = FiniteDifferenceModel(pe.finiteDifferenceOperator, pe.BCs, pe.fdEvolverFunc)

  pe.prices = copy(pe.intrinsicValues)

  rollback!(model, pe.prices.values, get_residual_time(pe), 0.0, pe.timeSteps)
  opt.results.value = value_at_center(pe.prices)
  opt.results.delta = first_derivative_at_center(pe.prices)
  opt.results.gamma = second_derivative_at_center(pe.prices)
  opt.results.theta = black_scholes_theta(pe.process, opt.results.value, opt.results.delta, opt.results.gamma)
  return pe, opt
end

function setup_args!(pe::FDMultiPeriodEngine, opt::VanillaOption)
  n = length(opt.exercise.dates)

  pe.stoppingTimes = zeros(n)
  for i in eachindex(pe.stoppingTimes)
    pe.stoppingTimes[i] = get_time(pe.process, opt.exercise.dates[i])
  end

  return pe
end

get_dividend_time(pe::FDMultiPeriodEngine, idx::Int) = pe.stoppingTimes[idx]

function execute_intermediate_step!(pe::FDBermudanEngine, ::Int)
  for i in eachindex(pe.intrinsicValues.grid)
    pe.prices.values[i] = max(pe.prices.values[i], pe.intrinsicValues.values[i])
  end

  return pe
end

function _calculate!(pe::FDMultiPeriodEngine, opt::VanillaOption)
  setup_args!(pe, opt)
  pe.exerciseDate = opt.exercise.dates[end]
  payoff = opt.payoff

  dateNumber = length(pe.stoppingTimes)
  lastDateIsResTime = false
  firstIndex = 0
  lastIndex = dateNumber
  firstDateIsZero = false
  firstNonZeroDate = get_residual_time(pe)

  dateTolerance = 1e-6

  resid_time = get_residual_time(pe)

  if dateNumber > 0
    get_dividend_time(pe, 1) >= 0 || error("first date cannot be negative")
    if get_dividend_time(pe, 1) < resid_time * dateTolerance
      firstDateIsZero = true
      firstIndex = 1
      if dateNumber >= 2
        firstNonZeroDate = get_dividend_time(pe, 2)
      end
    end

    if abs(get_dividend_time(pe, lastIndex) - resid_time) < dateTolerance
      lastDateIsResTime = true
      lastIndex = dateNumber - 1
    end

    if ~firstDateIsZero
      firstNonZeroDate = get_dividend_time(pe, 1)
    end

    if dateNumber >= 2
      issorted(pe.stoppingTimes) || error("dates must be in increasing order")
    end
  end

  dt = resid_time / (pe.timeStepPerPeriod * (dateNumber + 1.0))

  # ensure that dt is always smaller than the first non-zero date
  if firstNonZeroDate <= dt
    dt = firstNonZeroDate / 2.0
  end

  set_grid_limits!(pe, opt)
  initialize_initial_condition!(pe, opt)
  initialize_operator!(pe)
  initialize_boundary_conditions!(pe)
  model = FiniteDifferenceModel(pe.finiteDifferenceOperator, pe.BCs, pe.fdEvolverFunc)

  pe.prices = copy(pe.intrinsicValues)

  if lastDateIsResTime
    execute_intermediate_step!(pe, dateNumber - 1)
  end

  j = lastIndex
  while true
    if j == dateNumber
      beginDate = get_residual_time(pe)
    else
      beginDate = get_dividend_time(pe, j+1)
    end

    if j >= 1
      endDate = get_dividend_time(pe, j)
    else
      endDate = dt
    end

    rollback!(model, pe.prices.values, beginDate, endDate, pe.timeStepPerPeriod)

    if j >= 1
      execute_intermediate_step!(pe, j)
    end

    j -= 1
    if j < firstIndex
      break
    end
  end

  rollback!(model, pe.prices.values, dt, 0.0, 1)

  if firstDateIsZero
    execute_intermediate_step!(pe, 0)
  end

  opt.results.value = value_at_center(pe.prices)
  opt.results.delta = first_derivative_at_center(pe.prices)
  opt.results.gamma = second_derivative_at_center(pe.prices)

  return pe, opt
end

initialize_step_condition(pe::FDAmericanEngine) = build_AmericanStepCondition(pe.intrinsicValues.values)

function _calculate!(pe::FDStepConditionEngine, opt::VanillaOption)
  pe.exerciseDate = opt.exercise.dates[end]
  payoff = opt.payoff

  set_grid_limits!(pe, opt)
  initialize_initial_condition!(pe, opt)
  initialize_operator!(pe)
  initialize_boundary_conditions!(pe)
  stepCond = initialize_step_condition(pe)

  pe.prices = copy(pe.intrinsicValues)
  pe.controlPrices = copy(pe.intrinsicValues)
  pe.controlOperator = copy(pe.finiteDifferenceOperator)
  pe.controlBCs[1] = pe.BCs[1]
  pe.controlBCs[2] = pe.BCs[2]

  operatorSet = TridiagonalOperator[pe.finiteDifferenceOperator, pe.controlOperator]
  arraySet = Vector{Float64}[x for x in (pe.prices.values, pe.controlPrices.values)]
  bcSet = hcat(pe.BCs, pe.controlBCs)
  conditionSet = StepConditionSet(StepCondition[stepCond, FdmNullStepCondition()])

  model = FiniteDifferenceModel(operatorSet, bcSet, pe.fdEvolverFunc)

  rollback!(model, arraySet, get_residual_time(pe), 0.0, pe.timeSteps, conditionSet)

  pe.prices.values = arraySet[1]
  pe.controlPrices.values = arraySet[2]

  variance = black_variance(pe.process.blackVolatility, opt.exercise.dates[end], payoff.strike)

  dividendDiscount = discount(pe.process.dividendYield, opt.exercise.dates[end])

  riskFreeDiscount = discount(pe.process.riskFreeRate, opt.exercise.dates[end])

  spot = state_variable(pe.process).value
  forwardPrice = spot *  dividendDiscount / riskFreeDiscount

  black = BlackCalculator(payoff, forwardPrice, sqrt(variance), riskFreeDiscount)

  opt.results.value = value_at_center(pe.prices) - value_at_center(pe.controlPrices) + value(black)
  opt.results.delta = first_derivative_at_center(pe.prices) - first_derivative_at_center(pe.controlPrices) + delta(black, spot)
  opt.results.gamma = second_derivative_at_center(pe.prices) - second_derivative_at_center(pe.prices) + gamma(black, spot)

  return pe, opt
end


text:
@everywhere using PowerModels # Uses Julia package for Steady-State Power Network Optimization.

##### Generic Helper Functions #####

function remove_comment(string)
    return split(string, "/")[1] # Split str into an array of substrings on occurrences of the delimiter (eg ,)
end

##### GOC Initialization File Parser (.ini) #####

function parse_goc_files(ini_file; scenario_id="") # Function to parse the input files
    files, scenario_id = find_goc_files(ini_file, scenario_id=scenario_id) # Calls the function named "find_goc_files"
    return parse_goc_files(files["con"], files["inl"], files["raw"], files["rop"], ini_file=ini_file, scenario_id=scenario_id) # Returns the output of function "parse_goc_files"
end

function find_goc_files(ini_file; scenario_id="")
    files = Dict(
        "rop" => "x",
        "raw" => "x",
        "con" => "x",
        "inl" => "x"
    )

    if !endswith(ini_file, ".ini") # Checks for the format of file
        warn(LOGGER, "given init file does not end with .ini, $(ini_file)") # Gives warning if the file format is not correct
    end

    open(ini_file) do io # Opens the ini_file
        for line in readlines(io)
            line = strip(line) # Strip fucntion is used to trim both leading & trailing white space.
            #println(line)
            if startswith(line, "[INPUTS]")
                # do nothing
            elseif startswith(line, "ROP")
                files["rop"] = strip(split(line,"=")[2]) # Split returns a list of the words in the string, using sep as the delimiter string
            elseif startswith(line, "RAW")
                files["raw"] = strip(split(line,"=")[2])
            elseif startswith(line, "CON")
                files["con"] = strip(split(line,"=")[2])
            elseif startswith(line, "INL")
                files["inl"] = strip(split(line,"=")[2])
            else
                warn(LOGGER, "unknown input given in ini file: $(line)") # Warning if none of the conditions satisfied
            end
        end
    end

    #println(files)

    ini_dir = dirname(ini_file)

    #println(ini_dir)
    scenario_dirs = [file for file in readdir(ini_dir) if isdir(joinpath(ini_dir, file))]
    scenario_dirs = sort(scenario_dirs) # Sort the scenario_dirs
    #println(scenario_dirs)

    if length(scenario_id) == 0 # Checks the length of scenario_id if it is equal to zero
        scenario_id = scenario_dirs[1] # Assigns the value of scenario_dirs[1] to scenario_id
        info(LOGGER, "no scenario specified, selected directory \"$(scenario_id)\"") # If scenario_id is zero then no scenario is specified.
    else
        if !(scenario_id in scenario_dirs) # Checks whether scenario_id is in scenario_dirs
            error(LOGGER, "$(scenario_id) not found in $(scenario_dirs)") # Error message if scenario_id not found.
        end
    end

    for (id, path) in files
        if path == "."
            files[id] = ini_dir # Assigns ini_dir to files[id]
        elseif path == "x"
            files[id] = joinpath(ini_dir, scenario_id) # Join path components into a full path.
        else
            error(LOGGER, "unknown file path directive $(path) for file $(id)") # Error (unknown file path) if none of the above conditions fulfilled
        end
    end

    files["raw"] = joinpath(files["raw"], "case.raw")
    files["rop"] = joinpath(files["rop"], "case.rop")
    files["inl"] = joinpath(files["inl"], "case.inl")
    files["con"] = joinpath(files["con"], "case.con")

    return files, scenario_id
end

@everywhere function parse_goc_files(con_file, inl_file, raw_file, rop_file; ini_file="", scenario_id="none") # Function that performs the parsing of goc files
    files = Dict(
        "rop" => rop_file,
        "raw" => raw_file,
        "con" => con_file,
        "inl" => inl_file
    )

    info(LOGGER, "Parsing Files") # To notify the current situation of the program
    info(LOGGER, "  raw: $(files["raw"])") # raw: raw_file
    info(LOGGER, "  rop: $(files["rop"])") # rop: rop_file
    info(LOGGER, "  inl: $(files["inl"])") # inl: inl_file
    info(LOGGER, "  con: $(files["con"])") # con: con_file

    info(LOGGER, "skipping power models data warnings") # Informs what things are skipped
    pm_logger_level = getlevel(getlogger(PowerModels)) # "getlogger" returns the current root/global logger (to record events)
    setlevel!(getlogger(PowerModels), "error")
    network_model = PowerModels.parse_file(files["raw"], import_all=true) # Parse the raw file
    setlevel!(getlogger(PowerModels), pm_logger_level)

    gen_cost = parse_rop_file(files["rop"]) # parsing the rop file
    response = parse_inl_file(files["inl"]) # parsing the inl file
    contingencies = parse_con_file(files["con"]) # parsing the con file

    return (ini_file=ini_file, scenario=scenario_id, network=network_model, cost=gen_cost, response=response, contingencies=contingencies, files=files) # Returns these values
end

function parse_goc_opf_files(ini_file; scenario_id="") # Function to parse the opf files
    files = Dict(
        "rop" => "x",
        "raw" => "x",
    )

    if !endswith(ini_file, ".ini") # Chcks the file format
        warn(LOGGER, "given init file does not end with .ini, $(ini_file)") # Gives warning due to wrong file format
    end

    open(ini_file) do io # Open the ini_file and allow read/write permissions
        for line in readlines(io) # For loop to read lines
            line = strip(line) # Removes leading and trailing characters from line
            #println(line)
            if startswith(line, "[INPUTS]") # Returns true if line starts with INPUTS.
                # do nothing
            elseif startswith(line, "ROP") # Returns true if line starts with ROP.
                files["rop"] = strip(split(line,"=")[2]) # Split: splits the characters on a new line
            elseif startswith(line, "RAW") # Returns true if line starts with RAW.
                files["raw"] = strip(split(line,"=")[2]) # Split: splits the characters on a new line
            else
                warn(LOGGER, "unknown input given in ini file: $(line)") # Warning if none of the above conditions satisfied
            end
        end
    end

    #println(files)

    ini_dir = dirname(ini_file) # Assigns the directory name of ini_file to ini_dir

    #println(ini_dir)
    scenario_dirs = [file for file in readdir(ini_dir) if isdir(joinpath(ini_dir, file))] # Reads the ini_dir and if directory is present then joins the path with file
    scenario_dirs = sort(scenario_dirs) # To sort the scenario_dirs
    #println(scenario_dirs)

    if length(scenario_id) == 0 # Checks the length of scenario_id if it is equal to zero
        scenario_id = scenario_dirs[1] # Assigns the value of scenario_dirs[1] to scenario_id
        info(LOGGER, "no scenario specified, selected directory \"$(scenario_id)\"") # If scenario_id is zero then no scenario is specified
    else
        if !(scenario_id in scenario_dirs) # Checks whether scenario_id is in scenario_dirs
            error(LOGGER, "$(scenario_id) not found in $(scenario_dirs)") # Error message if scenario_id not found.
        end
    end

    for (id, path) in files
        if path == "."
            files[id] = ini_dir  # Assigns ini_dir to files[id]
        elseif path == "x"
            files[id] = joinpath(ini_dir, scenario_id) # Join path components into a full path.
        else
            error(LOGGER, "unknown file path directive $(path) for file $(id)") # Error (unknown file path) if none of the above conditions fulfilled
        end
    end

    files["raw"] = joinpath(files["raw"], "case.raw")
    files["rop"] = joinpath(files["rop"], "case.rop")

    info(LOGGER, "Parsing Files")
    info(LOGGER, "  raw: $(files["raw"])")
    info(LOGGER, "  rop: $(files["rop"])")

    network_model = PowerModels.parse_file(files["raw"], import_all=true) # Parsing the raw file to know details about network model
    gen_cost = parse_rop_file(files["rop"]) # Parsing the rop file to get gen_cost

    return (ini_file=ini_file, scenario=scenario_id, network=network_model, cost=gen_cost, files=files) # Returns the specified values
end

##### Unit Inertia and Governor Response Data File Parser (.inl) #####

@everywhere function parse_inl_file(file::String) # Function to parse the .inl file
    open(file) do io # Open the file for read/write access
        return parse_inl_file(io) # Return the parsed file
    end
end

@everywhere function parse_inl_file(io::IO) # Function including the details for parsing the inl file
    inl_list = []
    for line in readlines(io)
        #line = remove_comment(line)

        if startswith(strip(line), "0") # Checks if the line starts with zero
            debug(LOGGER, "inl file sentinel found")
            break
        end
        line_parts = split(line, ",") # Split the line based on commas
        @assert length(line_parts) >= 7 # Assert macro allows the users to optionally specify their own error message, instead of just printing the failed expression

        inl_data = Dict(
            "i"    => parse(Int, line_parts[1]),
            "id"   => strip(line_parts[2]),
            "h"    => strip(line_parts[3]),
            "pmax" => strip(line_parts[4]),
            "pmin" => strip(line_parts[5]),
            "r"    => parse(Float64, line_parts[6]),
            "d"    => strip(line_parts[7])
        )

        @assert inl_data["r"] >= 0.0 # Assert macro allows the users to optionally specify their own error message, instead of just printing the failed expression

        #println(inl_data)
        push!(inl_list, inl_data) # Inserts inl_data at the end of inl_list
    end
    return inl_list
end

##### Generator Cost Data File Parser (.rop) #####

@everywhere rop_sections = [
    "mod" => "Modification Code",
    "bus_vm" => "Bus Voltage Attributes",
    "shunt_adj" => "Adjustable Bus Shunts",
    "load" => "Bus Loads",
    "load_adj" => "Adjustable Bus Load Tables",
    "gen" => "Generator Dispatch Units",
    "disptbl" => "Active Power Dispatch Tables",
    "gen_reserve" => "Generator Reserve Units",
    "qg" => "Generation Reactive Capability",
    "branch_x" => "Adjustable Branch Reactance",
    "ctbl" => "Piecewise Linear Cost Curve Tables",
    "pwc" => "Piecewise Quadratic Cost Curve Tables",
    "pec" => "Polynomial & Exponential Cost Curve Tables",
    "reserve" => "Period Reserves",
    "branch_flow" => "Branch Flows",
    "int_flow" => "Interface Flows",
    "lin_const" => "Linear Constraint Dependencies",
    "dc_const" => "Two Terminal DC Line Constraint Dependencies",
]

@everywhere function parse_rop_file(file::String) # Calls the function named parse_rop_file
    open(file) do io # Opens the file for read/write operations
        return parse_rop_file(io)
    end
end

@everywhere function parse_rop_file(io::IO) # Function including the details for parsing the rop file
    active_section_idx = 1
    active_section = rop_sections[active_section_idx]

    section_data = Dict()
    section_data[active_section.first] = []

    line_idx = 1
    lines = readlines(io) # Read from io
    while line_idx < length(lines) # Compare the length of lines and line index
        #line = remove_comment(lines[line_idx])
        line = lines[line_idx]
        if startswith(strip(line), "0") # Checks if the stripped line starts with zero
            debug(LOGGER, "finished reading rop section $(active_section.second) with $(length(section_data[active_section.first])) items")
            active_section_idx += 1 # active_section_idx = active_section_idx + 1
            if active_section_idx > length(rop_sections) # If statement to make comparison b/w rop_sections & active_section_idx
                debug(LOGGER, "finished reading known rop sections")
                break
            end
            active_section = rop_sections[active_section_idx]
            section_data[active_section.first] = []
            line_idx += 1
            continue
        end

        if active_section.first == "gen" # Checks if the first entity in active_section is gen.
            push!(section_data[active_section.first], _parse_rop_gen(line)) # Inserts _parse_rop_gen(line) into the section_data file
        elseif active_section.first == "disptbl" # Checks if the first entity in active_section is disptbl.
            push!(section_data[active_section.first], _parse_rop_pg(line)) # Inserts _parse_rop_gen(line) into the section_data file
        elseif active_section.first == "ctbl" # Checks if the first entity in active_section is ctbl.
            pwl_line_parts = split(line, ",") # Splits the line based on commas
            @assert length(pwl_line_parts) >= 3

            num_pwl_lines = parse(Int, pwl_line_parts[3])
            @assert num_pwl_lines > 0 # Assert macro allows the users to optionally specify their own error message, instead of just printing the failed expression

            pwl_point_lines = lines[line_idx+1:line_idx+num_pwl_lines]
            #pwl_point_lines = remove_comment.(pwl_point_lines)
            push!(section_data[active_section.first], _parse_rop_pwl(pwl_line_parts, pwl_point_lines)) # Inserts _parse_rop_pwl into the section_data file
            line_idx += num_pwl_lines # line_idx = line_idx + num_pwl_lines
        else
            info(LOGGER, "skipping data line: $(line)") # If none of the above statements satisfied the this info appears.
        end
        line_idx += 1 # line_idx = line_idx + 1
    end
    return section_data # Returns the content of section_data
end

@everywhere function _parse_rop_gen(line) # Calls the function named _parse_rop_gen
    line_parts = split(line, ",") # Splits the line based on comma
    @assert length(line_parts) >= 4 # Assert macro allows the users to optionally specify their own error message, instead of just printing the failed expression

    data = Dict(
        "bus"     => parse(Int, line_parts[1]),
        "genid"   => strip(line_parts[2]),
        "disp"    => strip(line_parts[3]),
        "disptbl" => parse(Int, line_parts[4]),
    )

    @assert data["disptbl"] >= 0 # Assert macro allows the users to optionally specify their own error message, instead of just printing the failed expression

    return data
end

@everywhere function _parse_rop_pg(line)
    line_parts = split(line, ",")
    @assert length(line_parts) >= 7

    data = Dict(
        "tbl"      => parse(Int, line_parts[1]),
        "pmax"     => strip(line_parts[2]),
        "pmin"     => strip(line_parts[3]),
        "fuelcost" => strip(line_parts[4]),
        "ctyp"     => strip(line_parts[5]),
        "status"   => strip(line_parts[6]),
        "ctbl"     => parse(Int, line_parts[7]),
    )

    @assert data["tbl"] >= 0
    @assert data["ctbl"] >= 0

    return data
end

@everywhere function _parse_rop_pwl(pwl_parts, point_lines)
    @assert length(pwl_parts) >= 2

    points = []

    for point_line in point_lines
        point_line_parts = split(point_line, ",") # Splits the point_line based on commas
        @assert length(point_line_parts) >= 2 # Assert macro allows to optionally specify our own error message, instead of just printing the failed expression
        x = parse(Float64, point_line_parts[1])
        y = parse(Float64, point_line_parts[2])

        push!(points, (x=x, y=y))
    end

    data = Dict(
        "ltbl"   =>  parse(Int, pwl_parts[1]),
        "label"  => strip(pwl_parts[2]),
        "points" => points
    )

    @assert data["ltbl"] >= 0

    return data
end

##### Contingency Description Data File (.con) #####

# OPEN BRANCH FROM BUS *I TO BUS *J CIRCUIT *1CKT
@everywhere branch_contigency_structure = [
    1 => "OPEN",
    2 => "BRANCH",
    3 => "FROM",
    4 => "BUS",
    #5 => "I",
    6 => "TO",
    7 => "BUS",
    #8 => "J",
    9 => "CIRCUIT",
    #10 => "CKT"
]

#=
# OPEN BRANCH FROM BUS *I TO BUS *J CIRCUIT *1CKT
@everywhere branch_contigency_structure_alt = [
    1 => "OPEN",
    2 => "BRANCH",
    3 => "FROM",
    4 => "BUS",
    #5 => "I",
    6 => "TO",
    7 => "BUS",
    #8 => "J",
    9 => "CKT",
    #10 => "CKT"
]
=#

# REMOVE UNIT *ID FROM BUS *I
# Generator contigency details
@everywhere generator_contigency_structure = [
    1 => "REMOVE",
    2 => "UNIT",
    #3 => "ID",
    4 => "FROM",
    5 => "BUS",
    #6 => "I"
]

@everywhere function parse_con_file(file::String) # Calls the parse_con_file
    open(file) do io # Open the file for read/write purposes
        return parse_con_file(io) # Returns the parsed file
    end
end

@everywhere function parse_con_file(io::IO) # Function that describe details regarding how to parse the contingency file
    con_lists = []

    tokens = []

    for line in readlines(io) # Reads from the io file
        #line_tokens = split(strip(remove_comment(line)))
        line_tokens = split(strip(line))
        #println(line_tokens)
        append!(tokens, line_tokens) # Inserts line_tokens at the end of tokens
    end

    #println(tokens)

    token_idx = 1
    while token_idx <= length(tokens) # Infinite loop until token_idx > length(tokens)
        token = tokens[token_idx] # Assigns tokens[token_idx] to token
        if token == "END"
            debug(LOGGER, "end of contingency file found") # Notify the end of contingency file
            break
        elseif token == "CONTINGENCY" # Checks if token is equal to CONTINGENCY
            # start reading contingencies

            contingency_name = tokens[token_idx+1] # Increments the token index to read the next contingency
            debug(LOGGER, "reading contingency $(contingency_name)") # Notify the contingency name

            token_idx += 2 # token_idx = token_idx + 2
            token = tokens[token_idx]
            remaining_tokens = length(tokens) - token_idx # Finds the remaining tokens

            if token == "OPEN" # branch contingency case
                # OPEN BRANCH FROM BUS *I TO BUS *J CIRCUIT *1CKT

                @assert remaining_tokens >= 9
                branch_tokens = tokens[token_idx:token_idx+9] # Informs how many tokens are there for a branch
                #println(branch_tokens)

                #if !all(branch_tokens[idx] == val for (idx, val) in branch_contigency_structure) && !all(branch_tokens[idx] == val for (idx, val) in branch_contigency_structure_alt)
                if any(branch_tokens[idx] != val for (idx, val) in branch_contigency_structure) # If there is contradiction b/w branch contingency structure and branch tokens
                    error(LOGGER, "incorrect branch contingency structure: $(branch_tokens)")
                end

                bus_i = parse(Int, branch_tokens[5])
                @assert bus_i >= 0

                bus_j = parse(Int, branch_tokens[8])
                @assert bus_j >= 0

                ckt = branch_tokens[10]

                branch_contingency = Dict(
                    "label" => contingency_name,
                    "component" => "branch",
                    "action" => "open",
                    "i" => bus_i,
                    "j" => bus_j,
                    "ckt" => ckt,
                )

                push!(con_lists, branch_contingency) # Inserts branch contingency in contingency list

                token_idx += 9 # token_idx = token_idx + 9
            elseif token == "REMOVE"
                # REMOVE UNIT *ID FROM BUS *I

                @assert remaining_tokens >= 5
                generator_tokens = tokens[token_idx:token_idx+5]
                #println(generator_tokens)

                if any(generator_tokens[idx] != val for (idx, val) in generator_contigency_structure) # checks if there is contradiction b/w  generator contingency structure and generator tokens
                    error(LOGGER, "incorrect generator contingency structure: $(generator_tokens)")
                end

                gen_id = generator_tokens[3]

                bus_i = parse(Int, generator_tokens[6])
                @assert bus_i >= 0

                generator_contingency = Dict(
                    "label" => contingency_name,
                    "component" => "generator",
                    "action" => "remove",
                    "id" => gen_id,
                    "i" => bus_i,
                )

                push!(con_lists, generator_contingency) # Inserts generator contingencies into the contingency list

                token_idx += 5 # token_idx = token_idx + 5
            elseif token == "END"
                warn(LOGGER, "no action provided for contingency $(contingency_name)") # If no action is mentioned in the file for a contingency
                token_idx -= 1
            else
                warn(LOGGER, "unrecognized token $(token)") # For unexpected tokens
            end

            token_idx += 1
            token = tokens[token_idx]
            if token != "END"
                error(LOGGER, "expected END token at end of CONTINGENCY, got $(token)")
            end
        else
            warn(LOGGER, "unrecognized token $(token)")
        end
        token_idx += 1
    end

    return con_lists
end

@everywhere function parse_solution1_file(file::String) # Calls function named parse_solution1_file
    open(file) do io # Opens the file for read/write operations
        return parse_solution1_file(io)
    end
end

@everywhere function parse_solution1_file(io::IO) # Function containing details how to parse the solutin1 file
    bus_data_list = []
    gen_data_list = []

    lines = readlines(io) # Reads the input

    # skip bus list header section
    idx = 1

    separator_count = 0
    skip_next = false

    while idx <= length(lines)
        line = lines[idx]
        if length(strip(line)) == 0
            warn(LOGGER, "skipping blank line in solution1 file ($(idx))")
        elseif skip_next
            skip_next = false
        elseif startswith(strip(line), "--")
            separator_count += 1
            skip_next = true
        else
            if separator_count == 1
                parts = split(line, ",")
                @assert length(parts) >= 4
                bus_data = (
                    bus = parse(Int, parts[1]),
                    vm = parse(Float64, parts[2]),
                    va = parse(Float64, parts[3]),
                    bcs = parse(Float64, parts[4])
                )
                push!(bus_data_list, bus_data) # Inserts bus data at the end of bus_data_list
            elseif separator_count == 2
                parts = split(line, ",")
                @assert length(parts) >= 4
                gen_data = (
                    bus = parse(Int, parts[1]),
                    id = strip(strip(parts[2]), ['\'', ' ']),
                    pg = parse(Float64, parts[3]),
                    qg = parse(Float64, parts[4])
                )
                push!(gen_data_list, gen_data) # Inserts the gen_data at the end of gen_data_list
            else
                warn(LOGGER, "skipping line in solution1 file ($(idx)): $(line)")
            end
        end
        idx += 1
    end

    return (bus=bus_data_list, gen=gen_data_list) # Returns bus and gen data
end


text:
using Test
using JDF, TimeZones, DataFrames

a = some_elm(ZonedDateTime)

ar = Vector{ZonedDateTime}(undef, 10)
ar .= a

meta = open("plsdel.io", "w") do io
    compress_then_write(ar, io)
end

ar_loaded = open("plsdel.io") do io
    column_loader(meta.type, io, meta)
end

@test ar == ar_loaded

df = DataFrame(ar = ar)

JDF.save(df, "df.jdf")

df_loaded = DataFrame(JDF.load("df.jdf"))

@test df == df_loaded

ar_timezone = [ar.timezone for ar in ar]
meta = open("plsdel.io", "w") do io
    compress_then_write(ar_timezone, io)
end

ar_timezone_copy = open("plsdel.io") do io
    column_loader(meta.type, io, meta)
end

@test ar_timezone == ar_timezone_copy

meta = open("plsdel.io", "w") do io
    compress_then_write(a.timezone.transitions, io)
end

a_timezone_transitions_loaded = open("plsdel.io") do io
    column_loader(meta.type, io, meta)
end

@test a.timezone.transitions == a_timezone_transitions_loaded

rm("plsdel.io", force = true)
rm("df.jdf", force = true, recursive = true)

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
using OpenStreetMapX, Graphs, PyCall

# This code assumes that folium has benn installed

function plot_map(m::MapData, filename::AbstractString; tiles="Stamen Toner" )
    MAP_BOUNDS = [ ( m.bounds.min_y, m.bounds.min_x), ( m.bounds.max_y, m.bounds.max_x) ]
    flm = pyimport("folium")
    m_plot = flm.Map(tiles=tiles)
    for e in edges(m.g)
		info = "Edge from: $(e.src) to $(e.dst)
[information from the
.e
 and
.w
 fields] "
        flm.PolyLine(  (latlon(m,e.src), latlon(m,e.dst)),
            color="brown", weight=4, opacity=1).add_to(m_plot)
    end

    for n in keys(m.nodes)
        lla = LLA(m.nodes[n],m.bounds)
        info = "Node: $(n)\n
Lattitude: $(lla.lat)\n
Longitude: $(lla.lon)
[information from the
.node
 field] "
        flm.Circle(
            (lla.lat, lla.lon),
            popup=info,
            tooltip=info,
            radius=10,
            color="orange",
            weight=3,
            fill=true,
            fill_color="orange"
          ).add_to(m_plot)
    end

	for nn in keys(m.n)
		n = m.n[nn]
        lla = LLA(m.nodes[n],m.bounds)
        info =     "Graph: $nn
Node: $(n)\n
Lattitude: $(lla.lat)\n
Longitude: $(lla.lon)

		[The node identifiers are hold in the
.n
 field and location in the
.nodes
 field]"
		flm.Rectangle(
            [(lla.lat-0.00014, lla.lon-0.0002), (lla.lat+0.00014, lla.lon+0.0002)],
            popup=info,
            tooltip=info,
            color="green",
            weight=1.5,
            fill=false,
            fill_opacity=0.2,
            fill_color="green",
        ).add_to(m_plot)
    end

    MAP_BOUNDS = [( m.bounds.min_y, m.bounds.min_x),( m.bounds.max_y, m.bounds.max_x)]
    flm.Rectangle(MAP_BOUNDS, color="black",weight=4).add_to(m_plot)
    m_plot.fit_bounds(MAP_BOUNDS)
    m_plot.save(filename)
end

pth = joinpath(dirname(pathof(OpenStreetMapX)),"..","test","data","reno_east3.osm")

m2 =  OpenStreetMapX.get_map_data(pth,use_cache = false, trim_to_connected_graph=true);

plot_map(m2, "mymap.html")


text:
function load_net2(dir::AbstractString)
    nodeFfiles = readdir(dir)
    filter!(x->ismatch(r"\w+.txt$",x),nodeFfiles)
    node_ids_byname = Dict{AbstractString,Integer}()
    node_ids_bynum = Dict{Integer,AbstractString}()
    node_tables = Dict{AbstractString,Array{AbstractString,1}}()

    for name in	nodeFfiles
        f = open(joinpath(dir,name))
        try
            node_id = split(name,['.','_'])[1:2]
            node_tables[node_id[2]] = readlines(f)
            node_number = parse(Int,node_id[1])
            node_ids_byname[node_id[2]] = node_number
            node_ids_bynum[node_number] = node_id[2]
        finally
            close(f)
        end
    end

    netstate = Vector{BNode}()
    for (node_num,node_name) in sort(collect(node_ids_bynum))
        push!(netstate,BNode(node_name))
    end
    net = Net(netstate,DiGraph(length(netstate)))
    for sigma in net.state
        test_pattern = match(r"\d+", strip(node_tables[sigma.id][1]))
        sigma.def_state = test_pattern != nothing ? parse(Int,test_pattern.match) : error("Invalid value: First line of node file $(sigma.id) must contain a numeric value to set default state")
        # set_node!(sigma)
        test_pattern = match(r"\d+", strip(node_tables[sigma.id][2]))
        max_states = test_pattern != nothing ? parse(Int,test_pattern.match) : error("Invalid value: Second line of node file $(sigma.id) must contain the maximum number of possible states")
        sigma.state_range = collect(range(0, max_states)) # All possible states is given by a sequence of consecutive integer numbers starting from 0
        push!(sigma.thresholds,(max_states-1)/max_states*sigma.state_range[2:end]...)
        reg_ids = map(x->split(x,['_','.']),matchall(r"(\d+_)?[[:alnum:]]+(\.\d+)?",strip(node_tables[sigma.id][3])))
        for reg_id in reg_ids
            if length(reg_id) == 1
                reg_id_name = reg_id[1]
            elseif length(reg_id) == 2
                reg_id_name = reg_id[2]
                push!(sigma.tau_in,0)
            elseif length(reg_id) == 3
                reg_id_name = reg_id[2]
                push!(sigma.tau_in,parse(Int,reg_id[3]))
            else
                error("Format error for list of regulators for node $(sigma.id)")
            end
            index_reg = findfirst(x -> x.id==reg_id_name,net.state)
            index_reg!=0 || error("Unregistered node in list of regulators: There is no file for node $reg_id_name")
            push!(sigma.sigma_in, index_reg)
            add_edge!(net.graph,findfirst(x -> x.id==reg_id_name,net.state),findfirst(x -> x.id==sigma.id,net.state))
        end
        for reg_phrase in node_tables[sigma.id][4:end]
            r = matchall(r"\d+",reg_phrase)
            if(!isempty(r))
                sigma.F[tuple([parse(Int,ri) for ri in r[1:(end-1)]]...)] = parse(Int,r[end])
            end
        end
    end
    for n in vertices(net.graph)
        push!(net.state[n].sigma_out,out_neighbors(net.graph,n)...)
    end
    mem_size = maximum(vcat([sigma.tau_in for sigma in net.state]...)) + 1
    for (nsigma,sigma) in enumerate(net.state)
        # sigma.state = fill(sigma.def_state,maximum([sigmak.tau_in[findfirst(sigmak.sigma_in,nsigma)] for sigmak in net.state[sigma.sigma_out]])+1)
        sigma.state = fill(sigma.def_state,mem_size)
        sigma.sigma_in_state = [net.state[sigmak].def_state for sigmak in sigma.sigma_in]
        sigma.sigma_in_state_mu = Vector{Float16}(length(sigma.sigma_in))
        # set_node!(sigma)
    end
    return net
end

function fill_tables!(netstate::Vector{BNode},p_scheme::Symbol,p::AbstractFloat)
    if p_scheme==:semifix
        for sigma in netstate
            itp = product((netstate[sigmak].state_range for sigmak in sigma.sigma_in)...)
            table_size = length(itp)
            n_ones = round(Int,p*table_size)
            table_output = trues(n_ones)
            push!(table_output,bitrand(table_size-n_ones)...)
            shuffle!(table_output)
            sigma.F = Dict(reg_phrase => table_output[i] for (i,reg_phrase) in enumerate(itp))
        end
    elseif p_scheme==:prob
        for sigma in netstate
            itp = product((netstate[sigmak].state_range for sigmak in sigma.sigma_in)...)
            table_output = rand(length(itp)) .< p
            sigma.F = Dict(reg_phrase => table_output[i] for (i,reg_phrase) in enumerate(itp))
        end
    elseif p_scheme==:fix
        for sigma in netstate
            itp = product((netstate[sigmak].state_range for sigmak in sigma.sigma_in)...)
            table_size = length(itp)
            table_output = falses(table_size)
            n_ones = round(Int,p*table_size)
            table_output[1:n_ones] = true
            shuffle!(table_output)
            sigma.F = Dict(reg_phrase => table_output[i] for (i,reg_phrase) in enumerate(itp))
        end
    else
        error("Invalid scheme name: Choose from :fix,:prob,:semifix")
    end
    return nothing
end

function fill_taus!(net_graph::DiGraph,netstate::Vector{BNode},tau_scheme::Symbol;tau_params...)
    tau = 0
    taus = [0]
    found = false
    if tau_scheme==:fix
        for (id,value) in tau_params
            if id==:tau
                tau = value
                found = true
                break
            end
        end
        found || error("Expected tau parameter")
        for (nsigma,sigma) in enumerate(netstate)
            push!(sigma.tau_in,fill(tau,length(sigma.sigma_in))...)
            push!(sigma.sigma_in_state,zeros(Int8,length(sigma.sigma_in))...)
            sigma.sigma_in_state_mu = Vector{Float16}(length(sigma.sigma_in))
            # set_node!(sigma)
        end
    elseif tau_scheme==:random
        for (id,value) in tau_params
            if id==:taus
                taus = value
                found = true
                break
            end
        end
        found || error("Expected taus parameter")
        for (nsigma,sigma) in enumerate(netstate)
            push!(sigma.tau_in,rand(taus,length(sigma.sigma_in))...)
            push!(sigma.sigma_in_state,zeros(Int8,length(sigma.sigma_in))...)
            sigma.sigma_in_state_mu = Vector{Float16}(length(sigma.sigma_in))
            # set_node!(sigma)
        end
    else
        error("Invalid scheme name: Choose from :fix,:random")
    end
    mem_size = maximum(vcat([sigma.tau_in for sigma in netstate]...)) + 1
    for sigma in netstate
        sigma.state = zeros(Int8,mem_size)
    end
    return nothing
end

function generate_RBN(N::Integer,K::Integer;p_scheme::Symbol=:fix,tau_scheme::Symbol=:fix,params...)
    p = 0.5
    tau = 0
    taus = [0]
    for (id,value) in params
        if id==:p
            p = value
        elseif id==:tau
            tau = value
        elseif id==:taus
            taus = value
        end
    end
    net_graph = random_regular_digraph(N, K; dir=:in, seed=-1)
    netstate = Vector{BNode}(N)
    findfirst([:fix,:prob,:semifix],p_scheme)!=0 || error("Invalid scheme name: Choose from :fix,:prob,:semifix")
    for n in vertices(net_graph)
        netstate[n] = BNode(join(["n",string(n)]))
    end
    println("Succesfully added ",length(netstate)," nodes")
    for n in vertices(net_graph)
        isempty(in_neighbors(net_graph,n))||push!(netstate[n].sigma_in, in_neighbors(net_graph,n)...)
        isempty(out_neighbors(net_graph,n))||push!(netstate[n].sigma_out, out_neighbors(net_graph,n)...)
    end
    fill_tables!(netstate,p_scheme,p)
    fill_taus!(net_graph,netstate,tau_scheme,tau=tau,taus=taus)
    return Net(netstate,net_graph)
end

function knock_out!(netstate::Vector{BNode},nodeKO_id::AbstractString)
    index_KO = findfirst(x->x.id==nodeKO_id,netstate)
    if index_KO!=0
        nodeKO = netstate[index_KO]
        init_node!(nodeKO,nodeKO.def_state)
        nodeKO.thresholds = [nodeKO.thresholds[findfirst(sort([nodeKO.def_state,nodeKO.thresholds...]),nodeKO.def_state)]]
        nodeKO.state_range = [nodeKO.def_state]
        for k in keys(nodeKO.F)
            nodeKO.F[k] = nodeKO.def_state
        end
    else
        error("Node id invalid: $nodeKO_id does not exist in the network")
    end
end

function knock_down!(netstate::Vector{BNode},nodeKD_id::AbstractString)
    index_KD = findfirst(x->x.id==nodeKD_id,netstate)
    if index_KD!=0
        nodeKD = netstate[index_KD]
        init_node!(nodeKD,nodeKO.def_state)
        for k in keys(nodeKD.F)
            nodeKD.F[k] = nodeKD.def_state
        end
    else
        error("Node id invalid: $nodeKD_id does not exist in the network")
    end
end

function constitutive_state!(netstate::Vector{BNode},const_node_id::AbstractString,const_state::Integer)
    const_node = filter(x->x.id==const_node_id,netstate)
    if(!isempty(const_node))
        const_node = const_node[1]
        [const_node.F[k] = const_state for k in keys(const_node.F)]
        init_node!(const_node,const_state)
        #const_node.state_range = [const_state]
    else
        error("Node id invalid: $nodeCA_id does not exist in the network")
    end
end


text:
"""
    MalmquistDEAModel
An data structure representing a malmquist productivity DEA model.
"""
struct MalmquistDEAModel <: AbstractProductivityDEAModel
    n::Int64
    m::Int64
    s::Int64
    periods::Int64
    orient::Symbol
    rts::Symbol
    refperiod::Symbol
    dmunames::Union{Vector{String},Nothing}
    Prod::Matrix
    EC::Matrix
    TC::Matrix
end

"""
    malmquist(X, Y)
Compute the Malmquist productivity index using data envelopment analysis for inputs X and outputs Y.

# Optional Arguments
- `orient=:Input`: chooses between input oriented radial model `:Input` or output oriented radial model `:Output`.
- `refperiod=:Geomean`: chooses reference period for technological change: `:Base`, `:Comparison` or `:Geomean`.
- `rts=:CRS`: chooses constant returns to scale. For variable returns to scale choose `:VRS`.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = Array{Float64,3}(undef, 5, 1, 2);

julia> X[:, :, 1] = [2; 3; 5; 4; 4];

julia> X[:, :, 2] = [1; 2; 4; 3; 4];

julia> Y = Array{Float64,3}(undef, 5, 1, 2);

julia> Y[:, :, 1] = [1; 4; 6; 3; 5];

julia> Y[:, :, 2] = [1; 4; 6; 3; 3];

julia> malmquist(X, Y)
Mamlmquist DEA Model
DMUs = 5; Inputs = 1; Outputs = 1; Time periods = 2
Orientation = Input; Returns to Scale = CRS
Referene period = Geomean
─────────────────────────
         M        EC   TC
─────────────────────────
1  2.0      1.33333   1.5
2  1.5      1.0       1.5
3  1.25     0.833333  1.5
4  1.33333  0.888889  1.5
5  0.6      0.4       1.5
─────────────────────────
M  = Malmquist Productivity Index
EC = Efficiency Change
TC = Technological Change
```
"""
function malmquist(X::Array{Float64,3}, Y::Array{Float64,3};
    orient::Symbol = :Input, rts::Symbol = :CRS, refperiod::Symbol = :Geomean,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::MalmquistDEAModel

    # Check parameters
    nx, m, Tx = size(X)
    ny, s, Ty = size(Y)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if Tx != Ty
        throw(DimensionMismatch("number of time periods in X and Y ($Tx, $Ty) are not equal"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(GLPK.Optimizer)
    end

    # Compute efficiency for each DMU
    n = nx
    T = Tx

    Prod  = zeros(n, T - 1)
    EC    = zeros(n, T - 1)
    TC    = zeros(n, T - 1)

    # For each time period
    for t = 1:T - 1

        # Compute efficiency in t
        efft = efficiency(dea(X[:, :, t], Y[:, :, t], orient = orient, rts = rts, slack = false, optimizer = optimizer))

        # Compute efficiency in t + 1
        efft1 = efficiency(dea(X[:, :, t + 1], Y[:, :, t + 1], orient = orient, rts = rts, slack = false, optimizer = optimizer))

        # Compute efficiency in t +1, with reference t
        efft1_reft = efficiency(dea(X[:, :, t + 1], Y[:, :, t + 1], orient = orient, rts = rts, slack = false,
                                    Xref = X[:, :, t], Yref = Y[:, :, t],
                                    optimizer = optimizer))

        # Compute efficiency in t, with reference t + 1
        efft_reft1 = efficiency(dea(X[:, :, t], Y[:, :, t], orient = orient, rts = rts, slack = false,
                                    Xref = X[:, :, t + 1], Yref = Y[:, :, t + 1],
                                    optimizer = optimizer))

        # Inver if output oriented
        if orient == :Output
            efft       .= 1 ./ efft
            efft1      .= 1 ./ efft1
            efft1_reft .= 1 ./ efft1_reft
            efft_reft1 .= 1 ./ efft_reft1
        end

        # Technical Efficiency Change
        EC[:, t] .= efft1 ./ efft

        # Technological Change
        if refperiod == :Base
            TC[:, t] .= efft1_reft ./ efft1
        elseif refperiod == :Comparison
            TC[:, t] .= efft ./ efft_reft1
        elseif refperiod == :Geomean
            TC[:, t] .= sqrt.( (efft ./ efft_reft1) .* (efft1_reft ./ efft1) )
        else
            throw(ArgumentError("`refperiod` must be :Base, :Comparison or :Geomean"));
        end

        # Mamlmquist Index
        Prod[:, t] .= EC[:, t] .* TC[:, t]

    end

    return MalmquistDEAModel(n, m, s, T, orient, rts, refperiod, names, Prod, EC, TC)

end

function Base.show(io::IO, x::MalmquistDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    dmunames = names(x)

    periods = nperiods(x)
    Prod = prodchange(x, :Prod)
    EC   = prodchange(x, :EC)
    TC   = prodchange(x, :TC)

    if !compact
        print(io, "Mamlmquist DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "; Time periods = ", periods)
        print(io, "\n")
        print(io, "Orientation = ", string(x.orient))
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        print(io, "Referene period = ", string(x.refperiod))
        print(io, "\n")
        if periods == 2
            show(io, CoefTable(hcat(Prod, EC, TC), ["M", "EC", "TC"], dmunames))
        end
        if periods > 2
            show(io, CoefTable(hcat(Prod, EC, TC),
                               vcat(["M$i" for i in 1:2], ["EC$i" for i in 1:2], ["TC$i" for i in 1:2 ]),
                               dmunames))
        end
        print(io, "\n")
        print(io, "M  = Malmquist Productivity Index \n")
        print(io, "EC = Efficiency Change \n")
        print(io, "TC = Technological Change")
    end

end