text:
RealVector = Union{Array{Float64},Array{Real},Array{Int}}
import Distributions.estimate
rhoxb(x::Real, b::Real) = 2*b*b + 2.5 - sqrt(4*b^4 + 6*b*b+2.25 - x*x - x/b)
function multiply!(des::RealVector, x::RealVector, y::Real, n::Int=length(x))
    for i in 1:n
        @inbounds des[i] = x[i]*y
    end
end
multiply!(x::RealVector, y::Real) = multiply!(x, x, y)
function divide!(des::RealVector, x::RealVector, y::Real, n::Int=length(x))
    for i in 1:n
        @inbounds des[i] = x[i]/y
    end
end
divide!(x::RealVector, y::Real) = divide!(x, x, y)
function minus!(des::RealVector, y::Float64, x::RealVector, n::Int64=length(x))
   for i in 1:n
       @inbounds des[i] = y - x[i]
   end
   nothing
end
function add!(x::Vector{Float64}, y::Float64, n::Int64=length(x))
   for i in 1:n
       @inbounds x[i] = x[i] + y
   end
   nothing
end

function factorial2(x::Int)
    k = x;
    prod = 1;
    while (k > 0)
        prod *= k;
        k-=2;
    end
    return(prod);
end

function abs2!(des::RealVector, x::RealVector, n::Int64=length(x))
   for i in 1:n
       @inbounds des[i] = abs2(x[i])
   end
   nothing
end
import Base.abs2
function abs2(x::RealVector)
    return abs(x'*x)
end
function abs2(x::Array{Float64,2})
    return [x[i,:]'*x[i,:] for i=1:size(x,1)]
end
# # Second order gaussiankernel
# function gaussiankernel(x::Real, xdata::RealVector, h::Real, w::Vector, n::Int)
#     h1= 1.0/h
#     tmp = log(h) + log2π/2
#     for ind in 1:n
#         @inbounds w[ind]=-0.5*abs2((x - xdata[ind])*h1) - tmp
#     end
#     w .= exp.(w)
#
#     nothing
# end
#
# # MultiVariate prototype
# function gaussiankernel(x::RealVector, xdata::RealVector, h::RealVector, w::Vector, n::Int)
#     h1= 1.0/h
#     tmp = log(prod(h)) + log2π/2
#     x1 =x;
#     xdata1 = xdata;
#     for d = 1: size(h,1)
#         @inbounds x1[d] = x[d] / h[d]
#         @inbounds xdata1[:,d] = xdata[:,d] / h[d]
#         d+=1
#     end
#
#     for ind in 1:n
#         @inbounds w[ind]=-0.5*abs2(x1 - xdata1[ind,:]) - tmp
#     end
#     # add!(w, tmp, n)
#     w .= exp.(w)
#
#     nothing
# end

# Kern Function
#____________________
kern_epan_dict = Dict(2=> u->ifelse(abs2(u)>=1.0, 0.0, 0.75*(1-abs2(u))),
                      4=> u->ifelse(abs2(u)>=1.0, 0.0, 15/32*(7*(abs2(u)^2)-10*abs2(u)+3)),
                      6=> u->ifelse(abs2(u)>=1.0, 0.0, 105/256*(-33*(abs2(u)^3)+63*(abs2(u)^2)-35*abs2(u)+5)));
kern_biw_dict = Dict();
kern_triw_dict = Dict();
kern_gaussian_dict = Dict(2=> u-> (1 /sqrt(2*π)) * exp(-abs2(u)/2),
                    4=> u-> ifelse((3 - abs2(u))<=0, 0.0, 1 / (2*sqrt(2*π)))*(3 - abs2(u))*exp(-abs2(u)/2),
                    6=> u-> ifelse((15-10*abs2(u)+abs2(u)^2)<=0,0.0 ,1 / (8*sqrt(2*π)))*(15-10*abs2(u)+abs2(u)^2)*exp(-abs2(u)/2) );
kern_dict = Dict(:epan=>kern_epan_dict,:biw=>kern_biw_dict,
                :triw=>kern_triw_dict,:gaussian=>kern_gaussian_dict)
# Bandwidth
#____________________
κ_epan_dict = Dict(2=>1/5,4=>-1/21,6=>5/429);
κ_biw_dict  = Dict(2=>1/7,4=>-1/33,6=>1/143);
κ_triw_dict = Dict(2=>1/9,4=>-3/143,6=>1/221);
κ_gaussian_dict = Dict(2=>1.0,4=>-3.0,6=>15.0 );
κ_dict = Dict(:epan=>κ_epan_dict,:biw=>κ_biw_dict,:triw=>κ_triw_dict,:gaussian=>κ_gaussian_dict)

R_epan_dict = Dict(2=>3/5,4=>5/4,6=>1575/832);
R_biw_dict  = Dict(2=>5/7,4=>805/572,6=>29295/14144);
R_triw_dict = Dict(2=>350/429,4=>3780/2431,6=>301455/134368);
R_gaussian_dict = Dict(2=>1.0/(2 * sqrt(π)),4=>27 / (32 * sqrt(π)),6=>2265/(2048 * sqrt(π) ));
R_dict = Dict(:epan=>R_epan_dict,:biw=>R_biw_dict,:triw=>R_triw_dict,:gaussian=>R_gaussian_dict)

function bw_constant(q::Int,ν::Int,k_type::Symbol)
# Hansen, B. E. (2009). Lecture Notes on Nonparametrics.
# Retrieved from https://www.ssc.wisc.edu/~bhansen/718/NonParametrics1.pdf
    if !(ν in [2,4,6])    error("Kernel order has to be 2,4,6") end
    # if !(k_type in [:epan,:biw,:triw,:gaussian])    error("Kernel order has to be :epan,:biw,:triw,:gaussian")
    κ = κ_dict[k_type][ν];
    R = R_dict[k_type][ν];
    if q > 1
        Cₖ = ( (π^(q/2))*(2^(q+ν-1))*((factorial(ν))^2)*(R^q));
        Cₖ /= (ν*(κ^2)*(factorial2(2*ν-1)+(q-1)*(factorial2(ν-1)^2)));
        Cₖ = Cₖ^(1/(2*ν + q));
    elseif q == 1
         Cₖ = π^(1/2)*(factorial(ν)^3)*R;
         Cₖ /= (2*ν*factorial(2*ν)*κ^2);
         Cₖ = Cₖ^(1/(2*ν+1));
         Cₖ *= 2;
    else
        error();
    end
    return(Cₖ)
end

function bw_compute(σ::Union{Real,RealVector},n::Int, q::Int,ν::Int,k_type::Symbol)
    Cₖ = bw_constant(q,ν,k_type);
    bw = σ / Cₖ;
    bw *= n^(-1/(2*ν+q));
    return(bw);
end
# Kernel Type
#____________________

mutable struct Kernel
    order::Int
    type::Symbol
    kern::Function
    function Kernel(k_type::Symbol,ν::Int)
        # Constructor for the kernel function
        if !(ν in [2,4,6]) error("Kernel order has to be 2,4,6") end
        kern=kern_dict[k_type][ν];
        return(new(ν,k_type,kern));
    end
    # TODO: Can add self-defined kernel function
end

# Fourth order Epanechnikov ekernel
function compute_w(x::Real, xdata::RealVector, h::Real, w::Vector, n::Int, kern::Function)
    ind = 1
    ind_end = 1+n
    @inbounds while ind < ind_end
        u = (x - xdata[ind]) / h
        w[ind] = kern(u);
        ind += 1
    end
    multiply!(w, 1 / h)
    nothing
end
function compute_w(x::RealVector,xdata::RealVector,h::RealVector,w::Vector,n::Int,kern::Function)
    ind = 1
    d = 1
    ind_end = 1+n
    u_all = (x' .- xdata) ./ h';
    @inbounds while ind < ind_end
        u = u_all[ind,:];
        w[ind] =  kern(u);
        ind += 1
    end
    multiply!(w, 1 /prod(h))
    nothing
end

# Fourth order Epanechnikov ekernel
ekernel4 = (x,xdata,h,w,n)->compute_w(x,xdata,h,w,n,kern_dict[:epan][4])
ekernel2 = (x,xdata,h,w,n)->compute_w(x,xdata,h,w,n,kern_dict[:epan][2])
# ApproxFn
#____________________

mutable struct ApproxFn
    xdata::Union{Real,RealVector}
    y::RealVector #Realized value
    n::Int #Number of observation
    h::Union{Real,RealVector} #bandwidth
    q::Int #Number of regressors
    kern::Kernel #Kernel type
    function ApproxFn(xdata::Union{Real,RealVector},y::RealVector,k_type::Symbol,ν::Int)
        n = size(xdata,1);
        q = size(xdata,2);
        if length(y) != n
            error("The dimension of x,y must match");
        end
        if (q>1)
            σ = sqrt.(diag(cov(xdata)));
        elseif (q==1)
            σ = sqrt(cov(xdata));
        end
        h=bw_compute(σ,n,q,ν,k_type);
        if (q==1)
            h=h[1];
        end
        # self = new(x,y,n,"SVR",SVR(kernel="rbf",degree=4,gamma="scale"));
        # self = new(x,y,n,"KNN",KNeighborsRegressor());
        # fit!(self.model,x,y);
        kern=Kernel(k_type,ν);
        self=new(xdata,y,n,h,q,kern);
        return(self)
    end

    function (self::ApproxFn)(x::Union{Real,RealVector})
        if ((typeof(x)<:Real) & (self.q==1))
            return(forecast(x,self,0));
        elseif (typeof(x)<:Vector)& (self.q>1) &(size(x)[1]==self.q)
            return(forecast(x,self,0));
        elseif (typeof(x)<:Vector)&(self.q==1)
            forecast_y = zeros(size(x)[1]);
            for i = 1:size(x)[1]
                forecast_y[i] = forecast(x[i],self,0)[1];
            end
            return(forecast_y);
        elseif (self.q>1)&(size(x)[2]==self.q)
            forecast_y = zeros(size(x)[1]);
            for i = 1:size(x)[1]
                forecast_y[i] = forecast(x[i,:],self,0);
            end
            return(forecast_y);
        else
            error("Dimension of input incorrect")
        end
    end
end

function find_nn(neighbour_dist::Array{Float64,1},n::Int)
    neighbour_dist_sort = unique(sort(neighbour_dist))[1:n]; #Sort from smallest to
    # largest
    index=zeros(n);
    for i = 1:n
         for j = 1:length(neighbour_dist)
             if (neighbour_dist_sort[i]==neighbour_dist[j])
                index[i]=j;
            end
         end
    end
    return(convert(Array{Int,1},index));
end

function apply_row(fun,xdata)
    n_rows=size(xdata,1);
    n_data_col=size(xdata,2);
    # n_cols=size(xdata,2);
    if n_data_col > 1
        n_cols = length(fun(xdata[1,:]));
        if n_cols > 1
            val = zeros(n_rows,n_cols);
            for i = 1:n_rows
                val[i,:]=fun(xdata[i,:])
            end
        else
            val = zeros(n_rows);
            for i = 1:n_rows
                val[i]=fun(xdata[i,:])
            end
        end
    else
        n_cols = length(fun(xdata[1]));
        if n_cols > 1
            val = zeros(n_rows,n_cols);
            for i = 1:n_rows
                val[i,:]=fun(xdata[i])
            end
        else
            val = zeros(n_rows);
            for i = 1:n_rows
                val[i]=fun(xdata[i])
            end
        end
    end
    return(val);
end

function estimatenn(x::Union{Real,RealVector},af::ApproxFn,n::Int)
    # Use 10 nearest neighbour, I don't know, may be should do more
    dist2=y->abs2(y.- x);
    # n=10;
    neighbour_dist = apply_row(dist2,af.xdata);
    nn_index=find_nn(neighbour_dist,n);
    w=zeros(af.n);
    for i in nn_index
        w[i]=1/n;
    end
    w_diag = diagm(0=>w);
    # β_kernel = inv(af.xdata'*w_diag * af.xdata) * (af.xdata' * w_diag * af.y);
    return(w_diag);
end
# The ApproxFn
# The nonparametric function such that Y = ApproxFn(X)
function estimate(x::Union{Real,RealVector},af::ApproxFn)
    w = zeros(af.n);
    compute_w(x,af.xdata,af.h,w,af.n,af.kern.kern);
    # if (sum(w)==0)
    # @show sum(w.>0);
    if (sum(w.>0) < 5)
        nn = max(1,sum(w.>0));
        return(estimatenn(x,af,nn));
    else
        w_diag = diagm(0=>w);
        # β_kernel = inv(af.xdata'*w_diag * af.xdata) * (af.xdata' * w_diag * af.y);
        return(w_diag);
    end
end

function forecast(x::Union{Real,RealVector},af::ApproxFn,p::Int)
    w_diag = estimate(x,af);
    xdata = ones(af.n);
    diff_pop=apply_row(y->y-x,af.xdata);
    for pp = 1:p
        xdata=hcat(x,diff_pop.^p);
    end
    β_kernel = inv(xdata' * w_diag * xdata) * (xdata' * w_diag * af.y);
    return(β_kernel[1]);
end

function forecastnn(x::Union{Real,RealVector},af::ApproxFn)
    # KNN method to forecast
    # Use 10 nearest neighbour, I don't know, may be should do more
    w_diag = estimatenn(x,af);
    x = ones(af.n);
    diff_pop=apply_row(y->y-x,af.xdata);
    for pp = 1:p
        x=hcat(x,diff_pop.^p);
    end
    β_kernel = inv(x' * w_diag * x) * (x' * w_diag * af.y);
    return(β_kernel[1]);
end

function UpdateData(self::ApproxFn,xdata,y)
    self.xdata = xdata;
    self.y     = y;
    # fit!(self.model,self.x,self.y);
    return(self)
end

function UpdateVal(self::ApproxFn,y)
    self.y     = y;
    # fit!(self.model,self.x,self.y);
    return(self)
end

function forecast_fit(Kern::Kernel)
    yfit = zeros(Kern.n)
    for i = 1:Kern.n
        x_i = Kern.xdata[i,:];
        yfit[i] = Kern.forecast(x_i);
    end
    return(yfit);
end


text:
test_folder = "conjunction"
cd(test_folder)

openBF.runSimulation("$test_folder.yml", verbose=true, out_files=false)

p = readdlm("$test_folder"*"_results/d1_P.last")

@test isapprox(minimum(p[:,end]), -1771.0, atol=1)
@test isapprox(mean(p[:,end]), 623.0, atol=1)
@test isapprox(maximum(p[:,end]), 6529.0, atol=1)

rm("$test_folder"*"_results/", recursive=true)
cd("..")

ext: jl


text:
#******************************************************************************#
#                               Introduction                                   #
#******************************************************************************#

# PyPlot is package for Julia that able you to use matplotlib library.

# First step to start plotting in Julia using Pyplot is install PyPlot package.
# You should do it only once.

ENV["PYTHON"] = ""
Pkg.add("PyPlot")

# And start using PyPlot module.
using PyPlot

# Now you can working.

# To do plot use plot():

x = linspace(-2*π, 2*π, 2000) # Generating 2000 sclars from -2π to 2π
y = sin(x)
plot(x,y)

# If you need two lines in one plot:

x = linspace(-2*π, 2*π, 2000) # Generating 2000 sclars from -2π to 2π
y = sin(x)
z = cos(x)
plot(x,y)
plot(x,z)

# But it's only raw chart. Let's describe it:

plot(x,y, color = "green", linewidth = "2.0")
plot(x,z, color = "red", linewidth = "2.0", "--")
title("Sine wave")
ylabel("Amplitude")
xlabel("Time")
grid(true)
text(0, 0, "0")
legend(["sin(x)","cos(x)"], loc = "lower left", fontsize = 12, borderpad = 0.3,
    handlelength=3)

# You can do subplots using subplot(nrows,ncols, plot_number):

sinh = sinh(x)
cosh = cosh(x)

subplot(221)
    plot(x,y)
    legend(["sin(x)"], fontsize = 8,  loc = "best" )
subplot(222)
    plot(x,z)
    legend(["cos(x)"], fontsize = 8,  loc = "best"  )
subplot(223)
    plot(x,sinh)
    legend(["sinh(x)"], fontsize = 8,  loc = "best"  )
subplot(224)
    plot(x,cosh)
    legend(["cosh(x)"], fontsize = 8,  loc = "best"  )

#******************************************************************************#
#                               Scatter plots                                  #
#******************************************************************************#

n = 1024
X = rand(1,1,n)
Y = rand(1,1,n)
colors = rand(1,1,n)

scatter(X,Y, s = 50 ,c = colors)
scatter(X,Y, s = 50 ,c = colors, marker = "^", alpha = .5)

#******************************************************************************#
#                                  Bar plots                                   #
#******************************************************************************#

n = 15
X = range(0,n)
Y1 = (1-X/n) * rand(1:n)
Y2 = (1-X/n) * rand(1:n)
axes([0.025,0.025,0.95,0.95])
bar(X, +Y1, facecolor="#9999ff", edgecolor="white")
bar(X, -Y2, facecolor="#ff9999", edgecolor="white")
margins(0.05) # 5% margins

# Code below generating describe under or above each bar.

for (x,y) in zip(X,Y1)
    text(x+0.4, y+0.05, "$(round(y,2))", ha="center", va= "bottom")
end

for (x,y) in zip(X,Y2)
    text(x+0.4, -y-0.05, "$(round(y,2))",  ha="center", va= "top")
end

#******************************************************************************#
#                               Contour plots                                  #
#******************************************************************************#

f(x,y) = (1-x/2+x^5+y^3)*exp(-x^2-y^2)
n = 256
x = linspace(-3,3,n)
y = linspace(-3,3,n)
z = zeros(n,n)
for i in 1:n
    for j in 1:n
        z[i,j] = f(x[j],y[i])
    end
end
X = repmat(x',n,1)
Y = repmat(y,1,n)

contourf(X, Y, z, 8, alpha=.75, cmap="jet")

# Easier way to get contour plots is to use Gadfly package:
Pkg.add("Gadfly")

using Gadfly

plot((x,y) -> (1-x/2+x^5+y^3)*exp(-x^2-y^2), -4, 4, -4, 4)

# You can use more parameters. To get more information about Gadfly contour
# plots check this topic: https://github.com/dcjones/Gadfly.jl/issues/293

plot((x,y) -> (1-x/2+x^5+y^3)*exp(-x^2-y^2), -4, 4, -4, 4,
    Stat.contour(samples = 200, levels = linspace(-1,1,25)))

# In this example we used levels = linspace(-1,1,25) to show 25 lines from -1
# to 1 include 0. If you want to see more lines you can precise it in
# levels parameter.

#******************************************************************************#
#                                Pie charts                                    #
#******************************************************************************#

# Back to using PyPlot

using PyPlot

labels = ["PyPlot";"Plots";"Gadfly"]
colors = ["yellowgreen";"lightskyblue";"lightcoral"]
sizes = NaN*zeros(3)
explode = zeros(length(sizes))
explode[2] = 0.1 # Move slice 2 out by 0.1
sizes = [80,15,10]
font = Dict("fontname"=>"Sans","weight"=>"semibold")
axes([0.025,0.025,0.95,0.95])
pie(sizes,
    labels=labels,
    shadow=true,
    startangle=90,
    explode=explode,
    colors=colors,
    autopct="%1.1f%%",
    textprops=font
)

#******************************************************************************#
#                                   3D Plots                                   #
#******************************************************************************#

# Example no. 1

f(x) = (x[1]^2 + x[2]^2)

n = 100
x = linspace(-1, 1, n)
y = linspace(-1,1,n)

xgrid = repmat(x',n,1)
ygrid = repmat(y,1,n)

z = zeros(n,n)

for i in 1:n
    for j in 1:n
        z[i:i,j:j] = f([x[i],y[j]])
    end
end

plot_wireframe(xgrid,ygrid,z)

# Example no. 2

f(x,y) = sin(x)*cos(y)

n = 100
x = linspace(-2, 2, n)
y = linspace(-1, 1,n)

xgrid = repmat(x',n,1)
ygrid = repmat(y,1,n)

z = zeros(n,n)
for i in 1:n
    for j in 1:n
        z[i,j] = f(x[j],y[i])
    end
end

plot_surface(xgrid,ygrid,z, rstride=1, cstride=1, cmap="coolwarm",
                       linewidth=0)

# More 3D plots examples you can find there:
# http://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html

#******************************************************************************#
#                               Plots package                                  #
#******************************************************************************#

# You can use simpliest package - Plots.
Pkg.add("Plots")

using Plots

x = linspace(-2*π, 2*π, 2000) # Generating 2000 scalars from -2π to 2π
y = sin(x)
z = cos(x)
plot(x,y)
plot!(x,z)

# You can find more about plots there:
# https://en.wikibooks.org/wiki/Introducing_Julia/Plotting.

# To master skill in PyPlot you should visit this link:
# http://matplotlib.org/contents.html
# If you prefer Gadfly, visit manual page:
# http://gadflyjl.org/stable/

max_stars_repo_name: plon-io/Plotting-in-Julia

max_issues_repo_name: plon-io/Plotting-in-Julia

max_forks_repo_name: plon-io/Plotting-in-Julia


text:
using Nemo

# Globals for the spaces and symbols associated with the field of fractions
# of integer polynomials used as entries in the matrix space for the Burau
# representation.
# We represent the polynomials in `Z[t, t′]` as the field of fractions of
# `Z[t]` for computational simplicity. It can be trivially shown the two are
# isomorphic so, most importantly perhaps, the property remains that the
# trivial braid is represented by the identity matrix.
BurauR, Buraux = PolynomialRing(ZZ, "t")
BurauFF = FractionField(BurauR)

"""
Unreduced Burau representation for a single generator. (Helper)
"""
function unreducedburaugen(i, isinv, n, MS)
    m = MS()
    for j = 1:i-1
        m[j,j] = one(BurauFF)
    end
    if isinv
        m[i,i+1] = one(BurauR)
        m[i+1,i] = one(BurauFF) * inv(BurauFF(Buraux))
        m[i+1,i+1] = one(BurauFF) - (one(BurauFF) * inv(BurauFF(Buraux)))
    else
        m[i,i] = one(BurauFF) - BurauFF(Buraux)
        m[i+1,i] = one(BurauFF)
        m[i,i+1] = BurauFF(Buraux)
    end
    for j = i+2:n
        m[j,j] = one(BurauFF)
    end
    m
end

"""
Get the unreduced Burau representation of `b0`. This will come in the form of
a matrix with entries in the field of fractions of `Z[t]`, the ring of
polynomials with integer coefficients. If it is equal to the standard identity
matrix, then `b0` is ambient isomorphic to the trivial braid.
"""
function unreducedburau(b0::BraidElement)
    b = unprettybraid(b0)
    n = b.group.n
    MS = MatrixSpace(BurauFF, n, n)
    mats = []
    for i = 1:length(b.crossings)
        x = b.crossings[i]
        push!(mats, unreducedburaugen(x[1], x[2] < 0, n, MS))
    end
    m = MS()
    for i = 1:n
        m[i,i] = one(BurauFF)
    end
    for x in mats
        m = m * x
    end
    m
end


text:
function comm!(commAB, A, B)
	N = div(size(A)[1], 4)
	for j = 1: N
		indice = 4*(j - 1)
		comm2ls!(commAB, A, B, indice)
	end
	return nothing
end

function comm2ls!(commAB, A, B, indice)
	i1, i2, i3, i4 = indice + 1, indice + 2, indice + 3, indice + 4
	commAB[i1] = (A[i2]*B[i3] - A[i3]*B[i2])*2
	commAB[i2] = (A[i3]*B[i1] - A[i1]*B[i3])*2
	commAB[i3] = (A[i1]*B[i2] - A[i2]*B[i1])*2
	commAB[i4] = 0
	return nothing
end
function comm2ls!(commAB, A, B)
	commAB[1] = (A[2]*B[3] - A[3]*B[2])*2
	commAB[2] = (A[3]*B[1] - A[1]*B[3])*2
	commAB[3] = (A[1]*B[2] - A[2]*B[1])*2
	commAB[4] = 0
	return nothing
end
#Decompose a generic 2x2 Hermitian matrix M in the Pauli basis.
function pauli_decomp(M)
	c = zeros(4)
	c[1] = real(M[1, 2])
	c[2] = imag(M[2, 1])
	c[3] = real(M[1, 1] - M[2, 2])/2
	c[4] = real(M[1, 1] + M[2, 2])/2
	return c
end

function tensor2matrix(tensor)
	N = size(tensor)[3]
	matrix = zeros(Complex{Float64}, 2*N, 2*N)
	for n=1:N
		matrix[2*(n - 1) + 1: 2*n, 2*(n - 1) + 1: 2*n] .= tensor[:, :, n]
	end
	return matrix
end

function vector2matrix(v::Array{<:Number,1}, σ)
	#matrix = zeros(Complex{Float64}, 2, 2)
	matrix = v[1].*σ[1] .+ v[2].*σ[2] .+ v[3].*σ[3] .+ v[4].*σ[4]
	return matrix
end
function vector2matrix2(v, σ)
	#matrix = zeros(Complex{Float64}, 2, 2)
	matrix = v[1]*σ[1] + v[2]*σ[2] + v[3]*σ[3] + v[4]*σ[4]
	return matrix
end

function vector2matrix!(matrix, v, σ)
	matrix .= v[1].*σ[1] .+ v[2].*σ[2] .+ v[3].*σ[3] .+ v[4].*σ[4]
	return nothing
end

function vector2matrix(v::Array{<:Number,2}, σ)
	N = size(v)[1]
	matrix = zeros(Complex{Float64}, 2*N, 2*N)
	vj = zeros(Complex{Float64}, 4)
	for j=1:N
		vj .= v[:,j]
		matrix[2*(j - 1) + 1: 2*j, 2*(j - 1) + 1: 2*j] .= vector2matrix(vj , σ)
	end
	return matrix
end

function vector2matrix!(matrix, v::Array{<:Number,2}, σ)
	N = size(v)[1]
	vj = zeros(Complex{Float64}, 4)
	for j=1:N
		vj .= v[:,j]
		matrix[2*(j - 1) + 1: 2*j, 2*(j - 1) + 1: 2*j] .= vector2matrix(vj , σ)
	end
	return nothing
end

function fidelity(U, U_ideal, dims)
	U_ideal_dag = adjoint(U_ideal)
	M = U_ideal_dag[1: dims, 1: dims]*U[1: dims, 1: dims]
	M_dag = adjoint(M)
	fidelity = (tr(M*M_dag) + abs(tr(M))^2)/(dims*(dims + 1))
	return real(fidelity)
end

function fidelity2(U, U_ideal, dims)
	U_ideal_dag = adjoint(U_ideal)
	M = U_ideal_dag[1: dims, 1: dims]*U[1: dims, 1: dims]
	fidelity = abs(tr(M))/dims
	return fidelity
end

function exp_matrix(v)
	exp = zeros(Complex{Float64}, 4)
	a = norm(v)
	if a > 100.0*eps()
		n = v/a
		sina, cosa = sin(a), cos(a)
		exp[1] = 1im*sina*n[1]
		exp[2] = 1im*sina*n[2]
		exp[3] = 1im*sina*n[3]
		exp[4] = cosa
	else
		exp[4] = 1
	end
	return exp
end

function exp_matrix(v, N)
	exp = zeros(Complex{Float64}, 4, N)
	vj = zeros(Complex{Float64}, 4)
	index = levels .+ 1
	for j=1:N
		vj .= v[:, j]
		exp[:, j] .= exp_matrix(vj)
	end
	return exp
end
#Calculate y.σ = exp[i v.σ] x.σ exp[-i v.σ].
#a.σ = a[1]*σ1 + a[2]*σ2 + a[3]*σ3 + a[4]*σ4
function unitary_tranformation(x, v)
	a = norm(v)
	if a > 100*eps()
		n = SA[v[1]/a, v[2]/a, v[3]/a]
		cosa, sina = cos(a), sin(a)
		y = x*cosa*cosa
		y13 = view(y, 1: 3)
		x13 = view(x, 1: 3)
		p = 2*dot(n, x13)
		y13 .= y13 .- sin(2*a).*cross(n, x13) .+ (sina*sina).*(p.*n .- x13)
	else
		y = x
	end
	return y
end
function unitary_tranformation!(x, v)
	a = norm(v)
	if a > 100*eps()
		n = (v[1]/a, v[2]/a, v[3]/a)
		cosa, sina = cos(a), sin(a)
		cosa2, sina2 = cosa*cosa, sina*sina
		sin2a = 2*sina*cosa#sin(2.0*a)

		nXx = (n[2]*x[3] - n[3]*x[2], n[3]*x[1] - n[1]*x[3], n[1]*x[2] - n[2]*x[1])
		p = 2*(n[1]*x[1] + n[2]*x[2] + n[3]*x[3])
		x[1] = x[1]*cosa2 - sin2a*nXx[1] + sina2*(p*n[1] - x[1])
		x[2] = x[2]*cosa2 - sin2a*nXx[2] + sina2*(p*n[2] - x[2])
		x[3] = x[3]*cosa2 - sin2a*nXx[3] + sina2*(p*n[3] - x[3])
		x[4] = x[4]*cosa2
	end
	return nothing
end
#Calculate the unitary tranformation in the Hilbert space formed by the outer
#product between the two level system space and the cavity space:
#<n| y_n |n> = exp[i v[n,:].σ] x[n,:].σ exp[-i v[n,:].σ]
function unitary_tranformation(x, v, N::Int)
	u = zeros(4, N)
	for j = 1: N
		xj = view(x, :, j)
		vj = view(v, :, j)
		u[:,j] .= unitary_tranformation(xj, vj)
	end
	return u
end
function unitary_tranformation!(x, v, N::Int)
	for j = 1: N
		xj = view(x, :, j)
		vj = view(v, :, j)
		unitary_tranformation!(xj, vj)
	end
	return nothing
end

#Returns the angular frequency vector of a DFT
function fftfreq(n::Integer, dt)
	if n%2 == 0
		v1 = collect(0: 1: div(n, 2) - 1)
		v2 = -collect(div(n, 2): -1: 1)
		v = 2*pi*[v1; v2]/(dt*n)
	else
		v1 = collect(0: 1: div(n - 1,2))
		v2 = -collect(div(n - 1, 2): -1: 1)
		v = 2*pi*[v1; v2]/(dt*n)
	end
	return v
end

max_stars_repo_name: thalesfr/MagnusCorrection

max_issues_repo_name: thalesfr/MagnusCorrection

max_forks_repo_name: thalesfr/MagnusCorrection


text:
import GL
import GLFW

GLFW.Init()
GLFW.OpenWindowHint(GLFW.OPENGL_VERSION_MAJOR, 4)
GLFW.OpenWindowHint(GLFW.OPENGL_VERSION_MINOR, 2)
GLFW.OpenWindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
GLFW.OpenWindowHint(GLFW.OPENGL_FORWARD_COMPAT, 1)
GLFW.OpenWindow(0, 0, 0, 0, 0, 0, 0, 0, GLFW.WINDOW)
GLFW.SetWindowTitle("Hello Cube")

println("Vendor:   ", GL.GetString(GL.VENDOR))
println("Renderer: ", GL.GetString(GL.RENDERER))
println("Version:  ", GL.GetString(GL.VERSION))
println("GLSL:     ", GL.GetString(GL.SHADING_LANGUAGE_VERSION))

while GLFW.GetWindowParam(GLFW.OPENED)
	GLFW.SwapBuffers()
end

GLFW.Terminate()


text:
# Data structure for the simulation start

using Parameters

@with_kw struct MDinput

  temp       :: Float64 = 298.15
  total_time :: Float64 = 10
  dt         :: Float64 = 0.01
  nsteps     :: Int64   = trunc(Int64,total_time / dt)
  iprint     :: Int64   = 10

end

export MDinput


text:
module CoDa

export alpha, alphaTransform, Helmert
export ilr, clr

"""
    Helmert(n)

Calculate ``n×n`` Helmert matrix
"""
function Helmert(n)
    H = zeros(n,n)
  for i = 1:n, j = 1:n
      if i==1
        H[i,j] = 1.0/sqrt(n)
      elseif j < i
        H[i,j] = 1.0/sqrt(i*(i-1))
      elseif i==j
        H[i,j] = (1.0-i)/sqrt(i*(i-1))
      end
  end
    return -H
end

"""
    alphaTransform(X::Matrix, α::Float64)

Apply alpha-transform with Helmert product to ``X``
"""
function alphaTransform(X::Matrix{Float64}, α::Float64)
    D = size(X,2)
    H = Helmert(D)[2:end,:]
    Xᵃ = X.^α
    Xs = sum(Xᵃ, dims=2)
    U = (D*Xᵃ./Xs .- 1)./α
    return U*transpose(H)
end

"""
    alpha(X::Matrix, α::Float64)

Apply alpha-transform without Helmert product
"""
function alpha(X::Matrix{Float64}, α::Float64)
    D = size(X,2)
#    H = Helmert(D)[2:end,:]
    Xᵃ = X.^α
    Xs = sum(Xᵃ,dims=2)
    U = (D*Xᵃ./Xs .- 1)./α
    return U
end

"""
    ilr(X::Matrix)

Isometric logratio transform
"""
function ilr(X::Matrix{Float64})
    D = size(X,2)
    H = Helmert(D)[2:end,:]
    return log.(X) * transpose(H)
end

"""
    clr(X::Matrix)

Centered logratio transform
"""
function clr(X::Matrix{Float64})
    D = size(X,2)
    ratio = X ./ (prod(X,2).^(1/D))
    return log.(ratio)
end

end


text:
################################################################################
#################### Scaled Grid Interpolation #################################
################################################################################
# The following code generates three objects:
    # 1) A "TypeScale" type that transforms ranges into a type of scaled ranges
    # 2) A "TypeRange" type that foroms scaled ranges of a given type.
    #   These ranges behave "just as" normal ranges
    # 3) A ScaledInterpolations function that wraps "interpolations.jl" to be used
    #   with TypeRange ranges
# Types of Scales: parameters bounds for grid (a,b) and curvature (θ)
    # p = PolyScale(a,b,θ)
        # Maps a number x∈[0,1] to a polynomial grid with curvature θ: p(x)∈[a,b]
        # Usage: p(x), where x is a number, or p.(x) where x is vector/range
        # p(x) = a + (b - a) * x^θ
    # p = InvPolyScale(a,b,θ)
        # Maps a number y∈[a,b] from a polynomial grid with curvature θ to the unit interval: p(x)∈[0,1]
        # Usage: p(y), where x is a number, or p.(y) where y is vector/range
        # p(y) = ((y-a) / (b-a) )^(1/θ)
    # p = ExpScale(a,b,θ)
        # Maps a number x∈[0,1] to a exponential grid with curvature θ: p(x)∈[a,b]
        # Usage: p(x), where x is a number, or p.(x) where x is vector/range
        # p(x) = a + (b - a) * ( (exp(θ*x) - 1)/(exp(θ) - 1) )
    # p = LogScale(a,b,θ)
        # Maps a number y∈[a,b] from an exponential grid with curvature θ to the unit interval: p(x)∈[0,1]
        # Usage: p(y), where x is a number, or p.(y) where y is vector/range
        # p(y) = log( (y-a)*(exp(θ) - 1)/(b - a) + 1 ) / θ
# Types of Scaled Ranges
    # grid = PolyRange(a,b,θ=t,N=n) = PolyScale(a,b,t)(range(0,1,length=n))
    # grid =  ExpRange(a,b,θ=t,N=n) =  ExpScale(a,b,t)(range(0,1,length=n))
# ScaledInterpolations
    # itp = ScaledInterpolations(grid,f_grid,Itp_Type)
        # grid is either a PolyRange or ExpRange type
        # f_grid is a vector with the values of the function at the grid nodes
        # Itp_Type is a type from the interpolations.jl package
            # Examples: BSpline(Cubic(Line(OnGrid()))) // shortcut: CubicSpline
            #           BSpline(Cubic(Flat(OnGrid())))
            #           FritschButlandMonotonicInterpolation() // shortcut: MonotoneSpline
            #           BSpline(Linear()) // shortcut: LinearInterp
# ScaledInterpolations with multiple dimensions
    # itp_md = ScaledInterpolations( (grid_1,...,grid_N) , (f_grid_1,...,f_grid_N) , (Type_1,...,Type_N)  )
# Jacob Adenbaum, 2020

#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# Define Scaled types
    #-------------------------------------------------------------------------------
    # Polynomial Scaling
    abstract type Scale{T}       end
    abstract type ScaledRange{T} end

    for Foo in [:PolyScale, :InvPolyScale]

        @eval struct $Foo{T,TF} <: Scale{T}
            a::T
            b::T
            θ::TF
        end

        @eval function $Foo(a, b, θ = 1)
            T = reduce(promote_type, map(typeof, (a,b)))
            $Foo(convert(T, a), convert(T, b), θ)
        end
    end

    (p::PolyScale)(x)           = p.a + (p.b - p.a) * x^p.θ
    (p::InvPolyScale)(y)        = ((y-p.a) / (p.b-p.a) )^(1/p.θ)
    Base.inv(p::PolyScale)      = InvPolyScale(p.a, p.b, p.θ)
    Base.inv(p::InvPolyScale)   = PolyScale(p.a, p.b, p.θ)

    #-------------------------------------------------------------------------------
    # Exponential Scaling
    for Foo in [:ExpScale, :LogScale]

        @eval struct $Foo{T, TF} <: Scale{T}
            a::T
            b::T
            θ::TF
            et::T
            s::T
        end

        @eval function $Foo(a, b, θ = 1)
            s = (exp(θ) - 1)/(b - a)
            et= exp(θ)
            T = reduce(promote_type, map(typeof, (a,b,θ, et,s)))
            $Foo(convert(T, a), convert(T, b), θ, et,s)
        end
    end

    function (p::ExpScale)(x)
        t  = (exp(p.θ * x) - 1)/(p.et - 1)
        return p.a + (p.b - p.a) * t
    end

    function (p::LogScale)(y)
        θ  = p.θ
        log( (y-p.a)*p.s + 1 ) / p.θ
    end

    Base.inv(p::ExpScale)       = LogScale(p.a, p.b, p.θ)
    Base.inv(p::LogScale)       = ExpScale(p.a, p.b, p.θ)

    #-------------------------------------------------------------------------------
    # Finalize Types
    scales = Dict(:PolyRange => :PolyScale, :ExpRange => :ExpScale)

    Base.inv(::Type{ExpScale}) = LogScale
    Base.inv(::Type{LogScale}) = ExpScale
    Base.inv(::Type{PolyScale})= InvPolyScale
    Base.inv(::Type{InvPolyScale}) = PolyScale

#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# Define Scaled Ranges types
    #-------------------------------------------------------------------------------
    # Polynomial Range
    for Typ in [:PolyRange] @eval begin
        struct $Typ{T,TF} <: AbstractRange{T}
            a::T
            b::T
            θ::TF
            N::Int
            vals::Vector{Float64}
            scale::$(scales[Typ]){T,TF}
            invscale::inv($(scales[Typ])){T,TF}
        end

        function $Typ(a,b; θ = 1, N = 100)
            T       = reduce(promote_type, map(typeof, (a,b,θ, 1.0)))
            Ta      = convert(T, a)
            Tb      = convert(T, b)

            vals    = domscale($Typ)(Ta,Tb,θ).(LinRange(0, 1, N))
            scale   = $(scales[Typ])(Ta, Tb, θ)
            invscale= inv($(scales[Typ])(Ta, Tb, θ))
            $Typ(convert(T, a), convert(T, b), θ, N, vals, scale, invscale)
        end

        domscale(p::$Typ)      = p.scale
        domscale(::Type{$Typ}) = $(scales[Typ])

        function rangescale(r::$Typ, x)
            r.invscale(x)*(length(r)-1) + 1
        end

        Base.getindex(p::$Typ, i::Int) = p.vals[i]
        Base.length(p::$Typ)  = p.N
        Base.show(io::IO, p::$Typ) = begin
            Typ = $Typ
            print(io, "$Typ($(p.a), $(p.b), $(p.θ), $(p.N))")
        end

        function Base.searchsortedfirst(p::$Typ, x, args...)
            searchsortedfirst(p.vals, x, args...)
        end

    end end

    #-------------------------------------------------------------------------------
    # Exponential Range
    for Typ in [:ExpRange] @eval begin
        struct $Typ{T,TF} <: AbstractRange{T}
            a::T
            b::T
            θ::TF
            N::Int
            vals::Vector{Float64}
            scale::$(scales[Typ]){T,TF}
            invscale::inv($(scales[Typ])){T,TF}
        end

        function $Typ(a,b; θ = 1, N = 100)
            vals    = domscale($Typ)(a,b,θ).(LinRange(0, 1, N))
            scale   = $(scales[Typ])(a, b, θ)
            invscale= inv($(scales[Typ])(a, b, θ))
            T       = reduce(promote_type, map(typeof, (a,b,θ,scale.s)))
            $Typ(convert(T, a), convert(T, b), θ, N, vals, scale, invscale)
        end

        domscale(p::$Typ)      = p.scale
        domscale(::Type{$Typ}) = $(scales[Typ])

        function rangescale(r::$Typ, x)
            r.invscale(x)*(length(r)-1) + 1
        end

        Base.getindex(p::$Typ, i::Int) = p.vals[i]
        Base.length(p::$Typ)  = p.N
        Base.show(io::IO, p::$Typ) = begin
            Typ = $Typ
            print(io, "$Typ($(p.a), $(p.b), $(p.θ), $(p.N))")
        end

        function Base.searchsortedfirst(p::$Typ, x, args...)
            searchsortedfirst(p.vals, x, args...)
        end
    end end

    #-------------------------------------------------------------------------------
    # Finalize Ranges
    function rangescale(r::AbstractRange, x)
        ((x - first(r)))/(last(r) - first(r)) * (length(r) -1) + 1
    end

    rangescale(r::UnitRange, x::Int) = x

#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# Scaled Interpolation Wraper
    mutable struct ScaledInterpolations{T, N, IT <: Interpolations.AbstractInterpolation,
                                RT <: NTuple{N, AbstractRange},
                                VT <: AbstractArray{T, N}, AT}
        r::RT       # Interpolation grid (ranges)
        v::VT       # Interpolation values
        itp::IT     # Interpolation object
        args::AT    # Arguments to pass to construct interpolant
    end

    function ScaledInterpolations(r, v, args...)
        itp = extrapolate(interpolate(v, args...), Interpolations.Flat())
        return ScaledInterpolations(tuplefy(r), v, itp, args)
    end

    tuplefy(x::Tuple) = x
    tuplefy(x)        = tuple(x)

    """
    ```
    fit!(s::ScaledInterpolations, [v])
    ```
    Re-fit the interpolation after the underlying data has updated.  Works well if
    you use a view into the original array.  If an array of values is passed, copies
    the values to the scaled interpolations value array and then refits the
    interpolation.
    """
    function fit!(s::ScaledInterpolations)
        s.itp = extrapolate(interpolate(s.v, s.args...), Interpolations.Flat())
    end

    function fit!(s::Array{T,N}) where {T <: ScaledInterpolations, N}
        for v in s
            fit!(v)
        end
    end

    dim(::Type{ScaledInterpolations{T, N, IT, RT, VT, AT}}) where {T, N, IT, RT, VT, AT} = N
    Base.eltype(::Type{ScaledInterpolations{T, N, IT, RT, VT, AT}}) where {T, N, IT, RT, VT, AT} = T
    dim(itp::ScaledInterpolations) = dim(typeof(itp))
    Base.eltype(itp::ScaledInterpolations) = eltype(typeof(itp))

    @generated function (s::ScaledInterpolations)(x::Vararg)
        Ns = dim(s)
        Nx = length(x)
        Ns == Nx || begin
            return quote
                Ns = $Ns
                Nx = $Nx
                throw(ArgumentError("Must have $Ns arguments.  You passed $Nx"))
            end
        end

        ex = Expr(:call, :(s.itp))
        for i = 1:Ns
            push!(ex.args, :(rangescale(s.r[$i], x[$i])))
        end
        ex
    end

#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# Scaled Interpolation Shortcuts for Itp_Type
    const CubicSpline    = BSpline(Cubic(Natural(OnGrid())))
    const MonotoneSpline = FritschButlandMonotonicInterpolation()
    const LinearInterp   = BSpline(Linear())
    const LinearInterp2  = BSpline(Linear())

    function ScaledInterpolations(grid, values, type::FritschButlandMonotonicInterpolation)
        return interpolate(collect(grid), values, type)
    end

max_stars_repo_name: py-r-hans/AdvMacro

max_issues_repo_name: py-r-hans/AdvMacro

max_forks_repo_name: py-r-hans/AdvMacro


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

"""
    apply_archive_rule(analyzer_arn, rule_name)
    apply_archive_rule(analyzer_arn, rule_name, params::Dict{String,<:Any})

Retroactively applies the archive rule to existing findings that meet the archive rule
criteria.

# Arguments
- `analyzer_arn`: The Amazon resource name (ARN) of the analyzer.
- `rule_name`: The name of the rule to apply.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: A client token.
"""
apply_archive_rule(analyzerArn, ruleName; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("PUT", "/archive-rule", Dict{String, Any}("analyzerArn"=>analyzerArn, "ruleName"=>ruleName, "clientToken"=>string(uuid4())); aws_config=aws_config)
apply_archive_rule(analyzerArn, ruleName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("PUT", "/archive-rule", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn, "ruleName"=>ruleName, "clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    create_access_preview(analyzer_arn, configurations)
    create_access_preview(analyzer_arn, configurations, params::Dict{String,<:Any})

Creates an access preview that allows you to preview Access Analyzer findings for your
resource before deploying resource permissions.

# Arguments
- `analyzer_arn`: The ARN of the account analyzer used to generate the access preview. You
  can only create an access preview for analyzers with an Account type and Active status.
- `configurations`: Access control configuration for your resource that is used to generate
  the access preview. The access preview includes findings for external access allowed to the
  resource with the proposed access control configuration. The configuration must contain
  exactly one element.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: A client token.
"""
create_access_preview(analyzerArn, configurations; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("PUT", "/access-preview", Dict{String, Any}("analyzerArn"=>analyzerArn, "configurations"=>configurations, "clientToken"=>string(uuid4())); aws_config=aws_config)
create_access_preview(analyzerArn, configurations, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("PUT", "/access-preview", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn, "configurations"=>configurations, "clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    create_analyzer(analyzer_name, type)
    create_analyzer(analyzer_name, type, params::Dict{String,<:Any})

Creates an analyzer for your account.

# Arguments
- `analyzer_name`: The name of the analyzer to create.
- `type`: The type of analyzer to create. Only ACCOUNT and ORGANIZATION analyzers are
  supported. You can create only one analyzer per account per Region. You can create up to 5
  analyzers per organization per Region.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"archiveRules"`: Specifies the archive rules to add for the analyzer. Archive rules
  automatically archive findings that meet the criteria you define for the rule.
- `"clientToken"`: A client token.
- `"tags"`: The tags to apply to the analyzer.
"""
create_analyzer(analyzerName, type; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("PUT", "/analyzer", Dict{String, Any}("analyzerName"=>analyzerName, "type"=>type, "clientToken"=>string(uuid4())); aws_config=aws_config)
create_analyzer(analyzerName, type, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("PUT", "/analyzer", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerName"=>analyzerName, "type"=>type, "clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    create_archive_rule(analyzer_name, filter, rule_name)
    create_archive_rule(analyzer_name, filter, rule_name, params::Dict{String,<:Any})

Creates an archive rule for the specified analyzer. Archive rules automatically archive new
findings that meet the criteria you define when you create the rule.

# Arguments
- `analyzer_name`: The name of the created analyzer.
- `filter`: The criteria for the rule.
- `rule_name`: The name of the rule to create.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: A client token.
"""
create_archive_rule(analyzerName, filter, ruleName; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("PUT", "/analyzer/$(analyzerName)/archive-rule", Dict{String, Any}("filter"=>filter, "ruleName"=>ruleName, "clientToken"=>string(uuid4())); aws_config=aws_config)
create_archive_rule(analyzerName, filter, ruleName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("PUT", "/analyzer/$(analyzerName)/archive-rule", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("filter"=>filter, "ruleName"=>ruleName, "clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    delete_analyzer(analyzer_name)
    delete_analyzer(analyzer_name, params::Dict{String,<:Any})

Deletes the specified analyzer. When you delete an analyzer, Access Analyzer is disabled
for the account or organization in the current or specific Region. All findings that were
generated by the analyzer are deleted. You cannot undo this action.

# Arguments
- `analyzer_name`: The name of the analyzer to delete.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: A client token.
"""
delete_analyzer(analyzerName; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("DELETE", "/analyzer/$(analyzerName)", Dict{String, Any}("clientToken"=>string(uuid4())); aws_config=aws_config)
delete_analyzer(analyzerName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("DELETE", "/analyzer/$(analyzerName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    delete_archive_rule(analyzer_name, rule_name)
    delete_archive_rule(analyzer_name, rule_name, params::Dict{String,<:Any})

Deletes the specified archive rule.

# Arguments
- `analyzer_name`: The name of the analyzer that associated with the archive rule to delete.
- `rule_name`: The name of the rule to delete.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: A client token.
"""
delete_archive_rule(analyzerName, ruleName; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("DELETE", "/analyzer/$(analyzerName)/archive-rule/$(ruleName)", Dict{String, Any}("clientToken"=>string(uuid4())); aws_config=aws_config)
delete_archive_rule(analyzerName, ruleName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("DELETE", "/analyzer/$(analyzerName)/archive-rule/$(ruleName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    get_access_preview(access_preview_id, analyzer_arn)
    get_access_preview(access_preview_id, analyzer_arn, params::Dict{String,<:Any})

Retrieves information about an access preview for the specified analyzer.

# Arguments
- `access_preview_id`: The unique ID for the access preview.
- `analyzer_arn`: The ARN of the analyzer used to generate the access preview.

"""
get_access_preview(accessPreviewId, analyzerArn; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/access-preview/$(accessPreviewId)", Dict{String, Any}("analyzerArn"=>analyzerArn); aws_config=aws_config)
get_access_preview(accessPreviewId, analyzerArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/access-preview/$(accessPreviewId)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn), params)); aws_config=aws_config)

"""
    get_analyzed_resource(analyzer_arn, resource_arn)
    get_analyzed_resource(analyzer_arn, resource_arn, params::Dict{String,<:Any})

Retrieves information about a resource that was analyzed.

# Arguments
- `analyzer_arn`: The ARN of the analyzer to retrieve information from.
- `resource_arn`: The ARN of the resource to retrieve information about.

"""
get_analyzed_resource(analyzerArn, resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzed-resource", Dict{String, Any}("analyzerArn"=>analyzerArn, "resourceArn"=>resourceArn); aws_config=aws_config)
get_analyzed_resource(analyzerArn, resourceArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzed-resource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn, "resourceArn"=>resourceArn), params)); aws_config=aws_config)

"""
    get_analyzer(analyzer_name)
    get_analyzer(analyzer_name, params::Dict{String,<:Any})

Retrieves information about the specified analyzer.

# Arguments
- `analyzer_name`: The name of the analyzer retrieved.

"""
get_analyzer(analyzerName; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzer/$(analyzerName)"; aws_config=aws_config)
get_analyzer(analyzerName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzer/$(analyzerName)", params; aws_config=aws_config)

"""
    get_archive_rule(analyzer_name, rule_name)
    get_archive_rule(analyzer_name, rule_name, params::Dict{String,<:Any})

Retrieves information about an archive rule. To learn about filter keys that you can use to
create an archive rule, see Access Analyzer filter keys in the IAM User Guide.

# Arguments
- `analyzer_name`: The name of the analyzer to retrieve rules from.
- `rule_name`: The name of the rule to retrieve.

"""
get_archive_rule(analyzerName, ruleName; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzer/$(analyzerName)/archive-rule/$(ruleName)"; aws_config=aws_config)
get_archive_rule(analyzerName, ruleName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzer/$(analyzerName)/archive-rule/$(ruleName)", params; aws_config=aws_config)

"""
    get_finding(analyzer_arn, id)
    get_finding(analyzer_arn, id, params::Dict{String,<:Any})

Retrieves information about the specified finding.

# Arguments
- `analyzer_arn`: The ARN of the analyzer that generated the finding.
- `id`: The ID of the finding to retrieve.

"""
get_finding(analyzerArn, id; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/finding/$(id)", Dict{String, Any}("analyzerArn"=>analyzerArn); aws_config=aws_config)
get_finding(analyzerArn, id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/finding/$(id)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn), params)); aws_config=aws_config)

"""
    list_access_preview_findings(access_preview_id, analyzer_arn)
    list_access_preview_findings(access_preview_id, analyzer_arn, params::Dict{String,<:Any})

Retrieves a list of access preview findings generated by the specified access preview.

# Arguments
- `access_preview_id`: The unique ID for the access preview.
- `analyzer_arn`: The ARN of the analyzer used to generate the access.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"filter"`: Criteria to filter the returned findings.
- `"maxResults"`: The maximum number of results to return in the response.
- `"nextToken"`: A token used for pagination of results returned.
"""
list_access_preview_findings(accessPreviewId, analyzerArn; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("POST", "/access-preview/$(accessPreviewId)", Dict{String, Any}("analyzerArn"=>analyzerArn); aws_config=aws_config)
list_access_preview_findings(accessPreviewId, analyzerArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("POST", "/access-preview/$(accessPreviewId)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn), params)); aws_config=aws_config)

"""
    list_access_previews(analyzer_arn)
    list_access_previews(analyzer_arn, params::Dict{String,<:Any})

Retrieves a list of access previews for the specified analyzer.

# Arguments
- `analyzer_arn`: The ARN of the analyzer used to generate the access preview.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return in the response.
- `"nextToken"`: A token used for pagination of results returned.
"""
list_access_previews(analyzerArn; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/access-preview", Dict{String, Any}("analyzerArn"=>analyzerArn); aws_config=aws_config)
list_access_previews(analyzerArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/access-preview", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn), params)); aws_config=aws_config)

"""
    list_analyzed_resources(analyzer_arn)
    list_analyzed_resources(analyzer_arn, params::Dict{String,<:Any})

Retrieves a list of resources of the specified type that have been analyzed by the
specified analyzer..

# Arguments
- `analyzer_arn`: The ARN of the analyzer to retrieve a list of analyzed resources from.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return in the response.
- `"nextToken"`: A token used for pagination of results returned.
- `"resourceType"`: The type of resource.
"""
list_analyzed_resources(analyzerArn; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("POST", "/analyzed-resource", Dict{String, Any}("analyzerArn"=>analyzerArn); aws_config=aws_config)
list_analyzed_resources(analyzerArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("POST", "/analyzed-resource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn), params)); aws_config=aws_config)

"""
    list_analyzers()
    list_analyzers(params::Dict{String,<:Any})

Retrieves a list of analyzers.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return in the response.
- `"nextToken"`: A token used for pagination of results returned.
- `"type"`: The type of analyzer.
"""
list_analyzers(; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzer"; aws_config=aws_config)
list_analyzers(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzer", params; aws_config=aws_config)

"""
    list_archive_rules(analyzer_name)
    list_archive_rules(analyzer_name, params::Dict{String,<:Any})

Retrieves a list of archive rules created for the specified analyzer.

# Arguments
- `analyzer_name`: The name of the analyzer to retrieve rules from.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return in the request.
- `"nextToken"`: A token used for pagination of results returned.
"""
list_archive_rules(analyzerName; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzer/$(analyzerName)/archive-rule"; aws_config=aws_config)
list_archive_rules(analyzerName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzer/$(analyzerName)/archive-rule", params; aws_config=aws_config)

"""
    list_findings(analyzer_arn)
    list_findings(analyzer_arn, params::Dict{String,<:Any})

Retrieves a list of findings generated by the specified analyzer. To learn about filter
keys that you can use to retrieve a list of findings, see Access Analyzer filter keys in
the IAM User Guide.

# Arguments
- `analyzer_arn`: The ARN of the analyzer to retrieve findings from.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"filter"`: A filter to match for the findings to return.
- `"maxResults"`: The maximum number of results to return in the response.
- `"nextToken"`: A token used for pagination of results returned.
- `"sort"`: The sort order for the findings returned.
"""
list_findings(analyzerArn; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("POST", "/finding", Dict{String, Any}("analyzerArn"=>analyzerArn); aws_config=aws_config)
list_findings(analyzerArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("POST", "/finding", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn), params)); aws_config=aws_config)

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Retrieves a list of tags applied to the specified resource.

# Arguments
- `resource_arn`: The ARN of the resource to retrieve tags from.

"""
list_tags_for_resource(resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("GET", "/tags/$(resourceArn)", params; aws_config=aws_config)

"""
    start_resource_scan(analyzer_arn, resource_arn)
    start_resource_scan(analyzer_arn, resource_arn, params::Dict{String,<:Any})

Immediately starts a scan of the policies applied to the specified resource.

# Arguments
- `analyzer_arn`: The ARN of the analyzer to use to scan the policies applied to the
  specified resource.
- `resource_arn`: The ARN of the resource to scan.

"""
start_resource_scan(analyzerArn, resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("POST", "/resource/scan", Dict{String, Any}("analyzerArn"=>analyzerArn, "resourceArn"=>resourceArn); aws_config=aws_config)
start_resource_scan(analyzerArn, resourceArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("POST", "/resource/scan", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn, "resourceArn"=>resourceArn), params)); aws_config=aws_config)

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Adds a tag to the specified resource.

# Arguments
- `resource_arn`: The ARN of the resource to add the tag to.
- `tags`: The tags to add to the resource.

"""
tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes a tag from the specified resource.

# Arguments
- `resource_arn`: The ARN of the resource to remove the tag from.
- `tag_keys`: The key for the tag to add.

"""
untag_resource(resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)

"""
    update_archive_rule(analyzer_name, filter, rule_name)
    update_archive_rule(analyzer_name, filter, rule_name, params::Dict{String,<:Any})

Updates the criteria and values for the specified archive rule.

# Arguments
- `analyzer_name`: The name of the analyzer to update the archive rules for.
- `filter`: A filter to match for the rules to update. Only rules that match the filter are
  updated.
- `rule_name`: The name of the rule to update.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: A client token.
"""
update_archive_rule(analyzerName, filter, ruleName; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("PUT", "/analyzer/$(analyzerName)/archive-rule/$(ruleName)", Dict{String, Any}("filter"=>filter, "clientToken"=>string(uuid4())); aws_config=aws_config)
update_archive_rule(analyzerName, filter, ruleName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("PUT", "/analyzer/$(analyzerName)/archive-rule/$(ruleName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("filter"=>filter, "clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    update_findings(analyzer_arn, status)
    update_findings(analyzer_arn, status, params::Dict{String,<:Any})

Updates the status for the specified findings.

# Arguments
- `analyzer_arn`: The ARN of the analyzer that generated the findings to update.
- `status`: The state represents the action to take to update the finding Status. Use
  ARCHIVE to change an Active finding to an Archived finding. Use ACTIVE to change an
  Archived finding to an Active finding.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: A client token.
- `"ids"`: The IDs of the findings to update.
- `"resourceArn"`: The ARN of the resource identified in the finding.
"""
update_findings(analyzerArn, status; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("PUT", "/finding", Dict{String, Any}("analyzerArn"=>analyzerArn, "status"=>status, "clientToken"=>string(uuid4())); aws_config=aws_config)
update_findings(analyzerArn, status, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = accessanalyzer("PUT", "/finding", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn, "status"=>status, "clientToken"=>string(uuid4())), params)); aws_config=aws_config)


text:
##

import ACE, ACEbase
using LinearAlgebra, StaticArrays, Test, Printf, ACE.Testing
using ACE: evaluate, evaluate_d, evaluate_ed, Trig1pBasis, Trig

##

@info(" --------- Trig Tests --------- ")
@info("Test conversion to angle")
for ntest = 1:30
   θ = 2*π*rand() - π
   rr = [cos(θ), sin(θ)]
   print_tf(@test( ACE._theta(rr) ≈ θ ))
end
println()

##

@info("test gradient of coordinate transform")
for ntest = 1:30
   rr = randn(2) # [cos(θ), sin(θ)]
   G = rr -> ACE._theta(rr)
   dG = rr -> ACE._theta_ed(rr)[2] |> Vector
   print_tf(@test( all(
            ACEbase.Testing.fdtest(G, dG, rr; verbose=false)
         )))
end
println()

##

@info("test evaluation of trig basis")

bE = Trig(10)

for ntest = 1:30
   L = bE.L
   θ = π * rand()
   E = [ exp(im * l * θ) for l = -L:L ]
   rr = [cos(θ), sin(θ)]
   print_tf( @test evaluate(bE, rr) ≈ E )
end
println()

##

@info("test jacobian of trig basis ")
θ = π * rand()
rr = [cos(θ), sin(θ)]
E, dE = evaluate_ed(bE, rr)

for ntest = 1:30
   rr = randn(2)
   U = randn(length(E))
   F = rr -> ACE.contract(evaluate(bE, rr), U)
   dF = rr -> Vector(ACE.contract(evaluate_ed(bE, rr)[2], U))
   print_tf(@test( all(
          ACEbase.Testing.fdtest(F, dF, rr; verbose=false)
      )))
end
println()

##


text:
abstract type Field end
abstract type DiscreteField end

struct ScalarField <: Field end
struct VectorField <: Field end
struct MatrixField <: Field end

struct DiscretizedField <: Field

end


text:
using VisDyn
using Test

@testset "VisDyn.jl" begin
    # Write your tests here.
    plotTrajectory!(plot(),1:10)
end


text:
using HTTP

HTTP.serve() do request            # 1
  try                              # 2
    return HTTP.Response("Still lots of ToDos!")
  catch e
    return HTTP.Response(404, "Error: $e")
  end
end

# ::HTTP.Request


text:
# Functions for handling and executing actions on a state

const no_op = Action(Compound(Symbol("--"), []), @julog(true), @julog(and()))

"Get preconditions of an action as a list."
function get_preconditions(act::Action, args::Vector{<:Term};
                           converter::Function=flatten_conjs)
    subst = Subst(var => val for (var, val) in zip(act.args, args))
    precond = substitute(act.precond, subst)
    return converter(precond)
end

get_preconditions(act::Action; converter::Function=flatten_conjs) =
    converter(act.precond)

get_preconditions(act::Term, domain::Domain; kwargs...) =
    get_preconditions(domain.actions[act.name], get_args(act); kwargs...)

"Get effect term of an action with variables substituted by arguments."
function get_effect(act::Action, args::Vector{<:Term})
    subst = Subst(var => val for (var, val) in zip(act.args, args))
    return substitute(act.effect, subst)
end

get_effect(act::Term, domain::Domain) =
    get_effect(domain.actions[act.name], get_args(act))

"Global flag for whether to cache available actions."
const _use_available_cache = Ref(true)
"Cache of available actions for a given domain and state."
const _available_action_cache = Dict{UInt,Dict{UInt,Vector{Term}}}()

"Globally enable or disable available action cache."
use_available_action_cache!(val::Bool=true) = _use_available_cache[] = val

"Clear cache of available actions."
clear_available_action_cache!() =
    (map(empty!, values(_available_action_cache));
     empty!(_available_action_cache))
clear_available_action_cache!(domain::Domain) =
    empty!(_available_action_cache[objectid(domain)])

"""
    available(state, domain; use_cache)

Return the list of available actions in a given `state` and `domain`.
If `use_cache` is true, memoize the results in a global cache.
"""
function available(state::State, domain::Domain;
                   use_cache::Bool=_use_available_cache[])
    if use_cache # Look up actions in cache
        cache = get!(_available_action_cache, objectid(domain),
                     Dict{UInt,Vector{Term}}())
        state_hash = hash(state)
        if haskey(cache, state_hash) return copy(cache[state_hash]) end
    end
    # Ground all action definitions with arguments
    actions = Term[]
    for act in values(domain.actions)
        typecond = (@julog($ty(:v)) for (v, ty) in zip(act.args, act.types))
        # Include type conditions when necessary for correctness
        p = act.precond
        if has_fluent(p, domain) || has_axiom(p, domain) || has_quantifier(p)
            conds = prepend!(flatten_conjs(p), typecond)
        elseif domain.requirements[:typing]
            conds = append!(flatten_conjs(p), typecond)
        else
            conds = flatten_conjs(p)
        end
        # Find all substitutions that satisfy preconditions
        sat, subst = satisfy(conds, state, domain; mode=:all)
        if !sat continue end
        for s in subst
            args = [s[v] for v in act.args if v in keys(s)]
            if any(!is_ground(a) for a in args) continue end
            term = isempty(args) ? Const(act.name) : Compound(act.name, args)
            push!(actions, term)
        end
    end
    if use_cache cache[state_hash] = copy(actions) end
    return actions
end

"""
    available(act::Action, args, state, domain=nothing)
    available(act::Term, state, domain)

Check whether `act` can be executed in the given `state` and `domain`.
"""
function available(act::Action, args::Vector{<:Term}, state::State,
                   domain::Union{Domain,Nothing}=nothing)
    if any(!is_ground(a) for a in args)
       error("Not all arguments are ground.")
    end
    subst = Subst(var => val for (var, val) in zip(act.args, args))
    # Construct type conditions of the form "type(val)"
    typecond = (all(ty == :object for ty in act.types) ? Term[] :
               [@julog($ty(:v)) for (v, ty) in zip(args, act.types)])
    # Check whether preconditions hold
    precond = substitute(act.precond, subst)
    conds = has_fluent(precond, state) || has_quantifier(precond) ?
        [typecond; precond] : [precond; typecond]
    sat, _ = satisfy(conds, state, domain)
    return sat
end

available(act::Term, state::State, domain::Domain) =
    available(domain.actions[act.name], act.args, state, domain)

"Global flag for whether to cache relevant actions."
const _use_relevant_cache = Ref(true)
"Cache of relevant actions for a given domain and state."
const _relevant_action_cache = Dict{UInt,Dict{UInt,Vector{Term}}}()

"Globally enable or disable available action cache."
use_relevant_action_cache!(val::Bool=true) = _use_relevant_cache[] = val

"Clear cache of relevant actions."
clear_relevant_action_cache!() =
    (map(empty!, values(_relevant_action_cache));
     empty!(_relevant_action_cache))
clear_relevant_action_cache!(domain::Domain) =
    empty!(_relevant_action_cache[objectid(domain)])

"""
    relevant(state, domain; strict=false, use_cache)

Return the list of actions relevant to achieving a `state` in a given `domain`.
If `strict` is true, check that all added facts are true in `state`.
If `use_cache` is true, memoize the results in a global cache.
"""
function relevant(state::State, domain::Domain;
                  strict::Bool=false, use_cache::Bool=_use_relevant_cache[])
    if use_cache # Look up actions in cache
        cache = get!(_relevant_action_cache, hash(strict, objectid(domain)),
                     Dict{UInt,Vector{Term}}())
        state_hash = hash(state)
        if haskey(cache, state_hash) return copy(cache[state_hash]) end
    end
    actions = Term[]
    for act in values(domain.actions)
        # Compute postconditions from the action's effect
        diff = effect_diff(act.effect)
        addcond = strict ? diff.add : [Compound(:or, diff.add)]
        delcond = [@julog(not(:t)) for t in diff.del]
        typecond = [@julog($ty(:v)) for (v, ty) in zip(act.args, act.types)]
        # Include type conditions when necessary for correctness
        if any(has_fluent(c, domain) ||
               has_quantifier(c) for c in [addcond; delcond])
            conds = [typecond; addcond; delcond]
        else
            conds = [addcond; typecond; delcond]
        end
        # Find all substitutions that satisfy the postconditions
        sat, subst = satisfy(conds, state, domain; mode=:all)
        if !sat continue end
        for s in subst
            args = [get(s, var, var) for var in act.args]
            if any(!is_ground(a) for a in args) continue end
            term = isempty(args) ? Const(act.name) : Compound(act.name, args)
            push!(actions, term)
        end
    end
    if use_cache cache[state_hash] = copy(actions) end
    return actions
end

"""
    relevant(act::Action, args, state, domain=nothing; strict=false)
    relevant(act::Term, state, domain; strict=false)

Check whether `act` is relevant (can lead to) a `state` in the given `domain`.
If `strict` is true, check that all added facts are true in `state`.
"""
function relevant(act::Action, args::Vector{<:Term}, state::State,
                  domain::Union{Domain,Nothing}=nothing; strict::Bool=false)
   if any(!is_ground(a) for a in args)
       error("Not all arguments are ground.")
   end
   subst = Subst(var => val for (var, val) in zip(act.args, args))
   # Compute postconditions from the action's effect
   diff = effect_diff(substitute(act.effect, subst))
   postcond = Term[strict ? diff.add : Compound(:or, diff.add);
                   [@julog(not(:t)) for t in diff.del]]
   # Construct type conditions of the form "type(val)"
   typecond = (all(ty == :object for ty in act.types) ? Term[] :
               [@julog($ty(:v)) for (v, ty) in zip(args, act.types)])
   # Check whether postconditions hold
   sat, _ = satisfy([postcond; typecond], state, domain)
   return sat
end

relevant(act::Term, state::State, domain::Domain; kwargs...) =
    relevant(domain.actions[act.name], act.args, state, domain; kwargs...)

"""
    execute(act::Action, args, state, domain=nothing; kwargs...)
    execute(act::Term, state, domain; kwargs...)

Execute the action `act` on the given `state`. If `act` is an `Action`
definition, `args` must be supplied for the action's parameters. If `act` is
a `Term`, parameters will be extracted from its arguments, but `domain` must
be supplied.

Returns the resulting state by default, but this can modified by keyword
arguments `as_dist=true` (returns a distribution), `as_diff=true` (returns a
`Diff` between states), or both (returns a distribution over `Diff`s).
By default, also `check`s that action preconditions hold, failing with
an error if `fail_mode=:error`, or treating the action as a `no_op` if
`fail_mode=:no_op`.
"""
function execute(act::Action, args::Vector{<:Term}, state::State,
                 domain::Union{Domain,Nothing}=nothing;
                 as_dist::Bool=false, as_diff::Bool=false,
                 check::Bool=true, fail_mode::Symbol=:error)
    # Check whether references resolve and preconditions hold
    if check && !available(act, args, state, domain)
        if fail_mode == :no_op return as_diff ? no_effect(as_dist) : state end
        error("Precondition $(act.precond) does not hold.") # Error by default
    end
    # Substitute arguments and preconditions
    # TODO : Check for non-ground terms outside of quantified formulae
    subst = Subst(var => val for (var, val) in zip(act.args, args))
    effect = substitute(act.effect, subst)
    # Compute effects in the appropriate form
    if as_dist
        # Compute categorical distribution over differences
        diff = effect_dist(effect, state, domain)
    else
        # Sample a possible difference
        diff = effect_diff(effect, state, domain)
    end
    # Return either the difference or the updated state
    return as_diff ? diff : update(state, diff)
end

function execute(act::Term, state::State, domain::Domain; options...)
    if act.name in keys(domain.actions)
        act_def, act_args = domain.actions[act.name], get_args(act)
        execute(act_def, act_args, state, domain; options...)
    elseif act.name == Symbol("--")
        execute(no_op, Term[], state, domain; options...)
    else
        error("Unknown action: $act")
    end
end

function execute(actions::Vector{<:Term}, state::State, domain::Domain;
                 as_dist::Bool=false, as_diff::Bool=false, options...)
    state = copy(state)
    for act in actions
        diff = execute(domain.actions[act.name], get_args(act), state, domain;
                       as_dist=as_dist, as_diff=true, options...)
        update!(state, diff)
    end
    # Return either the difference or the final state
    return as_diff ? diff : state
end

"Execute a list of actions in sequence on a state."
execseq(actions::Vector{<:Term}, state::State, domain::Domain; options...) =
    execute(actions, state, domain; options...)

function execute(actions::Set{<:Term}, state::State, domain::Domain;
                 as_dist::Bool=false, as_diff::Bool=false, options...)
    diffs = [execute(domain.actions[act.name], get_args(act), state, domain;
                     as_dist=as_dist, as_diff=true, options...)
             for act in actions]
    diff = combine(diffs...)
    # Return either the difference or the updated state
    return as_diff ? diff : update(state, diff)
end

"Execute a set of actions in parallel on a state."
execpar(actions::Set{<:Term}, state::State, domain::Domain; options...) =
    execute(actions, state, domain; options...)
execpar(actions::Vector{<:Term}, state::State, domain::Domain; options...) =
    execute(Set(actions), state, domain; options...)

"""
    regress(act::Action, args, state, domain=nothing; kwargs...)
    regress(act::Term, state, domain; kwargs...)

Return the predecessor(s) that would lead to `state` upon executing `act`.
If `act` is an `Action` definition, `args` must be supplied for the action's
parameters. If `act` is a `Term`, parameters will be extracted from its
arguments, but `domain` must be supplied.
"""
function regress(act::Action, args::Vector{<:Term}, state::State,
                 domain::Union{Domain,Nothing}=nothing; as_diff::Bool=false,
                 check::Bool=true, fail_mode::Symbol=:error)
    # Check whether action is relevant
    if check && !relevant(act, args, state, domain)
        if fail_mode == :no_op return as_diff ? Diff() : state end
        error("Effect $(act.effect) is not relevant.") # Error by default
    end
    subst = Subst(var => val for (var, val) in zip(act.args, args))
    precond = substitute(act.precond, subst)
    effect = substitute(act.effect, subst)
    # Compute regression difference as Precond - Additions
    # TODO: Handle conditional effects, disjunctive preconditions, etc.
    pre_diff = precond_diff(precond, state)
    eff_diff = effect_diff(effect, state)
    append!(pre_diff.del, eff_diff.add)
    return as_diff ? pre_diff : update(state, pre_diff)
end

function regress(act::Term, state::State, domain::Domain; options...)
    if act.name in keys(domain.actions)
        act_def, act_args = domain.actions[act.name], get_args(act)
        regress(act_def, act_args, state, domain; options...)
    elseif act.name == Symbol("--")
        regress(no_op, Term[], state, domain; options...)
    else
        error("Unknown action: $act")
    end
end


text:
using MatrixMarket
using HSL
using LinearOperators

K = MatrixMarket.mmread("K_0.mtx")
rhs = readdlm("rhs_0.rhs")

# Version 1: simple
# x = ma97_solve(K, rhs)

# Version 2: intermediate
# LBL = ma97_factorize(K)
# x = ma97_solve(LBL, b)  # or x = LBL \ b

# Version 3: finer control
LBL = Ma97(K)
ma97_factorize(LBL)
x = ma97_solve(LBL, rhs)  # or x = LBL \ b

# Only the lower triangle appears in K.
# We use a linear operator to compute products with the symmetric matrix.
H = opHermitian(K)
@printf("Relative residual: %7.1e\n", norm(H * x[:] - rhs) / norm(rhs))


text:
@test_throws SystemError("mkdir", 2, "/foo/bar/app") matte_example("hello_world", "/foo/bar/app")
new_matte_app("test_new_app", activate = false) == "test_new_app"
@test basename(pwd()) == "test_new_app"
@test isfile("app.jl")
cd("../")
rm("test_new_app", recursive = true)

matte_example("hello_world", "hello_world", activate = false) == "hello_world"
@test basename(pwd()) == "hello_world"
@test isfile("app.jl")
cd("../")
rm("hello_world", recursive = true)


text:
@snn_kw struct IZParameter{FT=Float32}
    a::FT = 0.01
    b::FT = 0.2
    c::FT = -65
    d::FT = 2
end

@snn_kw mutable struct IZ{VFT=Vector{Float32},VBT=Vector{Bool}}
    param::IZParameter = IZParameter()
    N::Int32 = 100
    v::VFT = fill(-65.0, N)
    u::VFT = param.b * v
    fire::VBT = zeros(Bool, N)
    I::VFT = zeros(N)
    records::Dict = Dict()
end

"""
[Izhikevich Neuron](https://www.izhikevich.org/publications/spikes.htm)
"""
IZ

function integrate!(p::IZ, param::IZParameter, dt::Float32)
    @unpack N, v, u, fire, I = p
    @unpack a, b, c, d = param
    @inbounds for i = 1:N
        v[i] += 0.5f0dt * (0.04f0v[i]^2 + 5f0v[i] + 140f0 - u[i] + I[i])
        v[i] += 0.5f0dt * (0.04f0v[i]^2 + 5f0v[i] + 140f0 - u[i] + I[i])
        u[i] += dt * (a * (b * v[i] - u[i]))
    end

    @inbounds for i = 1:N
        fire[i] = v[i] > 0f0#30f0
        v[i] = ifelse(fire[i], c, v[i])
        u[i] += ifelse(fire[i], d, 0f0)
    end

end
#=
@inbounds for m = 1:N
    vT = v[m]+ 0.5f0dt * (k*(v[m] - vr)*(v[m] - vt)-u[m] + Iext[m])/C;
    v[m+1] = vT + 0.5f0dt  * (k*(v[m] - vr)*(v[m] - vt)-u[m] + Iext[m])/C;
    u[m+1] = u[m] + dt * a*(b*(v[m]-vr)-u[m]);
    u[m+1] = u[m] + dt * a*(b*(v[m+1]-vr)-u[m]);
    if v[m+1]>= vPeak:# % a spike is fired!
        v[m] = vPeak;# % padding the spike amplitude
        v[m+1] = c;# % membrane voltage reset
        u[m+1] = u[m+1] + d;# % recovery variable update
    end
end
=#


text:
using Test
using JPLEphemeris

using AstroBase

const AU = astronomical_unit(Float64, KiloMeter())[1]

const SPK_URL = Dict{Int, String}(
    430 => "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de430.bsp",
    405 => "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/a_old_versions/de405.bsp",
)
const TEST_URL = Dict{Int, String}(
    430 => "ftp://ssd.jpl.nasa.gov/pub/eph/planets/test-data/430/testpo.430",
    405 => "ftp://ssd.jpl.nasa.gov/pub/eph/planets/test-data/testpo.405",
)

path = abspath(joinpath(dirname(@__FILE__), "..", "deps"))

for denum in (430, 405)
    if !isdir(path)
        mkdir(path)
    end
    if !isfile("$path/de$denum.bsp")
        download(SPK_URL[denum], "$path/de$denum.bsp")
    end
    if !isfile("$path/testpo.$denum")
        download(TEST_URL[denum], "$path/testpo.$denum")
    end
end

@testset "JPLEphemeris" begin
    include("basic.jl")
    include("kernels.jl")
end


text:
# Load in the dependencies
using InfiniteOpt, JuMP, MathOptInterface, Distributions, Random,
FastGaussQuadrature, DataStructures

# load the test module
using Test

# Define convenient aliases
const IC = InfiniteOpt.Collections
const MT = InfiniteOpt.MeasureToolbox
const IOTO = InfiniteOpt.TranscriptionOpt
const JuMPC = JuMP.Containers
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
const MOIUC = MathOptInterface.Utilities.CleverDicts
const FGQ = FastGaussQuadrature
const IOMT = InfiniteOpt.MeasureToolbox

# Load in testing utilities
include("utilities.jl")

# Run unit tests
println("----------------------------------------------------------------------------")
println("---------------------------------UNIT TESTS---------------------------------")
println("----------------------------------------------------------------------------")
@time @testset "Collections" begin
    include("Collections/VectorTuple.jl")
    include("Collections/DualDict.jl")
end
println("")
@time @testset "Datatypes" begin include("datatypes.jl") end
println("")
@time @testset "Utilities" begin include("utility_tests.jl") end
println("")
@time @testset "Infinite Set Methods" begin include("infinite_sets.jl") end
println("")
@time @testset "General Variable Methods" begin
    include("general_variables.jl")
end
println("")
@time @testset "Optimizer Setup Methods" begin include("optimizer_setup.jl") end
println("")
@time @testset "Macro Utilities" begin include("macro_utilities.jl") end
println("")
@time @testset "Parameter Methods" begin
   @testset "Scalar" begin include("scalar_parameters.jl") end
   @testset "Array" begin include("array_parameters.jl") end
end
println("")
@time @testset "Variable Methods" begin
   @testset "Infinite Variables" begin include("infinite_variables.jl") end
   @testset "Point Variables" begin include("point_variables.jl") end
   @testset "Hold Variables" begin include("hold_variables.jl") end
   @testset "Info Constraints" begin include("variable_info.jl") end
   @testset "Reduced Variables" begin include("reduced_variables.jl") end
end
println("")
@time @testset "Derivative Methods" begin include("derivatives.jl") end
println("")
@time @testset "Operators" begin include("operators.jl") end
println("")
@time @testset "Expression Methods" begin include("expressions.jl") end
println("")
@time @testset "Macro Expressions" begin include("macro_expressions.jl") end
println("")
@time @testset "Measure Methods" begin include("measures.jl") end
println("")
@time @testset "Measure Toolbox Methods" begin
    @testset "Integrals" begin include("MeasureToolbox/integrals.jl") end
    @testset "Expectations" begin include("MeasureToolbox/expectations.jl") end
    @testset "Support Sums" begin include("MeasureToolbox/support_sums.jl") end
end
println("")
@time @testset "Objective Methods" begin include("objective.jl") end
println("")
@time @testset "Constraint Methods" begin include("constraints.jl") end
println("")
@time @testset "Printing Methods" begin include("show.jl") end
println("")
@time @testset "Deletion Methods" begin include("deletion.jl") end
println("")
@time @testset "Expansion Methods" begin include("measure_expansions.jl") end
println("")
@time @testset "Derivative Evaluation" begin include("derivative_evaluation.jl") end
println("")
@time @testset "TranscriptionOpt" begin
    @testset "Model" begin include("TranscriptionOpt/model.jl") end
    @testset "Measures" begin include("TranscriptionOpt/measure.jl") end
    @testset "Transcribe" begin include("TranscriptionOpt/transcribe.jl") end
    @testset "Optimize" begin include("TranscriptionOpt/optimize.jl") end
end
println("")
@time @testset "Solution Methods" begin include("optimizer.jl") end
println("")
@time @testset "Solution Queries" begin include("results.jl") end
println("")
@time @testset "Extensions" begin include("extensions.jl") end
println("")
println("----------------------------------------------------------------------------")
println("-----------------------------TESTING COMPLETE!------------------------------")
println("----------------------------------------------------------------------------")


text:
export set_inputs, get_outputs, process

function reset(ind::NEATInd)
    for n in ind.neurons
        n.processed = false
    end
end

function set_inputs(ind::NEATInd, inputs::Array{Float64})
    # assume length(inputs) == ind.n_in
    for i in eachindex(inputs)
        ind.neurons[i].output = inputs[i]
        ind.neurons[i].processed = true
    end
end

function get_outputs(ind::NEATInd)
    outputs = zeros(ind.n_out)
    for i in 1:ind.n_out
        outputs[i] = ind.neurons[ind.n_in + i].output
    end
    outputs
end

function activation(x::Float64)
    1 / (1 + exp(-x))
end

function process_neuron(ind::NEATInd, neuron::Int)
    # TODO: recurrent nodes based on genetic order?
    for c in ind.connections
        if c.enabled && c.out_node == neuron
            if !ind.neurons[c.in_node].processed && c.in_node < c.out_node
                process_neuron(ind, c.in_node)
            end
            ind.neurons[neuron].input += ind.neurons[c.in_node].output * c.weight
        end
    end
    # TODO: bias neurons
    ind.neurons[neuron].output = activation(ind.neurons[neuron].input)
    ind.neurons[neuron].processed = true
end

function run_neurons(ind::NEATInd)
    for i in 1:ind.n_hidden
        process_neuron(ind, ind.n_in + ind.n_out + i)
    end
    for i in 1:ind.n_out
        process_neuron(ind, ind.n_in + i)
    end
end

function process(ind::NEATInd, inputs::Array{Float64})
    reset(ind)
    set_inputs(ind, inputs)
    run_neurons(ind)
    get_outputs(ind)
end


text:
using NCDatasets
using Formatting

fn1 = "clim_LENS_B1850C5CN_005_gx3v7_TEMP.nc"
fn2 = "clim_LENS_B1850C5CN_005_gx3v7_zMLMML_TEMP.nc"
varname = "TEMP"

idx_i = 45   # longitude idx
idx_j = 10   # latitude  idx

Dataset(fn1, "r") do ds

    global Nx = ds.dim["Nx"]
    global Ny = ds.dim["Nx"]
    global Nz_old = ds.dim["Nz"]

    global lat = ds["lat"][:]
    global lon = ds["lon"][:]
    global zs_old = (ds["z_w_top"][:] + ds["z_w_bot"][:] ) / 2 / 100.0

    global data_old = replace(ds[varname][idx_i, idx_j, :, 1], missing=>NaN)

end

Dataset(fn2, "r") do ds

    global Nz_new = ds.dim["Nz"]

    global zs_new = ( ds["zs"][1:end-1] +  ds["zs"][2:end] ) / 2

    global data_new = replace(ds[varname][idx_i, idx_j, :, 1], missing=>NaN)

end

print("Loading PyPlot...")
using PyPlot
println("done.")

plt[:figure]()
plt[:plot](zs_old, data_old, "b-",  label="old data")
plt[:plot](zs_new, data_new, "r--", label="new data")
plt[:legend]()

plt[:title](format("(lat, lon) = ({:.2f}, {:.2f})", lat[idx_i, idx_j], lon[idx_i, idx_j]))
plt[:show]()

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM