text:
# IJulia does not use IPython magic syntax, because there are better
# alternatives in Julia, whereas "magics" are inherently limited to
# running in Jupyter cells only.   The alternatives in Julia, which
# can work in any Julia code regardless of whether it is in a notebook
# cell, are:
#     * To interact with the running kernel, we can use IJulia.foo(...),
#       since the running IJulia kernel is loaded into Main.
#     * For implementing special syntax, Julia has @macros and
#       stringmacros"....".
# However, to ease the transition for users coming from other Jupyter kernels,
# we implement some magics that tell the user how to do the corresponding
# thing in Julia.   These are executed whenever a cell beginning with %
# is encountered.

# regex to find the presence of magics in a code-cell string
const magics_regex = r"^\s*(%%?[A-Za-z_][A-Za-z_0-9]*\b)\s*(.*)\s*$"m

# given the cell contents "code" (which start with %foo) output help
# that indicates the Julia equivalent, if possible.
function magics_help(code::AbstractString)
    for m in eachmatch(magics_regex, code)
        display(get(magic_help, m[1], generic_magic_help)(m[1], m[2]))
    end
end

using Markdown

const magic_help_string = """
    Julia does not use the IPython `%magic` syntax.   To interact
    with the IJulia kernel, use `IJulia.somefunction(...)`, for
    example.  Julia macros, string macros, and functions can be used to
    accomplish most of the other functionalities of IPython magics."""
generic_magic_help(magic::AbstractString, args::AbstractString) =
    Base.Markdown.parse("""
        Unrecognized magic `$magic`.

        $magic_help_string""")

lsmagic_help(magic::AbstractString, args::AbstractString) =
    Base.Markdown.parse("""
    $magic_help_string

    The Julia analogues of many IPython magics are printed if
    you try to input the IPython magic in a code cell.""")

alias_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The equivalent of `%alias` magic in Julia is simply to write
    a function that [runs a command object](http://docs.julialang.org/en/latest/manual/running-external-programs/).

    For example, the equivalent of `%alias bracket echo "Input in brackets: <%l>"`
    in IPython would be the Julia function

        bracket(l) = run(`echo "Input in brackets: <$l>"`)

    which you can then run with e.g. `bracket("hello world")`."""

function cd_magic_help(magic::AbstractString, args::AbstractString)
    if magic == "%cd" && !occursin(r"\s*-", args)
        return md"""The equivalent of `%cd 'dir'` in IPython is `cd("dir")` in Julia."""
    else
        return md"""
    The equivalent of `%cd 'dir'` in IPython is `cd("dir")` in Julia.

    Julia does not keep a history of visited directories, so it
    has no built-in equivalent of IPython's `%cd -`, `%dhist`,
    `%pushd`, and `%popd`.

    However, if you are interested, it would be easy to write a
    Julia function `mycd(x)` that worked like IPython's `%cd`
    (accepting a string, `-`, `-n`, and so on, and keeping a stack
    of directory history).  This might make a nice little Julia package
    if you are sufficiently motivated."""
    end
end

debug_magic_help(magic::AbstractString, args::AbstractString) =
md"""The Julia interactive debugger is provided by the [Gallium](https://github.com/Keno/Gallium.jl) package."""

edit_magic_help(magic::AbstractString, args::AbstractString) = md"""
    An analogue of IPython's `%edit` magic is provided by
    the `edit` function and the `@edit` macro in Julia, as
    described [in the Julia manual](http://docs.julialang.org/en/latest/stdlib/base/#Base.edit).

    For example, instead of `%edit -n line 'filename'` in IPython, you would do
    `edit("filename", line)` followed by `include("filename")` if you
    want to execute the file when you are done editing.  The analogue
    of `%edit` (which creates a temporary file) is `n=tempname(); edit(n)`
    followed by `include(n)`.

    If you have a function `foo()`, then `%edit foo` in IPython opens
    the file where `foo` was defined.  The analogue of this in Julia is
    to do `@edit foo()`   Note that you need to supply sample arguments
    to `foo` (which is not actually called), so that Julia knows *which*
    method of `foo` you want to edit.

    If `s` is a string variable, IPython's `%edit s` loads its contents
    into an editor.  In Julia, you would have to explicitly write
    to a temporary file with `n=tempname(); write(n, s); edit(n)` followed
    by `include(n)` to execute it."""

env_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The equivalent of `%env` in IPython, which lists all environment
    variables, is `ENV` in Julia, and the equivalent of
    `%env var=val` is `ENV["var"]=val`.

    (The `ENV` built-in global variable in Julia provides a dictionary-like
    interface to the system environment variables.)
"""

gui_magic_help(magic::AbstractString, args::AbstractString) = md"""
    There is an analogue of IPython's `%gui` in the
    [PyCall package](https://github.com/stevengj/PyCall.jl)
    for calling Python from Julia.   If you have done `using PyCall`
    to load PyCall, then the analogue of `%gui wx` is `pygui_start(:wx)`,
    and the analogue of `%gui` is `pygui_stop()`.  See the PyCall
    documentation for more information.

    Other Julia packages for GUI toolkits, e.g. the
    [Tk package](https://github.com/JuliaGraphics/Tk.jl) or
    the [Gtk package](https://github.com/JuliaGraphics/Gtk.jl), also
    provide their own event-loop integration."""

history_magic_help(magic::AbstractString, args::AbstractString) = md"""
    An analogue of the `%history` or `%hist` magic of IPython, which provides
    access to the input history, is given by IJulia.history(). It is
    based on the global variable `In` in IJulia. `In` is a dictionary
    mapping cell numbers to the inputs. However, IJulia does not currently
    keep any other history, e.g. it discards input cells that you overwrite.
"""

load_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of `%load filename` in IPython is `IJulia.load("filename")`
    in IJulia (to load code from `filename` into the current frontend)

    The analogue of `%load url` is `IJulia.load(download("url"))`.
"""

paste_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of `%paste` in IPython is `IJulia.load_string(clipboard())`
    in IJulia (to load code from the clipboard into the current frontend).
"""

matplotlib_magic_help(magic::AbstractString, args::AbstractString) =
Base.Markdown.parse("""
    The analogue of IPython's `$magic` in Julia is to use
    the [PyPlot package](https://github.com/stevengj/PyPlot.jl),
    which gives a Julia interface to Matplotlib including inline
    plots in IJulia notebooks.   (The equivalent of `numpy` is already
    loaded by default in Julia.)

    Given PyPlot, the analogue of `$magic inline` is `using PyPlot`,
    since PyPlot defaults to inline plots in IJulia.

    To enable separate GUI windows in PyPlot, analogous to `$magic`,
    do `using PyPlot; pygui(true)`.   To specify a particular gui
    backend, analogous to `$magic gui`, you can either do
    `using PyPlot; pygui(:gui); using PyPlot; pygui(true)` (where
    `gui` is `wx`, `qt`, `tk`, or `gtk`), or you can do
    `ENV["MPLBACKEND"]=backend; using PyPlot; pygui(true)` (where
    `backend` is the name of a Matplotlib backend, like `tkagg`).

    For more options, see the PyPlot documentation.""")

pdef_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%pdef somefunction` in Julia is
    `methods(somefunction)`, which prints out all of the possible
    calling signatures of `somefunction`.  (Note that, unlike Python,
    Julia is based on multiple dispatch, so a given function will
    often be callable with many different argument signatures.)"""

pdoc_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%pdoc object` is `?object` in IJulia."""

pfile_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%pfile somefunction` (also `%psource`) is
    roughly `@less somefunction(somearguments...)` in Julia.   The
    reason that you need to supply sample arguments (the function is
    not actually evaluated) is because Julia functions can
    have multiple definitions for different argument types.
    `@less somefunction(somearguments...)` will print (or run
    through the `less` pager) the `somefunction` definition
    corresponding to the arguments you supply."""

file_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%file somefunction` is
    roughly `methods(somefunction)` in Julia.   This lists
    all the methods of `sumfunction` along with the locations where
    they are defined."""

precision_magic_help(magic::AbstractString, args::AbstractString) = md"""
    There currently is no way to globally set the output precision
    in Julia, analogous to `%precision` in IPython, as discussed
    in [Julia issue #6493](https://github.com/JuliaLang/julia/issues/6493)."""

prun_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%prun statement` in Julia is
    `@profile statement`, which runs the
    [Julia profiler](http://docs.julialang.org/en/latest/manual/profile/).
    The analogue of `%%prun ...code...` is
    ```
    @profile begin
        ...code...
    end
    ```
    Note, however, that you should put all performance-critical
    code into a function, avoiding global variables, before
    doing performance measurements in Julia; see the
    [performance tips in the Julia manual](http://docs.julialang.org/en/latest/manual/performance-tips/).

    See also the
    [ProfileView package](https://github.com/timholy/ProfileView.jl) for
    richer graphical display of profiling output."""

psearch_magic_help(magic::AbstractString, args::AbstractString) = md"""
    A rough analogue of IPython's `%psearch PATTERN` in Julia might be
    `filter(s -> occursin(r"PATTERN", string(s)), names(Base))`, which
    searches all the symbols defined in the `Base` module for a given
    regular-expression pattern `PATTERN`."""

pwd_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%pwd` is `pwd()` in Julia."""

qtconsole_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%qtconsole` is `IJulia.qtconsole()` in Julia."""

recall_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%recall n` in IJulia is
    `IJulia.load_string(In[n])`, and the analogue of
    `%recall n-m` is `IJulia.load_string(join([get(In,i,"") for i in n:m],"\n"))`."""

run_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%run file` is `include("file")` in Julia.

    The analogue of `%run -t file`, to print timing information, is
    `@time include("file")` in Julia.  The analogue of `%run -t -N n file`
    is `@time for i in 1:n; include("file"); end` in Julia.

    For running under a debugger, see the
    [Gallium Julia debugger](https://github.com/Keno/Gallium.jl).  To
    run other IJulia notebooks, see the
    [NBInclude package](https://github.com/stevengj/NBInclude.jl)."""

save_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%save filename n1-n2 n3-n4` is
    ```
    open("filename","w") do io
        for i in [n1:n2; n3:n4]
            println(get(In,i,""))
        end
    end
    ```
    in IJulia."""

sc_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%sc shell command` is `; shell command` in IJulia, or
        read(`shell command`, String)
    to capture the output as a string in Julia."""

set_env_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%set_env var val` is `ENV["var"]=val` in Julia."""

sx_magic_help(magic::AbstractString, args::AbstractString) =
Base.Markdown.parse("""
    The analogue of IPython's `$magic shell command` is
        split(read(`shell command`, String),'\n')
    in Julia.""")

time_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%time statement` (also `%timeit`) in Julia is
    `@time statement`.  The analogue of `%%time ...code...` is
    ```
    @time begin
        ...code...
    end
    ```
    Note, however, that you should put all performance-critical
    code into a function, avoiding global variables, before
    doing performance measurements in Julia; see the
    [performance tips in the Julia manual](http://docs.julialang.org/en/latest/manual/performance-tips/).

    The `@time` macro prints the timing results, and returns the
    value of evaluating the expression.  To instead return the time
    (in seconds), use `@elapsed statement`.

    For more extensive benchmarking tools, including the ability
    to collect statistics from multiple runs, see the
    [BenchmarkTools package](https://github.com/JuliaCI/BenchmarkTools.jl)."""

who_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%who` and `%whos` is `whos()` in Julia.

    You can also use `whos(r"PATTERN")` to find variables matching
    a given regular expression `PATTERN`."""

html_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%%html` is
    ```
    HTML(""\"
    ...html text...
    ""\")
    ```
    in Julia."""

javascript_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%%javascript ...code...` in Julia can be
    constructed by first evaluating
    ```
    macro javascript_str(s) display("text/javascript", s); end
    ```
    to define the `javascript"...."` [string macro](http://docs.julialang.org/en/latest/manual/strings/#non-standard-string-literals)
    in Julia.  Subsequently, you can simply do:
    ```
    javascript""\"
    ...code...
    ""\"
    ```
    to execute the script in an IJulia notebook."""

latex_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%%latex` is
    ```
    display("text/latex", ""\"
    ...latex text...
    ""\")
    ```
    in Julia.  Note, however, that `$` and `\` in the LaTeX text
    needs to be escaped as `\$` and `\\` so that they aren't
    interpreted by Julia.  See, however, the
    [LaTeXStrings package](https://github.com/stevengj/LaTeXStrings.jl)
    for easier input of LaTeX text as
    L""\"
    ...latex text...
    ""\"
    without requiring any extra backslashes.
"""

function pipe_magic_help(magic::AbstractString, args::AbstractString)
    cmd = magic[3:end] # magic is "%%cmd"
    if cmd == "script"
        arglist = split(args)
        cmd = isempty(arglist) ? "someprogram" : arglist[end]
        magic = "%%script $cmd"
    end
    Base.Markdown.parse("""
    The analogue of IPython's `$magic ...code...` in Julia can be
    constructed by first evaluating
    ```
    macro $(cmd)_str(s) open(`$cmd`,"w",stdout) do io; print(io, s); end; end
    ```
    to define the `$cmd"...."` [string macro](http://docs.julialang.org/en/latest/manual/strings/#non-standard-string-literals)
    in Julia.  Subsequently, you can simply do:
    ```
    $cmd""\"
    ...code...
    ""\"
    ```
    to evaluate the code in `$cmd` (outputting to `stdout`).""")
end

svg_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%%svg` is
    ```
    display("image/svg+xml", ""\"
    ...svg text...
    ""\")
    ```
    in Julia. To be even nicer, you can define
    ```
    macro svg_str(s) display("image/svg+xml", s); end
    ```
    to define the `svg"...."` [string macro](http://docs.julialang.org/en/latest/manual/strings/#non-standard-string-literals)
    in Julia.  Subsequently, you can simply do:
    ```
    svg""\"
    ...svg text...
    ""\"
    ```
    to display the SVG image.  Using a custom string macro like this
    has the advantage that you don't need to escape `$` and `\` if they
    appear in the SVG code."""

writefile_magic_help(magic::AbstractString, args::AbstractString) = md"""
    The analogue of IPython's `%%writefile filename` is
    `write("filename", In[IJulia.n])`.

    (`IJulia.n` is the index of the current code cell.  Of
    course, you can also use `In[N]` for some other `N` to output
    the contents of a different input cell.)"""

# map from magic to helpfunction(magic, magicargument)
const magic_help = Dict{String, Function}(
    "%alias" => alias_magic_help,
    "%cd" => cd_magic_help,
    "%dhist" => cd_magic_help,
    "%dirs" => cd_magic_help,
    "%popd" => cd_magic_help,
    "%pushd" => cd_magic_help,
    "%debug" => debug_magic_help,
    "%pdb" => debug_magic_help,
    "%edit" => edit_magic_help,
    "%env" => env_magic_help,
    "%gui" => gui_magic_help,
    "%hist" => history_magic_help,
    "%history" => history_magic_help,
    "%load" => load_magic_help,
    "%loadpy" => load_magic_help,
    "%paste" => paste_magic_help,
    "%lsmagic" => lsmagic_help,
    "%matplotlib" => matplotlib_magic_help,
    "%pylab" => matplotlib_magic_help,
    "%pdef" => pdef_magic_help,
    "%pdoc" => pdoc_magic_help,
    "%pinfo" => pdoc_magic_help,
    "%pinfo2" => pdoc_magic_help,
    "%pfile" => pfile_magic_help,
    "%file" => file_magic_help,
    "%psource" => pfile_magic_help,
    "%precision" => precision_magic_help,
    "%prun" => prun_magic_help,
    "%%prun" => prun_magic_help,
    "%psearch" => psearch_magic_help,
    "%pwd" => pwd_magic_help,
    "%qtconsole" => qtconsole_magic_help,
    "%recall" => recall_magic_help,
    "%run" => run_magic_help,
    "%save" => save_magic_help,
    "%sc" => sc_magic_help,
    "%set_env" => set_env_magic_help,
    "%sx" => sx_magic_help,
    "%system" => sx_magic_help,
    "%time" => time_magic_help,
    "%%time" => time_magic_help,
    "%timeit" => time_magic_help,
    "%%timeit" => time_magic_help,
    "%who" => who_magic_help,
    "%who_ls" => who_magic_help,
    "%whos" => who_magic_help,
    "%%html" => html_magic_help,
    "%%javascript" => javascript_magic_help,
    "%%latex" => latex_magic_help,
    "%%bash" => pipe_magic_help,
    "%%perl" => pipe_magic_help,
    "%%python" => pipe_magic_help,
    "%%python2" => pipe_magic_help,
    "%%python3" => pipe_magic_help,
    "%%ruby" => pipe_magic_help,
    "%%script" => pipe_magic_help,
    "%%sh" => pipe_magic_help,
    "%%svg" => svg_magic_help,
    "%%writefile" => writefile_magic_help,
)


text:
# script for Fig.12, Fig.13, Fig.14, Fig.15, Fig.16

using MultiscaleGraphSignalTransforms, JLD, Plots, LightGraphs, Distances
gr(dpi = 200)

## Build weighted toronto street network graph
G = loadgraph(joinpath(@__DIR__, "../datasets", "new_toronto_graph.lgz")); N = nv(G)
X = load(joinpath(@__DIR__, "../datasets", "new_toronto.jld"), "xy")
dist_X = pairwise(Euclidean(), X; dims = 1)
W = 1.0 .* adjacency_matrix(G)
Weight = zeros(N, N); Weight[W .> 0] = 1 ./ dist_X[W .> 0]; Weight = W .* Weight
L = Matrix(Diagonal(sum(Weight; dims = 1)[:]) - Weight)
𝛌, 𝚽 = eigen(L);
sgn = (maximum(𝚽, dims = 1)[:] .> -minimum(𝚽, dims = 1)[:]) .* 2 .- 1; 𝚽 = 𝚽 .* sgn';
Q = incidence_matrix(G; oriented = true)
edge_weight = 1 ./ sqrt.(sum((Q' * X).^2, dims = 2)[:])

## Build Dual Graph by DAG metric
distDAG = eigDAG_Distance(𝚽, Q, N; edge_weight = edge_weight) #52.375477 seconds
Gstar_Sig = dualgraph(distDAG)
G_Sig = GraphSig(W, xy = X); G_Sig = Adj2InvEuc(G_Sig)
GP_dual = partition_tree_fiedler(Gstar_Sig; swapRegion = false)
GP_primal = pairclustering(𝚽, GP_dual)

@time PC_NGWP = pc_ngwp(𝚽, GP_dual, GP_primal) #119.590168 seconds (12.14 M allocations: 158.035 GiB, 13.89% gc time)
@time VM_NGWP = vm_ngwp(𝚽, GP_dual) #1315.821724 seconds (3.05 M allocations: 495.010 GiB, 7.04% gc time)

#################### Fig. 12(a) a smooth spatial distribution of the street intersections graph signal
f = zeros(N); for i in 1:N; f[i] = length(findall(dist_X[:,i] .< 1/minimum(edge_weight))); end #fneighbor
G_Sig.f = reshape(f, (N, 1))
gplot(W, X; width=1); signal_plt = scatter_gplot!(X; marker = f, plotOrder = :s2l, ms = 3)
# savefig(signal_plt, joinpath(@__DIR__, "../paperfigs/Toronto_fdensity.png"))

#################### Fig. 12(b) spatial distribution signal relative l2 approximation error by various methods
DVEC = getall_expansioncoeffs(G_Sig, GP_dual, VM_NGWP, PC_NGWP, 𝚽)
approx_error_plot(DVEC);
plt = plot!(xguidefontsize = 14, yguidefontsize = 14, legendfontsize = 10)
# savefig(plt, joinpath(@__DIR__, "../paperfigs/Toronto_fdensity_DAG_approx.png"))

#################### Fig. 13 fdensity 16 most important VM-NGWP vectors (ignore the DC vector)
dmatrix_VM = ngwp_analysis(G_Sig, VM_NGWP)
dvec_vm_ngwp, BS_vm_ngwp = ngwp_bestbasis(dmatrix_VM, GP_dual)
important_idx = sortperm(dvec_vm_ngwp[:].^2; rev = true)
for i in 2:17
    dr, dc = BS_vm_ngwp.levlist[important_idx[i]]
    w = VM_NGWP[dr, dc, :]
    println("(j, k, l) = ", NGWP_jkl(GP_dual, dr, dc))
    gplot(W, X; width=1)
    scatter_gplot!(X; marker = w, plotOrder = :s2l, ms = 3)
    plt = plot!(cbar = false, clims = (-0.075,0.075))
    # savefig(plt, joinpath(@__DIR__,
    #     "../paperfigs/Toronto_fdensity_DAG_VM_NGW_important_basis_vector$(lpad(i,2,"0")).png"))
end

#################### Fig. 14(a) pedestrian volume graph signal
fp = load(joinpath(@__DIR__, "../datasets", "new_toronto.jld"), "fp")
G_Sig.f = reshape(fp, (N, 1))
gplot(W, X; width=1); signal_plt = scatter_gplot!(X; marker = fp, plotOrder = :s2l, ms = 3)
# savefig(signal_plt, joinpath(@__DIR__, "../paperfigs/Toronto_fp.png"))

#################### Fig. 14(b) pedestrian signal relative l2 approximation error by various methods
DVEC = getall_expansioncoeffs(G_Sig, GP_dual, VM_NGWP, PC_NGWP, 𝚽)
approx_error_plot(DVEC)
plt = plot!(xguidefontsize = 14, yguidefontsize = 14, legendfontsize = 10)
# savefig(plt, joinpath(@__DIR__, "../paperfigs/Toronto_fp_DAG_approx.png"))

#################### Fig. 15 pedestrian signal 16 most important VM-NGWP vectors
dmatrix_VM = ngwp_analysis(G_Sig, VM_NGWP)
dvec_vm_ngwp, BS_vm_ngwp = ngwp_bestbasis(dmatrix_VM, GP_dual)
important_idx = sortperm(dvec_vm_ngwp[:].^2; rev = true)
for i in 1:16
    dr, dc = BS_vm_ngwp.levlist[important_idx[i]]
    w = VM_NGWP[dr, dc, :]
    println("(j, k, l) = ", NGWP_jkl(GP_dual, dr, dc))
    gplot(W, X; width=1)
    scatter_gplot!(X; marker = w, plotOrder = :s2l, ms = 3)
    plt = plot!(cbar = false, clims = (-0.075,0.075))
    # savefig(plt, joinpath(@__DIR__,
    #     "../paperfigs/Toronto_fp_DAG_VM_NGW_important_basis_vector$(lpad(i,2,"0")).png"))
end

#################### Fig. 16 pedestrian signal 16 most important PC-NGWP vectors
dmatrix_PC = ngwp_analysis(G_Sig, PC_NGWP)
dvec_pc_ngwp, BS_pc_ngwp = ngwp_bestbasis(dmatrix_PC, GP_dual)
important_idx = sortperm(dvec_pc_ngwp[:].^2; rev = true)
for i in 1:16
    dr, dc = BS_pc_ngwp.levlist[important_idx[i]]
    w = PC_NGWP[dr, dc, :]
    println("(j, k, l) = ", NGWP_jkl(GP_dual, dr, dc))
    sgn = (maximum(w) > -minimum(w)) * 2 - 1
    gplot(W, X; width=1)
    scatter_gplot!(X; marker = sgn .* w, plotOrder = :s2l, ms = 3)
    plt = plot!(cbar = false, clims = (-0.075,0.075))
    # savefig(plt, joinpath(@__DIR__,
    #     "../paperfigs/Toronto_fp_DAG_PC_NGW_important_basis_vector$(lpad(i,2,"0")).png"))
end

ext: jl


text:
import Statistics
import DataFrames

# TODO: Fitness algorithm
function get_fitness(M::Array{Float64, 2}; iterations::Int, product_codes::Array{String, 1}, country_codes::Array{String, 1})

	# Check inputs:
	@assert size(product_codes, 1) == size(M, 2) "`product_codes` needs to be same length as width of `M`."
	@assert size(country_codes, 1) == size(M, 1) "`country_codes` needs to be same length as the height of `M`."

	# Define initial conditions:
	F_current = ones(size(M, 1), 1)
	Q_current = ones(size(M, 2), 1)

	# Run algorithm
	for i in 1:iterations
		F = M * Q_current
		Q = 1 ./ (M' * (1 ./ F_current))
		F_current = F / Statistics.mean(F)
		Q_current = Q / Statistics.mean(Q)
	end

	# Return results:
	complexity_tbl = DataFrames.DataFrame(product_code = product_codes, complexity = [Q_current...])
	fitness_tbl = DataFrames.DataFrame(country_code = country_codes, fitness = [F_current...])

	return Dict("Q" => complexity_tbl, "F" => fitness_tbl)
end

# TODO: HH algorithm

max_stars_repo_name: sbpost/compost

max_issues_repo_name: sbpost/compost

max_forks_repo_name: sbpost/compost


text:
@testset "Function tests" begin

    @testset "Logs" begin

        var_states = Dict{Int,Vector{Int}}()
        var_changes = Dict{Int,Vector{Tuple{Symbol,Int,Int,Int}}}()
        activity = Dict{Int,Float64}()
        children = CS.TreeLogNode{Int}[]
        l1 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            0,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            children,
        )
        l2 = CS.TreeLogNode(
            0,
            :Closed,
            true,
            0,
            0,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            children,
        )
        @test !CS.same_logs(l1, l2)

        # different children order
        children = CS.TreeLogNode{Int}[]
        tln1 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            1,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            CS.TreeLogNode{Int}[],
        )
        tln2 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            2,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            CS.TreeLogNode{Int}[],
        )
        push!(children, tln1)
        push!(children, tln2)
        l1 = CS.TreeLogNode(
            0,
            :Closed,
            false,
            0,
            0,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            children[[1, 2]],
        )
        l2 = CS.TreeLogNode(
            0,
            :Closed,
            false,
            0,
            0,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            children[[2, 1]],
        )
        @test !CS.same_logs(l1, l2)

        # missing one child
        children = CS.TreeLogNode{Float64}[]
        tln1 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0.0,
            1,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            CS.TreeLogNode{Float64}[],
        )
        tln2 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0.0,
            2,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            CS.TreeLogNode{Float64}[],
        )
        push!(children, tln1)
        push!(children, tln2)
        l1 = CS.TreeLogNode(
            0,
            :Closed,
            true,
            0.0,
            0,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            children[[1, 2]],
        )
        l2 = CS.TreeLogNode(
            0,
            :Closed,
            true,
            0.0,
            0,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            children[[1]],
        )
        @test !CS.same_logs(l1, l2)

        # 2 layer children different struture
        children1 = CS.TreeLogNode{Int}[]
        tln13 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            3,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            CS.TreeLogNode{Int}[],
        )
        tln11 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            1,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            [tln13],
        )
        tln12 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            2,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            CS.TreeLogNode{Int}[],
        )
        push!(children1, tln11)
        push!(children1, tln12)
        children2 = CS.TreeLogNode{Int}[]
        tln23 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            3,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            CS.TreeLogNode{Int}[],
        )
        tln21 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            1,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            CS.TreeLogNode{Int}[],
        )
        tln22 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            2,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            [tln23],
        )
        push!(children2, tln21)
        push!(children2, tln22)
        l1 = CS.TreeLogNode(
            0,
            :Closed,
            true,
            0,
            0,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            children1,
        )
        l2 = CS.TreeLogNode(
            0,
            :Closed,
            true,
            0,
            0,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            children2,
        )
        @test !CS.same_logs(l1, l2)

        # 2 layer children different position
        children1 = CS.TreeLogNode{Int}[]
        tln13 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            3,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            CS.TreeLogNode{Int}[],
        )
        tln14 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            4,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            CS.TreeLogNode{Int}[],
        )
        tln11 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            1,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            [tln13, tln14],
        )
        tln12 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            2,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            CS.TreeLogNode{Int}[],
        )
        push!(children1, tln11)
        push!(children1, tln12)
        children2 = CS.TreeLogNode{Int}[]
        tln23 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            3,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            CS.TreeLogNode{Int}[],
        )
        tln24 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            4,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            CS.TreeLogNode{Int}[],
        )
        tln21 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            1,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            [tln24, tln23],
        )
        tln22 = CS.TreeLogNode(
            0,
            :Open,
            true,
            0,
            2,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            CS.TreeLogNode{Int}[],
        )
        push!(children2, tln21)
        push!(children2, tln22)
        l1 = CS.TreeLogNode(
            0,
            :Closed,
            true,
            0,
            0,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            children1,
        )
        l2 = CS.TreeLogNode(
            0,
            :Closed,
            true,
            0,
            0,
            0,
            0,
            0,
            var_states,
            var_changes,
            activity,
            children2,
        )
        @test !CS.same_logs(l1, l2)
    end

    function create_table_row(table::CS.TableSetup, vec::Vector)
        row = Vector{CS.TableEntry}(undef, length(vec))
        c = 1
        for (col, val) in zip(table.cols, vec)
            row[c] = CS.TableEntry(col.id, convert(col.type, val))
            c += 1
        end
        return row
    end

    @testset "TableLogging" begin
        table = CS.TableSetup(
            [
                CS.TableCol(:open_nodes, "#Open", Int, 10, :center),
                CS.TableCol(:closed_nodes, "#Closed", Int, 10, :center),
                CS.TableCol(:incumbent, "Incumbent", Float64, 20, :center),
                CS.TableCol(:best_bound, "Best Bound", Float64, 20, :center),
                CS.TableCol(:duration, "Time [s]", Float64, 10, :center),
            ],
            Dict(:min_diff_duration => 5.0),
        )
        table_header = CS.get_header(table)
        lines = split(table_header, "\n")
        @test length(lines) == 2
        @test length(lines[2]) == sum([c.width + 2 for c in table.cols])
        @test occursin("#Open", lines[1])
        @test occursin("#Closed", lines[1])

        table_row = create_table_row(table, [0, 0, 1.0, 1.0, 0.203])
        line = CS.get_row(table, table_row)
        @test length(line) == sum([c.width + 2 for c in table.cols])
        line_split = split(line, r"\s+")
        # +2 for first and last empty
        @test length(line_split) == length(table.cols) + 2
        @test line_split[2] == "0"
        # only precision 2
        @test line_split[end - 1] == "0.20"

        table_row = create_table_row(table, [1000000000000, 1000000000000, 1.0, 1.0, 0.203])
        line = CS.get_row(table, table_row)
        @test length(line) == sum([c.width + 2 for c in table.cols])
        line_split = split(line, r"\s+")
        # +2 for first and last empty
        @test length(line_split) == length(table.cols) + 2
        @test line_split[2] == ">>"
        @test line_split[3] == ">>"

        # duration too long
        table_row = create_table_row(table, [1, 2, 1.0, 1.0, 10000000000.203])
        line = CS.get_row(table, table_row)
        @test length(line) == sum([c.width + 2 for c in table.cols])
        line_split = split(line, r"\s+")
        # +2 for first and last empty
        @test length(line_split) == length(table.cols) + 2
        @test line_split[end - 1] == ">>"

        # better precision for bound
        table_row = create_table_row(table, [1, 2, 1.0, 0.000004, 0.203])
        line = CS.get_row(table, table_row)
        @test length(line) == sum([c.width + 2 for c in table.cols])
        line_split = split(line, r"\s+")
        # +2 for first and last empty
        @test length(line_split) == length(table.cols) + 2
        @test line_split[end - 2] == "0.000004"

        @assert CS.push_to_table!(
            table;
            open_nodes = 1,
            closed_nodes = 1,
            incumbent = 1.0,
            best_bound = 1.0,
            duration = 0.1,
        )
        # don't print because the duration difference is less than 5sec
        @assert !CS.push_to_table!(
            table;
            open_nodes = 1,
            closed_nodes = 1,
            incumbent = 1.0,
            best_bound = 1.0,
            duration = 0.2,
        )

        # Incumbent precision too high
        table = CS.TableSetup([
            CS.TableCol(:open_nodes, "#Open", Int, 10, :center),
            CS.TableCol(:closed_nodes, "#Closed", Int, 10, :left),
            CS.TableCol(:incumbent, "Incumbent", Float64, 10, :center), # will get increased to 11
            CS.TableCol(:best_bound, "Best Bound", Float64, 10, :center),
            CS.TableCol(:duration, "Time [s]", Float64, 10, :right),
        ])
        table_row = create_table_row(table, [1, 2, 100000000.02, 0.000004, 0.203])
        line = CS.get_row(table, table_row)
        println("line: $line")
        @test length(line) == sum([c.width + 2 for c in table.cols])
        line_split = split(line, r"\s+")
        # +2 for first and last empty
        @test length(line_split) == length(table.cols) + 2
        @test line_split[4] == "100000000.0"
        @assert CS.push_to_table!(
            table;
            open_nodes = 1,
            closed_nodes = 1,
            incumbent = 1.0,
            best_bound = 1.0,
            duration = 0.1,
        )
        @assert CS.push_to_table!(
            table;
            open_nodes = 1,
            closed_nodes = 1,
            incumbent = 1.0,
            best_bound = 1.0,
            duration = 0.2,
        )

        # need to increase size for #Open
        table = CS.TableSetup([
            CS.TableCol("Open", Int),
            CS.TableCol("#Closed", Int, 10),
            CS.TableCol(:incumbent, "Incumbent", Float64, 20),
            CS.TableCol(:best_bound, "Best Bound", Float64, 20, :center),
            CS.TableCol(:duration, "Time [s]", Float64, 10, :center),
        ])
        table_header = CS.get_header(table)
        lines = split(table_header, "\n")
        @test length(lines) == 2
        @test length(lines[2]) == sum([c.width + 2 for c in table.cols])
        @test length(lines[2]) > sum([1, 10, 20, 20, 10])
        @test occursin("Open", lines[1])
        @test occursin("Closed", lines[1])
    end

    @testset "Traverse" begin
        com = CS.ConstraintSolverModel()
        com.options.traverse_strategy = :DFS
        com.traverse_strategy = Val(:DFS)
        com.sense = MOI.MIN_SENSE
        com.backtrack_vec = Vector{CS.BacktrackObj{Float64}}()
        bounds = [0.4, 0.15, 0.15, 0.1, 0.1]
        depths = [3, 2, 1, 2, 2]
        bo = CS.BacktrackObj(com)
        for i in 1:length(bounds)
            bo.status = :Open
            bo.idx = i
            bo.best_bound = bounds[i]
            bo.depth = depths[i]
            push!(com.backtrack_vec, bo)
            CS.add2priorityqueue(com, com.backtrack_vec[end])
        end
        order = [1, 4, 5, 2, 3]
        for i in 1:length(bounds)
            found, bo = CS.get_next_node(com, com.backtrack_vec, true)
            @test found
            @test bo.idx == order[i]
            CS.close_node!(com, bo.idx)
        end

        com = CS.ConstraintSolverModel()
        com.traverse_strategy = Val(:BFS)
        com.options.traverse_strategy = :BFS
        com.sense = MOI.MIN_SENSE
        com.backtrack_vec = Vector{CS.BacktrackObj{Float64}}()
        bounds = [0.4, 0.15, 0.15, 0.1, 0.1]
        depths = [3, 2, 1, 2, 2]
        bo = CS.BacktrackObj(com)
        for i in 1:length(bounds)
            bo.status = :Open
            bo.idx = i
            bo.best_bound = bounds[i]
            bo.depth = depths[i]
            push!(com.backtrack_vec, bo)
            CS.add2priorityqueue(com, com.backtrack_vec[end])
        end

        order = [4, 5, 2, 3, 1]
        for i in 1:length(bounds)
            found, bo = CS.get_next_node(com, com.backtrack_vec, true)
            @test found
            @test bo.idx == order[i]
            bo.status = :Closed
        end
    end
end