text:
using Pkg # hideall
Pkg.activate("_literate/ISL-lab-8/Project.toml")
Pkg.update()
macro OUTPUT()
    return isdefined(Main, :Franklin) ? Franklin.OUT_PATH[] : "/tmp/"
end;

# ## Getting started

using MLJ
import RDatasets: dataset
using PrettyPrinting
import DataFrames: DataFrame, select, Not
MLJ.color_off() # hide
DTC = @load DecisionTreeClassifier pkg=DecisionTree

carseats = dataset("ISLR", "Carseats")

first(carseats, 3) |> pretty

# We encode a new variable `High` based on whether the sales are higher or lower than 8 and add that column to the dataframe:

High = ifelse.(carseats.Sales .<= 8, "No", "Yes") |> categorical;
carseats[!, :High] = High;

# Let's now train a basic decision tree classifier for `High` given the other features after one-hot-encoding the categorical features:

X = select(carseats, Not([:Sales, :High]))
y = carseats.High;

# ### Decision Tree Classifier

HotTreeClf = OneHotEncoder() |> DTC()

mdl = HotTreeClf
mach = machine(mdl, X, y)
fit!(mach);

# Note `|>` is syntactic sugar for creating a `Pipeline` model from component model instances or model types.
# Note also that the machine `mach` is trained on the whole data.

ypred = predict_mode(mach, X)
misclassification_rate(ypred, y)

# That's right... it gets it perfectly; this tends to be classic behaviour for a DTC to overfit the data it's trained on.
# Let's see if it generalises:

train, test = partition(eachindex(y), 0.5, shuffle=true, rng=333)
fit!(mach, rows=train)
ypred = predict_mode(mach, rows=test)
misclassification_rate(ypred, y[test])

# Not really...
#
# ### Tuning a DTC
#
# Let's try to do a bit of tuning

r_mpi = range(mdl, :(decision_tree_classifier.max_depth), lower=1, upper=10)
r_msl = range(mdl, :(decision_tree_classifier.min_samples_leaf), lower=1, upper=50)

tm = TunedModel(
    model=mdl,
    ranges=[r_mpi, r_msl],
    tuning=Grid(resolution=8),
    resampling=CV(nfolds=5, rng=112),
    operation=predict_mode,
    measure=misclassification_rate
)
mtm = machine(tm, X, y)
fit!(mtm, rows=train)

ypred = predict_mode(mtm, rows=test)
misclassification_rate(ypred, y[test])

# We can inspect the parameters of the best model

fitted_params(mtm).best_model.decision_tree_classifier

# ### Decision Tree Regressor
#

DTR = @load DecisionTreeRegressor pkg=DecisionTree

boston = dataset("MASS", "Boston")

y, X = unpack(boston, ==(:MedV))

train, test = partition(eachindex(y), 0.5, shuffle=true, rng=551);

scitype(X)

# Let's recode the Count as Continuous and then fit a DTR

X = coerce(X, autotype(X, rules=(:discrete_to_continuous,)))

dtr_model = DTR()
dtr = machine(dtr_model, X, y)

fit!(dtr, rows=train)

ypred = MLJ.predict(dtr, rows=test)
round(rms(ypred, y[test]), sigdigits=3)

# Again we can try tuning this a bit, since it's the same idea as before, let's just try to adjust the depth of the tree:

r_depth = range(dtr_model, :max_depth, lower=2, upper=20)

tm = TunedModel(
    model=dtr_model,
    ranges=[r_depth],
    tuning=Grid(resolution=10),
    resampling=CV(nfolds=5, rng=1254),
    measure=rms
)

mtm = machine(tm, X, y)

fit!(mtm, rows=train)

ypred = MLJ.predict(mtm, rows=test)
round(rms(ypred, y[test]), sigdigits=3)

# ## Random Forest
#
# **Note**: the package [`DecisionTree.jl`](https://github.com/bensadeghi/DecisionTree.jl) also has a RandomForest model but it is not yet interfaced with in MLJ.

RFR = @load RandomForestRegressor pkg=ScikitLearn

rf_mdl = RFR()
rf = machine(rf_mdl, X, y)
fit!(rf, rows=train)

ypred = MLJ.predict(rf, rows=test)
round(rms(ypred, y[test]), sigdigits=3)

# ## Gradient Boosting Machine

XGBR = @load XGBoostRegressor

xgb_mdl = XGBR(num_round=10, max_depth=10)
xgb = machine(xgb_mdl, X, y)
fit!(xgb, rows=train)

ypred = MLJ.predict(xgb, rows=test)
round(rms(ypred, y[test]), sigdigits=3)

# Again we could do some tuning for this.

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
export SI

"Constants written in the same system of physical units."
struct UnitSystem
    name::Symbol
    "Boltzmann constant."
    k::Float
    "Atomic mass unit."
    AMU::Float

    UnitSystem(; name::Symbol, k::Number, AMU::Number) = new(name, k, AMU)
end

"Constants written in SI units"
const SI = UnitSystem(
    name = :SI,
    k = 1.38064852e-23, # [m^2 kg/(s^2 K)]
    AMU = 1.66054e-27,  # [kg]
)

ext: jl

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
using BNPRegressionGGA2021
using Random
using Test

const DGSBP = BNPRegressionGGA2021

@testset "Skeleton - constructor" begin
    N1, N0, D0 = 2, 10, 2
    Random.seed!(1)
    y0 = randn(N0)
    y1 = randn(N1)
    X0 = randn(N0, D0)
    X1 = zeros(N1, D0)
    skl = DGSBP.Skeleton(; y0, X0, y1, X1)
end

@testset "Skeleton - accessors" begin
    N1, N0, D0 = 2, 10, 2
    Random.seed!(1)
    y0 = randn(N0)
    y1 = randn(N1)
    X0 = randn(N0, D0)
    X1 = zeros(N1, D0)
    skl = DGSBP.Skeleton(; y0, X0, y1, X1)
    @test DGSBP.ygrid(skl) === y1
    @test DGSBP.Xgrid(skl) === X1
    @test DGSBP.Ngrid(skl) == N1
    @test DGSBP.y(skl) === y0
    @test DGSBP.X(skl) === X0
    @test DGSBP.N(skl) == N0
    @test DGSBP.d(skl) == ones(N0)
    @test DGSBP.r(skl) == ones(N0)
    @test DGSBP.f(skl) == zeros(N1)
    @test DGSBP.rmax(skl) == 1
end

@testset "DGSBPNormal - constructor" begin
    N1, N0, D0 = 2, 10, 2
    Random.seed!(1)
    y0 = randn(N0)
    y1 = randn(N1)
    X0 = randn(N0, D0)
    X1 = zeros(N1, D0)
    skl = DGSBP.Skeleton(; y0, X0, y1, X1)
    m = DGSBP.DGSBPNormal(; skl)
end

@testset "DGSBPNormal - methods" begin
    N1, N0, D0 = 2, 10, 2
    Random.seed!(1)
    y0 = randn(N0)
    y1 = randn(N1)
    X0 = randn(N0, D0)
    X1 = zeros(N1, D0)
    skl = DGSBP.Skeleton(; y0, X0, y1, X1)
    m = DGSBP.DGSBPNormal(; skl)
    @test DGSBP.skeleton(m) === skl
    @test DGSBP.kernel_pdf(m, 0.0, 1) ≈ 1 / √(2 * π)
    @test_throws MethodError DGSBP.kernel_pdf(m, 0, 1)
end

# @testset "Skeleton.jl" begin
#     N1, N0, D0 = 2, 10, 2
#     Random.seed!(1)
#     y0 = randn(N0)
#     X0 = randn(N0, D0)
#     y1 = LinRange(-3, 3, N1) |> collect
#     X1 = zeros(N1, D0)

#     # sampler = BNPRegressionGGA2021.Sampler(; y0, X0, y1, X1)
#     # BNPRegressionGGA2021.m(sampler)
#     # BNPRegressionGGA2021.update_suffstats!(sampler)
#     # BNPRegressionGGA2021.update_χ!(rng, sampler)
#     # BNPRegressionGGA2021.update_d!(rng, sampler)
#     # BNPRegressionGGA2021.update_r!(rng, sampler)
#     # BNPRegressionGGA2021.update_ϕ!(sampler)
#     # BNPRegressionGGA2021.update_β!(rng, sampler)
#     # BNPRegressionGGA2021.update_f!(sampler)
#     # sampler = BNPRegressionGGA2021.Sampler(; y0, X0, y1, X1)
#     # BNPRegressionGGA2021.step!(rng, sampler)
#     # chain = BNPRegressionGGA2021.sample(rng, sampler)
# end


text:
module TestSSQueue

using ConcurrentCollections
using ConcurrentCollections.Implementations: NodeIterator, SSQNode, check_invariance, ==′
using ProgressLogging: @logprogress, @withprogress
using Test
using ..TestDLCRQ: check_concurrent_push_pop!
using ..Utils: ⊏

function test_push_pop_once_int()
    q = DualLinkedQueue{Int}()
    push!(q, 111)
    @test popfirst!(q) == 111
end

function test_push_pop_once_any()
    q = DualLinkedQueue()
    push!(q, 111)
    @test popfirst!(q) == 111
end

function test_error()
    q = DualLinkedQueue{Int}()
    t = @task popfirst!(q)
    yield(t)
    msg = "Interrupting DLQ @$(time_ns())"
    schedule(t, ErrorException(msg); error = true)
    err = try
        Some(wait(t))
    catch e
        e
    end
    @test err isa TaskFailedException
    @test occursin(msg, sprint(showerror, err))
    push!(q, 111)  # this dequeues but ignores the interrupted waiter
    @test popfirst!(q) == 111
end

function test_concurrent_push_pop(ntrials = 100)
    @withprogress name = "concurrent push-pop" begin
        @testset for trial in 1:ntrials
            @logprogress (trial - 1) / ntrials
            q = DualLinkedQueue{Int}()
            check_concurrent_push_pop!(q; nitems = 2^13)
        end
    end
end

function test_iter()
    q = DualLinkedQueue{Int}()
    @test eltype(q) === Int
    push!(q, 111)
    push!(q, 222)
    @test collect(q) ==′ [111, 222]
end

function test_nodeiterator()
    q = DualLinkedQueue{Int}()
    push!(q, 111)
    push!(q, 222)
    itr = NodeIterator(q)
    @test eltype(itr) === SSQNode{Int}
    nodes = collect(itr)
    @test nodes[1].value == 111
    @test nodes[2].value == 222
    @test check_invariance(q)
end

function test_print()
    q = DualLinkedQueue{Int}()
    push!(q, 111)
    push!(q, 222)
    str = sprint(show, "text/plain", q)
    @test "DualLinkedQueue" ⊏ str
end

end  # module


text:
"""
	contourf(cmd0::String="", arg1=nothing; kwargs...)

Performs Delaunay triangulation on x,y[,z] data, i.e., it find how the points should be connected
to give the most equilateral triangulation possible.

Full option list at [`triangulate`]($(GMTdoc)triangulate.html)

Parameters
----------
- **A** | **annot** :: [Type => Str | Number]       ``Arg = [-|[+]annot_int][labelinfo]``

    *annot_int* is annotation interval in data units; it is ignored if contour levels are given in a file.
    ($(GMTdoc)contour.html#a)
- $(GMT.opt_B)
- **C** | **cont** | **contour** | **contours** | **levels** :: [Type => Str | Number | GMTcpt]  ``Arg = [+]cont_int``

    Contours to be drawn may be specified in one of three possible ways.
    ($(GMTdoc)grdcontour.html#c)
- **E** | **index** :: [Type => Str | Mx3 array]

    Give name of file with network information. Each record must contain triplets of node
    numbers for a triangle.
    ($(GMTdoc)contour.html#e)
- **G** | **labels** :: [Type => Str]

    Controls the placement of labels along the quoted lines.
    ($(GMTdoc)contour.html#g)
- $(GMT.opt_J)
- $(GMT.opt_P)
- **Q** | **cut** :: [Type => Str | Number]         ``Arg = [cut[unit]][+z]]``

    Do not draw contours with less than cut number of points.
    ($(GMTdoc)contour.html#q)
- $(GMT.opt_R)
- **S** | **skip** :: [Type => Str | []]            ``Arg = [p|t]``

    Skip all input xyz points that fall outside the region (Used when input data is a table).
    ($(GMTdoc)contour.html#s)
- **S** | **smooth** :: [Type => Number]

    Used to resample the contour lines at roughly every (gridbox_size/smoothfactor) interval.
    (Used when input data is a grid)
    ($(GMTdoc)grdcontour.html#s)
- **T** | **ticks** :: [Type => Str]                 ``Arg = [+|-][+a][+dgap[/length]][+l[labels]]``

    Draw tick marks pointing in the downward direction every *gap* along the innermost closed contours.
    ($(GMTdoc)contour.html#t)
- $(GMT.opt_U)
- $(GMT.opt_V)
- **W** | **pen** :: [Type => Str | Number]

    Sets the attributes for the particular line.
    ($(GMTdoc)contour.html#w)
- **Z** | **xyz** | **triplets** :: [Type => Bool]

    ($(GMTdoc)triangulate.html#z)
- $(GMT.opt_bi)
- $(GMT.opt_bo)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_swap_xy)

Examples
--------

```julia
    G = GMT.peaks();
    C = makecpt(T=(-7,9,2));

    contourf(G, show=1)
    contourf(G, C=[-2, 0, 2, 5], show=1)
    contourf(G, C, contour=[-2, 0, 2, 5], show=1)
    contourf(G, C, annot=[-2, 0, 2, 5], show=1)
    contourf(G, C, annot=2, show=1)
    contourf(G, C, contour=1, annot=[-2, 0, 2, 5], show=1)
    contourf(G, C, annot=:none, show=1)

    d = [0 2 5; 1 4 5; 2 0.5 5; 3 3 9; 4 4.5 5; 4.2 1.2 5; 6 3 1; 8 1 5; 9 4.5 5];
    contourf(d, limits=(-0.5,9.5,0,5), pen=0.25, labels=(line=(:min,:max),), show=1)
```
"""
function contourf(cmd0::String="", arg1=nothing, arg2=nothing; first=true, kwargs...)

	d = KW(kwargs)
	dict_auto_add!(d)					# The ternary module may send options via another channel
	CPT_arg = (isa(arg1, GMTcpt)) ? arg1 : (isa(arg2, GMTcpt) ? arg2 : nothing)		# Fish a CPT, if any.

	CPT = nothing;		C_contours = "";	C_int = 0;
	if ((val = find_in_dict(d, [:C :cont :contour :contours :levels])[1]) !== nothing)
		if (isa(val, GMTcpt))
			CPT = val
		elseif (isa(val, Array{<:Number}))
			if (CPT_arg === nothing)	# No CPT yet, compute one
				CPT = makecpt(T=arg2str(val,','), M=true, par=(COLOR_BACKGROUND="white", COLOR_FOREGROUND="white"))
			else						# We already have one CPT so VAL is interpreted as contours values
				C_contours = arg2str(val,',')
				if (!occursin(",", C_contours))  C_contours *= ","  end
			end
		else
			t = string(val)
			if     (occursin(".cpt", t)) CPT = gmtread(t)
			elseif ((x = parse(Float64, t)) != NaN)  C_int = x
			else   error("Bad levels option")
			end
		end
	end

	if (cmd0 != "")		# Then it must be a file name. Of what?
		if     ((val = find_in_dict(d, [:grd :grid])[1]) !== nothing) X = gmtread(cmd0, grd=true)
		elseif ((val = find_in_dict(d, [:data :dataset :table])[1]) !== nothing) X = gmtread(cmd0, data=true)
		elseif ((val = find_in_dict(d, [:ogr])[1]) !== nothing)  X = gmtread(cmd0, ogr=true)
		else                                                     X = gmtread(cmd0, ignore_grd=1) # Try luck with guessing
		end
		if (X !== nothing)		# X === nothing when a recognized grid was passed by name
			if (arg1 !== nothing)  arg2 = arg1  end		# Then arg1 MUST be a CPT
			arg1 = X
			cmd0 = ""
		end
	end

	if (isa(arg1, GMTgrid) || cmd0 != "")
		# All of these are grdcontour options also, so must fish them, store and use later.
		opt_A = find_in_dict(d, [:A :annot])[1]			# Will either have something or nothing
		opt_G = find_in_dict(d, [:G :labels])[1]
		opt_T = find_in_dict(d, [:T :ticks])[1]
		opt_L = find_in_dict(d, [:L :range])[1]
		opt_Q = find_in_dict(d, [:Q :cut])[1]
		opt_S = find_in_dict(d, [:S :smooth])[1]
	end
	opt_W = find_in_dict(d, [:W :pen])[1]

	if (isa(arg1, GMTgrid) || cmd0 != "")
		#isa(arg2, GMTcpt) ? d[:N] = arg2 : (isa(arg1, GMTcpt) ? d[:N] = arg1 : d[:N] = true)
		if (CPT === nothing && CPT_arg === nothing)
			if (cmd0 != "")
				info = grdinfo(cmd0 * " -C")
				C_inc, min, max = gen_contour_vals(info[1].data[5:6], C_int)
			else
				C_inc, min, max = gen_contour_vals(arg1, C_int)
			end
			CPT = makecpt(T=(min, max, C_inc))
		elseif (CPT === nothing && CPT_arg !== nothing)
			CPT = CPT_arg
		end

		show = false;	fmt = nothing;	savefig = ""
		if ((val = find_in_dict(d, [:show])[1]) !== nothing)  show = true  end
		if ((val = find_in_dict(d, [:fmt])[1])  !== nothing)  fmt  = val   end
		if ((val = find_in_dict(d, [:savefig :figname :name])[1]) !== nothing)  savefig = val  end

		d[:C] = CPT;	d[:Q] = "s";	done = false#	d[:W] = "c"
		if ((C_int == 0) && (opt_A == "none" || opt_A == :none))	# Then no need grdcontour
			if (show)  d[:show] = true  end
			if (fmt !== nothing) d[:fmt] = fmt  end
			if (savefig != "")   d[:savefig] = savefig  end
			if (opt_W !== nothing)  d[:W] = opt_W  end
			d[:W] = "c";  done = true
		end
		grdview(cmd0, arg1; first=first, d...)
		del_from_dict(d, [[:C :z], [:Q], [:W]])
		if (done)  return  end

		if (C_int > 0 || C_contours != "" || opt_A !== nothing)
			if     (C_int > 0)        d[:C] = string(C_int)
			elseif (C_contours != "") d[:C] = C_contours
			end
			if (opt_A !== nothing)    d[:A] = opt_A  end
		else
			d[:C] = CPT;
		end

		if (show)  d[:show] = true  end
		if (fmt !== nothing) d[:fmt] = fmt  end
		if (savefig != "")   d[:savefig] = savefig  end
		if (opt_G !== nothing)  d[:G] = opt_G  end
		if (opt_L !== nothing)  d[:L] = opt_L  end
		if (opt_Q !== nothing)  d[:Q] = opt_Q  end
		if (opt_T !== nothing)  d[:T] = opt_T  end
		if (opt_S !== nothing)  d[:S] = opt_S  end
		if (opt_W !== nothing)  d[:W] = opt_W  end
		grdcontour(cmd0, arg1; first=false, d...)
	else
		if (isa(CPT_arg, GMTcpt))
			d[:C] = CPT_arg;
		elseif (isa(CPT, GMTcpt))
			d[:C] = CPT;
		else
			D = gmtinfo(arg1, C=1)
			C_inc, min, max = gen_contour_vals(D[1].data[5:6], C_int)
			d[:C] = makecpt(T=(min, max, C_inc))
		end
		d[:I] = true
		(C_int != 0 && opt_W === nothing) && (opt_W = "0.25p")
		(opt_W !== nothing) && (d[:W] = opt_W)
		contour("", arg1; first=first, d...)
	end

end

# ---------------------------------------------------------------------------------------------------
function gen_contour_vals(in, C_int)
	# If C_int (contour interval) != 0 use it otherwise guess one
	if (isa(in, GMTgrid))  low = in.range[5];	high = in.range[6]
	else                   low = in[1];		high = in[2]
	end
	range = high - low
	(C_int == 0) && (C_int = auto_contour_interval(range))
	min = floor(low / C_int) * C_int;
	if (min > low)  min -= C_int  end
	max = ceil(high / C_int) * C_int;
	if (max < high)  max += C_int  end
	return C_int, min, max
end

# ---------------------------------------------------------------------------------------------------
function auto_contour_interval(range)
	# Do the same as GMT C
	x  = 10 ^ (floor(log10(range)) - 1.0)
	nx = div(range, x, RoundNearest)
	x *= (nx > 40) ? 5 : (nx > 20 ? 2 : 1)
	return x		# Contour intervals
end

# ---------------------------------------------------------------------------------------------------
contourf!(cmd0::String="", arg1=nothing, arg2=nothing; kw...) = contourf(cmd0, arg1, arg2; first=false, kw...)
contourf(arg1, arg2=nothing; kw...) = contourf("", arg1, arg2; first=true, kw...)
contourf!(arg1, arg2=nothing; kw...) = contourf("", arg1, arg2; first=false, kw...)


text:
using ClangSA

fixeddb = ClangSA.createFixedCompileDatabase(".", ["cc",
  string("-I", joinpath(JULIA_HOME,"../../src")),
  string("-I", joinpath(JULIA_HOME,"../../src/support")),
  string("-I", joinpath(JULIA_HOME,"../include/julia")),
  ])

ClangSA.runTool(fixeddb, ["GCPushPop.cpp"], ClangSA.getAnalysisAction(), verify = true)

#ClangSA.runTool(fixeddb, [path], ClangSA.getAnalysisAction())


text:
using Pkg

Pkg.activate(@__DIR__)
Pkg.instantiate()

using PackageCompiler

cd(@__DIR__)

create_sysimage(:RoME, sysimage_path="RoMESysimage.so", precompile_execution_file="precompile_triggers.jl")

## to use RoME with the newly created sysimage, start julia with:
# julia -J RoMESysimage.so


text:
using Touchstone
TS = Touchstone
using Test

include( "TestTouchstone.jl" )
include( "TestParseTouchstone.jl" )
include( "TestWriteTouchstone.jl" )
include( "TestFilesTouchstone.jl" )


text:
# PyNeuralPose2Pose2

using Random, Statistics
using DistributedFactorGraphs, TransformUtils
@everywhere using Random, Statistics
@everywhere using DistributedFactorGraphs, TransformUtils

@everywhere begin
import Base: convert
import IncrementalInference: getSample

struct PyNeuralPose2Pose2{P,D<:Vector,M<:SamplableBelief} <: FunctorPairwise
  predictFnc::P
  joyVelData::D
  naiveModel::M
  naiveFrac::Float64
  Zij::Pose2Pose2
  specialSampler::Function # special keyword field name used to invoke 'specialSampler' logic
end

function sampleNeuralPose2(nfb::PyNeuralPose2Pose2,
                          N::Int,
                          fmd::FactorMetadata,
                          Xi::DFGVariable,
                          Xj::DFGVariable)::Tuple
 #

  # calculate naive model and Predictive fraction of samples, respectively
  Nn = round(Int, nfb.naiveFrac*N)
  # calculate desired number of predicted values
  Np = N - Nn
  len = length(nfb.joyVelData) # expect this to be only 25 at developmenttime, likely to change

  # model samples (all for theta at this time)
  smpls_mAll = rand(nfb.naiveModel, N)

  # sample predictive fraction
  iT, jT = getTimestamp(Xi), getTimestamp(Xj)
  iPts, jPts = (getKDE(Xi) |> getPoints), (getKDE(Xj) |> getPoints)
  @assert size(jPts,2) == size(iPts,2) "sampleNeuralPose2 can currently only evaluate equal population   size variables"

  # calculate an average velocity component
  DT = jT - iT
  DXY = (@view jPts[1:2,:]) - (@view jPts[1:2,:])
  # rotate delta position from world to local iX frame
  for i in 1:size(iPts,2)
    DXY[1:2,i] .= TransformUtils.R(iPts[3,i])'*DXY[1:2,i]
  end
  # replace delta (velocity) values for this sampling
  mVXY = Statistics.mean(DXY, dims=2)
  # divide time to get velocity
  mVXY ./= 1e-3*DT.value
  mVXY[1] = isnan(mVXY[1]) ? 0.0 : mVXY[1]
  mVXY[2] = isnan(mVXY[2]) ? 0.0 : mVXY[2]

  for i in 1:len
    nfb.joyVelData[i][3:4] = mVXY
  end
  # and predict
      # A = [rand(4) for i in 1:25]
  smpls_pAll = nfb.predictFnc(nfb.joyVelData)

  # number of predictors to choose from, and choose random subset
  Npreds = size(smpls_pAll,1)
  allPreds = 1:Npreds |> collect
  # randomly select particles for prediction (with possible duplicates forwhen Np > size(iPts,2))
  Npp = Np < Npreds ? Np : Npreds
  Nnn = N - Npp
  selPreds = @view shuffle!(allPreds)[1:Npp] # TODO better in-place
  smpls_p = smpls_pAll[selPreds,:] # sample per row??
  smpls_p[:,3] .= smpls_mAll[3,Nnn+1:N] # use naive delta theta at this time

  # naive fraction
  smpls_n = @view smpls_mAll[:, 1:Nnn] # sample per column
  # join and shuffle predicted odo values
  shfSmpl = shuffle!(1:N |> collect)
  smpls = hcat(smpls_n, smpls_p')[:,shfSmpl]

  # @show Statistics.mean(smpls, dims=2)
  # @show maximum(smpls, dims=2)
  # @show minimum(smpls, dims=2)

  return (smpls, )
end

# Convenience function to help call the right constuctor
PyNeuralPose2Pose2(nn::P,
                   jvd::D,
                   md::M,
                   naiveFrac::Float64=0.4,
                   ss::Function=sampleNeuralPose2) where {P, M <: SamplableBelief, D <: Vector} = PyNeuralPose2Pose2{P,D,M}(nn,jvd,md,naiveFrac,Pose2Pose2(MvNormal(zeros(3),diagm(ones(3)))),ss )
#

function (nfb::PyNeuralPose2Pose2)(
            res::AbstractArray{<:Real},
            userdata::FactorMetadata,
            idx::Int,
            meas::Tuple{AbstractArray{<:Real},},
            Xi::AbstractArray{<:Real,2},
            Xj::AbstractArray{<:Real,2}  )
  #
  nfb.Zij(res,userdata,idx,meas,Xi,Xj)
  nothing
end

## packing converters

struct PackedPyNeuralPose2Pose2 <: IncrementalInference.PackedInferenceType
  joyVelData::Vector{Vector{Float64}}
  naiveModel::String
  naiveFrac::Float64
end

function convert(::Type{PyNeuralPose2Pose2}, d::PackedPyNeuralPose2Pose2)
  PyNeuralPose2Pose2(PyTFOdoPredictorPoint2,d.joyVelData,extractdistribution(d.naiveModel),d.naiveFrac)
end

function convert(::Type{PackedPyNeuralPose2Pose2}, d::PyNeuralPose2Pose2)
  PackedPyNeuralPose2Pose2(d.joyVelData, string(d.naiveModel), d.naiveFrac)
end

end # everywhere

#


text:
function build_ZYB!( branches, num_branches, num_dyns )  #=

builds the matrices Z, Y, B1 and B2 which are submatrices of the branch constitutive equations (BCEs) of the sparse tableau formulation (STF). They are built from the information of branches given as the arguments of this function. The first argument "branches" is an array of the struct Branch containing the information of a branch. The second argument "num_branches" and the third argument "num_dyns" are the number of branches and that of dynamic branches respectively, where an m-phase branch is counted as m branches. This function will return the matrices Z, Y, B1 and B2.

Author: Taku Noda
Started on May 10, 2019
Update on June 29, 2021

Copyright 2021 Central Research Institute of Electric Power Industry

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

=#

# -- Internally used functions -----------------------------------------

function set_meye!( A, ipos, jpos, nph )
for iph = 0:( nph - 1 )
  A[ ipos + iph, jpos + iph ] = -1.0
end
end

function calc_slope0( x, y )
N = length( x )
if x[ 1 ] > 0.0
  a = ( y[ 2 ] - y[ 1 ] )/( x[ 2 ] - x[ 1 ] )
  op1 = 0
elseif x[ N ] < 0.0
  a = ( y[ N ] - y[ N - 1 ] )/( x[ N ] - x[ N - 1 ] )
  op1 = N
else
  for n = 2:N
    if x[ n ] >= 0.0
      a = ( y[ n ] - y[ n - 1 ] )/( x[ n ] - x[ n - 1 ] )
      op1 = n - 1
      break
    end
  end
end
return a, op1
end

# -- Main body starts here ---------------------------------------------

Z = spzeros( num_branches, num_branches )
Y = spzeros( num_branches, num_branches )
B1 = spzeros( num_branches, num_dyns )
B2 = spzeros( num_dyns, 2*num_branches )

for branch in branches
  kpos = branch.p_vars
  nph = branch.nph
  mpos = kpos + nph - 1
  if branch.type == "R"
    Z[ kpos:mpos, kpos:mpos ] = eval( branch.val )
    set_meye!( Y, kpos, kpos, nph )
  elseif branch.type == "L"
    Z[ kpos:mpos, kpos:mpos ] = eval( branch.val )
    ppos = branch.p_de
    set_meye!( B1, kpos, ppos, nph )
    set_meye!( B2, ppos, num_branches + kpos, nph )
  elseif branch.type == "G"
    Y[ kpos:mpos, kpos:mpos ] = eval( branch.val )
    set_meye!( Z, kpos, kpos, nph )
  elseif branch.type == "C"
    Y[ kpos:mpos, kpos:mpos ] = eval( branch.val )
    ppos = branch.p_de
    set_meye!( B1, kpos, ppos, nph )
    set_meye!( B2, ppos, kpos, nph )
  elseif branch.type == "NR"
    Z[ kpos, kpos ], op1 = calc_slope0( branch.plx, branch.ply )
    ( branch.op )[ 1 ] = op1
    Y[ kpos, kpos ] = -1.0
  elseif branch.type == "NL"
    Z[ kpos, kpos ], op1 = calc_slope0( branch.plx, branch.ply )
    ( branch.op )[ 1 ] = op1
    ppos = branch.p_de
    B1[ kpos, ppos ] = -1
    B2[ ppos, num_branches + kpos ] = -1
  elseif branch.type == "NC"
    Y[ kpos, kpos ], op1 = calc_slope0( branch.plx, branch.ply )
    ( branch.op )[ 1 ] = op1
    ppos = branch.p_de
    B1[ kpos, ppos ] = -1
    B2[ ppos, kpos ] = -1
  elseif branch.type == "E"
    Y[ kpos, kpos ] = 1.0
  elseif branch.type == "J"
    Z[ kpos, kpos ] = 1.0
  end
end

return Z, Y, B1, B2

end

max_stars_repo_name: takunoda/omegasim

max_issues_repo_name: takunoda/omegasim

max_forks_repo_name: takunoda/omegasim


text:
# Autogenerated wrapper script for nauty_jll for x86_64-linux-musl
export NRswitchg, addedgeg, amtog, biplabg, blisstog, catg, checks6, complg, converseg, copyg, countg, cubhamg, deledgeg, delptg, directg, dreadnaut, dretodot, dretog, genbg, genbgL, geng, genquarticg, genrang, genspecialg, gentourng, gentreeg, hamheuristic, labelg, libnauty, libnautyA1, libnautyL0, libnautyL1, libnautyS0, libnautyS1, libnautyW0, libnautyW1, linegraphg, listg, multig, newedgeg, pickg, planarg, ranlabg, shortg, showg, subdivideg, sumlines, twohamg, vcolg, watercluster2

using GMP_jll
JLLWrappers.@generate_wrapper_header("nauty")
JLLWrappers.@declare_executable_product(NRswitchg)
JLLWrappers.@declare_executable_product(addedgeg)
JLLWrappers.@declare_executable_product(amtog)
JLLWrappers.@declare_executable_product(biplabg)
JLLWrappers.@declare_executable_product(blisstog)
JLLWrappers.@declare_executable_product(catg)
JLLWrappers.@declare_executable_product(checks6)
JLLWrappers.@declare_executable_product(complg)
JLLWrappers.@declare_executable_product(converseg)
JLLWrappers.@declare_executable_product(copyg)
JLLWrappers.@declare_executable_product(countg)
JLLWrappers.@declare_executable_product(cubhamg)
JLLWrappers.@declare_executable_product(deledgeg)
JLLWrappers.@declare_executable_product(delptg)
JLLWrappers.@declare_executable_product(directg)
JLLWrappers.@declare_executable_product(dreadnaut)
JLLWrappers.@declare_executable_product(dretodot)
JLLWrappers.@declare_executable_product(dretog)
JLLWrappers.@declare_executable_product(genbg)
JLLWrappers.@declare_executable_product(genbgL)
JLLWrappers.@declare_executable_product(geng)
JLLWrappers.@declare_executable_product(genquarticg)
JLLWrappers.@declare_executable_product(genrang)
JLLWrappers.@declare_executable_product(genspecialg)
JLLWrappers.@declare_executable_product(gentourng)
JLLWrappers.@declare_executable_product(gentreeg)
JLLWrappers.@declare_executable_product(hamheuristic)
JLLWrappers.@declare_executable_product(labelg)
JLLWrappers.@declare_library_product(libnauty, "libnauty.so.2")
JLLWrappers.@declare_library_product(libnautyA1, "libnautyA1.so.2")
JLLWrappers.@declare_library_product(libnautyL0, "libnautyL0.so.2")
JLLWrappers.@declare_library_product(libnautyL1, "libnautyL1.so.2")
JLLWrappers.@declare_library_product(libnautyS0, "libnautyS0.so.2")
JLLWrappers.@declare_library_product(libnautyS1, "libnautyS1.so.2")
JLLWrappers.@declare_library_product(libnautyW0, "libnautyW0.so.2")
JLLWrappers.@declare_library_product(libnautyW1, "libnautyW1.so.2")
JLLWrappers.@declare_executable_product(linegraphg)
JLLWrappers.@declare_executable_product(listg)
JLLWrappers.@declare_executable_product(multig)
JLLWrappers.@declare_executable_product(newedgeg)
JLLWrappers.@declare_executable_product(pickg)
JLLWrappers.@declare_executable_product(planarg)
JLLWrappers.@declare_executable_product(ranlabg)
JLLWrappers.@declare_executable_product(shortg)
JLLWrappers.@declare_executable_product(showg)
JLLWrappers.@declare_executable_product(subdivideg)
JLLWrappers.@declare_executable_product(sumlines)
JLLWrappers.@declare_executable_product(twohamg)
JLLWrappers.@declare_executable_product(vcolg)
JLLWrappers.@declare_executable_product(watercluster2)
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll)
    JLLWrappers.@init_executable_product(
        NRswitchg,
        "bin/NRswitchg",
    )

    JLLWrappers.@init_executable_product(
        addedgeg,
        "bin/addedgeg",
    )

    JLLWrappers.@init_executable_product(
        amtog,
        "bin/amtog",
    )

    JLLWrappers.@init_executable_product(
        biplabg,
        "bin/biplabg",
    )

    JLLWrappers.@init_executable_product(
        blisstog,
        "bin/blisstog",
    )

    JLLWrappers.@init_executable_product(
        catg,
        "bin/catg",
    )

    JLLWrappers.@init_executable_product(
        checks6,
        "bin/checks6",
    )

    JLLWrappers.@init_executable_product(
        complg,
        "bin/complg",
    )

    JLLWrappers.@init_executable_product(
        converseg,
        "bin/converseg",
    )

    JLLWrappers.@init_executable_product(
        copyg,
        "bin/copyg",
    )

    JLLWrappers.@init_executable_product(
        countg,
        "bin/countg",
    )

    JLLWrappers.@init_executable_product(
        cubhamg,
        "bin/cubhamg",
    )

    JLLWrappers.@init_executable_product(
        deledgeg,
        "bin/deledgeg",
    )

    JLLWrappers.@init_executable_product(
        delptg,
        "bin/delptg",
    )

    JLLWrappers.@init_executable_product(
        directg,
        "bin/directg",
    )

    JLLWrappers.@init_executable_product(
        dreadnaut,
        "bin/dreadnaut",
    )

    JLLWrappers.@init_executable_product(
        dretodot,
        "bin/dretodot",
    )

    JLLWrappers.@init_executable_product(
        dretog,
        "bin/dretog",
    )

    JLLWrappers.@init_executable_product(
        genbg,
        "bin/genbg",
    )

    JLLWrappers.@init_executable_product(
        genbgL,
        "bin/genbgL",
    )

    JLLWrappers.@init_executable_product(
        geng,
        "bin/geng",
    )

    JLLWrappers.@init_executable_product(
        genquarticg,
        "bin/genquarticg",
    )

    JLLWrappers.@init_executable_product(
        genrang,
        "bin/genrang",
    )

    JLLWrappers.@init_executable_product(
        genspecialg,
        "bin/genspecialg",
    )

    JLLWrappers.@init_executable_product(
        gentourng,
        "bin/gentourng",
    )

    JLLWrappers.@init_executable_product(
        gentreeg,
        "bin/gentreeg",
    )

    JLLWrappers.@init_executable_product(
        hamheuristic,
        "bin/hamheuristic",
    )

    JLLWrappers.@init_executable_product(
        labelg,
        "bin/labelg",
    )

    JLLWrappers.@init_library_product(
        libnauty,
        "lib/libnauty.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnautyA1,
        "lib/libnautyA1.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnautyL0,
        "lib/libnautyL0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnautyL1,
        "lib/libnautyL1.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnautyS0,
        "lib/libnautyS0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnautyS1,
        "lib/libnautyS1.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnautyW0,
        "lib/libnautyW0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnautyW1,
        "lib/libnautyW1.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        linegraphg,
        "bin/linegraphg",
    )

    JLLWrappers.@init_executable_product(
        listg,
        "bin/listg",
    )

    JLLWrappers.@init_executable_product(
        multig,
        "bin/multig",
    )

    JLLWrappers.@init_executable_product(
        newedgeg,
        "bin/newedgeg",
    )

    JLLWrappers.@init_executable_product(
        pickg,
        "bin/pickg",
    )

    JLLWrappers.@init_executable_product(
        planarg,
        "bin/planarg",
    )

    JLLWrappers.@init_executable_product(
        ranlabg,
        "bin/ranlabg",
    )

    JLLWrappers.@init_executable_product(
        shortg,
        "bin/shortg",
    )

    JLLWrappers.@init_executable_product(
        showg,
        "bin/showg",
    )

    JLLWrappers.@init_executable_product(
        subdivideg,
        "bin/subdivideg",
    )

    JLLWrappers.@init_executable_product(
        sumlines,
        "bin/sumlines",
    )

    JLLWrappers.@init_executable_product(
        twohamg,
        "bin/twohamg",
    )

    JLLWrappers.@init_executable_product(
        vcolg,
        "bin/vcolg",
    )

    JLLWrappers.@init_executable_product(
        watercluster2,
        "bin/watercluster2",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# ---
# title: 480. Sliding Window Median
# id: problem480
# author: Indigo
# date: 2021-02-03
# difficulty: Hard
# categories: Sliding Window
# link: <https://leetcode.com/problems/sliding-window-median/description/>
# hidden: true
# ---
#
# Median is the middle value in an ordered integer list. If the size of the list
# is even, there is no middle value. So the median is the mean of the two middle
# value.
#
# Examples:
#
# `[2,3,4]` , the median is `3`
#
# `[2,3]`, the median is `(2 + 3) / 2 = 2.5`
#
# Given an array _nums_ , there is a sliding window of size _k_ which is moving
# from the very left of the array to the very right. You can only see the _k_
# numbers in the window. Each time the sliding window moves right by one
# position. Your job is to output the median array for each window in the
# original array.
#
# For example,
# Given _nums_ = `[1,3,-1,-3,5,3,6,7]`, and _k_ = 3.
#
#
#
#     Window position                Median
#     ---------------               -----
#     [1  3  -1] -3  5  3  6  7       1
#      1 [3  -1  -3] 5  3  6  7       -1
#      1  3 [-1  -3  5] 3  6  7       -1
#      1  3  -1 [-3  5  3] 6  7       3
#      1  3  -1  -3 [5  3  6] 7       5
#      1  3  -1  -3  5 [3  6  7]      6
#
#
# Therefore, return the median sliding window as `[1,-1,-1,3,5,6]`.
#
# **Note:**
# You may assume `k` is always valid, ie: `k` is always smaller than input
# array's size for non-empty array.
# Answers within `10^-5` of the actual value will be accepted as correct.
#
#
## @lc code=start
using LeetCode

function median_sliding_window(nums::Vector{Int}, k::Int)
    tree = AVLTree{Tuple{Int, Int}}()
    res = Float64[]
    for i in 1:length(nums)
        push!(tree, (nums[i], i))
        if i ≥ k
            push!(res, (tree[k ÷ 2 + 1][1] + tree[(k - 1) ÷ 2 + 1][1]) / 2)
            delete!(tree, (nums[i - k + 1], i - k + 1))
        end
    end
    res
end
## @lc code=end


text:
"""
    DependenciesParser.jl

    This package provides a quick way to access dependency information
"""
module DependenciesParser
    using Base.Iterators: flatten
    using Pkg: METADATA_compatible_uuid
    using Pkg.Operations: load_package_data_raw, deps_graph, simplify_graph!, resolve
    using Pkg.TOML: parsefile
    using Pkg.Types: Context, Fixed, Requires, UUID, uuid_julia, VersionRange, VersionSpec

    """
        All package names in the General registry
    """
    const data =
        readdir.(joinpath.(homedir(), ".julia/registries/General", string.('A':'Z'))) |>
        flatten |>
        collect |>
        (x -> filter!(x -> ~any(x ∈ ["julia", ".DS_Store"]), x))
    function find_repo(name)
        dir = joinpath(homedir(), ".julia/registries/General", uppercase(name[1:1]), name)
        toml = parsefile(joinpath(dir, "Package.toml"))
        repo = replace(toml["repo"], r"\.git$" => "")
        try
            request("GET", repo)
            true
        catch
            false
        end
    end
    # Identify deleted repositories
    # using HTTP: request
    # available = Vector{Bool}(undef, length(data))
    # @time for idx ∈ eachindex(available)
    #     println(idx)
    #     available[idx] = find_repo(data[idx])
    # end
    # Based on a cache solution from above on 2019-01-25
    """
        Packages that no longer exist (repositories have been deleted)
    """
    const deleted_repo =
        ["Arduino", "ChainRecursive", "Chunks", "CombinatorialBandits", "ControlCore",
         "DotOverloading", "DynamicalBilliardsPlotting", "GLUT", "GetC", "HTSLIB",
         "KeyedTables", "LazyCall", "LazyContext", "LazyQuery", "LibGit2", "NumberedLines",
         "OpenGL", "OrthogonalPolynomials", "Parts", "React", "RecurUnroll",
         "RequirementVersions", "SDL", "SessionHacker", "Sparrow", "StringArrays",
         "TypedBools", "ValuedTuples", "ZippedArrays"]
    filter!(pkg -> pkg ∉ deleted_repo, data)
    const deps = Dict{UUID,Dict{VersionRange,Dict{String,UUID}}}()
    const compat = Dict{UUID,Dict{VersionRange,Dict{String,VersionSpec}}}()
    const uuid_to_name =
        Dict(UUID("ade2ca70-3891-5945-98fb-dc099432e06a")=>"Dates",
             UUID("7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee")=>"FileWatching",
             UUID("ea8e919c-243c-51af-8825-aaa63cd721ce")=>"SHA",
             UUID("6462fe0b-24de-5631-8697-dd941f90decc")=>"Sockets",
             UUID("8bb1440f-4735-579b-a4ab-409b98df4dab")=>"DelimitedFiles",
             UUID("9a3f8284-a2c9-5f02-9a11-845980a1fd5c")=>"Random",
             UUID("8ba89e20-285c-5b6f-9357-94700520ee1b")=>"Distributed",
             UUID("37e2e46d-f89d-539d-b4ee-838fcccc9c8e")=>"LinearAlgebra",
             UUID("a63ad114-7e13-5084-954f-fe012c677804")=>"Mmap",
             UUID("2f01184e-e22b-5df5-ae63-d93ebab69eaf")=>"SparseArrays",
             UUID("cf7118a7-6976-5b1a-9a39-7adc72f591a4")=>"UUIDs",
             UUID("56ddb016-857b-54e1-b83d-db4d58db5568")=>"Logging",
             UUID("10745b16-79ce-11e8-11f9-7d13ad32a3b2")=>"Statistics",
             UUID("8bf52ea8-c179-5cab-976a-9e18b702a9bc")=>"CRC32c",
             UUID("44cfe95a-1eb2-52ea-b672-e2afdf69b78f")=>"Pkg",
             UUID("9e88b42a-f829-5b0c-bbe9-9e923198166b")=>"Serialization",
             UUID("d6f4376e-aef5-505a-96c1-9c027394607a")=>"Markdown",
             UUID("4607b0f0-06f3-5cda-b6b1-a6196a1729e9")=>"SuiteSparse",
             UUID("3fa0cd96-eef1-5676-8a61-b3b8758bbffb")=>"REPL",
             UUID("8f399da3-3557-5675-b5ff-fb832c97cbdb")=>"Libdl",
             UUID("2a0f44e3-6c83-55bd-87e4-b1978d98bd5f")=>"Base64",
             UUID("9abbd945-dff8-562f-b5e8-e1ebf5ef1b79")=>"Profile",
             UUID("b77e0a4c-d291-57a0-90e8-8db25a27a240")=>"InteractiveUtils",
             UUID("1a1011a3-84de-559e-8e89-a11a2f7dc383")=>"SharedArrays",
             UUID("de0858da-6303-5e67-8744-51eddeeeb8d7")=>"Printf",
             UUID("9fa8497b-333b-5362-9e8d-4d0656e87820")=>"Future",
             UUID("76f85450-5226-5b5a-8eaa-529ad045b433")=>"LibGit2",
             UUID("8dfed614-e22c-5e08-85e1-65c5234f0b40")=>"Test",
             UUID("4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5")=>"Unicode",
             uuid_julia=>"julia")
    const versions = Dict{UUID,Set{VersionNumber}}()
    for name ∈ data
        dir = joinpath(homedir(), ".julia/registries/General", uppercase(name[1:1]), name)
        uuid = UUID(parsefile(joinpath(dir, "Package.toml"))["uuid"])
        uuid_to_name[uuid] = name
        versions[uuid] = Set(VersionNumber.(keys(parsefile(joinpath(dir, "Versions.toml")))))
        deps[uuid] = load_package_data_raw(UUID, joinpath(dir, "Deps.toml"))
        compat[uuid] = load_package_data_raw(VersionSpec, joinpath(dir, "Compat.toml"))
    end
    """
        installable(pkg::AbstractString,
                    julia::VersionNumber = VERSION;
                    direct::Bool = false)::Tuple{Bool,Vector{String}}

        Return whether the package is installable and the dependencies for the solved version.
        If direct, only direct dependencies are returned.
    """
    function installable(pkg::AbstractString,
                         julia = VERSION::VersionNumber;
                         direct::Bool = false)
        uuid = METADATA_compatible_uuid(pkg)
        try
            graph = deps_graph(Context(),
                               uuid_to_name,
                               Requires(uuid => VersionSpec()),
                               Dict(uuid_julia => Fixed(julia)))
            simplify_graph!(graph)
            sol = get.(Ref(uuid_to_name),
                       filter(!isequal(uuid), keys(resolve(graph))),
                       nothing) |>
                  sort!
            if direct
                secondary = reduce((x,y) -> vcat(last(x), last(y)),
                                   installable.(sol)) |>
                            unique!
                sort!(filter!(dep -> dep ∉ secondary, sol))
            end
            return true, sol::Vector{String}
        catch
            return false, Vector{String}([pkg])
        end
    end
    # Code to get all installable packages
    # status = Vector{Tuple{Bool,Vector{String}}}()
    # @time for (idx, pkg) ∈ enumerate(data)
        # println(idx)
        # push!(status, installable(pkg))
    # end
    # data[first.(status)]
    # __init__() = append!(data, alldeps())
    export installable
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: cloudsearch_domain
using AWS.Compat
using AWS.UUIDs

"""
    Search()

Retrieves a list of documents that match the specified search criteria. How you specify the search criteria depends on which query parser you use. Amazon CloudSearch supports four query parsers:  simple: search all text and text-array fields for the specified string. Search for phrases, individual terms, and prefixes.  structured: search specific fields, construct compound queries using Boolean operators, and use advanced features such as term boosting and proximity searching. lucene: specify search criteria using the Apache Lucene query parser syntax. dismax: specify search criteria using the simplified subset of the Apache Lucene query parser syntax defined by the DisMax query parser.  For more information, see Searching Your Data in the Amazon CloudSearch Developer Guide. The endpoint for submitting Search requests is domain-specific. You submit search requests to a domain's search endpoint. To get the search endpoint for your domain, use the Amazon CloudSearch configuration service DescribeDomains action. A domain's endpoints are also displayed on the domain dashboard in the Amazon CloudSearch console.

# Required Parameters
- `q`: Specifies the search criteria for the request. How you specify the search criteria depends on the query parser used for the request and the parser options specified in the queryOptions parameter. By default, the simple query parser is used to process requests. To use the structured, lucene, or dismax query parser, you must also specify the queryParser parameter.  For more information about specifying search criteria, see Searching Your Data in the Amazon CloudSearch Developer Guide.

# Optional Parameters
- `cursor`: Retrieves a cursor value you can use to page through large result sets. Use the size parameter to control the number of hits to include in each response. You can specify either the cursor or start parameter in a request; they are mutually exclusive. To get the first cursor, set the cursor value to initial. In subsequent requests, specify the cursor value returned in the hits section of the response.  For more information, see Paginating Results in the Amazon CloudSearch Developer Guide.
- `expr`: Defines one or more numeric expressions that can be used to sort results or specify search or filter criteria. You can also specify expressions as return fields.  You specify the expressions in JSON using the form {\"EXPRESSIONNAME\":\"EXPRESSION\"}. You can define and use multiple expressions in a search request. For example:  {\"expression1\":\"_score*rating\", \"expression2\":\"(1/rank)*year\"}   For information about the variables, operators, and functions you can use in expressions, see Writing Expressions in the Amazon CloudSearch Developer Guide.
- `facet`: Specifies one or more fields for which to get facet information, and options that control how the facet information is returned. Each specified field must be facet-enabled in the domain configuration. The fields and options are specified in JSON using the form {\"FIELD\":{\"OPTION\":VALUE,\"OPTION:\"STRING\"},\"FIELD\":{\"OPTION\":VALUE,\"OPTION\":\"STRING\"}}. You can specify the following faceting options:   buckets specifies an array of the facet values or ranges to count. Ranges are specified using the same syntax that you use to search for a range of values. For more information, see  Searching for a Range of Values in the Amazon CloudSearch Developer Guide. Buckets are returned in the order they are specified in the request. The sort and size options are not valid if you specify buckets.   size specifies the maximum number of facets to include in the results. By default, Amazon CloudSearch returns counts for the top 10. The size parameter is only valid when you specify the sort option; it cannot be used in conjunction with buckets.   sort specifies how you want to sort the facets in the results: bucket or count. Specify bucket to sort alphabetically or numerically by facet value (in ascending order). Specify count to sort by the facet counts computed for each facet value (in descending order). To retrieve facet counts for particular values or ranges of values, use the buckets option instead of sort.    If no facet options are specified, facet counts are computed for all field values, the facets are sorted by facet count, and the top 10 facets are returned in the results. To count particular buckets of values, use the buckets option. For example, the following request uses the buckets option to calculate and return facet counts by decade.  {\"year\":{\"buckets\":[\"[1970,1979]\",\"[1980,1989]\",\"[1990,1999]\",\"[2000,2009]\",\"[2010,}\"]}}  To sort facets by facet count, use the count option. For example, the following request sets the sort option to count to sort the facet values by facet count, with the facet values that have the most matching documents listed first. Setting the size option to 3 returns only the top three facet values.  {\"year\":{\"sort\":\"count\",\"size\":3}}  To sort the facets by value, use the bucket option. For example, the following request sets the sort option to bucket to sort the facet values numerically by year, with earliest year listed first.   {\"year\":{\"sort\":\"bucket\"}}  For more information, see Getting and Using Facet Information in the Amazon CloudSearch Developer Guide.
- `fq`: Specifies a structured query that filters the results of a search without affecting how the results are scored and sorted. You use filterQuery in conjunction with the query parameter to filter the documents that match the constraints specified in the query parameter. Specifying a filter controls only which matching documents are included in the results, it has no effect on how they are scored and sorted. The filterQuery parameter supports the full structured query syntax.  For more information about using filters, see Filtering Matching Documents in the Amazon CloudSearch Developer Guide.
- `highlight`: Retrieves highlights for matches in the specified text or text-array fields. Each specified field must be highlight enabled in the domain configuration. The fields and options are specified in JSON using the form {\"FIELD\":{\"OPTION\":VALUE,\"OPTION:\"STRING\"},\"FIELD\":{\"OPTION\":VALUE,\"OPTION\":\"STRING\"}}. You can specify the following highlight options:   format: specifies the format of the data in the text field: text or html. When data is returned as HTML, all non-alphanumeric characters are encoded. The default is html.   max_phrases: specifies the maximum number of occurrences of the search term(s) you want to highlight. By default, the first occurrence is highlighted.   pre_tag: specifies the string to prepend to an occurrence of a search term. The default for HTML highlights is &lt;em&gt;. The default for text highlights is *.   post_tag: specifies the string to append to an occurrence of a search term. The default for HTML highlights is &lt;/em&gt;. The default for text highlights is *.   If no highlight options are specified for a field, the returned field text is treated as HTML and the first match is highlighted with emphasis tags: &lt;em>search-term&lt;/em&gt;. For example, the following request retrieves highlights for the actors and title fields.  { \"actors\": {}, \"title\": {\"format\": \"text\",\"max_phrases\": 2,\"pre_tag\": \"\",\"post_tag\": \"\"} }
- `partial`: Enables partial results to be returned if one or more index partitions are unavailable. When your search index is partitioned across multiple search instances, by default Amazon CloudSearch only returns results if every partition can be queried. This means that the failure of a single search instance can result in 5xx (internal server) errors. When you enable partial results, Amazon CloudSearch returns whatever results are available and includes the percentage of documents searched in the search results (percent-searched). This enables you to more gracefully degrade your users' search experience. For example, rather than displaying no results, you could display the partial results and a message indicating that the results might be incomplete due to a temporary system outage.
- `q.options`: Configures options for the query parser specified in the queryParser parameter. You specify the options in JSON using the following form {\"OPTION1\":\"VALUE1\",\"OPTION2\":VALUE2\"...\"OPTIONN\":\"VALUEN\"}. The options you can configure vary according to which parser you use:  defaultOperator: The default operator used to combine individual terms in the search string. For example: defaultOperator: 'or'. For the dismax parser, you specify a percentage that represents the percentage of terms in the search string (rounded down) that must match, rather than a default operator. A value of 0% is the equivalent to OR, and a value of 100% is equivalent to AND. The percentage must be specified as a value in the range 0-100 followed by the percent (%) symbol. For example, defaultOperator: 50%. Valid values: and, or, a percentage in the range 0%-100% (dismax). Default: and (simple, structured, lucene) or 100 (dismax). Valid for: simple, structured, lucene, and dismax. fields: An array of the fields to search when no fields are specified in a search. If no fields are specified in a search and this option is not specified, all text and text-array fields are searched. You can specify a weight for each field to control the relative importance of each field when Amazon CloudSearch calculates relevance scores. To specify a field weight, append a caret (^) symbol and the weight to the field name. For example, to boost the importance of the title field over the description field you could specify: \"fields\":[\"title^5\",\"description\"]. Valid values: The name of any configured field and an optional numeric value greater than zero. Default: All text and text-array fields. Valid for: simple, structured, lucene, and dismax. operators: An array of the operators or special characters you want to disable for the simple query parser. If you disable the and, or, or not operators, the corresponding operators (+, |, -) have no special meaning and are dropped from the search string. Similarly, disabling prefix disables the wildcard operator (*) and disabling phrase disables the ability to search for phrases by enclosing phrases in double quotes. Disabling precedence disables the ability to control order of precedence using parentheses. Disabling near disables the ability to use the ~ operator to perform a sloppy phrase search. Disabling the fuzzy operator disables the ability to use the ~ operator to perform a fuzzy search. escape disables the ability to use a backslash () to escape special characters within the search string. Disabling whitespace is an advanced option that prevents the parser from tokenizing on whitespace, which can be useful for Vietnamese. (It prevents Vietnamese words from being split incorrectly.) For example, you could disable all operators other than the phrase operator to support just simple term and phrase queries: \"operators\":[\"and\",\"not\",\"or\", \"prefix\"]. Valid values: and, escape, fuzzy, near, not, or, phrase, precedence, prefix, whitespace. Default: All operators and special characters are enabled. Valid for: simple. phraseFields: An array of the text or text-array fields you want to use for phrase searches. When the terms in the search string appear in close proximity within a field, the field scores higher. You can specify a weight for each field to boost that score. The phraseSlop option controls how much the matches can deviate from the search string and still be boosted. To specify a field weight, append a caret (^) symbol and the weight to the field name. For example, to boost phrase matches in the title field over the abstract field, you could specify: \"phraseFields\":[\"title^3\", \"plot\"] Valid values: The name of any text or text-array field and an optional numeric value greater than zero. Default: No fields. If you don't specify any fields with phraseFields, proximity scoring is disabled even if phraseSlop is specified. Valid for: dismax. phraseSlop: An integer value that specifies how much matches can deviate from the search phrase and still be boosted according to the weights specified in the phraseFields option; for example, phraseSlop: 2. You must also specify phraseFields to enable proximity scoring. Valid values: positive integers. Default: 0. Valid for: dismax. explicitPhraseSlop: An integer value that specifies how much a match can deviate from the search phrase when the phrase is enclosed in double quotes in the search string. (Phrases that exceed this proximity distance are not considered a match.) For example, to specify a slop of three for dismax phrase queries, you would specify \"explicitPhraseSlop\":3. Valid values: positive integers. Default: 0. Valid for: dismax. tieBreaker: When a term in the search string is found in a document's field, a score is calculated for that field based on how common the word is in that field compared to other documents. If the term occurs in multiple fields within a document, by default only the highest scoring field contributes to the document's overall score. You can specify a tieBreaker value to enable the matches in lower-scoring fields to contribute to the document's score. That way, if two documents have the same max field score for a particular term, the score for the document that has matches in more fields will be higher. The formula for calculating the score with a tieBreaker is (max field score) + (tieBreaker) * (sum of the scores for the rest of the matching fields). Set tieBreaker to 0 to disregard all but the highest scoring field (pure max): \"tieBreaker\":0. Set to 1 to sum the scores from all fields (pure sum): \"tieBreaker\":1. Valid values: 0.0 to 1.0. Default: 0.0. Valid for: dismax.
- `q.parser`: Specifies which query parser to use to process the request. If queryParser is not specified, Amazon CloudSearch uses the simple query parser.  Amazon CloudSearch supports four query parsers:   simple: perform simple searches of text and text-array fields. By default, the simple query parser searches all text and text-array fields. You can specify which fields to search by with the queryOptions parameter. If you prefix a search term with a plus sign (+) documents must contain the term to be considered a match. (This is the default, unless you configure the default operator with the queryOptions parameter.) You can use the - (NOT), | (OR), and * (wildcard) operators to exclude particular terms, find results that match any of the specified terms, or search for a prefix. To search for a phrase rather than individual terms, enclose the phrase in double quotes. For more information, see Searching for Text in the Amazon CloudSearch Developer Guide.   structured: perform advanced searches by combining multiple expressions to define the search criteria. You can also search within particular fields, search for values and ranges of values, and use advanced options such as term boosting, matchall, and near. For more information, see Constructing Compound Queries in the Amazon CloudSearch Developer Guide.   lucene: search using the Apache Lucene query parser syntax. For more information, see Apache Lucene Query Parser Syntax.   dismax: search using the simplified subset of the Apache Lucene query parser syntax defined by the DisMax query parser. For more information, see DisMax Query Parser Syntax.
- `return`: Specifies the field and expression values to include in the response. Multiple fields or expressions are specified as a comma-separated list. By default, a search response includes all return enabled fields (_all_fields). To return only the document IDs for the matching documents, specify _no_fields. To retrieve the relevance score calculated for each document, specify _score.
- `size`: Specifies the maximum number of search hits to include in the response.
- `sort`: Specifies the fields or custom expressions to use to sort the search results. Multiple fields or expressions are specified as a comma-separated list. You must specify the sort direction (asc or desc) for each field; for example, year desc,title asc. To use a field to sort results, the field must be sort-enabled in the domain configuration. Array type fields cannot be used for sorting. If no sort parameter is specified, results are sorted by their default relevance scores in descending order: _score desc. You can also sort by document ID (_id asc) and version (_version desc). For more information, see Sorting Results in the Amazon CloudSearch Developer Guide.
- `start`: Specifies the offset of the first search hit you want to return. Note that the result set is zero-based; the first result is at index 0. You can specify either the start or cursor parameter in a request, they are mutually exclusive.  For more information, see Paginating Results in the Amazon CloudSearch Developer Guide.
- `stats`: Specifies one or more fields for which to get statistics information. Each specified field must be facet-enabled in the domain configuration. The fields are specified in JSON using the form: {\"FIELD-A\":{},\"FIELD-B\":{}} There are currently no options supported for statistics.
"""
search(q; aws_config::AWSConfig=global_aws_config()) = cloudsearch_domain("GET", "/2013-01-01/search?format=sdk&pretty=true", Dict{String, Any}("q"=>q); aws_config=aws_config)
search(q, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudsearch_domain("GET", "/2013-01-01/search?format=sdk&pretty=true", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("q"=>q), args)); aws_config=aws_config)

"""
    Suggest()

Retrieves autocomplete suggestions for a partial query string. You can use suggestions enable you to display likely matches before users finish typing. In Amazon CloudSearch, suggestions are based on the contents of a particular text field. When you request suggestions, Amazon CloudSearch finds all of the documents whose values in the suggester field start with the specified query string. The beginning of the field must match the query string to be considered a match.  For more information about configuring suggesters and retrieving suggestions, see Getting Suggestions in the Amazon CloudSearch Developer Guide.  The endpoint for submitting Suggest requests is domain-specific. You submit suggest requests to a domain's search endpoint. To get the search endpoint for your domain, use the Amazon CloudSearch configuration service DescribeDomains action. A domain's endpoints are also displayed on the domain dashboard in the Amazon CloudSearch console.

# Required Parameters
- `q`: Specifies the string for which you want to get suggestions.
- `suggester`: Specifies the name of the suggester to use to find suggested matches.

# Optional Parameters
- `size`: Specifies the maximum number of suggestions to return.
"""
suggest(q, suggester; aws_config::AWSConfig=global_aws_config()) = cloudsearch_domain("GET", "/2013-01-01/suggest?format=sdk&pretty=true", Dict{String, Any}("q"=>q, "suggester"=>suggester); aws_config=aws_config)
suggest(q, suggester, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudsearch_domain("GET", "/2013-01-01/suggest?format=sdk&pretty=true", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("q"=>q, "suggester"=>suggester), args)); aws_config=aws_config)

"""
    UploadDocuments()

Posts a batch of documents to a search domain for indexing. A document batch is a collection of add and delete operations that represent the documents you want to add, update, or delete from your domain. Batches can be described in either JSON or XML. Each item that you want Amazon CloudSearch to return as a search result (such as a product) is represented as a document. Every document has a unique ID and one or more fields that contain the data that you want to search and return in results. Individual documents cannot contain more than 1 MB of data. The entire batch cannot exceed 5 MB. To get the best possible upload performance, group add and delete operations in batches that are close the 5 MB limit. Submitting a large volume of single-document batches can overload a domain's document service.  The endpoint for submitting UploadDocuments requests is domain-specific. To get the document endpoint for your domain, use the Amazon CloudSearch configuration service DescribeDomains action. A domain's endpoints are also displayed on the domain dashboard in the Amazon CloudSearch console.  For more information about formatting your data for Amazon CloudSearch, see Preparing Your Data in the Amazon CloudSearch Developer Guide. For more information about uploading data for indexing, see Uploading Data in the Amazon CloudSearch Developer Guide.

# Required Parameters
- `Content-Type`: The format of the batch you are uploading. Amazon CloudSearch supports two document batch formats:  application/json application/xml
- `documents`: A batch of documents formatted in JSON or HTML.

"""
upload_documents(Content_Type, documents; aws_config::AWSConfig=global_aws_config()) = cloudsearch_domain("POST", "/2013-01-01/documents/batch?format=sdk", Dict{String, Any}("documents"=>documents, "headers"=>Dict{String, Any}("Content-Type"=>Content_Type)); aws_config=aws_config)
upload_documents(Content_Type, documents, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = cloudsearch_domain("POST", "/2013-01-01/documents/batch?format=sdk", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("documents"=>documents, "headers"=>Dict{String, Any}("Content-Type"=>Content_Type)), args)); aws_config=aws_config)


text:
module FastRef

using LazySets, LazySets.Approximations
using Polyhedra, CDDLib

using LinearAlgebra
using Parameters
using Interpolations # only for PiecewiseLinear

import LazySets: dim, HalfSpace # necessary to avoid conflict with Polyhedra

using Requires

using Base.Threads

include("activation.jl")
include("network.jl")
include("problem.jl")
include("util.jl")

function __init__()
  @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" include("flux.jl")
end

export
    Solver,
    Network,
    AbstractActivation,
    #PolytopeComplement,
    #complement,
    # NOTE: not sure if exporting these is a good idea as far as namespace conflicts go:
    # ReLU,
    # Max,
    # Id,
    GeneralAct,
    PiecewiseLinear,
    LinearPieces,
    PiecewiseLinearActivation,
    Problem,
    Result,
    BasicResult,
    CounterExampleResult,
    AdversarialResult,
    ReachabilityResult,
    read_nnet,
    solve,
    forward_network,
    check_inclusion

export solve

include("reachability.jl")

include("maxSens.jl")
include("hullGrid.jl")
include("dimGrid.jl")
include("fastGrid.jl")

include("speGuid.jl")
include("hullTree.jl")
include("dimTree.jl")
include("fastTree.jl")

export MaxSens, HullGrid, DimGrid
export FastGrid
export SpeGuid, HullTree, DimTree
export FastTree

end


text:
ENV["GKSwstype"] = "100"  # set 'GR environment' to 'no output' (for Travis CI)
using Documenter, ReachabilityAnalysis

DocMeta.setdocmeta!(ReachabilityAnalysis, :DocTestSetup,
                   :(using ReachabilityAnalysis); recursive=true)

# generate Literate documentation
include("generate.jl")

makedocs(
    format = Documenter.HTML(prettyurls = haskey(ENV, "GITHUB_ACTIONS"),
                             collapselevel = 1,
                             assets = ["assets/juliareach-light.css"]),
    sitename = "ReachabilityAnalysis.jl",
    doctest = false,
    strict = false,
    pages = [
        "Home" => "index.md",
        "Manual" => Any["Basics" => "man/basics.md",
                        "Linear ODEs" => "man/linear.md",
                        "Nonlinear ODEs" => "man/nonlinear.md",
                        #"Semidiscrete PDEs" => "man/pde.md"
                        "Hybrid systems" => "man/hybrid.md",
                        "Clocked systems" => "man/clocked.md",
                        "Exploiting structure" => "man/structure.md",
                        "Parametric reachability" => "man/parametric.md",
                        "Backwards reachability" => "man/backwards.md",
                        "Benchmarks" => "man/benchmarks.md",
                        "Model library" => "man/library.md",
                        "Frequently Asked Questions" => "man/faq.md"],
        "Examples" => Any["Overview" => "man/examples_overview.md",
                          "Operational amplifier" => "models/OpAmp.md",
                          "Square wave oscillator" => "models/SquareWaveOscillator.md",
                          "Van der Pol oscillator" => "models/VanDerPol.md",
                          "Duffing oscillator" => "models/DuffingOscillator.md",
                          "Transmision line" => "models/TransmissionLine.md",
                          "Laub-Loomis" => "models/LaubLoomis.md",
                          "Building" => "models/Building.md",
                          "Production-Destruction" => "models/ProductionDestruction.md",
                          "Lotka-Volterra" => "models/LotkaVolterra.md",
                          "Brusselator" => "models/Brusselator.md",
                          "Structural Model" => "models/ISS.md",
                          "Lorenz system" => "models/Lorenz.md",
                          "Platoon" => "models/Platoon.md",
                          "Quadrotor" => "models/Quadrotor.md",
                          "SEIR model" => "models/SEIR.md",
                          "Spacecraft" => "models/Spacecraft.md"],
                          # "Electromechanic break" => "man/EMBrake.md"
        "Algorithms" => Any["ASB07" => "lib/algorithms/ASB07.md",
                            "BFFPSV18" => "lib/algorithms/BFFPSV18.md",
                            "BOX" => "lib/algorithms/BOX.md",
                            "CARLIN" => "lib/algorithms/CARLIN.md",
                            "GLGM06" => "lib/algorithms/GLGM06.md",
                            "INT" => "lib/algorithms/INT.md",
                            "LGG09" => "lib/algorithms/LGG09.md",
                            "ORBIT" => "lib/algorithms/ORBIT.md",
                            "QINT" => "lib/algorithms/QINT.md",
                            "TMJets" => "lib/algorithms/TMJets.md",
                            "VREP" => "lib/algorithms/VREP.md"],
        "API Reference" => Any["Reach-sets" => "lib/reachsets.md",
                               "Flowpipes" => "lib/flowpipes.md",
                               "Solutions" => "lib/solutions.md",
                               "Systems" => "lib/systems.md",
                               "Discretization" => "lib/discretize.md",
                               "Projections" => "lib/projections.md",
                               "Clustering" => "lib/clustering.md",
                               "Further set operations" => "lib/operations.md",
                               "Distributed computations" => "lib/distributed.md",
                               "Internal functions and macros" => "lib/internals.md"],
        "References" => "references.md",
        "About" => "about.md"
    ]
)

# Deploy built documentation from Travis.
deploydocs(
    repo = "github.com/JuliaReach/ReachabilityAnalysis.jl.git",
    push_preview = true,
)


text:
using Thebes, Luxor, Random

include(dirname(pathof(Thebes)) * "/../data/moreobjects.jl")

mygfunction(vertices, faces, labels) =
    Thebes.simplegfunction(vertices, faces, labels, action=:fill)

function frame(scene, framenumber)
    background("white")
    # flyaround
    theta = rescale(framenumber, 1, scene.framerange.stop, 0, 2pi)
    eyepoint(Point3D(150 + 80(cos(theta)), 150 + 80(sin(theta)), 40(sin(theta))))

    sethue("white")
    carpet(1000)
    drawcube(250)
    axes3D(250)

    setline(.5)
    setopacity(0.5)
    st = [-2, 1, 2]
    for x in st, y in st, z in st
      n = x * y * z
      object = make(tetrahedron)
      scaleby!(object, 10, 10, 10)
      moveby!(object, x * 25, y * 25, z * 25)
      d = distance(object.vertices[end], Point3D(0, 0, 0))
      if d < 1200
          pin(object, gfunction=mygfunction)
      end
    end
    sethue("black")
    text(string(eyepoint()), Point(0, boxheight(BoundingBox()/2) - 20))
end

width, height = (800, 800)

threedmovie = Movie(width, height, "threedmovie")

animate(threedmovie, [
    Scene(threedmovie, frame, 1:25)
    ],
    creategif=true,
    framerate=20,
    pathname="/tmp/threed-animation.gif")


text:
using LinearAlgebra
function Sparse_matrix_SINDYPI(Theta, dx, lambda, nvars, maxiter)

normLib = zeros(1, size(Theta,2))

Thetan = deepcopy(Theta)

    for norm_k=1:size(Theta,2)
            normLib[norm_k] = norm(Thetan[:,norm_k]);
            Thetan[:,norm_k] = Thetan[:,norm_k]/normLib[norm_k];
    end

   #Epsilon = Thetan\dx

#Below:Only to be used for ill conditioned matrices, like in the glycolysis example. Check on this
Epsilon = zeros(size(Thetan,2),1)
ldiv!(Epsilon, qr(Thetan), dx)

#=
   Epsilon_c = similar(Epsilon)
   Epsilon_c .= Epsilon
=#

    for i in 1:maxiter
        index_remove = abs.(Epsilon) .<= lambda
        Epsilon[index_remove] .= zero(eltype(Epsilon))

        for j in 1:nvars
            index_retain = @. !index_remove[:,j]
            Epsilon[index_retain, j] = Thetan[:, index_retain] \ dx[:,j]
        end
   #=
    if norm(Epsilon_c - Epsilon, 2) < eps()
        break
    else
        Epsilon_c .= Epsilon
    end
=#
    end

        for norm_k=1:length(Epsilon)
            Epsilon[norm_k,:] = Epsilon[norm_k,:]/normLib[norm_k];
        end

    return Epsilon
end

#=
A. Q factor

qrf.R

Ans = inv(qrf.R)*(qrf.Q)'*dx
=#


text:
##### Beginning of file

Test.@test( Base.VERSION >= VersionNumber("1.0") )

Test.@test( PredictMDFull.version() > VersionNumber(0) )

Test.@test(
    PredictMDFull.version() ==
        PredictMDFull.version(PredictMDFull)
    )

Test.@test(
    PredictMDFull.version() ==
        PredictMDFull.version(first(methods(PredictMDFull.eval)))
    )

Test.@test(
    PredictMDFull.version() ==
        PredictMDFull.version(PredictMDFull.eval)
    )

Test.@test(
    PredictMDFull.version() ==
        PredictMDFull.version(PredictMDFull.eval, (Any,))
    )

Test.@test( PredictMDFull.version(TestModuleA) == VersionNumber("1.2.3") )

Test.@test( PredictMDFull.version(TestModuleB) == VersionNumber("4.5.6") )

Test.@test_throws(
    ErrorException,
    PredictMDFull._TomlFile(joinpath(mktempdir(),"1","2","3","4")),
    )

##### End of file


text:
module OpenTelemetryInstrumentationBase

using OpenTelemetryAPI
import TOML

const PKG_VERSION =
    VersionNumber(TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"])

const INSTRUMENTATION_INFO = InstrumentationInfo(
    name = string(@__MODULE__),
    version = PKG_VERSION,
    schema_url = "https://oolong.dev/OpenTelemetry.jl/dev/OpenTelemetryInstrumentationBase/",
)

"""
## Schema

| Meter Name | Instrument Name     | Instrument Type   | Unit        | Description                                                                                |
|:---------- |:------------------- |:----------------- |:----------- |:---------------------------------------------------------------------------------------------------------------- |
| `Base.Sys` | `uptime`            | `ObservableGauge` | second      | The current system uptime in seconds.                                                                            |
| `Base.Sys` | `free_memory`       | `ObservableGauge` | bytes       | The total free memory in RAM in bytes.                                                                           |
| `Base.Sys` | `free_memory_ratio` | `ObservableGauge` | %           | The ratio of free memory in percentage.                                                                          |
| `Base.Sys` | `maxrss`            | `ObservableGauge` | bytes       | The maximum resident set size utilized in bytes.The current system uptime in seconds.                            |
| `Base`     | `jit_total_bytes`   | `ObservableGauge` | bytes       | The total amount (in bytes) allocated by the just-in-time compiler.                                              |
| `Base`     | `gc_live_bytes`     | `ObservableGauge` | bytes       | The total size of live objects after the last garbage collection, plus the number of bytes allocated since then. |
| `Base`     | `gc_time_ns`        | `ObservableGauge` | nanoseconds | Total time spend in garbage collection, in nanoseconds                                                           |
"""
function init(;
    meter_provider = global_meter_provider(),
    tracer_provider = global_tracer_provider(),
)
    m_sys = Meter(
        "Base.Sys";
        provider = meter_provider,
        instrumentation_info = INSTRUMENTATION_INFO,
    )

    ObservableGauge{Float64}(
        Sys.uptime,
        "uptime",
        m_sys;
        unit = "second",
        description = "The current system uptime in seconds.",
    )

    ObservableGauge{UInt64}(
        Sys.free_memory,
        "free_memory",
        m_sys;
        unit = "bytes",
        description = "The total free memory in RAM in bytes.",
    )

    ObservableGauge{Float64}(
        () -> Sys.free_memory() / Sys.total_memory() * 100,
        "free_memory_ratio",
        m_sys;
        unit = "%",
        description = "The ratio of free memory in percentage.",
    )

    ObservableGauge{UInt64}(
        Sys.maxrss,
        "maxrss",
        m_sys;
        unit = "bytes",
        description = "The maximum resident set size utilized in bytes.",
    )

    m = Meter(
        "Base";
        provider = meter_provider,
        instrumentation_info = INSTRUMENTATION_INFO,
    )

    ObservableGauge{Int64}(
        Base.jit_total_bytes,
        "jit_total_bytes",
        m;
        unit = "bytes",
        description = "The total amount (in bytes) allocated by the just-in-time compiler.",
    )

    ObservableGauge{Int64}(
        Base.gc_live_bytes,
        "gc_live_bytes",
        m;
        unit = "bytes",
        description = "The total size of live objects after the last garbage collection, plus the number of bytes allocated since then",
    )

    ObservableGauge{UInt64}(
        Base.gc_time_ns,
        "gc_time_ns",
        m;
        unit = "nanoseconds",
        description = "Total time spend in garbage collection, in nanoseconds",
    )
end

function __init__()
    init()
end

end # module


text:
# Co-location maps classificator
struct CoLocation <:  AbstractCoLocationMapClassificator
    mcr::AbstractMapClassificator
    CoLocation(mcr::AbstractMapClassificator) = return new(mcr);
end

# Structure for Co-location  Map Classification
struct CoLocationMapClassification <: AbstractMapClassification
    mcr::AbstractMapClassificator
    k::Int
    group::Vector{Int}
    ngroup::Vector{Int}
    grouplabs::Vector{String}
end

function mapclassify(mcr::AbstractCoLocationMapClassificator, x::Vector{Vector{T}} where T;
    lower::Symbol = :open, upper::Symbol = :closed)::CoLocationMapClassification

    l = length(x)
    n = length(x[1])
    lgroup = zeros(Int, n, l)
    lk = zeros(Int, l)
    comcr = mcr.mcr

    if isa(comcr, AbstractGraduatedMapClassificator)

        # Get categories for each variable
        for i in 1:l
            ni = length(x[i])
            (n == ni) || throw(DimensionMismatch("dimensions must match: x[1] has ($(n)), x[$(i)] has ($ni)"))

            mc = mapclassify(comcr, x[i], lower = lower, upper = upper)
            lgroup[:,i] = assignments(mc)
            lk[i] = length(mc)

            (lk[1] == lk[i]) || throw(DimensionMismatch("classifications categories must match: 1 has ($(lk[1])), $(i) has ($(lk[i]))"))
        end

        k = lk[1] + 1

        # Get common categories
        group = zeros(Int, n)
        ngroup = zeros(Int, k)

        for i in 1:n
            gi = lgroup[i,1]
            if all(lgroup[i,:] .== gi)
                group[i] = gi + 1
                ngroup[gi + 1] += 1
            else
                group[i] = 1
                ngroup[1] += 1
            end
        end

        grouplabs = vcat("Other", string.(1:k-1))
    elseif isa(comcr, AbstractUniqueMapClassificator)
        # Get all categories
        cats = unique(mapreduce(unique, vcat, x))
        sort!(cats)
        k = length(cats) + 1
        group = zeros(Int, n)
        ngroup = zeros(Int, k)
        grouplabs = vcat("Other", string.(cats))

        for i in 1:n
            xi = x[1][i]
            for j in 2:l
                # Count number of coincidences
                c = 1
                if x[j][i] == xi
                    c += 1
                end

                if c == l
                    g = findall(y -> y == xi, cats)[1]
                    group[i] = g + 1
                    ngroup[g + 1] += 1
                else
                    group[i] = 1
                    ngroup[1] += 1
                end
            end
        end

    end

    return CoLocationMapClassification(mcr, k, group, ngroup, grouplabs)
end

function maplabels(mc::CoLocationMapClassification; counts::Bool = true)::Vector{String}

    ngroup = mc.ngroup
    grouplabs = mc.grouplabs
    labels = grouplabs
    labels[1] = "Other"

    if counts
        labels = labels .* " (" .* string.(ngroup) .* ")"
    end

    labels
end

levels(mc::CoLocationMapClassification)::Vector{String} = return mc.grouplabs;


text:
using NeuralQuantum, Test
using NeuralQuantum: set_index!, trainable_first
using NeuralQuantum: out_similar, unsafe_get_batch
num_types = [Float32, Float64]
atol_types  = [1e-5, 1e-8]

machines = Dict()

ma = (T, N) -> RBMSplit(T, N, 2)
machines["RBMSplit"] = ma

ma = (T, N) -> NDM(T, N, 2, 3, NeuralQuantum.logℒ)
machines["NDM_softplus"] = ma

ma = (T, N) -> NDM(T, N, 2, 3, NeuralQuantum.logℒ2)
machines["NDM_cosh"] = ma

ma = (T, hilb) -> begin
    N = nsites(hilb)
    NDMSymm(T, N, 1, 2, translational_symm_table(HyperCube([N], periodic=true)), NeuralQuantum.logℒ2)
end
machines["NDMSymm_cosh"] = ma

ma = (T, N) -> RBM(T, N, 2, NeuralQuantum.logℒ)
machines["RBM_softplus"] = ma

ma = (T, N) -> RBM(T, N, 2, NeuralQuantum.logℒ2)
machines["RBM_cosh"] = ma

ma = (T, hi) -> begin
    N = nsites(hi)
    ch = Chain(Dense(T, N, N-1, af_softplus), Dense(T, N-1, N-2, af_softplus), WSum(T, N-2))
    return PureStateAnsatz(ch, hi)
end
machines["chain_pure"] = ma

ma = (T, hi) -> begin
    N = nsites(hi)
    ch = Chain(Dense(T, N, 2*N, af_softplus), sum_autobatch)
    return PureStateAnsatz(ch, N)
end
machines["chain_pure_softplus"] = ma

ma = (T, hi) -> begin
    N = nsites(hi)
    ch = Chain(Dense(T, N, 2*N, af_logcosh), sum_autobatch)
    return PureStateAnsatz(ch, N)
end
machines["chain_pure_cosh"] = ma

ma = (T, hi) -> begin
    N = nsites(hi)
    ch = Chain(DenseSplit(T, N, 2*N, af_softplus), sum_autobatch)
    return MixedStateAnsatz(ch, N)
end
machines["chain_mixed_softplus"] = ma

N = 4
T = Float32
b_sz = 3

@testset "test batched dispatch - values: $name" for name=keys(machines)
    name == "NDMSymm_cosh" && continue
    hilb = HomogeneousFock(N, 2)
    net = machines[name](T,hilb)
    if net isa NeuralQuantum.MatrixNet
        hilb = SuperOpSpace(hilb)
    end
    cnet = cached(net, b_sz)
    v  = state(T, hilb, net)
    vb = state(T, hilb, net, b_sz)

    rand!(vb, hilb)

    @test net(vb) ≈ cnet(vb)
    if (vb isa Tuple)
        @test net(vb...) ≈ cnet(vb...)
    end

    o = rand(Complex{T}, 1, b_sz)
    o2 = similar(o)

    oo  = logψ!(o, net, vb)
    oo2 = logψ!(o2, cnet, vb)
    @test oo === o
    @test oo2 === o2
    @test oo ≈ oo2
end

@testset "test cached dispatch - inplace gradients: $name" for name=keys(machines)
    hilb = HomogeneousFock(N, 2)
    net = machines[name](T,hilb)

    snet = cache(net)
    cnet = cached(net, b_sz)
    if net isa NeuralQuantum.MatrixNet
        hilb = SuperOpSpace(hilb)
    end

    v  = state(T, hilb, net)
    vb = state(T, hilb, net, b_sz)

    rand!(vb, hilb)

    g1 = grad_cache(net, b_sz)
    g2 = grad_cache(net, b_sz)
    g3 = grad_cache(net, b_sz)

    v1, gg1 = logψ_and_∇logψ!(g1, cnet, vb)
    v3 = similar(v1)
    vv3, gg3 = logψ_and_∇logψ!(g3, v3, cnet, vb)
    @test v1 ≈ v3
    @test gg1 ≈ gg3
    @test g3 === gg3
    @test v3 === vv3

    if vb isa Tuple
        v2, gg2 = logψ_and_∇logψ!(g2, cnet, vb...)
        @test v1 ≈ v2
        @test gg1 ≈ gg2
        @test g1 === gg1 && g2 === gg2

        vv3, gg3 = logψ_and_∇logψ!(g3, v3, cnet, vb...)
        @test v1 ≈ v3
        @test gg1 ≈ gg3
        @test g3 === gg3
        @test v3 === vv3
    end
end

@testset "test batched logψ: $name" for name=keys(machines)
    hilb = HomogeneousFock(N, 2)

    net = machines[name](T,hilb)

    cnet = cached(net)
    bnet = cached(net, b_sz)
    if net isa NeuralQuantum.MatrixNet
        hilb = SuperOpSpace(hilb)
    end

    v  = state(T, hilb, net)
    vb = state(T, hilb, net, b_sz)

    vals   = out_similar(bnet)
    vals_2 = out_similar(bnet)
    res    = Bool[]
    for batch_i=Iterators.partition(1:spacedimension(hilb), b_sz)
        for (i, h_i)=enumerate(batch_i)
            set!(unsafe_get_batch(vb, i), hilb, h_i)
        end
        logψ!(vals, bnet, vb)

        for (i, h_i)=enumerate(batch_i)
            set!(v, hilb, h_i)
            vals_2[i] = logψ(cnet, v)
        end
        push!(res, vals_2 ≈ vals)
    end
    @test all(res)

end

@testset "test batched ∇ψ: $name" for name=keys(machines)
    name == "NDMSymm_cosh" && continue

    hilb = HomogeneousFock(N, 2)

    net = machines[name](T,hilb)

    cnet = cached(net)
    bnet = cached(net, b_sz)
    if net isa NeuralQuantum.MatrixNet
        hilb = SuperOpSpace(hilb)
    end

    v  = state(T, hilb, net)
    vb = state(T, hilb, net, b_sz)

    vals   = out_similar(bnet)
    vals_2 = out_similar(bnet)

    gb = grad_cache(net, b_sz)
    g2 = grad_cache(net)
    gtmp = grad_cache(net)
    res    = Bool[]
    res_grad = Bool[]
    for batch_i=Iterators.partition(1:spacedimension(hilb), b_sz)
        for (i, h_i)=enumerate(batch_i)
            set!(unsafe_get_batch(vb, i), hilb, h_i)
        end
        logψ_and_∇logψ!(gb, vals, bnet, vb)

        for (i, h_i)=enumerate(batch_i)
            set!(v, hilb, h_i)
            vals_2[i] , _ = logψ_and_∇logψ!(g2, cnet, v)

            for (x,y)=zip(vec_data(gtmp), vec_data(gb))
                x .= y[:,i]
            end
            push!(res_grad, gtmp ≈ g2)
        end
        push!(res, vals_2 ≈ vals)
    end
    @test all(res)
    @test all(res_grad)
end


text:
module GinzburgLandau

# package code goes here

dct_nodes(N, L) = (0:N-1)*L/N + L/(2N)

AdamsMoulton = [1          0          0         0         0          0;
                1/2        1/2        0         0         0          0;
                5/12       8/12      -1/12      0         0          0;
                9/24       19/24     -5/24      1/24      0          0;
                251/720    646/720   -264/720   106/720  -19/720     0;
                475/1440   1427/1440 -798/1440  482/1440 -173/1440   27/1440]

AdamsBashforth=[1          0          0          0        0          0;
                3/2        -1/2       0          0        0          0;
                23/12      -16/12     5/12       0        0          0;
                55/24      -59/24     37/24     -9/24     0          0;
                1901/720   -2774/720  2616/720  -1274/720  251/720   0;
                4277/1440  -7923/1440 9982/1440 -7298/1440 2877/1440 -475/1440]

type GLECoeffs
    σ::Float64
    γ::Float64
    c₁::Float64
    c₃::Float64
    k::Vector{Complex{Float64}}
end

function GLECoeffs(σ, γc, c₁, c₃)
    k = zeros(Complex{Float64}, 3)
    γ = γc/c₁
    k[1] = σ
    k[2] = γ*(1 + im*c₁)
    k[3] = -σ*(1 - im*c₃)
    GLECoeffs(σ, γ, c₁, c₃, k)
end

import Base.getindex

getindex(coeffs::GLECoeffs, i) = coeffs.k[i]

type GLESolver1
    L::Float64
    N::Int
    Δt::Float64
    δ::Vector{Complex{Float64}}
    g::Vector{Complex{Float64}}
end

function GLESolver1(L, N, Δt, k::GLECoeffs)

    δ = zeros(Complex{Float64}, N)

    for i = 1:N
        w = π/L*(i-1)
        δ[i] = 1.0 - Δt*(k[1] - k[2]*w^2)
    end

    g = zeros(Complex{Float64}, N)
    GLESolver1(L, N, Δt, δ, g)
end

nmodes(slv::GLESolver1) = slv.N
nodes(slv::GLESolver1) = dct_nodes(slv.N, slv.L)
tstep(slv::GLESolver1, i) = (i-1)*slv.Δt

function gle_step!(slv::GLESolver1, k::GLECoeffs, a₀::AbstractVector{Complex{Float64}},
                   a₁::AbstractVector{Complex{Float64}})
    # Calcularo f
    g = slv.g
    dt = slv.Δt
    N = slv.N
    δ = slv.δ

    copy!(a₁, a₀)
    idct!(a₁)

    for i = 1:N
        g[i] = k[3] * abs(a₁[i])^2 * a₁[i]
    end

    dct!(g)
    for i = 1:N
        a₁[i] = ( a₀[i] + dt*g[i] ) / δ[i]
    end

    idct!(a₁)
    for i = 1:N
        g[i] = k[3] * abs(a₁[i])^2 * a₁[i]
    end
    dct!(g)

    for i = 1:N
        a₁[i] = ( a₀[i] + dt*g[i] ) / δ[i]
    end

    return nothing

end

type GLEAdamsBashforth
    L::Float64
    N::Int
    s::Int
    Δt::Float64
    f::Vector{Vector{Complex{Float64}}
    a₀::Vector{Complex{Float64}}
    idx::Vector{Integer}
    step::Int
    function GLEAdamsBashforth(L, N, s, Δt)
        f = Vector{Vector{Complex{Float64}}}(s)
        for i = 1:s
            f[i] = zeros(Complex{Float64}, N)
        end
        a₀ = zeros(Complex{Float64}, N)
        new(L, N, s, Δt, f, a₀, zeros(Int, s), 0)
    end

end
nmodes(slv::GLEAdamsBashforth) = slv.N
nodes(slv::GLEAdamsBashforth) = dct_nodes(slv.N, slv.L)
tstep(slv::GLEAdamsBashforth, i) = (i-1)*slv.Δt
nstep(slv::GLEAdamsBashforth) = slv.s

function init_gle!(slv::GLEAdamsBashforth, a::AbstractVector)
    s = nstep(slv)
    slv.step = 0
    fill!(slv.idx, 0)
    for i = 1:s
        fill!(slv.f[i], 0)
    end
    fill!(slv.a₀, a)

end

function simple_step!(slv::GLEAdamsBashforth, a::AbstractVector{Complex{Float64}})

end # module


text:
# Count the OpenMP changes to the BATS-R-US source code.
#
# Hongyang Zhou, hyzhou@umich.edu 05/09/2019

using DataFrames, Printf, Query

function extractOpenmpInfo(File::AbstractString; Dir::AbstractString=".",
   Debug::Bool=false)

nLine,nOpenmpLine = 0,0 # Initialize counters
for line in eachline(joinpath(Dir,File))
   if occursin(r"^\s*!\$omp",line) nOpenmpLine += 1 end
   # Exclude comment/empty lines
   if !(occursin(r"^\s*![^$]",line) || occursin(r"^\s*$",line))
      nLine += 1
   end
end

if Debug
   println("Searching in $File...")
   println("Total lines:",nLine)
   println("Total OpenMP lines:",nOpenmpLine)
end

df = DataFrame(Name=File, TotalLine=nLine, OpenmpLine=nOpenmpLine)

return df

end

function postprocess(DirIn::AbstractString=".")

dfFiles = DataFrame(Name=String[], TotalLine=Int32[], OpenmpLine=Int32[])

# Extract run time parameters and timings from runlog files
for (root, dirs, files) in walkdir(DirIn)
   !occursin(r"\/src\w*$",root) && continue

   println("Files in $root")
   for file in filter(x -> endswith(x,"f90"),files)
      df = extractOpenmpInfo(file,Dir=root)
      try
         append!(dfFiles,df)
      catch
         println("No info from file ",file)
      end
   end
end

return dfFiles

end

dfFiles = postprocess("/Users/hyzhou/SWMF/test/BATSRUS/srcBATL")
display(describe(dfFiles))
#showall(dfFiles)
ratio = sum(dfFiles[:OpenmpLine]) / sum(dfFiles[:TotalLine])
@printf("OpenMP ratio = %4.2f%%",ratio*100)


text:
module GaussianInterpolation

using CMBLensing
using FFTW
using Interpolations
using LinearAlgebra
using Match
using NearestNeighbors: KDTree, knn
using ProgressMeter
using Statistics
using StaticArrays
using SparseArrays
import Statistics: cov
import Base: |, eltype, getindex, +, -, *, length, similar, parent, iterate, show, eachindex
using LinearAlgebra
using PositiveFactorizations

export @unpack
export coords

axes_centered(A) = Tuple((ax .- offset for (ax, offset) in zip(axes(A), size(A) .÷ 2 .+ 1)))

function coords(m::AbstractMatrix)
    I = CartesianIndices(m)
    X = Matrix{Int32}(undef, length(I), 2)
    for (i, ind) in zip(LinearIndices(I), I)
        X[i, :] .= Tuple(ind)
    end
    X
end
coords(f::FlatMap) = coords(f.Ix)
coords(f::FlatFourier) = coords(Map(f).Ix)
coords(f::FieldTuple) = coords(first(f))

# Functions for filling covariance matrices / vectors
function cov!(out::AbstractVector, Kitp, X)
    for i in 1:length(out)
        out[i] = Kitp(X[i, 1], X[i, 2])
    end
    out
end
function cov!(out::AbstractMatrix, Kitp, X, Y=X)
    @assert size(out) == @. first(size((X, Y)))
    ni, nj = size(out)
    for j in 1:nj, i in 1:ni
        out[i,j] = Kitp(X[i, 1] - Y[j, 1], X[i, 2] - Y[j, 2])
    end
    out
end
function cov!(out::Symmetric, Kitp, X)
    n = first(size(X))
    @assert all(n .== size(out))
    for j in 1:n, i in 1:j,
        out.data[i,j] = Kitp(X[i, 1] - X[j, 1], X[i, 2] - X[j, 2])
    end
    out
end

# Kernels and Masks
export make_kernel, make_mask
make_kernel(L::LinOp) = make_kernel(diag(L))
make_kernel(f::FlatFourier) = make_kernel(fftshift(Map(f).Ix))
make_kernel(m::Matrix) = CubicSplineInterpolation(axes_centered(m), m, extrapolation_bc=NaN)
# make_kernel(m::Matrix) = LinearInterpolation(axes_centered(m), m, extrapolation_bc=0.0)

function parse_predict(predict::Tuple)
	Tuple((p isa Symbol ? p => p : p) for p in predict)
end
function kernel_dict(C::LinOp, predict::NTuple)
    names = (Symbol(p...) for p in predict)
    Dict(name => make_kernel(C[name]) for name in names)
end

make_mask(L::LazyBinaryOp)= make_mask(L.b)
make_mask(D::DiagOp{<:Field{Map}})= diag(D)
make_mask(D::DiagOp{<:FlatS02})= diag(D)[:I]

export sparse_gp
function sparse_gp(ds::DataSet; m=15)
	X = coords(ds.d)
	Xtree = KDTree(permutedims(float.(X)))
	skip_predicate = (i -> iszero(mask[i]) ? true : false)
    design_inds(i) = knn(Xtree, X[i, :], m, false, skip_predicate)[1]

	Kf, Kn = make_kernel.((ds.Cf, ds.Cn))
	mask = make_mask(ds.M)

    function compute_row!(i, J, V)
        js = design_inds(i)
        Xj = X[js, :] .- X[i:i, :]

        s = zeros(T, m)
        S = zeros(T, (m,m))
        N = similar(S)
        cov!(s, Kf, Xj)
        cov!(S, Kf, Xj)
        cov!(N, Kn, Xj)

        inds = (i-1)*m+1:i*m
        J[inds] = js
        V[inds] = (S + N) \ s
    end

    @unpack Nside, T = fieldinfo(ds.d)
    I = Vector{Int32}((m:m*Nside^2 + (m-1)) .÷ m)
    J = similar(I)
    V = similar(I,T)
    for i in 1:Nside^2
        compute_row!(i, J, V)
    end
    SparseMatGP(sparse(I,J,V,Nside^2,Nside^2))
end

mutable struct SparseMatGP
    M::SparseMatrixCSC
end
*(gp::SparseMatGP, f::FlatMap{P}) where P = FlatMap{P}(reshape(gp.M * f[:], size(f.Ix)))
*(gp::SparseMatGP, f::FlatFourier) = gp * Map(f)

export local_gp
function local_gp(ds::DataSet, predict::NTuple, m::Int; ϕ=nothing)
    @unpack Nside, T, Δx = fieldinfo(ds.d)

    # Kernels, Covariance Matrices, and Mask
    predict = parse_predict(predict)
    Ks, Kn = kernel_dict.((ds.Cf, ds.Cn), Ref(predict))
    s, S, N = zeros(T, m), zeros(T, (m,m)), zeros(T, (m,m))

    # Coordinates and Nearest Nearest Neighbors Tree
    X_predict = coords(ds.d)
    if ϕ !== nothing
        αy, αx = getindex.(Map.(∇ * ϕ ./ Δx ), :)
        X_obs = X_predict .+ [αx αy]
    else
        X_obs = X_predict
    end
    Xtree = KDTree(permutedims(float.(X_obs)))

    # Local GP Design
    if ds.M == 1
        skip_predicate = x -> false
    else
        mask = make_mask(ds.M)
        skip_predicate = (i -> iszero(mask[i]) ? true : false)
    end
    design_inds(i) = knn(Xtree, X_obs[i, :], m, true, skip_predicate)[1]

    function apply(f::FlatField{P, T, M}) where {P, T, M}
        f_in, f_out = Dict{Symbol, Matrix{T}}(), Dict{Symbol, Matrix{T}}()
        for (p_in, p_out) in predict
            f_in[p_in] = f[Symbol(string(p_in)*"x")]
            f_out[p_out] = similar(f_in[p_in])
        end
        for i in 1:Nside^2
            # calculate reference point coordinates
            js = design_inds(i)
            Xj = X_obs[js, :] .- X_predict[i:i, :]

            # loop over things we're predicting
            for (in, out) in predict
                for (m, kernel) in ((s, Ks), (S, Ks), (N, Kn))
                    cov!(m, kernel[Symbol(in, out)], Xj)
                end
                f_out[out][i] = s ⋅ ((S + N) \ f_in[in][js])
            end
        end
        (; (p_out => FlatMap{P,T}(f_out[p_out]) for p_out in last.(predict))...)
    end
    apply
end

export full_gp
function full_gp(ds::DataSet, predict::NTuple; ϕ=nothing)
    @unpack Nside, T, Δx = fieldinfo(ds.d)
    Npix = Nside^2

    # Kernels, Covariance Matrices, and Mask
    predict = parse_predict(predict)
    Ks, Kn = kernel_dict.((ds.Cf, ds.Cn), Ref(predict))
    S, N = zeros(T, (Npix, Npix)), zeros(T, (Npix, Npix))

    # Coordinates and Nearest Nearest Neighbors Tree
    X_predict = coords(ds.d)
    if ϕ !== nothing
        αy, αx = getindex.(Map.(∇ * ϕ ./ Δx ), :)
        X_obs = X_predict .+ [αx αy]
    else
        X_obs = X_predict
    end

    function apply(f::FlatField{P, T, M}) where {P, T, M}
        f_in, f_out = Dict{Symbol, Matrix{T}}(), Dict{Symbol, Matrix{T}}()
        for (p_in, p_out) in predict
            f_in[p_in] = f[Symbol(string(p_in)*"x")]
            f_out[p_out] = similar(f_in[p_in])
        end

        # loop over things we're predicting
        for (in, out) in predict
            cov!(S, Ks[Symbol(in, out)], X_obs)
            cov!(N, Kn[Symbol(in, out)], X_obs)
            vec(f_out[out]) .= S * (cholesky(Positive, S + N) / f_in[in][:])
        end

        (; (p_out => FlatMap{P,T}(f_out[p_out]) for p_out in last.(predict))...), S, N
    end
    apply
end

# module end
end


text:
"Returns the time-scale for advection on a regular grid across a single grid cell."
function cell_advection_timescale(u, v, w, grid)
    umax = maximum(abs, u)
    vmax = maximum(abs, v)
    wmax = maximum(abs, w)

    Δx = grid.Δx
    Δy = grid.Δy
    Δz = grid.Δz

    return min(Δx/umax, Δy/vmax, Δz/wmax)
end

cell_advection_timescale(model) =
    cell_advection_timescale(model.velocities.u.data.parent,
                             model.velocities.v.data.parent,
                             model.velocities.w.data.parent,
                             model.grid)


text:
export
    scale2log,
    normlog,
    sph2cart,
    cart2sph,
    mod_angle,
    mod_angle_deg,
    mod_angle_rad

# operation
scale2log(x) = 20log10(x)
normlog(vec) = vec .- maximum(vec)
sph2cart(θ, ϕ, r) = begin
    (
        x=r * sin(θ)cos(ϕ),
        y=r * sin(θ)sin(ϕ),
        z=r * cos(θ)
    )
end
cart2sph(x, y, z) = begin
    # XXX r are not calcuated, may cause problem
    r = sqrt(x^2 + y^2 + z^2)
    θ = acos(z / r)
    ϕ = atan(y, x)
    (
        θ=θ,
        ϕ=ϕ,
        r=r
    )
end
# convert to domain θ ∈[0°, 180°], ϕ∈ [-180°, 180°]
# bad implementation
function mod_angle(θ, ϕ)
    sign = x -> x >= 0 ? 1 : -1
    item1 = -sign.(mod2pi.(θ) .- π)
    θ_ = mod1.(θ .* item1, π)
    θ_[θ.==0] .= 0
    ϕ_ = item1 .* ϕ .|> mod2pi
    (θ_, ϕ_)
end
function mod_angle(θ::Number, ϕ)
    sign = x -> x >= 0 ? 1 : -1
    item1 = -sign(mod2pi(θ) - π)
    θ_ = (θ==0 ?  0.0 : mod1(θ * item1, π))
    ϕ_ = (item1 .* ϕ) .|> mod2pi
    (θ_, ϕ_)
end

# convert to domain θ ∈[0°, 180°], ϕ∈ [-180°, 180°]
function mod_angle_deg(θ::Number, ϕ)
    θ_res = 0.
    ϕ_res = 0.
    invert_ϕ = false
    # in domain
    if θ >= 0 && θ <= 180
        θ_res = θ
    elseif θ < 0 && θ >= -180
        invert_ϕ = true
        θ_res = -θ
    else
        throw("θ = $(θ) is out of range [-180,180]")
    end

    if ϕ >=-180 && ϕ <= 180
        ϕ_res = ϕ
    elseif ϕ > 180 && ϕ <=360
        ϕ_res = ϕ - 360
    else
        throw("ϕ = $(ϕ) is out of range [-180,360]")
    end

    if invert_ϕ
        if ϕ_res <= 0
            ϕ_res += 180
        else
            ϕ_res -= 180
        end
    end
    θ_res, ϕ_res
end

# not test
function mod_angle_rad(θ::Number, ϕ)
    θ_res = 0.
    ϕ_res = 0.
    invert_ϕ = false
    # in domain
    if θ >= 0 && θ <= pi
        θ_res = θ
    elseif θ < 0 && θ >= -pi
        invert_ϕ = true
        θ_res = -θ
    else
        throw("θ = $(θ) is out of range [-pi, pi]")
    end

    if ϕ >=-pi && ϕ <= pi
        ϕ_res = ϕ
    elseif ϕ > pi && ϕ <=2pi
        ϕ_res = ϕ - 2pi
    else
        throw("ϕ = $(ϕ) is out of range [-pi, 2pi]")
    end

    if invert_ϕ
        if ϕ_res <= 0
            ϕ_res += pi
        else
            ϕ_res -= pi
        end
    end
    θ_res, ϕ_res
end