text:
function basicThomasian(givenT, givenJ, givenF, givenK, givenRegions, givenD, givenS0, givenS, givenG)
	#%realistic is the implementation of the "On a more realistic lock
	#%contention model and its analysis" paper from ICDE 1994
	T_total = givenT; #% the requested TPS
	J = givenJ; #%  of transaction classes, e.g. in TPC-C, J=5
	f = givenF; #% frequency matrix for each class
	K = givenK; #% the number of locks requested in each transaction class
	regions = givenRegions; #% number of database regions, perhaps should be equal to the number of tables?!
	D = givenD; #% number of data items in the i'th DB region
	maxK = Max(K)
	S = givenS; #% S_jn is the processing of the n'th step of a transaction of class C_j
	S0 = givenS0; #% Barzan assumes that everything that precedes the first lock is constant time for all transactions
	g = givenG; #%g_jni is the probability that a tran type j access the i'the database region in its n'th step
	cores = Inf;

	R = 0;
	M_total = 0;
	Vp = 0;
	V = 0;
	Pcon = 0;
	L = 0;

	if length(f)!=J || length(K)!=J || length(D)!= regions || size(S,1) != J || size(S,2) != maxK
		error("Invalid input arguments");
	end

	for j=1:J
		for n=1:K[j]
			if abs(sum(g[j,n,:]) - 1.0) > 1e-7
				j
				n
				#%g[j,n,:]
				why=sum(g[j,n,:])
				error("The access pattern g is malformed");
			end
		end
	end

	K_total = dot(vec(f),vec(K));

	D_total = sum(D);

	U = zeros(convert(Int,J),convert(Int,maxK)); #% initialization
	for j=1:J
		U[j,1:K[j]] = 0;
	end

	deltaW=1000; W = 1000;
	iter = 0;

	ClientsNum = 160;

	while abs(deltaW)> 0.001 && !isnan(W) && W!=Inf && iter <=100 #% repetetive process
		iter = iter + 1;
		Scum = cumsum(S,2);
		Ucum = cumsum(U,2);

		B = sum(U',1); #%3.2: B_j=delay of tran type j due to blocking

		B_total = dot(vec(f),vec(B)); #%3.1: B_total=mean blocking time over all trans classes

		R = S0 + sum(S'+U',1); #%3.1: R_j=latency of tran type J

		beta = Div(B,R); #%3.12: beta_j=fraction of time trans of type J are blocked=fraction of trans type J that are blocked

		L = zeros(convert(Int,J),convert(Int,regions)); #%3.6: L_ji= of locks held by trans of type j in DBR i
		Lp = zeros(convert(Int,J),convert(Int,regions)); #%3.6: Lp_ji=locks held by blocked trans of type j in DBR i
		buf = zeros(size(U))
		for i=1:regions
			gJK = cumsum(g[:,:,i],2);
			gJKshifted = gJK - g[:,:,i];
			buf = U+0;
			for j=1:length(gJKshifted)
				if gJKshifted[j] == 0
					buf[j] = 0
				end
			end
			tempDot = zeros(size(buf,1),1)
			for j= 1:size(buf,1)
				for k=1:size(buf,2)
					tempDot[j,1]=tempDot[j,1]+buf[j,k]*gJKshifted[j,k]
				end
				#tempDot[j,1] = dot(vec(buf[j,:]),vec(gJKshifted[j,:]))
			end
			Lp[:,i] = (1./ R') .* tempDot
			tempDot2 = zeros(size(S,1),1)
			for j= 1:size(S,1)
				for k=1:size(S,2)
					tempDot2[j,1]=tempDot2[j,1]+S[j,k]*gJK[j,k]
				end
				#tempDot2[j,1] = dot(vec(S[j,:]),vec(gJK[j,:]))
			end
			L[:,i] = (1./ R') .* tempDot2 + Lp[:,i];
		end

		R_total = dot(vec(f), vec(R)); #%3.2: R=mean delay over all trans classes
		M_total = T_total * R_total; #%3.3: M_total=total number of transactions in the system
		if 1==0
			if (M_total>ClientsNum)
				M_total = ClientsNum;
				newR_total = M_total / T_total;
				scaleR = newR_total / R_total;
				U = U * scaleR;
				S0 = S0 * scaleR;
				S = S * scaleR;
				R = S0 + sum(S',1) + sum(U',1);
				R_total = dot(vec(f), vec(R));
			end
		end

	    	if((typeof(T_total) <: Array) && size(f) == size(T_total))
			T = f * T_total;
		elseif(typeof(T_total) <: Array)
			T = f .* T_total[1];
		elseif(typeof(T_total) <: Number)
			T = f .*  T_total;
		end

		#%3.3: T_j=TPS of trans type j
		M = T .* R; #%3.3: M_j=# of trans type j in the system
		#Julia: change from / to ./
		beta_total = dot(vec(beta), vec(M ./ M_total)); #%3.11: beta_total=fraction of time that trans are blocked in the system, or fraction of blocked trans
		N = zeros(convert(Int,J),convert(Int,J),convert(Int,regions)); #%3.5: N_jli= of exclusive locks held by transations of type l in DBR i, (when considering a trans of type j)
		for j=1:J
			for i=1:regions
				N[:,j,i] = M[j] * L[j,i];
				if M[j]>=1
					N[j,j,i] = (M[j]-1) * L[j,i];
				else
					N[j,j,i] = 0;
				end
			end
		end

		P = zeros(convert(Int,J),convert(Int,maxK),convert(Int,J),convert(Int,regions)); #%3.4: P_jnli=prob of a conflict upon n'th lock request of tran type j with a tran of type l in DBR i
		for j=1:J
			for n=1:K[j]
				for jj=1:J
					for i=1:regions
					    P[j,n,jj,i] = Div(g[j,n,i] * N[j,jj,i], D[i]);
					end
				end
			end
		end

		Pcon = zeros(1, convert(Int,J)); #%3.12: Pcon_j=prob of lock conflict PER LOCK request for transactions of type j
		for j=1:J
			Pcon[j] = Div(sum(P[j,:,:,:]) , K[j]);
		end

		Pcon_total = dot(vec(Div(f .* K , K_total)), vec(Pcon)); #%3.13: Pcon_total=prob of lock conflict per lock request

		rhop = Div(Lp, L); #%3.9: rhop_ji=prob of a lock conflict with a blocked trans of type j which requested a lock in DBR i

		Q = zeros(convert(Int,J), convert(Int,regions)); #%3.11: Q_ji=prob of lock conflict by trans of type j requesting a lock in DBR i
		for j=1:J
			for i=1:regions
				Q[j,i] = (1/K[j]) * sum(P[j,:,:,i]);
			end
		end

		rho = Div(sum(rhop .* Q) , sum(Q)); #%3.10: rho=prob that a requested lock is held by a blocked trans when there is a lock conflict

		H = zeros(1, convert(Int,regions)); #% the normalization constant
		for i=1:regions
			H[i] = sum (f .* sum((S.*cumsum(g[:,:,i],2))',1));
		end

		Vp = zeros(1, convert(Int,regions)); #%3.6: Vp_i= mean waiting time w.r.t. the active transactions in DBR i
		for i=1:regions
			for j=1:J
				tempJ = 0;
				for n=1:K[j]
					tempK = 0;
					for m=1:n
						if Ucum[j,n]==Inf
							buf = Inf;
							if Ucum[j,K[j]] < Ucum[j,n]
								error("Something went wrong with a running sum!");
							end
						else
							buf = Ucum[j,K[j]]-Ucum[j,n];
						end
						tempK = tempK + bmult(g[j,m,i], (Scum[j,K[j]]-Scum[j,n]+S[j,n] + buf) );
					end
					tempJ = tempJ + S[j,n] * tempK;
				end
				Vp[i] = Vp[i] + bmult(f[j], tempJ);
			end
		end
		Vp = Div(Vp, H);

		Vp_total = 0; #%3.7: Vp_total=mean waiting time w.r.t. active transactions in all DBRs
		for j=1:J
			tempI = 0;
			for i=1:regions
				tempI = tempI + bmult(Vp[i], sum(P[j,:,:,i]));
			end
			Vp_total = Vp_total + Div(bmult(f[j], tempI), Pcon[j]);
		end
		Vp_total = Div(Vp_total , K_total);

		V = Vp + Div(Vp_total* rho*(1+rho) , (2*(1-rho)^2)); #%3.8: V_i=blocking time when transactions encounter a lock conflict in DBR i

		old_W = W;
		W = 0; #%3.14: W=mean waiting time per lock conflict
		for j=1:J
			tempI = 0;
			for i=1:regions
				buf = sum(P[j,:,:,i]);
				if buf!=0;
					tempI = tempI + V[i]*buf;
				end
			end
			W = W + Div(bmult(f[j], tempI) , Pcon[j]);
		end
		W = Div(W , K_total);

		deltaW = W - old_W;

		#%finally
		for j=1:J
			for n=1:K[j]
				tempI = 0;
				for i=1:regions
					tempI = tempI + bmult(V[i], sum(P[j,n,:,i]));
				end
				U[j,n] = tempI+0;  #%3.2: U_jn=mean delay incurred by trans type j when they encounter a lock conflict at step n
			end
		end

		if typeof(M_total) <: Array
			if M_total[1] * (1-beta_total) > cores
				S0 = givenS0 * M_total[1] * (1-beta_total) /cores;
				S = givenS .* M_total[1] * (1-beta_total) /cores;
			end
		else
			if M_total * (1-beta_total) > cores
				S0 = givenS0 * M_total * (1-beta_total) /cores;
				S = givenS .* M_total * (1-beta_total) /cores;
			end
		end
	end

	#%fprintf(1,'iter deltaW  R W T V=');
	#%[iter deltaW R W T V]
	#%U

	#%R
	#%T
	#%M_total

	return R,T_total,M_total,Vp,V,W,U,Pcon,L

end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
using Distributions: pdf, Binomial, rand
import Plots
using Plots: plot, plot!, scatter!, histogram
using StatisticalRethinking: hpdi, PI
using Gen: @gen, binom, (~)

"""
Suppose the globe tossing data (Chapter 2) had turned out to be 4 water
and 11 land. Construct the posterior distribution, using grid approximation.
Use the same flat prior as in the book.
"""
function problem1(grid_length=20)

    # define grid
    p_grid = [n for n in range(0,1,grid_length)]
    @debug "p_grid" p_grid

    # define prior
    prior = fill(1.0, grid_length)
    @debug "prior" prior

    # Distribution
    binomials =  Binomial.(15, p_grid)
    @debug "binomials" binomials

    # compute likelihood at each value in grid
    likelihood = pdf.(binomials, 4)
    @debug "likelihood" likelihood

    # compute product of likelihood and prior
    unstd_posterior = likelihood .* prior
    @debug "unstd_posterior" unstd_posterior

    # standardize the posterior, so it sums to 1
    posterior = unstd_posterior ./ sum(unstd_posterior)
    @debug "posterior" posterior

    return posterior
end

"""
Now suppose the data are 4 water and 2 land. Compute the posterior
again, but this time use a prior that is zero below p = 0.5 and a constant
above p = 0.5. This corresponds to prior information that a majority of the
Earth’s surface is water.
"""
function problem1_plot(grid_length=20)
    xs = range(0,1,length=grid_length)
    data = problem1(grid_length)
    title = "HW1 Problem 1 (Gridsize=" * string(grid_length) * ")"
    xlabel="probability of water"
    ylabel="posterior probability"
    labels="y"
    plot(xs, data, title=title, labels=labels, xlabel=xlabel, ylabel=ylabel)
end

function problem2(grid_length=20)
     # define grid
    p_grid = [n for n in range(0,1,grid_length)]
    @debug "p_grid" p_grid

    # define prior
    prior = [g < 0.5 ? 0. : 1. for g in p_grid]
    @debug "prior" prior

    # Distribution
    binomials =  Binomial.(6, p_grid)
    @debug "binomials" binomials

    # compute likelihood at each value in grid
    likelihood = pdf.(binomials, 4)
    @debug "likelihood" likelihood

    # compute product of likelihood and prior
    unstd_posterior = likelihood .* prior
    @debug "unstd_posterior" unstd_posterior

    # standardize the posterior, so it sums to 1
    posterior = unstd_posterior ./ sum(unstd_posterior)
    @debug "posterior" posterior

    return posterior
end

function problem2_plot(grid_length=20)
    xs = range(0,1,length=grid_length)
    data = problem2(grid_length)
    title = "HW1 Problem 2 (Gridsize=" * string(grid_length) * ")"
    xlabel="probability of water"
    ylabel="posterior probability"
    labels="y"
    plot(xs, data, title=title, labels=labels, xlabel=xlabel, ylabel=ylabel)
end

"""
For the posterior distribution from 2, compute 89% percentile and HPDI
intervals. Compare the widths of these intervals. Which is wider? Why? If
you had only the information in the interval, what might you misunderstand
about the shape of the posterior distribution
"""
function problem3(grid_length=20)
    posterior = problem2(grid_length)
    return (posterior, PI(posterior;perc_prob=0.89), hpdi(posterior;alpha=0.11))
end

function problem3_plot(grid_length=20)
    posterior, percentile, highprobabilitydensityinterval = problem3(grid_length)
    xs = Vector(range(0,1,grid_length))
    plot(xs, posterior, layout=2)
    scatter!(percentile, posterior, fillrange=0., fillalpha=0.3, fillcolor=:blue)
end

"""
Suppose there is bias in sampling so that Land
is more likely than Water to be recorded. Specifically, assume that 1-in-5
(20%) of Water samples are accidentally recorded instead as ”Land”. First,
write a generative simulation of this sampling process. Assuming the true
proportion of Water is 0.70, what proportion does your simulation tend to
produce instead? Second, using a simulated sample of 20 tosses, compute
the unbiased posterior distribution of the true proportion of water.
"""
function problem4(samples=100, tosses=20, proportion_water=0.7, bias=0.8)
    biased_sample = zeros(tosses+1)
    biased_true = pdf.(Binomial(tosses, proportion_water * bias), 0:tosses)
    unbiased = pdf.(Binomial(tosses, proportion_water), 0:tosses)
    for _ in 1:samples
        sample = biased_water(tosses, proportion_water, bias)
        biased_sample[sample] += 1
    end

    biased_sample = biased_sample ./ samples

    return (biased_sample, biased_true, unbiased)
end

@gen function biased_water(tosses=20, proportion_water=0.7, bias=0.8)

    true_water ~ binom(tosses, proportion_water)
    recorded_water ~ binom(true_water, bias)

    return recorded_water
end

function problem4_plot(samples=100, tosses=20, proportion_water=0.7, bias=0.8)
    xs = range(0, tosses, step=1)
    @debug xs

    biased_sample, biased_true, unbiased = problem4(samples, tosses, proportion_water, bias)

    labels = "biased sample"
    xlabel="probability of water"
    ylabel="posterior probability"
    seriestype=:bar

    biased_sample_plot = plot(xs, biased_sample, seriestype=seriestype, color=:blue, title="Biased (" * string(samples) * " samples)", labels=labels, xlabel=xlabel, ylabel=ylabel)
    @debug biased_sample_plot

    biased_true_plot = plot(xs, biased_true, seriestype=seriestype, color=:green, title="Biased", labels=labels, xlabel=xlabel, ylabel=ylabel)
    @debug biased_true_plot

    unbiased_plot = plot(xs, unbiased, seriestype=seriestype, color=:red, title="Unbiased", labels=labels, xlabel=xlabel, ylabel=ylabel)
    @debug unbiased_plot

    overlay_plot = plot(xs, [biased_sample biased_true unbiased], seriestype=seriestype, seriesalpha=[0.33 0.33 0.33], color=[:blue :green :red], labels=["biased (sample)" "biased" "unbiased"], xlabel=xlabel, ylabel=ylabel)
    @debug overlay_plot

    plot(biased_sample_plot, biased_true_plot, unbiased_plot, overlay_plot, layout=4, legend=false)
end

ext: jl


text:
module Underscores

@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    read(path, String)
end Underscores

export @_

isquoted(ex) = ex isa Expr && ex.head in (:quote, :inert, :meta)

function _replacesyms(sym_map, ex)
    if ex isa Symbol
        return sym_map(ex)
    elseif ex isa Expr
        if isquoted(ex)
            return ex
        end
        args = map(e->_replacesyms(sym_map, e), ex.args)
        return Expr(ex.head, args...)
    else
        return ex
    end
end

function add_closures(ex, prefix, pattern)
    if ex isa Expr && (ex.head == :kw || ex.head == :parameters)
        return Expr(ex.head, map(e->add_closures(e,prefix,pattern), ex.args)...)
    end
    plain_nargs = false
    numbered_nargs = 0
    body = _replacesyms(ex) do sym
        m = match(pattern, string(sym))
        if m === nothing
            sym
        else
            argnum_str = m[1]
            if isempty(argnum_str)
                plain_nargs = true
                argnum = 1
            else
                if !isdigit(argnum_str[1])
                    argnum_str = map(c->c-'₀'+'0', argnum_str)
                end
                argnum = parse(Int, argnum_str)
                numbered_nargs = max(numbered_nargs, argnum)
            end
            Symbol(prefix, argnum)
        end
    end
    if plain_nargs && numbered_nargs > 0
        throw(ArgumentError("Cannot mix plain and numbered `$prefix` placeholders in `$ex`"))
    end
    nargs = max(plain_nargs, numbered_nargs)
    if nargs == 0
        return ex
    end
    argnames = map(i->Symbol(prefix,i), 1:nargs)
    return :(($(argnames...),) -> $body)
end

replace_(ex)  = add_closures(ex, "_", r"^_([0-9]*|[₀-₉]*)$")
replace__(ex) = add_closures(ex, "__", r"^__([0-9]*|[₀-₉]*)$")

const _square_bracket_ops = [:comprehension, :typed_comprehension, :generator,
                             :ref, :vcat, :typed_vcat, :hcat, :typed_hcat, :row]

_isoperator(x) = x isa Symbol && Base.isoperator(x)

function lower_inner(ex)
    if ex isa Expr
        if ex.head == :(=) ||
            (ex.head == :call && length(ex.args) > 1 && _isoperator(ex.args[1]))
            # Infix operators do not count as outermost function call
            return Expr(ex.head, ex.args[1],
                map(lower_inner, ex.args[2:end])...)
        elseif ex.head in _square_bracket_ops || ex.head == :if
            # Indexing & other square brackets not counted as outermost function
            return Expr(ex.head, map(lower_inner, ex.args)...)
        elseif ex.head == :. && length(ex.args) == 2 && ex.args[2] isa QuoteNode
            # Getproperty also doesn't count
            return Expr(ex.head, map(lower_inner, ex.args)...)
        elseif ex.head == :. && length(ex.args) == 2 &&
            ex.args[2] isa Expr && ex.args[2].head == :tuple
            # Broadcast calls treated as normal calls for underscore lowering
            return Expr(ex.head, replace_(ex.args[1]),
                                  Expr(:tuple, map(replace_, ex.args[2].args)...))
        else
            # For other syntax, replace _ in args individually
            return Expr(ex.head, map(replace_, ex.args)...)
        end
    else
        return ex
    end
end

# In principle this can be extended locally by a package for use within the
# package and for prototyping purposes. However note that this will interact
# badly with precompilation. (If it makes sense we could fix this per-package
# by storing a per-module _pipeline_ops in the module using @_.)
const _pipeline_ops = [:|>, :<|, :∘, :.|>, :.<|]

function lower_underscores(ex)
    if ex isa Expr
        if isquoted(ex)
            return ex
        elseif ex.head == :call && length(ex.args) > 1 &&
               ex.args[1] in _pipeline_ops
            # Special case for pipelining and composition operators
            return Expr(ex.head, ex.args[1],
                        map(lower_underscores, ex.args[2:end])...)
        elseif ex.head == :do
            error("@_ expansion for `do` syntax is reserved")
        else
            # For other syntax, replace __ over the entire expression
            return replace__(lower_inner(ex))
        end
    else
        return ex
    end
end

"""
    @_ func(ex1, [ex2 ...])

Convert `ex1,ex2,...` into anonymous functions when they have `_` placeholders,
and *pass them along* to `func`.

The detailed rules are:
1. Uses of the placeholder `_` expand to the single argument of an anonymous
   function which is passed to the outermost ordinary function call.
2. Numbered placeholders `_1,_2,...` (or `_₁,_₂,...`) may be used if you need
   more than one argument. Numbers indicate position in the argument list.
3. The double underscore placeholder `__` (and numbered versions `__1,__2,...`)
   expands the closure scope to the whole expression.
4. Piping and composition chains with `|>,<|,∘` are treated as a special case
   where the replacement recurses into sub-expressions.

These rules imply the following equivalences

| Expression                 |  Rules  | Meaning                        |
|:-------------------------- |:------- |:------------------------------ |
| `@_ map(_+1, a)`           | (1)     | `map(x->x+1, a)`               |
| `@_ map(_^_, a)`           | (1)     | `map(x->x^x, a)`               |
| `@_ map(_2/_1, a, b)`      | (1,2)   | `map((x,y)->y/x, a, b)`        |
| `@_ func(a,__,b)`          | (3)     | `x->func(a,x,b)`               |
| `@_ func(a,__2,b)`         | (3)     | `(x,y)->func(a,y,b)`           |
| `@_ data \\|> map(_.f,__)` | (1,3,4) | `data \\|> (d->map(x->x.f,d))` |

# Extended help

## Examples

`@_` can be used for simple mapping operations in cases where broadcasting
syntax is awkward. For example, to get the second last element of each array in
a collection:

```jldoctest
julia> @_ map(_[end-1],  [[1,2,3], [4,5]])
2-element Array{Int64,1}:
 2
 4
```

If you need to repeat an argument more than once, just use `_` multiple times:

```jldoctest
julia> @_ map(_^_,  [1,2,3])
3-element Array{Int64,1}:
  1
  4
 27
```

For manipulating tabular data `@_` provides convenient syntax which is
especially useful when combined with double underscore placeholders `__` and
piping syntax. Think of `__` as the table, and `_` as an individual row:

```jldoctest
julia> table = [(x="a", y=1),
                (x="b", y=2),
                (x="c", y=3)];

julia> @_ table |>
          filter(!startswith(_.x, "a"), __) |>
          map(_.y, __)
2-element Array{Int64,1}:
 2
 3
```

## Extraordinary functions

The scope of `_` as described in rule 1 depends on "ordinary" function call.
This excludes the following operations:

* Square brackets: In `map(_^2,__)[3]`, it is `map` which receives an anonymous
  function, as this happens before the indexing is lowered to `getindex(...,3)`.
  Comprehensions (`collect`) and explicit matrix constructions (`hvcat`) are
  treated similarly.

* Broadcasting, and field access: In `f.(_,xs)` and `f(_,x).y` the function `f`
  is the ordinary call, not the internal `broadcast` or `getproperty`.

* Infix operators: While `sum(_^2,x) / length(x)` can be written in prefix form
  `/(...,...)`, the convention of `@_` is not to view this as an ordinary call,
  and hence to pass the anonymous function to `sum` instead. This also applies to
  broadcasted operators, such as `map(_^2,x) ./ length(x)`.

* If statements, including the ternary operator. Note that this has higher
  precedence than pipes: `data |> (any(_.x<0, __) ? abs.(__) : __) |> step`
  needs these brackets.

The scope of `__` is unaffected by these concerns.

| Expression                       | Meaning                            |
|:-------------------------------- |:---------------------------------- |
| `@_ data \\|> map(_[2],__)[3]`   | `data \\|> (d->map(x->x[2],d)[3])` |
| `@_ [sum(_*_, z) for z in a]`    | `[sum(x->x*x, z) for z in a]`      |
| `@_ sum(1+_^2, data).re`         | `sum(x->1+x^2, data).re`           |
| `@_ sum(_^2,a) / length(a)`      | `sum(x->x^2,a) / length(a)`        |
| `@_ /(sum(_^2,a), length(a))`    | The same, infix form is canonical. |
| `@_ data \\|> filter(_>3,__).^2` | `data \\|> d->(filter(>(3),d).^2)` |
| `@_ any(_>3,xs) ? 0 : map(_,ys)` | `any(x->x>3,xs) ? 0 : map(y->y,ys)`|

"""
macro _(ex)
    esc(lower_underscores(ex))
end

end


text:
function estimatelearningcore(y,x,Choice,bstart,N,T,S,J,BigN,Delta,sig,Resid,resid,tresid,Psi1,abilsub,Csum)
    abil=zeros(N*S,J)

    idelta   = inv(Delta)
    vtemp2   = zeros(J,J)
    vabil    = zeros(N*S,J)
    vabilw   = zeros(N*S,J)
    for j=1:J
        tresid[:,j] = (sum(reshape(Resid[:,j],(T,N*S));dims=1))'
    end

    for i=1:S*N
        psit=Psi1[i,:]
        Psi = zeros(J,J)
        for j=1:J
            Psi[j,j] = psit[j]./sig[j]
        end

        vtemp=inv(idelta.+Psi)

        vectres = zeros(J)
        for j=1:J
            vectres[j] = tresid[i,j]./sig[j]
        end
        temp=(vtemp*vectres)'
        abil[i,:]=temp

        for j=1:J
            vabil[i,j] = vtemp[j,j] # later need PTypesl[i]*vtemp[j,j]
        end

        vtemp2.+=(vtemp.+temp'*temp) # later need PTypesl[i]*(vtemp+temp'*temp)
    end

    Delta=vtemp2./N

    vabilw    = deepcopy(vabil)
    Abil      = kron(abil,ones(T,1))
    Vabil     = kron(vabil,ones(T,1))

    sigdem = zeros(J)
    for j=1:J
        abilsub[j] = Abil[vec(Choice.==j),j]
        sigdem[j]  = sum((resid[j].-abilsub[j]).^2) # later need sum(PTypesub[j].*(resid[j].-abilsub[j]).^2)
    end

    sig=((sum(Csum.*vabilw;dims=1))'.+sigdem)./BigN

    for j=1:J
        flag                     = (y[:,j].!=999)
        bstart[:,j]              = x[flag,:,j]\(y[flag,j].-abilsub[j]) # later weight by PTypesub[j]
        resid[j]                 = y[flag,j].-x[flag,:,j]*bstart[:,j]
        Resid[vec(Choice.==j),j] = resid[j]
    end

    return bstart,Delta,sig,Resid,resid,tresid,abilsub
end


text:
import Base: rand,
             ∈,
             split

export Zonotope,
       order,
       linear_map,
       scale,
       ngens,
       reduce_order,
       constraints_list,
       generators

"""
    Zonotope{N<:Real} <: AbstractZonotope{N}

Type that represents a zonotope.

### Fields

- `center`     -- center of the zonotope
- `generators` -- matrix; each column is a generator of the zonotope

### Notes

Mathematically, a zonotope is defined as the set

```math
Z = \\left\\{ c + ∑_{i=1}^p ξ_i g_i,~~ ξ_i \\in [-1, 1]~~ ∀ i = 1,…, p \\right\\},
```
where ``c \\in \\mathbb{R}^n`` is its *center* and ``\\{g_i\\}_{i=1}^p``,
``g_i \\in \\mathbb{R}^n``, is the set of *generators*.
This characterization defines a zonotope as the finite Minkowski sum of line
segments.
Zonotopes can be equivalently described as the image of a unit infinity-norm
ball in ``\\mathbb{R}^n`` by an affine transformation.

- `Zonotope(center::AbstractVector{N},
            generators::AbstractMatrix{N}) where {N<:Real}`

- `Zonotope(center::AbstractVector{N},
            generators_list::AbstractVector{VN}
           ) where {N<:Real, VN<:AbstractVector{N}}`

The optional argument `remove_zero_generators` controls whether we remove zero
columns from the `generators` matrix.
This option is active by default.

### Examples

A two-dimensional zonotope with given center and set of generators:

```jldoctest zonotope_label
julia> Z = Zonotope([1.0, 0.0], [0.1 0.0; 0.0 0.1])
Zonotope{Float64}([1.0, 0.0], [0.1 0.0; 0.0 0.1])
julia> dim(Z)
2
```

Compute its vertices:

```jldoctest zonotope_label
julia> vertices_list(Z)
4-element Array{Array{Float64,1},1}:
 [1.1, 0.1]
 [0.9, 0.1]
 [1.1, -0.1]
 [0.9, -0.1]
```

Evaluate the support vector in a given direction:

```jldoctest zonotope_label
julia> σ([1., 1.], Z)
2-element Array{Float64,1}:
 1.1
 0.1
```

Alternative constructor: A zonotope in two dimensions with three generators:

```jldoctest
julia> Z = Zonotope(ones(2), [[1., 0.], [0., 1.], [1., 1.]])
Zonotope{Float64}([1.0, 1.0], [1.0 0.0 1.0; 0.0 1.0 1.0])
julia> Z.generators
2×3 Array{Float64,2}:
 1.0  0.0  1.0
 0.0  1.0  1.0
```
"""
struct Zonotope{N<:Real} <: AbstractZonotope{N}
    center::AbstractVector{N}
    generators::AbstractMatrix{N}

    function Zonotope(center::AbstractVector{N}, generators::AbstractMatrix{N};
                      remove_zero_generators::Bool=true) where {N<:Real}
        if remove_zero_generators
            generators = delete_zero_columns!(generators)
        end
        new{N}(center, generators)
    end
end

# constructor from center and list of generators
Zonotope(center::AbstractVector{N}, generators_list::AbstractVector{VN};
         remove_zero_generators::Bool=true
        ) where {N<:Real, VN<:AbstractVector{N}} =
    Zonotope(center, hcat(generators_list...);
             remove_zero_generators=remove_zero_generators)

# --- AbstractCentrallySymmetric interface functions ---

"""
    center(Z::Zonotope{N})::Vector{N} where {N<:Real}

Return the center of a zonotope.

### Input

- `Z` -- zonotope

### Output

The center of the zonotope.
"""
function center(Z::Zonotope{N})::Vector{N} where {N<:Real}
    return Z.center
end

# --- AbstractZonotope interface functions ---

"""
   genmat(Z::Zonotope)

Return the generator matrix of a zonotope.

### Input

- `Z` -- zonotope

### Output

A matrix where each column represents one generator of the zonotope `Z`.
"""
function genmat(Z::Zonotope)
    return Z.generators
end

"""
    generators(Z::Zonotope)

Return an iterator over the generators of a zonotope.

### Input

- `Z` -- zonotope

### Output

An iterator over the generators of `Z`.
"""
function generators(Z::Zonotope)
    return generators_fallback(Z)
end

# --- AbstractPolytope interface functions ---

"""
    vertices_list(Z::Zonotope{N})::Vector{Vector{N}} where {N<:Real}

Return the vertices of a zonotope.

### Input

- `Z` -- zonotope

### Output

List of vertices as a vector of vectors.

### Algorithm

If the zonotope has ``p`` generators, each of the ``2^p`` vertices is computed
by taking the sum of the center and a linear combination of generators, where
the combination factors are ``ξ_i ∈ \\{-1, 1\\}``.

### Notes

For high dimensions, it would be preferable to develop a `vertex_iterator`
approach.
"""
function vertices_list(Z::Zonotope{N})::Vector{Vector{N}} where {N<:Real}
    p = ngens(Z)
    vlist = Vector{Vector{N}}()
    sizehint!(vlist, 2^p)

    for ξi in Iterators.product([[1, -1] for i = 1:p]...)
        push!(vlist, Z.center .+ Z.generators * collect(ξi))
    end

    return vlist
end

# --- LazySet interface functions ---

"""
    ρ(d::AbstractVector{N}, Z::Zonotope{N}) where {N<:Real}

Return the support function of a zonotope in a given direction.

### Input

- `d` -- direction
- `Z` -- zonotope

### Output

The support function of the zonotope in the given direction.

### Algorithm

The support value is ``cᵀ d + ‖Gᵀ d‖₁`` where ``c`` is the center and ``G`` is
the generator matrix of `Z`.

"""
function ρ(d::AbstractVector{N}, Z::Zonotope{N}) where {N<:Real}
    return dot(center(Z), d) + sum(abs.(transpose(Z.generators) * d))
end

"""
    σ(d::AbstractVector{N}, Z::Zonotope{N}) where {N<:Real}

Return the support vector of a zonotope in a given direction.

### Input

- `d` -- direction
- `Z` -- zonotope

### Output

Support vector in the given direction.
If the direction has norm zero, the vertex with ``ξ_i = 1 \\ \\ ∀ i = 1,…, p``
is returned.
"""
function σ(d::AbstractVector{N}, Z::Zonotope{N}) where {N<:Real}
    return Z.center .+ Z.generators * sign_cadlag.(_At_mul_B(Z.generators, d))
end

"""
    ∈(x::AbstractVector{N}, Z::Zonotope{N};
      solver=GLPKSolverLP(method=:Simplex))::Bool where {N<:Real}

Check whether a given point is contained in a zonotope.

### Input

- `x`      -- point/vector
- `Z`      -- zonotope
- `solver` -- (optional, default: `GLPKSolverLP(method=:Simplex)`) the backend
              used to solve the linear program

### Output

`true` iff ``x ∈ Z``.

### Examples

```jldoctest
julia> Z = Zonotope([1.0, 0.0], [0.1 0.0; 0.0 0.1]);

julia> ∈([1.0, 0.2], Z)
false
julia> ∈([1.0, 0.1], Z)
true
```

### Algorithm

The membership problem is computed by stating and solving the following linear
program with the simplex method.
Let ``p`` and ``n`` be the number of generators and ambient dimension,
respectively.
We consider the minimization of ``x_0`` in the ``p+1``-dimensional space of
elements ``(x_0, ξ_1, …, ξ_p)`` constrained to ``0 ≤ x_0 ≤ ∞``,
``ξ_i ∈ [-1, 1]`` for all ``i = 1, …, p``, and such that ``x-c = Gξ`` holds.
If a feasible solution exists, the optimal value ``x_0 = 0`` is achieved.

### Notes

This function is parametric in the number type `N`. For exact arithmetic use
an appropriate backend, e.g. `solver=GLPKSolverLP(method=:Exact)`.
"""
function ∈(x::AbstractVector{N}, Z::Zonotope{N};
           solver=GLPKSolverLP(method=:Simplex))::Bool where {N<:Real}
    @assert length(x) == dim(Z)

    p, n = ngens(Z), dim(Z)
    # (n+1) x (p+1) matrix with block-diagonal blocks 1 and Z.generators
    A = [[one(N); zeros(N, p)]'; [zeros(N, n) Z.generators]]
    b = [zero(N); (x - Z.center)]
    lbounds = [zero(N); fill(-one(N), p)]
    ubounds = [N(Inf); ones(N, p)]
    sense = ['>'; fill('=', n)]
    obj = [one(N); zeros(N, p)]

    lp = linprog(obj, A, sense, b, lbounds, ubounds, solver)
    return (lp.status == :Optimal) # Infeasible or Unbounded => false
end

"""
    rand(::Type{Zonotope}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
        )::Zonotope{N}

Create a random zonotope.

### Input

- `Zonotope`       -- type for dispatch
- `N`              -- (optional, default: `Float64`) numeric type
- `dim`            -- (optional, default: 2) dimension
- `rng`            -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`           -- (optional, default: `nothing`) seed for reseeding
- `num_generators` -- (optional, default: `-1`) number of generators of the
                      zonotope (see comment below)

### Output

A random zonotope.

### Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1.

The number of generators can be controlled with the argument `num_generators`.
For a negative value we choose a random number in the range `dim:2*dim` (except
if `dim == 1`, in which case we only create a single generator).
"""
function rand(::Type{Zonotope};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing,
              num_generators::Int=-1
             )::Zonotope{N}
    rng = reseed(rng, seed)
    center = randn(rng, N, dim)
    if num_generators < 0
        num_generators = (dim == 1) ? 1 : rand(dim:2*dim)
    end
    generators = randn(rng, N, dim, num_generators)
    return Zonotope(center, generators)
end

# --- Zonotope functions ---

"""
    order(Z::Zonotope)::Rational

Return the order of a zonotope.

### Input

- `Z` -- zonotope

### Output

A rational number representing the order of the zonotope.

### Notes

The order of a zonotope is defined as the quotient of its number of generators
and its dimension.
"""
function order(Z::Zonotope)::Rational
    return ngens(Z) // dim(Z)
end

"""
    linear_map(M::AbstractMatrix{N}, Z::Zonotope{N}) where {N<:Real}

Concrete linear map of a zonotope.

### Input

- `M` -- matrix
- `Z` -- zonotope

### Output

The zonotope obtained by applying the linear map to the center and generators
of ``Z``.
"""
function linear_map(M::AbstractMatrix{N}, Z::Zonotope{N}) where {N<:Real}
    @assert dim(Z) == size(M, 2) "a linear map of size $(size(M)) cannot be " *
                                 "applied to a set of dimension $(dim(Z))"

    c = M * Z.center
    gi = M * Z.generators
    return Zonotope(c, gi)
end

"""
    scale(α::Real, Z::Zonotope)

Concrete scaling of a zonotope.

### Input

- `α` -- scalar
- `Z` -- zonotope

### Output

The zonotope obtained by applying the numerical scale to the center and
generators of ``Z``.
"""
function scale(α::Real, Z::Zonotope)
    c = α .* Z.center
    gi = α .* Z.generators
    return Zonotope(c, gi)
end

"""
    ngens(Z::Zonotope)::Int

Return the number of generators of a zonotope.

### Input

- `Z` -- zonotope

### Output

Integer representing the number of generators.
"""
ngens(Z::Zonotope)::Int = size(Z.generators, 2)

"""
    reduce_order(Z::Zonotope, r)::Zonotope

Reduce the order of a zonotope by overapproximating with a zonotope with less
generators.

### Input

- `Z` -- zonotope
- `r` -- desired order

### Output

A new zonotope with less generators, if possible.

### Algorithm

This function implements the algorithm described in A. Girard's
*Reachability of Uncertain Linear Systems Using Zonotopes*, HSCC. Vol. 5. 2005.

If the desired order is smaller than one, the zonotope is *not* reduced.
"""
function reduce_order(Z::Zonotope{N}, r)::Zonotope{N} where {N<:Real}
    c, G = Z.center, Z.generators
    d, p = dim(Z), ngens(Z)

    if r * d >= p || r < 1
        # do not reduce
        return Z
    end

    h = zeros(N, p)
    for i in 1:p
        h[i] = norm(G[:, i], 1) - norm(G[:, i], Inf)
    end
    ind = sortperm(h)

    m = p - floor(Int, d * (r - 1)) # subset of ngens that are reduced
    rg = G[:, ind[1:m]] # reduced generators

    # interval hull computation of reduced generators
    Gbox = Diagonal(sum(abs.(rg), dims=2)[:])
    if m < p
        Gnotred = G[:, ind[m+1:end]]
        Gred = [Gnotred Gbox]
    else
        Gred = Gbox
    end
    return Zonotope(c, Gred)
end

"""
    split(Z::Zonotope, j::Int)

Return two zonotopes obtained by splitting the given zonotope.

### Input

- `Z` -- zonotope
- `j` -- index of the generator to be split

### Output

The zonotope obtained by splitting `Z` into two zonotopes such that
their union is `Z` and their intersection is possibly non-empty.

### Algorithm

This function implements [Prop. 3, 1], that we state next. The zonotope
``Z = ⟨c, g^{(1, …, p)}⟩`` is split into:

```math
Z₁ = ⟨c - \\frac{1}{2}g^{(j)}, (g^{(1, …,j-1)}, \\frac{1}{2}g^{(j)}, g^{(j+1, …, p)})⟩ \\\\
Z₂ = ⟨c + \\frac{1}{2}g^{(j)}, (g^{(1, …,j-1)}, \\frac{1}{2}g^{(j)}, g^{(j+1, …, p)})⟩,
```
such that ``Z₁ ∪ Z₂ = Z`` and ``Z₁ ∩ Z₂ = Z^*``, where

```math
Z^* = ⟨c, (g^{(1,…,j-1)}, g^{(j+1,…, p)})⟩.
```

[1] *Althoff, M., Stursberg, O., & Buss, M. (2008). Reachability analysis of
nonlinear systems with uncertain parameters using conservative linearization.
In Proc. of the 47th IEEE Conference on Decision and Control.*
"""
function split(Z::Zonotope, j::Int)
    @assert 1 <= j <= ngens(Z) "cannot split a zonotope with $(ngens(Z)) generators along index $j"
    c, G = Z.center, Z.generators
    Gj = G[:, j]
    Gj_half = Gj / 2

    c₁ = c - Gj_half
    c₂ = c + Gj_half

    G₁ = copy(G)
    G₁[:, j] = Gj_half
    G₂ = copy(G₁)

    Z₁ = Zonotope(c₁, G₁)
    Z₂ = Zonotope(c₂, G₂)
    return Z₁, Z₂
end

"""
    constraints_list(P::Zonotope{N}) where {N<:Real}

Return the list of constraints defining a zonotope.

### Input

- `Z` -- zonotope

### Output

The list of constraints of the zonotope.

### Algorithm

This is the (inefficient) fallback implementation for rational numbers.
It first computes the vertices and then converts the corresponding polytope
to constraint representation.
"""
function constraints_list(Z::Zonotope{N}) where {N<:Real}
    return constraints_list(VPolytope(vertices_list(Z)))
end

"""
    constraints_list(Z::Zonotope{N}) where {N<:AbstractFloat}

Return the list of constraints defining a zonotope.

### Input

- `Z`               -- zonotope
- `check_full_rank` -- (optional; default: `true`) flag for checking whether the
                       generator matrix has full rank

### Output

The list of constraints of the zonotope.

### Notes

The algorithm assumes that no generator is redundant.
The result has ``2 \\binom{p}{n-1}`` (with ``p`` being the number of generators
and ``n`` being the ambient dimension) constraints, which is optimal under this
assumption.

If ``p < n`` or the generator matrix is not full rank, we fall back to the
(slower) computation based on the vertex representation.

### Algorithm

We follow the algorithm presented in *Althoff, Stursberg, Buss: Computing
Reachable Sets of Hybrid Systems Using a Combination of Zonotopes and Polytopes.
2009.*

The one-dimensional case is not covered by that algorithm; we manually handle
this case, assuming that there is only one generator.
"""
function constraints_list(Z::Zonotope{N}; check_full_rank::Bool=true
                         ) where {N<:AbstractFloat}
    G = genmat(Z)
    p = ngens(Z)
    n = dim(Z)

    # use fallback implementation if order < 1 or matrix is not full rank
    if p < n || (check_full_rank && rank(G) < n)
        return invoke(constraints_list, Tuple{Zonotope{<:Real}}, Z)
    end

    # special handling of 1D case
    if n == 1
        if p > 1
            error("1D-zonotope constraints currently only support a single " *
                  "generator")
        end

        c = Z.center[1]
        g = G[:, 1][1]
        constraints = [LinearConstraint([N(1)], c + g),
                       LinearConstraint([N(-1)], g - c)]
        return constraints
    end

    i = 0
    c = Z.center
    m = binomial(p, n - 1)
    constraints = Vector{LinearConstraint{N, Vector{N}}}(undef, 2 * m)
    for columns in StrictlyIncreasingIndices(p, n-1)
        i += 1
        c⁺ = cross_product(view(G, :, columns))
        normalize!(c⁺, 2)

        Δd = sum(abs.(transpose(G) * c⁺))

        d⁺ = dot(c⁺, c) + Δd
        c⁻ = -c⁺
        d⁻ = -d⁺ + 2 * Δd  # identical to dot(c⁻, c) + Δd

        constraints[i] = LinearConstraint(c⁺, d⁺)
        constraints[i + m] = LinearConstraint(c⁻, d⁻)
    end
    @assert i == m "expected 2*$m constraints, but only created 2*$i"
    return constraints
end

"""
    translate(Z::Zonotope{N}, v::AbstractVector{N}; share::Bool=false
             ) where {N<:Real}

Translate (i.e., shift) a zonotope by a given vector.

### Input

- `Z`     -- zonotope
- `v`     -- translation vector
- `share` -- (optional, default: `false`) flag for sharing unmodified parts of
             the original set representation

### Output

A translated zonotope.

### Notes

The generator matrix is shared with the original zonotope if `share == true`.

### Algorithm

We add the vector to the center of the zonotope.
"""
function translate(Z::Zonotope{N}, v::AbstractVector{N}; share::Bool=false
                  ) where {N<:Real}
    @assert length(v) == dim(Z) "cannot translate a $(dim(Z))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    c = center(Z) + v
    generators = share ? Z.generators : copy(Z.generators)
    return Zonotope(c, generators)
end


text:
@testset "Dicts can be diffed" begin
    d1 = Dict(
        :foo => "foo",
        :bar => "bar",
        :baz => Dict(
            :fizz => "fizz",
            :buzz => "buzz"
        )
    )

    @testset "One Changed" begin
        d = deepdiff(d1, Dict(
            :foo => "foo",
            :bar => "biz",
            :baz => Dict(
                :fizz => "fizz",
                :buzz => "buzz"
            )
        ))
        @test added(d) == Set()
        @test removed(d) == Set()
        @test changed(d) == Dict{Symbol, DeepDiffs.DeepDiff}(:bar => deepdiff("bar", "biz"))
    end

    @testset "One Removed" begin
        d = deepdiff(d1, Dict(
            :foo => "foo",
            :baz => Dict(
                :fizz => "fizz",
                :buzz => "buzz"
            )
        ))
        @test added(d) == Set()
        @test removed(d) == Set([:bar])
        @test changed(d) == Dict()
    end

    @testset "One Added" begin
        d = deepdiff(d1, Dict(
            :foo => "foo",
            :bar => "bar",
            :biz => "biz",
            :baz => Dict(
                :fizz => "fizz",
                :buzz => "buzz"
            )
        ))
        @test added(d) == Set([:biz])
        @test removed(d) == Set()
        @test changed(d) == Dict()
    end

    @testset "Inner Dict Modified" begin
        d = deepdiff(d1, Dict(
            :foo => "foo",
            :bar => "bar",
            :baz => Dict(
                :fizz => "fizz",
                :buzz => "bizzle"
            )
        ))
        @test added(d) == Set()
        @test removed(d) == Set()
        @test changed(d) == Dict{Symbol, DeepDiffs.DeepDiff}(:baz => deepdiff(
            Dict(
                :fizz => "fizz",
                :buzz => "buzz"
            ),
            Dict(
                :fizz => "fizz",
                :buzz => "bizzle"
            )
        ))
    end

    @testset "Totally Removed" begin
        d = deepdiff(d1, Dict())
        @test added(d) == Set()
        @test removed(d) == Set([:foo, :bar, :baz])
        @test changed(d) == Dict()
    end

    @testset "Totally added" begin
        d = deepdiff(Dict(), d1)
        @test added(d) == Set([:foo, :bar, :baz])
        @test removed(d) == Set()
        @test changed(d) == Dict()
    end

    @testset "missing/NaN" begin
        dnan = Dict(:d => NaN)
        d = deepdiff(dnan, dnan)
        @test added(d) == removed(d) == Set()
        @test changed(d) == Dict()

        dmis = Dict(:d => missing)
        d = deepdiff(dmis, dmis)
        @test added(d) == removed(d) == Set()
        @test changed(d) == Dict()

        dnank = Dict(NaN => true)
        d = deepdiff(dnan, dnan)
        @test added(d) == removed(d) == Set()
        @test changed(d) == Dict()

        dmisk = Dict(missing => true)
        d = deepdiff(dmis, dmis)
        @test added(d) == removed(d) == Set()
        @test changed(d) == Dict()

        d = DeepDiffs.deepdiff(dnank, dmisk)
        @test d == d

    end
end

max_stars_repo_name: ssfrr/StructureDiffs

max_issues_repo_name: ssfrr/StructureDiffs

max_forks_repo_name: ssfrr/StructureDiffs


text:
# High level API exports

export Devices, dB, gainrange

## KWArgs

abstract type KWArgs <: AbstractDict{Symbol, String}; end

mutable struct OwnedKWArgs <: KWArgs
    ptr::SoapySDRKwargs
    function OwnedKWArgs(kw::SoapySDRKwargs)
        this = new(kw)
        finalizer(SoapySDRKwargs_clear, this)
        this
    end
end
Base.unsafe_load(o::OwnedKWArgs) = o.ptr
function ptr(kw::OwnedKWArgs)
    return pointer_from_objref(kw)
end

mutable struct KWArgsList <: AbstractVector{KWArgs}
    ptr::Ptr{SoapySDRKwargs}
    length::Csize_t
    function KWArgsList(ptr::Ptr{SoapySDRKwargs}, length::Csize_t)
        this = new(ptr, length)
        finalizer(SoapySDRKwargsList_clear, this)
        this
    end
end
Base.size(kwl::KWArgsList) = (kwl.length,)

struct KWArgsListRef <: KWArgs
    list::KWArgsList
    idx::Int
    function Base.getindex(kwl::KWArgsList, i::Integer)
        checkbounds(kwl, i)
        new(kwl, i)
    end
end

function ptr(kwl::KWArgsListRef)
    @assert kwl.list.ptr !== C_NULL
    kwl.list.ptr + (kwl.idx-1)*sizeof(SoapySDRKwargs)
end

Base.unsafe_load(kw::KWArgs) = unsafe_load(ptr(kw))
Base.length(kw::KWArgs) = unsafe_load(kw).size
function _getindex(kw::KWArgs, i::Integer)
    1 <= i <= length(kw) || throw(BoundsError(kw, i))
    @GC.preserve kw begin
        return Symbol(unsafe_string(unsafe_load(unsafe_load(kw).keys, i))) =>
               unsafe_string(unsafe_load(unsafe_load(kw).vals, i))
    end
end

Base.iterate(kw::KWArgs, i=1) = i > length(kw) ? nothing : (_getindex(kw, i), i+1)

function SoapySDRKwargs_clear(kw::OwnedKWArgs)
    SoapySDRKwargs_clear(pointer(kw))
end

function SoapySDRKwargsList_clear(list::KWArgsList)
    SoapySDRKwargsList_clear(list.ptr, list.length)
    list.ptr = C_NULL
end

##

mutable struct StringList <: AbstractVector{String}
    strs::Ptr{Cstring}
    length::Csize_t
    function StringList(strs::Ptr{Cstring}, length::Integer)
        this = new(strs, Csize_t(length))
        finalizer(SoapySDRStrings_clear, this)
        this
    end
end
Base.size(s::StringList) = (s.length,)
function Base.getindex(s::StringList, i::Integer)
    checkbounds(s, i)
    unsafe_string(unsafe_load(s.strs, i))
end

SoapySDRStrings_clear(s::StringList) = @GC.preserve s SoapySDRStrings_clear(pointer_from_objref(s), s.length)
##

struct Devices
    kwargslist::KWArgsList
    Devices() = new(KWArgsList(SoapySDRDevice_enumerate()...))
end
Base.length(d::Devices) = length(d.kwargslist)

function Base.show(io::IO, d::Devices)
    if length(d) == 0
        print(io, "< No devices available >")
    end
    for (i, dev) in enumerate(d.kwargslist)
        print(io, "[$i] ")
        join(io, dev, ", ")
        println(io)
    end
end

##

mutable struct Device
    ptr::Ptr{SoapySDRDevice}
    function Device(ptr::Ptr{SoapySDRDevice})
        this = new(ptr)
        finalizer(SoapySDRDevice_unmake, this)
        return this
    end
end
SoapySDRDevice_unmake(d::Device) = SoapySDRDevice_unmake(d.ptr)
Base.cconvert(::Type{<:Ptr{SoapySDRDevice}}, d::Device) = d
Base.unsafe_convert(::Type{<:Ptr{SoapySDRDevice}}, d::Device) = d.ptr

function Base.show(io::IO, d::Device)
    print(io, "SoapySDR ", d.hardware, " device (driver: ", d.driver, ") w/ ", length(d.tx), " TX channels and ", length(d.rx), " RX channels")
end

function Base.getindex(d::Devices, i::Integer)
    Device(SoapySDRDevice_make(ptr(d.kwargslist[i])))
end
Base.iterate(d::Devices, state=1) = state > length(d) ? nothing : (d[state], state+1)

function Base.getproperty(d::Device, s::Symbol)
    if s === :info
        OwnedKWArgs(SoapySDRDevice_getHardwareInfo(d))
    elseif s === :driver
        Symbol(unsafe_string(SoapySDRDevice_getDriverKey(d)))
    elseif s === :hardware
        Symbol(unsafe_string(SoapySDRDevice_getHardwareKey(d)))
    elseif s === :tx
        ChannelList(d, Tx)
    elseif s === :rx
        ChannelList(d, Rx)
    else
        getfield(d, s)
    end
end

##

struct Channel
    device::Device
    direction::Direction
    idx::Int
end
Base.show(io::IO, c::Channel) =
    print(io, "Channel(", c.device.hardware, ", ", c.direction, ", ", c.idx, ")")

# Express everything in (kHz, MHz, GHz)
function pick_freq_unit(val::Quantity)
    iszero(val.val) && return val
    abs(val) >= 1.0u"GHz" ? uconvert(u"GHz", val) :
    abs(val) >= 1.0u"MHz" ? uconvert(u"MHz", val) :
                            uconvert(u"kHz", val)
end

# Print to 3 digits of precision, no decimal point
print_3_digit(io::IO, val::Quantity) =  print(io, Base.Ryu.writeshortest(round(val.val, sigdigits=3),
    #= plus =# false,
    #= space =# false,
    #= hash =# false,
    ))

function print_unit(io::IO, val::Quantity)
    print_3_digit(io, val)
    print(io, " ", unit(val))
end

function print_unit_interval(io::IO, min, max)
    if unit(min) == unit(max) || iszero(min.val)
        print_3_digit(io, min)
        print(io, "..")
        print_3_digit(io, max)
        print(io, " ", unit(max))
    else
        print_unit(io, min)
        print(io, " .. ")
        print_unit(io, max)
    end
end
print_unit_interval(io::IO, x::Interval{<:Any, Closed, Closed}) =
    print_unit_interval(io, minimum(x), maximum(x))

using Intervals: Closed
function print_hz_range(io::IO, x::Interval{<:Any, Closed, Closed})
    min, max = pick_freq_unit(minimum(x)), pick_freq_unit(maximum(x))
    print_unit_interval(io, min, max)
end

function Base.show(io::IO, ::MIME"text/plain", c::Channel)
    println(io, c.direction == Tx ? "TX" : "RX", " Channel #", c.idx + 1, " on ", c.device.hardware)
    if !get(io, :compact, false)
        print(io, "  Selected Antenna [")
            join(io, AntennaList(c), ", ")
            println(io, "]: ", c.antenna)
        print(io, "  Bandwidth [ ")
            join(io, map(x->sprint(print_hz_range, x), bandwidth_ranges(c)), ", ")
            println(io, " ]: ", pick_freq_unit(c.bandwidth))
        print(io, "  Frequency [ ")
                join(io, map(x->sprint(print_hz_range, x), frequency_ranges(c)), ", ")
                println(io, " ]: ", pick_freq_unit(c.frequency))
            for element in FrequencyComponentList(c)
                print(io, "    ", element, " [ ")
                join(io, map(x->sprint(print_hz_range, x), frequency_ranges(c, element)), ", ")
                println(io, " ]: ", pick_freq_unit(c[element]))
            end
        println(io, "  Gain ", gainrange(c), ": ", c.gain)
            for element in GainElementList(c)
                println(io, "    ", element, " ", gainrange(c,element), ": ", c[element])
            end
        print(io, "  Sample Rate [ ", )
            join(io, map(x->sprint(print_hz_range, x), sample_rate_ranges(c)), ", ")
            println(io, " ]: ", pick_freq_unit(c.sample_rate))
        c.dc_offset_mode !== missing && println(io, "  Automatic DC offset correction: ", c.dc_offset_mode ? "ON" : "OFF")
        c.dc_offset !== missing && println(io, "  DC offset correction: ", c.dc_offset)
        c.iq_balance_mode !== missing && println(io, "  Automatic IQ balance correction: ", c.iq_balance_mode ? "ON" : "OFF")
        c.iq_balance !== missing && println(io, "  IQ balance correction: ", c.iq_balance)
        c.frequency_correction !== missing && println(io, "  Frequency correction: ", c.frequency_correction, " ppm")
    end
end

struct ChannelList <: AbstractVector{Channel}
    device::Device
    direction::Direction
end

function Base.size(cl::ChannelList)
    (SoapySDRDevice_getNumChannels(cl.device, cl.direction),)
end

function Base.getindex(cl::ChannelList, i::Integer)
    checkbounds(cl, i)
    Channel(cl.device, cl.direction, i-1)
end

function Base.getproperty(c::Channel, s::Symbol)
    if s === :info
        return OwnedKWArgs(SoapySDRDevice_getChannelInfo(c.device.ptr, c.direction, c.idx))
    elseif s === :antenna
        return Antenna(Symbol(unsafe_string(SoapySDRDevice_getAntenna(c.device.ptr, c.direction, c.idx))))
    elseif s === :gain
        return SoapySDRDevice_getGain(c.device.ptr, c.direction, c.idx)*dB
    elseif s === :dc_offset_mode
        if !SoapySDRDevice_hasDCOffsetMode(c.device.ptr, c.direction, c.idx)
            return missing
        end
        return SoapySDRDevice_getDCOffsetMode(c.device.ptr, c.direction, c.idx)
    elseif s === :dc_offset
        if !SoapySDRDevice_hasDCOffset(c.device.ptr, c.direction, c.idx)
            return missing
        end
        return SoapySDRDevice_getDCOffset(c.device.ptr, c.direction, c.idx)
    elseif s === :iq_balance_mode
        if !SoapySDRDevice_hasIQBalanceMode(c.device.ptr, c.direction, c.idx)
            return missing
        end
        return SoapySDRDevice_getIQBalanceMode(c.device.ptr, c.direction, c.idx)
    elseif s === :iq_balance
        if !SoapySDRDevice_hasIQBalance(c.device.ptr, c.direction, c.idx)
            return missing
        end
        return SoapySDRDevice_getIQBalance(c.device.ptr, c.direction, c.idx)
    elseif s === :gain_mode
        if !SoapySDRDevice_hasGainMode(c.device.ptr, c.direction, c.idx)
            return missing
        end
        return SoapySDRDevice_getGainMode(c.device.ptr, c.direction, c.idx)
    elseif s === :frequency_correction
        if !SoapySDRDevice_hasFrequencyCorrection(c.device.ptr, c.direction, c.idx)
            return missing
        end
        # TODO: ppm unit?
        return SoapySDRDevice_getFrequencyCorrection(c.device.ptr, c.direction, c.idx)
    elseif s === :sample_rate
        return SoapySDRDevice_getSampleRate(c.device.ptr, c.direction, c.idx) * Hz
    elseif s === :bandwidth
        return SoapySDRDevice_getBandwidth(c.device.ptr, c.direction, c.idx) * Hz
    elseif s === :frequency
        return SoapySDRDevice_getFrequency(c.device.ptr, c.direction, c.idx) * Hz
    else
        return getfield(c, s)
    end
end

function Base.setproperty!(c::Channel, s::Symbol, v)
    if s === :frequency
        if isa(v, Quantity)
            SoapySDRDevice_setFrequency(c.device.ptr, c.direction, c.idx, uconvert(u"Hz", v).val, C_NULL)
        elseif isa(v, FreqSpec)
            error("TODO")
        else
            error("Frequency must be specified as either a Quantity or a FreqSpec!")
        end
    elseif s === :bandwidth
        SoapySDRDevice_setBandwidth(c.device.ptr, c.direction, c.idx, uconvert(u"Hz", v).val)
    elseif s === :gain
        SoapySDRDevice_setGain(c.device.ptr, c.direction, c.idx, uconvert(u"dB", v).val)
    elseif s === :sample_rate
        SoapySDRDevice_setSampleRate(c.device.ptr, c.direction, c.idx, uconvert(u"Hz", v).val)
    else
        return setfield!(c, s, v)
    end
end

## Antenna/GainElement/FrequencyComponent

abstract type Component; end
Base.print(io::IO, c::Component) = print(io, c.name)

struct ComponentList{T<:Component} <: AbstractVector{T}
    s::StringList
end
Base.size(l::ComponentList) = (length(l.s),)
Base.getindex(l::ComponentList{T}, i::Integer) where {T} = T(Symbol(l.s[i]))

for (e, f) in zip((:Antenna, :GainElement, :FrequencyComponent),
               (:SoapySDRDevice_listAntennas, :SoapySDRDevice_listGains, :SoapySDRDevice_listFrequencies))
    @eval begin
        struct $e <: Component; name::Symbol; end
        $f(channel::Channel) =
            $f(channel.device, channel.direction, channel.idx)
        $(Symbol(string(e, "List")))(channel::Channel) = ComponentList{$e}(StringList($f(channel)...))
    end
end

using Base: unsafe_convert
function Base.getindex(c::Channel, ge::GainElement)
    SoapySDRDevice_getGainElement(c.device, c.direction, c.idx, Cstring(unsafe_convert(Ptr{UInt8}, ge.name))) * dB
end

function Base.getindex(c::Channel, fe::FrequencyComponent)
    SoapySDRDevice_getFrequencyComponent(c.device, c.direction, c.idx, Cstring(unsafe_convert(Ptr{UInt8}, fe.name))) * Hz
end

## GainElement

function _gainrange(soapyr::SoapySDRRange)
    if soapyr.step == 0.0
        # Represents an interval rather than a range
        return (soapyr.minimum*dB)..(soapyr.maximum*dB)
    end
    return range(soapyr.minimum*dB; stop=soapyr.maximum*dB, step=soapyr.step*dB)
end

function gainrange(c::Channel)
    return _gainrange(SoapySDRDevice_getGainRange(c.device, c.direction, c.idx))
end

gainrange(c::Channel, ge::GainElement) =
    return _gainrange(SoapySDRDevice_getGainElementRange(c.device, c.direction, c.idx, Cstring(unsafe_convert(Ptr{UInt8}, ge.name))))

function Base.setindex!(c::Channel, gain::typeof(1.0dB), ge::GainElement)
    SoapySDRDevice_setGainElement(c.device, c.direction, c.idx, Cstring(unsafe_convert(Ptr{UInt8}, ge.name, gain.val)))
    return gain
end

function _hzrange(soapyr::SoapySDRRange)
    if soapyr.step == 0.0
        # Represents an interval rather than a range
        return (soapyr.minimum*Hz)..(soapyr.maximum*Hz)
    end
    return range(soapyr.minimum*Hz; stop=soapyr.maximum*Hz, step=soapyr.step*Hz)
end

function bandwidth_ranges(c::Channel)
    (ptr, len) = SoapySDRDevice_getBandwidthRange(c.device.ptr, c.direction, c.idx)
    arr = map(_hzrange, unsafe_wrap(Array, Ptr{SoapySDRRange}(ptr), (len,)))
    SoapySDR_free(ptr)
    arr
end

function frequency_ranges(c::Channel)
    (ptr, len) = SoapySDRDevice_getFrequencyRange(c.device.ptr, c.direction, c.idx)
    arr = map(_hzrange, unsafe_wrap(Array, Ptr{SoapySDRRange}(ptr), (len,)))
    SoapySDR_free(ptr)
    arr
end

function frequency_ranges(c::Channel, fe::FrequencyComponent)
    (ptr, len) = SoapySDRDevice_getFrequencyRangeComponent(c.device.ptr, c.direction, c.idx, Cstring(unsafe_convert(Ptr{UInt8}, fe.name)))
    arr = map(_hzrange, unsafe_wrap(Array, Ptr{SoapySDRRange}(ptr), (len,)))
    SoapySDR_free(ptr)
    arr
end

function sample_rate_ranges(c::Channel)
    (ptr, len) = SoapySDRDevice_getSampleRateRange(c.device.ptr, c.direction, c.idx)
    arr = map(_hzrange, unsafe_wrap(Array, Ptr{SoapySDRRange}(ptr), (len,)))
    SoapySDR_free(ptr)
    arr
end

### Frequency Setting
struct FreqSpec{T}
    val::T
    kwargs::Dict{Any, String}
end

### Streams

struct StreamFormat
    T
end
function Base.print(io::IO, sf::StreamFormat)
    T = sf.T
    if T <: Complex
        print(io, 'C')
        T = real(T)
    end
    if T <: AbstractFloat
        print(io, 'F')
    elseif T <: Signed
        print(io, 'S')
    elseif T <: Unsigned
        print(io, 'U')
    else
        error("Unknown format")
    end
    print(io, 8sizeof(T))
end

function StreamFormat(s::String)
    complex = first(s) === 'C'
    complex && (s = s[2:end])
    t = first(s)
    nbits = parse(Int, s[2:end])
    if t == 'F'
        T = Dict(32 => Float32, 64 => Float64)[nbits]
    elseif t === 'S' || t === 'U'
        T = Dict(8 => UInt8, 16 => UInt16, 32 => UInt32, 64 => UInt64)[nbits]
        T === 'S' && (T = signed(T))
    end
    return StreamFormat{complex ? Complex{T} : T}()
end

function stream_formats(c::Channel)
    slist = StringList(SoapySDRDevice_getStreamFormats(c.device.ptr, c.direction, c.idx)...)
    map(StreamFormat, slist)
end

## Stream

mutable struct Stream{T}
    d::Device
    nchannels::Int
    ptr::Ptr{SoapySDRStream}
    function Stream{T}(d::Device, nchannels::Int, ptr::Ptr{SoapySDRStream}) where {T}
        this = new{T}(d, nchannels, ptr)
        finalizer(SoapySDRDevice_closeStream, this)
        return this
    end
end
Base.cconvert(::Type{<:Ptr{SoapySDRStream}}, s::Stream) = s
Base.unsafe_convert(::Type{<:Ptr{SoapySDRStream}}, s::Stream) = s.ptr
SoapySDRDevice_closeStream(s::Stream) = SoapySDRDevice_closeStream(s.d, s)

function Base.show(io::IO, s::Stream)
    print(io, "Stream on ", s.d.hardware)
end

function Stream(format::Union{StreamFormat, Type}, device::Device, direction::Direction; kwargs...)
    format = StreamFormat(format)
    isempty(kwargs) || error("TODO")
    Stream{T}(device, 1, SoapySDRDevice_setupStream(device, direction, string(format), C_NULL, 0, C_NULL))
end

function Stream(format::Union{StreamFormat, Type}, channels::Vector{Channel}; kwargs...)
    format = StreamFormat(format)
    isempty(kwargs) || error("TODO")
    isempty(channels) && error("Must specify at least one channel or use the device/direction constructor for automatic.")
    device = first(channels).device
    direction = first(channels).direction
    if !all(channels) do channel
                channel.device == device && channel.direction == direction
            end
        throw(ArgumentError("Channels must agree on device and direction"))
    end
    Stream{format.T}(device, length(channels), SoapySDRDevice_setupStream(device, direction, string(format), map(x->x.idx, channels), length(channels), C_NULL))
end

function _read!(s::Stream{T}, buffers::NTuple{N, Vector{T}}; timeout=nothing) where {N, T}
    timeout === nothing && (timeout = 0.1u"s") # Default from SoapySDR upstream
    buflen = length(first(buffers))
    @assert all(buffer->length(buffer) == buflen, buffers)
    @assert N == s.nchannels
    n, flags, timens = SoapySDRDevice_readStream(s.d, s, Ref(map(pointer, buffers)), buflen, uconvert(u"μs", timeout).val)
    timens = timens * u"ns"
    n, flags, timens
end

function Base.read!(s, buffers; kwargs...)
    _read!(s, buffers; kwargs...)[1]
end

struct SampleBuffer{N, T}
    bufs::NTuple{N, Vector{T}}
    flags::Cint
    timens::typeof(1u"ns")
end
Base.length(sb::SampleBuffer) = length(sb.bufs[1])

function Base.read(s::Stream{T}, n::Int; kwargs...) where {T}
    bufs = ntuple(_->Vector{T}(undef, n), s.nchannels)
    nread, flags, timens = _read!(s, bufs; kwargs...)
    if nread != n
        @info("assertion debugging", nread, n)
        @assert nread == n
    end
    SampleBuffer(bufs, flags, timens)
end

function activate!(s::Stream; flags = 0, timens = nothing, numElems=0)
    SoapySDRDevice_activateStream(s.d, s, flags, timens === nothing ? 0 : uconvert(u"ns", timens).val, numElems)
    nothing
end

function deactivate!(s::Stream; flags = 0, timens = nothing)
    SoapySDRDevice_deactivateStream(s.d, s, flags, timens === nothing ? 0 : uconvert(u"ns", timens).val)
    nothing
end

function Base.write(s::Stream{T}, buffers::NTuple{N, Vector{T}}; timeout = nothing) where {N, T}
    timeout === nothing && (timeout = 0.1u"s") # Default from SoapySDR upstream
    buflen = length(first(buffers))
    @assert all(buffer->length(buffer) == buflen, buffers)
    @assert N == s.nchannels
    SoapySDRDevice_writeStream(s.d, s, Ref(map(pointer, buffers)), buflen, 0, 0, uconvert(u"μs", timeout).val)
end


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

mutable struct OutboundNatRulePropertiesFormat <: SwaggerModel
    allocatedOutboundPorts::Any # spec type: Union{ Nothing, Int32 } # spec name: allocatedOutboundPorts
    frontendIPConfigurations::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: frontendIPConfigurations
    backendAddressPool::Any # spec type: Union{ Nothing, SubResource } # spec name: backendAddressPool
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function OutboundNatRulePropertiesFormat(;allocatedOutboundPorts=nothing, frontendIPConfigurations=nothing, backendAddressPool=nothing, provisioningState=nothing)
        o = new()
        validate_property(OutboundNatRulePropertiesFormat, Symbol("allocatedOutboundPorts"), allocatedOutboundPorts)
        setfield!(o, Symbol("allocatedOutboundPorts"), allocatedOutboundPorts)
        validate_property(OutboundNatRulePropertiesFormat, Symbol("frontendIPConfigurations"), frontendIPConfigurations)
        setfield!(o, Symbol("frontendIPConfigurations"), frontendIPConfigurations)
        validate_property(OutboundNatRulePropertiesFormat, Symbol("backendAddressPool"), backendAddressPool)
        setfield!(o, Symbol("backendAddressPool"), backendAddressPool)
        validate_property(OutboundNatRulePropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type OutboundNatRulePropertiesFormat

const _property_map_OutboundNatRulePropertiesFormat = Dict{Symbol,Symbol}(Symbol("allocatedOutboundPorts")=>Symbol("allocatedOutboundPorts"), Symbol("frontendIPConfigurations")=>Symbol("frontendIPConfigurations"), Symbol("backendAddressPool")=>Symbol("backendAddressPool"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_OutboundNatRulePropertiesFormat = Dict{Symbol,String}(Symbol("allocatedOutboundPorts")=>"Int32", Symbol("frontendIPConfigurations")=>"Vector{SubResource}", Symbol("backendAddressPool")=>"SubResource", Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ OutboundNatRulePropertiesFormat }) = collect(keys(_property_map_OutboundNatRulePropertiesFormat))
Swagger.property_type(::Type{ OutboundNatRulePropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_OutboundNatRulePropertiesFormat[name]))}
Swagger.field_name(::Type{ OutboundNatRulePropertiesFormat }, property_name::Symbol) =  _property_map_OutboundNatRulePropertiesFormat[property_name]

function check_required(o::OutboundNatRulePropertiesFormat)
    (getproperty(o, Symbol("backendAddressPool")) === nothing) && (return false)
    true
end

function validate_property(::Type{ OutboundNatRulePropertiesFormat }, name::Symbol, val)
end


text:
using Adapt

struct MultipleForcings{N, F}
    forcings :: F
end

Adapt.adapt_structure(to, mf::MultipleForcings) = MultipleForcings(adapt(to, mf.forcings))

Base.getindex(mf::MultipleForcings, i) = mf.forcings[i]

"""
    MultipleForcings(forcings)

Return a lightweight tuple-wrapper representing multiple user-defined `forcings`.
Each forcing in `forcings` is added to the specified field's tendency.
"""
function MultipleForcings(forcings)
    N = length(forcings)
    F = typeof(forcings)
    return MultipleForcings{N, F}(forcings)
end

function regularize_forcing(forcing_tuple::Tuple, field, field_name, model_field_names)
    forcings = Tuple(regularize_forcing(f, field, field_name, model_field_names)
                     for f in forcing_tuple)
    return MultipleForcings(forcings)
end

@inline (mf::MultipleForcings{1})(i, j, k, grid, clock, model_fields) = mf.forcings[1](i, j, k, grid, clock, model_fields)

@inline (mf::MultipleForcings{2})(i, j, k, grid, clock, model_fields) = mf.forcings[1](i, j, k, grid, clock, model_fields) +
                                                                        mf.forcings[2](i, j, k, grid, clock, model_fields)

@inline (mf::MultipleForcings{3})(i, j, k, grid, clock, model_fields) = mf.forcings[1](i, j, k, grid, clock, model_fields) +
                                                                        mf.forcings[2](i, j, k, grid, clock, model_fields) +
                                                                        mf.forcings[3](i, j, k, grid, clock, model_fields)

@inline (mf::MultipleForcings{4})(i, j, k, grid, clock, model_fields) = mf.forcings[1](i, j, k, grid, clock, model_fields) +
                                                                        mf.forcings[2](i, j, k, grid, clock, model_fields) +
                                                                        mf.forcings[3](i, j, k, grid, clock, model_fields) +
                                                                        mf.forcings[4](i, j, k, grid, clock, model_fields)

# The magic (which doesn't seem to work on GPU now)
@inline function (mf::MultipleForcings{N})(i, j, k, grid, clock, model_fields) where N
    total_forcing = zero(eltype(grid))
    forcings = mf.forcings
    ntuple(Val(N)) do n
        Base.@_inline_meta
        @inbounds begin
            nth_forcing = forcings[n]
            total_forcing += nth_forcing(i, j, k, grid, clock, model_fields)
        end
    end
    return total_forcing
end

Base.summary(mf::MultipleForcings) = string("MultipleForcings with ", length(mf.forcings), " forcing",
                                            ifelse(length(mf.forcings) > 1, "s", ""))

function Base.show(io::IO, mf::MultipleForcings)
    start = summary(mf) * ":"

    Nf = length(mf.forcings)
    if Nf > 1
        body = [string("├ ", prettysummary(f), '\n') for f in mf.forcings[1:end-1]]
    else
        body = []
    end

    push!(body, string("└ ", prettysummary(mf.forcings[end])))

    print(io, start, '\n', body...)

    return nothing
end


text:
@testset "normal distribution, integral normed to 1" begin
    @test 1 / √(2π)^2 / 1 * exp.(.- rr2((5,5)) / 2 .* 1^2) ≈ normal((5,5))
    @test 1 / √(2π)^2 / 1 * exp.(.- rr2((8,5)) / 2 .* 1^2) ≈ normal((8,5))
    @test 1 / √(2π)^2 / 1 * exp.(.- rr2((8,5)) / 2 .* 1^2) ≈ normal(zeros(8,5))
    @test eltype(normal(Float32, (8,5))) === Float32
end

@testset "Gaussian function, normed to 1" begin
    arr = 1 / √(2π)^2 / 1 * exp.(.- rr2((5,5)) / 2 .* 1^2)
    arr ./= maximum(arr)
    @test arr ≈ gaussian((5,5))
    @test eltype(gaussian(Float32, (8,5))) === Float32
end


text:
is_nan(cur_value::Number) = isnan(cur_value)

function is_nan(cur_value::SymEngine.Basic)
  isempty(free_symbols(cur_value)) || return false

  isnan(N(cur_value))
end

is_finite(cur_value::Number) = isfinite(cur_value)

function is_finite(cur_value::SymEngine.Basic)
  isempty(free_symbols(cur_value)) || return true
  is_nan(cur_value) && return false

  isfinite(N(cur_value))
end


text:
"""
json2volume(path::String)

Read a file `.json` of volume model.

# Example of a volume file json structure.

```
{
   "clip":true,
   "name":"name",
   "scale":{
      "x":1.,
      "y":1.,
      "z":1.
   },
   "position":{
   	  "x":0.,
	  "y":0.,
	  "z":0.
   },
   "rotation":{
      "x":0.,
      "y":0.,
      "z":0.
   },
   "permitExtraction":true
}
```
"""
function json2volume(path::String)
	dict = Dict{String,Any}[]

	open(path, "r") do f
	    dict = JSON.parse(f)  # parse and transform data
	end

	position = dict["position"]
	scale = dict["scale"]
	rotation = dict["rotation"]

	return Volume(  [scale["x"],scale["y"],scale["z"]],
					[position["x"],position["y"],position["z"]],
					[rotation["x"],rotation["y"],rotation["z"]]
				)
end

"""
	json2LARvolume(path::String)

Return LAR model of Potree volume tools.
"""
function json2LARvolume(path::String)
	@assert isfile(path) "json2LARvolume: $path not an existing file"

	volume = json2volume(path)
	return Common.volume2LARmodel(volume::Volume)
end

"""
	json2ucs(file::String)

Read a file `.json` of UCS.

# Example of a UCS file json structure.

```
{
	"id":"270b6b7a-d00c-46f6-ba58-1c76310558aa",
	"data":{
		"plane":{
			"A":0,
			"B":0,
			"C":0
		},
		"xAxis":{
			"x":0.21194956712662227,
			"y":-0.9771563332378362,
			"z":-0.015584652964510243
		},
		"yAxis":{
			"x":0.007193328934240251,
			"y":-0.014386657868480502,
			"z":0.9998706316790285
		},
		"zAxis":{
			"x":-0.9772541312338778,
			"y":-0.21203425310209192,
			"z":0.003979761017532582
		},
		"origin":{
			"x":2.250469923019409,
			"y":-6.521675497293472,
			"z":-1.4895449578762054
		}
	},
	"name":"XY_PIANO_FINESTRA"

}
```
"""
function json2ucs(file::String)
	dict = Dict{String,Any}[]

	open(file, "r") do f
	    dict = JSON.parse(f)  # parse and transform data
	end

	return dict
end

function ucs2matrix(file::String)
	dict = json2ucs(file)

	origin = dict["data"]["origin"]
	xAxis = dict["data"]["xAxis"]
	yAxis = dict["data"]["yAxis"]
	zAxis = dict["data"]["zAxis"]
	M =  zeros(4,4)

	M[1,1] = xAxis["x"]
	M[1,2] = xAxis["y"]
	M[1,3] = xAxis["z"]

	M[2,1] = yAxis["x"]
	M[2,2] = yAxis["y"]
	M[2,3] = yAxis["z"]

	M[3,1] = zAxis["x"]
	M[3,2] = zAxis["y"]
	M[3,3] = zAxis["z"]

	O = M[1:3, 1:3] * -[origin["x"], origin["y"], origin["z"]]
	M[1,4] = O[1]
	M[2,4] = O[2]
	M[3,4] = O[3]

	M[4,4] = 1.0
	return M
end

"""
measure JSON from Potree 1.7
"""
function seedPointsFromFile(path::String)
	dataset = Matrix[]
	dict=Dict{String,Any}[]
	open(path, "r") do f
	    dict = JSON.parse(f)  # parse and transform data
	end

	objects = dict["features"]

	for obj in objects
		geometry = obj["geometry"]
		if geometry["type"] == "Polygon"
			coords = []
			vec_coords = geometry["coordinates"][1]
			for element in vec_coords
				push!(coords,element)
			end
			mycoords = hcat(coords...)
			push!(dataset,mycoords)
		end
	end
	return dataset
end


text:
using .JSort

function sortfile(parameters::Parameters)
    labr = Vector{LaBrEvent}
    for file in parameters.datafiles
        println("Processing file $file")
        @time labr = processfile(file, parameters)
    end
    println("Saving.")
    @show length(labr)
    @time save(labr, parameters.savepath)
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
using ConvexTests, SCS

@info "Starting SCS tests"

do_tests("SCS", () -> SCS.Optimizer(verbose=0, eps=1e-6); exclude = [r"mip"], description = "Tests run with `eps=1e-6`.")


text:
using WignerFunctions
using Test
using Random

d222(β::AbstractFloat) = (1 + cos(β))^2 / 4
d221(β::AbstractFloat) = -sin(β) * (1 + cos(β)) / 2
d220(β::AbstractFloat) = sin(β) * sin(β) * √(3 / 8)
d221̄(β::AbstractFloat) = -sin(β) * (1 - cos(β)) / 2
d222̄(β::AbstractFloat) = (1 - cos(β))^2 / 4
d211(β::AbstractFloat) = (2cos(β) * cos(β) + cos(β) - 1) / 2
d210(β::AbstractFloat) = -sin(2β) * √(3 / 8)
d211̄(β::AbstractFloat) = (-2cos(β) * cos(β) + cos(β) + 1) / 2
d200(β::AbstractFloat) = (3cos(β) * cos(β) - 1) / 2
WignerFunctions.Index(::typeof(d222)) = WignerFunctions.Index((2, 2, 2))
WignerFunctions.Index(::typeof(d221)) = WignerFunctions.Index((2, 2, 1))
WignerFunctions.Index(::typeof(d220)) = WignerFunctions.Index((2, 2, 0))
WignerFunctions.Index(::typeof(d221̄)) = WignerFunctions.Index((2, 2, -1))
WignerFunctions.Index(::typeof(d222̄)) = WignerFunctions.Index((2, 2, -2))
WignerFunctions.Index(::typeof(d211)) = WignerFunctions.Index((2, 1, 1))
WignerFunctions.Index(::typeof(d210)) = WignerFunctions.Index((2, 1, 0))
WignerFunctions.Index(::typeof(d211̄)) = WignerFunctions.Index((2, 1, -1))
WignerFunctions.Index(::typeof(d200)) = WignerFunctions.Index((2, 0, 0))

@testset "Naive" begin
    naive = WignerFunctions.naive
    @testset "diagonal: β=-0, l=$l, m₁=$m₁, m₂=$m₂" for l ∈ 0:4, m₁ ∈ -l:l, m₂ ∈ -l:l
        @test naive(BigFloat(0), l, m₁, m₂) ≈ (m₁ == m₂ ? 1 : 0)
    end

    functions = [d222, d221, d220, d221̄, d222̄, d211, d210, d211̄, d211̄, d200]
    @testset "function $func" for func ∈ functions
        β = 2π*rand(10)
        l, m, m′ = WignerFunctions.Index(func)
        @test all(func.(β) .≈ naive.(BigFloat.(β), l, m, m′))
        @test all(func.(β) .≈ ((m - m′) % 2 == 0 ? 1 : -1) .* naive.(BigFloat.(β),  l, m′, m))
        @test all(naive.(BigFloat.(β), l, m′, m) .≈ naive.(BigFloat.(β),  l, -m, -m′))
    end
end

@testset "$method" for method ∈ (WignerFunctions.Risbo.nosym,)
    dmatrix = method(0)
    @testset "diagonal: β=-0, l=$l, m₁=$m₁, m₂=$m₂" for l ∈ 0:4, m₁ ∈ -l:l, m₂ ∈ -l:l
        @test dmatrix(l, m₁, m₂) ≈ (m₁ == m₂ ? 1 : 0)
    end

    l = rand(2:20, 50)
    m = rand.(range.(-l, l, step=1))
    m′ = rand.(range.(-l, l, step=1))
    β = 2π * rand(length(l))
    @testset "Symmetries" begin
        @testset "β=$(β[i]), l=$(l[i]), m=$(m[i]), m′=$(m′[i]) " for i ∈ 1:length(l)
            k, m₁, m₂ = l[i], m[i], m′[i]
            dmatrix = method(β[i])
            @test dmatrix(k, m₁, m₂) ≈ ((m₁ - m₂) % 2 == 0 ? 1 : -1) * dmatrix(k, m₂, m₁)
            @test dmatrix(k, m₁, m₂) ≈ ((m₁ - m₂) % 2 == 0 ? 1 : -1) * dmatrix(k, -m₁, -m₂)
            @test dmatrix(k, m₁, m₂) ≈ dmatrix(k, -m₂, -m₁)
        end
    end

    naive = WignerFunctions.naive
    oldprec = precision(BigFloat)
    setprecision(512)
    @testset "Against naive β=$(β[i]), l=$(l[i]), m=$(m[i]), m′=$(m′[i]) " for i ∈ 1:length(l)
        @test method(β[i])(l[i], m[i], m′[i]) ≈ float(naive(BigFloat(β[i]), l[i], m[i], m′[i]))
    end
    setprecision(oldprec)
end
nothing

@testset "Delta recurrence" begin
    recurrence = WignerFunctions.Trapani.delta_recurrence
    @test recurrence(0, 0, 0) == 1
    @test recurrence(3, 4, 0) == 0
    @test recurrence(3, 0, -4) == 0

    ls = rand(2:20, 50)
    ms = rand.(range.(-ls, ls, step=1))
    ms′ = rand.(range.(-ls, ls, step=1))
    oldprec = precision(BigFloat)
    setprecision(512)
    naive = (l, m₁, m₂) -> float(WignerFunctions.naive(BigFloat(pi)/2, l, m₁, m₂))
    @testset "Against naive l=$(ls[i]), m=$(ms[i]), m′=$(ms′[i]) " for i ∈ 1:length(ls)
        l, m₁, m₂ = ls[i], ms[i], ms′[i]
        @test recurrence(l, m₁, m₂) ≈ naive(l, m₁, m₂)
    end
    setprecision(oldprec)
end


text:
## -------------------------------------------------------------------
# krebs per ider per exp per src
let

    for (iJR, Data) in [
        (ChK.iJR904, ChK.KayserData),
        (ChN.iJR904, ChN.NanchenData),
        (ChF.iJR904, ChF.FolsomData),
        (ChH.iJR904, ChH.HeerdenData),
    ]

        src = nameof(Data)

        datfile = iJR.procdir("dat.bson")
        DAT = UJL.load_data(datfile; verbose = false)

        METHODS = [:FBA_Z_FIX_MIN_COST, :ME_MAX_POL]
        EXPS = DAT[:EXPS]
        EXCH_FLX_IDERS = iJR.load_krebs_iders()

        for exp in EXPS
            p = plot(;xlabel = "flx ider", ylabel = "flx val / exp glcex")
            glc_exp_val = abs(DAT[:exp, :flx, "GLC", exp])

            for method in METHODS
                model_vals = abs.(DAT[method, :flx, EXCH_FLX_IDERS, exp]) ./ glc_exp_val

                color = method_colors[method]
                scatter!(p, EXCH_FLX_IDERS, model_vals;
                    color, label = string(method),
                    m = 8
                )
                plot!(p, EXCH_FLX_IDERS, model_vals; label = "",
                    color, ls = :dash, alpha = 0.5, lw = 3
                )
            end

            plot!(;xrotation = 35, ylim = [0.0, Inf])
            mysavefig(p, "krebs"; src, exp)
        end
    end
end

## -------------------------------------------------------------------
# krebs avs
let

    fba_method = :FBA_Z_FIX_MIN_COST
    me_method = :ME_MAX_POL
    EXCH_FLX_IDERS = KiJR.krebs_iders

    for (iJR, Data) in [
        (ChK.iJR904, ChK.KayserData),
        (ChN.iJR904, ChN.NanchenData),
        (ChF.iJR904, ChF.FolsomData),
        (ChH.iJR904, ChH.HeerdenData),
    ]
        me_avs, me_stds = [], []
        fba_avs, fba_stds = [], []
        src = nameof(Data)
        for ider in EXCH_FLX_IDERS

            datfile = iJR.procdir("dat.bson")
            DAT = UJL.load_data(datfile; verbose = false)

            EXPS = DAT[:EXPS]

            glc_exp_vals = abs.(DAT[:exp, :flx, "GLC", EXPS])
            fba_vals = abs.(DAT[fba_method, :flx, ider, EXPS]) ./ glc_exp_vals
            me_vals = abs.(DAT[me_method, :flx, ider, EXPS]) ./ glc_exp_vals
            push!(fba_avs, mean(fba_vals)); push!(fba_stds, std(fba_vals))
            push!(me_avs, mean(me_vals)); push!(me_stds, std(me_vals))
        end

        p = plot(;title = string(src), xlabel = "flx ider", ylabel = "flx val / exp glcex")

        for (avs, stds, method) in [
            (me_avs, me_stds, me_method),
            (fba_avs, fba_stds, fba_method),
        ]
            color = method_colors[method]
            lb, ub = avs .- stds, avs .+ stds
            lb, ub = avs .- stds, avs .+ stds
            scatter!(p, EXCH_FLX_IDERS, [avs avs]; label = "",
                fillrange = [lb ub], fillalpha = 0.2, color,
                markeralpha = 0.0
            )
        end

        marker = (8, source_markers[Data])
        for (avs, method) in [
            (me_avs, me_method),
            (fba_avs, fba_method),
        ]
            color = method_colors[method]
            scatter!(p, EXCH_FLX_IDERS, avs; label = string(method),
                color, marker
            )
            plot!(p, EXCH_FLX_IDERS, avs; label = "",
                color, lw = 1, alpha = 0.3, ls = :dash
            )
        end

        plot!(p; xrotation = 35, ylim = [0.0, Inf])
        mysavefig(p, "krebs_ave_comparison"; src)
    end
end


text:
@testset "Test the Jacobian of the pressure with respect to the position for a set of point vortices" begin

	atol_AD  = 1000*eps()

	Ny = 15
	sensors = randn(ComplexF64, Ny)

	# Test for a set of vortices
	Ns = 10
	zs = randn(ComplexF64, Ns)
	Ss = randn(Ns)

	sys = Vortex.Point.(zs, Ss)
	U = 0.0*im
	freestream = Freestream(U)

	dpdzAD, dpdzstarAD = PotentialFlow.Elements.jacobian_position(x->pressure_AD(sensors, x, 0.0), sys)

	dpdz, dpdzstar = analytical_jacobian_position(sensors, sys, freestream, 0.0)

	@test isapprox(dpdz, dpdzAD, atol = atol_AD)
	@test isapprox(dpdzstar, dpdzstarAD, atol = atol_AD)
end

@testset "Test the Jacobian of the pressure with respect to the strength for a set of point vortices" begin

	atol_AD  = 1000*eps()

	Ny = 15
	sensors = randn(ComplexF64, Ny)

	# Test for a set of vortices
	Ns = 10
	zs = randn(ComplexF64, Ns)
	Ss = randn(Ns)

	sys = Vortex.Point.(zs, Ss)
	U = 0.0*im
	freestream = Freestream(U)

	dpdSAD = PotentialFlow.Elements.jacobian_strength(x->pressure_AD(sensors, x, 0.0), sys)

	dpdS, dpdSstar = analytical_jacobian_strength(sensors, sys, freestream, 0.0)

	# Here ∂p/∂ΓL = -i(∂p/∂SL - ∂p/∂S̄L) with SL = QL - iΓL, as S = Γ + iQ is used as the convention in PotentialFlow.jl
	@test isapprox(-im*(dpdS - dpdSstar), dpdSAD, atol = atol_AD)
end

@testset "Test the Jacobian of the pressure with respect to the position for a set of regularized vortices" begin

	atol_AD  = 100*eps()

	Ny = 15
	sensors = randn(ComplexF64, Ny)

	δ = 1e-1

	# Test for a set of vortices
	Ns = 10
	zs = randn(ComplexF64, Ns)
	Ss = randn(Ns)

	sys = Vortex.Blob.(zs, Ss, δ*ones(Ns))
	U = 0.0*im
	freestream = Freestream(U)

	dpdzAD, dpdzstarAD = PotentialFlow.Elements.jacobian_position(x->pressure_AD(sensors, x, 0.0), sys)

	dpdz, dpdzstar = analytical_jacobian_position(sensors, sys, freestream, 0.0)

	@test isapprox(dpdz, dpdzAD, atol = atol_AD)
	@test isapprox(dpdzstar, dpdzstarAD, atol = atol_AD)
end

@testset "Test the Jacobian of the pressure with respect to the strength for a set of regularized vortices" begin

	atol_AD  = 1000*eps()

	Ny = 15
	sensors = randn(ComplexF64, Ny)

	δ = 1e-1

	# Test for a set of vortices
	Ns = 10
	zs = randn(ComplexF64, Ns)
	Ss = randn(Ns)

	sys = Vortex.Blob.(zs, Ss, δ*ones(Ns))
	U = 0.0*im
	freestream = Freestream(U)

	dpdSAD = PotentialFlow.Elements.jacobian_strength(x->pressure_AD(sensors, x, 0.0), sys)

	dpdS, dpdSstar = analytical_jacobian_strength(sensors, sys, freestream, 0.0)

	# Here ∂p/∂ΓL = -i(∂p/∂SL - ∂p/∂S̄L) with SL = QL - iΓL, as S = Γ + iQ is used as the convention in PotentialFlow.jl
	@test isapprox(-im*(dpdS - dpdSstar), dpdSAD, atol = atol_AD)
end

@testset "Validate analytical_jacobian_pressure for a set of singularities" begin
    atol_AD = 1000*eps()
	Nv = 10
	zv = randn(Nv) + im *rand(Nv)
	Sv = randn(Nv)

	sensors = randn(ComplexF64, 15)

	sys = Vortex.Point.(zv, Sv)
	U = 0.0*im
	freestream = Freestream(U)

	∂p = analytical_jacobian_pressure(sensors, sys, freestream, 0.0)

	∂pAD = AD_symmetric_jacobian_pressure(sensors, sys, 0.0)

	@test isapprox(∂p, ∂pAD, atol = atol_AD)
end

@testset "Validate analytical_jacobian_pressure for a set of regularized singularities" begin
    atol_AD = 1000*eps()
	Nv = 10
	zv = randn(Nv) + im *rand(Nv)
	Sv = randn(Nv)
    δ = 1e-1

	sensors = randn(ComplexF64, 15)

	sys = Vortex.Blob.(zv, Sv, δ*ones(Nv))
	U = 0.0*im
	freestream = Freestream(U)

	∂p = analytical_jacobian_pressure(sensors, sys, freestream, 0.0)

	∂pAD = AD_symmetric_jacobian_pressure(sensors, sys, 0.0)

	@test isapprox(∂p, ∂pAD, atol = atol_AD)
end


text:
using Reexport

@reexport module SMeasure
export s_measure

using NearestNeighbors, Distances, Neighborhood, DelayEmbeddings

"""
    s_measure(x::AbstractDataset, y::AbstractDataset;
        K::Int = 2, metric = SqEuclidean(), tree_metric = Euclidean()) → Float64
    s_measure(x::AbstractVector, y::AbstractVector;
        K::Int = 2, metric = SqEuclidean(), tree_metric = Euclidean(),
        dx::Int = 2, my::Int = 2, τx::Int = 1, τy::Int = 1) → Float64
    s_measure(x::AbstractDataset, y::AbstractVector;
        K::Int = 2, metric = SqEuclidean(), tree_metric = Euclidean(),
        dx::Int = 2, τx::Int = 1) → Float64
    s_measure(x::AbstractDataset, y::AbstractVector;
        K::Int = 2, metric = SqEuclidean(), tree_metric = Euclidean(),
        dy::Int = 2, τy::Int = 1) → Float64

S-measure [^Grassberger1999][^Quiroga2000] for the directional dependence between
univariate/multivariate time series `x` and `y`.

Returns a scalar `s ∈ [0, 1]`, where `s = 0` indicates independence between `x` and `y`,
and higher values indicate synchronization between `x` and `y`, with complete
synchronization for `s = 1.0`.

## Input data

The algorithm is slightly modified from [1] to allow univariate time series as input.
In that case, these time series are embedded using the parameters `mx`/`τx` (dimension/lag
for time series `x`) or `my`/`τy` (dimension/lag for time series `y`). For custom embeddings,
do the embedding yourself and input `Dataset`s as both `x` and `y`.

- If `x` and `y` are `dx` and `dy`-dimensional [`Dataset`](@ref)s, respectively,
    then use `x` and `y` as is.
- If `x` and `y` are scalar time series, then create `dx` and `dy` dimensional embeddings,
    respectively, of both `x` and `y`, resulting in `N` different `m`-dimensional embedding points
    ``X = \\{x_1, x_2, \\ldots, x_N \\}`` and ``Y = \\{y_1, y_2, \\ldots, y_N \\}``.
    `τx` and `τy` control the embedding lags for `x` and `y`.
- If `x` is a scalar-valued vector and `y` is a [`Dataset`](@ref), or vice versa,
    then create an embedding of the scalar time series using parameters `dx`/`τx` or `dy`/`τy`.

In all three cases, input datasets are length-matched by eliminating points at the end of
the longest dataset (after the embedding step, if relevant) before analysis.

## Description

1. Let ``r_{i,j}`` and ``s_{i,j}`` be the indices of the `K`-th nearest neighbors
    of ``x_i`` and ``y_i``, respectively.

2. Compute the the mean squared Euclidean distance to the ``K`` nearest neighbors
    for each ``x_i``, using the indices ``r_{i, j}``.

```math
R_i^{(k)}(x) = \\dfrac{1}{k} \\sum_{i=1}^{k}(x_i, x_{r_{i,j}})^2
```

- Compute the y-conditioned mean squared Euclidean distance to the ``K`` nearest
    neighbors for each ``x_i``, now using the indices ``s_{i,j}``.

```math
R_i^{(k)}(x|y) = \\dfrac{1}{k} \\sum_{i=1}^{k}(x_i, x_{s_{i,j}})^2
```

- Define the following measure of independence, where ``0 \\leq S \\leq 1``, and
    low values indicate independence and values close to one occur for
    synchronized signals.

```math
S^{(k)}(x|y) = \\dfrac{1}{N} \\sum_{i=1}^{N} \\dfrac{R_i^{(k)}(x)}{R_i^{(k)}(x|y)}
```

## Example

```julia
using CausalityTools

# A two-dimensional Ulam lattice map
sys = ulam(2)

# Sample 1000 points after discarding 5000 transients
orbit = trajectory(sys, 1000, Ttr = 5000)
x, y = orbit[:, 1], orbit[:, 2]

# 4-dimensional embedding for `x`, 5-dimensional embedding for `y`
s_measure(x, y, dx = 4, τx = 3, dy = 5, τy = 1)
```

[^Quiroga2000]: Quian Quiroga, R., Arnhold, J. & Grassberger, P. [2000] “Learning driver-response relationships from synchronization patterns,” Phys. Rev. E61(5), 5142–5148.
[^Grassberger1999]: Arnhold, J., Grassberger, P., Lehnertz, K., & Elger, C. E. (1999). A robust method for detecting interdependences: application to intracranially recorded EEG. Physica D: Nonlinear Phenomena, 134(4), 419-430.
"""
function s_measure(x::AbstractDataset{D1, T}, y::AbstractDataset{D2, T}; K::Int = 3,
        metric = SqEuclidean(),
        tree_metric = Euclidean(),
        theiler_window::Int = 0 # only point itself excluded
        ) where {D1, D2, T}

    # Match length of datasets by excluding end points.
    lx = length(x); ly = length(y)
    lx > ly ? X = x[1:ly, :] : X = x
    ly > lx ? Y = y[1:lx, :] : Y = y
    N = length(X)

    # Pre-allocate vectors to hold indices and distances during loops
    dists_x = zeros(T, K)
    dists_x_cond_y = zeros(T, K)

    # Mean squared distances in X, and
    # mean squared distances in X conditioned on Y
    Rx = zeros(T, N)
    Rx_cond_y = zeros(T, N)

    # Search for the K nearest neighbors of each points in both X and Y
    treeX = searchstructure(KDTree, X, tree_metric)
    treeY = searchstructure(KDTree, Y, tree_metric)
    neighborhoodtype, theiler = NeighborNumber(K), Theiler(0)
    idxs_X = bulkisearch(treeX, X, neighborhoodtype, theiler)
    idxs_Y = bulkisearch(treeY, Y, neighborhoodtype, theiler)

    for n in 1:N
        pxₙ = X[n]

        for j = 1:K
            rₙⱼ = idxs_X[n][j] # nearest neighbor indices in X
            sₙⱼ = idxs_Y[n][j] # nearest neighbor indices in Y
            dists_x[j] = evaluate(metric, pxₙ, X[rₙⱼ])
            dists_x_cond_y[j] = evaluate(metric, pxₙ, X[sₙⱼ])
        end

        Rx[n] = sum(dists_x) / K
        Rx_cond_y[n] = sum(dists_x_cond_y) / K
    end

    return sum(Rx ./ Rx_cond_y) / N
end

function s_measure(x::AbstractVector{T}, y::AbstractVector{T}; K::Int = 3,
    dx::Int = 2, dy::Int = 2, τx = 1, τy = 1, metric = SqEuclidean(),
    tree_metric = Euclidean()) where T

    X = embed(x, dx, τx)
    Y = embed(y, dy, τy)
    lX, lY = length(X), length(Y)

    s_measure(
        # TODO: cut the last points of the shortest resulting embedding.
        lX > lY ? X[1:lY, :] : X,
        lY > lX ? Y[1:lX, :] : Y,
        K = K, metric = metric, tree_metric = tree_metric)
end

function s_measure(X::AbstractDataset{D}, y::AbstractVector{T}; K::Int = 3,
        dy::Int = 2, τy = 1,
        metric = SqEuclidean(), tree_metric = Euclidean()) where {D, T}

    Y = embed(y, dy, τy)
    s_measure(X[1:length(Y), :], Y, K = K, metric = metric, tree_metric = tree_metric)
end

function s_measure(x::AbstractVector{T}, Y::AbstractDataset{D}; K::Int = 3,
        dx::Int = 2, τx = 1, metric = SqEuclidean(),
        tree_metric = Euclidean()) where {D, T}
    X = embed(x, dx, τx)

    s_measure(X, Y[1:length(X), :], K = K, metric = metric, tree_metric = tree_metric)
end

end


text:
using Distributed
using DistributedArrays

DA = DistributedArrays

function prepare!(grid::DArray)
    m, n = size(grid)
    P = procs(grid)
    # Note: To simplify the code, we'll set the corners of the matrix twice.
    @sync begin
        for p in P[1,:]
            @async remotecall_fetch(p) do
                _, cols = localindices(grid)
                lgrid = localpart(grid)
                view(lgrid, 1, :) .= grid[m-1,cols]
                nothing
            end
        end
        for p in P[end,:]
            @async remotecall_fetch(p) do
                _, cols = localindices(grid)
                lgrid = localpart(grid)
                lm, ln = size(lgrid)
                view(lgrid, lm, :) .= grid[2,cols]
                nothing
            end
        end
    end
    # Now we are good to include the corners as well:
    @sync begin
        for p in P[:,1]
            @async remotecall_fetch(p) do
                rows, _ = localindices(grid)
                lgrid = localpart(grid)
                view(lgrid, :, 1) .= grid[rows,n-1]
                nothing
            end
        end
        for p in P[:,end]
            @async remotecall_fetch(p) do
                rows, _ = localindices(grid)
                lgrid = localpart(grid)
                lm, ln = size(lgrid)
                view(lgrid, :, ln) .= grid[rows,2]
                nothing
            end
        end
    end
    nothing
end

function step!(new::DArray, old::DArray, ::ProcParallel)
    prepare!(old)
    @sync for p in DA.procs(old)
        @async remotecall_fetch(p) do
            @assert localindices(new) == localindices(old)
            m, n = size(old)
            rows, cols = localindices(old)

            # start and length of the window to compute in:
            r = c = 1
            h = length(rows)
            w = length(cols)

            # check whether to cut of mirrored edges of the grid:
            if rows[1] == 1
                h -= 1
                r = 2
            end
            if rows[end] == m
                h -= 1
            end
            if cols[1] == 1
                w -= 1
                c = 2
            end
            if cols[end] == n
                w -= 1
            end

            newindices = (r:r+h-1, c:c+w-1)
            oldindices = (rows[newindices[1]], cols[newindices[2]])

            life_rule!(localpart(new), old, newindices, oldindices)
            nothing
        end
    end
end


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

mutable struct ApplicationGatewayPrivateLinkResourceListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{ApplicationGatewayPrivateLinkResource} } # spec name: value
    nextLink::Any # spec type: Union{ Nothing, String } # spec name: nextLink

    function ApplicationGatewayPrivateLinkResourceListResult(;value=nothing, nextLink=nothing)
        o = new()
        validate_property(ApplicationGatewayPrivateLinkResourceListResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(ApplicationGatewayPrivateLinkResourceListResult, Symbol("nextLink"), nextLink)
        setfield!(o, Symbol("nextLink"), nextLink)
        o
    end
end # type ApplicationGatewayPrivateLinkResourceListResult

const _property_map_ApplicationGatewayPrivateLinkResourceListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"), Symbol("nextLink")=>Symbol("nextLink"))
const _property_types_ApplicationGatewayPrivateLinkResourceListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{ApplicationGatewayPrivateLinkResource}", Symbol("nextLink")=>"String")
Base.propertynames(::Type{ ApplicationGatewayPrivateLinkResourceListResult }) = collect(keys(_property_map_ApplicationGatewayPrivateLinkResourceListResult))
Swagger.property_type(::Type{ ApplicationGatewayPrivateLinkResourceListResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayPrivateLinkResourceListResult[name]))}
Swagger.field_name(::Type{ ApplicationGatewayPrivateLinkResourceListResult }, property_name::Symbol) =  _property_map_ApplicationGatewayPrivateLinkResourceListResult[property_name]

function check_required(o::ApplicationGatewayPrivateLinkResourceListResult)
    true
end

function validate_property(::Type{ ApplicationGatewayPrivateLinkResourceListResult }, name::Symbol, val)
end


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

# helper function to print a large iterator
# in multiple lines with a given tabulation
function io_lines(itr, tab="  ")
  vec = collect(itr)
  N = length(vec)
  I, J = N > 10 ? (5, N-4) : (N, N+1)
  lines = [["$(tab)└─$(vec[i])" for i in 1:I]
           (N > 10 ? ["$(tab)⋮"] : [])
           ["$(tab)└─$(vec[i])" for i in J:N]]
  join(lines, "\n")
end


text:
module PauliStrings

using Combinatorics
using SparseArrays
using LinearAlgebra

export PauliString
export coef, set_coef!, conj, conj!, commutePauli
export adH, simplify_stringList, ∥, show

include("pauliStringStruct.jl")
include("convertPauliString.jl")

# to do: some routines for operator evolution
#include("operatorEvo.jl")

end # module

max_stars_repo_name: Gawatz/PauliStrings

max_issues_repo_name: Gawatz/PauliStrings

max_forks_repo_name: Gawatz/PauliStrings


text:
######################
# 1: Attributes of ToricVarieties
######################

@doc Markdown.doc"""
    dim( v::AbstractNormalToricVariety )

Computes the dimension of the normal toric variety `v`.
"""
function dim( v::AbstractNormalToricVariety )
    return v.polymakeNTV.FAN_DIM::Int
end
export dim

@doc Markdown.doc"""
    dim_of_torusfactor( v::AbstractNormalToricVariety )

Computes the dimension of the torus factor of the normal toric variety `v`.
"""
function dim_of_torusfactor( v::AbstractNormalToricVariety )

    if has_torusfactor( v ) == false
        return 0
    end

    dimension_of_fan = v.polymakeNTV.FAN_DIM::Int
    ambient_dimension = v.polymakeNTV.FAN_AMBIENT_DIM::Int
    return ambient_dimension - dimension_of_fan
end
export dim_of_torusfactor

@doc Markdown.doc"""
    euler_characteristic( v::AbstractNormalToricVariety )

Computes the Euler characteristic of the normal toric variety `v`.
"""
function euler_characteristic( v::AbstractNormalToricVariety )
    f_vector = Vector{Int}(v.polymakeNTV.F_VECTOR)
    return f_vector[ dim( v ) ]
end
export euler_characteristic

"""
    nef_cone( v::NormalToricVariety )

Computes the nef cone of the normal toric variety `v`.

# Examples
```jldoctest
julia> pp = projective_space(2)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2

julia> nef = nef_cone(pp)
A polyhedral cone in ambient dimension 1

julia> dim(nef)
1
```
"""
nef_cone( v::NormalToricVariety ) = Cone( v.polymakeNTV.NEF_CONE )
export nef_cone

"""
    mori_cone( v::NormalToricVariety )

Computes the mori cone of the normal toric variety `v`.

# Examples
```jldoctest
julia> pp = projective_space(2)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2

julia> mori = mori_cone(pp)
A polyhedral cone in ambient dimension 1

julia> dim(mori)
1
```
"""
mori_cone( v::NormalToricVariety ) = Cone( v.polymakeNTV.MORI_CONE )
export mori_cone

######################
# 2: Methods of ToricVarieties
######################

@doc Markdown.doc"""
    ith_betti_number( v::AbstractNormalToricVariety, i::Int )

Compute the i-th Betti number of the normal toric variety `v`.
"""
function ith_betti_number( v::AbstractNormalToricVariety, i::Int )
    if isodd(i)
        return 0
    end
    k = div(i, 2)
    f_vector = Vector{Int}(v.polymakeNTV.F_VECTOR)
    pushfirst!(f_vector, 1)
    betti_number = sum( (-1)^(i-k) * binomial(i,k) * f_vector[ dim( v ) - i + 1 ] for i=k:dim( v ))
    return betti_number

end
export ith_betti_number

@doc Markdown.doc"""
    toric_ideal_binomial_generators(antv::AffineNormalToricVariety)

Get the exponent vectors corresponding to the generators of the toric ideal
associated to the affine normal toric variety `antv`.

# Examples
Take the cyclic quotient singularity corresponding to the pair of integers
`(2,5)`.
```jldoctest
julia> C = Oscar.positive_hull([-2 5; 1 0])
A polyhedral cone in ambient dimension 2

julia> antv = AffineNormalToricVariety(C)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2

julia> toric_ideal_binomial_generators(antv)
pm::Matrix<long>
-1 -1 2 1
-1 0 3 -1
0 -1 -1 2
```
"""
function toric_ideal_binomial_generators(antv::AffineNormalToricVariety)
    pmntv = pm_ntv(antv)
    result = pmntv.TORIC_IDEAL.BINOMIAL_GENERATORS
    return result
end
export toric_ideal_binomial_generators