text:
# This file contains tests for the time series generator

# A module which contains the code to
# perform tests on .Gen module from Scats
module TestGen

using Scats.Internal.Prec      # Precision module from Scats
using Scats: Internal.Gen, API # API and .Gen module from Scats
using Test                     # A package to perform tests

# Load the general utilities
include("utilities.jl")

# Print the header
println("\e[1;32mRUNNING\e[0m: gen_test.jl")

"""
Read generator parameters from a file
"""
macro gen_read(file)
    return esc(:(s.Gen.read!($file)))
end

"""
Generate an example of a file containing generator parameters
"""
macro gen_example(file)
    return esc(:(s.Gen.example($file)))
end

"""
Reset the values of the generator parameters
"""
macro gen_reset()
    return esc(:(s.Gen.reset!()))
end

"""
Test the values of generator parameters
"""
macro test_values()
    return esc(quote
        @test s.Gen.N == 230
        @test s.Gen.Δt == 1
        @test s.Gen.q == 0.01
        @test s.Gen.α == 0.1
        @test s.Gen.β == 0.05
        @test s.Gen.r == 1
        @test s.Gen.A == [1]
        @test s.Gen.ν == [0.1]
        @test s.Gen.ϕ == [0]
        @test s.Gen.γ == 0.5
    end)
end

"""
Create a temporary file, write exemplary
generator parameters in it and read them
"""
macro read_example_params()
    return esc(quote
        @file
        @gen_example file
        @gen_read file
    end)
end

"""
Test the end of file exception on a file
"""
macro test_eof_exception()
    return esc(quote
        try
            @gen_read file
        catch e
            @test e isa Gen.ScatsGenEOF
            @test sprint(showerror, e) == retrieve_messages(
                Gen.ScatsGenEOF, file,
            )
        end
    end)
end

@construct_exceptions "Gen.ScatsGen"

# Test exceptions related to file status
@testset "Check file status" begin

    @instantiate

    # Test when a wrong path has been specified
    try
        @gen_read "Wrong file path!"
    catch e
        @test e isa Gen.ScatsGenNotAFile
        @test sprint(showerror, e) == retrieve_messages(
            Gen.ScatsGenNotAFile, "Wrong file path!"
        )
    end

    @file
    @test_eof_exception

    # Test the end of file exception on files with different number of lines
    for i in 1:28

        # Generate valid data on input lines
        if !(i in range(2, 29, step = 3))
            println(io, "Line #", i)
        elseif i == 2 || i == 17
            println(io, IT(1))
        else
            println(io, RT(1))
        end

        # Update the file immediately
        flush(io)

        @test_eof_exception

    end

end

# Test different ways to read parameters
@testset "Read `good` parameters" begin

    @instantiate

    @read_example_params

    @test_values
    @gen_reset

    # Read the generator parameters (another way)
    s.Gen.read!(file)

    @test_values
    @gen_reset

    # Read the generator parameters (another way)
    s.Gen(file)

    @test_values

end

# println(a)

# Test the creation of an example of generator parameters
@testset "Check creation of an example" begin

    @instantiate
    @read_example_params
    @test_values

    # Create a temporary directory
    dir = mktempdir()

    # Test the exception being thrown when a directory has been passed
    try
        @gen_example dir
    catch e
        @test e isa Gen.ScatsGenIsADir
        @test sprint(showerror, e) == retrieve_messages(
            Gen.ScatsGenIsADir, dir,
        )
    end

    # Write generator parameters (another way)
    s.gen_example(file)

    @gen_read file
    @test_values

    # Test the exception being thrown when a directory has been passed (another way)
    try
        s.gen_example(dir)
    catch e
        @test e isa Gen.ScatsGenIsADir
        @test sprint(showerror, e) == retrieve_messages(
            Gen.ScatsGenIsADir, dir,
        )
    end

end

# Test all exceptions related to invalid input
@testset "Read `bad` parameters" begin

    @instantiate

    # Construct an ordered list of exceptions
    exceptions = construct_exceptions(
        "WR", "N", "Δt", "q", "α", "β", "r", "A", "ν", "ϕ", "γ",
    )

    # Get an ordered list of expected messages
    messages = retrieve_messages(exceptions, "gen")

    # Create a temporary file to contain valid parameters
    good_file, _ = mktemp()

    @gen_example good_file

    # Corrupt a file on a specific line
    @inline function break_a_line!(ln::Int)

        @file

        open(good_file) do good_io
            k = 0
            # Print lines of the `good` file into the `bad` file
            for line in eachline(good_io)

                k += 1

                # Corrupt a line
                if k == ln
                    line = "Hello there!\n"
                end

                # Print a line
                println(io, line)

            end
        end

        # Close the stream
        close(io)

        # Save the temporary file as `gen`
        cp(file, "gen", force = true)

    end

    # Test all cases of corrupted input
    for i in 1:10

        # Calculate the number of a line to corrupt
        n = 2 + (i - 1) * 3

        # Corrupt a specific line
        break_a_line!(n)

        # Test an exception
        try
            @gen_read "gen"
        catch e
            @test e isa exceptions[i]
            @test sprint(showerror, e) == messages[i]
        end

    end

    # Remove the `gen` file
    isfile("gen") && rm("gen")

end

# Test generation of time series
@testset "Check generation" begin

    @instantiate
    @read_example_params

    # Generate time series
    s.gen!()

    # Test scalars
    @test s.Input.N != 0
    @test s.Input.Δt != 0
    @test s.Input.q != 0

    # Test arrays
    @test s.Input.t[1] == 0
    for t in s.Input.t[2:s.Input.N]
        @test t != 0
    end
    for x in s.Input.x[1:s.Input.N]
        @test x != 0
    end

end

# Test values resetting
@testset "Check resetting" begin

    @instantiate
    @read_example_params

    @gen_reset

    # Test values
    @test s.Gen.N == 0
    @test s.Gen.Δt == 0
    @test s.Gen.q == 0
    @test s.Gen.α == 0
    @test s.Gen.β == 0
    @test s.Gen.r == 0
    @test s.Gen.A == []
    @test s.Gen.ν == []
    @test s.Gen.ϕ == []
    @test s.Gen.γ == 0

end

end


text:
@testset "Central difference method" begin
    # Ref. Example 9.1 pp. 773-774
    prob, U = example_9_1_Bathe()
    alg = CentralDifference(Δt=0.28)

    # solution
    sol = solve(prob, alg, NSTEPS=12) |> displacements

    @test abs(sol[13][1] - 1.02) < 5e-3
    @test abs(sol[13][2] - 2.60) < 5e-3
end

@testset "Houbolt method" begin
    # Ref. Example 9.2 pp. 776-777
    prob, U = example_9_1_Bathe()
    alg = Houbolt(Δt=0.28)

    # solution
    sol = solve(prob, alg, NSTEPS=12) |> displacements

    @test abs(sol[13][1] - 1.72) < 5e-3
    @test abs(sol[13][2] - 2.28) < 5e-3
end

@testset "Newmark method" begin
    # Ref. Example 9.3 pp. 778-779
    prob, U = example_9_1_Bathe()
    alg = Newmark(Δt=0.28, α=0.25, δ=0.5)

    # solution
    sol = solve(prob, alg, NSTEPS=12) |> displacements

    @test abs(sol[13][1] - 1.40) < 5e-3
    @test abs(sol[13][2] - 2.31) < 5e-3
end

@testset "Bathe method" begin
    # Ref. Example 9.4 pp. 781-782
    prob, U = example_9_1_Bathe()
    alg = Bathe(Δt=0.28)

    # solution
    sol = solve(prob, alg, NSTEPS=12) |> displacements

    @test abs(sol[13][1] - 1.28) < 5e-3
    @test abs(sol[13][2] - 2.40) < 5e-3
end

@testset "Bathe method with forcing term" begin
    ivp, NSTEPS, Δt = harmonic_oscillator_forced()

    alg = Bathe(Δt = Δt)
    sol = solve(ivp, alg, NSTEPS=NSTEPS)

    # TODO add analytic solution for test
end

@testset "BackwardEuler method" begin
    prob, U = heat_transfer()

    testΔt = 0.0001
    alg = BackwardEuler(Δt=testΔt)

    sol = solve(prob, alg, NSTEPS=12) |> displacements

    # relative error at node 1 verification
    analyticVal = U(testΔt)
    @test (abs(sol[13][1] - analyticVal) / abs(analyticVal)) < 5e-3
end

@testset "Nonlinear Newmark method: Pendulum" begin

    # load model
    prob = pendulum_nonlinear()

    # solve using Trapezoidal scheme
    α = 1/4
    δ = 1/2
    Δt = 0.01
    alg = Newmark(Δt=Δt, α=α, δ=δ)

    sol = solve(prob, alg, T=3.0)

    # FIXME review tests
    @test sol.U[end] ≈ [1.9113595927872666, 1.410419996390988]
    @test sol.U′[end] ≈ [1.0007721533056833, 3.2500475651804632]
    @test sol.U′′[end] ≈ [-8.286971439792033, -7.2537920159991955]
end

@testset "Nonlinear Newmark method: Von Mises Truss" begin

    # load model
    prob = von_mises_truss()

    # algorithm parameters
    α = 1/4
    δ = 1/2
    Δt = 0.0025
    alg = Newmark(Δt=Δt, α=α, δ=δ)

    # solve using Trapezoidal scheme
    sol = solve(prob, alg, T=2.0)

    # FIXME update tests
    sol.U[end] ≈ [0.001888137400895257, -0.17208311613153432]
    sol.U′[end] ≈ [-0.07080458802627602, 0.14912664021827504]
    sol.U′′[end] ≈ [6.0851727632872485, -6.18000018114007]
end

@testset "Nonlinear Central difference method: Von Mises Truss" begin

    # load model
    prob = von_mises_truss()

    # algorithm parameters
    Δt = 0.000025
    alg = CentralDifference(Δt=Δt)

    # solve
    sol = solve(prob, alg, T=2.0)

    # FIXME update tests
end

ext: jl


text:
module ThreeDP3

greet() = print("Hello World!")

using Gen

import GLRenderer as GL
import GenDirectionalStats
import MiniGSG as S
import Rotations as R
import PoseComposition as PC
import PoseComposition: Pose
import StaticArrays:SVector, @SVector, StaticVector, SMatrix
import Colors
import PyCall
import LightGraphs as LG
import MetaGraphs as MG

include("utils/bbox.jl")
include("utils/clustering.jl")
include("utils/data_loader.jl")
include("utils/icp.jl")
include("utils/occlusions.jl")
include("utils/pose.jl")
include("utils/pose_utils.jl")
include("utils/scale_camera.jl")
include("utils/utils.jl")

include("distributions/distributions.jl")

include("contact/contact.jl")
include("graphs/graphs.jl")
include("shape/shape.jl")

include("model/sg_model.jl")
include("model/sg_model_utils.jl")

include("inference/inference.jl")

end # module


text:
function part1(input)
    board, position = parse_input(input, part = 1)
    return play(board, position, UInt8.(max.(board.area, 0)))
end

function part2(input)
    board, position = parse_input(input, part = 2)
    return play(board, position, UInt8.(max.(board.area, 0)))
end

struct Board
    area::Vector{Int}
    distances::Matrix{Int}
    moves::Vector{Vector{Vector{Int}}}
end

function parse_input(input; part)
    raw_board = collect.(readlines(input))
    if part == 2
        insert!(raw_board, 4, collect("  #D#C#B#A#"))
        insert!(raw_board, 5, collect("  #D#B#A#C#"))
    end

    n = 2 * part
    area = vcat(0, 0, -1, 0, -1, 0, -1, 0, -1, 0, 0,
                fill(1, n), fill(2, n), fill(3, n), fill(4, n))
    neighbors = [[2],
                 [1, 3],
                 [2, 4, 12],
                 [3, 5],
                 [4, 6, 12 + n],
                 [5, 7],
                 [6, 8, 12 + 2 * n],
                 [7, 9],
                 [8, 10, 12 + 3 * n],
                 [9, 11],
                 [10]]
    for i = 1:4
        for j = 1:n
            j == 1 && push!(neighbors, [1 + 2 * i, 13 + n * (i - 1)])
            1 < j < n && push!(neighbors, [10 + n * (i - 1) + j, 12 + n * (i - 1) + j])
            j == n && push!(neighbors, [10 + n * (i - 1) + j])
        end
    end

    N = length(area)
    distances = zeros(Int, N, N)
    moves = Vector{Vector{Int}}[]
    for i = 1:length(area)
        if area[i] < 0
            push!(moves, Int[])
            continue
        end
        visited = falses(N)
        dist = zeros(N)
        push!(moves, Vector{Int}[])
        dfs!(area, neighbors, visited, Int[], dist, moves[i], i, 0)
        if area[i] == 0
            filter!(x -> area[last(x)] > 0, moves[i])
        else
            filter!(x -> area[last(x)] != area[first(x)], moves[i])
        end
        distances[i, :] .= dist
    end

    position = vcat(zeros(UInt8, 11))
    for i = 1:4
        for j = 1:n
            push!(position, 1 + raw_board[2 + j][2 + 2 * i] - UInt8('A'))
        end
    end

    return Board(area, distances, moves), position
end

function dfs!(area, neighbors, visited, path, dist, moves, vertex, depth)
    visited[vertex] = true
    dist[vertex] = depth
    if depth > 0
        push!(path, vertex)
        if area[vertex] >= 0
            push!(moves, copy(path))
        end
    end
    for neighbor in neighbors[vertex]
        if !visited[neighbor]
            dfs!(area, neighbors, visited, path, dist, moves, neighbor, depth + 1)
        end
    end
    if depth > 0
        pop!(path)
    end
end

mutable struct PriorityQueue
    queue::Dict{Int, Vector{Vector{UInt8}}}
    min_cost::Int
end

PriorityQueue() = PriorityQueue(Dict{Int, Vector{Vector{UInt8}}}(), 0)

Base.isempty(queue::PriorityQueue) = isempty(queue.queue)

function Base.push!(queue::PriorityQueue, item, cost)
    @assert cost >= queue.min_cost
    push!(get!(() -> Vector{UInt8}[], queue.queue, cost), item)
end

function Base.pop!(queue::PriorityQueue)
    while !haskey(queue.queue, queue.min_cost) || isempty(queue.queue[queue.min_cost])
        if haskey(queue.queue, queue.min_cost)
            delete!(queue.queue, queue.min_cost)
        end
        queue.min_cost += 1
    end
    return pop!(queue.queue[queue.min_cost])
end

function play(board, position, goal)
    costs = Dict(copy(position) => 0)
    queue = PriorityQueue()
    push!(queue, position, 0)
    best_cost = typemax(Int)
    while !isempty(queue)
        position = pop!(queue)
        if costs[position] < queue.min_cost
            continue
        end
        if position == goal
            return costs[position]
        end
        for i in 1:length(position)
            color = position[i]
            color == 0 && continue
            for move in board.moves[i]
                destination = last(move)
                board.area[destination] ∉ (0, color) && continue
                any(position[j] != 0 for j in move) && continue
                if board.area[destination] > 0
                    if !valid_destination(board, position, destination, color)
                        continue
                    end
                end
                new_position = copy(position)
                new_position[i] = 0
                new_position[destination] = color
                cost = costs[position] + move_cost(color, board.distances[i, destination])
                cost > best_cost && continue
                if new_position == goal
                    best_cost = min(cost, best_cost)
                end
                if !haskey(costs, new_position)
                    costs[new_position] = cost
                    if cost < best_cost || new_position == goal
                        push!(queue, new_position, cost)
                    end
                else
                    if cost < costs[new_position]
                        costs[new_position] = cost
                        push!(queue, new_position, cost)
                    end
                end
            end
        end
    end
end

move_cost(color, distance) = 10^(color - 1) * distance

function valid_destination(board, position, destination, color)
    while true
        destination += 1
        if destination > length(board.area) || board.area[destination] != color
            return true
        end
        if position[destination] != color
            return false
        end
    end
end


text:
import LibGEOS
import GeoJSON
import Base.convert
using DataFrames
using Distributions
using Printf
import CSV
import JSON
import Proj4
import Proj4: transform

import GeoInterface
using GeoInterface: features, coordinates, geometry, properties

const SchDistr = Int

const SALE_PRICE = Symbol("SALE_PRICE")
const SQFT = Symbol("GROSS_SQUARE_FEET")
const BUILDING_CLASS_AT_TIME_OF_SALE = Symbol("BUILDING_CLASS_AT_TIME_OF_SALE")
const BUILDING_CLASS_CATEGORY = Symbol("BUILDING_CLASS_CATEGORY")
const TAX_CLASS_AT_TIME_OF_SALE = Symbol("TAX_CLASS_AT_TIME_OF_SALE")

function transform(src::Proj4.Projection, dest::Proj4.Projection, multicoords::AbstractVector{V} where V<:AbstractVector)
    return transform.(Ref(src), Ref(dest), multicoords)
end
function transform(src, dist, geom::G) where G<:GeoInterface.AbstractGeometry
    G(transform(src, dist, coordinates(geom)))
end
function transform_epsg(coords; epsg_from::Int, epsg_to::Int)
    proj_from = Proj4.Projection(Proj4.epsg[epsg_from])
    proj_to = Proj4.Projection(Proj4.epsg[epsg_to])
    return transform(proj_from, proj_to, coords)
end

function read_distr_shapes(; filedir="nysd_16c", filename="nysd.json", data_dir="NYC_data", epsg_from::Int=4326, epsg_to::Int=2263)
    nysd_json = GeoJSON.read(read(joinpath(data_dir, filedir, filename)))
    schdistr_shape_dict = Dict{SchDistr, GeoRDD.RegionType}()
    for feature in features(nysd_json)
        schdistr = convert(SchDistr, properties(feature)["SchoolDist"])
        projected_geom = transform_epsg(geometry(feature); epsg_from=epsg_from, epsg_to=epsg_to)
        shape = convert(GeoRDD.RegionType, projected_geom)
        schdistr_shape_dict[schdistr] = shape
    end
    return schdistr_shape_dict
end

function read_processed_sales(; filename="NYC_sales.csv", filedir="processed", data_dir="NYC_data")
    NYC_sales=CSV.read(joinpath(data_dir, filedir, filename),
                       types=Dict("TAX CLASS AT PRESENT" => Union{Missings.Missing, String},
                                  "TAX CLASS AT TIME OF SALE" => Union{Missings.Missing, String}),
                       copycols=true
                      )
    nyc_schdistrs = CategoricalVector(NYC_sales[!, :SchDistr])
    # schd_strings = [ismissing(sd)?missing:dec(sd,2) for sd in  nyc_schdistrs]
    # str_schdistrs = CategoricalVector(schd_strings)
    NYC_sales[!, :SchDistr] = nyc_schdistrs
    # categorical variables
    categorical!(NYC_sales, [
        :BOROUGH,
        :BUILDING_CLASS_CATEGORY,
        :BUILDING_CLASS_AT_TIME_OF_SALE,
        :ZIP_CODE,
        :TAX_CLASS_AT_PRESENT,
        :TAX_CLASS_AT_TIME_OF_SALE,
        :NEIGHBORHOOD,
        ])
    sort!(NYC_sales, :SchDistr)
    return NYC_sales
end

const DWELLINGS_DICT = Dict(
    "01  ONE FAMILY DWELLINGS"=>true,
    "02  TWO FAMILY DWELLINGS"=>true,
    "03  THREE FAMILY DWELLINGS"=>true,
    "04  TAX CLASS 1 CONDOS"=>false,
    "05  TAX CLASS 1 VACANT LAND"=>false,
    "06  TAX CLASS 1 - OTHER"=>false,
    "07  RENTALS - WALKUP APARTMENTS"=>false,
    "08  RENTALS - ELEVATOR APARTMENTS"=>false,
    "09  COOPS - WALKUP APARTMENTS"=>false,
    "10  COOPS - ELEVATOR APARTMENTS"=>false,
    "11  SPECIAL CONDO BILLING LOTS"=>false,
    "11A CONDO-RENTALS"=>false,
    "12  CONDOS - WALKUP APARTMENTS"=>false, # why are these duplicated?
    "13  CONDOS - ELEVATOR APARTMENTS"=>false,
    "14  RENTALS - 4-10 UNIT"=>false,
    "15  CONDOS - 2-10 UNIT RESIDENTIAL"=>false,
    "16  CONDOS - 2-10 UNIT WITH COMMERCIAL UNIT"=>false,
    "17  CONDO COOPS"=>false,
    "18  TAX CLASS 3 - UNTILITY PROPERTIES"=>false,
    "21  OFFICE BUILDINGS"=>false,
    "22  STORE BUILDINGS"=>false,
    "23  LOFT BUILDINGS"=>false,
    "24  TAX CLASS 4 - UTILITY BUREAU PROPERTIES"=>false,
    "25  LUXURY HOTELS"=>false,
    "26  OTHER HOTELS"=>false,
    "27  FACTORIES"=>false,
    "28  COMMERCIAL CONDOS"=>false,
    "29  COMMERCIAL GARAGES"=>false,
    "30  WAREHOUSES"=>false,
    "31  COMMERCIAL VACANT LAND"=>false,
    "32  HOSPITAL AND HEALTH FACILITIES"=>false,
    "33  EDUCATIONAL FACILITIES"=>false,
    "34  THEATRES"=>false,
    "35  INDOOR PUBLIC AND CULTURAL FACILITIES"=>false,
    "36  OUTDOOR RECREATIONAL FACILITIES"=>false,
    "37  RELIGIOUS FACILITIES"=>false,
    "38  ASYLUMS AND HOMES"=>false,
    "39  TRANSPORTATION FACILITIES"=>false,
    "40  SELECTED GOVERNMENTAL FACILITIES"=>false,
    "41  TAX CLASS 4 - OTHER"=>false,
    "42  CONDO CULTURAL/MEDICAL/EDUCATIONAL/ETC"=>false,
    "43  CONDO OFFICE BUILDINGS"=>false,
    "44  CONDO PARKING"=>false,
    "45  CONDO HOTELS"=>false,
    "46  CONDO STORE BUILDINGS"=>false,
    "47  CONDO NON-BUSINESS STORAGE"=>false,
    "48  CONDO TERRACES/GARDENS/CABANAS"=>false,
    "49  CONDO WAREHOUSES/FACTORY/INDUS"=>false,
    )
const RESIDENTIAL_DICT = Dict(
    "01  ONE FAMILY DWELLINGS"=>true,
    "02  TWO FAMILY DWELLINGS"=>true,
    "03  THREE FAMILY DWELLINGS"=>true,
    "04  TAX CLASS 1 CONDOS"=>true,
    "05  TAX CLASS 1 VACANT LAND"=>false,
    "06  TAX CLASS 1 - OTHER"=>false,
    "07  RENTALS - WALKUP APARTMENTS"=>false,
    "08  RENTALS - ELEVATOR APARTMENTS"=>false,
    "09  COOPS - WALKUP APARTMENTS"=>true,
    "10  COOPS - ELEVATOR APARTMENTS"=>true,
    "11  SPECIAL CONDO BILLING LOTS"=>false,
    "11A CONDO-RENTALS"=>false,
    "12  CONDOS - WALKUP APARTMENTS"=>true, # why are these duplicated?
    "13  CONDOS - ELEVATOR APARTMENTS"=>true,
    "14  RENTALS - 4-10 UNIT"=>false,
    "15  CONDOS - 2-10 UNIT RESIDENTIAL"=>true,
    "16  CONDOS - 2-10 UNIT WITH COMMERCIAL UNIT"=>false,
    "17  CONDO COOPS"=>true,
    "18  TAX CLASS 3 - UNTILITY PROPERTIES"=>false,
    "21  OFFICE BUILDINGS"=>false,
    "22  STORE BUILDINGS"=>false,
    "23  LOFT BUILDINGS"=>false,
    "24  TAX CLASS 4 - UTILITY BUREAU PROPERTIES"=>false,
    "25  LUXURY HOTELS"=>false,
    "26  OTHER HOTELS"=>false,
    "27  FACTORIES"=>false,
    "28  COMMERCIAL CONDOS"=>false,
    "29  COMMERCIAL GARAGES"=>false,
    "30  WAREHOUSES"=>false,
    "31  COMMERCIAL VACANT LAND"=>false,
    "32  HOSPITAL AND HEALTH FACILITIES"=>false,
    "33  EDUCATIONAL FACILITIES"=>false,
    "34  THEATRES"=>false,
    "35  INDOOR PUBLIC AND CULTURAL FACILITIES"=>false,
    "36  OUTDOOR RECREATIONAL FACILITIES"=>false,
    "37  RELIGIOUS FACILITIES"=>false,
    "38  ASYLUMS AND HOMES"=>false,
    "39  TRANSPORTATION FACILITIES"=>false,
    "40  SELECTED GOVERNMENTAL FACILITIES"=>false,
    "41  TAX CLASS 4 - OTHER"=>false,
    "42  CONDO CULTURAL/MEDICAL/EDUCATIONAL/ETC"=>false,
    "43  CONDO OFFICE BUILDINGS"=>false,
    "44  CONDO PARKING"=>false,
    "45  CONDO HOTELS"=>false,
    "46  CONDO STORE BUILDINGS"=>false,
    "47  CONDO NON-BUSINESS STORAGE"=>false,
    "48  CONDO TERRACES/GARDENS/CABANAS"=>false,
    "49  CONDO WAREHOUSES/FACTORY/INDUS"=>false,
    )

function filter_sales(NYC_sales::DataFrame)

    logSalePricePerSQFT = log.(NYC_sales[!,SALE_PRICE]) .- log.(NYC_sales[!,SQFT])
    believable = zeros(Bool, size(NYC_sales,1))
    removed_reason = Dict(
        :residential => 0,
        :noprice => 0,
        :nosqft => 0,
        :geocode => 0,
        :toosmall => 0,
        :outlier => 0,
       )
    for i in 1:size(NYC_sales,1)
        if ismissing(NYC_sales[i,BUILDING_CLASS_AT_TIME_OF_SALE])
            # remove sales with missing covariates
            removed_reason[:residential] += 1
        elseif ismissing(NYC_sales[i,BUILDING_CLASS_CATEGORY])
            # remove sales with missing covariates
            removed_reason[:residential] += 1
        elseif !DWELLINGS_DICT[NYC_sales[i,BUILDING_CLASS_CATEGORY]]
            # remove non-residential properties
            # in fact, remove things that aren't dwellings (houses?)
            removed_reason[:residential] += 1
        elseif ismissing(NYC_sales[i,SALE_PRICE])
            # remove sales without a sale price
            removed_reason[:noprice] += 1
        elseif ismissing(NYC_sales[i,SQFT])
            # remove sales with missing GROSS SQUARE FEET information
            removed_reason[:nosqft] += 1
        elseif ismissing(NYC_sales[i,:XCoord])
            # remove properties with failed geocoding
            removed_reason[:geocode] += 1
        elseif ismissing(NYC_sales[i,:YCoord])
            # remove properties with failed geocoding
            removed_reason[:geocode] += 1
        elseif ismissing(NYC_sales[i,:SchDistr])
            # remove properties without a school district
            removed_reason[:geocode] += 1
        # elseif ismissing(NYC_sales[i,TAX_CLASS_AT_TIME_OF_SALE])
            # # remove sales with missing covariates
            # removed_reason[7] += 1
        elseif NYC_sales[i,SQFT] < 100.0
            # remove properties smaller to 100sqft (seems unlikely to be real)
            removed_reason[:toosmall] += 1
        elseif logSalePricePerSQFT[i] < 3.0
            # that's too cheap (remove outliers)
            removed_reason[:outlier] += 1
        elseif logSalePricePerSQFT[i] > 8.0
            # that's too expensive (remove outliers)
            removed_reason[:outlier] += 1
        else
            # otherwise keep
            believable[i] = true
        end
    end

    filtered = copy(NYC_sales[believable,:])
    filtered[!,:logSalePricePerSQFT] = convert.(Float64, logSalePricePerSQFT[believable])
    for col in (:XCoord, :YCoord)
        filtered[!,col] = convert.(Float64, filtered[!,col])
    end
    return Dict(
        :filtered => filtered,
        :believable => believable,
        :removed_reason => removed_reason,
        )
end

# function sales_dicts(NYC_sales::DataFrame)
    # schdistrs = sort(NYC_sales[:SchDistr].pool.levels)
    # schdistr_indices = Dict{SchDistr,Vector{Int}}()
    # schdistrs_col = NYC_sales[:SchDistr]
    # for distr in schdistrs
        # indices = find(schdistrs_col.refs .== find(schdistrs_col.pool.index .== distr)[1])
        # nobsv_schdistr = length(indices)
        # @printf("District %d has %d sales\n", distr, nobsv_schdistr)
        # schdistr_indices[distr] = indices
    # end
    # Y_dict=Dict{SchDistr, Vector{Float64}}()
    # X_dict=Dict{SchDistr, Array{Float64,2}}()
    # for distr in schdistrs
        # Y_dict[distr] = NYC_sales[schdistr_indices[distr], :logSalePricePerSQFT]
        # X_dict[distr] = NYC_sales[schdistr_indices[distr],[:XCoord, :YCoord]]
    # end
    # return Dict(
            # :schdistr_indices => schdistr_indices,
            # :schdistrs => schdistrs,
            # :X_dict => X_dict,
            # :Y_dict => Y_dict
            # )
# end

function table_τpair(τpost_pairs::Dict{Tuple{NYC.SchDistr,NYC.SchDistr},Normal}, include_distr::Vector{Int})
    τ_pair_nested = Dict{NYC.SchDistr, Dict{NYC.SchDistr, Normal}}()
    for (distrA, distrB) in keys(τpost_pairs)
        τpost = τpost_pairs[distrA, distrB]
        if distrA ∉ include_distr
            continue
        end
        if distrB ∉ include_distr
            continue
        end
        if distrA ∉ keys(τ_pair_nested)
            τ_pair_nested[distrA] = Dict{NYC.SchDistr, Normal}()
        end
        τ_pair_nested[distrA][distrB] = τpost
    end
    ;

    for distrA in sort(collect(keys(τ_pair_nested)))
        if isempty(τ_pair_nested[distrA])
            continue
        end
        @printf("\\( \\mathbf{%2d} \\)", distrA)
        for distrB in sort(collect(keys(τ_pair_nested[distrA])))
            τ = τ_pair_nested[distrA][distrB]
            @printf("& \\( \\mathbf{%2d:}~%+.2f \\pm %.2f \\)", distrB, -mean(τ), std(τ))
        end
        print("\\\\ \n")
    end
end


text:
# ---
# title: 752. Open the Lock
# id: problem752
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Breadth-first Search
# link: <https://leetcode.com/problems/open-the-lock/description/>
# hidden: true
# ---
#
# You have a lock in front of you with 4 circular wheels. Each wheel has 10
# slots: `'0', '1', '2', '3', '4', '5', '6', '7', '8', '9'`. The wheels can
# rotate freely and wrap around: for example we can turn `'9'` to be `'0'`, or
# `'0'` to be `'9'`. Each move consists of turning one wheel one slot.
#
# The lock initially starts at `'0000'`, a string representing the state of the
# 4 wheels.
#
# You are given a list of `deadends` dead ends, meaning if the lock displays any
# of these codes, the wheels of the lock will stop turning and you will be
# unable to open it.
#
# Given a `target` representing the value of the wheels that will unlock the
# lock, return the minimum total number of turns required to open the lock, or
# -1 if it is impossible.
#
#
#
# **Example 1:**
#
#
#
#     Input: deadends = ["0201","0101","0102","1212","2002"], target = "0202"
#     Output: 6
#     Explanation:
#     A sequence of valid moves would be "0000" -> "1000" -> "1100" -> "1200" -> "1201" -> "1202" -> "0202".
#     Note that a sequence like "0000" -> "0001" -> "0002" -> "0102" -> "0202" would be invalid,
#     because the wheels of the lock become stuck after the display becomes the dead end "0102".
#
#
# **Example 2:**
#
#
#
#     Input: deadends = ["8888"], target = "0009"
#     Output: 1
#     Explanation:
#     We can turn the last wheel in reverse to move from "0000" -> "0009".
#
#
# **Example 3:**
#
#
#
#     Input: deadends = ["8887","8889","8878","8898","8788","8988","7888","9888"], target = "8888"
#     Output: -1
#     Explanation:
#     We can't reach the target without getting stuck.
#
#
# **Example 4:**
#
#
#
#     Input: deadends = ["0000"], target = "8888"
#     Output: -1
#
#
#
#
# **Constraints:**
#
#   * `1 <= deadends.length <= 500`
#   * `deadends[i].length == 4`
#   * `target.length == 4`
#   * target **will not be** in the list `deadends`.
#   * `target` and `deadends[i]` consist of digits only.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
POMDPs.solve(solver::IPFTSolver, mdp::Union{POMDP,MDP}) = IPFTPlanner(solver, mdp)

"""
Delete existing decision tree.
"""
function clear_tree!(p::IPFTPlanner)
    p.tree = nothing
end

MCTS.init_Q(q::Array{Float64}, mdp::Union{MDP,POMDP}, s, a) = q

"""
Construct an IPFT tree and choose the best action.
"""
POMDPs.action(p::IPFTPlanner, s) = first(action_info(p, s))

"""
Construct an IPFT tree and choose the best action. Also output some information.
"""
function POMDPModelTools.action_info(p::IPFTPlanner, s; tree_in_info=false)
    local a::actiontype(p.mdp)
    info = Dict{Symbol, Any}()
    try
        if isterminal(p.mdp, s)
            error("""
                  IPFT cannot handle terminal states. action was called with
                  s = $s
                  """)
        end

        S = statetype(p.mdp)
        @assert s isa S "s:$(typeof(s)), S=$S"
        A = actiontype(p.mdp)
        if p.solver.keep_tree
            if p.tree == nothing
                tree = IPFTree{S,A}(p.solver.n_iterations)
                p.tree = tree
            else
                tree = p.tree
            end
            if haskey(tree.s_lookup, s)
                snode = tree.s_lookup[s]
            else
                snode = insert_state_node!(tree, s, true)
            end
        else
            tree = IPFTree{S,A}(p.solver.n_iterations)
            p.tree = tree
            snode = insert_state_node!(tree, s, p.solver.check_repeat_state)
        end

        nquery = 0
        start_us = CPUtime_us()
        for j = 1:p.solver.n_iterations
            nquery += 1
            @debug "Start simulation $j"
            simulate(p, snode, s, p.solver.depth)
            if CPUtime_us() - start_us >= p.solver.max_time * 1e6
                break
            end
        end
        info[:search_time_us] = CPUtime_us() - start_us
        info[:tree_queries] = nquery
        if p.solver.tree_in_info || tree_in_info
            info[:tree] = tree
        end

        sanode = select_best(p.solver.final_criterion, IPFTStateNode(tree,snode), p.rng)
        # XXX some publications say to choose action that has been visited the most
        a = tree.a_labels[sanode] # choose action with highest approximate value
    catch ex
        a = convert(actiontype(p.mdp), default_action(p.solver.default_action, p.mdp, s, ex))
        info[:exception] = ex
    end

    return a, info
end

"""
Return the reward for one iteration of IPFT.
"""
function POMDPModelTools.simulate(p::IPFTPlanner{P, S, A, SE, NA, RNG}, snode::Int, s::S, d::Int) where {P, S, A, SE, NA, RNG}
    sol = p.solver
    @debug "Run simulate on depth $(sol.depth - d)"
    tree = p.tree
    # s = tree.s_labels[snode]
    if d == 0
        return [0.0, 0.0]
    end
    # action progressive widening
    if sol.enable_action_pw
        if length(tree.children[snode]) <= sol.k_action*tree.total_n[snode]^sol.alpha_action # criterion for new action generation
            a = next_action(p.next_action, p.mdp, s, IPFTStateNode(tree, snode)) # action generation step
            if !sol.check_repeat_action || !haskey(tree.a_lookup, (snode, a))
                n0 = init_N(sol.init_N, p.mdp, s, a)
                insert_action_node!(tree, snode, a, n0,
                                    init_Q(sol.init_Q, p.mdp, s, a),
                                    sol.check_repeat_action
                                   )
                tree.total_n[snode] += n0
            end
        end
    elseif isempty(tree.children[snode])
        for a in actions(p.mdp, s)
            n0 = init_N(sol.init_N, p.mdp, s, a)
            insert_action_node!(tree, snode, a, n0,
                                init_Q(sol.init_Q, p.mdp, s, a),
                                false)
            tree.total_n[snode] += n0
        end
    end
    sanode = select_best(sol.criterion, IPFTStateNode(tree,snode), p.rng)
    a = tree.a_labels[sanode]

    # state progressive widening
    new_node = false
    if tree.n_a_children[sanode] <= sol.k_observation*tree.n[sanode]^sol.alpha_observation
        sp, r = generate_sr(p.mdp, s, a, p.rng)

        if sol.check_repeat_state && haskey(tree.s_lookup, sp)
            spnode = tree.s_lookup[sp]
        else
            spnode = insert_state_node!(tree, sp, sol.keep_tree || sol.check_repeat_state)
            new_node = true
        end

        push!(tree.transitions[sanode], (spnode, sp, r))

        if !sol.check_repeat_state
            tree.n_a_children[sanode] += 1
        elseif !((sanode,spnode) in tree.unique_transitions)
            push!(tree.unique_transitions, (sanode,spnode))
            tree.n_a_children[sanode] += 1
        end
    else
        spnode, sp, r = rand(p.rng, tree.transitions[sanode])
        if !sol.reuse_states
            sp, r = generate_sr(p.mdp, s, a, p.rng)
        end
    end

    if isterminal(p.mdp, sp, a)
        q = r
    elseif new_node
        q = r + discount(p.mdp)*estimate_value(p.solved_estimate, p.mdp, sp, d-1)
    else
        q = r + discount(p.mdp)*simulate(p, spnode, sp, d-1)
    end

    tree.n[sanode] += 1
    tree.total_n[snode] += 1

    tree.q[sanode] += (q - tree.q[sanode])/tree.n[sanode]

    return q
end


text:
using GraphBenchmarks
using Base.Test

include(joinpath(Pkg.dir("GraphBenchmarks", "src", "examples"), "lg.jl"))
#include(joinpath(Pkg.dir("GraphBenchmarks", "src", "examples"), "stinger.jl"))
# write your own tests here
const benchmarkspecs = (Graph500(),)#(StingerBenchmarks(), Graph500())
for benchmarkspec in benchmarkspecs
    info("Running $benchmarkspec benchmarks")
    #= #Uncomment to run the StingerWrapper tests
    info("Running SerialBFS for StingerGraph")
    @test isa(benchmark(benchmarkspec, StingerGraph, SerialBFS(), Kronecker(10, 16)), BenchmarkTools.Trial)
    info("Running LevelSynchronousBFS for StingerGraph")
    @test isa(benchmark(benchmarkspec, StingerGraph, LevelSynchronousBFS(), Kronecker(10, 16)), BenchmarkTools.Trial)
    =#
    info("Running BFS for LightGraphs Graph")
    @test isa(benchmark(benchmarkspec, LGGraph{LightGraphs.Graph}, BFS(), Kronecker(10, 16)), BenchmarkTools.Trial)
    info("Running BFS for LightGraphs DiGraph")
    @test isa(benchmark(benchmarkspec, LGGraph{LightGraphs.DiGraph}, BFS(), Kronecker(10, 16)), BenchmarkTools.Trial)
end


text:
module TerminalExtensions

#
# None of these functions are operating system specific because they
# might be connect via e.g. SSH to a different client operating system
#

import REPL

const DCS = "\eP"
const ST  = "\e\\"

function readDCS(io::IO)
    while bytesavailable(stdin) >= 2
        c1 = read(io,UInt8)
        c1 == 0x90 && return true
        if c1 == UInt8('\e')
            read(io,UInt8) == UInt8('P') && return true
        end
    end
    return false
end

function readST(io::IO)
    c1 = read(io,UInt8)
    c1 == 0x90 && return true
    c1 != UInt8('\e') && return false
    read(io,UInt8) != UInt8('\\') && return false
    return true
end

module iTerm2

    import REPL: display
    import Base64

    struct InlineDisplay <: AbstractDisplay
        io::IO
    end
    InlineDisplay() = InlineDisplay(stdout)

    function set_mark()
        "\033]50;SetMark\007"
    end

    # Runs after interactively edited command but before execution
    function preexec()
        "\033]133;C\007"
    end

    function remotehost_and_currentdir()
        return string("\033]1337;RemoteHost=", ENV["USER"], "@",
                      read(`hostname -f`, String), "\007",
                      "\033]1337;CurrentDir=", pwd(), "\007")
    end

    function prompt_prefix(last_success = true)
        return string("\033]133;D;$(convert(Int, last_success))\007",remotehost_and_currentdir(),"\033]133;A\007")
    end

    function prompt_suffix()
        return "\033]133;B\007"
    end

    function shell_version_number()
        return "\033]1337;ShellIntegrationVersion=1;shell=julia\007"
    end

    function prepare_display_file(io::IO=stdout;filename="Unnamed file", size=nothing, width=nothing, height=nothing, preserveAspectRatio::Bool=true, inline::Bool=false)
        q = "\e]1337;File="
        options = String[]
        filename != "Unnamed file" && push!(options,"name=" * Base64.base64encode(filename))
        size !== nothing && push!(options,"size=" * dec(size))
        height !== nothing && push!(options,"height=" * height)
        width !== nothing && push!(options,"width=" * width)
        preserveAspectRatio !== true && push!(options,"preserveAspectRatio=0")
        inline !== false && push!(options,"inline=1")
        q *= join(options,';')
        q *= ":"
        write(io,q)
    end

    function display_file(data::Vector{UInt8}; io::IO=stdout, kwargs...)
        prepare_display_file(io;kwargs...)
        write(io, Base64.base64encode(data))
        write(io,'\a')
    end

    # Incomplete list. Will be extended as necessity comes up
    const iterm2_mimes = ["image/png", "image/gif", "image/jpeg", "application/pdf", "application/eps"]

    for mime in iterm2_mimes
        @eval begin
            function display(d::InlineDisplay, m::MIME{Symbol($mime)}, x)
                buf = IOBuffer()
                show(Base64.Base64EncodePipe(buf),m,x)
                prepare_display_file(d.io;filename="image",inline=true)
                write(d.io, take!(buf))
                write(d.io,'\a')
            end
        end
    end

    function display(d::InlineDisplay,x)
        for m in iterm2_mimes
            if showable(m,x)
                return display(d,m,x)
            end
        end
        throw(MethodError(display, (d,x)))
    end

end

function __init__()
    if !(isinteractive() && isdefined(Base, :active_repl))
        return
    end
    begin
        term = REPL.Terminals.TTYTerminal("xterm",stdin,stdout,stderr)
        REPL.Terminals.raw!(term,true)
        Base.start_reading(stdin)

        # Detect iTerm support
        print(stdout, "\e[1337n\e[5n")
        readuntil(stdin, "\e")
        itermname = ""
        c = read(stdin, Char)
        c1 = read(stdin, Char)
        if c == '[' && c1 != '0'
            itermname = string(c1, readuntil(stdin, "\e")[1:end-2])
            read(stdin, Char); read(stdin, Char)
        end
        # Read the rest of the \e[5n query
        read(stdin, Char)

        if startswith(itermname, "ITERM2")
            # Inform iTerm of the shell integration version and that we're julia
            write(stdout, iTerm2.shell_version_number())

            pushdisplay(iTerm2.InlineDisplay())
            repl = Base.active_repl

            if !isdefined(repl,:interface)
                repl.interface = REPL.setup_interface(repl)
            end

            let waserror = false
                prefix = repl.interface.modes[1].prompt_prefix
                repl.interface.modes[1].prompt_prefix = function ()
                    (TerminalExtensions.iTerm2.prompt_prefix(waserror) * (isa(prefix,Function) ? prefix() : prefix))
                end
                suffix = repl.interface.modes[1].prompt_suffix
                repl.interface.modes[1].prompt_suffix = function ()
                    ((isa(suffix,Function) ? suffix() : suffix) * TerminalExtensions.iTerm2.prompt_suffix())
                end
                for mode in repl.interface.modes
                    if isdefined(mode,:on_done)
                        of = mode.on_done
                        mode.on_done = function (args...)
                            print(stdout,TerminalExtensions.iTerm2.preexec())
                            of(args...)
                            waserror = repl.waserror
                        end
                    end
                end
            end
        end
    end
end

export iTerm2

end # module


text:
abstract type AbstractArduino end

Base.isopen(a::AbstractArduino) = isopen(a.sp)
Base.close(a::AbstractArduino) = isopen(a) && close(a.sp)
function Base.open(a::AbstractArduino)
    close(a)
    a.sp = LibSerialPort.open(a.port, baudrate[])
    empty!(a.pwm.listeners)
    connect!(a)
end
restart(a::AbstractArduino) = (close(a); open(a))


text:
module Scilabs

using Compat
using DataStructures
using Distributions
using GLM
using Formatting
using MAT
using Glob
using Match
using Pipe

# using Compose
# using Memoize

export Web, writemime, load_md_table, glob, @pipe, @ls, @_, gpath, load_csv_table
export simplify_test_data
export OrderedDict, disp
export sort, filter

export FN,FN21,FN21,FN22,FN32,FN42,FN52,FN53,FN54

## Notebook Web Content
type Web
   content::String

   Web(s::String) = new(s)

   Web(s) = new(stringmime("text/html",s))
end

import Base.writemime

function writemime(io::IO, ::MIME"text/html", x::Web)
    write(io, x.content);
end

function disp(x)
	display(Web(x))
end

disp(x) = display(Web(x))

# Function
function FN(x::Real, n=2, f=2)
    abs(x) < 1/10^n || abs(x) > 10^n ? fmt(FormatSpec("$(n).$(f)e"), x) : fmt(FormatSpec("$(n).$(f)f"), x)
end

FN(n::Int=2, f::Int=2) = x-> FN(x, n, f)
FN21(x) = FN(x, 2,1)
FN22(x) = FN(x, 2,2)
FN32(x) = FN(x, 3,2)
FN42(x) = FN(x, 4,2)
FN52(x) = FN(x, 5,2)
FN53(x) = FN(x, 5,3)
FN54(x) = FN(x, 5,4)

## Path Helpers
function gpath(xs; path=".")
    res = glob(xs, path)
    isempty(res) && error("File not found for pattern: `$xs` in path: `$(path |> realpath)`")
    res |> first
end

macro ls()
    :( readdir() )
end

import Base.filter

function contains(s::String)
    function containsCurry(xs)
        return contains(xs, s)
    end
    return containsCurry
end

function filter(foo::Function)
    function filterCurry(xs)
        return filter(foo, xs)
    end
    return filterCurry
end

import Base.sort
function sort(foo::Function)
    function sortCurry(xs)
        return sort(foo, xs)
    end
    return sortCurry
end

macro _(xs)
    :( _ -> $(xs) )
end

## Array indexes
import Base.(.==)
import Base.∪

≈(x,y::Real) = isapprox(x, y)
≈(x,yt::Tuple) = isapprox(x, yt[1], atol=yt[2])
≈(xs::Array,y::Int) = find(x-> trunc(x) == y, xs)
≈(xs::Array,y::Real) = find(x->isapprox(x, y), xs)
≈(xs::Array,yt::FloatRange) = find(x->isapprox(x, yt.start, atol=yt.step), xs)

±(x,y) = range(x,y,0)
.==(arr::Array,y::FloatRange) = [ isapprox(x,y.start, atol=y.step) for x in arr ]

⊂(x::String, xs::String) = contains(xs, x)
⊂(x::String, name::Symbol) = xs -> contains(xs, x)

export ⊂, ±, .==, ≈, ∪, ⊂

## Dictionary Handling ##

markdownsplitter(x::String) = @pipe x |> strip(_) |> strip(_,'|') |> split(_,'|') |> map(strip,_)
dictFromArray(fields) = x -> OrderedDict(zip(fields,x))
arrayFromDictKeyArray(d::Dict, keys::Array, exclude...) = [ d[k] for k in filter(x-> !(x in exclude),keys) ]

⊂(d::Dict, keys::Array) = Dict( [ k=>v for (k,v) in filter((x,y)-> in(x,Set(keys)),d) ] )
∪(u::Dict, v::Dict) = Dict( concat(u ⊂ collect(keys(v)),v) )

## Matrix Helpers

simplify_mat(d::Dict{ASCIIString, Any}) = Dict{Symbol, Any}( [ symbol(string(k)) => simplify_mat(v) for (k,v) in d ] )
simplify_mat(arr::Array{Float64,2}) = arr[:] # slicing???
simplify_mat(o::Any) = o

## Markdown Helpers
function simplify_test_data(testdata::Dict)
    Dict( [ symbol(string(k)) => simplify_mat(v) for (k,v) in testdata])
end

function load_md_table_raw(testfile; f::Function=(x)->x)
    tests = open("$testfile") do f readlines(f) end
    headers = tests[1] |> markdownsplitter |> xs->map(lowercase, xs) |> xs->map(symbol, xs)
    @show headers

    @pipe( tests[3:end]
        |> map(markdownsplitter,_)
        |> filter( l -> any(  x -> !isempty(x), l), _)
        |> map(f,_)
        |> map(dictFromArray(headers),_)
    )
end

function load_md_table(testfile)

    function parser(x)
        asint(i) = try parse(Int, i) catch () end
        asfloat(i) = try parse(Float, i) catch () end
        y = ()
        if (isempty(y)) y = asint(x) end
        if (isempty(y)) y = asfloat(x) end
        if (isempty(y)) y = string(x) end
        return y
    end

    load_md_table_raw(testfile, f= x->map(parser,x))
end

function load_csv_table(testfile; usesymbols=false, postProcessor=x->x)
    dataMat = readcsv(testfile)

    if usesymbols
        headers = Symbol[ symbol(lowercase(replace(d, " ", ""))) for d in dataMat[1,:] ]
    else
        headers = String[ string(d) for d in dataMat[1,:] ]
    end

    makeRow = dictFromArray(headers)
    @show(headers, size(dataMat))

    dat = @compat Associative{eltype(headers), Any}[]
    for i in 2:size(dataMat)[1]
        row = dataMat[i,:]
        push!(dat, makeRow(row |> postProcessor))
    end

    return dat
end

include("Fitting.jl")
include("PyPlotHelpers.jl")

end


text:
module Nernst

using LinearAlgebra

EXP(x)  = cospi(x)+im*sinpi(x)

    include("AnomalousNernst_1st.jl")
    include("AnomalousHall.jl")
    include("berry.jl")
    include("derivative.jl")
    include("fastplot.jl")
    include("nernst_2nd.jl")
end


text:
# ---
# title: 834. Sum of Distances in Tree
# id: problem834
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/sum-of-distances-in-tree/description/>
# hidden: true
# ---
#
# An undirected, connected tree with `N` nodes labelled `0...N-1` and `N-1`
# `edges` are given.
#
# The `i`th edge connects nodes `edges[i][0] `and` edges[i][1]` together.
#
# Return a list `ans`, where `ans[i]` is the sum of the distances between node
# `i` and all other nodes.
#
# **Example 1:**
#
#
#
#     Input: N = 6, edges = [[0,1],[0,2],[2,3],[2,4],[2,5]]
#     Output: [8,12,6,10,10,10]
#     Explanation:
#     Here is a diagram of the given tree:
#       0
#      / \
#     1   2
#        /|\
#       3 4 5
#     We can see that dist(0,1) + dist(0,2) + dist(0,3) + dist(0,4) + dist(0,5)
#     equals 1 + 1 + 2 + 2 + 2 = 8.  Hence, answer[0] = 8, and so on.
#
#
# Note: `1 <= N <= 10000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

@test big(2.0)^big(3) == 8

for T in [Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128, BigInt]
    @test T(2)^big(3.0) == 8
    @test big(2.0)^T(3) == 8
end

# issue 15659
@test (setprecision(53) do; big(1/3); end) < 1//3


text:
using VegaLite
using CSV
using DataFrames

# This file is
data = CSV.read(joinpath("data", "oss", "working", "for_figure.csv"), DataFrame)

plt = data |>
    @vlplot(:line,
            x = {"year", axis = {title = "Year", format=:c}},
            y = {"repos", axis = {title = "Repositories Contributed To", format = "s"}},
            color = {"sector", axis = {title = "Series", format = "n"}, legend = {orient = :bottom, title = nothing}},
            height = 275,
            width = 325,
            title = "Federal Government Contribution to OSS on GitHub")
plt |> save(joinpath("repos.svg"))

plt = data |>
    @vlplot(:line,
            x = {"year", axis = {title = "Year", format=:c}},
            y = {"additions", axis = {title = "Lines Added", format = "s"}},
            color = {"sector", axis = {title = nothing, format = "n"}, legend = {orient = :bottom, title = nothing}},
            height = 275,
            width = 325,
            title = "Federal Government Contribution to OSS on GitHub")
plt |> save(joinpath("additions.svg"))
@vlplot(
    :bar,
    data={
        values=[
            {task="A",start=1,stop=3},
            {task="B",start=3,stop=8},
            {task="C",start=8,stop=10}
        ]
    },
    y="task:o",
    x="start:q",
    x2="stop:q"
)


text:
const ADSUPPORT = (:ForwardDiff, :Zygote)
const ADAVAILABLE = Dict{Module, Function}()

Hamiltonian(metric::M, ℓπ::T, m::Module) where {M<:AbstractMetric,T} = ADAVAILABLE[m](metric, ℓπ)

function Hamiltonian(metric::AbstractMetric, ℓπ)
    available = collect(keys(ADAVAILABLE))
    if length(available) == 0
        support_list_str = join(ADSUPPORT, " or ")
        error("MethodError: no method matching Hamiltonian(metric::AbstractMetric, ℓπ) because no backend is loaded. Please load an AD package ($support_list_str) first.")
    elseif length(available) > 1
        available_list_str = join(keys(ADAVAILABLE), " and ")
        constructor_list_str = join(map(m -> "Hamiltonian(metric, ℓπ, $m)", available), "\n  ")
        error("MethodError: Hamiltonian(metric::AbstractMetric, ℓπ) is ambiguous because multiple AD pakcages are available ($available_list_str). Please use AD explictly. Candidates:\n  $constructor_list_str")
    else
        return Hamiltonian(metric, ℓπ, first(available))
    end
end


text:
# Copied from
# https://github.com/kul-forbes/ProximalOperators.jl/tree/master/benchmark
using ArgParse
using PkgBenchmark
using Markdown

function displayresult(result)
    md = sprint(export_markdown, result)
    md = replace(md, ":x:" => "❌")
    md = replace(md, ":white_check_mark:" => "✅")
    display(Markdown.parse(md))
end

function printnewsection(name)
    println()
    println()
    println()
    printstyled("▃" ^ displaysize(stdout)[2]; color=:blue)
    println()
    printstyled(name; bold=true)
    println()
    println()
end

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--target"
            help = "the branch/commit/tag to use as target"
            default = "HEAD"
        "--baseline"
            help = "the branch/commit/tag to use as baseline"
            default = "main"
        "--retune"
            help = "force re-tuning (ignore existing tuning data)"
            action = :store_true
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()

    mkconfig(; kwargs...) =
        BenchmarkConfig(
            env = Dict(
                "JULIA_NUM_THREADS" => "1",
            );
            kwargs...
        )

    target = parsed_args["target"]
    group_target = benchmarkpkg(
        dirname(@__DIR__),
        mkconfig(id = target),
        resultfile = joinpath(@__DIR__, "result-$(target).json"),
        retune = parsed_args["retune"],
    )

    baseline = parsed_args["baseline"]
    group_baseline = benchmarkpkg(
        dirname(@__DIR__),
        mkconfig(id = baseline),
        resultfile = joinpath(@__DIR__, "result-$(baseline).json"),
    )

    printnewsection("Target result")
    displayresult(group_target)

    printnewsection("Baseline result")
    displayresult(group_baseline)

    judgement = judge(group_target, group_baseline)

    printnewsection("Judgement result")
    displayresult(judgement)
end

main()

# # This file was adapted from Transducers.jl
# # which is available under an MIT license (see LICENSE).
# using PkgBenchmark
# include("pprinthelper.jl")

# mkconfig(; kwargs...) =
#     BenchmarkConfig(
#         env = Dict(
#             "JULIA_NUM_THREADS" => "1",
#         );
#         kwargs...
#     )

# group_target = benchmarkpkg(
#     dirname(@__DIR__),
#     mkconfig(),
# )

# # @show readdir("benchmark")
# # out = run(`git diff`)
# # @show read(out, String)
# # out = run(`cat benchmark/Project.toml`)
# # @show read(out, String)

# group_baseline = benchmarkpkg(
#     dirname(@__DIR__),
#     mkconfig(id = "baseline"),
# )

# judgement = judge(group_target, group_baseline)

# displayresult(judgement)

# printnewsection("Target result")
# displayresult(group_target)

# printnewsection("Baseline result")
# displayresult(group_baseline)


text:
######### Stan program example  ###########

using StanSample

eightschools ="
data {
  int<lower=0> J; // number of schools
  real y[J]; // estimated treatment effects
  real<lower=0> sigma[J]; // s.e. of effect estimates
}
parameters {
  real mu;
  real<lower=0> tau;
  real eta[J];
}
transformed parameters {
  real theta[J];
  for (j in 1:J)
    theta[j] <- mu + tau * eta[j];
}
model {
  eta ~ normal(0, 1);
  y ~ normal(theta, sigma);
}
"

schools8data = Dict("J" => 8,
  "y" => [28,  8, -3,  7, -1,  1, 18, 12],
  "sigma" => [15, 10, 16, 11,  9, 11, 10, 18],
  "tau" => 25
)

sm = SampleModel("schools8", eightschools)

(sample_file, log_file) = stan_sample(sm, data=schools8data)

if !(sample_file == nothing)
  chns = read_samples(sm)

  chn = set_section(chns, Dict(
    :parameters => ["mu", "tau"],
    :thetas => ["theta.$i" for i in 1:8],
    :etas => ["eta.$i" for i in 1:8],
    :internals => ["lp__", "accept_stat__", "stepsize__", "treedepth__", "n_leapfrog__",
      "divergent__", "energy__"]
    )
  )

  describe(chn)

end


text:
using Test

import Refactoring: econvert, variables, assignements

@testset "extract method" begin

eval_type = s -> @eval Main typeof($s)

ex = Meta.parse("x = 1")
ex = econvert(ex)
@test variables(ex) == [:x]
@test assignements(ex) == [:x]

ex = Meta.parse("x,y = 1,2")
ex = econvert(ex)
@test variables(ex) == [:x,:y]
@test assignements(ex) == [:x,:y]

ex = Meta.parse("f(x) = x")
ex = econvert(ex)
@test variables(ex) == []
@test assignements(ex) == [:f]

ex = Meta.parse("f(x=1, y=rand(N))")
ex = econvert(ex)
@test variables(ex) == [:N]
@test assignements(ex) == []

ex = quote
    x,y,z,k = 1,c,3,4
    g = x*d
end
ex = econvert(ex)
@test variables(ex) == [:x,:y,:z,:k,:c,:g,:d]
@test assignements(ex) == [:x,:y,:z,:k,:g]
@test unassigned_variables(ex, eval_type) == [:c,:d]

@test unassigned_variables(quote
    x,y,z,k = 1,2,3,4
    g = x*y
end, eval_type) == Symbol[]

@test unassigned_variables(quote
    x.a = b
end, eval_type) == [:x,:b]

if !@isdefined Geom
    @eval Main module Geom end
end

ex = quote
    using Gadfly
    mu = K1
    const alpha = 1
    const beta = (alpha+1)*mu
    for i=1:10
        x = rand(100, f(K2))
    end
    plot(x = K3, y = pdf(InverseGamma(K4,beta),K5), Geom.line)
    ind, K7 = K6.asd
    x = [K8[K7] for i=1:K6]
end
@test unassigned_variables(ex, eval_type) == [:K1,:K2,:K3,:K4,:K5,:K6,:K8]

ex = quote
    pmin = pinit
    para = zeros(length(pmin),Niter)
    lik, lr = Float64[], Float64[]
    res = DiffBase.GradientResult(pmin)
    w = ones(length(pmin)); #w[1] *= 0.1
end
@test unassigned_variables(ex, eval_type) == [:pinit, :Niter]

ex = quote
    x = f(Int, Base, sin)
    y = cos(x)
end
@test unassigned_variables(ex, eval_type) == [:sin]

# comprehensions
ex = quote
    x = [i for i=1:10]
end
@test unassigned_variables(ex, eval_type) == []

ex = quote
    x = [i for i=1:N]
end
@test unassigned_variables(ex, eval_type) == [:N]

ex = quote
    x = [f(i,k) for i=1:length(N)]
end
@test unassigned_variables(ex, eval_type) == [:k,:N]

ex = quote
    f(x,y) = x*y + a
    ex = f(1,2)
end
@test unassigned_variables(ex, eval_type) == [:a]

ex = quote
    sel  = df[Symbol("col")] .== "x"
end
@test unassigned_variables(ex, eval_type) == [:df]

ex = quote
    sel  = df[Symbol(S)] .== "x"
end
@test unassigned_variables(ex, eval_type) == [:df, :S]

ex = quote
    sel  = df[Symbol(S)] .== x
end
@test unassigned_variables(ex, eval_type) == [:df, :S, :x]

ex = quote
    for y in umks
        for umk in y
            u1, u2 = split_umk(umk)
            d_ = Mod.distance(u1,u2)
            P[length(u1),length(u2)] += 1
            if d_ > 3
                push!(j,(u1,u2))
            end
            push!(d, d_)
        end
    end
end
@test unassigned_variables(ex, eval_type) == [:umks, :P, :j, :d]

ex =  quote
    for var_idx = 1:length(pv_fp)
        x[i] = 1
    end
end
@test unassigned_variables(ex, eval_type) == [:pv_fp, :i]

m = extract_method("x = sin(y)")
@test m ==
"
function (y)
    x = sin(y)
end
"

end


text:
#=
# 1D advection matching the example in the book: Nodal Discontinuous Galerkin Method.
=#
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("advection1d");

@useLog("advection1dlog")

@domain(1)
@solver(DG)
@functionSpace(LEGENDRE, 2)
@stepper(LSRK4)

@mesh(LINEMESH, 10, 2, [0,2])

@variable(u)
@testSymbol(v)

T = 1;
@timeInterval(T)
@initial(u, "sin(x)")

@boundary(u, 1, DIRICHLET, "-sin(2*pi*t)")
@boundary(u, 2, NO_BC)

@coefficient(a, 2*pi)
@coefficient(beta, 0.5)

###################################################
# Uncomment one of the following two groups of code(or both)

### To generate the Julia code and export it to a set of files, use these two lines.
@weakForm(u, "Dt(u*v) + a*grad(u)*v + surface(-a*normal()*u*v + a*normal()*ave(u)*v - 0.5*(1-beta)*a*jump(u)*v)")
@exportCode("advec1d")
###

### To import previously generated or modified code, use this line
#@importCode("advec1d")
###

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

solve(u);

# # Uncomment to plot
# using Plots
# pyplot();
# display(plot(Femshop.grid_data.allnodes[1,:], u.values[1,:], marker=:circle, reuse=false))

@finalize()


text:
using .CUDA

change_vector_eltype(S0::Type{<:CUDA.CuVector}, T) = S0.name.wrapper{T, 1, CUDA.Mem.DeviceBuffer}

convert_mat(M::CUDA.CUSPARSE.CuSparseMatrixCSC, T) = CUDA.CUSPARSE.CuSparseMatrixCSC(
  convert(CUDA.CuArray{Int, 1, CUDA.Mem.DeviceBuffer}, M.colPtr),
  convert(CUDA.CuArray{Int, 1, CUDA.Mem.DeviceBuffer}, M.rowVal),
  convert(CUDA.CuArray{T, 1, CUDA.Mem.DeviceBuffer}, M.nzVal),
  M.dims,
)
convert_mat(M::CUDA.CUSPARSE.CuSparseMatrixCSR, T) = CUDA.CUSPARSE.CuSparseMatrixCSR(
  convert(CUDA.CuArray{Int, 1, CUDA.Mem.DeviceBuffer}, M.rowPtr),
  convert(CUDA.CuArray{Int, 1, CUDA.Mem.DeviceBuffer}, M.colVal),
  convert(CUDA.CuArray{T, 1, CUDA.Mem.DeviceBuffer}, M.nzVal),
  M.dims,
)
convert_mat(M::CUDA.CuMatrix, T) = convert(typeof(M).name.wrapper{T, 2, CUDA.Mem.DeviceBuffer}, M)

function sparse_dropzeros(rows, cols, vals::CuVector, nrows, ncols)
  CPUvals = Vector(vals)
  M = sparse(rows, cols, CPUvals, ncols, nrows)
  dropzeros!(M)
  MGPU = CUDA.CUSPARSE.CuSparseMatrixCSR(M)
  return MGPU
end

function get_diag_Q_dense(Q::CUDA.CUSPARSE.CuSparseMatrixCSR{T}) where {T <: Real}
  n = size(Q, 1)
  diagval = CUDA.zeros(T, n)
  fill_diag_Q_dense!(Q.rowPtr, Q.colVal, Q.nzVal, diagval, n)
  return diagval
end

function fill_diag_Q_dense!(
  Q_rowPtr,
  Q_colVal,
  Q_nzVal::CUDA.CuVector{T},
  diagval::CUDA.CuVector{T},
  n,
) where {T <: Real}
  function kernel(Q_rowPtr, Q_colVal, Q_nzVal, diagval, n)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    k = Q_rowPtr[index + 1] - 1
    if k > 0
      i = Q_colVal[k]
      if index == i
        diagval[index] = Q_nzVal[k]
      end
    end
    return nothing
  end

  threads = min(n, 256)
  blocks = ceil(Int, n / threads)
  @cuda name = "diagq" threads = threads blocks = blocks kernel(
    Q_rowPtr,
    Q_colVal,
    Q_nzVal,
    diagval,
    n,
  )
end

function check_bounds(x::T, lvar, uvar) where {T <: Real}
  ϵ = T(1.0e-4)
  if lvar >= x
    x = lvar + ϵ
  end
  if x >= uvar
    x = uvar - ϵ
  end
  if (lvar < x < uvar) == false
    x = (lvar + uvar) / 2
  end
  return x
end

# starting points
function update_rngbounds!(x, irng, lvar, uvar, ϵ) where {T <: Real}
  @views broadcast!(check_bounds, x[irng], x[irng], lvar[irng], uvar[irng])
end

# α computation (in iterations.jl)
function ddir(dir_vi::T, vi::T) where {T <: Real}
  if dir_vi < zero(T)
    α_new = -vi * T(0.999) / dir_vi
    return α_new
  end
  return one(T)
end

function compute_α_dual_gpu(v, dir_v, store_v)
  map!(ddir, store_v, dir_v, v)
  return minimum(store_v)
end

function pdir_l(dir_vi::T, lvari::T, vi::T) where {T <: Real}
  if dir_vi < zero(T)
    α_new = (lvari - vi) * T(0.999) / dir_vi
    return α_new
  end
  return one(T)
end

function pdir_u(dir_vi::T, uvari::T, vi::T) where {T <: Real}
  if dir_vi > zero(T)
    α_new = (uvari - vi) * T(0.999) / dir_vi
    return α_new
  end
  return one(T)
end

function compute_α_primal_gpu(v, dir_v, lvar, uvar, store_v)
  map!(pdir_l, store_v, dir_v, lvar, v)
  α_l = minimum(store_v)
  map!(pdir_u, store_v, dir_v, uvar, v)
  α_u = minimum(store_v)
  return min(α_l, α_u)
end

@inline function compute_αs_gpu(
  x::CuVector,
  s_l::CuVector,
  s_u::CuVector,
  lvar::CuVector,
  uvar::CuVector,
  Δxy::CuVector,
  Δs_l::CuVector,
  Δs_u::CuVector,
  nvar,
  store_vpri::CuVector,
  store_vdual_l::CuVector,
  store_vdual_u::CuVector,
)
  α_pri = @views compute_α_primal_gpu(x, Δxy[1:nvar], lvar, uvar, store_vpri)
  α_dual_l = compute_α_dual_gpu(s_l, Δs_l, store_vdual_l)
  α_dual_u = compute_α_dual_gpu(s_u, Δs_u, store_vdual_u)
  return α_pri, min(α_dual_l, α_dual_u)
end

# dot gpu with views
function dual_obj_gpu(
  b,
  y,
  xTQx_2,
  s_l,
  s_u,
  lvar,
  uvar,
  c0,
  ilow,
  iupp,
  store_vdual_l,
  store_vdual_u,
)
  store_vdual_l .= @views lvar[ilow]
  store_vdual_u .= @views uvar[iupp]
  dual_obj = dot(b, y) - xTQx_2 + dot(s_l, store_vdual_l) - dot(s_u, store_vdual_u) + c0
  return dual_obj
end

import LinearAlgebra.rdiv!, LinearAlgebra.ldiv!

function LinearAlgebra.rdiv!(M::CUDA.AbstractGPUArray{T}, D::Diagonal) where {T}
  D.diag .= one(T) ./ D.diag
  rmul!(M, D)
  D.diag .= one(T) ./ D.diag
  return M
end

function LinearAlgebra.ldiv!(D::Diagonal, M::CUDA.AbstractGPUArray{T}) where {T}
  D.diag .= one(T) ./ D.diag
  lmul!(D, M)
  D.diag .= one(T) ./ D.diag
  return M
end


text:
function workerPrinterCheck()
    f = @spawn myid()
    println("====")
    println(fetch(f))
end

function returnArrayOne()
    return [myid(), myid() + 1]
end

function changeToOne(cutest)
    cutest[1] = 8
    cutest[2] = 9
    return nothing
end

function check_leaks()
    if length(DistributedArrays.refs) > 0
        sleep(1)  # allow time for any cleanup to complete and test again
        length(DistributedArrays.refs) > 0 && @warn("Probable leak of ", length(DistributedArrays.refs), " darrays")
    end
end

@inline function updateLocalGhost(loc_ghost_holder, dist_globaldata)
    localkeys = keys(loc_ghost_holder[1])
    # println(localkeys)
    for iter in localkeys
        #Dict To Array Equality
        loc_ghost_holder[1][iter] = dist_globaldata[loc_ghost_holder[1][iter].globalID]
    end
    return nothing
end

@inline function updateLocalGhostQ(loc_ghost_holder, loc_globaldata_mutable, dist_q)
    localkeys = keys(loc_ghost_holder[1])
    # println(localkeys)
    for iter in localkeys
        #Dict To Array Equality
        @. loc_ghost_holder[1][iter].q = dist_q[loc_ghost_holder[1][iter].globalID].q
        @. loc_globaldata_mutable[9:12, iter] = loc_ghost_holder[1][iter].q
    end
    return nothing
end

@inline function updateLocalGhostDQ(loc_ghost_holder, loc_globaldata_mutable, dist_dq)
    localkeys = keys(loc_ghost_holder[1])
    # println(localkeys)
    for iter in localkeys
        #Dict To Array Equality
        @. loc_ghost_holder[1][iter].dq = dist_dq[loc_ghost_holder[1][iter].globalID].dq
        @. loc_globaldata_mutable[13:16, iter] = loc_ghost_holder[1][iter].dq[1]
        @. loc_globaldata_mutable[17:20, iter] = loc_ghost_holder[1][iter].dq[2]
    end
    return nothing
end

function getInitialPrimitive(configData)
    rho_inf = configData["core"]["rho_inf"]::Float64
    mach = configData["core"]["mach"]::Float64
    machcos::Float64 = mach * cos(calculateTheta(configData))
    machsin::Float64 = mach * sin(calculateTheta(configData))
    pr_inf = configData["core"]["pr_inf"]::Float64
    primal = [rho_inf, machcos, machsin, pr_inf]
    return primal
end

function getInitialPrimitive2(configData)
    dataman = open("prim_soln_clean")
    data = read(dataman, String)
    data1 = split(data, "\n")
    finaldata = Array{Array{Float64,1},1}(undef, 0)
    for (idx,itm) in enumerate(data1)
        # try
        da = split(itm)
        da1 = parse.(Float64, da)
        push!(finaldata, da1)
    end
    close(dataman)
    return finaldata
end

function matchInitialQ(loc_q, loc_globaldata)
    local_size = length(loc_q)
    for idx in 1:local_size
        loc_q[idx].q = loc_globaldata[idx].q
    end
end

function placeNormals(loc_globaldata, globaldata, loc_ghost_holder, interior, wall, outer)
    local_size = length(loc_globaldata)
    # println(local_size, " is the local size")
    updateLocalGhost(loc_ghost_holder, globaldata)
    for idx in 1:local_size
        flag = loc_globaldata[idx].flag_1
        if flag == wall || flag == outer
            currpt = getxy(loc_globaldata[idx])
            leftpt = loc_globaldata[idx].left
            rightpt = loc_globaldata[idx].right

            if leftpt <= local_size
                leftpt = getxy(loc_globaldata[leftpt])
            else
                leftpt = getxy(loc_ghost_holder[1][leftpt])
            end

            if rightpt <= local_size
                rightpt = getxy(loc_globaldata[rightpt])
            else
                rightpt = getxy(loc_ghost_holder[1][rightpt])
            end

            normals = calculateNormals(leftpt, rightpt, currpt[1], currpt[2])
            setNormals(loc_globaldata[idx], normals)
        elseif flag == interior
            setNormals(loc_globaldata[idx], (0,1))
        else
            @warn "Illegal Point Type"
        end
    end
    return nothing
end

function calculateNormals(left, right, mx, my)
    lx = left[1]
    ly = left[2]

    rx = right[1]
    ry = right[2]

    nx1 = my - ly
    nx2 = ry - my

    ny1 = mx - lx
    ny2 = rx - mx

    nx = 0.5*(nx1 + nx2)
    ny = 0.5*(ny1 + ny2)

    det = hypot(nx, ny)

    nx = -nx/det
    ny = ny/det

    return (nx,ny)
end

function calculateConnectivity(loc_globaldata, globaldata, loc_ghost_holder)
    local_size = length(loc_globaldata)

    for idx in 1:local_size
        ptInterest = loc_globaldata[idx]
        currx = ptInterest.x
        curry = ptInterest.y
        nx = ptInterest.nx
        ny = ptInterest.ny

        flag = ptInterest.flag_1

        xpos_conn,xneg_conn,ypos_conn,yneg_conn = Array{Int32,1}(undef, 0),Array{Int32,1}(undef, 0),Array{Int32,1}(undef, 0),Array{Int32,1}(undef, 0)

        tx = ny
        ty = -nx

        for itm in ptInterest.conn
            if itm <= local_size
                itmx = loc_globaldata[itm].x
                itmy = loc_globaldata[itm].y
            else
                itmx = loc_ghost_holder[1][itm].x
                itmy = loc_ghost_holder[1][itm].y
            end

            delx = itmx - currx
            dely = itmy - curry

            dels = delx*tx + dely*ty
            deln = delx*nx + dely*ny
            if dels <= 0.0
                push!(xpos_conn, itm)
            end
            if dels >= 0.0
                push!(xneg_conn, itm)
            end
            if flag == 1
                if deln <= 0.0
                    push!(ypos_conn, itm)
                end
                if deln >= 0.0
                    push!(yneg_conn, itm)
                end
            elseif flag == 0
                push!(yneg_conn, itm)
            elseif flag == 2
                push!(ypos_conn, itm)
            end
        end
        setConnectivity(loc_globaldata[idx], (xpos_conn, xneg_conn, ypos_conn, yneg_conn))
    end
    return nothing
end

function fpi_solver(iter_store, ghost_holder, dist_globaldata, dist_q, dist_dq, dist_globaldata_mutable, configData, res_old, res_new, numPoints)
    # println(IOContext(stdout, :compact => false), globaldata[3].prim)
    # print(" 111\n")
    str = CuStream()
    power = configData["core"]["power"]
    cfl = configData["core"]["cfl"]::Float64
    Mach::Float64 = configData["core"]["mach"]::Float64
    gamma::Float64 = configData["core"]["gamma"]::Float64
    pr_inf::Float64 = configData["core"]["pr_inf"]::Float64
    rho_inf::Float64 = configData["core"]["rho_inf"]::Float64
    theta = calculateTheta(configData)
    threadsperblock = parse(Int, ARGS[3])

    for iter in 1:iter_store
        if iter == 1
            println("Starting FuncDelta")
        end
    # @sync

        @sync for ip in procs(dist_globaldata)
            @spawnat ip begin
                updateLocalGhost(ghost_holder[:L], dist_globaldata)
            end
        end

        @sync for ip in procs(dist_globaldata)
            @spawnat ip begin
                # println(length(localpart(globaldata)))
                func_delta(dist_globaldata[:L], ghost_holder[:L], cfl, numPoints)
            end
        end

        for rk in 1:4
        #    # if iter == 1
        #        # println("Starting QVar")
        #    # end
            @sync for ip in procs(dist_globaldata)
                @spawnat ip begin
                    q_variables(dist_globaldata[:L], dist_q[:L], dist_globaldata_mutable[:L])
                end
            end

            @sync for ip in procs(dist_globaldata)
                @spawnat ip begin
                    updateLocalGhostQ(ghost_holder[:L], dist_globaldata_mutable[:L], dist_q)
                end
            end

            @sync for ip in procs(dist_globaldata)
                @spawnat ip begin
                    q_var_derivatives(dist_globaldata[:L], dist_dq[:L], ghost_holder[:L], dist_globaldata_mutable[:L], power)
                end
            end
        # println(IOContext(stdout, :compact => false), dist_globaldata[3])
        #    # if iter == 1
        #        # println("Starting Calflux")
        #    # end

            @sync for ip in procs(dist_globaldata)
                @spawnat ip begin
                    updateLocalGhostDQ(ghost_holder[:L], dist_globaldata_mutable[:L], dist_dq)
                end
            end

            # @sync for ip in procs(dist_globaldata)
            #     @spawnat ip begin
            #         cal_flux_residual(dist_globaldata[:L], ghost_holder[:L], dist_globaldata_mutable[:L], configData)
            #     end
            # end

            @sync for ip in procs(dist_globaldata)
                @spawnat ip begin
                    gpuLocDataRest = CuArray(dist_globaldata_mutable[:L])
                    locBlocksPerGrid = Int(ceil(gpuLocNumPoints/threadsperblock))
                    @cuda blocks= locBlocksPerGrid threads= threadsperblock cal_flux_residual_kernel(gpuLocDataConn, gpuLocDataFixedPoint, gpuLocDataRest, gpuConfigData, gpuLocNumPoints)
                    dist_globaldata_mutable[:L] = Array(gpuLocDataRest)
                end
            end
            synchronize()
        #    # println(IOContext(stdout, :compact => false), globaldata[3].prim)
        #    # println(IOContext(stdout, :compact => false), globaldata[3].prim)
        #    # residue = 0
        #    # if iter == 1
        #        # println("Starting StateUpdate")
        #    # end
            @sync for ip in procs(dist_globaldata)
                @spawnat ip begin
                    state_update(dist_globaldata[:L], dist_globaldata_mutable[:L], Mach, gamma, pr_inf, rho_inf, theta, iter, res_old[:L], res_new[:L], rk, numPoints)
                end
            end
        end
        residue = 0
        if iter <= 2
            # res_old[1] = res_new[1]
            residue = 0
        else
            residue = log10(sqrt(sum(res_new))/sqrt(sum(res_old)))
        end
        println("Iteration Number ", iter, " ", residue)
    end

    # println(IOContext(stdout, :compact => false), globaldata[3].prim)
    # residue = res_old
    return nothing
end

@inline function q_variables(loc_globaldata, loc_q, loc_globaldata_mutable)
    for (idx, itm) in enumerate(loc_globaldata)
        rho = itm.prim[1]
        u1 = itm.prim[2]
        u2 = itm.prim[3]
        pr = itm.prim[4]

        beta = 0.5 * (rho / pr)
        loc_globaldata_mutable[9, idx] = log(rho) + log(beta) * 2.5 - (beta * ((u1 * u1) + (u2 * u2)))
        two_times_beta = 2.0 * beta

        loc_globaldata_mutable[10, idx] = (two_times_beta * u1)
        loc_globaldata_mutable[11, idx] = (two_times_beta * u2)
        loc_globaldata_mutable[12, idx] = -two_times_beta
        loc_globaldata[idx].q = @views loc_globaldata_mutable[9:12, idx]
        loc_q[idx].q = loc_globaldata[idx].q
    end
    return nothing
end

function q_var_derivatives(loc_globaldata, loc_dq, loc_ghost_holder, loc_globaldata_mutable, power)
    sum_delx_delq = zeros(Float64, 4)
    sum_dely_delq = zeros(Float64, 4)
    dist_length = length(loc_globaldata)

    for (idx, itm) in enumerate(loc_globaldata)
        x_i = itm.x
        y_i = itm.y
        sum_delx_sqr = zero(Float64)
        sum_dely_sqr = zero(Float64)
        sum_delx_dely = zero(Float64)
        sum_delx_delq = fill!(sum_delx_delq, 0.0)
        sum_dely_delq = fill!(sum_dely_delq, 0.0)
        for i in 1:4
            loc_globaldata[idx].max_q[i] = loc_globaldata[idx].q[i]
            loc_globaldata_mutable[20 + i, idx] = loc_globaldata[idx].q[i]
            loc_globaldata[idx].min_q[i] = loc_globaldata[idx].q[i]
            loc_globaldata_mutable[24 + i, idx] = loc_globaldata[idx].q[i]
        end
        for conn in itm.conn
            if conn <= dist_length
                globaldata_conn = loc_globaldata[conn]
            else
                # update_ghost_q_variables(loc_ghost_holder[1][conn])
                globaldata_conn = loc_ghost_holder[1][conn]
            end
            x_k = globaldata_conn.x
            y_k = globaldata_conn.y
            delx = x_k - x_i
            dely = y_k - y_i
            dist = hypot(delx, dely)
            weights = dist ^ power
            sum_delx_sqr += ((delx * delx) * weights)
            sum_dely_sqr += ((dely * dely) * weights)
            sum_delx_dely += ((delx * dely) * weights)

            for i in 1:4
                sum_delx_delq[i] += (weights * delx * (loc_globaldata_mutable[8+i, conn] - loc_globaldata_mutable[8+i, idx]))
                sum_dely_delq[i] += (weights * dely * (loc_globaldata_mutable[8+i, conn] - loc_globaldata_mutable[8+i, idx]))
                if loc_globaldata_mutable[20 + i, idx] < loc_globaldata_mutable[8+i, conn]
                    loc_globaldata[idx].max_q[i] = loc_globaldata_mutable[8+i, conn]
                    loc_globaldata_mutable[20 + i, idx] = loc_globaldata_mutable[8+i, conn]
                end
                if loc_globaldata_mutable[24 + i, idx] > loc_globaldata_mutable[8+i, conn]
                    loc_globaldata[idx].min_q[i] = loc_globaldata_mutable[8+i, conn]
                    loc_globaldata_mutable[24 + i, idx] = loc_globaldata_mutable[8+i, conn]
                end
            end
        end
        det = (sum_delx_sqr * sum_dely_sqr) - (sum_delx_dely * sum_delx_dely)
        one_by_det = 1.0 / det
        @. loc_globaldata_mutable[13:16, idx] = one_by_det * (sum_delx_delq * sum_dely_sqr - sum_dely_delq * sum_delx_dely)
        loc_globaldata[idx].dq[1] = @views loc_globaldata_mutable[13:16, idx]
        @. loc_globaldata_mutable[17:20, idx] = one_by_det * (sum_dely_delq * sum_delx_sqr - sum_delx_delq * sum_delx_dely)
        loc_globaldata[idx].dq[2] = @views loc_globaldata_mutable[17:20, idx]
        loc_dq[idx].dq[1] = loc_globaldata[idx].dq[1]
        loc_dq[idx].dq[2] = loc_globaldata[idx].dq[2]
        # globaldata[idx].dq = [tempsumx, tempsumy]
    end
    # println(IOContext(stdout, :compact => false), globaldata[3].dq)
    # println(IOContext(stdout, :compact => false), globaldata[3].max_q)
    # println(IOContext(stdout, :compact => false), globaldata[3].min_q)
    return nothing
end

# @inline function update_ghost_q_variables(ghost_point)
#     rho = ghost_point.prim[1]
#     u1 = ghost_point.prim[2]
#     u2 = ghost_point.prim[3]
#     pr = ghost_point.prim[4]

#     beta = 0.5 * (rho / pr)
#     ghost_point.q[1] = log(rho) + log(beta) * 2.5 - (beta * ((u1 * u1) + (u2 * u2)))
#     two_times_beta = 2.0 * beta

#     ghost_point.q[2] = (two_times_beta * u1)
#     ghost_point.q[3] = (two_times_beta * u2)
#     ghost_point.q[4] = -two_times_beta
#     return nothing
# end

# function updateLocalGhostNew(loc_ghost_holder, dist_globaldata)
#     localkeys = keys(loc_ghost_holder[1])
#     # println(localkeys)
#     for iter in localkeys
#         #Dict To Array Equality
#         loc_ghost_holder[1][iter].prim = dist_globaldata[loc_ghost_holder[1][iter].globalID].prim
#         loc_ghost_holder[1][iter].q = dist_globaldata[loc_ghost_holder[1][iter].globalID].q
#     end
#     return nothing
# end

# function updateGhost(ghost_holder, dist_globaldata)
#     for i in 1:nworkers()
#         localkeys = keys(ghost_holder[i])
#         for iter in localkeys
#             ghost_holder[i][iter] = dist_globaldata[ghost_holder[i][iter].globalID]
#         end
#     end
# end

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
Match() = Match(-1, -1, Inf)

HalfMatch() = HalfMatch(-1, Inf)

function Normalizer(dfp::DfPair)
    dfs1, dfs2 = sample.([dfp.df1, dfp.df2])
    raw_dist_arr = [col_dist.(dfs1[!, col], dfs2[!, col]) for col in names(dfs1)]
    means, stds = (mean.(raw_dist_arr), std.(raw_dist_arr))
    nonzero_stds = @. ifelse(stds > 0, stds, 1)
    multiplier =
        [
            (arr .- mu) ./ rho
            for (arr, mu, rho) in zip(raw_dist_arr, means, nonzero_stds)
        ] |>
        mean |>
        std
    Normalizer(means, nonzero_stds * multiplier)
end

function SortedJoiner(arr::Array{Int, 2})
    sinds = [sortperm(arr[:, i]) for i in axes(arr, 2)]
    SortedJoiner(arr, sinds, [view(arr, sinds[i], i) for i in axes(arr, 2)])
end

SortedJoiner(df::DataFrame) = SortedJoiner(hcat(df[:, 1], df[:, 2]))

function RelPair(rels::Array, colmap::Tuple{String, String})
    RelPair(SortedJoiner.(rels)...,colmap)
end

CorefProposal() = CorefProposal(-1, fill(0.0, 5))  # Latent variable count

function CorefSpace(longside::Int, shortside::Int, preplen::Int)
    latent_match_vars = Dict(
        k=> [[CorefProposal() for _ in  1:size(df)[1]]
        for df in [dfp.df1, dfp.df2]] for (k, dfp) in crs.esp_dict
    );
    prepspace = fill(HalfMatch(),(shortside, preplen))
    free_i = fill(false, shortside)
    free_j = fill(false, longside)
    CorefSpace(latent_match_vars, prepspace, free_i, free_j, [])
end

function CorefSystem(esp_dict::Dict{String, DfPair}, relps::Array{RelPair})
    normalizers = Dict(k=>Normalizer(dfp) for (k, dfp) in esp_dict)
    CorefSystem(esp_dict, normalizers, relps, Dict(k=>Set() for k in keys(esp_dict)))
end

function CorefSystem(dir_path::String)
    esp_base, relp_base = [joinpath(dir_path, s) for s in ["esp", "relp"]]
    esp_dict = Dict(
        esp_p => DfPair(load_dfs(esp_base, esp_p)...) for esp_p in readdir(esp_base)
    )
    relps = [
        RelPair(
            SortedJoiner.(load_dfs(relp_base, relp_p))...,
            JSON.parsefile(joinpath(relp_base, relp_p, "colmap.json")) |> Tuple
            )
        for relp_p in readdir(relp_base)
    ]
    CorefSystem(esp_dict,relps)
end

function SpacePair(maxh1::Int, maxh2::Int, maxw::Int)
    SpacePair([
        EntityRefSpace(
            Array{Int,2}(undef, x, maxw),
            0,
            0,
            Array{String}(undef, maxw))
             for x in [maxh1, maxh2]]...
                )
end

function CorefResolver(
    crs::CorefSystem,
    shortside::Int,
    longside::Int,
    max_depth::Int,
    preplen::Int,
    )
    CorefResolver(
        crs,
        CorefSpace(longside, shortside, preplen),
        Tuple(SpacePair(shortside, longside, max_depth) for _ in 1:2),
        true,
    )
end

InitParams(entity_type::String) = InitParams(entity_type, typemax(Int), typemax(Int))


text:
@testset "test_complex_systems" begin

C_1 = zpk([], [-1+im], 1.0+1im)
C_2 = zpk([-1+im], [], 1.0+1im)

# Basic arittmetic
@test C_1 + C_1 == zpk([], [-1+im], 2+2im)
@test (1.0+im)*C_1 == zpk([], [-1+im], 2.0im)
@test C_2 + C_2 == zpk([-1+im], [], 2+2im)

@test C_1 * C_1 == zpk([], [-1+im,-1+im], 2.0im)
@test C_2 * C_2 ≈ zpk([-1+im,-1+im], [], 2.0im)

@test im*ss(1) == ss(im)

@test pole(zpk([], [-1+im,-1+im,0], 2.0im)) == [-1+im,-1+im,0]
@test tzero(zpk([-1+im,-1+im,0], [-2], 2.0im)) == [-1+im,-1+im,0]

@test zpk( tf([1.0, 1+im], [1.0, 2+im]) ) == zpk( [-1-im], [-2-im], 1.0+0im)

@test minreal(zpk([-1+im], [-1+im,-1+im],1+0im)) == zpk([], [-1+im],1+0im)
@test minreal(zpk([-1+im, -1+im], [-1+im],1+1im)) == zpk([-1+im], [], 1+1im)

@test_throws AssertionError zpk([-1+im], [-1+im,-1+im],1) #  Given the type of k this should be a real-coefficient system, but poles and zeros don't come in conjugate pairs

@test zpk([-2+im], [-1+im],1+0im)*zpk([], [-1+im],1+0im) == zpk([-2+im], [-1+im, -1+im], 1+0im)
@test zpk([], [-2], 2) + zpk([], [-1], 1) == zpk([-4/3], [-2,-1], 3)

@test tf(zpk([-2+im], [-1+im],1+0im)) == tf([1, 2-im], [1, 1-im])

@test 1 / ( tf("s") + 1 + im ) == tf([1], [1, 1+im])

end


text:
import StatsFuns.betacdf

seed!(12345)

function testBeta(α, β, low, high, h, N)
  function cdf(x::Float64)
    return betacdf(α, β, x)
  end
  function sampler()
    return sampleBeta(α, β)
  end
  return testGOFContinuous(cdf, sampler, low, high, h, N)
end

@test testBeta(3.0, 2.0, 0.05, 0.95, 0.01, 2^21) > 0.01
@test testBeta(0.1, 0.2, 0.05, 0.95, 0.01, 2^21) > 0.01
@test testBeta(0.8, 12.3, 0.01, 0.55, 0.01, 2^21) > 0.01
@test testBeta(10.0, 25.0, 0.1, 0.6, 0.01, 2^21) > 0.01


text:
using GJRand
using Test

@testset "GJRand.jl" begin
    # Write your tests here.
end


text:
# Second iteration of explicitly exloring the marginalization of discrete variables onto the continuous state space.
# Although this code is still excessive and messy, this is a significant feature expansion; and future versions will
# generalize the marginalization process to allow for implicit hypothesis exploration.  The messy explicit version of code is
# intended to help develop the required generalistic unit tests.  The unit tests will then be validated, frozen and used to
# confirm that future "algebraic" marginalization (implicit) versions operate correctly.
# FYI, the complexity of general multihypothesis convolutions can be deceiving, however, note that the coding
# complexity is contained for each indivual factor at a time.  Global Bayes tree inference then creates the symphony
# of non-Gaussian (multimodal) posterior beliefs from the entire factor graph.
#
# 2018/6/01 @dehann

"""
    $(SIGNATURES)

This function explicitly codes that marginalization of a discrete categorical selection variable for ambiguous data association situations.  Improved implementations should implicitly induce the same behaviour through summation (integration) when marginalizing any number of discrete variables.  This function populates `allelements` with particle indices associated with particular multihypothesis selection while `activehypo` simultaneously contains the hypothesis index and factor graph variables associated with that hypothesis selection.  The return value `certainidx` are the hypotheses that are not in question.

```
# `allelements` example BearingRange [:x1, 0.5:l1a, 0.5:l1b]
# sfidx = (1=:x1,2=:l1a,3=:l1b)
if solvefor :x1, then allelem = [mhidx.==:l1a; mhidx.==l1b]
if solvefor :l1a, then allelem = [mhidx.==:l1a] and ARR[solvefor][:,mhidx.==:l1b]=ARR[:l1b][:,mhidx.==:l1b]
if solvefor :l1b, then allelem = [mhidx.==:l1b] and ARR[solvefor][:,mhidx.==:l1a]=ARR[:l1a][:,mhidx.==:l1a]
if solvefor 1, then allelem = [mhidx.==2; mhidx.==3]
if solvefor 2, then allelem = [mhidx.==2] and ARR[solvefor][:,mhidx.==3]=ARR[3][:,mhidx.==3]
if solvefor 3, then allelem = [mhidx.==3] and ARR[solvefor][:,mhidx.==2]=ARR[2][:,mhidx.==2]

# `activehypo` in example mh=[1.0;0.5;0.5]
sfidx=1, mhidx=2:  ah = [1;2]
sfidx=1, mhidx=3:  ah = [1;3]
sfidx=2, mhidx=2:  ah = [1;2]
sfidx=2, mhidx=3:  2 should take a value from 3
sfidx=3, mhidx=2:  3 should take a value from 2
sfidx=3, mhidx=3:  ah = [1;3]

# `activehypo` in example mh=[1.0;0.33;0.33;0.34]
sfidx=1, mhidx=2:  ah = [1;2]
sfidx=1, mhidx=3:  ah = [1;3]
sfidx=1, mhidx=4:  ah = [1;4]

sfidx=2, mhidx=2:  ah = [1;2]
sfidx=2, mhidx=3:  2 should take a value from 3
sfidx=2, mhidx=4:  2 should take a value from 4

sfidx=3, mhidx=2:  3 should take a value from 2
sfidx=3, mhidx=3:  ah = [1;3]
sfidx=3, mhidx=4:  3 should take a value from 4

sfidx=4, mhidx=2:  4 should take a value from 2
sfidx=4, mhidx=3:  4 should take a value from 3
sfidx=4, mhidx=4:  ah = [1;4]
```
"""
function assembleHypothesesElements!(
            mh::Categorical,
            maxlen::Int,
            sfidx::Int,
            lenXi::Int  )
  #
  allelements = []
  activehypo = []
  mhidx = Int[]

  allidx = 1:maxlen
  allmhp = 1:length(mh.p)
  certainidx = allmhp[mh.p .< 1e-10]

  # prep mmultihypothesis selection values
  mhidx = rand(mh, maxlen) # selection of which hypothesis is correct

  pidx = 0
  for pval in mh.p
    pidx += 1
    iterarr = allidx[mhidx .== pidx]
    iterah = Int[]
    if pval >= 1e-10 && mh.p[sfidx] < 1e-10
      iterah = sort([sfidx;pidx]) # TODO -- currently only support binary factors in multihypo mode
    elseif pval < 1e-10 && mh.p[sfidx] < 1e-10
      iterarr = Int[]
      # info("assembleHypothesesElements! -- not processing pval < 1e-10 && mh.p[sfidx] < 1e-10")
    elseif pval < 1e-10 && mh.p[sfidx] >= 1e-10
      iterah = sort([sfidx;pidx]) # TODO -- currently only support binary factors in multihypo mode
    elseif pval >= 1e-10 && mh.p[sfidx] >= 1e-10
      iterah = allmhp[mh.p .> 1e-10]
    else
      error("Unknown hypothesis case, got sfidx=$(sfidx) with mh.p=$(mh.p), pidx=$(pidx)")
    end
    push!(allelements, iterarr)
    push!(activehypo, (pidx,iterah))
  end

  return certainidx, allelements, activehypo, mhidx
end

function assembleHypothesesElements!(
            mh::Void,
            maxlen::Int,
            sfidx::Int,
            lenXi::Int  )
  #
  allelements = []
  activehypo = []
  mhidx = Int[]

  allidx = 1:maxlen
  certainidx = 1:lenXi
  doneall = false
  for i in certainidx
    if !doneall
      push!(allelements, allidx)
      push!(activehypo, (i,certainidx))
      doneall = true
    else
      push!(allelements, Int[])
      push!(activehypo, (i,Int[]))
    end
  end
  return certainidx, allelements, activehypo, mhidx # certainidx = allhp
end

#


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    MappingMethod

A method for mapping spatial data onto a domain.
"""
abstract type MappingMethod end

"""
    map(sdata, sdomain, targetvars, method)

Map the `targetvars` in `sdata` to `sdomain` with mapping `method`.
"""
Base.map(data, domain, targetvars, ::MappingMethod)

# ----------------
# IMPLEMENTATIONS
# ----------------
include("mappings/nearest.jl")
include("mappings/copy.jl")


text:
"""
    to_vec(x)

Transform `x` into a `Vector`, and return the vector, and a closure which inverts the
transformation.
"""
function to_vec(x::Real)
    function Real_from_vec(x_vec)
        return first(x_vec)
    end
    return [x], Real_from_vec
end

function to_vec(z::Complex)
    function Complex_from_vec(z_vec)
        return Complex(z_vec[1], z_vec[2])
    end
    return [real(z), imag(z)], Complex_from_vec
end

# Base case -- if x is already a Vector{<:Real} there's no conversion necessary.
to_vec(x::Vector{<:Real}) = (x, identity)

function to_vec(x::AbstractVector)
    x_vecs_and_backs = map(to_vec, x)
    x_vecs, backs = first.(x_vecs_and_backs), last.(x_vecs_and_backs)
    function Vector_from_vec(x_vec)
        sz = cumsum(map(length, x_vecs))
        x_Vec = [backs[n](x_vec[sz[n] - length(x_vecs[n]) + 1:sz[n]]) for n in eachindex(x)]
        return oftype(x, x_Vec)
    end
    return vcat(x_vecs...), Vector_from_vec
end

function to_vec(x::AbstractArray)

    x_vec, from_vec = to_vec(vec(x))

    function Array_from_vec(x_vec)
        return oftype(x, reshape(from_vec(x_vec), size(x)))
    end

    return x_vec, Array_from_vec
end

# Some specific subtypes of AbstractArray.

function to_vec(x::T) where {T<:LinearAlgebra.AbstractTriangular}
    x_vec, back = to_vec(Matrix(x))
    function AbstractTriangular_from_vec(x_vec)
        return T(reshape(back(x_vec), size(x)))
    end
    return x_vec, AbstractTriangular_from_vec
end

function to_vec(x::Symmetric)
    function Symmetric_from_vec(x_vec)
        return Symmetric(reshape(x_vec, size(x)))
    end
    return vec(Matrix(x)), Symmetric_from_vec
end

function to_vec(X::Diagonal)
    function Diagonal_from_vec(x_vec)
        return Diagonal(reshape(x_vec, size(X)...))
    end
    return vec(Matrix(X)), Diagonal_from_vec
end

function to_vec(X::Transpose)
    function Transpose_from_vec(x_vec)
        return Transpose(permutedims(reshape(x_vec, size(X))))
    end
    return vec(Matrix(X)), Transpose_from_vec
end

function to_vec(X::Adjoint)
    function Adjoint_from_vec(x_vec)
        return Adjoint(conj!(permutedims(reshape(x_vec, size(X)))))
    end
    return vec(Matrix(X)), Adjoint_from_vec
end

# Non-array data structures

function to_vec(x::Tuple)
    x_vecs, x_backs = zip(map(to_vec, x)...)
    sz = cumsum(collect(map(length, x_vecs)))
    function Tuple_from_vec(v)
        return ntuple(length(x)) do n
            return x_backs[n](v[sz[n] - length(x_vecs[n]) + 1:sz[n]])
        end
    end
    return reduce(vcat, x_vecs), Tuple_from_vec
end

# Convert to a vector-of-vectors to make use of existing functionality.
function to_vec(d::Dict)
    d_vec, back = to_vec(collect(values(d)))
    function Dict_from_vec(v)
        v_vec_vec = back(v)
        return Dict(key => v_vec_vec[n] for (n, key) in enumerate(keys(d)))
    end
    return d_vec, Dict_from_vec
end


text:
module StreamingTimeSeries

import Base.Dates.TimeType

export
    EMAFeature,
    EMVFeature,
    TimeSinceFeature,
    LastValueFeature,
    DecayFeature,
    LagSumFeature,
    LagLastValueFeature,
    update!,
    valueat

# see EMA_lin in http://www.eckner.com/papers/ts_alg.pdf
type EMAFeature{T<:TimeType}
    decayRate::Float64
    lastTime::T
    lastValue::Float64
    value::Float64
end
EMAFeature{T<:TimeType}(decayRate::Float64, lastTime::T) = EMAFeature(decayRate, lastTime, 0.0, 0.0)
function update!{T<:TimeType}(feature::EMAFeature{T}, time::T, value::Float64)
    @assert time >= feature.lastTime
    if time == feature.lastTime
        return feature.value
    end

    dt = Dates.value(time - feature.lastTime)/60000 # convert from milliseconds to minutes
    a = dt / feature.decayRate
    u = exp(-a)
    v = (1 - u) / a
    feature.value = u*feature.value + (v - u)*feature.lastValue + (1.0 - v)*value
    feature.lastTime = time
    feature.lastValue = value

    feature.value
end
function valueat{T<:TimeType}(feature::EMAFeature{T}, time::T)
    @assert time >= feature.lastTime
    feature.value
end

type EMVFeature{T<:TimeType}
    ema::EMAFeature{T}
    ema2::EMAFeature{T}
    value::Float64
end
EMVFeature{T<:TimeType}(decayRate::Float64, lastTime::T) = EMVFeature(EMAFeature(decayRate, lastTime), EMAFeature(decayRate, lastTime), 0.0)
function update!{T<:TimeType}(feature::EMVFeature{T}, time::T, value::Float64)
    @assert time >= feature.ema.lastTime
    if time == feature.ema.lastTime
        return feature.value
    end

    update!(feature.ema, time, value)
    update!(feature.ema2, time, value^2)

    feature.value = feature.ema2.value - feature.ema.value^2
end
function valueat{T<:TimeType}(feature::EMVFeature{T}, time::T)
    @assert time >= feature.ema.lastTime
    feature.value
end

# simply keep track of how far we are from the reference time point
# note that we collapse all negative values to 0.0 since otherwise they
# would provide information about the future.
type TimeSinceFeature{T<:TimeType}
    referenceTime::T
end
function update!{T<:TimeType}(feature::TimeSinceFeature{T}, time::T)
    feature.referenceTime = time
end
function update!{T<:TimeType}(feature::TimeSinceFeature{T}, time::T, ignored::Float64)
    update!(feature, time)
end
function valueat{T<:TimeType}(feature::TimeSinceFeature{T}, time::T)
    numMin = Dates.value(time - feature.referenceTime)/60000 # convert from milliseconds to minutes
    numMin < 0.0 ? 0.0 : numMin
end

# simply keep track of the last value seen
type LastValueFeature{T<:TimeType}
    time::T
    value::Float64
end
LastValueFeature{T<:TimeType}(time::T) = LastValueFeature(time, 0.0)
function update!{T<:TimeType}(feature::LastValueFeature{T}, time::T, value::Float64)
    feature.time = time
    feature.value = value
end
function valueat{T<:TimeType}(feature::LastValueFeature{T}, time::T)
    time < feature.time ? 0.0 : feature.value
end

# Sum all the values and let them decay exponentially
type DecayFeature{T<:TimeType}
    decayRate::Float64
    lastTime::T
    value::Float64
end
DecayFeature{T<:TimeType}(halfLife::Float64, lastTime::T) = DecayFeature(exp(log(0.5)/halfLife), lastTime, 0.0)
function update!{T<:TimeType}(feature::DecayFeature{T}, time::T, value::Float64)
    @assert time >= feature.lastTime

    dt = Dates.value(time - feature.lastTime)/60000 # convert from milliseconds to minutes
    feature.lastTime = time
    feature.value *= feature.decayRate^dt
    feature.value += value
end
function valueat{T<:TimeType}(feature::DecayFeature{T}, time::T)
    @assert time >= feature.lastTime
    dt = Dates.value(time - feature.lastTime)/60000 # convert from milliseconds to minutes
    feature.value * feature.decayRate^dt
end

"Sum all the values in a window, the starting point is excluded and the end point included."
type LagSumFeature{T<:TimeType}
    lagStart::Dates.TimePeriod # exclusive
    lagEnd::Dates.TimePeriod # inclusive
    bufferTimes::Array{T,1}
    bufferValues::Array{Float64,1}
    windowTimes::Array{T,1}
    windowValues::Array{Float64,1}
    lastTime::T
    value::Float64
end
LagSumFeature{T<:TimeType}(lagStart::Dates.TimePeriod, lagEnd::Dates.TimePeriod, lastTime::T) = LagSumFeature(lagStart, lagEnd, T[], Float64[], T[], Float64[], lastTime, 0.0)
function update!{T<:TimeType}(feature::LagSumFeature{T}, time::T, value::Float64)
    @assert time >= feature.lastTime
    feature.lastTime = time
    unshift!(feature.bufferTimes, time)
    unshift!(feature.bufferValues, value)
end
function valueat{T<:TimeType}(feature::LagSumFeature{T}, time::T)
    @assert time >= feature.lastTime

    # move data from the buffer to the window
    for i in length(feature.bufferTimes):-1:1
        if time - feature.bufferTimes[i] >= feature.lagEnd
            unshift!(feature.windowTimes, pop!(feature.bufferTimes))
            unshift!(feature.windowValues, pop!(feature.bufferValues))
            feature.value += feature.windowValues[1]
        else
            break
        end
    end

    # move data out of the window
    for i in length(feature.windowTimes):-1:1
        if time - feature.windowTimes[i] >= feature.lagStart
            pop!(feature.windowTimes)
            feature.value -= pop!(feature.windowValues)
        else
            break
        end
    end
    feature.lastTime = time

    feature.value
end

"Report the last value with event time <= a given lag point."
type LagLastValueFeature{T<:TimeType}
    lag::Dates.TimePeriod # inclusive
    bufferTimes::Array{T,1}
    bufferValues::Array{Float64,1}
    lastTime::T
    value::Float64
end
LagLastValueFeature{T<:TimeType}(lag::Dates.TimePeriod, lastTime::T) = LagLastValueFeature(lag, T[], Float64[], lastTime, 0.0)
function update!{T<:TimeType}(feature::LagLastValueFeature{T}, time::T, value::Float64)
    @assert time >= feature.lastTime
    feature.lastTime = time
    unshift!(feature.bufferTimes, time)
    unshift!(feature.bufferValues, value)
end
function valueat{T<:TimeType}(feature::LagLastValueFeature{T}, time::T)
    @assert time >= feature.lastTime

    # move data from the buffer until we get to the current time point
    for i in length(feature.bufferTimes):-1:1
        if time - feature.bufferTimes[i] >= feature.lag
            pop!(feature.bufferTimes)
            feature.value = pop!(feature.bufferValues)
        else
            break
        end
    end

    feature.value
end

end # module