text:
function add_point_to_map!(m::MapData, point::Tuple{Float64,Float64}, external_id::Int64 = -999, verbose::Bool = false)
    """
    Description:
        Converts a pair Latitude-Longitude of coordinates `point` to a node on a map `m`.
    Returns:
        It updates m.n and m.nodes and m.v.
    """
    point_LLA = OpenStreetMapX.LLA(point[1],point[2])
    add_point_to_map!(m, point_LLA, external_id, verbose)
end

function add_point_to_map!(m::MapData, point::OpenStreetMapX.LLA, external_id::Int64 = -999, verbose::Bool = false)
    """
    Description:
        Converts a pair Latitude-Longitude of coordinates `point` to a node on a map `m`.
    Returns:
        It updates m.n and m.nodes and m.v.
    """
    point_ENU = OpenStreetMapX.ENU(point, m.bounds)

    # generate node id for the point
    unq = false
    new_node_id = 0
    while unq == false
        new_node_id = -abs(rand(Int32))
        unq = ifelse(new_node_id ∉  m.n, true, false)
    end

    # Add to nodes
    push!(m.n, new_node_id)
    m.nodes[new_node_id] = point_ENU

    # Add to vertices
    m.v[new_node_id] = m.v.count + 1

    # Add to graph
    add_vertex!(m.g)

    if verbose == true
        println("New node and vertice created.")
        println("Node id: ", new_node_id)
        println("Vertice id: ", m.v.count)
    end

    if external_id != -999
        return new_node_id => external_id
    end
end

function find_route_id_by_stop_id(stop_id::Int64)
    """
    Description:
        Gives route_id to which stop_id is assigned.
    Returns:
        route_id which is associated with given stop_id.
    """
    for (nr,route) in enumerate(keys.(values(routes_dict)))
        if stop_id in route
           return collect(keys(routes_dict))[nr]
        end
    end
end

function describe_node(node_id::Int)
	"""
	Takes node_id from the final graph and if this node is a TTC node the description will be given.
	"""
    stop_my_id = m.v[node_id] - orig_map_nodes_num;
    stop_org_id = trip_longest_df[coalesce.(trip_longest_df[!, "my_stop_id"] .== stop_my_id, false),:][!, "stop_id"][1]
    stop_name = stops[coalesce.(stops[!, "stop_id"] .== stop_org_id, false),:][!, "stop_name"][1]
    route_id = find_route_id_by_stop_id(stop_my_id)
    route_name = routes_filtered[coalesce.(routes_filtered[!, "route_id"] .== parse(Int64,split(route_id,"_")[1]), false),:][!, "route_long_name"][1]
    route_number = routes_filtered[coalesce.(routes_filtered[!, "route_id"] .== parse(Int64,split(route_id,"_")[1]), false),:][!, "route_short_name"][1]
    println("It is stop '$stop_name' on the route '$route_name ($route_number) [$route_id]'.")
end

function nearest_1_node_allowed_list(nodes::Dict{Int,T}, loc::T, node_list::Array{Int64,1}, max_dist::Float64) where T<:(Union{OpenStreetMapX.ENU,OpenStreetMapX.ECEF})
	"""
	Find the nearest one node to a given location `loc`.
	"""
    ind_dist = OrderedDict{Int64, Float64}()

    for ind in node_list
        dist = OpenStreetMapX.distance(nodes[ind], loc)
        if dist <= max_dist
            ind_dist[ind] = dist
        end
    end

    dict_sorted_by_distance = sort(ind_dist; byvalue=true)
    keys_sorted_dict_by_distance = collect(keys(dict_sorted_by_distance))

    selected_ind = Dict{Int64,Float64}()

	@assert length(dict_sorted_by_distance) > 0 "There is no even one point!"
    selected_ind[keys_sorted_dict_by_distance[1]] = dict_sorted_by_distance[keys_sorted_dict_by_distance[1]]

    return selected_ind
end

function nearest_n_nodes_allowed_list(nodes::Dict{Int,T}, loc::T, node_list::Array{Int64,1}, n_nearest::Int64, max_dist::Float64) where T<:(Union{OpenStreetMapX.ENU,OpenStreetMapX.ECEF})
    """
	Find the nearest N nodes to a given location `loc`
	"""
	ind_dist = OrderedDict{Int64, Float64}()

    for ind in node_list
        dist = OpenStreetMapX.distance(nodes[ind], loc)
        if dist <= max_dist
            ind_dist[ind] = dist
        end
    end

    keys_sorted_dict_by_distance = collect(keys(sort(ind_dist; byvalue=true)))
    selected_ind = Int64[]
    @assert length(keys_sorted_dict_by_distance) > 0 "There is no even one point!"
	k = 0
	while k < n_nearest
	   push!(selected_ind,keys_sorted_dict_by_distance[k+1])
	   k += 1
	end

    return selected_ind
end

function nodes_within_range_with_allowed_list(nodes::Dict{Int,T}, loc::T, node_list::Array{Int64,1}, range::Float64 = 100.0) where T<:(Union{OpenStreetMapX.ENU,OpenStreetMapX.ECEF})
    """
	Find nodes within some range to a given location. Only nodes in the `node_list` will be checked.
	"""
	if range == Inf
        return node_list
    end
    indices = Int[]
    for ind in node_list
        dist = OpenStreetMapX.distance(nodes[ind], loc)
        if dist < range
            push!(indices, ind)
        end
    end
    return indices
end

## Examples ##
# nodes_within_range_with_allowed_list(m.nodes, m.nodes[394502528], m.n,  10.0)
# nearest_n_nodes_allowed_list(m.nodes, m.nodes[394502528], [i for i in m.n if i ∉ collect(keys(node_id_stop_id_dict))], 1, Inf )
# nearest_1_node_allowed_list(m.nodes, m.nodes[394502528], [i for i in m.n if i ∉ collect(keys(node_id_stop_id_dict))], Inf )

function add_edges_and_weights_to_map!(;m::MapData
        , TTC_graph::SimpleWeightedDiGraph{Int64,Float64}
        , stop_id_to_node_id_dict::Dict{Int64,Int32}
        , node_id_stop_id_dict::Dict{Int32,Int64}
        , wait_time::Int64
        , orig_map_nodes_num::Int64
		, agents_speed_in_sec_per_m::Float64)
    """
    Description:
        Adds information from the TTC graph to the main city graph (e.g. central_toronto).
		Edges are connections between stops.
		Weights are numbers of seconds between two stops.

		NOTE! Assumption: Healthy human walks at a pace of 100m in 80 sec.
		NOTE! Assumption: Each TTC node must be connected at least with one city node.In general each TTC node is connected with all city nodes within 100m range.
		NOTE! Assumption: Each agent assume that expected waiting time is half of TTC frequency. For example, if for a TTC car A it takes 10m to get from stop X to stop Y an agent will assume that he/she will wait for 5m at stop Y.
    Returns:
        It modifies the map m.
    """

    edges = adjacency_matrix(TTC_graph)

    matrix_rows = [i[1] for i in findall(!iszero, m.w)];
    matrix_cols = [i[2] for i in findall(!iszero, m.w)];
    matrix_weights = [i * agents_speed_in_sec_per_m for i in m.w.nzval]; # The average, healthy human walks at a pace of 100m in 80 sec. Here we assume 100m=80s, so 1m = 0.8s.

    n_added = 0

    for edge in findall(!iszero, edges) # only non-zero edges
        egde_node_a = stop_id_to_node_id_dict[edge[1]]
        edge_node_b = stop_id_to_node_id_dict[edge[2]]

        # NOTE! edges contain also reversed nodes ==> e.g. (a --> b) and (b --> a) with the same weights.
        # That us why we do push only once. Our TTC_graph already contains reversed connections.

        # add new edges to m.e
        push!(m.e, (egde_node_a, edge_node_b)) # m.e is an object without any weights

        # add new edges to m.w
        push!(matrix_rows, m.v[egde_node_a])
        push!(matrix_cols, m.v[edge_node_b])

        # add new weights to m.w
        push!(matrix_weights, edges[(edge[1],edge[2])])

        # add edges to m.g
        add_edge!(m.g, (m.v[egde_node_a],m.v[edge_node_b]))

        # add classes to m.class
        push!(m.class, 6)

        n_added += 1
    end

    # Integrating TTC and map
    # In purpose to integrate TTC map and OpenStreetMap we decided to connect TTC nodes with all not TTC nodes, which
    # are in 100m range. If there is no nodes in 100m, we then connect TTC node to the nearest city node.

    for ttc_node in  collect(keys(node_id_stop_id_dict))
        # [i for i in m.n if i > 0] --> it helps us to select only non-TTC node keys.
        ttc_neigbours_nodes = nodes_within_range_with_allowed_list(m.nodes, m.nodes[ttc_node], [i for i in m.n if i ∉ collect(keys(node_id_stop_id_dict))],  100.0)
        time_from_ttc = agents_speed_in_sec_per_m * 100 #80 sec to go 100m.
        if length(ttc_neigbours_nodes) == 0
            # if there is no nodes in a range of 100m, then we connect TTC node to the nearest one city node (both directions!).
            node_dist = nearest_1_node_allowed_list(m.nodes, m.nodes[ttc_node], [i for i in m.n if i ∉ collect(keys(node_id_stop_id_dict))], Inf)
            ttc_neigbours_nodes = collect(keys(node_dist))[1]
            time_from_ttc = collect(values(node_dist))[1] * agents_speed_in_sec_per_m
        end
        for neigbour_node in ttc_neigbours_nodes
            # NOTE! Here we need to make connections for both sides ==> (a --> b) and (b --> a). That us why we push two times.

            # add new edges to m.e
            push!(m.e, (ttc_node, neigbour_node))
            push!(m.e, (neigbour_node, ttc_node))

            # add new edges to m.w
            push!(matrix_rows, m.v[ttc_node])
            push!(matrix_cols, m.v[neigbour_node])

            push!(matrix_rows, m.v[neigbour_node])
            push!(matrix_cols, m.v[ttc_node])

            # add new weights to m.w
            push!(matrix_weights, time_from_ttc)
            push!(matrix_weights, (time_from_ttc + wait_time)) # To every arrow from city to TTC we need to add WAIT_TIME as time of waiting for TTC.

            # add edges to m.g
            add_edge!(m.g, (m.v[ttc_node], m.v[neigbour_node]))
            add_edge!(m.g, (m.v[neigbour_node], m.v[ttc_node]))

            # add classes to m.class
            push!(m.class, 6)
            push!(m.class, 6)

            n_added += 2

        end
    end

    # New weights structure for the map
    m.w = SparseArrays.sparse(matrix_rows, matrix_cols, matrix_weights, m.v.count, m.v.count)

    # check_correctness_of_map_udpate
	@assert adjacency_matrix(m.g)[1+orig_map_nodes_num,2+orig_map_nodes_num] > 0.0 "SOMETHING WENT WRONG 1+orig_map_nodes_num-->2+orig_map_nodes_num should have a non-zero weight!"
	@assert adjacency_matrix(m.g)[2+orig_map_nodes_num,1+orig_map_nodes_num] == 0.0 "SOMETHING WENT WRONG 2+orig_map_nodes_num-->1+orig_map_nodes_num should have zero weight!"
	@assert adjacency_matrix(m.g)[616+orig_map_nodes_num,615+orig_map_nodes_num] > 0.0 "SOMETHING WENT WRONG 616+orig_map_nodes_num-->615+orig_map_nodes_num should have a non-zero weight!"
    @assert adjacency_matrix(m.g)[615+orig_map_nodes_num,616+orig_map_nodes_num] == 0.0 "SOMETHING WENT WRONG 615+orig_map_nodes_num-->616+orig_map_nodes_num should have zero weight!"

	println("Expected TTC car waiting time: ", wait_time, " sec")
	println("Assumed agent`s speend in seconds per m: ", agents_speed_in_sec_per_m, "sec/meter")
    println("TTC info was added to the map.Total of $n_added edges and weights added to the map")
end

function adding_TTC_points_to_map!(routes_dict::OrderedDict{Symbol, OrderedDict{Int,Tuple{Float64,Float64}}}, m::MapData)
	"""
	Description:
		Adds TTC stops as vertices to the map graph (e.g. central_toronto.osm).
		It also keeps the information about added points in the dictionaries.
	Returns:
		1) Modifies the map
		2) node_id --> stop_id and stop_id --> node_id dictionaries
	"""
    node_id_stop_id_dict = Dict{Int32, Int64}()
    for route in collect(values(routes_dict))
        for (stop_id, stop_LLA) in sort(route)
           new_point_dict = add_point_to_map!(m, stop_LLA, stop_id, false)
           push!(node_id_stop_id_dict, new_point_dict)
        end
    end
    # reversing the dict so we map previous node keys to map node keys
    stop_id_to_node_id_dict= Dict(value => key for (key, value) in node_id_stop_id_dict);
    println("TTC POINTS WERE ADDED TO THE MAP")
    return node_id_stop_id_dict, stop_id_to_node_id_dict
end

function create_routes_dicts(routes_dict, secs_between_stops, orig_map_nodes_num)
	"""
	Description:
		Stops on every route which was added to the main map (e.g. central_toronto.osm) and seconds it takes to travel between two stops on the route.
		orig_map_nodes - number of nodes in the original graph (e.g. central_toronto.osm).
		NOTE! TTC nodes are added directly after original nodes, so that is why we need (stop_id + orig_map_nodes_num).
	Returns:
		Two dictionaries, where one shows all stops on each route and another one shows the time it takes to get from one stop to another on each route.
	"""
    # All points on a road
    routes_path = OrderedDict{Symbol, OrderedDict{Int64,Int64}}()
    routes_distances = OrderedDict{Symbol, OrderedDict{Int64,Float64}}()
    for route in keys(routes_dict)
        stops_sel = OrderedDict{Int64,Int64}()
        distances  = OrderedDict{Int64,Float64}()
        for stop in 1:size(collect(keys(routes_dict[route])))[1]-1
            stop_a = collect(keys(routes_dict[route]))[stop]
            stop_b = collect(keys(routes_dict[route]))[stop+1]
            stops_sel[stop_a+orig_map_nodes_num] = stop_b+orig_map_nodes_num
            distances[stop_a+orig_map_nodes_num] = secs_between_stops[(stop_a,stop_b)]
        end
        routes_path[route] = stops_sel
        routes_distances[route] = distances
    end
    return routes_path, routes_distances
end

println("Complete map with TTC functions were read.")

max_stars_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_issues_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_forks_repo_name: NykPol/EpidemicInUrbanNetworkToronto


text:
#   Unit "libSVM.jl" of the PosDefManifoldML Package for Julia language
#   v 0.2.1 - last update 18th of October 2019
#
#   MIT License
#   Copyright (c) 2019,
#   Anton Andreev, CNRS, Grenoble, France:

# ? CONTENTS :
#   This unit implements a wrapper to libSVM. It projects data to tangent space
#   and it applies SVM classification using Julia's SVM wrapper.

abstract type SVMmodel<:TSmodel end

"""
**Examples**:
```
# Note: creating models with the default creator is possible,
# but not useful in general.

using PosDefManifoldML

# generate data
PTr, PTe, yTr, yTe=gen2ClassData(10, 30, 40, 60, 80)

# create and train an SVM model with default parameters for tangent space calculation and SVM
model=fit(SVM(), PTr, yTr)

# predict using this model
yPred=predict(model, PTe, :l)

# calculate prediction error
predictErr(yTe, yPred)

You can supply parameters for both tangent space calculaton and SVM:

model=fit(SVM(), PTr, yTr, cost = 1.0)

```
"""

mutable struct SVM <: SVMmodel
    	metric        :: Metric
		svmtype       :: Type
		kernel        :: Kernel.KERNEL
		meanISR
		featDim
		# LIBSVM args
		epsilon
		cost
		gamma
		svmModel #used to store the training model from the SVM library
    function SVM( metric :: Metric=Fisher;
				  svmtype = SVC,
                  kernel  = Kernel.RadialBasis,
				  meanISR = nothing,
				  featDim = nothing,
				  epsilon = nothing,
				  cost    = nothing,
				  gamma   = nothing,
				  svmModel = nothing)
	   	 			 new(metric, svmtype, kernel, meanISR, featDim,
					     epsilon, cost, gamma, svmModel)
    end
end

"""
```
mutable struct SVM <: SVMmodel
    	metric        :: Metric
		svmtype       :: Type
		kernel        :: Kernel.KERNEL
		meanISR
		epsilon
		cost
		gamma
		svmModel
end
```

ATTENTION: This class is not used to set the parameters, just to store them. You need to set the parameterts in the FIT function.
Example:
```
model=fit(SVM(), PTr, yTr, cost = 1.0)
```

`.metric`, of type
[Metric](https://marco-congedo.github.io/PosDefManifold.jl/dev/MainModule/#Metric::Enumerated-type-1),
is the metric that will be adopted to compute the mean used
as base-point for tangent space projection. By default the
Fisher metric is adopted. See [mdm.jl](@ref)
for the available metrics. If the data used to train the model
are not positive definite matrices, but Euclidean feature vectors,
the `.metric` field has no use.

`.meanISR` is optionally passed to the [`fit`](@ref)
function. By default it is computed thereby.
If the data used to train the model
are not positive definite matrices, but Euclidean feature vectors,
the `.meanISR` field has no use and is set to `nothing`.

The following are parameters that are passed to the LIBSVM package:

`svmtype::Type=SVC`: Type of SVM to train `SVC` (for C-SVM), `NuSVC`
    `OneClassSVM`, `EpsilonSVR` or `NuSVR`. Defaults to `OneClassSVM` if
    `y` is not used

`kernel::Kernels.KERNEL=Kernel.RadialBasis`: Model kernel `Linear`, `polynomial`,
    `RadialBasis`, `Sigmoid` or `Precomputed`

`gamma::Float64=1.0/size(X, 1)` : γ for kernels

`cost::Float64=1.0`: cost parameter C of C-SVC, epsilon-SVR, and nu-SVR

`epsilon::Float64=0.1`: epsilon in loss function of epsilon-SVR
"""

function fit(model  :: SVMmodel,
               𝐏Tr  :: Union{ℍVector, Matrix{Float64}},
               yTr  :: IntVector;
		   # Tnagent space parameters
		   w        :: Union{Symbol, Tuple, Vector} = [],
           meanISR  :: Union{ℍ, Nothing} = nothing,
		   vecRange :: UnitRange = 𝐏Tr isa ℍVector ? (1:size(𝐏Tr[1], 2)) : (1:size(𝐏Tr, 2)),
		   # SVM paramters
		   svmtype  :: Type = SVC,
		   kernel   :: Kernel.KERNEL = Kernel.RadialBasis,
		   epsilon  :: Float64 = 0.001,
		   cost     :: Float64 = 1.0,
		   gamma    :: Float64 = 1/_getDim(𝐏Tr, vecRange),
		   # Generic parametes
           verbose  :: Bool = true,
		         ⏩  :: Bool = true,
          parallel  :: Bool=false)

    #println(defaultFont, "Start")
    ⌚=now() # get the time in ms

    # output model
    ℳ=deepcopy(model)

    # checks
    # 𝐏Tr isa ℍVector ? nObs=length(𝐏Tr) : nObs=size(𝐏Tr, 1)

	# check w argument and get weights
    (w=_getWeights(w, yTr, "fit ("*_modelStr(ℳ)*" model)")) == nothing && return

    # projection onto the tangent space
    if 𝐏Tr isa ℍVector
        verbose && println(greyFont, "Projecting data onto the tangent space...")
        if meanISR==nothing
            (X, G⁻½)=tsMap(ℳ.metric, 𝐏Tr; w=w, vecRange=vecRange, ⏩=⏩)
            ℳ.meanISR = G⁻½
        else
            X=tsMap(ℳ.metric, 𝐏Tr; w=w, vecRange=vecRange, meanISR=meanISR, ⏩=⏩)
            ℳ.meanISR = meanISR
        end
    else
        X=𝐏Tr
    end

	ℳ.svmtype = svmtype
	ℳ.kernel = kernel
    ℳ.gamma = gamma
	ℳ.epsilon = epsilon
	ℳ.cost = cost
	ℳ.featDim = size(X, 2)

	#convert data to LIBSVM format; first dim is features, second dim is observations
	instances = X'

    verbose && println(defaultFont, "Fitting SVM model...")
    model = svmtrain(instances, yTr; svmtype = ℳ.svmtype, kernel = ℳ.kernel, epsilon = ℳ.epsilon, cost=ℳ.cost, gamma = ℳ.gamma);

    ℳ.svmModel = model

    verbose && println(defaultFont, "Done in ", now()-⌚,".")
    return ℳ
end

function predict(model   :: SVMmodel,
                 𝐏Te     :: Union{ℍVector, Matrix{Float64}},
                 what    :: Symbol = :labels,
                vecRange :: UnitRange = 𝐏Te isa ℍVector ? (1:size(𝐏Te[1], 2)) : (1:size(𝐏Te, 2)),
                 checks  :: Bool = true,
                 verbose :: Bool = true,
                  ⏩     :: Bool = true)

    ⌚=now()

    # checks
    if checks
        if !_whatIsValid(what, "predict ("*_modelStr(model)*")") return end
    end

    # projection onto the tangent space
    if 𝐏Te isa ℍVector
        verbose && println(greyFont, "Projecting data onto the tangent space...")
        X=tsMap(model.metric, 𝐏Te; meanISR=model.meanISR, ⏩=⏩, vecRange=vecRange)
    else X=𝐏Te[:, vecRange] end

    #convert data to LIBSVM format first dim is features, second dim is observations
    instances = X'

	(predicted_labels, decision_values) = svmpredict(model.svmModel, instances;)
    🃏 = predicted_labels

    verbose && println(defaultFont, "Done in ", now()-⌚,".")
    verbose && println(titleFont, "\nPredicted ",_what2Str(what),":", defaultFont)
    return 🃏
end

# ++++++++++++++++++++  Show override  +++++++++++++++++++ # (REPL output)
function Base.show(io::IO, ::MIME{Symbol("text/plain")}, M::SVM)
    println(io, titleFont, "\n↯ SVM LIBSVM machine learning model")
    println(io, defaultFont, "  ", _modelStr(M))
    println(io, separatorFont, "⭒  ⭒    ⭒       ⭒          ⭒", defaultFont)
    println(io, "type    : PD Tangent Space model")
    println(io, "features: tangent vectors of length $(M.featDim)")
    println(io, "classes : 2")
    println(io, "fields  : ")
	println(io, separatorFont," .metric      ", defaultFont, string(M.metric))

	if 		M.svmtype==SVC s="SVC"
	elseif  M.svmtype==C-SVM s="C-SVM"
	elseif  M.svmtype==EpsilonSVR s="EpsilonSVR"
	elseif  M.svmtype==OneClassSVM s="OneClassSVM"
	elseif  M.svmtype==NuSVR s="NuSVR"
	else    s = "Warning: the SVM type is unknown"
	end
	println(io, separatorFont," .svmtype     ", defaultFont, "$s")

	println(io, separatorFont," .kernel      ", defaultFont, "$(string(M.kernel))")

	#`kernel::Kernels.KERNEL=Kernel.RadialBasis`: Model kernel `Linear`, `polynomial`,
	#    `RadialBasis`, `Sigmoid` or `Precomputed`

	if M.meanISR == nothing
        println(io, greyFont, " .meanISR      not created")
    else
        n=size(M.meanISR, 1)
        println(io, separatorFont," .meanISR     ", defaultFont, "$(n)x$(n) Hermitian matrix")
    end

	M.epsilon==nothing ? println(io, "       not created ") :
						 println(io, separatorFont," .epsilon     ", defaultFont, "$(round(M.epsilon, digits=9))")

    M.cost==nothing ?    println(io, "       not created ") :
						 println(io, separatorFont," .cost        ", defaultFont, "$(round(M.cost, digits=3))")

    M.gamma==nothing ?   println(io, "       not created ") :
 						 println(io, separatorFont," .gamma       ", defaultFont, "$(round(M.gamma, digits=5))")

	println(io, separatorFont," .svmModel ", defaultFont, "   (LIBSVM model)")
end

ext: jl


text:
using RegressionAndOtherStories

helicopters = CSV.read(ros_data("Helicopters", "helicopters.txt"),
    DataFrame; delim=" ")
f = open(ros_datadir("Helicopters", "helicopters.csv"), "w")
CSV.write(f, helicopters[:, 1:4])
close(f)

helicopters2 = CSV.read(ros_datadir("Helicopters", "helicopters.csv"), DataFrame)


text:
function count(context::NamedTuple,selections::NamedTuple)
    request = context.request
	try
		response = HTTP.Response(200)
        HTTP.Messages.setheader(response,"Content-type" => "application/json")

        body = IOBuffer(response.body, write=true)
        t = @elapsed write( body, JSON.json(count(selections)) )

		@info """$(request.method) $(replace(request.target, r"\?.+" => "")) | $t seconds"""
		return response

	catch exception
		printstyled("ERROR $(replace(request.target, r"\?.+" => "")) | $exception\n",color=:red)
		return HTTP.Response(500)
	end
end

function count(selections::NamedTuple)
	return [ (names = decode(code, selections.names.keys ), count = count) for (code, count) ∈ countmap(selections.codes[selections.rows]) ]
end

max_stars_repo_name: gszep/flow-explorer

max_issues_repo_name: gszep/flow-explorer

max_forks_repo_name: gszep/flow-explorer


text:
let

    METHODS = [
        :FBA_MAX_Z_MIN_COST,
        :FBA_Z_FIX_MAX_VG_MIN_COST,
        # :FBA_Z_FIX_MAX_VG_MAX_COST,
        :FBA_Z_FIX_MIN_VG_MIN_COST,
        # :FBA_Z_FIX_MIN_VG_MAX_COST,
        # :FBA_Z_FIX_ADJ_MAX_VG_MIN_COST,
        # :FBA_Z_FIX_ADJ_MAX_VG_MAX_COST,
        # :FBA_Z_FIX_ADJ_MIN_VG_MIN_COST,
        # :FBA_Z_FIX_ADJ_MIN_VG_MAX_COST,
        # :FBA_Z_FIX_MAX_ATP_MIN_COST,
        :FBA_Z_FIX_MAX_ATP_MAX_COST,
        :ME_MAX_POL,
        # :ME_MAX_POL_B0
    ]

    fontsize = 10
    commom_params = (;
        # dpi = 1000,
        # thickness_scaling = 1.0,
        # xlabel = "exp",
        # ylabel = "model",
        titlefont = 15,
        # aspectatio = 1.0,
        # xguidefontsize = fontsize,
        # yguidefontsize = fontsize
        # xaxis = false,
        # yaxis = false
    )

    avDs = [:exch, 0.05, 0.1, 0.2, 0.3, 0.4, :tot]
    ps = Plots.Plot[]
    for avD in avDs, method in METHODS
        title = string(method, "\n ------- \nD = ", avD)
        if avD == :exch
            k = (method, :exchs_flx_cor, :all_iders)
            p = plot!(deepcopy(PS[k...]); title, commom_params...)
        else
            # corrs
            k = (method, :inner_flx_cor, :all_subs, avD)
            p = plot!(deepcopy(PS[k...]); title, commom_params...)
        end
        push!(ps, p)
    end
    layout = (length(avDs), length(METHODS))
    mysavefig(ps, "inner_corr_callage"; layout)
    return

    # comparizon
    # push!(ps, PS[:inner_flx_comparizon, :all_subs, avD])

    # MSE
    push!(ps, PS[:inner_flx_MSE, :all_subs, avD])

    layout = (2, 4)
    margin = -3
    mysavefig(ps, "inner_callage"; avD, layout, margin)
end


text:
__precompile__()
module TKPappKitchenSink

using ArgCheck
using DocStringExtensions
using PGFPlotsX
using StaticArrays

include("utilities.jl")
include("math.jl")
include("plots.jl")

end # module


text:
# import Base.copy

export
    Point,
    set!

mutable struct Point{T <: AbstractFloat}
    x::T
    y::T

    function Point{T}() where {T <: AbstractFloat}
        new(0.0, 0.0)
    end

    function Point{T}(x::T, y::T) where {T <: AbstractFloat}
        new(x, y)
    end
end

# copy specializations
copy(p::Point) = Point{Float64}(p.x, p.y)

# setters/getters set_point
function set!(p::Point{T}, x::T, y::T) where {T <: AbstractFloat}
    p.x = x
    p.y = y;
end

function set!(p::Point{T}, from::Point{T}) where {T <: AbstractFloat}
    p.x = from.x
    p.y = from.y;
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
# This is based on the runtests.jl file of MPI.jl.

using InteractiveUtils
using MPI: MPI, mpiexec

# Load test packages to trigger precompilation
using PencilArrays

# These tests can be run in serial mode
include("permutations.jl")

test_files = [
    "reductions.jl",
    "broadcast.jl",
    "arrays.jl",
    "pencils.jl",
    "array_interface.jl",
    "io.jl",
]

Nproc = let N = get(ENV, "JULIA_MPI_TEST_NPROCS", nothing)
    N === nothing ? clamp(Sys.CPU_THREADS, 4, 6) : parse(Int, N)
end

println()
versioninfo()
println("\n", MPI.MPI_LIBRARY_VERSION_STRING, "\n")

for fname in test_files
    @info "Running $fname with $Nproc processes..."
    mpiexec() do cmd
        run(`$cmd -n $Nproc $(Base.julia_cmd()) $fname`)
    end
    println()
end


text:
function get_growth_rate(file, meta_data_file, stan_file)
    df = read_OD(file, meta_data_file)

    s = open(stan_file) do file
        read(file, String)
    end
    growth = Array{Float64, 1}[]
    rate = Array{Float64, 1}[]
    bw_model = Stanmodel(
      name="bw",
      nchains=4,
      num_warmup=1000,
      num_samples=1000,
      thin=10,
      model=s,
      printsummary=false
    )
    max_rates = zeros(Float64, 48)
    for (j, well) in enumerate(unique(df[!, :well]))
        t = df[df.well .== well, Symbol("time_[s]")]
        N = df[df.well .== well, :OD]
        strain = df[df.well .== well, :strain] |> unique
        if strain[1] != "blank"
            bw_data = Dict(
                "N"=>length(t),
                "x"=>t,
                "y"=>N,
                "N_predict"=>72,
                "x_predict"=>range(minimum(t), stop=maximum(t), length=72),
            )

            _, bw_chains, bw_names = stan(bw_model, bw_data, summary=false)

            d = [collect_params_from_chain(bw_names, bw_chains[:,:,i]) for i in 1:4]
            chain_maximums = map(x->maximum(mean(x["g_predict"], dims=1)), d)
            max_rates[j] = maximum(chain_maximums)
            println(string("Well ", well, " done"))

            chain_means_growth = map(x->mean(x["y_predict"], dims=1),d)
            mean_growths = mean.([[chain_means_growth[i][j] for i in 1:4] for j in 1:72])
            chain_means_rates = map(x->mean(x["g_predict"], dims=1),d)
            mean_rates = mean.([[chain_means_rates[i][j] for i in 1:4] for j in 1:72])

            push!(growth, mean_growths)
            push!(rate, mean_rates)

        end
    end
    return max_rates, growth, rate
end


text:
@testset "generate indices" begin
    @test DeconvOptim.generate_indices(5, 2, 1, 5) == (Any[:i1, :(i2 + 1), :i3, :i4, :i5], Any[:i1, :(i2 + 5), :i3, :i4, :i5])
end

@testset "generate_laplace" begin
    x = DeconvOptim.generate_laplace(2, [1, 2], [4 , 5], debug=true)
    @test x==Any[:(res = abs2((4 * arr[i1 + 1, i2] + 4 * arr[i1 + -1, i2]) + (5 * arr[i1, i2 + 1] + 5 * arr[i1, i2 + -1]) + -(18* arr[i1, i2])))]
    x = DeconvOptim.generate_laplace(2, [1, 2], [1 , 1], debug=true)
    @test x==Any[:(res = abs2((1 * arr[i1 + 1, i2] + 1 * arr[i1 + -1, i2]) + (1 * arr[i1, i2 + 1] + 1 * arr[i1, i2 + -1]) + -(4 * arr[i1, i2])))]

end

@testset "Tikhonov" begin
    x = [1,2,3,1,3,1,12.0,2,2,3,2.0]
    reg = Tikhonov(num_dims=1, sum_dims=[1], weights=[1])
    @test 756 ≈ reg(x)

    reg = Tikhonov(num_dims=1, mode="spatial_grad_square")
    @test 188 ≈ reg(x)

    reg = Tikhonov(num_dims=1, mode="identity")
    @test 190 ≈ reg(x)

end

@testset "Good's roughness" begin

    x = generate_GR(5, [1,2], [4, 5], 1, -1, debug=true)

    @test x == Any[:(res = -2.0 * arr[i1, i2, i3, i4, i5] * (4 * (arr[i1 + 1, i2, i3, i4, i5] + arr[i1 + -1, i2, i3, i4, i5]) + 5 * (arr[i1, i2 + 1, i3, i4, i5] + arr[i1, i2 + -1, i3, i4, i5]) + -18 * arr[i1, i2, i3, i4, i5]))]

    x = [1,2,3,1,3,1,12.0,2,2,3,2.0]
    reg = GR(num_dims=1, sum_dims=[1], weights=[1])
    @test 22.71233466779126 ≈ reg(x)

end

@testset "TV" begin
    x = [1,2,3,1,3,1,12.0,2,2,3,2.0]
    reg = TV(num_dims=1, sum_dims=[1], weights=[1])
    @test 31.00010002845424 ≈ reg(x)

    @test TV_cuda(num_dims=2)(x) ≈ reg(x)
    @test TV_cuda(num_dims=3)(x) ≈ reg(x)

    x = generate_TV(4, [1,2], [5, 7], 1, -1, debug=true)
    @test x == Any[:(res = sqrt(5 * abs2(arr[i1 + 1, i2, i3, i4] - arr[i1 + -1, i2, i3, i4]) + 7 * abs2(arr[i1, i2 + 1, i3, i4] - arr[i1, i2 + -1, i3, i4]) + 1.0f-8))]

end


text:
function J1d(J::Matrix{Float64},Q::Float64,A::Float64,beta::Float64,rho::Float64,ts::Float64,
    xs::Float64,Qs::Float64,As::Float64)
    if A < 0.
        error("A negative. A = $A.")
    end
    J[1,2] = ts*Qs/(As*xs);
    J[2,1] = -Qs*ts/(As*xs)*Q^2/(A^2) + As^1.5*ts/(Qs*xs)*beta*sqrt(A)/(2.*rho);
    J[2,2] = Qs*ts/(As*xs)*2.*Q/A;
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
using Test
using CommutativeRings

@testset "generic" begin include("generic.jl") end
@testset "typevars" begin include("typevars.jl") end
@testset "ZZ" begin include("zz.jl") end
@testset "QQ" begin include("qq.jl") end
@testset "ZZmod" begin include("zzmod.jl") end
@testset "univarpoly" begin include("univarpolynom.jl") end
@testset "multivarpoly" begin include("multivarpolynom.jl") end
@testset "ideal" begin include("ideal.jl") end
@testset "fraction" begin include("fraction.jl") end
@testset "quotient" begin include("quotient.jl") end
@testset "factorization" begin include("factorization.jl") end
@testset "intfactors" begin include("intfactorization.jl") end
@testset "galoisfields" begin include("galoisfields.jl") end
@testset "numbertheory" begin include("numbertheoretical.jl") end
@testset "enumerations" begin include("enumerations.jl") end
@testset "linearalgebra" begin include("linearalgebra.jl") end
#@testset "gcd" begin include("gcd.jl") end


text:
module test_julia_abstractset

using Test

@test union((1,3,5,5,5)) == union([1,3,5,5,5]) == [1,3,5]

end # module test_julia_abstractset


text:
baremodule ConcurrentCollections

export
    #
    ConcurrentDict,
    ConcurrentQueue,
    ConcurrentStack,
    Delete,
    DualLinkedConcurrentRingQueue,
    DualLinkedQueue,
    Keep,
    LinkedConcurrentRingQueue,
    WorkStealingDeque,
    length_upper_bound,
    length_upper_bound,
    modify!,
    maybeget,
    maybepop!,
    maybepopfirst!

import Base

struct Keep{Value}
    value::Value
end

struct Delete{Value}
    value::Value
end

abstract type ConcurrentDict{Key,Value} <: Base.AbstractDict{Key,Value} end

function modify! end
function maybepop! end
function maybepopfirst! end
function maybeget end
function length_lower_bound end
function length_upper_bound end

module Implementations

using Base:
    #
    HasEltype,
    IteratorEltype,
    RefValue,
    aligned_sizeof

using ..ConcurrentCollections:
    ConcurrentCollections,
    ConcurrentDict,
    Delete,
    Keep,
    length_lower_bound,
    length_upper_bound,
    modify!,
    maybeget,
    maybepop!,
    maybepopfirst!

include("UnsafeAtomics.jl")
using .UnsafeAtomics: acq_rel, acquire, monotonic, release, seq_cst, unordered

include("utils.jl")
include("cache.jl")
include("atomicsutils.jl")
include("promise.jl")
include("dict.jl")
include("workstealing.jl")
include("msqueue.jl")
include("stack.jl")
include("lcrq.jl")
include("dlcrq.jl")
include("ssqueue.jl")
include("misc.jl")

end  # module Implementations

using .Implementations:
    ConcurrentQueue,
    ConcurrentStack,
    DualLinkedConcurrentRingQueue,
    DualLinkedQueue,
    LinkedConcurrentRingQueue,
    WorkStealingDeque

Implementations.define_docstrings()

end  # baremodule ConcurrentCollections


text:
@testset "typing" begin
    g = rand(Typing_sampler(5, 10), 7)
    @test g isa Vector{Vector{Vector{Int}}}
    @test length(g) == 7
    @test length(g[3]) == 5
    @test Set(vcat(vcat.(g...)...)) == Set(1:10)

    g = rand(Typing_sampler(3, .4, [.1, .2, .3], ['a', 'b', 'c']), 30)
    @test g isa Vector{Vector{Vector{Char}}}
    @test length(g) == 30
    @test all(length.(g) .== 3)
    @test Set(vcat(vcat.(g...)...)) == Set("abc")
end

max_stars_repo_name: LilithHafner/FBD

max_issues_repo_name: LilithHafner/FBD

max_forks_repo_name: LilithHafner/FBD


text:
module Constants

export μ₀, ϵ₀, c₀, η₀, twopi, fourpi, tol, dbmin, tdigits

"Permeability of free space [H/m]"
const μ₀ = 1.25663706212*1e-6

"Speed of light [m/s]"
const c₀ = 299792458.0

"Permeability of free space [F/m]"
const ϵ₀ = 1/(μ₀*c₀^2)

"Intrinsic impedance of free space [F/m]"
const η₀ = sqrt(μ₀/ϵ₀)

const twopi = 2π
const fourpi = 4π
const tol = 1e-4

"Minimum modal attenuation"
const dbmin = 30.0

"Min. elect. length for a layer to not be included in a `GBLOCK`"
const min_elength = twopi * 0.1

"Number of digits to use in printing elapsed time"
const tdigits = 4

end


text:
####
# UI Oracle tests
####
facts("suppFcnTest UI") do
	srand(8675309); data = rand(200, 2)
	w = UIOracle(data, [0, 0], [1, 1], .1, .2)

	zstar, ustar = suppFcn([1, 1], w, :Max)
	@fact zstar --> roughly(1.8772156628616448, 1e-8)
	@fact ustar[1] --> roughly(0.9385081140320037, 1e-8)
	@fact ustar[2] --> roughly(0.9387075488296409, 1e-8)

	zstar, ustar = suppFcn([1, -1], w, :Max)
	@fact zstar --> roughly(0.8932039312511011, 1e-8)
	@fact ustar[1] --> roughly(0.9492096614168394, 1e-8)
	@fact ustar[2] --> roughly(0.056005730165738316, 1e-8)

	zstar, ustar = suppFcn([-1, 1], w, :Max)
	@fact zstar --> roughly(0.8767942540553981, 1e-8)
	@fact ustar[1] --> roughly(0.06360683171372103, 1e-8)
	@fact ustar[2] --> roughly(0.9404010857691192, 1e-8)

	#extra tests for Min and 0 entries
	zstar, ustar = suppFcn([-1, 1], w, :Min)
	@fact zstar --> roughly(-0.893203928668029, 1e-8)
	@fact ustar[1] --> roughly(0.9492096601291875, 1e-8)
	@fact ustar[2] --> roughly(0.056005731461158545, 1e-8)

	zstar, ustar = suppFcn([0, 1], w, :Min)
	@fact zstar --> roughly(0.001943343938739852, 1e-8)
	@fact ustar[1] --> roughly(0.5715291128670335, 1e-8)
	@fact ustar[2] --> roughly(0.001943343938739852, 1e-8)

	zstar, ustar = suppFcn([0, 1], w, :Max)
	@fact zstar --> roughly(0.9981094005216729, 1e-8)
	@fact ustar[1] --> roughly(0.5715291128670335, 1e-8)
	@fact ustar[2] --> roughly(0.9981094005216729, 1e-8)

end

facts("portTest UI") do
	srand(8675309)
	data = rand(500, 2)
	w = UIOracle(data, [0., 0.], [1., 1.], .1, .2)
	portTest(w, 0.07068906144892116, [0.5039487187951147, 0.4960512812048853])
end


text:
# Quantum trajectory simulation of spin model for two-level atoms as presented in
# Marco T. Manzoni, Darrick E. Chang and James S. Douglas,
# Nature Communications 8, 1743 (2017)

using MAT
using MatrixProductStates
BLAS.set_num_threads(2)

# spin model parameters
const na = 50 # number of atoms
const gam_1d = 1.0*ones(na) # coupling of probe to each atom
const gam_eg = 0.0 # eg spontaneous decay rate
const del_p = 0.0 # detuning of pump beam from eg transition del_p = ωp - ωeg
const k_wg = 0.2*pi # waveguide wavevector
const rj = collect(1.0:na) # atom positions
const k_in = 0.5*pi # pump beam wavevector
const f_amp = 0.0 # probe amplitude
const sigma_t = 6.0 # time width of gaussian input pulse
const t_peak = 10.0 # time peak of gaussian input pulse
f(t) = f_amp*(pi*sigma_t^2/4)^(-1/4)*exp(-2(t - t_peak)^2/sigma_t^2) # time envelope of pump beam

# simulation parameters
const setrand = 1 # random number generator seed
const dt = 0.01 # time step
const t_fin = 50.0 # final time
const d = 2 # dimension of atom Hilbert space
const d_max = 50 # maximum bond dimension
const measure_int = 10 # number of time steps per observable measurement
const path_data = "/home/jdouglas/data/" # save directory
const base_filename = string(path_data,"SM_N", na, "_D", d_max,
	"_Tf", t_fin, "_dt", dt, "_traj", setrand)

# local spin operators
const hgg = [1 0; 0 0]
const hee = [0 0; 0 1]
const hge = [0 1; 0 0]
const heg = hge'
const id = eye(2)

# types used in mps representation
const TN = Complex{Float64}
const TA = Array

function time_evolve()

    # construct left-canonical initial MPS ground state where the
    #A, dims = mpsgroundstate(TN, TA, na, d_max, d)
    A, dims = mpsproductstate(TN, TA, na, d_max, d, [0, 1])

    # create jump and measurement operators
    jumpleft = makempo(TN, TA, [jj->id jj->im*sqrt(gam_1d[jj]/2)*exp(im*k_wg*rj[jj])*hge;
                                jj->0 jj->id], na, d)
    jumpright = makempo(TN, TA, [jj->id jj->im*sqrt(gam_1d[jj]/2)*exp(-im*k_wg*rj[jj])*hge;
                                 jj->0 jj->id], na, d)
    jumpright[1][1, :, 2, :] = f(0.0)*id + im*sqrt(gam_1d[1]/2)*exp(-im*k_wg*rj[1])*hge
    ir_mpo = applyMPOtoMPO(jumpright, conj_mpo(jumpright))
    ir2_mpo = applyMPOtoMPO(applyMPOtoMPO(jumpright, ir_mpo), conj_mpo(jumpright))

    # step times and measurement times
    t = 0.0:dt:t_fin
    t_m = t[1:measure_int:end]
    tstep = length(t) - 1

    # initialize output measures
    p_coh = zeros(tstep)
    p_sum = zeros(tstep)
    p_acc = zeros(tstep, na + 3)
    e_pop = zeros(na, length(t_m))
    epop = zeros(na)
    I_r = zeros(TN, length(t_m))
    I2_r = zeros(TN, length(t_m))
    times = zeros(tstep, 4)
    t_r = Float64[]
    t_l = Float64[]
    t_eg = Float64[]
    pos_eg = Float64[]
    cont = zeros(Int64, 3)

    # initial measurements
    I_r[1] = scal_op_prod(A, ir_mpo, A)
    I2_r[1] = scal_op_prod(A, ir2_mpo, A)
    measure_excitations!((@view e_pop[:, 1]), A)

    # temporary arrays for time evolution
    envop = build_env(TN, TA, dims, dims, ones(na + 1)*4)
    envop_jump = build_env(TN, TA, dims, dims, ones(na + 1)*2)
    A_coh = similar(A)
    A_r = similar(A)
    A_l = similar(A)
    A_coh .= copy.(A)
    A_r .= copy.(A)
    A_l .= copy.(A)

    # random jump variables
    srand(setrand)
    r1 = rand(tstep)
    r2 = rand(tstep)

    # create linear evolution MPO at initial time
    expH = sm_hamiltonian(TN, TA, 1, dt, t[1])

    for i = 1:tstep
        time_in = time()

        # update time-evolution operator
        update_h!(expH, 1, dt, t[i])

        # coherent time evolution probability
        compress_var_apply_H!(A_coh, envop, A, A, expH, 1)
        p_coh[i] = normalize_lo!(A_coh)
        times[i,1] = time() - time_in

        # coherent evolution or jumps
        if r1[i] < p_coh[i]
            A .= copy.(A_coh)
        else
            # update waveguide jump operator
            jumpright[1][1, :, 2, :] = f(t[i])*id +
                    im*sqrt(gam_1d[1]/2)*exp(-im*k_wg*rj[1])*hge

            # jump probabilities for decay into the waveguide
            compress_var_apply_H!(A_r, envop_jump, A, A, jumpright, 1)
            p_r = normalize_lo!(A_r)
            compress_var_apply_H!(A_l, envop_jump, A, A, jumpleft, 1)
            p_l = normalize_lo!(A_l)

            # jump probabilities decay outside the waveguide
            measure_excitations!(epop, A)
            p_acc[i, 1] = 0.0
            for j = 1:na
                p_acc[i, j + 1] = p_acc[i, j] + dt*gam_eg*epop[j]
            end
            p_acc[i, na + 2] = p_acc[i, na + 1] + dt*p_r
            p_acc[i, na + 3] = p_acc[i, na + 2] + dt*p_l
            p_sum[i] = p_acc[i, end]
            println(p_sum[i] + p_coh[i])
            println(t[i])
            p_acc[i, :] = p_acc[i, :]/p_acc[i, end]

            # jumps according to above probabilities
            if r2[i] < p_acc[i, na + 1]
                for j = 1:na
                    if r2[i] < p_acc[i, j + 1]
                        apply_site_operator!(A[j], hge)
                        leftorthonormalizationQR!(A)
                        push!(t_eg, t[i])
                        push!(pos_eg, j)
                        cont[1] += 1
                        break
                    end
                end
            elseif r2[i] < p_acc[i, na + 2]
                A .= copy.(A_r)
                push!(t_r, t[i])
                cont[2] += 1
            else
                A .= copy.(A_l)
                push!(t_l, t[i])
                cont[3] += 1
            end
            times[i, 2] = time() - times[i, 1] - time_in

        end

        if rem(i, measure_int) == 0

            mind = div(i, measure_int) + 1
            # Atomic state populations and cavity occupation number
            measure_excitations!((@view e_pop[:, mind]), A)
            times[i, 3] = time() - times[i, 2]  - times[i, 1] - time_in

            # Output right field update
            jumpright[1][1, :, 2, :] = f(t[i + 1])*id +
                    im*sqrt(gam_1d[1]/2)*exp(-im*k_wg*rj[1])*hge
            ir_mpo[1] =  apply_site_MPOtoMPO(jumpright[1], conj_site_mpo(jumpright[1]))
            ir2_mpo[1] =  apply_site_MPOtoMPO(
                apply_site_MPOtoMPO(jumpright[1], ir_mpo[1]),
                conj_site_mpo(jumpright[1]))

            # right field output observables
            I_r[mind] = scal_op_prod(A, ir_mpo, A)
            I2_r[mind] = scal_op_prod(A, ir2_mpo, A)
            times[i,4] = time() - times[i, 3] - times[i, 2] - times[i, 1] - time_in

        end

        # write temporary file
        if rem(i,measure_int*10) == 0
            write_data_file(string(base_filename, "_temp.mat"),
                            t_m, e_pop, p_coh, p_acc, p_sum, I_r,
                            I2_r, cont, t_r, t_l, t_eg, A, times)
        end

    end

    # saving final data
    write_data_file(string(base_filename, ".mat"),
                           t_m, e_pop, p_coh, p_acc, p_sum, I_r,
                           I2_r, cont, t_r, t_l, t_eg, A, times)

end

# calculate local expectation values, assuming mps is already left normalized
function measure_excitations!(e_pop, A)

    # initialize right environment
    env = ones(1, 1)
    # sweep left taking expectation value
    for n = na:(-1):1
        Aenv = prod_LR(A[n], env)
        e_pop[n] = real(local_exp(Aenv, A[n], hee))
        env = update_renv(A[n], Aenv)
    end

end

# creates the spin model hamiltonian when con = 0 and delt = 1
# with con = 1 and delt = dt creates the linear time evolution operator 1 - im*dt*H
function sm_hamiltonian(::Type{TN}, ::Type{TA}, con, delt, time) where {TN, TA}

    H = Array{TA{TN, 4}, 1}(na + 1)

    drj = diff(rj)
    ph = exp.(im*k_wg*drj)
    cp = sqrt.(delt*gam_1d/2)

    H[1] = zeros(TN, 1, 2, 4, 2)
    H[1][1, :, 1, :] = id
    H[1][1, :, 2, :] = -cp[1]*ph[1]*heg
    H[1][1, :, 3, :] = -cp[1]*ph[1]*hge
    H[1][1, :, 4, :] = local_h(1, con, delt, time)

    H[na] = zeros(TN, 4, 2, 1, 2)
    H[na][1, :, 1, :] = id
    H[na][1, :, 1, :] = local_h(na, con, delt, time)
    H[na][2, :, 1, :] = cp[na]*hge
    H[na][3, :, 1, :] = cp[na]*heg
    H[na][4, :, 1, :] = id

    for jj = 2:(na - 1)
        H[jj] = zeros(TN, 4, 2, 4, 2)
        H[jj][1, :, 1, :] = id
        H[jj][1, :, 2, :] = -cp[jj]*ph[jj]*heg
        H[jj][1, :, 3, :] = -cp[jj]*ph[jj]*hge
        H[jj][1, :, 4, :] = local_h(jj, con, delt, time)
        H[jj][2, :, 4, :] = cp[jj]*hge
        H[jj][3, :, 4, :] = cp[jj]*heg
        H[jj][2, :, 2, :] = ph[jj]*id
        H[jj][3, :, 3, :] = ph[jj]*id
        H[jj][4, :, 4, :] = id
    end

    return H

end

# update of H created using vit_hamiltonian() for time dependent input pulse
function update_h!(H, con, delt, time)

    H[1][1, :, 4, :] = local_h(1, con, delt, time)

    H[na][1, :, 1, :] = local_h(na, con, delt, time)

    for jj = 2:(na - 1)
        H[jj][1, :, 4, :] = local_h(jj, con, delt, time)
    end

end

# local site Hamiltonian
function local_h(jj, con, delt, time)

    con*id/na + delt*((im*del_p - (gam_eg + gam_1d[jj])/2)*hee +
        im*sqrt(gam_1d[jj]/2)*f(time)*exp(im*k_in*rj[jj])*heg -
        (1/(2na))*abs2(f(time))*id)

end

function write_data_file(filename, t_m, e_pop, p_coh, p_acc,
                         p_sum, I_r, I2_r, cont, t_r, t_l, t_eg,
                         A, times)

    file = matopen(string(filename), "w")
    write(file, "TN",string(TN))
    write(file, "TA",string(TA))
    write(file, "na", na)
    write(file, "gam_1d", gam_1d)
    write(file, "gam_eg", gam_eg)
    write(file, "del_p", del_p)
    write(file, "k_wg", k_wg)
    write(file, "rj", rj)
    write(file, "k_in", k_in)
    write(file, "f_amp", f_amp)
    write(file, "t_peak", t_peak)
    write(file, "sigma_t",sigma_t)
    write(file, "f", f.(t_m))
    write(file, "dt", dt)
    write(file, "t_fin", t_fin)
    write(file, "d_max", d_max)
    write(file, "measure_int", measure_int)

    write(file, "t_m", collect(t_m))
    write(file, "e_pop", e_pop)
    write(file, "p_coh", p_coh)
    write(file, "p_acc", p_acc)
    write(file, "p_sum", p_sum)
    write(file, "I_r", I_r)
    write(file, "I2_r", I2_r)
    write(file, "cont", cont)
    write(file, "t_r", t_r)
    write(file, "t_l", t_l)
    write(file, "t_eg", t_eg)
    write(file, "A", collect.(A))
    write(file, "times", times)
    close(file)

end

time_evolve()


text:
function MainWindow()
  t = importOrders()
  win = @Window("test", 1920, 1080)
  (sells, buys) = groupby(t, :bid)
  listViewSells = DataFrameView(sells)
  listViewBuys = DataFrameView(buys)
  f = @Frame("orders")
  panes = @Paned(true)
  panes[1] = listViewSells
  panes[2] = listViewBuys
  push!(f, panes)
  juliaEval = @Window("Julia")
  evalDummy = @Frame()
  Gtk.GAccessor.attached_to(juliaEval, evalDummy)
  tabs = @Notebook()
  push!(tabs, f, "Orders")
  push!(tabs, evalDummy, "Julia")
  push!(win, tabs)
  showall(win)
end

max_stars_repo_name: barcharcraz/eve-julia-gui

max_issues_repo_name: barcharcraz/eve-julia-gui

max_forks_repo_name: barcharcraz/eve-julia-gui


text:
const RestrictionMatrix = BandedMatrix{<:Int, <:FillArrays.Ones}

const RestrictionTuple{B<:Basis} = Tuple{B, <:RestrictionMatrix}
const AdjointRestrictionTuple{B<:Basis} = Tuple{<:Adjoint{<:Any,<:RestrictionMatrix}, <:QuasiAdjoint{<:Any,B}}

const RestrictedBasis{B<:Basis} = Mul{<:Any,<:RestrictionTuple{B}}
const AdjointRestrictedBasis{B<:Basis} = Mul{<:Any,<:AdjointRestrictionTuple{B}}

const RestrictedQuasiArray{T,N,B<:Basis} = MulQuasiArray{T,N,<:RestrictionTuple{B}}
const AdjointRestrictedQuasiArray{T,N,B<:Basis} = MulQuasiArray{T,N,<:AdjointRestrictionTuple{B}}

const BasisOrRestricted{B<:Basis} = Union{B,RestrictedBasis{<:B},<:RestrictedQuasiArray{<:Any,<:Any,<:B}}
const AdjointBasisOrRestricted{B<:Basis} = Union{<:QuasiAdjoint{<:Any,B},AdjointRestrictedBasis{<:B},<:AdjointRestrictedQuasiArray{<:Any,<:Any,<:B}}


text:
function τ_lag(rule::Function, T::Array{Float64,1}, X::Array{Float64,1})::Float64

    # filter -
    idx = findfirst(rule, X)

    # return -
    return T[idx]
end

function τ_lag(rule::Function, T::Array{Float64,2}, X::Array{Float64,2})::Array{Float64,1}

    # initialize -
    time_lag_array = Array{Float64,1}()

    # size -
    (_, C) = size(X)
    for i ∈ 1:C

        # compute -
        τ = τ_lag(rule, T[:, i], X[:, i])

        # filter -
        push!(time_lag_array, τ)
    end

    # return -
    return time_lag_array
end

function max_FIIa(X::Array{Float64,2})::Array{Float64,1}

    # initialize -
    peak_array = Array{Float64,1}()

    # size -
    (R, C) = size(X)
    for i ∈ 1:C

        # get sim run -
        sim_col = X[:, i]

        # what is the max value?
        max_value = maximum(sim_col)

        # capture -
        push!(peak_array, max_value)
    end

    # return -
    return peak_array
end

function peak_time(T::Array{Float64,1}, X::Array{Float64,1})::Float64

    # find the index of the max value -
    idx = argmax(X)
    return T[idx]
end

function peak_time(T::Array{Float64,2}, X::Array{Float64,2})::Array{Float64,1}

    # initialize -
    peak_time_array = Array{Float64,1}()
    (R, C) = size(X)
    for i ∈ 1:C
        value = peak_time(T[:, i], X[:, i])
        push!(peak_time_array, value)
    end

    # return -
    return peak_time_array
end

function auc(T::Array{Float64,1}, X::Array{Float64,1})::Float64
    return integrate(T, X)
end

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
# tests for functors

using Base.Test
using NumericExtensions

## basic operators

@test evaluate(Negate(),   2) == -2
@test evaluate(Add(),      2, 3) == 5
@test evaluate(Subtract(), 2, 3) == -1
@test evaluate(Multiply(), 2, 3) == 6
@test evaluate(Divide(),   2, 3) == 2 / 3
@test evaluate(Modulo(),   5, 3) == 2
@test evaluate(Pow(),      2, 3) == 8

@test evaluate(Greater(),      2, 3) == false
@test evaluate(Greater(),      2, 2) == false
@test evaluate(GreaterEqual(), 2, 3) == false
@test evaluate(GreaterEqual(), 2, 2) == true

@test evaluate(Less(),      3, 2) == false
@test evaluate(Less(),      3, 2) == false
@test evaluate(LessEqual(), 3, 2) == false
@test evaluate(LessEqual(), 2, 2) == true

@test evaluate(Equal(), 3, 2)    == false
@test evaluate(Equal(), 2, 2)    == true
@test evaluate(NotEqual(), 3, 2) == true
@test evaluate(NotEqual(), 2, 2) == false

@test evaluate(Not(), true)  == false
@test evaluate(Not(), false) == true
@test evaluate(And(), true, false) == false
@test evaluate(Or(),  true, false) == true

@test evaluate(BitwiseNot(), 1) == ~1
@test evaluate(BitwiseAnd(), 3, 6) == 2
@test evaluate(BitwiseOr(),  3, 6) == 7
@test evaluate(BitwiseXor(), 3, 6) == 5

# elementary math

x = 0.8
y = 0.6
z = 1.2

@test evaluate(MinFun(), x, y) == min(x, y)
@test evaluate(MaxFun(), x, y) == max(x, y)

@test evaluate(AbsFun(), -0.5) == 0.5
@test evaluate(Abs2Fun(), x) == abs2(x)
@test evaluate(SqrFun(),  x) == sqr(x)
@test evaluate(SqrtFun(), x) == sqrt(x)
@test evaluate(CbrtFun(), x) == cbrt(x)
@test evaluate(RcpFun(),  x) == rcp(x)
@test evaluate(RsqrtFun(), x) == rsqrt(x)
@test evaluate(RcbrtFun(), x) == rcbrt(x)

@test evaluate(HypotFun(), x, y) == hypot(x, y)

@test evaluate(IsfiniteFun(), 1.0) == true
@test evaluate(IsfiniteFun(), Inf) == false
@test evaluate(IsfiniteFun(), NaN) == false

@test evaluate(IsinfFun(), 1.0) == false
@test evaluate(IsinfFun(), Inf) == true
@test evaluate(IsinfFun(), NaN) == false

@test evaluate(IsnanFun(), 1.0) == false
@test evaluate(IsnanFun(), Inf) == false
@test evaluate(IsnanFun(), NaN) == true

@test evaluate(FloorFun(), x) == floor(x)
@test evaluate(CeilFun(), x)  == ceil(x)
@test evaluate(RoundFun(), x) == round(x)
@test evaluate(TruncFun(), x) == trunc(x)

@test evaluate(IfloorFun(), x) == ifloor(x)
@test evaluate(IceilFun(), x)  == iceil(x)
@test evaluate(IroundFun(), x) == iround(x)
@test evaluate(ItruncFun(), x) == itrunc(x)

@test evaluate(ModFun(), x, y) == mod(x, y)
@test evaluate(RemFun(), x, y) == rem(x, y)
@test evaluate(FldFun(), x, y) == fld(x, y)
@test evaluate(DivFun(), 8, 6) == 1

@test evaluate(ExpFun(), x)   == exp(x)
@test evaluate(Exp2Fun(), x)  == exp2(x)
@test evaluate(Exp10Fun(), x) == exp10(x)
@test evaluate(Expm1Fun(), x) == expm1(x)

@test evaluate(LogFun(), x)   == log(x)
@test evaluate(Log2Fun(), x)  == log2(x)
@test evaluate(Log10Fun(), x) == log10(x)
@test evaluate(Log1pFun(), x) == log1p(x)

@test evaluate(SinFun(), x) == sin(x)
@test evaluate(CosFun(), x) == cos(x)
@test evaluate(TanFun(), x) == tan(x)
@test evaluate(AsinFun(), x) == asin(x)
@test evaluate(AcosFun(), x) == acos(x)
@test evaluate(AtanFun(), x) == atan(x)
@test evaluate(Atan2Fun(), x, y) == atan2(x, y)

# ternary functors

@test evaluate(FMA(), x, y, z) == (x + y * z)
@test evaluate(IfelseFun(), true, y, z) == y
@test evaluate(IfelseFun(), false, y, z) == z

# Test unary functors

for e in [
    (Negate, -),
    (AbsFun, abs), (Abs2Fun, abs2), (SqrtFun, sqrt), (CbrtFun, cbrt),
    (FloorFun, floor), (CeilFun, ceil), (RoundFun, round), (TruncFun, trunc),
    (ExpFun, exp), (Exp2Fun, exp2), (Exp10Fun, exp10), (Expm1Fun, expm1),
    (LogFun, log), (Log2Fun, log2), (Log10Fun, log10), (Log1pFun, log1p),
    (SinFun, sin), (CosFun, cos), (TanFun, tan), (AsinFun, asin), (AcosFun, acos), (AtanFun, atan),
    (SinhFun, sinh), (CoshFun, cosh), (TanhFun, tanh), (AsinhFun, asinh), (AcoshFun, acosh), (AtanhFun, atanh),
    (ErfFun, erf), (ErfcFun, erfc), (GammaFun, gamma), (LgammaFun, lgamma), (DigammaFun, digamma),
    (IsfiniteFun, isfinite), (IsinfFun, isinf), (IsnanFun, isnan)]

    T = e[1]
    f = e[2]

    @test evaluate(T(), 1.0) === f(1.0)
    @test evaluate(T(), 1) === f(1)
end

# Test binary functors

for e in [
    (Add, +),
    (Subtract, -),
    (Multiply, *),
    (Divide, /),
    (Pow, ^),
    (MaxFun, max),
    (MinFun, min),
    (HypotFun, hypot),
    (Atan2Fun, atan2),
    (Greater, >),
    (GreaterEqual, >=),
    (Less, <),
    (LessEqual, <=),
    (Equal, ==),
    (NotEqual, !=)]

    T = e[1]
    f = e[2]

    @test evaluate(T(), 2.0, 3.0) === f(2.0, 3.0)
    @test evaluate(T(), 2, 3) === f(2, 3)
    @test evaluate(T(), 2.0, 3) === f(2.0, 3)
    @test evaluate(T(), 2, 3.0) === f(2, 3.0)
end


text:
const LengthType=typeof(1.0u"m");

"""
    MonoLightField2D

`MonoLightField2D` is used to represent a 2-dimensional monochromatic light field.

# Examples

```jldoctest; setup = :(using OpticalPropagation)
julia> using Unitful

julia> MonoLightField2D([1 2; 3 4],wavelength=632.8u"nm",size=(1u"mm",1u"mm"))
MonoLightField2D:
    wavelength: 6.328e-7 m
    physical size: (0.001 m, 0.001 m)
    complex amplitude data:
2×2 Matrix{ComplexF64}:
 1.0+0.0im  2.0+0.0im
 3.0+0.0im  4.0+0.0im
```
"""
struct MonoLightField2D <: AbstractArray{ComplexF64,2}
    data::Array{ComplexF64,2}
    wavelength::LengthType
    size::Tuple{LengthType,LengthType}

    MonoLightField2D(data::AbstractMatrix; wavelength::Unitful.Length, size::Tuple{<: Unitful.Length,<: Unitful.Length}) =
        new(ComplexF64.(data), float(wavelength), float.(size))
end

function MonoLightField2D(x::MonoLightField2D; kwargs...)
    dist_data=get(kwargs,:data,x.data)
    wavelen=get(kwargs,:wavelength,x.wavelength)
    phy_size=get(kwargs,:size,x.size)
    MonoLightField2D(dist_data,wavelength=wavelen,size=phy_size)
end

function Base.show(io::IO, x::MonoLightField2D)
    print(io, "MonoLightField2D(")
    show(io, x.data)
    print(io, ", wavelength=")
    show(io, x.wavelength)
    print(io, ", size=")
    show(io, x.size)
    print(io, ")")
end

function Base.show(io::IO, mime::MIME"text/plain", x::MonoLightField2D)
    print(io,"MonoLightField2D:")
    print(io,"\n    wavelength: ")
    show(io, mime, x.wavelength)
    print(io, "\n    physical size: ")
    show(io, mime, x.size)
    print(io,"\n    complex amplitude data:\n")
    show(io, mime, x.data)
end

function +(a::MonoLightField2D,b::MonoLightField2D)
    if ! (a.wavelength≈b.wavelength)
        throw(DimensionMismatch("Wavelengths $(a.wavelength) and $(b.wavelength) are not equal."))
    elseif ! all(a.size .≈ b.size)
        throw(DimensionMismatch("Physical size $(a.size) and $(b.size) are not equal."))
    else
        MonoLightField2D(a.data .+ b.data,
            wavelength=a.wavelength,
            size=a.size)
    end
end

function -(a::MonoLightField2D,b::MonoLightField2D)
    if ! (a.wavelength≈b.wavelength)
        throw(DimensionMismatch("Wavelengths $(a.wavelength) and $(b.wavelength) are not equal."))
    elseif ! all(a.size .≈ b.size)
        throw(DimensionMismatch("Physical size $(a.size) and $(b.size) are not equal."))
    else
        MonoLightField2D(a.data .- b.data,
            wavelength=a.wavelength,
            size=a.size)
    end
end

*(r::Number, x::MonoLightField2D) =
    MonoLightField2D(r .* x.data,
        wavelength = x.wavelength,
        size = x.size)

*(x::MonoLightField2D, r::Number) = r*x

/(x::MonoLightField2D, r::Number) = (1/r)*x

\(r::Number, x::MonoLightField2D) = (1/r)*x

function Base.isapprox(a::MonoLightField2D,b::MonoLightField2D; kwargs...)
    isapprox(a.wavelength,b.wavelength,kwargs...) &&
      all(isapprox.(a.size,b.size,kwargs...)) &&
      isapprox(a.data,b.data,kwargs...)
end

Base.size(a::MonoLightField2D) = size(a.data)
Base.getindex(a::MonoLightField2D, inds::Vararg{Int,2}) = a.data[inds...]
Base.setindex!(a::MonoLightField2D, val, inds::Vararg{Int,2}) = a.data[inds...] = val

Base.BroadcastStyle(::Type{MonoLightField2D}) = Broadcast.ArrayStyle{MonoLightField2D}()
function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{MonoLightField2D}}, ::Type{ElType}) where ElType
    # Scan the inputs for the MonoLightField2D:
    A = find_monolightfield2d(bc)
    MonoLightField2D(A, data = similar(Array{ElType}, axes(bc)))
end

find_monolightfield2d(bc::Base.Broadcast.Broadcasted) = find_monolightfield2d(bc.args)
find_monolightfield2d(::Tuple{}) = nothing
find_monolightfield2d(args::Tuple) = find_monolightfield2d(find_monolightfield2d(args[1]), Base.tail(args))
find_monolightfield2d(x) = x
find_monolightfield2d(a::MonoLightField2D, rest) = a
find_monolightfield2d(::Any, rest) = find_monolightfield2d(rest)

@recipe function plot(light::MonoLightField2D; datafunction=abs, unit = u"mm")
    (m, n) = size(light.data)
    x = ((-1/2):(1/m):(1/2-1/m))*uconvert(Unitful.NoUnits, light.size[1]/unit)
    y = ((-1/2):(1/n):(1/2-1/n))*uconvert(Unitful.NoUnits, light.size[2]/unit)
    z = datafunction.(light.data)
    seriestype --> :heatmap
    aspect_ratio --> :equal
    xguide --> unit
    yguide --> unit
    size --> (m, n)
    x, y, z
end


text:
# An experimental d3 backend for compose.

export D3

# In d3 svg we use millimeters as the base measure.
type D3Measure <: NativeMeasure
    value::Float64
end

*(u::Float64, v::D3Measure) = D3Measure(u * v.value)
+(u::D3Measure, v::D3Measure) = D3Measure(u.value + v.value)
-(u::D3Measure, v::D3Measure) = D3Measure(u.value - v.value)
convert(::Type{Float64}, u::D3Measure) = u.value
convert(::Type{D3Measure}, u::Float64) = D3Measure(u)
function convert(::Type{SimpleMeasure{MillimeterUnit}}, u::D3Measure)
    SimpleMeasure{MillimeterUnit}(u.value)
end

function show(io::IO, u::D3Measure)
    @printf(io, "%smm", svg_fmt_float(u.value))
end

type D3 <: Backend
    # Image size in millimeters
    width::D3Measure
    height::D3Measure

    # Output stream
    out::IO

    # Should f be closed when the backend is finished?
    close_stream::Bool

    # Current level of indentation.
    indentation::Int

    # Code to be inserted after the end of a property group.
    hooks::Vector{String}

    # Keep track of which properties that are push are empty to we can avoid
    # printiing them.
    empty_properties::Vector{Bool}

    # Form count used to generate unique classes for DataForms.
    dataform_count::Int

    # A counter to assign unique ids to clip paths.
    clippath_count::Int

    # Store datasets for serialization.
    data::Dict{Uint64, (Any, Int)}

    # Emit the graphic on finish when writing to a buffer.
    emit_on_finish::Bool

    function D3(f::IO,
                width::MeasureOrNumber,
                height::MeasureOrNumber,
                emit_on_finish::Bool=true)
        img = new()
        img.width  = native_measure(width, img)
        img.height = native_measure(height, img)
        img.out = f
        img.close_stream = false
        img.indentation = 0
        img.hooks = Array(String, 0)
        img.empty_properties = Array(Bool, 0)
        img.dataform_count = 0
        img.clippath_count = 0
        img.data = Dict{Uint64, (Any, Int)}()
        img.emit_on_finish = emit_on_finish

        width_value = svg_fmt_float(width.value)
        height_value = svg_fmt_float(height.value)
        write(img.out,
              """
              function draw_with_data(data, parent_id) {
                var g = d3.select(parent_id)
                          .append("svg")
                            .attr("width", "$(width_value)mm")
                            .attr("height", "$(height_value)mm")
                            .attr("viewBox", "0 0 $(width_value) $(height_value)")
                            .attr("stroke-width", "0.5")
                            .attr("style", "stroke:black;fill:black");
                g.append("defs");
                var t = {"scale": 1.0};
              """)
        img
    end

    function D3(filename::String,
                width::MeasureOrNumber,
                height::MeasureOrNumber)
        f = open(filename, "w")
        img = D3(f, width, height)
        img.close_stream = true
        img
    end

    function D3(width::MeasureOrNumber, height::MeasureOrNumber,
                emit_on_finish::Bool=true)
        img = D3(IOString(), width, height, emit_on_finish=emit_on_finish)
        img.close_stream = false
        img
    end
end

json(c::ColorValue) = repr("#$(hex(c))")

function write_data(img::D3, d::AbstractArray)
    write(img.out, "  [")
    n = length(d)
    for (i, x) in enumerate(d)
        write(img.out, json(x))
        if i < n
            write(img.out, ",")
        end
    end
    write(img.out, "]")
end

function write_data(img::D3)
    write(img.out, "var data = [\n")
    datapairs = Array(Tuple, length(img.data))
    for (i, (_, (d, j))) in enumerate(img.data)
        datapairs[i] = (j, d)
    end
    sort!(datapairs, by=jd -> jd[1])

    for (i, (_, d)) in enumerate(datapairs)
        write_data(img, d)
        if i < length(img.data)
            write(img.out, ",\n")
        end
    end
    write(img.out, "];\n\n")
end

function finish(img::D3)
    write(img.out, "}\n\n")
    write_data(img)
    write(img.out,
        """
        var draw = function(parent_id) {
            draw_with_data(data, parent_id);
        };
        """)
    if img.close_stream
        close(img.out)
    end

    if img.emit_on_finish && typeof(img.out) == IOString
        seek(img.out, 0)
        emit(Emitable("application/javascript", readall(img.out)))
        close(img.out)
    end
end

# Generate a unique class for a data form
function next_dataform_class(backend::D3)
    class = @sprintf("form%d", backend.dataform_count)
    backend.dataform_count += 1
    class
end

# Index of the the given array in the serialized data.
function data_idx(backend::D3, d::AbstractArray)
    id = object_id(d)
    if !haskey(backend.data, d)
        backend.data[id] = (d, length(backend.data))
    end
    backend.data[id][2]
end

# Conversion to D3 units (i.e. millimeters)
function native_measure(u::SimpleMeasure{PixelUnit},
                        img::D3)
    native_measure(convert(SimpleMeasure{MillimeterUnit}, u), img)
end

function native_measure(u::SimpleMeasure{MillimeterUnit},
                        img::D3)
    D3Measure(u.value)
end

function root_box(img::D3)
    NativeBoundingBox(
        D3Measure(0.0),
        D3Measure(0.0),
        img.width,
        img.height)
end

native_zero(::D3) = D3Measure(0.0)

function indent(img::D3)
    for i in 1:img.indentation
        write(img.out, "  ")
    end
end

# Nop catchall
function draw(img::D3, form::FormPrimitive)
end

function draw(img::D3, form::Polygon)
    n = length(form.points)
    if n <= 1; return; end

    indent(img)
    write(img.out, "g.append(\"svg:path\")\n")
    indent(img)
    write(img.out, "   .attr(\"d\", \"")
    print_svg_path(img.out, form.points)
    write(img.out, " z\");\n")
end

function draw(img::D3, form::Lines)
    n = length(form.points)
    if n <= 1; return; end

    indent(img)
    write(img.out, "g.append(\"svg:path\")\n")
    indent(img)
    write(img.out, "   .attr(\"d\", \"")
    print_svg_path(img.out, form.points)
    write(img.out, "\");\n")
end

# Converting text with tango markup to d3 code. Messy business.

function pango_tag_to_tspan(out, tag::String, indent::String)
    print(out, ".append(\"tspan\")\n")
    if tag == "sup"
        print(out, indent, "  ", ".attr(\"baseline-shift\", \"super\")\n")
    elseif tag == "sub"
        print(out, indent, "  ", ".attr(\"baseline-shift\", \"sub\")\n")
    end
end

function pango_to_d3(out, text::String, indentnum=0)
    print(out, Base.repeat(" ", indentnum))
    print(out, ".call(function(text) {\n")

    indentnum += 2; indent = Base.repeat(" ", indentnum)
    print(out, indent, "text");
    indentnum += 4; indent = Base.repeat(" ", indentnum)

    opentagpat = r"<\s*([^>\s]*)\s*>"
    openmat = match(opentagpat, text)
    if !is(openmat, nothing)

        # find end tag
        closemat = match(Regex(string(".*(</\\s*", openmat.captures[1], "\\s*>)")), text)
        if closemat === nothing
            error("No matching tag for <$(openmat.captures[1])>")
        end

        tag = openmat.captures[1]
        i0 = openmat.offset
        i1 = i0 + length(openmat.match)
        j0 = closemat.offsets[1]
        j1 = j0 + length(closemat.captures[1])

        if i0 > 1
            print(out, ".append(\"tspan\").text(\"$(escape_string(text[1:i0-1]))\")\n")
            print(out, indent)
        end

        pango_tag_to_tspan(out, tag, indent)
        pango_to_d3(out, text[i1:j0-1], indentnum)

        if j1 < length(text)
            print(out, indent, ".append(\"tspan\").text(\"$(escape_string(text[j1:end]))\")\n")
        end
        print(out, indent, ";\n");
    else
        print(out, ".text(\"$(escape_string(text))\");\n")
    end

    indentnum -= 6; indent = Base.repeat(" ", indentnum)
    print(out, indent, "})\n");
end

function draw(img::D3, form::Text)
    indent(img)
    write(img.out, "g.append(\"svg:text\")\n")
    indent(img)
    @printf(img.out, "   .attr(\"x\", %s)\n", svg_fmt_float(form.pos.x.value))
    indent(img)
    @printf(img.out, "   .attr(\"y\", %s)\n", svg_fmt_float(form.pos.y.value))

    if is(form.halign, hcenter)
        indent(img)
        print(img.out, "   .attr(\"text-anchor\", \"middle\")\n")
    elseif is(form.halign, hright)
        indent(img)
        print(img.out, "   .attr(\"text-anchor\", \"end\")\n")
    end

    if is(form.valign, vcenter)
        indent(img)
        print(img.out, "   .style(\"dominant-baseline\", \"central\")\n")
    elseif is(form.valign, vtop)
        indent(img)
        print(img.out, "   .style(\"dominant-baseline\", \"text-before-edge\")\n")
    end

    if !isidentity(form.t)
        indent(img)
        @printf(img.out, "   .attr(\"transform\", \"rotate(%s, %s, %s)\")\n",
                svg_fmt_float(radians2degrees(atan2(form.t.M[2,1], form.t.M[1,1]))),
                svg_fmt_float(form.pos.x.value),
                svg_fmt_float(form.pos.y.value))
    end

    indent(img)
    pango_to_d3(img.out, form.value)
    print(img.out, ";\n")
    #@printf(img.out, "   .text(\"%s\");\n",
            #escape_string(pango_to_svg(form.value)))
end

# Nop catchall
function push_property(img::D3, property::Property)
    p = property
    hasproperties = false;
    clip = nothing
    while !is(p, empty_property)
        if !is(typeof(p.primitive), D3Hook)
            hasproperties = true
        end
        if typeof(p.primitive) === Clip
                clip = p.primitive
        end
        p = p.next
    end

    # Nothing to do
    if property === empty_property
        push!(img.empty_properties, true)
        return
    end
    push!(img.empty_properties, false)

    indent(img)
    write(img.out, "(function (g) {\n")
    img.indentation += 1
    indent(img)

    if !is(clip, nothing)
        clipid = @sprintf("clippath%d", img.clippath_count)
        img.clippath_count += 1
        write(img.out,
            """
                d3.select("defs")
                  .append("svg:clipPath")
                    .attr("id", parent_id + "_$(clipid)")
                    .append("svg:path")
                      .attr("d", " """)
        print_svg_path(img.out, clip.points)
        write(img.out, " z\");")
    end

    if hasproperties
        write(img.out, "g")
    end
    p = property
    hook = ""
    while !is(p, empty_property)
        if typeof(p.primitive) === D3Hook
            hook = string(hook, "\n", p.primitive.code)
        else
            apply_property(img, p.primitive)
        end
        if !is(p.next , empty_property)
            write(img.out, "\n ")
            indent(img)
        end
        p = p.next
    end
    if hasproperties
        write(img.out, ";\n");
    end
    push!(img.hooks, hook)
end

function pop_property(img::D3)
    if pop!(img.empty_properties)
        return
    end

    hook = pop!(img.hooks)
    write(img.out, hook)
    img.indentation -= 1;
    indent(img)
    write(img.out, "}(g.append(\"g\")));\n")
end

function apply_property(img::D3, p::PropertyPrimitive)

end

function apply_property(img::D3, p::Stroke)
    @printf(img.out, ".attr(\"stroke\", \"%s\")", svg_fmt_color(p.value))
end

function apply_property(img::D3, p::Fill)
    @printf(img.out, ".attr(\"fill\", \"%s\")", svg_fmt_color(p.value))
end

function apply_property(img::D3, p::LineWidth)
    @printf(img.out, ".attr(\"stroke-width\", \"%s\")",
            svg_fmt_float(p.value.value))
end

function apply_property(img::D3, p::Font)
    @printf(img.out, ".attr(\"font-family\", \"%s\")",
            escape_string(p.family))
end

function apply_property(img::D3, p::FontSize)
    @printf(img.out, ".attr(\"font-size\", \"%s\")",
            svg_fmt_float(p.value.value))
end

function apply_property(img::D3, p::Clip)
    clipid = @sprintf("clippath%d", img.clippath_count - 1)
    @printf(img.out, ".attr(\"clip-path\", \"url(#\" + parent_id + \"_%s)\")", clipid)
end

function apply_property(img::D3, p::SVGID)
    @printf(img.out, ".attr(\"id\", \"%s\"\)", escape_string(p.value))
end

function apply_property(img::D3, p::SVGClass)
    @printf(img.out, ".attr(\"class\", \"%s\"\)", escape_string(p.value))
end

function apply_property(img::D3, p::SVGAttribute)
    @printf(img.out, ".attr(\"%s\", \"%s\")", p.attribute, p.value)
end

function apply_property(img::D3, p::D3Embed)
    print(img.out, p.code)
end


text:
# this file defines iterators used for looping over a grid
struct UpdateFlags
    nodes::Bool
    coords::Bool
    celldofs::Bool
end

UpdateFlags(; nodes::Bool=true, coords::Bool=true, celldofs::Bool=true) =
    UpdateFlags(nodes, coords, celldofs)

"""
    CellIterator(grid::Grid)
    CellIterator(grid::DofHandler)

Return a `CellIterator` to conveniently loop over all the cells in a grid.

# Examples
```julia
for cell in CellIterator(grid)
    coords = getcoordinates(cell) # get the coordinates
    dofs = celldofs(cell)         # get the dofs for this cell
    reinit!(cv, cell)             # reinit! the FE-base with a CellIterator
end
```
"""
struct CellIterator{dim,N,T,M}
    flags::UpdateFlags
    grid::Grid{dim,N,T,M}
    current_cellid::ScalarWrapper{Int}
    nodes::Vector{Int}
    coords::Vector{Vec{dim,T}}
    dh::DofHandler{dim,N,T,M}
    celldofs::Vector{Int}

    function CellIterator{dim,N,T,M}(dh::DofHandler{dim,N,T,M}, flags::UpdateFlags) where {dim,N,T,M}
        cell = ScalarWrapper(0)
        nodes = zeros(Int, N)
        coords = zeros(Vec{dim,T}, N)
        n = ndofs_per_cell(dh)
        celldofs = zeros(Int, n)
        return new{dim,N,T,M}(flags, dh.grid, cell, nodes, coords, dh, celldofs)
    end

    function CellIterator{dim,N,T,M}(grid::Grid{dim,N,T,M}, flags::UpdateFlags) where {dim,N,T,M}
        cell = ScalarWrapper(0)
        nodes = zeros(Int, N)
        coords = zeros(Vec{dim,T}, N)
        return new{dim,N,T,M}(flags, grid, cell, nodes, coords)
    end
end

CellIterator(grid::Grid{dim,N,T,M},     flags::UpdateFlags=UpdateFlags()) where {dim,N,T,M} =
    CellIterator{dim,N,T,M}(grid, flags)
CellIterator(dh::DofHandler{dim,N,T,M}, flags::UpdateFlags=UpdateFlags()) where {dim,N,T,M} =
    CellIterator{dim,N,T,M}(dh, flags)

# iterator interface
function Base.iterate(ci::CellIterator, state = 1)
    if state > getncells(ci.grid)
        return nothing
    else
        return (reinit!(ci, state), state+1)
    end
end
Base.length(ci::CellIterator)  = getncells(ci.grid)

Base.IteratorSize(::Type{T})   where {T<:CellIterator} = Base.HasLength() # this is default in Base
Base.IteratorEltype(::Type{T}) where {T<:CellIterator} = Base.HasEltype() # this is default in Base
Base.eltype(::Type{T})         where {T<:CellIterator} = T

# utility
@inline getnodes(ci::CellIterator) = ci.nodes
@inline getcoordinates(ci::CellIterator) = ci.coords
@inline nfaces(ci::CellIterator) = nfaces(eltype(ci.grid.cells))
@inline onboundary(ci::CellIterator, face::Int) = ci.grid.boundary_matrix[face, ci.current_cellid[]]
@inline cellid(ci::CellIterator) = ci.current_cellid[]
@inline celldofs!(v::Vector, ci::CellIterator) = celldofs!(v, ci.dh, ci.current_cellid[])
@inline celldofs(ci::CellIterator) = ci.celldofs

function reinit!(ci::CellIterator{dim,N}, i::Int) where {dim,N}
    nodeids = ci.grid.cells[i].nodes
    ci.current_cellid[] = i
    @inbounds for j in 1:N
        nodeid = nodeids[j]
        ci.flags.nodes  && (ci.nodes[j] = nodeid)
        ci.flags.coords && (ci.coords[j] = ci.grid.nodes[nodeid].x)
    end
    if isdefined(ci, :dh) && ci.flags.celldofs # update celldofs
        celldofs!(ci.celldofs, ci)
    end
    return ci
end

@inline reinit!(cv::CellValues{dim,T}, ci::CellIterator{dim,N,T}) where {dim,N,T} = reinit!(cv, ci.coords)
@inline reinit!(fv::FaceValues{dim,T}, ci::CellIterator{dim,N,T}, face::Int) where {dim,N,T} = reinit!(fv, ci.coords, face)


text:
############################################################################################
function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::U,
    x::AbstractArray
) where {T<:AbstractFloat,F<:FlattenTypes,U<:UnflattenTypes}
## Obtain flatten/unflatten constructor for each element
    x_vecs_and_backs = map(x) do xᵢ
        construct_flatten(T, flattentype, unflattentype, xᵢ)
    end
    _flatten, _unflatten = first.(x_vecs_and_backs), last.(x_vecs_and_backs)
## Run flatten once to obtain dimensions of individual elements
    x_vecs = map(iter -> _flatten[iter](x[iter]), eachindex(x))
    lngth = map(length, x_vecs)
    sz = cumsum(lngth)
## Return flatten/unflatten for AbstractArray
    function construct_flatten_AbstractArr(x::AbstractArray)
        #!NOTE: reduce(vcat) seems to be faster for Array than tuple, thus speedincrease
        reduce(vcat,
            map(iter -> _flatten[iter](x[iter]), eachindex(x))
            )
    end
    function construct_unflatten_AbstractArr(v::Union{R,AbstractVector{R}}) where {R<:Real}
        return [
            _unflatten[n](view(v, (sz[n] - lngth[n] + 1):sz[n])) for n in eachindex(x)
        ]
    end
    return construct_flatten_AbstractArr, construct_unflatten_AbstractArr
end

############################################################################################
function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::U,
    constraint::AbstractArray,
    x::AbstractArray
) where {T<:AbstractFloat,F<:FlattenTypes,U<:UnflattenTypes}
## Obtain flatten/unflatten constructor for each element
    x_vecs_and_backs = map(x, constraint) do xᵢ, constraintᵢ
        construct_flatten(T, flattentype, unflattentype, constraintᵢ, xᵢ)
    end
    _flatten, _unflatten = first.(x_vecs_and_backs), last.(x_vecs_and_backs)
## Run flatten once to obtain dimensions of individual elements
    x_vecs = map(iter -> _flatten[iter](x[iter]), eachindex(x))
    lngth = map(length, x_vecs)
    sz = cumsum(lngth)
## Return flatten/unflatten for AbstractArray
    function construct_flatten_AbstractArr(x::AbstractArray)
        #!NOTE: reduce(vcat) seems to be faster for Array than tuple, thus speedincrease
        reduce(vcat,
            map(iter -> _flatten[iter](x[iter]), eachindex(x))
            )
    end
    function construct_unflatten_AbstractArr(v::Union{R,AbstractVector{R}}) where {R<:Real}
        return [
            _unflatten[n](view(v, (sz[n] - lngth[n] + 1):sz[n])) for n in eachindex(x)
        ]
    end
    return construct_flatten_AbstractArr, construct_unflatten_AbstractArr
end

############################################################################################
#!NOTE: If we map over collections, need to extend this to params.jl functions:
# _check, construct_transform, constrain, unconstrain, log_abs_det_jac

############################################################################################
function _check(
    _rng::Random.AbstractRNG,
    constraint::AbstractArray,
    x::AbstractArray
)
    return all(
            map(constraint, x) do cᵢ, xᵢ
                _check(_rng, cᵢ, xᵢ)
            end
    )
end

############################################################################################
function construct_transform(
    constraint::AbstractArray,
    x::AbstractArray
)
## Obtain transform/inversetransform constructor for each element
    x_transforms = map(x, constraint) do xᵢ, constraintᵢ
        construct_transform(constraintᵢ, xᵢ)
    end
    _transform, _inversetransform = first.(x_transforms), last.(x_transforms)
## Return flatten/unflatten for AbstractArray
    return _transform, _inversetransform
end

############################################################################################
function constrain(
    b⁻¹ᵥ::AbstractArray,
    θᵤ::AbstractArray
)
    return map(constrain, b⁻¹ᵥ, θᵤ)
end

############################################################################################
function unconstrain(
    bᵥ::AbstractArray,
    θ::AbstractArray
)
    return map(unconstrain, bᵥ, θ)
end

############################################################################################
function log_abs_det_jac(
    bᵥ::AbstractArray,
    θ::AbstractArray
)
    return sum(map(log_abs_det_jac, bᵥ, θ))
end

############################################################################################
#Export
export
    construct_flatten,
    construct_transform,
    constrain,
    unconstrain,
    log_abs_det_jac


text:
function rrule(::typeof(eigen), X::AbstractMatrix{<:Real};k=1)
    F = eigen(X)
    function eigen_pullback(Ȳ::Composite{<:Eigen})
        ∂X = @thunk(eigen_rev(F, Ȳ.values, Ȳ.vectors,k))
        return (NO_FIELDS, ∂X)
    end
    return F, eigen_pullback
end

function rrule(::typeof(getproperty), F::T, x::Symbol) where T <: Eigen
    function getproperty_eigen_pullback(Ȳ)
        C = Composite{T}
        ∂F = if x === :values
            C(values=Ȳ,)
        elseif x === :vectors
            C(vectors=Ȳ,)
        end
        return NO_FIELDS, ∂F, DoesNotExist()
    end
    return getproperty(F, x), getproperty_eigen_pullback
end

function eigen_rev(ΛV::Eigen,Λ̄,V̄,k)

    Λ = ΛV.values
    V = ΛV.vectors
    A = V*diagm(Λ)*V'

    Ā = zeros(size(A))
    tempĀ = zeros(size(A))
    # eigen(A).values are in descending order, this current implementation
    # assumes that the input matrix positive definite
    for j = length(Λ):-1:1
        tempĀ = (I-V[:,j]*V[:,j]') ./ norm(A*V[:,j])
        for i = 1:k-1
            tempĀ += A^i * (I-V[:,j]*V[:,j]') ./ (norm(A*V[:,j])^(i+1))
        end
        tempĀ *= V̄[:,j]*V[:,j]'
        # x = A ./ norm(A*V[:,j])
        # p = fill(I, (k - 1))
        # evalpoly(x, p)
        # v = V[:,j]
        # w = (v̄ .-v .* dot(v, v̄)) ./ norm(A*v)
        # tempĀ = (tempĀ * w) * v'
        Ā += tempĀ
        A = A - A*V[:,j]*V[:,j]'
    end
    return Ā
end

# using BLAS

function frule(
    (_, ΔA),
    ::typeof(eigen!),
    A::LinearAlgebra.RealHermSymComplexHerm{<:BLAS.BlasReal,<:StridedMatrix};
    sortby::Union{Function,Nothing}=nothing,
)
    F = eigen!(A; sortby=sortby)
    ΔA isa AbstractZero && return F, ΔA
    λ, U = F.values, F.vectors
    tmp = U' * ΔA
    ∂K = mul!(ΔA.data, tmp, U)
    ∂Kdiag = @view ∂K[diagind(∂K)]
    ∂λ = real.(∂Kdiag)
    ∂K ./= λ' .- λ
    fill!(∂Kdiag, 0)
    ∂U = mul!(tmp, U, ∂K)
    _eigen_norm_phase_fwd!(∂U, A, U)
    ∂F = Composite{typeof(F)}(values = ∂λ, vectors = ∂U)
    return F, ∂F
end

function rrule(
    ::typeof(eigen),
    A::LinearAlgebra.RealHermSymComplexHerm;
    sortby::Union{Function,Nothing}=nothing,
)
    F = eigen(A; sortby=sortby)
    function eigen_pullback(ΔF::Composite{<:Eigen})
        λ, U = F.values, F.vectors
        Δλ, ΔU = ΔF.values, ΔF.vectors
        ∂A = eigen_rev(A, λ, U, Δλ, ΔU)
        return NO_FIELDS, ∂A
    end
    eigen_pullback(ΔF::AbstractZero) = (NO_FIELDS, ΔF)
    return F, eigen_pullback
end

function eigen_rev(A::LinearAlgebra.RealHermSymComplexHerm, λ, U, ∂λ, ∂U)
    if ∂U isa AbstractZero
        ∂λ isa AbstractZero && return (NO_FIELDS, ∂λ + ∂U)
        ∂K = Diagonal(∂λ)
        ∂A = U * ∂K * U'
    else
        ∂U = copyto!(similar(∂U), ∂U)
        _eigen_norm_phase_rev!(∂U, A, U)
        ∂K = U' * ∂U
        ∂K ./= λ' .- λ
        ∂K[diagind(∂K)] = ∂λ
        ∂A = mul!(∂K, U * ∂K, U')
    end
    return _hermitrize!(∂A, A)
end

_eigen_norm_phase_fwd!(∂V, ::LinearAlgebra.RealHermSym, V) = ∂V
function _eigen_norm_phase_fwd!(∂V, A::Hermitian, V)
    k = A.uplo === 'U' ? size(A, 1) : 1
    @inbounds for i in axes(V, 2)
        v = @view V[:, i]
        vₖ, ∂vₖ = real(v[k]), ∂V[k, i]
        ∂v .-= v .* (imag(∂vₖ) / ifelse(iszero(vₖ), one(vₖ), vₖ))
    end
    return ∂V
end

_eigen_norm_phase_rev!(∂V, ::LinearAlgebra.RealHermSym, V) = ∂V
function _eigen_norm_phase_rev!(∂V, A::Hermitian, V)
    k = A.uplo === 'U' ? size(A, 1) : 1
    @inbounds for i in axes(V, 2)
        v, ∂v = @views V[:, i], ∂V[:, i]
        vₖ = real(v[k])
        ∂c = dot(v, ∂v)
        ∂v[k] -= im * (imag(∂c) / ifelse(iszero(vₖ), one(vₖ), vₖ))
    end
    return ∂V
end

#####
##### `eigvals!`/`eigvals`
#####

function frule(
    (_, ΔA),
    ::typeof(eigvals!),
    A::LinearAlgebra.RealHermSymComplexHerm{<:BLAS.BlasReal,<:StridedMatrix};
    sortby::Union{Function,Nothing}=nothing,
)
    ΔA isa AbstractZero && return eigvals!(A; sortby=sortby), ΔA
    F = eigen!(A; sortby=sortby)
    λ, U = F.values, F.vectors
    tmp = ΔA * U
    # diag(U' * tmp) without computing matrix product
    ∂λ = similar(λ)
    @inbounds for i in eachindex(λ)
        ∂λ[i] = @views real(dot(U[:, i], tmp[:, i]))
    end
    return λ, ∂λ
end

function rrule(
    ::typeof(eigvals),
    A::LinearAlgebra.RealHermSymComplexHerm;
    sortby::Union{Function,Nothing}=nothing,
)
    F = eigen(A; sortby=sortby)
    λ = F.values
    function eigvals_pullback(Δλ)
        U = F.vectors
        ∂A = _hermitrize!(U * Diagonal(Δλ) * U', A)
        return NO_FIELDS, ∂A
    end
    eigvals_pullback(Δλ::AbstractZero) = (NO_FIELDS, Δλ)
    return λ, eigvals_pullback
end

# in-place hermitrize matrix, optionally wrapping like A
function _hermitrize!(A)
    A .= (A .+ A') ./ 2
    return A
end
function _hermitrize!(∂A, A)
    _hermitrize!(∂A)
    return _symhermtype(A)(∂A, Symbol(A.uplo))
end

Zygote.@adjoint enumerate(xs) = enumerate(xs), diys -> (map(last, diys),)
_ndims(::Base.HasShape{d}) where {d} = d
_ndims(x) = Base.IteratorSize(x) isa Base.HasShape ? _ndims(Base.IteratorSize(x)) : 1
Zygote.@adjoint function Iterators.product(xs...)
                    d = 1
                    Iterators.product(xs...), dy -> ntuple(length(xs)) do n
                        nd = _ndims(xs[n])
                        dims = ntuple(i -> i<d ? i : i+nd, ndims(dy)-nd)
                        d += nd
                        func = sum(y->y[n], dy; dims=dims)
                        ax = axes(xs[n])
                        reshape(func, ax)
                    end
                end

function sum2(op,arr)
    return sum(op,arr)
end

function sum2adj( Δ, op, arr )
    n = length(arr)
    g = x->Δ*Zygote.gradient(op,x)[1]
    return ( nothing, map(g,arr))
end

Zygote.@adjoint function sum2(op,arr)
    return sum2(op,arr),Δ->sum2adj(Δ,op,arr)
end

ChainRulesCore.refresh_rules()
Zygote.refresh()


text:
struct RectilinearGrid{TX, TY, TZ, R, A1, A3} <: AbstractGrid{TX, TY, TZ}
    # corrdinates at cell centers, unit: meter
    xC::R
    yC::R
    zC::A1
    # corrdinates at cell faces, unit: meter
    xF::R
    yF::R
    zF::A1
    # grid spacing, unit: meter
    Δx::Float64
    Δy::Float64
    dzC::A1
    dzF::A1
    # number of grid points
    Nx::Int
    Ny::Int
    Nz::Int
    # number of halo points
    Hx::Int
    Hy::Int
    Hz::Int
    # landmask to indicate where is the land
    landmask::A3
end

"""
    RectilinearGrid(;size, x, y, z,
                     topology = (Periodic, Periodic, Bounded),
                     landmask = nothing,
                     halo = (2, 2, 2))
Creats a `RectilinearGrid` struct with `size = (Nx, Ny, Nz)` grid points.
    `x` and `y` directions must be regular spaced,
    `z` direction can be vertically stretched or regular spaced.

Keyword Arguments (Required)
============================
- `size` : A tuple prescribing the number of grid points.
                `size` is a 3-tuple no matter for 3D, 2D, or 1D model.
- `x` and `y` : A 2-tuple that specify the start and end points of the domain.
- `z` : is either a (1) 1D array that specifies the locations of cell faces in z direction,
                or (2) 2-tuples that specify the start and end points of the domain.
                Vertical indexing starts from surface and use negative numbers for depth.

Keyword Arguments (Optional)
============================
- `topology` : A 3-tuple specifying the topology of the domain.
                The topology can be either Periodic or Bounded in each direction.
- `landmask` : a 3-dimentional array to indicate where the land is.
- `halo` : A tuple of integers that specifies the size of the halo region of cells
                surrounding the physical interior for each direction.
                `halo` is a 3-tuple no matter for 3D, 2D, or 1D model.
                At least 2 halo points are needed for DST3FL advection scheme.
"""
function RectilinearGrid(;size, x, y, z,
                         topology = (Periodic, Periodic, Bounded),
                         landmask = nothing,
                         halo = (2, 2, 2))
    Nx, Ny, Nz = size
    Hx, Hy, Hz = halo
    TX, TY, TZ = validate_topology(topology)

    @assert isa(x, Tuple{<:Number, <:Number})
    @assert isa(y, Tuple{<:Number, <:Number})
    x₁, x₂ = x
    y₁, y₂ = y
    @assert x₁ < x₂
    @assert y₁ < y₂
    Δx = (x₂ - x₁) / Nx
    Δy = (y₂ - y₁) / Ny

    xF = range(-Hx * Δx, (Nx + Hx - 1) * Δx, length = Nx + 2 * Hx)
    yF = range(-Hy * Δy, (Ny + Hy - 1) * Δy, length = Ny + 2 * Hy)

    xC = range((0.5 - Hx) * Δx, (Nx + Hx - 0.5) * Δx, length = Nx + 2 * Hx)
    yC = range((0.5 - Hy) * Δy, (Ny + Hy - 0.5) * Δy, length = Ny + 2 * Hy)

    if isa(z, Tuple{<:Number, <:Number})
        z₁, z₂ = z
        z = collect(range(z₁, z₂, length = Nz+1))
    elseif isa(z, AbstractVector)
        z₁, z₂ = z[1], z[end]
    else
        throw(ArgumentError("z must be a 2-tuple or 1D array"))
    end

    @assert z₁ > z₂
    @assert Base.length(z) == Nz + 1

    zF = zeros(Nz+2Hz)
    zC = zeros(Nz+2Hz)
    dzF = zeros(Nz+2Hz)
    dzC = zeros(Nz+2Hz)

    zF[1+Hz:Nz+Hz+1] .= Float64.(z)
    zC[1+Hz:Nz+Hz] .= (zF[1+Hz:Nz+Hz] .+ zF[2+Hz:Nz+Hz+1]) ./ 2
    dzF[1+Hz:Nz+Hz] .= zF[1+Hz:Nz+Hz] .- zF[2+Hz:Nz+Hz+1]
    dzC[1+Hz:Nz+Hz-1] .= zC[1+Hz:Nz+Hz-1] .- zC[2+Hz:Nz+Hz]

    ##### fill halos
    @views @. dzF[1:Hz] = dzF[Hz+1]
    @views @. dzC[1:Hz] = dzF[Hz+1]
    @views @. dzF[Nz+Hz+1:Nz+2*Hz] = dzF[Nz+Hz]
    @views @. dzC[Nz+Hz:Nz+2*Hz] = dzF[Nz+Hz]

    @views @. zF[1:Hz] = zF[1+Hz:2*Hz] + dzF[1]*Hz
    @views @. zF[Nz+Hz+2:Nz+2*Hz] = zF[Nz+2:Nz+Hz] - dzF[end]*Hz
    for i in 1:Hz
        zC[Nz+Hz+i] = zC[Nz+Hz+i-1] - dzF[end]
    end
    for i in Hz:-1:1
        zC[i] = zC[i+1] + dzF[1]
    end

    landmask = landmask_validation(landmask, Nx, Ny, Nz, Hx, Hy, Hz, TX, TY)

    return RectilinearGrid{TX, TY, TZ, typeof(xF), typeof(zF), typeof(landmask)}(
        xC, yC, zC, xF, yF, zF, Δx, Δy, dzC, dzF, Nx, Ny, Nz, Hx, Hy, Hz, landmask)
end

function show(io::IO, g::RectilinearGrid{TX, TY, TZ}) where {TX, TY, TZ}
    xL, xR = g.xF[g.Hx+1], g.xF[g.Hx+1+g.Nx]
    yL, yR = g.yF[g.Hy+1], g.yF[g.Hy+1+g.Ny]
    zL, zR = g.zF[g.Hz+1], g.zF[g.Hz+1+g.Nz]
    dzF_min = minimum(g.dzF)
    dzF_max = maximum(g.dzF)
    print(io, "RegularRectilinearGrid{$TX, $TY, $TZ}\n",
              "domain: x ∈ [$xL, $xR], y ∈ [$yL, $yR], z ∈ [$zL, $zR]\n",
              "topology (Tx, Ty, Tz):     ", (TX, TY, TZ), '\n',
              "resolution (Nx, Ny, Nz):   ", (g.Nx, g.Ny, g.Nz), '\n',
              "halo size (Hx, Hy, Hz):    ", (g.Hx, g.Hy, g.Hz), '\n',
              "grid spacing (Δx, Δy, Δz): ", g.Δx, ", ", g.Δy, ", [min=", dzF_min, ", max=", dzF_max,"])")
end

function short_show(grid::RectilinearGrid{TX, TY, TZ}) where {TX, TY, TZ}
    return "RegularRectilinearGrid{$TX, $TY, $TZ}(Nx=$(grid.Nx), Ny=$(grid.Ny), Nz=$(grid.Nz))"
end

function validate_topology(topology)
    for T in topology
        if !isa(T(), AbstractTopology)
            e = "$T is not a valid topology!" *
                "Valid topologies are: Periodic and Bounded."
            throw(ArgumentError(e))
        end
    end

        return topology
end


text:
type SolverParams
    numbeats::Int64
    numbeatstotal::Int64
    numsteps::Int64
    nstart::Int64
    JL::Int64
    acols::Vector{Int8}
    qcols::Vector{Int8}
    acolspre::Vector{Int8}
    qcolspre::Vector{Int8}
    acolscor::Vector{Int8}
    qcolscor::Vector{Int8}
    colsint::Vector{Int8}
    CFL::Float64
    h::Float64
    tshift::Float64
    rho::Float64
    mu::Float64
    nu::Float64
    diffusioncoeff::Float64
    abserror::Float64
    maxiter::Int16
    epsJ::Float64
    epsN::Float64
    maxval::Float64
    totaliter::Int64
    th::Float64

    function SolverParams()
        this = new()
        this.CFL = 1.0;
        this.JL = 5;
        this.acols = [1:this.JL;];
        this.qcols = [this.JL+1:2*this.JL;];
        this.acolspre = [2:this.JL-1;];
        this.qcolspre = [this.JL+2:2*this.JL-1;];
        this.acolscor = [1:this.JL-2;];
        this.qcolscor = [this.JL-1:2*this.JL-4;];
        this.colsint = [2:this.JL-1;];
        this.nstart = 0;
        this.numbeats = 0;
        this.numbeatstotal = 0;
        this.tshift = 0;
        this.rho = 1060;
        this.mu = 0.004;
        this.nu = 0.5;
        this.diffusioncoeff = 22.0;
        this.abserror = 1e-10;
        this.maxiter = 200;
        this.epsJ = 1e-11;
        this.epsN = 1e-8;
        this.maxval = 1e7;
        this.totaliter = 0;
        this.th = 0.8;
        return this
    end
end