text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "PTSCOTCH"
version = v"6.1.3"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://gitlab.inria.fr/scotch/scotch/-/archive/v$(version)/scotch-v$(version).tar.gz","4e54f056199e6c23d46581d448fcfe2285987e5554a0aa527f7931684ef2809e"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/scotch*
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/native_build.patch"
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/Makefile.patch"
if [[ "${target}" == *apple* || "${target}" == *freebsd* ]]; then
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/OSX_FreeBSD.patch"
fi
if [[ "${target}" == *mingw* ]]; then
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/Windows.patch"
fi
cd src
make ptscotch
make ptesmumps
cp ../lib/libpt* ${libdir}
cp ../include/p* ${includedir}
install_license ../LICENSE_en.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; exclude=Sys.iswindows)

# The products that we will ensure are always built
products = [
    LibraryProduct("libptscotcherr", :libptscotcherr),
    LibraryProduct("libptscotcherrexit", :libptscotcherrexit),
    LibraryProduct("libptscotchparmetis", :libptscotchparmetis),
    LibraryProduct("libptscotch", :libptscotch, dont_dlopen=true),
    LibraryProduct("libptesmumps", :libptesmumps)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a")),
    Dependency("MPICH_jll"),
    Dependency("SCOTCH_jll", v"6.1.3")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"9.1.0", julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
import Base.exp

abstract type Function end

mutable struct Variable
    data::Array
    grad::Union{Array,Nothing}
    creator::Union{Function,Nothing}
    Variable(data) = new(data, nothing, nothing)
end
function backward(v::Variable)::Nothing
    if isnothing(v.grad)
        v.grad = ones(size(v.data))
    end
    funcs::Array{Function} = [v.creator]
    while !isempty(funcs)
        f::Function = pop!(funcs) # No need to check for nothingness at the first time
        gys::Array{Array} = [output.grad for output in f.outputs]
        gxs::Union{Tuple{Vararg{Array}},Array} = backward(f, gys...)
        if !isa(gxs, Tuple)
            gxs = (gxs,)
        end
        for (x, gx) in zip(f.inputs, gxs)
            if isnothing(x.grad)
                x.grad = gx
            else
                x.grad += gx
            end
            if !isnothing(x.creator)
                push!(funcs, x.creator)
            end
        end
    end
end
function cleargrad(v::Variable)
    v.grad = nothing
end

function _function(
    f::Function, forward::Core.Function, inputs::Variable...
)::Union{Array{Variable},Variable}
    xs::Array{Array} = [x.data for x in inputs]
    ys::Union{Tuple{Vararg{Array}},Array} = forward(xs...)
    if !isa(ys, Tuple)
        ys = (ys,)
    end
    outputs = [Variable(y) for y in ys]
    for output in outputs
        output.creator = f
    end
    f.inputs = inputs
    f.outputs = outputs
    return length(outputs) > 1 ? outputs : outputs[1]
end

mutable struct Add <: Function
    inputs::Union{Tuple{Vararg{Variable}},Nothing}
    outputs::Union{Array{Variable},Nothing}
    Add() = new()
end
(f::Add)(inputs::Variable...)::Variable = _function(f, (x0, x1) -> x0 + x1, inputs...)
add(x0::Variable, x1::Variable)::Variable = Add()(x0, x1)
backward(::Add, gy::Array)::Tuple{Array,Array} = gy, gy

mutable struct Square <: Function
    inputs::Union{Tuple{Vararg{Variable}},Nothing}
    outputs::Union{Array{Variable},Nothing}
    Square() = new()
end
(f::Square)(input::Variable)::Variable = _function(f, x -> x .^ 2, input)
square(x::Variable)::Variable = Square()(x)
backward(f::Square, gy::Array)::Array = 2 * f.inputs[1].data .* gy

ext: jl


text:
##
## L1 regularization
##
@with_kw mutable struct regularizer_l1 <: Regularizer
    λ::Float64
end

function wholespace_manifold(reg::regularizer_l1, x)
    return l1Manifold(ones(size(x)))
end

## 0th order
g(reg::regularizer_l1, x) = reg.λ * norm(x, 1)

## 1st order
function prox_αg!(reg::regularizer_l1, res, x, α)
    res .= softthresh.(x, reg.λ * α)
    return l1Manifold(map(x -> abs(x) > 0, res))
end
function prox_αg(reg::regularizer_l1, x, α)
    res = softthresh.(x, reg.λ * α)
    return res, l1Manifold(map(x -> abs(x) > 0, res))
end

∇M_g!(reg::regularizer_l1, ::l1Manifold, ∇M_g, x) = (@. ∇M_g = reg.λ * sign.(x))
function ∇M_g(reg::regularizer_l1, M::l1Manifold, x)
    ∇M_g = zeros(representation_size(M))
    return ∇M_g!(reg, M, ∇M_g, x)
end

## 2nd order
∇²M_g_ξ(::regularizer_l1, M::l1Manifold, x, ξ) = zeros(size(ξ))
∇²M_g_ξ!(::regularizer_l1, M::l1Manifold, res, x, ξ) = (res .= 0)

#
### Optimality status
#
function model_g_subgradient!(model, regularizer::regularizer_l1, M, x)
    n = size(x, 1)
    λ = regularizer.λ

    d = manifold_dimension(M)

    @variable(model, ḡ_normal[1:n-d])
    @constraint(model, -λ .<= ḡ_normal .<= λ)

    return ḡ_normal
end

function build_subgradient_from_normalcomp(regularizer::regularizer_l1, M, x, ḡ_normal)
    n = size(x, 1)
    ḡ = Vector(undef, n)
    λ = regularizer.λ
    ind_norm = 1
    for i in 1:n
        if M.nnz_coords[i]
            ḡ[i] = sign(x[i]) * λ
        else
            ḡ[i] = ḡ_normal[ind_norm]
            ind_norm += 1
        end
    end
    return ḡ
end

function build_normalcomp_from_subgradient(::regularizer_l1, M, x, ḡ)
    return ḡ[.!M.nnz_coords]
end


text:
# using BinDeps
# @BinDeps.setup

# libsmilejl = library_dependency("libsmilejl.so", os=:Unix)

# provides(Binaries,
# 	URI("http://dl.bintray.com/tawheeler/generic/libsmilejl.so/libsmilejl.so"),
# 	libsmilejl, os=:Unix, unpacked_dir="libsmilejl")

# @BinDeps.install [:libsmilejl => :libsmilejl]

downloads = joinpath(dirname(@__FILE__), "downloads")
println(downloads)

if !isdir(downloads)
	mkdir(downloads)
end

destfile = joinpath(downloads, "libsmilejl.so")
println(destfile)

download("https://bintray.com/artifact/download/tawheeler/generic/generic/libsmilejl.so", destfile)


text:
module LB

include("module_common.jl")

const _libcloud_lb_types = PyCall.PyNULL()
const _libcloud_lb_providers = PyCall.PyNULL()
const _libcloud_lb_base = PyCall.PyNULL()

immutable LBDriver
    o::PyObject

    function LBDriver(provider::Compat.String, args...; kwargs...)
        cls = _libcloud_lb_providers[:get_driver](provider)
        LBDriver(cls(args...; kwargs...))
    end
    LBDriver(driver::PyObject) = new(driver)
end
show(io::IO, c::LBDriver) = print(io, c.o[:__str__]())

immutable LoadBalancer
    o::PyObject

    id::Nullable{Compat.String}
    name::Compat.String
    state::Int
    ip::Compat.String
    port::Nullable{Int}
    driver::LBDriver
    extra::Dict

    function LoadBalancer(o::PyObject)
        new(o, o[:id], o[:name], o[:state], o[:ip], o[:port], LBDriver(o[:driver]), o[:extra])
    end
end
PyObject(o::LoadBalancer) = o.o
convert(::Type{LoadBalancer}, o::PyObject) =  LoadBalancer(o)
show(io::IO, o::LoadBalancer) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_lb_base, "LoadBalancer") LoadBalancer

immutable Member
    o::PyObject

    id::Nullable{Compat.String}
    ip::Nullable{Compat.String}
    port::Nullable{Compat.String}
    balancer::Nullable{LoadBalancer}
    extra::Dict

    function Member(o::PyObject)
        balancer = (o[:balancer] == nothing) ? Nullable{LoadBalancer}() : Nullable(LoadBalancer(o[:balancer]))
        new(o, o[:id], o[:ip], o[:port], balancer, o[:extra])
    end
end
PyObject(o::Member) = o.o
convert(::Type{Member}, o::PyObject) =  Member(o)
show(io::IO, o::Member) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_lb_base, "Member") Member

const _lbdriver_fns = [
    :list_protocols, :list_balancers, :create_balancer, :destroy_balancer,
    :get_balancer, :update_balancer, :balancer_attach_compute_node, :balancer_attach_member,
    :balancer_detach_member, :balancer_list_members, :list_supported_algorithms
]

const _lb_fns = [
    :attach_compute_node, :attach_member, :detach_member, :list_members, :destroy
]

for f in union(Set(_lbdriver_fns), Set(_base_driver_fns))
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_lb_base, "Driver", $sf) $(f)(lb::LBDriver, args...; kwargs...) = lb.o[$(sf)](args..., kwargs...)
end

for f in _lb_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_lb_base, "LoadBalancer", $sf) $(f)(lb::LoadBalancer, args...; kwargs...) = lb.o[$(sf)](args..., kwargs...)
end

# Initialize load balancer
function __init__()
    copy!(_libcloud_lb_types, pyimport("libcloud.loadbalancer.types"))
    copy!(_libcloud_lb_providers, pyimport("libcloud.loadbalancer.providers"))
    copy!(_libcloud_lb_base, pyimport("libcloud.loadbalancer.base"))

    global const LBProvider = pywrap(_libcloud_lb_types[:Provider])
    global const State = pywrap(_libcloud_lb_types[:State])
    global const MemberCondition = pywrap(_libcloud_lb_types[:MemberCondition])
    global const Algorithm = pywrap(_libcloud_lb_base[:Algorithm])
    global const DefaultAlgorithm = Algorithm.ROUND_ROBIN

    _map_types(_libcloud_lb_base, (Member, LoadBalancer))
end

# types
export LBProvider, LBDriver, State, MemberCondition, Algorithm, DefaultAlgorithm, LoadBalancer, Member
# lb driver functions
export list_protocols, list_balancers, create_balancer, destroy_balancer,
       get_balancer, update_balancer, balancer_attach_compute_node, balancer_attach_member,
       balancer_detach_member, balancer_list_members, list_supported_algorithms
# lb functons
export attach_compute_node, attach_member, detach_member, list_members, destroy

end # module LoadBalancer


text:
using Test
using SparseIR

@testset "kernel.jl" begin
    @testset "accuracy with K = $K" for K in (
        LogisticKernel(9.0),
        RegularizedBoseKernel(8.0),
        LogisticKernel(120_000.0),
        RegularizedBoseKernel(127_500.0),
        SparseIR.get_symmetrized(LogisticKernel(40_000.0), -1),
        SparseIR.get_symmetrized(RegularizedBoseKernel(35_000.0), -1),
    )
        T = Float32
        T_x = Float64

        rule = convert(SparseIR.Rule{T}, SparseIR.legendre(10))
        hints = SparseIR.sve_hints(K, eps(T_x))
        gauss_x = SparseIR.piecewise(rule, SparseIR.segments_x(hints))
        gauss_y = SparseIR.piecewise(rule, SparseIR.segments_y(hints))
        ϵ = eps(T)
        tiny = floatmin(T) / ϵ

        result = SparseIR.matrix_from_gauss(K, gauss_x, gauss_y)
        result_x = SparseIR.matrix_from_gauss(
            K, convert(SparseIR.Rule{T_x}, gauss_x), convert(SparseIR.Rule{T_x}, gauss_y)
        )
        magn = maximum(abs, result_x)

        @test result ≈ result_x atol = 2magn * ϵ rtol = 0
        reldiff = @. ifelse(abs(result) < tiny, 1, result / result_x)
        @test all(x -> isapprox(x, 1, atol=100ϵ, rtol=0), reldiff)
    end

    @testset "singularity with Λ = $Λ" for Λ in (10, 42, 10_000), x in 2rand(10) .- 1
        K = RegularizedBoseKernel(Λ)
        @test K(x, 0.0) ≈ 1 / Λ
    end
end


text:
@time include(joinpath("..", "required.jl"))

## Create glossary dataframe

# Landcover data
lcgloss = [
    "lc1" "bare" "Bare land cover (%)"
    "lc2" "crops" "Crops land cover (%)"
    "lc3" "grass" "Grass land cover (%)"
    "lc4" "moss" "Moss land cover (%)"
    "lc5" "shrub" "Shrub land cover (%)"
    "lc6" "snow" "Snow land cover (%)"
    "lc7" "tree" "Tree land cover (%)"
    "lc8" "urban" "Urban land cover (%)"
    "lc9" "water_permanent" "Permanent water land cover (%)"
    "lc10" "water_seasonal" "Seasonal water land cover (%)"
]
lcdf = DataFrame(;
    variable=lcgloss[:, 1],
    type="landcover",
    full_name=lcgloss[:, 2],
    description=lcgloss[:, 3],
)
# Climate data
wcgloss = [
    "wc1" "temperature_annual_mean" "Annual Mean Temperature (°C)"
    "wc2" "temperature_diurnal_range" "Mean Diurnal Range (Mean of monthly (max temp - min temp)) (°C)"
    "wc3" "temperature_isothermality" "Isothermality (BIO2/BIO7) (* 100) (°C)"
    "wc4" "temperature_seasonality" "Temperature Seasonality (standard deviation *100) (°C)"
    "wc5" "temperature_max_warmest_month" "Max Temperature of Warmest Month (°C)"
    "wc6" "temperature_min_coldest_month" "Min Temperature of Coldest Month (°C)"
    "wc7" "temperature_annual_range" "Temperature Annual Range (BIO5-BIO6) (°C)"
    "wc8" "temperature_mean_wettest_quarter" "Mean Temperature of Wettest Quarter (°C)"
    "wc9" "temperature_mean_driest_quarter" "Mean Temperature of Driest Quarter (°C)"
    "wc10" "temperature_mean_warmest_quarter" "Mean Temperature of Warmest Quarter (°C)"
    "wc11" "temperature_mean_coldest_quarter" "Mean Temperature of Coldest Quarter (°C)"
    "wc12" "precipitation_annual_mean" "Annual Precipitation (mm)"
    "wc13" "precipitation_wettest_month" "Precipitation of Wettest Month (mm)"
    "wc14" "precipitation_driest_month" "Precipitation of Driest Month (mm)"
    "wc15" "precipitation_seasonality" "Precipitation Seasonality (Coefficient of Variation) (mm)"
    "wc16" "precipitation_wettest_quarter" "Precipitation of Wettest Quarter (mm)"
    "wc17" "precipitation_driest_quarter" "Precipitation of Driest Quarter (mm)"
    "wc18" "precipitation_warmest_quarter" "Precipitation of Warmest Quarter (mm)"
    "wc19" "precipitation_coldest_quarter" "Precipitation of Warmest Quarter (mm)"
]
wcdf = DataFrame(;
    variable=wcgloss[:, 1],
    type="climate",
    full_name=wcgloss[:, 2],
    description=wcgloss[:, 3],
)

# Species data
@load joinpath("data", "jld2", "spenames.jld2") spenames specommon
spenames
specommon
spdf = DataFrame(;
    variable=string.("sp", eachindex(spenames)),
    type="species",
    full_name=spenames,
    description=specommon,
)

# Create glossary
glossdf = vcat(lcdf, wcdf, spdf)

# Export to CSV
CSV.write(joinpath("data", "proc", "glossary.csv"), glossdf)

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
### A Pluto.jl notebook ###
# v0.19.4

#> [frontmatter]
#> title = "ExoFinder.jl"
#> description = "Let's find some worlds!"

using Markdown
using InteractiveUtils

# ╔═╡ f19b358c-8506-11ec-252c-c39dcd644d06
begin
	import Pkg
	Pkg.activate(Base.current_project())

	using AstroImages, PlutoUI, Plots
	using MarkdownLiteral: @mdx
end

# ╔═╡ 2d5365c9-c7d9-4d14-b074-f77b39baec02
using AstroImages.AstroAngles

# ╔═╡ cc2aaba1-4dfe-414f-8ace-71075e2f3bf5
using HTTP

# ╔═╡ ad40234c-474c-489e-b328-73fa5ecfb6e2
using HTTP.URIs

# ╔═╡ acdf51db-09e0-4e4c-b529-2db8030ea57c
using CSV, DataFramesMeta

# ╔═╡ 3845b39a-a637-4d2b-b2b9-f4ac0294f0e9
@mdx """
# ExoFinder 🪐

In this notebook, we will build a simple map showing the locations of known exoplanets in the night sky. The final product will look like this:
"""

# ╔═╡ e58691c1-98ed-4e75-a5bb-e03102e62def
md"""
This demo uses the following packages:
"""

# ╔═╡ 3ea72aeb-356c-4fd1-b228-f0cf0cba2df7
with_terminal() do
	Pkg.status()
end

# ╔═╡ 6207daa3-fde6-4535-9ef3-d1e4a762a14d
@mdx """
## 🗺 Sky map

First, let's set up a way to create a coordinate-aware plot of a patch of the sky. We will accomplish this with the handy [AstroImages.jl](https://github.com/JuliaAstro/AstroImages.jl) package. First, we download a sample image (stored in the [FITS](https://en.wikipedia.org/wiki/FITS) file format):

!!! warning "TODO"
	Update link with Will's [PR](https://github.com/JuliaAstro/AstroImages.jl/pull/30)
"""

# ╔═╡ bb2425be-3e6a-456d-bad2-e665dc7408aa
fname = download("http://data.astropy.org/tutorials/FITS-images/HorseHead.fits")

# ╔═╡ a5d80eed-fb12-4980-82db-800a1c4dba49
img = load(fname);

# ╔═╡ 5fcff0be-3d80-4423-a239-2a00aa376db3
let
	implot(img; wcsticks=true, framestyle=:grid)
	X, Y = world_to_pix(img, [hms"5 41 00"deg, dms"-2 28"deg])
	scatter!([X], [Y]; ms=10, msw=0, color=:cyan)
end

# ╔═╡ 8830d13c-04e7-4333-babc-10bb267993fe
@mdx """
And then just plot it!
"""

# ╔═╡ 8762e1a5-505b-4eed-8863-35171e8cc8e1
implot(img)

# ╔═╡ 50666f3e-b5a7-4fab-86a6-979c8da62693
@mdx """
!!! Question
	What was the reasoning for separating `plot` from `implot` again?
"""

# ╔═╡ 67e9e89e-7442-4d96-b8e2-fadb900a8cc3
@mdx """
In the background, AstroImages.jl calls the generic plotting package [Plots.jl](https://docs.juliaplots.org/stable/), and automatically adds a colorbar, tick marks, labels, and grid lines for us. The relevant transformations from pixel space to world coordinates (in this case RA and Dec) are also automatically applied based on the [WCS](https://www.atnf.csiro.au/people/mcalabre/WCS/) information stored in the image of our fits file.

For more information about this package and its different capabilities, be sure to check out the [demo AstroImages.jl notebook]()!

!!! warning "TODO"
	Add a demo notebook from Will?
"""

# ╔═╡ d4552736-159b-4c6a-9a97-ef57024333f7
@mdx """
### Coordinate conversion aside

Now that we have a nice way to visualize astronomical data with coordinate-aware transformations, let's compose this with regular plot commands from Plots.jl. If we turn off the axes transformations, we recover the typical pixel space information stored in our FITs file:
"""

# ╔═╡ 6fd6108a-dc97-42ed-a722-bc74c19a0bde
implot(img; wcsticks=false, framestyle=:grid)

# ╔═╡ b8e06c27-5504-41d1-b07c-c5259ff2681f
@mdx """
By default, the origin is located in the bottom-left corner and indices increase to the right and up. Since `implot` returns a standard Plots.jl figure object, we can plot a test point on top of it [in the usual way](https://docs.juliaplots.org/latest/basics/):
"""

# ╔═╡ 12b9a88f-2f9a-4744-ba69-496545a0dac0
let
	p = implot(img; wcsticks=false, framestyle=:grid)
	scatter!(p, [500], [400]; ms=10, msw=0, color=:cyan)
end

# ╔═╡ e0a44c34-69f8-4037-82a8-60b196b87382
@mdx """
!!! Question
	`wcticks = false` seems to override `grid`. Is there a way to plot in pixel space while still having the grid lines enabled?
"""

# ╔═╡ aa60557e-faf0-4058-96fd-128ca093c773
@mdx """
We can go the other direction (from world cordinate space to pixel space) in the following two stages:

1. First, we use [AstroAngles.jl](https://github.com/JuliaAstro/AstroAngles.jl) to conveniently convert RA and Dec formatted input to their equivalent values in degrees for us. This package supports a wide range of formats, which are shown [here](https://github.com/JuliaAstro/AstroAngles.jl#usage).
1. Next, we pass these values to the `world_to_pix` function exported by `AstroImages.jl` to make the WCS transformations in our final plot.

!!! Question
	Are degrees used internally for all cases?
"""

# ╔═╡ ba4e4e4e-7f0a-4590-981c-619b53fd0bec
let
	implot(img; wcsticks=true, framestyle=:grid)
	X, Y = world_to_pix(img, [hms"5 41 00"deg, dms"-2 28 00"deg])
	scatter!([X], [Y]; ms=10, msw=0, color=:cyan)
end

# ╔═╡ f62bae29-e31a-416e-b844-9720a5ef57f2
@mdx """
!!! Question
	Is it better to do `using AstroAngles` explicitly here?
"""

# ╔═╡ a6368cc9-b7c6-4ffd-a9ba-5ec33be3cb2b
@mdx """
Now that we have a way to plot a single arbitray point, let's extend this to a whole collection of known exoplanet coordinates.
"""

# ╔═╡ d03ffae4-5a15-448f-a47b-e850049efe80
@mdx """
## 🎯 Exoplanet locations

The [NASA Exoplanet Archive](https://exoplanetarchive.ipac.caltech.edu/) stores an updated list of known exoplanets, along with additional information about its host star and orbital parameters. As of this writing, there are 5,000+ confirmed detections, and fortunately their is an API to query all of this information!

The archive provides a Table Access Protocol [(TAP)](https://exoplanetarchive.ipac.caltech.edu/docs/TAP/usingTAP.html) service to query the data using an astronomy specific extension of SQL known as Astronomical Data Query Language [(ADQL)](https://www.ivoa.net/documents/ADQL/). This essentially boils down to pasting a query into a url, and then pulling it down with a `GET` request, which we accomplish with [HTTP.jl](https://juliaweb.github.io/HTTP.jl/stable/):
"""

# ╔═╡ 1319c8bf-ea90-469a-8433-5c3b66b1af07
q = """
select top 10 hostname, pl_name, tic_id
from pscomppars
"""

# ╔═╡ ef38432c-0ec4-46b7-9444-9321180729d9
query = "query=" * escapeuri(q) * "&format=csv"

# ╔═╡ 4d073b62-c25b-4cd9-be89-d87144f2bfdb
@mdx """
!!! note "Note"
	We use the `escapeuri` function exported by [URIs.jl](https://docs.juliahub.com/URIs/eec2u/1.3.0/#Tutorial) to convert the spaces and other special characters into valid [URI](https://en.wikipedia.org/wiki/Uniform_Resource_Identifier) characters.
"""

# ╔═╡ 4581038a-fb53-49c7-a85f-60eb153b6f25
df = CSV.read(
	HTTP.get("https://exoplanetarchive.ipac.caltech.edu/TAP/sync"; query).body,
	DataFrame,
)

# ╔═╡ 403c435d-4d54-49d6-a50f-9f5362ae96d9
@mdx """
We now have the data in a convenient table format (provided by [DataFrames.jl](https://dataframes.juliadata.org/stable/)) for our queried exoplanets.
"""

# ╔═╡ 6db3d93a-0c35-4c27-ace7-1fb44966d864
@mdx """
### Data extraction

Whooo
"""

# ╔═╡ 127338cb-b917-4e2d-8ba1-3ed045c799a4
@mdx """
# Notebook setup 📦
"""

# ╔═╡ fcceea3e-db8f-4853-af49-240d66d54377
TableOfContents()

# ╔═╡ Cell order:
# ╟─3845b39a-a637-4d2b-b2b9-f4ac0294f0e9
# ╟─5fcff0be-3d80-4423-a239-2a00aa376db3
# ╟─e58691c1-98ed-4e75-a5bb-e03102e62def
# ╟─3ea72aeb-356c-4fd1-b228-f0cf0cba2df7
# ╟─6207daa3-fde6-4535-9ef3-d1e4a762a14d
# ╠═bb2425be-3e6a-456d-bad2-e665dc7408aa
# ╠═a5d80eed-fb12-4980-82db-800a1c4dba49
# ╟─8830d13c-04e7-4333-babc-10bb267993fe
# ╠═8762e1a5-505b-4eed-8863-35171e8cc8e1
# ╟─50666f3e-b5a7-4fab-86a6-979c8da62693
# ╟─67e9e89e-7442-4d96-b8e2-fadb900a8cc3
# ╟─d4552736-159b-4c6a-9a97-ef57024333f7
# ╠═6fd6108a-dc97-42ed-a722-bc74c19a0bde
# ╟─b8e06c27-5504-41d1-b07c-c5259ff2681f
# ╠═12b9a88f-2f9a-4744-ba69-496545a0dac0
# ╟─e0a44c34-69f8-4037-82a8-60b196b87382
# ╟─aa60557e-faf0-4058-96fd-128ca093c773
# ╠═ba4e4e4e-7f0a-4590-981c-619b53fd0bec
# ╠═2d5365c9-c7d9-4d14-b074-f77b39baec02
# ╟─f62bae29-e31a-416e-b844-9720a5ef57f2
# ╟─a6368cc9-b7c6-4ffd-a9ba-5ec33be3cb2b
# ╟─d03ffae4-5a15-448f-a47b-e850049efe80
# ╠═cc2aaba1-4dfe-414f-8ace-71075e2f3bf5
# ╠═ad40234c-474c-489e-b328-73fa5ecfb6e2
# ╠═1319c8bf-ea90-469a-8433-5c3b66b1af07
# ╠═ef38432c-0ec4-46b7-9444-9321180729d9
# ╟─4d073b62-c25b-4cd9-be89-d87144f2bfdb
# ╠═4581038a-fb53-49c7-a85f-60eb153b6f25
# ╟─403c435d-4d54-49d6-a50f-9f5362ae96d9
# ╟─6db3d93a-0c35-4c27-ace7-1fb44966d864
# ╠═acdf51db-09e0-4e4c-b529-2db8030ea57c
# ╟─127338cb-b917-4e2d-8ba1-3ed045c799a4
# ╠═fcceea3e-db8f-4853-af49-240d66d54377
# ╠═f19b358c-8506-11ec-252c-c39dcd644d06


text:
function getDiffusionKernel{T <: FloatingPoint}(:: Type{T})
        nType = T == Float64 ? "double" : "float"
        nOne = T == Float64 ? "1.0" : "1.0f"

        return "
        #if defined(cl_khr_fp64)  // Khronos extension available?
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
        #elif defined(cl_amd_fp64)  // AMD extension available?
        #pragma OPENCL EXTENSION cl_amd_fp64 : enable
        #endif
        #define number $nType
        #define number8 $(nType)8

        #define Conc(x,y) a[y*D2 + x]
        #define Pot(x,y) b[y*D2 + x]
        #define Flow(x,y) f[y*D2 + x]

        #define P_move(x,y) out[y*D2 + x]

        __kernel void flow(
                      __global const number *b,
                      __global number8 *f,
                      const int D1,
                      const int D2) {

            int i = get_global_id(0);
            int j = get_global_id(1);
            int west = j-1;
            int east = j+1;
            int north = i+1;
            int south = i-1;

            if(j == 0)
                west = D2 - 1;
            if(j == D2 - 1)
                east = 0;
            if(i == D1 - 1)
                north = 0;
            if(i == 0)
                south = D1 - 1;

            const number p = Pot(i,j);

            number8 ge;
            ge.s0 = Pot(south,west);
            ge.s1 = Pot(south,j   );
            ge.s2 = Pot(south,east);
            ge.s3 = Pot(i    ,east);
            ge.s4 = Pot(north,east);
            ge.s5 = Pot(north,j   );
            ge.s6 = Pot(north,west);
            ge.s7 = Pot(i    ,west);

            ge = ge - p;

            ge = -1 * ge / (1-exp(ge));

            //Check for NaN
            const number8 ones = (number8)($nOne);
            ge = select(ge, ones, isnan(ge));
            Flow(i,j) = ge / 8;

        }

        static number sum(const number8 n) {
            return n.s0 + n.s1 +n.s2 + n.s3 + n.s4 + n.s5 + n.s6 + n.s7;
        }

        __kernel void diffusion(
                      __global const number *a,
                      __global const number8 *f,
                      __global number *out,
                      const int D1,
                      const int D2) {

            int i = get_global_id(0);
            int j = get_global_id(1);
            int west = j-1;
            int east = j+1;
            int north = i+1;
            int south = i-1;

            if(j == 0)
                west = D2 - 1;
            if(j == D2 - 1)
                east = 0;
            if(i == D1 - 1)
                north = 0;
            if(i == 0)
                south = D1 - 1;

            // Get the flow out of cell ij

            number outflow =  Conc(i,j) * sum(Flow(i,j));

            // Calculate the inflow based on the outflow from other cells into this on.

            number inflow  =    Conc(south, west) * Flow(south, west).s4 +
                                Conc(south, j   ) * Flow(south, j   ).s5 +
                                Conc(south, east) * Flow(south, east).s6 +
                                Conc(i    , east) * Flow(i    , east).s7 +
                                Conc(north, east) * Flow(north, east).s0 +
                                Conc(north, j   ) * Flow(north, j   ).s1 +
                                Conc(north, west) * Flow(north, west).s2 +
                                Conc(i    , west) * Flow(i    , west).s3 ;

            // Inflow - outflow = change

            P_move(i,j) = inflow - outflow;
    }
"
end

function diffusionCL!{T <: FloatingPoint}(
    a_buff, b_buff,
    out_buff,
    d1 :: Int64, d2 :: Int64,
    ctx, queue, program, ::Type{T})

    flow_buff = Buffer(T, ctx, :rw, d1 * d2 * 8)

    kd = Kernel(program, "diffusion")
    kf = Kernel(program, "flow")

    call(queue, kf, (d1,d2), nothing, b_buff, flow_buff, int32(d1), int32(d2))
    call(queue, kd, (d1,d2), nothing, a_buff, flow_buff, out_buff, int32(d1), int32(d2))
end


text:
"""
module
code adapted from https://discourse.julialang.org/t/custom-subject-in-rocket-jl-for-mouse-events-from-glfw/65133/3
it is design to help processing data from
    -GLFW.SetCursorPosCallback(window, (_, x, y) -> println("cursor: x, y")) and  for example : cursor: 29.0, 469.0  types   Float64  Float64
    -GLFW.SetMouseButtonCallback(window, (_, button, action, mods) -> println("button action"))  for example types MOUSE_BUTTON_1 PRESS   GLFW.MouseButton  GLFW.Action
The main function is to mark the interaction of the mouse to be saved in appropriate mask and be rendered onto the screen
so we modify the data that is the basis of the mouse interaction mask  and we pass the data on so appropriate part of the texture would be modified to be displayed on screen

"""
module ReactOnMouseClickAndDrag
using Logging, Parameters, Rocket,Setfield, GLFW, ModernGL,  ..ForDisplayStructs, ..TextureManag,  ..OpenGLDisplayUtils
using  Dates, Parameters,  ..DataStructs,  ..StructsManag
export registerMouseClickFunctions
export reactToMouseDrag

"""
configuting Rocket on Subscribe so we get custom handler of input as we see we still need to define actor
"""
# function Rocket.on_subscribe!(handler::MouseCallbackSubscribable, actor::SyncActor{Any, ActorWithOpenGlObjects})
#     return subscribe!(handler.subject, actor)
# end

function Rocket.on_subscribe!(handler::MouseCallbackSubscribable, actor::SyncActor{Any, ActorWithOpenGlObjects})

    return subscribe!(handler.subject, actor)
end

"""
we define how handler should act on the subject - observable so it will pass event onto subject - here we have 2 events that we want to be ready for - mouse button press
example of possible inputs that we would be intrested in
for example : cursor: 29.0, 469.0  types   Float64  Float64
for example  MOUSE_BUTTON_1 PRESS    types GLFW.MouseButton  GLFW.Action
             MOUSE_BUTTON_1 RELEASE  types GLFW.MouseButton  GLFW.Action
We get two overloads so we will be able to respond with single handler to both mouse click and mouse position
Enum GLFW.Action:
RELEASE = 0
PRESS = 1
REPEAT = 2
Enum GLFW.MouseButton:
MOUSE_BUTTON_1 = 0
MOUSE_BUTTON_2 = 1

experiments show that max x,y in window is both 600 if window width and height is 600
 so in order to specify weather we are over aour quad we need to know how big is primary quad -
  defaoul it is occupying 100% of y axis and first left 80% of x axis
  hence we can calculate max height to equal the height of the window
"""
function (handler::MouseCallbackSubscribable)( a, x::Float64, y::Float64)
  point = CartesianIndex(Int(x),Int(y))
  handler.lastCoordinate = point

  if  (handler.isLeftButtonDown && x>=handler.xmin && x<=handler.xmax && y>=handler.ymin && y<= handler.ymax )
    push!(handler.coordinatesStoreForLeftClicks,point) # putting coordinate to the list it will be processed when context will be ready
        if(!handler.isBusy[])

            #sending mouse position only if all conditions are met
            next!(handler.subject,MouseStruct(
                isLeftButtonDown=handler.isLeftButtonDown
                ,isRightButtonDown=handler.isRightButtonDown
               ,lastCoordinates=  handler.coordinatesStoreForLeftClicks))#sending mouse position only if all conditions are met
            handler.coordinatesStoreForLeftClicks = [point]
        end#if
    #now we also get intrested about  right clicks becouse they help to controll
    elseif(handler.isRightButtonDown)
        next!(handler.subject,MouseStruct(
            isLeftButtonDown=handler.isLeftButtonDown
            ,isRightButtonDown=handler.isRightButtonDown
           ,lastCoordinates= [point] ))#sending mouse position only if all conditions are met

   end#if

end #handler

function (handler::MouseCallbackSubscribable)(a, button::GLFW.MouseButton, action::GLFW.Action,m)

    res=(button==GLFW.MOUSE_BUTTON_1 &&  action==GLFW.PRESS)#so it will stop either as we relese left mouse or click right
    handler.isLeftButtonDown = res
    handler.isRightButtonDown = (button==GLFW.MOUSE_BUTTON_2 &&  action==GLFW.PRESS)

    if(res)
        handler.coordinatesStoreForLeftClicks = [handler.lastCoordinate]
    end #if
end #second handler

"""
we pass coordinate of cursor only when isLeftButtonDown is true and we make it true
if left button is presed down - we make it true if the left button is pressed over image and false if mouse get out of the window or we get information about button release
imageWidth adn imageHeight are the dimensions of textures that we use to display
"""
function registerMouseClickFunctions(window::GLFW.Window
                                    ,stopListening::Base.Threads.Atomic{Bool}
                                    ,calcD::CalcDimsStruct
                                    ,isBusy::Base.Threads.Atomic{Bool}     )

 stopListening[]=true # stoping event listening loop to free the GLFW context

 # calculating dimensions of quad becouse it do not occupy whole window, and we want to react only to those mouse positions that are on main image quad
  mouseButtonSubs = MouseCallbackSubscribable(isLeftButtonDown=false
                                        ,isRightButtonDown=false
                                        ,xmin=Int32(calcD.windowWidthCorr)
                                        ,ymin=Int32(calcD.windowHeightCorr)
                                        ,xmax=Int32(calcD.avWindWidtForMain-calcD.windowWidthCorr)
                                        ,ymax=Int32(calcD.avWindHeightForMain-calcD.windowHeightCorr)
                                        ,coordinatesStoreForLeftClicks=[]
                                        ,lastCoordinate=CartesianIndex(1,1)
                                        ,isBusy=isBusy
                                        ,subject=Subject(MouseStruct  ,scheduler = AsyncScheduler()))

GLFW.SetCursorPosCallback(window, (a, x, y) -> mouseButtonSubs(a,x, y ) )# and  for example : cursor: 29.0, 469.0  types   Float64  Float64
GLFW.SetMouseButtonCallback(window, (a, button, action, mods) ->mouseButtonSubs(a,button, action,mods )) # for example types MOUSE_BUTTON_1 PRESS   GLFW.MouseButton  GLFW.Action

stopListening[]=false # reactivate event listening loop

return mouseButtonSubs

end #registerMouseScrollFunctions

 """
we use mouse coordinate to modify the texture that is currently active for modifications
    - we take information about texture currently active for modifications from variables stored in actor
    from texture specification we take also its id and its properties ...
"""
function reactToMouseDrag(mousestr::MouseStruct, actor::SyncActor{Any, ActorWithOpenGlObjects})
    obj = actor.actor.mainForDisplayObjects
    obj.stopListening[]=true #free GLFW context
    actor.actor.isBusy[]=true# mark that OpenGL is busy
    textureList = actor.actor.textureToModifyVec
    mouseCoords= mousestr.lastCoordinates
    if (!isempty(textureList)  && mousestr.isLeftButtonDown && textureList[1].isEditable)
        texture= textureList[1]
        calcDim =  actor.actor.calcDimsStruct

        # two dimensional coordinates on plane of intrest (current slice)
       mappedCoords =  translateMouseToTexture(texture.strokeWidth
                                                ,mouseCoords
                                                ,actor.actor.calcDimsStruct)

       twoDimDat= actor.actor.currentlyDispDat|> # accessing currently displayed data
       (singSl)-> singSl.listOfDataAndImageNames[singSl.nameIndexes[texture.name]] #accessing the texture data we want to modify

       toSet =  convert( parameter_type(texture) ,actor.actor.valueForMasToSet.value)
       sliceDat= modSlice!(twoDimDat, mappedCoords, convert(twoDimDat.type, toSet )) # modifying data associated with texture

     #  updateTexture(twoDimDat.type,sliceDat, texture,0,0,calcDim.imageTextureWidth,calcDim.imageTextureHeight  )

     singleSliceDat = setproperties(actor.actor.currentlyDispDat,( listOfDataAndImageNames=[ sliceDat ] ))

   updateImagesDisplayed(singleSliceDat
                        ,actor.actor.mainForDisplayObjects
                        ,actor.actor.textDispObj
                        ,actor.actor.calcDimsStruct
                        ,actor.actor.valueForMasToSet  )

        #to enable undoing we just set the point we modified back to 0
        addToforUndoVector(actor, ()-> begin
        modSlice!(twoDimDat, mappedCoords,  convert(twoDimDat.type, 0 ))|> # modifying data associated with texture
        (sliceDat)-> updateTexture(twoDimDat.type,sliceDat.dat, texture,0,0,calcDim.imageTextureWidth,calcDim.imageTextureHeight  )
         basicRender(obj.window) end)

    elseif( mousestr.isRightButtonDown  )
        #we save data about right click position in order to change the slicing plane accordingly
        mappedCoords =  translateMouseToTexture(Int32(1)
                    ,mouseCoords
                    ,actor.actor.calcDimsStruct)
        mappedCorrd= mappedCoords
        if(!isempty(mappedCorrd))
            cartMapped= cartTwoToThree(actor.actor.onScrollData.dataToScrollDims
            ,actor.actor.currentDisplayedSlice
        , mappedCoords[1]   )

            actor.actor.lastRecordedMousePosition =  cartMapped
        end#if
     end #if
    actor.actor.isBusy[]=false # we can do sth in opengl
    obj.stopListening[]=false # reactivete event listening loop
end#..ReactToScroll

"""
given list of cartesian coordinates and some window/ image characteristics - it translates mouse positions
to cartesian coordinates of the texture
strokeWidth - the property connected to the texture marking how thick should be the brush
mouseCoords - list of coordinates of mouse positions while left button remains pressed
calcDims - set of values usefull for calculating mouse position
return vector of translated cartesian coordinates
"""
function translateMouseToTexture(strokeWidth::Int32
                                ,mouseCoords::Vector{CartesianIndex{2}}
                                ,calcD::CalcDimsStruct )::Vector{CartesianIndex{2}}

    filteredList= map(c->CartesianIndex( getNewX(c[1],calcD),  getNewY(c[2] ,calcD)) ,mouseCoords)  |>
    (x)->filter(it->it[1]>0 && it[2]>0 ,x)       # we do not want to try access it in point 0 as julia is 1 indexed
   if(!isempty(filteredList))
    return map(point-> addStrokeWidth(point,Int64(strokeWidth) ) , filteredList) |>  # adding some points around the point of choice so will be better visible
    (matrix)->reduce(vcat, matrix)|># when we added some oints around we got list of lists so now we need to flatten it out
    unique|> # we want only unique elements
    uniq-> filter(it-> it[1]>0 && it[1]<calcD.imageTextureWidth && it[2]>0 && it[2]<calcD.imageTextureHeight  ,uniq)     # as we add new points they may end up getting outside the texture; we need to filter those out
   end #if
   #if we are here we do not have anything meaningfull else to return
   return Vector{CartesianIndex{2}}()
end #translateMouseToTexture

"""
adding the width to the stroke so we will be able to controll how thicly we are painting ...
"""
function addStrokeWidth(point::CartesianIndex{2},strokeW::Int64  )
    return CartesianIndices((-strokeW:strokeW,-strokeW:strokeW)) |> # set of cartesian indices that we will filter ot later
    list->list.+point|> # making coordinates around point of intrest
    added-> filter(x-> ( abs(point[1]- x[1]) + abs(x[2] -point[2]))<strokeW  ,added )# filtering to distant points
end#addStrokeWidth

xx = CartesianIndex(2,2)
xx[1]
"""
helper function for translateMouseToTexture
"""
function getNewX(x::Int,calcD::CalcDimsStruct)::Int
  # first we subtract windowWidthCorr as in window the image do not need to start at the begining  of the window
   return Int64(floor( ((x- calcD.windowWidthCorr)/(calcD.correCtedWindowQuadWidth))*calcD.imageTextureWidth))
end#getNewX

"""
helper function for translateMouseToTexture
"""
function getNewY(y::Int,calcD::CalcDimsStruct)::Int
    Int64(floor(  ((calcD.correCtedWindowQuadHeight-y+calcD.windowHeightCorr)/calcD.correCtedWindowQuadHeight)*calcD.imageTextureHeight)  )

end#getNewY

end #ReactOnMouseClickAndDrag


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct PolicyDefinitionsApi <: SwaggerApi
    client::Swagger.Client
end

"""
Creates or updates a policy definition in a subscription.
This operation creates or updates a policy definition in the given subscription with the given name.
Param: policyDefinitionName::String (required)
Param: parameters::PolicyDefinition (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PolicyDefinition
"""
function _swaggerinternal_policyDefinitionsCreateOrUpdate(_api::PolicyDefinitionsApi, policyDefinitionName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", PolicyDefinition, "/subscriptions/{subscriptionId}/providers/Microsoft.Authorization/policyDefinitions/{policyDefinitionName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "policyDefinitionName", policyDefinitionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyDefinitionsCreateOrUpdate(_api::PolicyDefinitionsApi, policyDefinitionName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsCreateOrUpdate(_api, policyDefinitionName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyDefinitionsCreateOrUpdate(_api::PolicyDefinitionsApi, response_stream::Channel, policyDefinitionName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsCreateOrUpdate(_api, policyDefinitionName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Creates or updates a policy definition in a management group.
This operation creates or updates a policy definition in the given management group with the given name.
Param: policyDefinitionName::String (required)
Param: parameters::PolicyDefinition (required)
Param: api_version::String (required)
Param: managementGroupId::String (required)
Return: PolicyDefinition
"""
function _swaggerinternal_policyDefinitionsCreateOrUpdateAtManagementGroup(_api::PolicyDefinitionsApi, policyDefinitionName::String, parameters, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", PolicyDefinition, "/providers/Microsoft.Management/managementgroups/{managementGroupId}/providers/Microsoft.Authorization/policyDefinitions/{policyDefinitionName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "policyDefinitionName", policyDefinitionName)  # type String
    Swagger.set_param(_ctx.path, "managementGroupId", managementGroupId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyDefinitionsCreateOrUpdateAtManagementGroup(_api::PolicyDefinitionsApi, policyDefinitionName::String, parameters, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsCreateOrUpdateAtManagementGroup(_api, policyDefinitionName, parameters, api_version, managementGroupId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyDefinitionsCreateOrUpdateAtManagementGroup(_api::PolicyDefinitionsApi, response_stream::Channel, policyDefinitionName::String, parameters, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsCreateOrUpdateAtManagementGroup(_api, policyDefinitionName, parameters, api_version, managementGroupId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Deletes a policy definition in a subscription.
This operation deletes the policy definition in the given subscription with the given name.
Param: policyDefinitionName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_policyDefinitionsDelete(_api::PolicyDefinitionsApi, policyDefinitionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/providers/Microsoft.Authorization/policyDefinitions/{policyDefinitionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "policyDefinitionName", policyDefinitionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyDefinitionsDelete(_api::PolicyDefinitionsApi, policyDefinitionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsDelete(_api, policyDefinitionName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyDefinitionsDelete(_api::PolicyDefinitionsApi, response_stream::Channel, policyDefinitionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsDelete(_api, policyDefinitionName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Deletes a policy definition in a management group.
This operation deletes the policy definition in the given management group with the given name.
Param: policyDefinitionName::String (required)
Param: api_version::String (required)
Param: managementGroupId::String (required)
Return: Nothing
"""
function _swaggerinternal_policyDefinitionsDeleteAtManagementGroup(_api::PolicyDefinitionsApi, policyDefinitionName::String, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/providers/Microsoft.Management/managementgroups/{managementGroupId}/providers/Microsoft.Authorization/policyDefinitions/{policyDefinitionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "policyDefinitionName", policyDefinitionName)  # type String
    Swagger.set_param(_ctx.path, "managementGroupId", managementGroupId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyDefinitionsDeleteAtManagementGroup(_api::PolicyDefinitionsApi, policyDefinitionName::String, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsDeleteAtManagementGroup(_api, policyDefinitionName, api_version, managementGroupId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyDefinitionsDeleteAtManagementGroup(_api::PolicyDefinitionsApi, response_stream::Channel, policyDefinitionName::String, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsDeleteAtManagementGroup(_api, policyDefinitionName, api_version, managementGroupId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieves a policy definition in a subscription.
This operation retrieves the policy definition in the given subscription with the given name.
Param: policyDefinitionName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PolicyDefinition
"""
function _swaggerinternal_policyDefinitionsGet(_api::PolicyDefinitionsApi, policyDefinitionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicyDefinition, "/subscriptions/{subscriptionId}/providers/Microsoft.Authorization/policyDefinitions/{policyDefinitionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "policyDefinitionName", policyDefinitionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyDefinitionsGet(_api::PolicyDefinitionsApi, policyDefinitionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsGet(_api, policyDefinitionName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyDefinitionsGet(_api::PolicyDefinitionsApi, response_stream::Channel, policyDefinitionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsGet(_api, policyDefinitionName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieve a policy definition in a management group.
This operation retrieves the policy definition in the given management group with the given name.
Param: policyDefinitionName::String (required)
Param: api_version::String (required)
Param: managementGroupId::String (required)
Return: PolicyDefinition
"""
function _swaggerinternal_policyDefinitionsGetAtManagementGroup(_api::PolicyDefinitionsApi, policyDefinitionName::String, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicyDefinition, "/providers/Microsoft.Management/managementgroups/{managementGroupId}/providers/Microsoft.Authorization/policyDefinitions/{policyDefinitionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "policyDefinitionName", policyDefinitionName)  # type String
    Swagger.set_param(_ctx.path, "managementGroupId", managementGroupId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyDefinitionsGetAtManagementGroup(_api::PolicyDefinitionsApi, policyDefinitionName::String, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsGetAtManagementGroup(_api, policyDefinitionName, api_version, managementGroupId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyDefinitionsGetAtManagementGroup(_api::PolicyDefinitionsApi, response_stream::Channel, policyDefinitionName::String, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsGetAtManagementGroup(_api, policyDefinitionName, api_version, managementGroupId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieves a built-in policy definition.
This operation retrieves the built-in policy definition with the given name.
Param: policyDefinitionName::String (required)
Param: api_version::String (required)
Return: PolicyDefinition
"""
function _swaggerinternal_policyDefinitionsGetBuiltIn(_api::PolicyDefinitionsApi, policyDefinitionName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicyDefinition, "/providers/Microsoft.Authorization/policyDefinitions/{policyDefinitionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "policyDefinitionName", policyDefinitionName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyDefinitionsGetBuiltIn(_api::PolicyDefinitionsApi, policyDefinitionName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsGetBuiltIn(_api, policyDefinitionName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyDefinitionsGetBuiltIn(_api::PolicyDefinitionsApi, response_stream::Channel, policyDefinitionName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsGetBuiltIn(_api, policyDefinitionName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieves policy definitions in a subscription
This operation retrieves a list of all the policy definitions in a given subscription.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PolicyDefinitionListResult
"""
function _swaggerinternal_policyDefinitionsList(_api::PolicyDefinitionsApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicyDefinitionListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Authorization/policyDefinitions", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyDefinitionsList(_api::PolicyDefinitionsApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsList(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyDefinitionsList(_api::PolicyDefinitionsApi, response_stream::Channel, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsList(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieve built-in policy definitions
This operation retrieves a list of all the built-in policy definitions.
Param: api_version::String (required)
Return: PolicyDefinitionListResult
"""
function _swaggerinternal_policyDefinitionsListBuiltIn(_api::PolicyDefinitionsApi, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicyDefinitionListResult, "/providers/Microsoft.Authorization/policyDefinitions", ["azure_auth"])
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyDefinitionsListBuiltIn(_api::PolicyDefinitionsApi, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsListBuiltIn(_api, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyDefinitionsListBuiltIn(_api::PolicyDefinitionsApi, response_stream::Channel, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsListBuiltIn(_api, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieve policy definitions in a management group
This operation retrieves a list of all the policy definitions in a given management group.
Param: api_version::String (required)
Param: managementGroupId::String (required)
Return: PolicyDefinitionListResult
"""
function _swaggerinternal_policyDefinitionsListByManagementGroup(_api::PolicyDefinitionsApi, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicyDefinitionListResult, "/providers/Microsoft.Management/managementgroups/{managementGroupId}/providers/Microsoft.Authorization/policyDefinitions", ["azure_auth"])
    Swagger.set_param(_ctx.path, "managementGroupId", managementGroupId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyDefinitionsListByManagementGroup(_api::PolicyDefinitionsApi, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsListByManagementGroup(_api, api_version, managementGroupId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyDefinitionsListByManagementGroup(_api::PolicyDefinitionsApi, response_stream::Channel, api_version::String, managementGroupId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyDefinitionsListByManagementGroup(_api, api_version, managementGroupId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export policyDefinitionsCreateOrUpdate, policyDefinitionsCreateOrUpdateAtManagementGroup, policyDefinitionsDelete, policyDefinitionsDeleteAtManagementGroup, policyDefinitionsGet, policyDefinitionsGetAtManagementGroup, policyDefinitionsGetBuiltIn, policyDefinitionsList, policyDefinitionsListBuiltIn, policyDefinitionsListByManagementGroup


text:
using MosekTools, PyPlot, LaTeXStrings, MarkovBounds
# use Catalyst to define the reaction network under consideration
# in this case, we investigate Schlögel's system
rn = @reaction_network begin
    0.01*2, 2X --> X
    1.0, ∅ --> X
end
X = species(rn)[1]
# specify initial molecular counts in system
x0 = Dict(X => 1)
x_scale =  Dict(X => 10.0)

# compute bounds on mean molecular counts of the product counts at steady state
# for different truncation orders 2:5
orders = 2:2:6
sol = Dict()
for m in orders
    sol[m] = stationary_mean(rn, [X], x0, m, Mosek.Optimizer, auto_scaling = false, scales = x_scale)
end

fig, ax = subplots()
ax.plot(orders, [sol[m].bounds[X][1] for m in orders], label="lower bound", color="blue")
ax.plot(orders, [sol[m].bounds[X][2] for m in orders], label="upper bound", color="red")
ax.legend()
ax.set_xlabel("truncation order")
ax.set_ylabel(L"⟨ X_{\infty} ⟩")
display(fig)

# same for variance
sol = Dict()
for m in orders
    sol[m] = stationary_variance(rn, [X], x0, m, Mosek.Optimizer, auto_scaling = false, scales = x_scale)
end
fig, ax = subplots()
ax.plot(orders, [sol[m].bounds[X] for m in orders], label="upper bound", color="red")
ax.legend()
ax.set_xlabel("truncation order")
ax.set_ylabel(L"V(X_{\infty})")
display(fig)

# transient moments
nT = 10
traj = []
Tf_range = [0.0, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 5.0, 7.5,
            10.0, 12.5, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0]
order = 4
for Tf in Tf_range
    trange = range(Tf/nT, stop=Tf, length=nT)
    sol = transient_mean(rn, [X], x0, order, trange, Mosek.Optimizer, scales = x_scale)
    push!(traj, sol)
end
fig, ax = subplots()
ax.plot(Tf_range, [sol.bounds[X][1] for sol in traj], color="red", label="lower bound")
ax.plot(Tf_range, [sol.bounds[X][2] for sol in traj], color="blue", label="upper bound")
ax.legend()
ax.set_xlabel("time")
ax.set_ylabel(L"⟨ X ⟩")
display(fig)

# variance
traj = []
for Tf in Tf_range
    trange = range(Tf/nT, stop=Tf, length=nT)
    sol = transient_variance(rn, [X], x0, order, trange, Mosek.Optimizer, scales =  x_scale)
    push!(traj, sol)
end
fig, ax = subplots()
ax.plot(Tf_range, [sol.bounds[X][1] for sol in traj], color="red", label="lower bound")
ax.legend()
ax.set_xlabel("time")
ax.set_ylabel(L"Var(X)")
display(fig)


text:
"""
    ConfMat(labels)
    ConfMat(labels, y, ŷ)

# Example
```julia
julia> labels = [:dog, :cat]
2-element Vector{Symbol}:
 :dog
 :cat

 julia> y = rand(labels, 100);

 julia> ŷ = rand(labels, 100);

 julia> cm = ConfMat(labels, y, ŷ)
 ConfMat{2,Symbol}
 Actual class ╲ Predicted class │ :dog  :cat
 ───────────────────────────────┼───────────
 :dog                           │   30    24
 :cat                           │   21    25

```
"""
struct ConfMat{N, T}
    labels
    counts
end

function ConfMat(labels::Vector, counts::NamedMatrix)
    N = length(labels)
    T = eltype(labels)
    ConfMat{N, T}(labels, counts)
end

function ConfMat(labels::Vector)
    #labels = vec(labels)
    n = length(labels)
    c = zeros(Int, n, n)
    c = NamedArray(c, (labels, labels), ("Actual class", "Predicted class"))
    ConfMat(labels, c)
end

function ConfMat(labels::T, y::T, ŷ::T) where {E, T <: Vector{E}}
    cm = ConfMat(labels)
    update!(cm, y, ŷ)
end

function update!(cm::ConfMat{2,}, y::OneHotMatrix, ŷ::AbstractMatrix;
                 threshold = 0.5)
    y = onecold(y, cm.labels)
    y = labidx(y, cm.labels)
    ŷ = (@. Int(ŷ > threshold) + 1) |> vec
    _update!(cm, y, ŷ)
end

function update!(cm::ConfMat{N,}, y::OneHotMatrix, ŷ::AbstractMatrix) where N
    y = onecold(y, cm.labels)
    ŷ = onecold(ŷ, cm.labels)
    update!(cm, y, ŷ)
end

function update!(cm::ConfMat{N, T}, a::Vector{T}, p::Vector{T}) where {N, T}
    a = labidx(a, cm.labels)
    p = labidx(p, cm.labels)
    _update!(cm, a, p)
end

function _update!(cm::ConfMat, y::Vector{Int}, ŷ::Vector{Int})
    length(ŷ) == length(y) || throw(DimensionMismatch(
        "Number of actual values and number of predicted values mismatch."))
    for i in eachindex(ŷ)
        @inbounds r = y[i]
        @inbounds c = ŷ[i]
        @inbounds cm.counts[r, c] += 1
    end
    cm
end

"""
    labidx(v, l)

"""
function labidx(v, l)
    i = indexin(v, l)
    all(.!isnothing.(i)) || throw(ArgumentError(
        "Unknown labels in supplied data."))
    Vector{Int}(i)
end

summary(cm::ConfMat{N, T}) where {N, T} =
    string("ConfMat{", N, ",", T, "}")

Base.print(cm::ConfMat) = print(cm.counts)

function Base.show(io::IO, cm::ConfMat)
    println(io, summary(cm))
    #show(io, cm.labels)
    tmp = IOBuffer()
    show(tmp, cm.counts)
    #println(io)
    sh = String(take!(tmp))
    for l ∈ split(sh, "\n")[2:end]
        println(io, l)
    end
end


text:
using HDFS

host = isempty(ARGS) ? "localhost" : ARGS[1]
port = isempty(ARGS) ? 9000 : int(ARGS[2])

# connect
fs = hdfs_connect(host, port)
println("connected to hdfs:://$(host):$(port)")

# print capacity, used
capacity = hdfs_get_capacity(fs)
used = hdfs_get_used(fs)
println("$(used) used from total $(capacity) bytes")
@assert 0 < used < capacity

blk_sz = hdfs_get_default_block_size(fs)
println("default block sz: $(blk_sz)")
@assert blk_sz > 0

# get current working directory
pwd1 = pwd(fs)
println("current directory: $(pwd1)")
@assert !isempty(pwd1)

# make directory
const TEST_DIR = "test_dir"
const TEST_FILE = "test_file"
hdfs_exists(fs, TEST_DIR) && hdfs_delete(fs, TEST_DIR)
mkdir(fs, TEST_DIR)
println("created new directory $(TEST_DIR)")
@assert true == isdir(fs, TEST_DIR)
cd(fs, TEST_DIR)
pwd2 = pwd(fs)
println("current directory: $(pwd2)")
@assert pwd2 == string(pwd1, "/", TEST_DIR)

@assert false == hdfs_exists(fs, TEST_FILE)                         # try exists on a file
fi = open(fs, TEST_FILE, "w")                                  # create the file
println("created new file $(TEST_FILE)")
@assert 600 == write(fi, repeat("hello ", 100))                # write to the file
flush(fi)                                                      # flush writes
@assert 0 == close(fi)                                         # close the file
println("wrote 600 bytes to file")
@assert false == isdir(fs, TEST_FILE)

@assert true == hdfs_exists(fs, TEST_FILE)                          # try exists on a file
fi = open(fs, TEST_FILE, "r")                                  # read from the file
avlb = nb_available(fi)
println("bytes available: $(avlb)")
finfo = stat(fs, TEST_FILE)
println("file size bytes: $(finfo.size)")
@assert 600 == finfo.size
@assert 0 <= avlb <= finfo.size

println("reading...")
buff = read(fi, Array(Uint8, 5))
readstr = bytestring(buff[1:5])
@assert "hello" == readstr

@assert 5 == position(fi)
seek(fi, 2)

buff = read(fi, Array(Uint8, 5))
readstr = bytestring(buff[1:5])
@assert "llo h" == readstr

# test read at end of file
seek(fi, finfo.size-10)
(buff,bytes) = hdfs_read(fi, 50)
@assert 10 == bytes

println("closing...")
@assert 0 == close(fi)

blocks = hdfs_blocks(fs, TEST_FILE, 0, finfo.size)
@assert length(blocks) == 1

const NEW_TEST_FILE = string(TEST_FILE, ".new")
const NEW_TEST_FILE_2 = string(TEST_FILE, ".new2")
cp(fs, TEST_FILE, fs, NEW_TEST_FILE)
@assert true == hdfs_exists(fs, NEW_TEST_FILE)
mv(fs, NEW_TEST_FILE, fs, NEW_TEST_FILE_2)
@assert true == hdfs_exists(fs, NEW_TEST_FILE_2)
@assert false == hdfs_exists(fs, NEW_TEST_FILE)

println("listing the directory...")
entries = readdir(fs, ".")
@assert 2 == length(entries)
for fname in entries
    println("\t$(fname)")
end

println("cleaning up...")
cd(fs, "..")
rm(fs, TEST_DIR)

println("passed")


text:
# Tests to simulate multivariate traits

## Get an ultrametric network ('net' already in global scope)
net = readTopology("(((Ag:5,(#H1:1::0.056,((Ak:2,(E:1,#H2:1::0.004):1):1,(M:2)#H2:1::0.996):1):1):1,(((((Az:1,Ag2:1):1,As:2):1)#H1:1::0.944,Ap:4):1,Ar:5):1):1,(P:4,20:4):3,165:7);");

@testset "Simulate data and check means and dimensions" begin

# one single trait
pars = ParamsMultiBM([5.0], ones(1,1))
sim = simulate(net, pars)
@test sim[:Tips, :Exp] == [5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0]

trait_dim = 3

## Simulate a MBD
Random.seed!(17920921); # fix the seed

μ = randn(trait_dim)
Σ = randn(trait_dim, trait_dim)
Σ = Σ * Σ' # needs to be positive definite

@test_throws ErrorException ParamsMultiBM(μ[1:2], Σ)

pars = ParamsMultiBM(μ, Σ); # params of a MBD
@test_logs show(devnull, pars)

sim = simulate(net, pars); # simulate according to a BM
@test_logs show(devnull, sim)
@test_throws ErrorException sim[:Tips, :Broken]

# Extract simulated values
traitsTips = sim[:Tips];
traitsNodes = sim[:InternalNodes];

# Check dimensions
@test size(traitsTips) == (trait_dim, net.numTaxa)
@test size(traitsNodes) == (trait_dim, net.numNodes - net.numTaxa)

# Check means (no shifts)
@test sim[:All, :Exp] ≈ μ * ones(net.numNodes)'

end

###############################################################################
## Test of distibution (fixed root)
###############################################################################
@testset "Simulate test distribution (fixed root)" begin
trait_dim = 3

## Generate some values
Random.seed!(18480224); # fix the seed
μ = randn(trait_dim)
Σ = randn(trait_dim, trait_dim)
Σ = Σ * Σ' # needs to be positive definite
pars = ParamsMultiBM(μ, Σ); # params of a MBD

N = 10000
S = length(tipLabels(net));
μ_sim = zeros(trait_dim, S)
Σ_sim = zeros(trait_dim * S, trait_dim * S)
for i = 1:N
    tips = simulate(net, pars)[:Tips]
    μ_sim .+= tips
    v_sim = vec(tips)
    Σ_sim .+= v_sim * v_sim'
end
μ_sim ./= N
Σ_sim ./= N
Σ_sim = Σ_sim - vec(μ_sim) * vec(μ_sim)'

## Check means
μ_true = μ * ones(S)'
@test isapprox(μ_sim, μ_true, atol=0.2)

## Check covariance
Ψ = Matrix(vcv(net))
Σ_true = kron(Ψ, Σ)
@test isapprox(Σ_sim, Σ_true, atol=3.0) # norm L2 of 36x36 matrix

end

###############################################################################
## Test of distibution (random root)
###############################################################################
@testset "Simulate test distribution (random root)" begin
trait_dim = 3

## Generate some values
Random.seed!(24452384); # fix the seed
μ = randn(trait_dim)
Σ = randn(trait_dim, trait_dim)
Σ = Σ * Σ' # needs to be positive definite
pars = ParamsMultiBM(μ, Σ); # params of a MBD

Σ_root = randn(trait_dim, trait_dim)
Σ_root = Σ_root * Σ_root' / 10.0 # needs to be positive definite (and not too big to reduce variance in test)
pars.varRoot = Σ_root
pars.randomRoot = true

@test_logs show(devnull, pars)
show(devnull, pars)

N = 10000
S = length(tipLabels(net));
μ_sim = zeros(trait_dim, S)
Σ_sim = zeros(trait_dim * S, trait_dim * S)
for i = 1:N
    tips = simulate(net, pars)[:Tips]
    μ_sim .+= tips
    v_sim = vec(tips)
    Σ_sim .+= v_sim * v_sim'
end
μ_sim ./= N
Σ_sim ./= N
Σ_sim = Σ_sim - vec(μ_sim) * vec(μ_sim)'

## Check means
@test isapprox(μ_sim, μ * ones(S)', atol=0.2)

## Check covariance
Ψ = Matrix(vcv(net))
Σ_true = kron(Ψ, Σ) + kron(ones(S, S), pars.varRoot)
@test isapprox(Σ_sim, Σ_true, atol=3.6) # norm L2 of 36x36 matrix

end

################################################################################
## With shifts
################################################################################

net = readTopology("(A:2.5,((B:1,#H1:0.5::0.4):1,(C:1,(D:0.5)#H1:0.5::0.6):1):0.5);")

@testset "Simulate with Shifts" begin
trait_dim = 3

Random.seed!(275698234545); # fix the seed

## Test construction function
@test_throws ErrorException ShiftNet(net.edge[7], [1.0, 2.0],  net) # can't put a shift on hybrid branch
@test_throws ErrorException ShiftNet(net.node[6], [1.0, 2.0],  net) # can't put a shift on hybrid branch
@test_throws ErrorException ShiftNet([net.node[7], net.node[6]], [1.0 2.0], net) # dimensions don't match
@test ShiftNet(net.edge[8], [1.0, 2.0],  net).shift ≈ ShiftNet([net.edge[8]], [1.0 2.0],  net).shift
@test ShiftNet(net.edge[8], [1.0, 2.0],  net).shift ≈ ShiftNet(net.node[7], [1.0, 2.0],  net).shift
@test ShiftNet(net.node[7], [1.0, 2.0],  net).shift ≈ ShiftNet([net.node[7]], [1.0 2.0],  net).shift

## Concatenate function
sh1 = ShiftNet(net.node[7], [1.0, 2.0],  net)*ShiftNet(net.node[9], [3.0, -1.5],  net)
@test_logs show(devnull, sh1)
@test sh1.shift ≈ ShiftNet([net.node[7], net.node[9]], [1.0 2.0; 3.0 -1.5],  net).shift
@test_throws ErrorException sh1*ShiftNet(net.edge[7], [4.0, 3.5],  net) # can't concatenate if the two affect the same edges
@test sh1.shift ≈ (sh1*ShiftNet([net.node[7]], [1.0 2.0],  net)).shift
@test_throws ErrorException sh1*ShiftNet(net.edge[8], [4.0, 3.5, 5.0],  net) # can't concatenate if the two affect the same edges

## Values and edge numbers functions
@test getShiftEdgeNumber(sh1) == [-1, 8]
@test getShiftValue(sh1) == [3.0 -1.5; 1.0 2.0]

## Hybrid shifts
@test shiftHybrid([4.5 2.0], net).shift ≈ ShiftNet(net.edge[6], [4.5, 2.0], net).shift
@test_throws ErrorException shiftHybrid([4.5 2.0; 3.0 5.0], net) # dimension mismatch

## Distributions
μ = randn(trait_dim)
Σ = randn(trait_dim, trait_dim)
Σ = Σ * Σ' # needs to be positive definite

@test ParamsMultiBM(μ, Σ, net).shift.shift ≈ ParamsMultiBM(μ, Σ, ShiftNet(net, trait_dim)).shift.shift
@test_throws ErrorException ParamsMultiBM(μ, Σ, ShiftNet(net, 1))

# Shift at root causes an error.
sh = ShiftNet(net.node[7], [1.0, 2.0, -1.0],  net)*ShiftNet(net.node[9], [3.0, -1.5, 4.2],  net)
pars = ParamsMultiBM(μ, Σ, sh)
@test_throws ErrorException simulate(net, pars)

# One shift, not at the root
sh = ShiftNet(net.node[7], [1.0, 2.0, -1.0],  net)
pars = ParamsMultiBM(μ, Σ, sh)

@test_logs show(devnull, pars)
show(devnull, pars)

N = 10000
S = length(tipLabels(net));
μ_sim = zeros(trait_dim, S)
Σ_sim = zeros(trait_dim * S, trait_dim * S)
for i = 1:N
    tips = simulate(net, pars)[:Tips]
    μ_sim .+= tips
    v_sim = vec(tips)
    Σ_sim .+= v_sim * v_sim'
end
μ_sim ./= N
Σ_sim ./= N
Σ_sim = Σ_sim - vec(μ_sim) * vec(μ_sim)'

## Check means
sim = simulate(net, pars)
μ_true = sim[:Tips, :Exp]
@test μ_true ≈ μ .+ [0 0 1 0.6; 0 0 2 1.2; 0 0 -1 -.6] # from shift on edge 8 only
@test isapprox(μ_sim, μ_true, atol=0.2)

## Check covariance
Ψ = Matrix(vcv(net))
Σ_true = kron(Ψ, Σ)
@test isapprox(Σ_sim, Σ_true, norm = x -> maximum(abs.(x)), atol=0.4)
end


text:
using Latexify
using WatchJuliaBurn

ord_keys = sort(collect(WatchJuliaBurn.func_to_emojis), by=x->string(x[1]))
function to_string(emojis)
    join([emoji isa Tuple ? string(emoji[1]) : string(emoji) for emoji in emojis], ", ")
end

function to_version(emojis)
    if any(x->x isa Tuple, emojis)
        return join([emoji isa Tuple ? string(emoji[2]) : "1" for emoji in emojis], ", ")
    else
        return ""
    end
end

ar = reduce(vcat, ['`' * string(key) * '`' to_string(emojis) to_version(emojis)]
    for (key, emojis) in ord_keys)
ar = vcat(["Function" "Emojis" "Julia Version"], ar)

md_ar = md(ar; latex=false)

dots = "https://raw.githubusercontent.com/JuliaLang/julia/master/doc/src/assets/julia.ico"
intro = """# ⌚<img src="$(dots)" height="26"/>🔥.jl

WatchJuliaBurn aims at destroying the look of your code by adding emojis like :smile: and kaomojis like c╯°□°ↄ╯ instead of your favorites Julia functions.
Don't hesitate to add your worst creations via PR.

## List of emojis

"""

outro = """
## Control Flow
You can now replace boring old try/catch/finally clauses with fancy monkey flow!
```julia
@🐒 begin
    🙈
    	💣()
    🙊(💥)
    	😥(💥)
    🙉
    	🍌()
    end
```
Parsing may behave weird when there are infix operators around the block. Try enclosing everything with parenthesis like `@🐒(begin ... end)` if that happens.

## REPL

You can use the [EmojiSymbols.jl](https://github.com/wookay/EmojiSymbols.jl) package to super-turbo-charge your REPL experience!
"""

open(joinpath(@__DIR__, "..", "README.md"), "w") do io
    write(io, intro * string(md_ar) * outro)
end


text:
using PrettyTables

# MRExperiment is summarized in terms of its loop counters.
@recipe function f(expt::MRExperiment)
    acq = expt.data
    nloop = length(acq[1].loop_counters)
    t = scanner_time(expt)
    t = t .- t[1]
    all_labels = counter_labels(expt)
    labels = []
    title --> "Loop counter summary"
    for i = 1:nloop
        c = [a.loop_counters[i] for a in acq]
        if any(c .!= 0)
            @series begin
                label --> "$(all_labels[i]) (index $i)"
                markershape --> :circle
                markersize --> 1
                t, c
            end
        end
    end
end

filter_zero_counters(counter_vals,col) = any(counter_vals[:,col] .!= 0)

"""
    summarize(expt::MRExperiment)

Print detailed textural summary of the MR experiment. `kwargs` are passed to
`PrettyTables.pretty_table()`, for example, to show loop counters which are all
zero, add the keyword argument `filters_col=()`.
"""
function summarize(io::IO, expt::MRExperiment; kwargs...)
    kwargs = merge((crop=:horizontal, filters_col=(filter_zero_counters,),), kwargs)
    counter_names = collect(fieldnames(LoopCounters))
    counter_vals = zeros(Int, length(expt.data), length(counter_names))
    for (i,acq) in enumerate(expt.data)
        counter_vals[i,:] = acq.loop_counters
    end
    pretty_table(io, counter_vals, counter_names; kwargs...)
end

function summarize(expt::MRExperiment; kwargs...)
    summarize(stdout, expt; kwargs...)
end

# Color map similar to the one which comes with Felix NMR
const RGBN0f8 = RGB{Colors.N0f8}

"""
A popular color map for multidimensional NMR data as defined in Felix NMR
"""
const felix_colors = [
    RGBN0f8(0.0,0.0,1.0)  ,
    RGBN0f8(0.063,0.0,1.0),
    RGBN0f8(0.247,0.0,1.0),
    RGBN0f8(0.749,0.0,1.0),
    RGBN0f8(1.0,0.0,0.749),
    RGBN0f8(1.0,0.0,0.498),
    RGBN0f8(1.0,0.0,0.247),
    RGBN0f8(1.0,0.0,0.0)  ,
    RGBN0f8(1.0,0.247,0.0),
    RGBN0f8(1.0,0.498,0.0),
    RGBN0f8(1.0,0.561,0.0),
    RGBN0f8(1.0,0.627,0.0),
    RGBN0f8(1.0,0.749,0.0),
    RGBN0f8(1.0,0.812,0.0),
    RGBN0f8(1.0,0.871,0.0),
    RGBN0f8(1.0,1.0,0.063),
    RGBN0f8(1.0,1.0,0.125),
    RGBN0f8(1.0,1.0,0.188),
    RGBN0f8(1.0,1.0,0.247),
    RGBN0f8(1.0,1.0,0.373),
    RGBN0f8(1.0,1.0,0.498),
    RGBN0f8(1.0,1.0,0.561),
    RGBN0f8(1.0,1.0,0.624),
    RGBN0f8(1.0,1.0,0.686),
    RGBN0f8(1.0,1.0,0.749),
    RGBN0f8(1.0,1.0,0.812),
    RGBN0f8(1.0,1.0,0.914),
    RGBN0f8(1.0,1.0,0.882),
]

# FIXME: Type piracy; should get this into AxisArrays instead!
using RecipesBase

@recipe function plot(a::AxisArray)
    ax1 = AxisArrays.axes(a,1)
    xlabel --> AxisArrays.axisname(ax1)
    if ndims(a) == 1
        ax1.val, a.data
    else
        ax2 = AxisArrays.axes(a,2)
        # Categorical axes print as a set of labelled series
        if axistrait(ax2) === Categorical
            for i in eachindex(ax2.val)
                @series begin
                    label --> "$(AxisArrays.axisname(ax2)) $(ax2.val[i])"
                    ax1.val, a.data[:,i]
                end
            end
        else
            # Other axes as a 2D array
            ylabel --> AxisArrays.axisname(ax2)
            seriestype --> :heatmap
            ax1.val, ax2.val, a.data
        end
    end
end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
#check balance level by level using a post-order treewalk (from deepest level up to root)

#flip jth digit of num and return
flipbit(num, j) = num ⊻ ( 1 << (j-1) )

const PERIODIC = true

function set_max_depth_AMR!(tree::Node{N,T,D}, max_depth::Int) where {N,T,D}
	delete_high_level_nodes!(tree, max_depth, tree.length[1])
end

function delete_high_level_nodes!(node::Node{N,T,D}, max_depth::Int, root_node_length::T) where {N,T,D}
	depth = log2(round(root_node_length / node.length[1]))
	if !isLeaf(node)
        #println("This is a node... ")
		if depth < max_depth
	    	@inbounds for i in 1:2^N
	        	#println("open this node")
				delete_high_level_nodes!(node.child[i], max_depth, root_node_length)
	    	end
		else
			#delete child nodes
			node.child = nothing
		end
    end
end

function get_min_node_length!(length::T, node::Node{N,T,D}) where{N,T,D}
	if isLeaf(node)
		return min(length, node.length[1])
	else
		@inbounds for i in 1:2^N
			length = get_min_node_length!(length, node.child[i])
		end
		return length
	end
end

function get_max_tree_depth(tree::Node{N,T,D}) where{N,T,D}
	len_min = 1000.
	len_min = get_min_node_length!(len_min, tree)
	return Int( log2( round( tree.length[1] / len_min) ) )
end

function balance_all_level!(tree::Node{N,T,D}) where{N,T,D}
	max_depth = get_max_tree_depth(tree)
	for i in max_depth:-1:1
		#println("balancing the tree at level ", i)
		balance!(tree, i, tree.length[1])
	end
end

function balance!(node::Node{N,T,D}, max_depth::Int, root_node_length::T) where{N,T,D}
	depth = Int(log2(round(root_node_length / node.length[1])))

    if isLeaf(node)
		#println("in a leaf node at level ", depth)
		if depth == max_depth #only do it at the deepest level
			if node.ID == 1 || node.ID == 2^N #we only need to check balance for two diagonal nodes
				#println("in a leaf node at level ", depth, "... check balance!")
				for i in 1:N
					#check both directions (0 & 1) along each dimension (i), which is N*2 directions in total
					ngb_node_0 = check_ngb_node_and_refine_if_necessary(node, i, 0, node.length[1])
					ngb_node_1 = check_ngb_node_and_refine_if_necessary(node, i, 1, node.length[1])
					@assert (ngb_node_0.length[1] / node.length[1]) < 2.5
					@assert (ngb_node_1.length[1] / node.length[1]) < 2.5
				end
			end
		end
    else
        #println("This is a node at level ", depth)
		#if depth >= max_depth and it's still a node,
		#then it's already checked during the previous iterations
		if depth < max_depth
	    	@inbounds for i in 1:2^N
	        	#println("open this node")
				balance!(node.child[i], max_depth, root_node_length)
	    	end
		end
    end
end

#search both directions for each dimension (dim)
function check_ngb_node_and_refine_if_necessary(node, dim, direction, length0)
    #if dim == 2 #find N & S ngbs
    #if direction == 1 #we're looking for the north ngb
    if node.parent == nothing
        #node is root! we're at the boundary! #TODO: peridoic B.C.?
        return nothing
    end

    #node is a south child, so it can give us a north ngb directly from its sibling (we need to distinguish between SW & SE)
    if getbit(node.ID-1, dim) != direction #meaning that node has a sibling along direction
        idx_ngb = flipbit(node.ID-1, dim) #flipbit = opposite direction in dim-th dimension
		#@show node.ID-1, dim, idx_ngb
        return node.parent.child[idx_ngb+1]
    end

    #Okay, node is a north child, so it can't give us a north ngb (it has no north siblings).
    #Therefore, we have to move up to the parent node. If that parent node is still not a south child,
    #then we keep moving up (to the parent's parent) until we find a parent who is a south child or until
    #we reach the root which means that we're at the boundary and there is no ngb along this direction
    uncle = check_ngb_node_and_refine_if_necessary(node.parent, dim, direction, length0) #an uncle is the sibling of a parent (duh)
    if uncle == nothing #root!
		if PERIODIC
			#@show "root!!!"
			uncle = node.parent #set uncle as node's parent and proceed (don't return!!!)
		else
			return uncle
		end
    end

    #if uncle is a leaf, then it's a least 2x bigger than node, so it doesn't matter if node is NW or NE
    #we check if the uncle needs to be refined
    if isLeaf(uncle)
        if uncle.length[1] / length0 > 2.5 #should be just 2 but use 2.5 to avoid round-off error
			#println("uncle.length = ", uncle.length[1], "  length0 = ", length0, "  refine the node!!!")
            split_node!(uncle)
        else
            return uncle #only differ by a factor of 2 (one level)
        end
    end

    #uncle is not a leaf, so we return the ngb cousin
    idx_ngb = flipbit(node.ID-1, dim) #flipbit = opposite direction in dim-th dimension
    return uncle.child[idx_ngb+1] #SW (bit flip at the dim-th index)
end

#for j in 1:2
#    for i in 1:2^2
#        @show getbit(i-1,j), bitstring((i-1))[end-8:end], bitstring(flipbit(i-1,j))[end-8:end]
#    end
#end

#NE = [1,1]
#NW = [0,1]
#SE = [1,0]
#SW = [0,0]

#dim = 1: x-direction , or E & W
#dim = 2: y-direction , or N & S

#direction = 0 or 1 (control E or W if dim == 1; N or S if dim ==2)

#if node is a north child, then its southern ngb is its sibling, so the balance is fine, so we don't have to check its southern ngb
#we only need to check its northern ngb

#=

#search both directions for each dimension (dim)
check_ngb_node_and_refine_if_necessary(node, dim, direction, length0)
    if dim == 2 #find N & S ngbs
        if direction == 1 #we're looking for the north ngb
            if node.parent == nothing
                #node is root! we're at the boundary! #TODO: peridoic B.C.?
                return nothing
            end

            #node is a south child, so it can give us a north ngb directly from its sibling (we need to distinguish between SW & SE)
            if node.ID[dim] != direction
            if node.ID == [0,0] #SW
                return node.parent.child[0,1] #NW (bit flip at the dim-th index)
            elseif node.ID == [1,0] #SE
                return node.parent.child[1,1] #NE (bit flip at the dim-th index)
            end
            end

            #Okay, node is a north child, so it can't give us a north ngb (it has no north siblings).
            #Therefore, we have to move up to the parent node. If that parent node is still not a south child,
            #then we keep moving up (to the parent's parent) until we find a parent who is a south child or until
            #we reach the root which means that we're at the boundary and there is no ngb along this direction
            uncle = check_ngb_node_and_refine_if_necessary(node.parent, dim, direction, length0) #an uncle is the sibling of a parent (duh)
            if uncle == nothing #root!
                return uncle
            end

            #if uncle is a leaf, then it's a least 2x bigger than node, so it doesn't matter if node is NW or NE
            #we check if the uncle needs to be refined
            if isLeaf(uncle)
                if uncle.length / length0 > 2.5 #should be just 2 but use 2.5 to avoid round-off error
                    split_node!(uncle)
                else
                    return uncle #only differ by a factor of 2 (one level)
                end
            end

            #uncle is not a leaf, so we return the ngb cousin
            if node.ID == [0,1] #NW
                return uncle.child[0,0] #SW (bit flip at the dim-th index)
            elseif node.ID == [1,1] #NE
                return uncle.child[1,0] #SE (bit flip at the dim-th index)
            end
        else  #direction == 0 #we're looking for the south ngb
            code...
        end #direction
    elseif dim == 1 #find E & W ngbs
        code...
    end #dim
end

=#


text:
function Base.show(io::IO,ct::Counter)
	print(io,"\n    ↳ In nodes           =  $(ct.node)")
	print(io,"\n    ↳ In virtual path    =  $(ct.virtual_path) ")
	print(io,"\n    ↳ In DDF             =  $(ct.ddf)   ")
	print(io,"\n    ↳ Number of branches =  $(ct.branch) " )
	print(io,"\n    ↳ Number of pruning  =  $(ct.prune) ")
end

function Base.show(io::IO, c::ConformationOutput)
	print(io,"\n ✔ Solver = $(c.solver)")
	print(io,"\n ✔ Number of solutions = $(c.number)\n")
	print(io,  " ✔ Solutions \n")
	n = length(c.molecules[1].atoms)
	for k=1:c.number
		print(io,"   ↳ Molecule $(k) with LDE = $(c.molecules[k].lde) \n")
		for i=1:5
			print(io,"    ( $(c.molecules[k].atoms[i].x) , $(c.molecules[k].atoms[i].y) ,  $(c.molecules[k].atoms[i].z) ) \n")
		end
			print(io,"    ⋮  \n")
		for i=n-2:n
			print(io,"    ( $(c.molecules[k].atoms[i].x) , $(c.molecules[k].atoms[i].y) ,  $(c.molecules[k].atoms[i].z) ) \n")
		end

	end
	print(io," ✔ Number of main operations [+-,*,/,√] ")
	print(io,"   $(c.nop) ")
end


text:
doc = """NTT parameter finder.

Usage:
  ntt-params.jl --modulus <q> [ --base <g> ]
  ntt-params.jl --length <n> [ --modulus <q> ] [ --base <g> ]
  ntt-params.jl -h | --help

Options:
  -h --help     Show this screen.
  -n --length   Length of input vector.
  -q --modulus  Modulo arithmetic.
  -g --base     Transform exponential base.
"""

using DocOpt, Primes

args = docopt(doc, version=v"2.0.0")
args_parse(T, key) = args[key] === nothing ? nothing : parse(T, args[key])
(g, q, n) = args_parse(BigInt, "<g>"), args_parse(BigInt, "<q>"), args_parse(Int64, "<n>")

include("find-ntt-methods.jl")
findntt = FindNTT(g, q, n)
findntt()


text:
# Sample postprocessing script for plotting the components of a vector from one snapshot.
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, PyPlot, LaTeXStrings, Printf

file = "/wrk/group/spacephysics/vlasiator/2D/ABA/bulk/bulk.0001000.vlsv"

axisunit = EARTH
colorscale = Linear
cmap = matplotlib.cm.RdBu

var = "B_vol"

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

meta = load(file)

pArgs = Vlasiator.set_args(meta, var, axisunit; normal=:none)

x, y = Vlasiator.get_axis(pArgs)
# [nT]
B = meta["B_vol"]
bx = @views B[1,:].*1e9
bx = reshape(bx, meta.ncells[1], meta.ncells[2])
by = @views B[2,:].*1e9
by = reshape(by, meta.ncells[1], meta.ncells[2])
bz = @views B[3,:].*1e9
bz = reshape(bz, meta.ncells[1], meta.ncells[2])

V = (bx, by, bz)
V_str = ("Bx", "By", "Bz")

# Symmetric range for diverging colormap
vmin = minimum(minimum.(V))
vmax = -vmin

cnorm1, cticks1 = Vlasiator.set_colorbar(colorscale, vmin, vmax)

fig, axs = subplots(1,3,
   figsize=(10, 4), sharex=true, sharey=true, constrained_layout=true)

c1 = axs[1].pcolormesh(x, y, V[1]'; norm=cnorm1, cmap)

for i in eachindex(axs)[2:end]
   axs[i].pcolormesh(x, y, V[i]'; norm=cnorm1, cmap)
end

for (i, ax) in enumerate(axs)
   ax.set_aspect("equal")
   ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
   ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
   ax.set_title(V_str[i])
   ax.set_xlabel(L"x [$R_E$]")
end

axs[1].set_ylabel(L"y [$R_E$]")

# One colorbar for all subplots
cb1 = colorbar(c1; ax=axs, ticks=cticks1, fraction=0.046, pad=0.04)
cb1.ax.set_ylabel("[nT]")
cb1.outline.set_linewidth(1.0)

str_title = @sprintf "B at t = %4.1fs" meta.time

fig.suptitle(str_title, fontsize="xx-large")

savefig("Bvector.png", bbox_inches="tight")


text:
export ZeroPad

"""
`ZeroPad([domainType::Type,] dim_in::Tuple, zp::Tuple)`

`ZeroPad(x::AbstractArray, zp::Tuple)`

Create a `LinearOperator` which, when multiplied to an array `x` of size `dim_in`, returns an expanded array `y` of size `dim_in .+ zp` where `y[1:dim_in[1], 1:dim_in[2] ... ] = x` and zero elsewhere.

```julia
julia> Z = ZeroPad((2,2),(0,2))
[I;0]  ℝ^(2, 2) -> ℝ^(2, 4)

julia> Z*ones(2,2)
2×4 Array{Float64,2}:
 1.0  1.0  0.0  0.0
 1.0  1.0  0.0  0.0

```

"""
struct ZeroPad{T,N} <: LinearOperator
	dim_in::NTuple{N,Int}
	zp::NTuple{N,Int}
end

# Constructors
#standard constructor
function ZeroPad(domainType::Type, dim_in::NTuple{N,Int}, zp::NTuple{M,Int}) where {N,M}
	M != N && error("dim_in and zp must have the same length")
	any([zp...].<0) && error("zero padding cannot be negative")
	ZeroPad{domainType,N}(dim_in,zp)
end

ZeroPad(dim_in::Tuple, zp::NTuple{N,Int})                   where {N} = ZeroPad(Float64, dim_in, zp)
ZeroPad(domainType::Type, dim_in::Tuple, zp::Vararg{Int,N}) where {N} = ZeroPad(domainType, dim_in, zp)
ZeroPad(dim_in::Tuple, zp::Vararg{Int,N})                   where {N} = ZeroPad(Float64, dim_in, zp)
ZeroPad(x::AbstractArray, zp::NTuple{N,Int})                where {N} = ZeroPad(eltype(x), size(x), zp)
ZeroPad(x::AbstractArray, zp::Vararg{Int,N})                where {N} = ZeroPad(eltype(x), size(x), zp)

# Mappings
@generated function mul!(y::AbstractArray{T,N},L::ZeroPad{T,N},b::AbstractArray{T,N}) where {T,N}

	# builds
	#for i1 =1:size(y,1), i2 =1:size(y,2)
	#	y[i1,i2] = i1 <= size(b,1) && i2 <= size(b,2)  ?  b[i1,i2] : 0.
	#end

	ex = "for "
	for i = 1:N ex *= "i$i =1:size(y,$i)," end

	ex = ex[1:end-1] #remove ,

	ex *= " y["
	for i = 1:N ex *= "i$i," end

	ex = ex[1:end-1] #remove ,

	ex *= "] = "
	for i = 1:N ex *= " i$i <= size(b,$i) &&" end

	ex = ex[1:end-2] #remove &&

	ex *= " ?  b["
	for i = 1:N ex *= "i$i," end

	ex = ex[1:end-1] #remove ,
	ex *= "] : 0. end"

	ex = Meta.parse(ex)

end

@generated function mul!(y::AbstractArray{T,N},L::AdjointOperator{ZeroPad{T,N}},b::AbstractArray{T,N}) where {T,N}

	#builds
	#for l = 1:size(y,1), m = 1:size(y,2)
	#	y[l,m] = b[l,m]
	#end

	ex = "for "
	for i = 1:N ex *= "i$i =1:size(y,$i)," end

	ex *= " y["
	idx = ""
	for i = 1:N idx *= "i$i," end

	idx = idx[1:end-1] #remove ,

	ex *= idx

	ex *= "] = b["
	ex *= idx
	ex *= "] end"

	ex = Meta.parse(ex)

end

# Properties

  domainType(L::ZeroPad{T}) where {T} = T
codomainType(L::ZeroPad{T}) where {T} = T

size(L::ZeroPad) = L.dim_in .+ L.zp, L.dim_in

fun_name(L::ZeroPad)       = "[I;0]"
is_AcA_diagonal(L::ZeroPad) = true
diag_AcA(L::ZeroPad) = 1

is_full_column_rank(L::ZeroPad) = true


text:
# ---- Include sub-files
include("objective.jl")
include("est.jl")
include("add_ens.jl")
include("iterate.jl")
include("update.jl")

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
using MaximumLikelihood
using Base.Test
using Distributions

@testset "multivariate normal likelihood" begin
    d = MultivariateNormal([1,2], 3)
    ℓ(θ) = logpdf(d, θ)
    ML = MaximumLikelihood.estimate_ML( ℓ, zeros(2))

    @test ML.θ ≈ [1.0, 2.0]
    @test full(ML.Σ) ≈ diagm([9.0, 9.0])
end


text:
@testset "CoreModel utilities" begin
    cp = test_LP()
    @test n_reactions(cp) == 3
    @test n_metabolites(cp) == 4
    @test n_coupling_constraints(cp) == 0

    cp2 = test_LP()
    @test isequal(cp, cp2)
    cp2.S[1] = 1
    @test !isequal(cp, cp2)
    @test isequal(cp, copy(cp))

    cp = test_coupledLP()
    @test n_coupling_constraints(cp) == 2000
    @test isequal(cp, copy(cp))
end


text:
using GPMTD
using Test

# @test hello("Julia") == "Hello, Julia"
# @test domath(2.0) ≈ 7.0

## test Distributions.Gamma
mean(Gamma(2.0, 0.5)) ≈ 1.0 # the second parameter is the shape

## Assure that Distances doesn't change
x = [0.0, 1.0, 1.0]
a = [[0.0, 1.0, 1.0] [1.0, 0.0, 0.0] [1.0, 0.0, 0.0]]
@test pairDistMat(x) ≈ a
@test pairDistMat( hcat(x,x) ) ≈ a .* sqrt(2.0)

## Covariance
x[3] = -1.0

mp = MaternParams(Inf, 1.0)
D = pairDistMat(x)

corr(0.0, mp)
corr(1.0, mp)

corrMat(D, mp)

v = 5.0
Cov = covMat(D, v, mp)

## Log-likelihood calculation
n = 50
σ2 = 2.0
y = 5.0*randn(n)
x = 1.0*randn(n)
D = pairDistMat(x)
v = 5.0
Cov = covMat(D, v, mp, tol=1.0e-10)

x2 = randn(Int(floor(n/2)))
Dxy = pairDistMat(x, x2)
CovXY = covMat(Dxy, v, mp, tol=1.0e-10)

# @benchmark llik_marg($y, $σ2, $Cov) ## winner
# @benchmark logpdf( MvNormal( $Cov + $σ2*I ), $y)

rpost_MvNNV(y, σ2, Cov)

x_old = reshape(randn(3), 3, 1)
x_new = 0.0
D = pairDistMat(x_old)
expanD(D, x_new, x_old)


text:
export produce_or_load, tagsave, @tagsave, safesave

"""
    produce_or_load([prefix="",] c, f; kwargs...) -> file
Let `s = savename(prefix, c, suffix)`.
If a file named `s` exists then load it and return it.

If the file does not exist then call `file = f(c)`, save `file` as
`s` and then return the `file`.
The function `f` must return a dictionary.
The macros [`@dict`](@ref) and [`@strdict`](@ref) can help with that.

## Keywords
* `tag = true` : Add the Git commit of the project in the saved file.
* `gitpath = projectdir()` : Path to search for a Git repo.
* `suffix = "bson"` : Used in `savename`.
* `force = false` : If `true` then don't check if file `s` exists and produce
  it and save it anyway.
* `kwargs...` : All other keywords are propagated to `savename`.

See also [`savename`](@ref) and [`tag!`](@ref).
"""
produce_or_load(c, f; kwargs...) = produce_or_load("", c, f; kwargs...)
function produce_or_load(prefix::String, c, f;
    tag::Bool = true, gitpath = projectdir(),
    suffix = "bson", force = false, kwargs...)

    s = savename(prefix, c, suffix; kwargs...)

    if !force && isfile(s)
        file = wload(s)
        return file
    else
        if force
            @info "Producing file $s now..."
        else
            @info "File $s does not exist. Producing it now..."
        end
        file = f(c)
        try
            mkpath(dirname(s))
            if tag
                tagsave(s, file, false, gitpath)
            else
            wsave(s, copy(file))
        end
        catch er
            @warn "Could not save file, got error $er. "*
            "\nReturning the file nontheless."
        end
        return file
    end
end

"""
    tagsave(file::String, d::Dict [, safe = false, gitpath = projectdir()])
First [`tag!`](@ref) dictionary `d` and then save `d` in `file`.
If `safe = true` save the file using [`safesave`](@ref).
"""
tagsave(file, d, p::String) = tagsave(file, d, false, p)
function tagsave(file, d, safe = false, gitpath = projectdir(), s = nothing)
    d2 = tag!(d, gitpath, s)
    if safe
        safesave(file, copy(d2))
    else
        wsave(file, copy(d2))
    end
    return d2
end

"""
    @tagsave(file::String, d::Dict [, safe = false, gitpath = projectdir()])
Same as [`tagsave`](@ref) but also add a field `script` that records
the local path of the script that called `@tagsave`, see [`@tag!`](@ref).
"""
macro tagsave(file, d, safe::Bool = false, gitpath = projectdir())
    s = QuoteNode(__source__)
    :(tagsave($(esc(file)), $(esc(d)), $(esc(safe)), $(esc(gitpath)), $s))
end

################################################################################
#                          Backup files before saving                          #
################################################################################

# Implementation inspired by behavior of GROMACS
"""
    safesave(filename, data)

Safely save `data` in `filename` by ensuring that no existing files
are overwritten. Do this by renaming already existing data with a backup-number
ending like `#1, #2, ...`. For example if `filename = test.bson`, the first
time you `safesave` it, the file is saved normally. The second time
the existing save is renamed to `test_#1.bson` and a new file `test.bson`
is then saved.

If a backup file already exists then its backup-number is incremented
(e.g. going from `#2` to `#3`). For example safesaving `test.bson` a third time
will rename the old `test_#1.bson` to `test_#2.bson`, rename the old
`test.bson` to `test_#1.bson` and then save a new `test.bson` with the latest
`data`.

See also [`tagsave`](@ref).
"""
function safesave(f, data)
    recursively_clear_path(f)
    wsave(f, data)
end

#take a path of a results file and increment its prefix backup number
function increment_backup_num(filepath)
    path, filename = splitdir(filepath)
    fname, suffix = splitext(filename)
    m = match(r"^(.*)_#([0-9]+)$", fname)
    if m == nothing
        return joinpath(path, "$(fname)_#1$(suffix)")
    end
    newnum = string(parse(Int, m.captures[2]) +1)
    return joinpath(path, "$(m.captures[1])_#$newnum$(suffix)")
end

#recursively move files to increased backup number
function recursively_clear_path(cur_path)
    isfile(cur_path) || return
    new_path=increment_backup_num(cur_path)
    if isfile(new_path)
        recursively_clear_path(new_path)
    end
    mv(cur_path, new_path)
end