text:
using Distributions
#using Distributions: check_args

export NegativeBinomialConvolution, pmf, cdf

immutable NegativeBinomialConvolution <: DiscreteUnivariateDistribution
    rs::Array{Float64,1}
    ps::Array{Float64,1}
    logrs::Array{Float64,1}
    logps::Array{Float64,1}
    logqs::Array{Float64,1}
    mean::Float64
    stddev::Float64

    """
    This is based on the paper:
    "On the sums of compound negative binomial and gamma random variables" by Vallaisamy, 2009
    http://projecteuclid.org/euclid.jap/1238592129
    """
    function NegativeBinomialConvolution(rs::Array{Float64,1}, ps::Array{Float64,1})
        # @check_args(NegativeBinomialConvolution, minimum(rs) > zero(rs[1]))
        # @Distributions.check_args(NegativeBinomialConvolution, minimum(ps) > zero(ps[1]))
        # @Distributions.check_args(NegativeBinomialConvolution, Base.maximum(ps) < one(ps[1]))

        # just add the means and variances to get the convolution mean and variance
        meanVal = sum([(ps[i]*rs[i])/(1-ps[i]) for i in 1:length(ps)])
        varVal = sum([(ps[i]*rs[i])/(1-ps[i])^2 for i in 1:length(ps)])

        new(rs, ps, log(rs), log(ps), log(1 - ps), meanVal, sqrt(varVal))
    end
    NegativeBinomialConvolution(rs::Array{Real,1}, ps::Array{Real,1}) = NegativeBinomialConvolution(convert(Array{Float64,1}, rs), convert(Array{Float64,1}, ps))
    NegativeBinomialConvolution(rs::Array{Int64,1}, ps::Array{Float64,1}) = NegativeBinomialConvolution(convert(Array{Float64,1}, rs), convert(Array{Float64,1}, ps))
end

@Distributions.distr_support NegativeBinomialConvolution 0 Inf

params(d::NegativeBinomialConvolution) = (d.rs, d.ps)

"Call f on all possible ways n non-negative integers can sum to equal s."
function allsums_apply(f::Function, n::Int, s::Int)
    if n == 1
        f([s])
    elseif n == 2
        for i in 0:s
            f(Int64[i,s-i])
        end
    elseif n == 3
        for i in 0:s
            for j in 0:s
                if i + j > s break end
                f(Int64[i,j,s-i-j])
            end
        end
    elseif n == 4
        for i in 0:s
            for j in 0:s
                if i + j > s break end
                for k in 0:s
                    if i + j + k > s break end
                    f(Int64[i,j,k,s-i-j-k])
                end
            end
        end
    elseif n == 5
        for i in 0:s
            for j in 0:s
                if i + j > s break end
                for k in 0:s
                    if i + j + k > s break end
                    for l in 0:s
                        if i + j + k + l > s break end
                        f(Int64[i,j,k,l,s-i-j-k-l])
                    end
                end
            end
        end
    else
        error("n = $n is not yet implemented.")
    end
end

"Compute the exact PMF at the given point with a sum."
function Distributions.pdf(d::NegativeBinomialConvolution, s::Int)
    exp(Distributions.logpdf(d::NegativeBinomialConvolution, s::Int))
end

function Distributions.logpdf(d::NegativeBinomialConvolution, s::Int)
    total = 0.0

    vals = Float64[]
    function inner_term(m)
        val = 0.0
        for j in 1:length(d.ps)
            binVal = lgamma(d.rs[j] + m[j]) - lgamma(m[j]+1) - lgamma(d.rs[j])
            val += binVal + d.rs[j]*d.logps[j] + m[j]*d.logqs[j]
        end
        push!(vals, val)
    end
    allsums_apply(inner_term, length(d.ps), s)

    maxVal = Base.maximum(vals)
    maxVal + log(sum(map(x->exp(x-maxVal), vals)))
end

function Distributions.cdf(d::NegativeBinomialConvolution, s::Int)
    total = 0.0
    for i in 0:s
        total += pdf(d, i)
    end
    total
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     cell_metadata_json: true
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.11.3
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# # Transport Streamfunction And Divergent Component
#
# Here we apply such [Helmholtz Decomposition](https://en.wikipedia.org/wiki/Helmholtz_decomposition) to vertically integrated transports defined over a Global Ocean model (`C-grid`).
#
# 1. read vertically integrated transport from file
# 2. compute its convergence and apply land mask
# 3. decompose into rotational and divergent components
# 4. derive streamfunction from the rotational component

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ### Read Transports & Grid From File
#
# 1. pre-requisites
# 2. read variables

# + {"slideshow": {"slide_type": "subslide"}, "cell_style": "center"}
using MeshArrays, Plots, Statistics
pth=dirname(pathof(MeshArrays))
include(joinpath(pth,"../examples/Plots.jl"))
include("helper_functions.jl")

pth=MeshArrays.GRID_LLC90
γ=GridSpec("LatLonCap",pth)
Γ=GridLoad(γ;option="full")
(Tx,Ty,τx,τy,η)=trsp_read("LatLonCap",pth);

# +
μ =Γ.hFacC[:,1]
μ[findall(μ.>0.0)].=1.0
μ[findall(μ.==0.0)].=NaN

lon=[i for i=-179.:2.0:179., j=-89.:2.0:89.]
lat=[j for i=-179.:2.0:179., j=-89.:2.0:89.]

#(f,i,j,w)=InterpolationFactors(Γ,vec(lon),vec(lat))
#λ=(lon=lon,lat=lat,f=f,i=i,j=j,w=w);
#df = DataFrame(f=λ.f[:], i=λ.i[:], j=λ.j[:], w=Float32.(λ.w[:]));
#CSV.write("interp_coeffs.csv", df)

df=DataFrame(CSV.File("interp_coeffs.csv"))
λ=(f=reshape(df.f,length(lon[:]),4), i=reshape(df.i,length(lon[:]),4),
    j=reshape(df.j,length(lon[:]),4), w=reshape(df.w,length(lon[:]),4));

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ### Helmholtz Decomposition
#
# 1. convergence & land mask
# 2. compute scalar potential
# 3. subtract divergent component
# 4. compute vector potential / streamfunction
#
# _note: masking avoids singularities related to isolated canyons_

# + {"slideshow": {"slide_type": "subslide"}}
#convergence & land mask
TrspCon=μ.*convergence(Tx,Ty)

#scalar potential
TrspPot=ScalarPotential(TrspCon)

#Divergent transport component
(TxD,TyD)=gradient(TrspPot,Γ)
TxD=TxD.*Γ.DXC
TyD=TyD.*Γ.DYC

#Rotational transport component
TxR = Tx-TxD
TyR = Ty-TyD

#vector Potential
TrspPsi=VectorPotential(TxR,TyR,Γ);

# + [markdown] {"slideshow": {"slide_type": "skip"}}
# ### Verify The Results Consistency
#
# _Next plot below normally shows that TrspCon-tmpCon is negligible compared with TrspCon_

# + {"slideshow": {"slide_type": "skip"}, "cell_style": "center"}
tmpCon=convergence(TxD,TyD)
tmp1=TrspCon[3]
tmp2=tmp1[findall(isfinite.(tmp1))]
errCon=1/sqrt(mean(tmp2.^2)).*(tmpCon[3]-TrspCon[3])
errCon[findall(isnan.(errCon))].=0.0
extrema(errCon)

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ### Map Out Streamfunction And Scalar Potential
#
# _In plots below we interpolate MeshArrays to a regular global grid_

# + {"slideshow": {"slide_type": "subslide"}, "cell_style": "center"}
TrspPsiI=Interpolate(1e-6*μ*TrspPsi,λ.f,λ.i,λ.j,λ.w)
contourf(vec(lon[:,1]),vec(lat[1,:]),
    TrspPsiI,title="Streamfunction",clims=(-50,50))

# + {"slideshow": {"slide_type": "subslide"}, "cell_style": "center"}
TrspPotI=Interpolate(1e-6*μ*TrspPot,λ.f,λ.i,λ.j,λ.w)
contourf(vec(lon[:,1]),vec(lat[1,:]),
        TrspPotI,title="Scalar Potential",clims=(-0.4,0.4))
# -

ext: jl

max_stars_repo_name: gaelforget/Notebooks

max_issues_repo_name: gaelforget/Notebooks

max_forks_repo_name: gaelforget/Notebooks


text:
using Econometrics, Plots, Statistics, Calculus, LinearAlgebra, DelimitedFiles

include("SVlib.jl") # library of functions for the SV model

# proposals for MCMC
# uniform random walk in one dimension
function proposal1(current, tuning)
    trial = copy(current)
    i = rand(1:size(trial,1))
    trial[i] += tuning[i].*randn()
    return trial
end
# MVN random walk, or occasional draw from prior
function proposal2(current, cholV)
    current + cholV'*randn(size(current))
end

function main()
# setup
θtrue = TrueParameters() # true param values, on param space
lb, ub = PriorSupport()

# generate a new data set, or load the example data
n = 500
burnin = 100
mcreps = 100
results = zeros(mcreps,22)
Threads.@threads for rep = 1:mcreps
    y = SVmodel(TrueParameters(),n, burnin)
    m0 =  auxstat(y)
    # Estimation by two-step MSM
    S = 10 # number of simulation reps
    shocks_u = randn(n+burnin,S) # used fixed shocks to control chatter
    shocks_e = randn(n+burnin,S)
    # define the moments
    ms = θ -> m0' .- auxstat(θ, shocks_e, shocks_u)  # moment contributions: sample stats minus simulated stats
    # first step
    W = eye(size(m0,1))
    m = θ -> vec(mean(ms(θ),dims=1)) # use average of moment contributions
    obj = θ -> InSupport(θ) ? m(θ)'*W*m(θ) : Inf
    θhat, objvalue, converged, details = samin(obj, θtrue, lb, ub; ns = 5, verbosity = 0, rt = 0.5)
    # Estimate the weight matrix using more draws
    W = inv((1.0 + 1.0/S)*cov(m0' .- auxstat(θhat, randn(n+burnin, 500), randn(n+burnin, 500))))
    # second step
    obj = θ -> InSupport(θ) ? m(θ)'*W*m(θ) : Inf
    θhat, objvalue, converged, details = samin(obj, θhat, lb, ub; ns = 5, verbosity = 0, rt = 0.5)
    # compute the estimated standard errors and CIs
    D = (Calculus.jacobian(m, vec(θhat), :central))
    se = sqrt.(diag(inv(D'*W*D)))
    # now do MCMC
    prior = θ -> Prior(θ)
    burnin = 100
    ChainLength = 1000
    # initial proposal moves one at a time
    tuning = 0.2/sqrt(12.0)*(ub-lb) # two tenths of a standard. dev. of prior
    proposal = θ -> proposal1(θ, tuning)
    obj = θ -> InSupport(θ) ? -m(θ)'*W*m(θ) : -Inf
    chain = mcmc(θhat, ChainLength, burnin, prior, obj, proposal, false)
    # now use a MVN random walk proposal
    Σ = NeweyWest(chain[:,1:3])
    tuning = 1.0
    mcmcreps = 3
    for j = 1:mcmcreps
        P = (cholesky(Σ)).U
        proposal = θ -> proposal2(θ,tuning*P)
        θinit = vec(mean(chain[:,1:3],dims=1))
        if j == mcreps
            ChainLength = 10000
        end
        chain = mcmc(θinit, ChainLength, 0, prior, obj, proposal, false)
        if j < mcreps
            accept = mean(chain[:,end])
            if accept > 0.35
                tuning *= 1.5
            elseif accept < 0.25
                tuning *= 0.5
            end
            θinit = vec(mean(chain[:,1:3],dims=1))
            Σ = 0.5*Σ + 0.5*cov(chain[:,1:3])
        end
    end
    q1 = quantile(chain[:,1],[0.005, 0.025, 0.05, 0.95, 0.975, 0.995])
    q2 = quantile(chain[:,2],[0.005, 0.025, 0.05, 0.95, 0.975, 0.995])
    q3 = quantile(chain[:,3],[0.005, 0.025, 0.05, 0.95, 0.975, 0.995])
    bayes_inci = [
            θtrue[1] >= q1[1] && θtrue[1] <= q1[6],
            θtrue[2] >= q2[1] && θtrue[2] <= q2[6],
            θtrue[3] >= q3[1] && θtrue[3] <= q3[6],
            θtrue[1] >= q1[2] && θtrue[1] <= q1[5],
            θtrue[2] >= q2[2] && θtrue[2] <= q2[5],
            θtrue[3] >= q3[2] && θtrue[3] <= q3[5],
            θtrue[1] >= q1[3] && θtrue[1] <= q1[4],
            θtrue[2] >= q2[3] && θtrue[2] <= q2[4],
            θtrue[3] >= q3[3] && θtrue[3] <= q3[4]
           ]
    extremum_inci = [
            θtrue[1] >= θhat[1] - 2.576*se[1] && θtrue[1] <= θhat[1] + 2.576*se[1],
            θtrue[2] >= θhat[2] - 2.576*se[2] && θtrue[2] <= θhat[2] + 2.576*se[2],
            θtrue[3] >= θhat[3] - 2.576*se[3] && θtrue[3] <= θhat[3] + 2.576*se[3],
            θtrue[1] >= θhat[1] - 1.96 *se[1] && θtrue[1] <= θhat[1] + 1.96* se[1],
            θtrue[2] >= θhat[2] - 1.96* se[2] && θtrue[2] <= θhat[2] + 1.96* se[2],
            θtrue[3] >= θhat[3] - 1.96* se[3] && θtrue[3] <= θhat[3] + 1.96* se[3],
            θtrue[1] >= θhat[1] - 1.645*se[1] && θtrue[1] <= θhat[1] + 1.645*se[1],
            θtrue[2] >= θhat[2] - 1.645*se[2] && θtrue[2] <= θhat[2] + 1.645*se[2],
            θtrue[3] >= θhat[3] - 1.645*se[3] && θtrue[3] <= θhat[3] + 1.645*se[3]
            ]
   results[rep,:] = vcat(converged, θhat, extremum_inci, bayes_inci)
end
return results
end
results = main()
extremum_inci = mean(results[:,5:13],dims=1)
bayes_inci = mean(results[:,14:22],dims=1)
println("99% CI, extremum")
prettyprint(extremum_inci[1:3]', ["ϕ"; "ρ"; "σ"])
println("99% CI, Bayes")
prettyprint(bayes_inci[1:3]', ["ϕ"; "ρ"; "σ"])
println("95% CI, extremum")
prettyprint(extremum_inci[4:6]', ["ϕ"; "ρ"; "σ"])
println("95% CI, Bayes")
prettyprint(bayes_inci[4:6]', ["ϕ"; "ρ"; "σ"])
println("90% CI, extremum")
prettyprint(extremum_inci[7:9]', ["ϕ"; "ρ"; "σ"])
println("90% CI, Bayes")
prettyprint(bayes_inci[7:9]', ["ϕ"; "ρ"; "σ"])

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
module MatchingMarkets

include("util.jl")
include("deferred_acceptance.jl")
include("matching_tools.jl")

export deferred_acceptance, SProposing, CProposing
export random_prefs, ReturnCaps

end # module


text:
module NonSymFEMDriver

using Test
using Gridap
using GridapPardiso

model = CartesianDiscreteModel(
  domain=(0.0,1.0,0.0,1.0), partition=(2,2))

const T = VectorValue{2,Float64}

order = 2
fespace1 = FESpace(
  reffe = :QLagrangian,
  conformity = :H1,
  valuetype = T,
  model = model,
  order = order,
  diritags = "boundary")

fespace2 = FESpace(
  reffe = :PLagrangian,
  conformity = :L2,
  valuetype = Float64,
  model = model,
  order = order-1)

fespace2 = ConstrainedFESpace(fespace2,[1,])

V = TestFESpace(fespace1)
Q = TestFESpace(fespace2)
Y = [V, Q]

U = TrialFESpace(fespace1, zero(T))
P = TrialFESpace(fespace2)
X = [U, P]

trian = Triangulation(model)
quad = CellQuadrature(trian,degree=4)

function a_Ω(y,x)
  v, q = y
  u, p = x
  inner(∇(v),∇(u)) - inner(div(v),p) + inner(q,div(u))
end

function b_Ω(y)
  v, q = y
  inner(v,(x)->2*x)
end

t_Ω = AffineFETerm(a_Ω,b_Ω,trian,quad)

op = LinearFEOperator(Y,X,t_Ω)

ls = PardisoSolver()
# This one works: ls = LUSolver()
solver = LinearFESolver(ls)

xh = solve(solver,op)

@show num_free_dofs(U)
@show num_free_dofs(P)

x = free_dofs(xh)
A = op.mat
b = op.vec

@show all(abs.(A - A') .< 1e-13)

r = A*x - b

tol = 1e-10
@test maximum(abs.(r))/maximum(abs.(b)) < tol

end # module


text:
function printarray(A, fpath)
    open(fpath, "w") do file
        println(file, string(eltype(A))*"[")
        for ii = 1:size(A, 1)
            println(file, join(A[ii, :], " ")*";")
        end
        println(file, "]")
    end
end

function data_glass(material::Symbol)
    if material == :SiO2
        Float64[
            2000.0 0.99993 1.503e-5;
            1952.0 0.99992 1.64e-5;
            1905.0 0.99992 1.78e-5;
            1860.0 0.99991 1.936e-5;
            1815.0 0.99991 6.3e-6;
            1772.0 0.99991 6.9e-6;
            1730.0 0.9999 7.56e-6;
            1688.0 0.9999 8.3e-6;
            1648.0 0.99989 9.08e-6;
            1609.0 0.99989 9.941e-6;
            1570.0 0.99988 1.09e-5;
            1533.0 0.99987 1.19e-5;
            1496.0 0.99987 1.31e-5;
            1460.0 0.99986 1.43e-5;
            1426.0 0.99985 1.56e-5;
            1392.0 0.99985 1.71e-5;
            1358.0 0.99984 1.87e-5;
            1326.0 0.99983 2.04e-5;
            1294.0 0.99982 2.24e-5;
            1263.0 0.99981 2.45e-5;
            1233.0 0.99981 2.67e-5;
            1204.0 0.9998 2.92e-5;
            1175.0 0.99979 3.19e-5;
            1147.0 0.99978 3.48e-5;
            1119.0 0.99976 3.81e-5;
            1093.0 0.99975 4.16e-5;
            1067.0 0.99974 4.54e-5;
            1041.0 0.99973 4.97e-5;
            1016.0 0.99971 5.42e-5;
            992.0 0.9997 5.92e-5;
            968.0 0.99968 6.47e-5;
            945.0 0.99967 7.06e-5;
            923.0 0.99965 7.69e-5;
            901.0 0.99964 8.36e-5;
            879.0 0.99962 9.12e-5;
            858.0 0.9996 9.928e-5;
            838.0 0.99958 0.000108;
            818.0 0.99956 0.000118;
            798.0 0.99954 0.000129;
            779.0 0.99951 0.000141;
            760.0 0.99949 0.000154;
            742.0 0.99946 0.000167;
            725.0 0.99944 0.000181;
            707.0 0.99941 0.000198;
            690.0 0.99938 0.000215;
            674.0 0.99935 0.000233;
            658.0 0.99932 0.000254;
            642.0 0.99928 0.000277;
            627.0 0.99925 0.000301;
            612.0 0.99921 0.000327;
            597.0 0.99917 0.000356;
            583.0 0.99913 0.000386;
            569.0 0.99909 0.00042;
            555.0 0.99904 0.000458;
            542.0 0.99899 0.000497;
            529.0 0.99894 0.000149;
            516.0 0.99889 0.000162;
            504.0 0.99884 0.000176;
            492.0 0.99878 0.000191;
            480.0 0.99872 0.000208;
            469.0 0.99866 0.000226;
            458.0 0.99859 0.000245;
            447.0 0.99852 0.000267;
            436.0 0.99845 0.00029;
            426.0 0.99837 0.000315;
            415.0 0.99828 0.000344;
            406.0 0.99821 0.000371;
            396.0 0.99812 0.000404;
            386.0 0.99802 0.00044;
            377.0 0.99792 0.000477;
            368.0 0.99782 0.000518;
            359.0 0.99771 0.000563;
            351.0 0.9976 0.000607;
            342.0 0.99747 0.000662;
            334.0 0.99735 0.000717;
            326.0 0.99722 0.000777;
            318.0 0.99708 0.000844;
            311.0 0.99695 0.000909;
            303.0 0.99678 0.000991;
            296.0 0.99663 0.00107;
            289.0 0.99646 0.00116;
            282.0 0.99629 0.00125;
            275.0 0.99609 0.00136;
            269.0 0.99592 0.00146;
            262.0 0.9957 0.00159;
            256.0 0.99549 0.00171;
            250.0 0.99527 0.00185;
            244.0 0.99504 0.002;
            238.0 0.99478 0.00216;
            232.0 0.99451 0.00234;
            227.0 0.99427 0.0025;
            221.0 0.99395 0.00272;
            216.0 0.99367 0.00292;
            211.0 0.99336 0.00314;
            206.0 0.99304 0.00337;
            201.0 0.99269 0.00363;
            196.0 0.99231 0.0039;
            191.0 0.9919 0.00421;
            187.0 0.99155 0.00447;
            182.0 0.99108 0.00482;
            178.0 0.99068 0.00511;
            174.0 0.99024 0.00544;
            170.0 0.9898 0.0058;
            166.0 0.9893 0.00619;
            162.0 0.9887 0.00662;
            158.0 0.9882 0.00699;
            154.0 0.9875 0.0073;
            151.2 0.9871 0.0073;
            147.6 0.9865 0.0076;
            144.2 0.986 0.0081;
            140.9 0.9855 0.0089;
            137.8 0.9853 0.0099;
            134.8 0.9854 0.0107;
            131.9 0.9858 0.0114;
            129.2 0.9868 0.0111;
            126.5 0.9872 0.0111;
            124.0 0.9874 0.0097;
            121.6 0.9865 0.0085;
            119.2 0.9851 0.0083;
            117.0 0.9841 0.0096;
            114.8 0.9848 0.0099;
            112.7 0.9844 0.009;
            110.7 0.9822 0.0093;
            108.8 0.9828 0.0159;
            106.9 0.9867 0.0106;
            105.1 0.9858 0.0075;
            103.3 0.9839 0.0065;
            101.6 0.9823 0.0068;
            99.99 0.9813 0.007;
            98.4 0.9803 0.0073;
            96.86 0.9794 0.0076;
            95.37 0.9789 0.0076;
            93.93 0.9778 0.0083;
            92.53 0.977 0.0087;
            91.17 0.9761 0.0093;
            89.84 0.9747 0.0103;
            82.66 0.9634 0.0199;
            77.49 0.9608 0.0226;
            72.93 0.9562 0.0244;
            68.88 0.9509 0.0281;
            65.26 0.9458 0.0335;
            61.99 0.9416 0.0374;
            59.04 0.9386 0.0437;
            56.36 0.9328 0.0462;
            53.91 0.9271 0.052;
            51.66 0.9222 0.0578;
            49.59 0.9164 0.065;
            47.69 0.9105 0.0726;
            45.92 0.9207 0.068;
            44.28 0.9175 0.075;
            42.75 0.9137 0.082;
            41.33 0.913 0.09;
            40.0 0.907 0.092;
            38.75 0.901 0.094;
            37.57 0.895 0.098;
            36.47 0.888 0.107;
            35.42 0.882 0.113;
            34.44 0.877 0.12;
            33.51 0.87 0.128;
            32.63 0.866 0.137;
            31.79 0.858 0.144;
            31.0 0.851 0.156;
            30.24 0.845 0.169;
            29.52 0.839 0.18;
            28.83 0.833 0.19;
            28.18 0.827 0.202;
            27.55 0.822 0.218;
            26.95 0.817 0.233;
            26.38 0.813 0.25;
            25.83 0.808 0.27;
            25.3 0.804 0.282;
            25.0 0.733 0.325;
            24.8 0.803 0.3;
            24.31 0.804 0.322;
            24.0 0.753 0.375;
            23.84 0.806 0.343;
            23.39 0.811 0.366;
            23.0 0.774 0.434;
            22.96 0.817 0.385;
            22.54 0.822 0.408;
            22.14 0.829 0.43;
            22.0 0.797 0.48;
            21.75 0.833 0.45;
            21.38 0.843 0.47;
            21.01 0.851 0.482;
            21.0 0.827 0.53;
            20.66 0.862 0.497;
            20.0 0.859 0.585;
            19.5 0.879 0.613;
            19.0 0.902 0.645;
            18.5 0.927 0.677;
            18.0 0.957 0.712;
            17.75 0.975 0.731;
            17.5 0.999 0.75;
            17.25 1.03 0.763;
            17.0 1.072 0.768;
            16.75 1.124 0.765;
            16.5 1.137 0.755;
            16.25 1.156 0.737;
            16.0 1.172 0.717;
            15.75 1.178 0.703;
            15.5 1.172 0.696;
            15.25 1.167 0.699;
            15.0 1.168 0.711;
            14.75 1.175 0.739;
            14.5 1.195 0.771;
            14.25 1.225 0.799;
            14.0 1.265 0.808;
            13.75 1.32 0.795;
            13.5 1.363 0.775;
            13.25 1.371 0.755;
            13.0 1.368 0.747;
            12.75 1.372 0.766;
            12.5 1.383 0.793;
            12.25 1.41 0.824;
            12.0 1.475 0.861;
            11.8 1.554 0.874;
            11.6 1.635 0.859;
            11.4 1.716 0.81;
            11.2 1.766 0.718;
            11.0 1.739 0.569;
            10.9 1.687 0.565;
            10.8 1.587 0.618;
            10.7 1.513 0.725;
            10.6 1.492 0.914;
            10.5 1.567 1.11;
            10.45 1.645 1.136;
            10.4 1.772 1.13;
            10.35 1.919 1.045;
            10.3 2.048 0.925;
            10.25 2.152 0.81;
            10.2 2.24 0.715;
            10.1 2.332 0.46;
            10.0 2.33 0.323;
            9.9 2.292 0.236;
            9.8 2.243 0.168;
            9.7 2.19 0.119;
            9.6 2.14 0.077;
            9.5 2.092 0.0561;
            9.4 2.047 0.043;
            9.3 2.006 0.0339;
            9.2 1.969 0.0271;
            9.1 1.935 0.0228;
            9.0 1.904 0.0189;
            8.9 1.876 0.0156;
            8.8 1.85 0.0132;
            8.7 1.825 0.0109;
            8.6 1.803 0.00838;
            8.5 1.783 0.00557;
            8.4 1.764 0.00317;
            8.3 1.747 0.0014;
            8.2 1.73 0.000463;
            8.1 1.716 0.000122;
            8.0 1.702 3.2e-5;
            7.8 1.676 4.7e-6;
            7.6 1.653 0.0;
            7.4 1.633 0.0;
            7.2 1.616 0.0;
            7.0 1.6 0.0;
            6.75 1.582 0.0;
            6.5 1.567 0.0;
            6.25 1.554 0.0;
            6.0 1.543 0.0;
            5.7976 1.53429 0.0;
            5.7819 1.53371 0.0;
            5.468 1.52276 0.0;
            5.3858 1.52009 0.0;
            5.2131 1.51474 0.0;
            5.1674 1.51338 0.0;
            4.9939 1.50841 0.0;
            4.6751 1.50004 0.0;
            4.594 1.49805 0.0;
            4.504 1.49592 0.0;
            4.4226 1.49404 0.0;
            4.2848 1.49099 0.0;
            4.1784 1.48873 0.0;
            4.1034 1.48719 0.0;
            3.7542 1.48053 0.0;
            3.7105 1.47976 0.0;
            3.6427 1.47858 0.0;
            3.577 1.47746 0.0;
            3.434 1.47512 0.0;
            3.3967 1.47453 0.0;
            3.064 1.46961 0.0;
            2.8448 1.46669 0.0;
            2.6503 1.46429 0.0;
            2.5504 1.46313 0.0;
            2.4379 1.46187 0.0;
            2.2705 1.46008 0.0;
            2.1489 1.45885 0.0;
            2.1411 1.45877 0.0;
            2.1102 1.45847 0.0;
            2.1041 1.45841 0.0;
            1.9257 1.45671 0.0;
            1.8892 1.45637 0.0;
            1.8566 1.45608 0.0;
            1.7549 1.45515 0.0;
            1.455 1.45248 0.0;
            1.3863 1.45185 0.0;
            1.2228 1.45025 0.0;
            1.1449 1.44941 0.0;
            1.0985 1.44888 0.0;
            0.91018 1.44621 0.0;
            0.88874 1.44584 0.0;
            0.84372 1.44497 0.0;
            0.81061 1.44427 0.0;
            0.74663 1.44267 0.0;
            0.7376 1.44241 0.0;
            0.73225 1.44226 0.0;
            0.72543 1.44205 0.0;
            0.68384 1.44069 0.0;
            0.62934 1.43851 0.0;
            0.60243 1.43722 0.0;
            0.57598 1.43576 0.0;
            0.53317 1.43292 0.0;
            0.50868 1.43095 0.0;
            0.38221 1.41314 0.0;
            0.37953 1.41253 0.0;
            0.37542 1.41155 0.0;
            0.3623 1.40819 0.0;
            0.35354 1.40568 0.0;
            0.34863 1.40418 0.0;
            0.33449 1.39936 0.0;
            0.3224 1.395 0.0;
            0.2976 1.383 0.0;
            0.2728 1.365 0.0;
            0.248 1.342 0.0;
            0.2232 1.306 0.0;
            0.2108 1.278 0.0;
            0.1984 1.239 0.0;
            0.1922 1.212 0.0;
            0.186 1.175 0.0;
            0.1835 1.158 0.0;
            0.181 1.135 0.0;
            0.1785 1.107 0.0;
            0.1761 1.084 0.0;
            0.1736 1.053 0.0;
            0.1711 1.014 0.0;
            0.1686 0.9702 0.0;
            0.1674 0.9488 0.0;
            0.1661 0.9175 0.0;
            0.1649 0.8897 0.0;
            0.1637 0.86 0.0;
            0.1624 0.8213 0.0;
            0.1612 0.7719 0.0;
            0.1599 0.7037 0.0;
            0.1587 0.6232 0.0;
            0.1575 0.5456 0.0;
            0.1562 0.4677 0.0;
            0.155 0.4113 0.0;
            0.1537 0.3931 0.0;
            0.1525 0.402 0.0;
            0.1513 0.4329 0.0;
            0.15 0.453 0.0;
            0.1488 0.46 0.0;
            0.1475 0.473 0.0;
            0.1463 0.4746 0.0;
            0.1451 0.4656 0.0;
            0.1438 0.4563 0.0;
            0.1426 0.4309 0.0;
            0.1413 0.3915 0.0;
            0.1401 0.3563 0.0;
            0.1389 0.3705 0.0;
            0.1376 0.5846 0.0;
            0.1364 1.043 0.0;
            0.1351 1.616 0.0;
            0.1333 2.25 0.0;
            0.1302 2.76 0.0;
            0.1271 2.839 0.0;
            0.124 2.694 0.0;
            0.1209 2.448 0.0;
            0.1178 2.224 0.0;
            0.1147 2.038 0.0;
            0.1116 1.869 0.0;
            0.1091 1.784 0.0;
            0.1066 1.69 0.0;
            0.1054 1.652 0.0;
            0.1041 1.619 0.0;
            0.1029 1.615 0.0;
            0.1017 1.658 0.0;
            0.1004 1.701 0.0;
            0.09919 1.753 0.0;
            0.09795 1.789 0.0;
            0.09671 1.811 0.0;
            0.09547 1.81 0.0;
            0.09423 1.779 0.0;
            0.09299 1.756 0.0;
            0.08989 1.698 0.0;
            0.08679 1.643 0.0;
            0.08369 1.598 0.0;
            0.08059 1.555 0.0;
            0.07749 1.502 0.0;
            0.07439 1.45 0.0;
            0.07191 1.401 0.0;
            0.06943 1.337 0.0;
            0.06695 1.235 0.0;
            0.06534 1.161 0.0;
            0.06447 1.05 0.0;
            0.06323 0.8857 0.0;
            0.06199 0.6616 0.0;
            ]
    else
        error("No lookup table for glass $material")
    end
end

"Lookup tables for complex refractive indices of metals."
function data_metal(material::Symbol)
    # Below: 0.127: W. S. M. Werner, K. Glantschnig, C. Ambrosch-Draxl.
    # Optical constants and inelastic electron-scattering data for 17 elemental metals,
    # J. Phys Chem Ref. Data 38, 1013-1092 (2009)
    # Above 0.206: S. Babar and J. H. Weaver.
    # Optical constants of Cu, Ag, and Au revisited, Appl. Opt. 54, 477-481 (2015)
    # Below 0.206: K. Stahrenberg, Th. Herrmann, K. Wilmers, N. Esser, W. Richter, and M. J. G. Lee.
    # Optical properties of copper and silver in the energy range 2.5-9.0 eV,
    # Phys Rev. B 64, 115111 (2001) (Numerical data kindly provided by Prof. Dr. Norbert Esser)
    # above 12.4: H.-J. Hagemann, W. Gudat, and C. Kunz.
    # Optical constants from the far infrared to the x-ray region: Mg, Al, Cu, Ag, Au, Bi, C, and Al2O3,
    # J. Opt. Soc. Am. 65, 742-744 (1975)
    dat = (Ag = Float64[
        0.00236 0.99667 0.00774;
        0.00316 0.99770 0.00530;
        0.00447 0.99198 0.00237;
        0.00676 0.9813 0.00291;
        0.01140 0.9172 0.00818;
        0.01355 0.860 0.0423;
        0.01714 0.863 0.293;
        0.02430 0.928 0.310;
        0.02563 0.927 0.346;
        0.03038 0.895 0.302;
        0.03270 0.894 0.304;
        0.03568 0.872 0.351;
        0.03793 0.896 0.427;
        0.04059 0.953 0.451;
        0.04478 0.889 0.425;
        0.04607 0.848 0.446;
        0.04895 0.833 0.507;
        0.05391 1.04 0.746;
        0.05843 1.09 0.658;
        0.06163 1.06 0.564;
        0.06714 1.03 0.586;
        0.07185 0.982 0.596;
        0.07359 1.13 0.751;
        0.07437 1.17 0.721;
        0.08345 1.28 0.750;
        0.08780 1.37 0.688;
        0.09198 1.32 0.657;
        0.09321 1.27 0.625;
        0.09888 1.32 0.586;
        0.10258 1.23 0.571;
        0.10482 1.27 0.577;
        0.10667 1.30 0.565;
        0.10857 1.30 0.550;
        0.11350 1.33 0.518;
        0.12007 1.27 0.490;
        0.12157 1.21 0.515;
        0.12782 1.10047 0.72025;
        0.12808 1.10131 0.72092;
        0.12835 1.10177 0.72135;
        0.12861 1.10169 0.72134;
        0.12888 1.10100 0.72080;
        0.12915 1.09967 0.71975;
        0.12942 1.09762 0.71822;
        0.12969 1.09481 0.71619;
        0.12996 1.09145 0.71374;
        0.13024 1.08796 0.71109;
        0.13051 1.08474 0.70849;
        0.13078 1.08194 0.70609;
        0.13106 1.07956 0.70392;
        0.13134 1.07761 0.70203;
        0.13162 1.07600 0.70043;
        0.13190 1.07450 0.69911;
        0.13218 1.07305 0.69803;
        0.13246 1.07165 0.69717;
        0.13275 1.07019 0.69646;
        0.13303 1.06863 0.69586;
        0.13332 1.06700 0.69539;
        0.13360 1.06519 0.69501;
        0.13389 1.06311 0.69468;
        0.13418 1.06069 0.69431;
        0.13447 1.05788 0.69384;
        0.13477 1.05470 0.69326;
        0.13506 1.05129 0.69264;
        0.13535 1.04786 0.69208;
        0.13565 1.04452 0.69169;
        0.13595 1.04136 0.69159;
        0.13625 1.03835 0.69186;
        0.13655 1.03546 0.69254;
        0.13685 1.03258 0.69357;
        0.13715 1.02959 0.69486;
        0.13745 1.02649 0.69631;
        0.13776 1.02333 0.69786;
        0.13807 1.02008 0.69935;
        0.13838 1.01667 0.70070;
        0.13868 1.01312 0.70195;
        0.13900 1.00959 0.70326;
        0.13931 1.00630 0.70479;
        0.13962 1.00327 0.70659;
        0.13994 1.00043 0.70860;
        0.14025 0.99765 0.71068;
        0.14057 0.99478 0.71266;
        0.14089 0.99170 0.71446;
        0.14121 0.98833 0.71602;
        0.14153 0.98477 0.71742;
        0.14186 0.98106 0.71871;
        0.14218 0.97724 0.71992;
        0.14251 0.97339 0.72109;
        0.14284 0.96943 0.72222;
        0.14317 0.96532 0.72333;
        0.14350 0.96113 0.72452;
        0.14383 0.95706 0.72591;
        0.14417 0.95333 0.72764;
        0.14450 0.95000 0.72972;
        0.14484 0.94698 0.73206;
        0.14518 0.94414 0.73458;
        0.14552 0.94137 0.73720;
        0.14586 0.93855 0.73982;
        0.14621 0.93563 0.74239;
        0.14655 0.93257 0.74486;
        0.14690 0.92931 0.74718;
        0.14725 0.92578 0.74931;
        0.14760 0.92200 0.75125;
        0.14795 0.91810 0.75314;
        0.14831 0.91416 0.75509;
        0.14866 0.91017 0.75713;
        0.14902 0.90608 0.75929;
        0.14938 0.90197 0.76166;
        0.14974 0.89796 0.76437;
        0.15010 0.89411 0.76744;
        0.15047 0.89043 0.77090;
        0.15083 0.88696 0.77474;
        0.15120 0.88370 0.77897;
        0.15157 0.88061 0.78354;
        0.15194 0.87762 0.78837;
        0.15231 0.87469 0.79342;
        0.15269 0.87180 0.79864;
        0.15307 0.86889 0.80393;
        0.15345 0.86587 0.80921;
        0.15383 0.86267 0.81443;
        0.15421 0.85933 0.81961;
        0.15459 0.85590 0.82477;
        0.15498 0.85246 0.82997;
        0.15537 0.84911 0.83528;
        0.15576 0.84597 0.84077;
        0.15615 0.84313 0.84652;
        0.15655 0.84073 0.85258;
        0.15694 0.83879 0.85897;
        0.15734 0.83725 0.86563;
        0.15774 0.83597 0.87247;
        0.15814 0.83480 0.87934;
        0.15855 0.83363 0.88619;
        0.15895 0.83246 0.89296;
        0.15936 0.83130 0.89969;
        0.15977 0.83023 0.90640;
        0.16019 0.82930 0.91312;
        0.16060 0.82853 0.91987;
        0.16102 0.82797 0.92667;
        0.16144 0.82762 0.93353;
        0.16186 0.82750 0.94043;
        0.16228 0.82762 0.94737;
        0.16271 0.82798 0.95434;
        0.16314 0.82857 0.96136;
        0.16357 0.82933 0.96839;
        0.16400 0.83021 0.97545;
        0.16444 0.83118 0.98251;
        0.16487 0.83223 0.98958;
        0.16531 0.83088 0.99959;
        0.16575 0.83127 1.01158;
        0.16620 0.82727 1.02560;
        0.16665 0.82808 1.03261;
        0.16709 0.82886 1.03820;
        0.16755 0.82962 1.04289;
        0.16800 0.83035 1.04727;
        0.16846 0.83108 1.05195;
        0.16892 0.83186 1.05693;
        0.16938 0.83273 1.06215;
        0.16984 0.83368 1.06808;
        0.17031 0.83476 1.07462;
        0.17078 0.83595 1.08125;
        0.17125 0.83723 1.08787;
        0.17172 0.83860 1.09421;
        0.17220 0.84003 1.10000;
        0.17268 0.84149 1.10533;
        0.17316 0.84298 1.11071;
        0.17365 0.84453 1.11637;
        0.17413 0.84615 1.12232;
        0.17463 0.84786 1.12868;
        0.17512 0.84966 1.13554;
        0.17561 0.85155 1.14216;
        0.17611 0.85353 1.14846;
        0.17662 0.85559 1.15476;
        0.17712 0.85774 1.16112;
        0.17763 0.85997 1.16755;
        0.17814 0.86228 1.17387;
        0.17865 0.86465 1.18014;
        0.17917 0.86708 1.18656;
        0.17969 0.86952 1.19291;
        0.18021 0.87199 1.19907;
        0.18073 0.87449 1.20500;
        0.18126 0.87706 1.21088;
        0.18179 0.87969 1.21643;
        0.18233 0.88237 1.22161;
        0.18287 0.88512 1.22669;
        0.18341 0.88797 1.23201;
        0.18395 0.89094 1.23792;
        0.18450 0.89402 1.24420;
        0.18505 0.89720 1.25057;
        0.18560 0.90045 1.25677;
        0.18616 0.90374 1.26250;
        0.18672 0.90705 1.26788;
        0.18729 0.91045 1.27341;
        0.18785 0.91394 1.27929;
        0.18843 0.91750 1.28535;
        0.18900 0.92111 1.29122;
        0.18958 0.92475 1.29654;
        0.19016 0.92845 1.30169;
        0.19074 0.93225 1.30693;
        0.19133 0.93614 1.31227;
        0.19193 0.94009 1.31770;
        0.19252 0.94411 1.32307;
        0.19312 0.94815 1.32826;
        0.19373 0.95224 1.33342;
        0.19433 0.95637 1.33861;
        0.19494 0.96053 1.34372;
        0.19556 0.96469 1.34858;
        0.19618 0.96886 1.35307;
        0.19680 0.97303 1.35722;
        0.19743 0.97727 1.36115;
        0.19806 0.98163 1.36507;
        0.19869 0.98618 1.36935;
        0.19933 0.99098 1.37427;
        0.19997 0.99599 1.37936;
        0.20062 1.00116 1.38449;
        0.20127 1.00650 1.38965;
        0.20193 1.01200 1.39469;
        0.20259 1.01760 1.39932;
        0.20325 1.02330 1.40363;
        0.20392 1.02910 1.40767;
        0.2066 1.079 1.247;
        0.2101 1.101 1.258;
        0.2138 1.121 1.267;
        0.2175 1.140 1.272;
        0.2214 1.157 1.280;
        0.2254 1.169 1.283;
        0.2296 1.181 1.287;
        0.2339 1.190 1.294;
        0.2384 1.198 1.306;
        0.2431 1.215 1.317;
        0.2480 1.233 1.326;
        0.2530 1.254 1.328;
        0.2583 1.270 1.327;
        0.2638 1.294 1.329;
        0.2695 1.321 1.310;
        0.2755 1.348 1.291;
        0.2818 1.380 1.247;
        0.2883 1.402 1.181;
        0.2952 1.427 1.079;
        0.3024 1.410 0.894;
        0.3100 1.265 0.656;
        0.3179 0.859 0.466;
        0.3263 0.307 0.651;
        0.3351 0.187 1.042;
        0.3444 0.142 1.269;
        0.3542 0.106 1.467;
        0.3647 0.076 1.645;
        0.3757 0.061 1.815;
        0.3875 0.050 1.981;
        0.3999 0.054 2.138;
        0.4133 0.050 2.292;
        0.4275 0.051 2.448;
        0.4428 0.052 2.604;
        0.4592 0.052 2.765;
        0.4769 0.053 2.930;
        0.4959 0.052 3.105;
        0.5166 0.052 3.288;
        0.5391 0.05168 3.483;
        0.5636 0.05009 3.694;
        0.5904 0.04977 3.918;
        0.6199 0.04803 4.164;
        0.6525 0.04964 4.432;
        0.6888 0.05080 4.725;
        0.7293 0.05151 5.048;
        0.7749 0.05269 5.409;
        0.8266 0.05504 5.814;
        0.8856 0.05896 6.276;
        0.9537 0.06542 6.802;
        1.033 0.07286 7.412;
        1.127 0.08427 8.128;
        1.240 0.09799 8.981;
        1.305 0.1056 9.472;
        1.378 0.1143 10.018;
        1.459 0.1256 10.63;
        1.550 0.1388 11.31;
        1.653 0.1543 12.08;
        1.771 0.1735 12.97;
        1.907 0.1988 13.98;
        2.066 0.2315 15.16;
        2.254 0.2724 16.56;
        2.480 0.3256 18.23;
        2.755 0.3970 20.26;
        3.100 0.4960 22.80;
        3.542 0.6387 26.05;
        4.133 0.8566 30.37;
        4.959 1.228 36.37;
        6.199 1.851 45.26;
        8.266 3.227 59.73;
        12.4000 5.079 86.53;
        2.480E+01 3.670E+01 1.73E+02;
        4.959E+01 1.180E+02 3.06E+02;
        1.240E+02 3.090E+02 5.06E+02;
        2.480E+02 5.310E+02 6.89E+02
        ],
    # Below 0.14: H.-J. Hagemann, W. Gudat, and C. Kunz.
    # Optical constants from the far infrared to the x-ray region: Mg, Al, Cu, Ag, Au, Bi, C, and Al2O3,
    # J. Opt. Soc. Am. 65, 742-744 (1975)
    # K. M. McPeak, S. V. Jayanti, S. J. P. Kress, S. Meyer, S. Iotti, A. Rossinelli, and D. J. Norris.
    # Plasmonic films can easily be better: Rules and recipes, ACS Photonics 2, 326-333 (2015)
    # Above 2.0: M. A. Ordal, R. J. Bell, R. W. Alexander, L. A. Newquist, M. R. Querry.
    # Optical properties of Al, Fe, Ti, Ta, W, and Mo at submillimeter wavelengths,
    # Appl. Opt. 27, 1203-1209 (1988)
    Al = Float64[
        1.033E-05 9.990E-01 8.24E-12;
        1.240E-05 1.001E+00 9.88E-12;
        1.550E-05 1.000E+00 2.47E-11;
        2.066E-05 1.001E+00 8.09E-11;
        2.480E-05 1.001E+00 1.50E-10;
        3.100E-05 1.002E+00 3.77E-10;
        4.133E-05 1.002E+00 1.06E-09;
        6.199E-05 1.001E+00 4.68E-09;
        8.266E-05 1.000E+00 1.46E-08;
        1.127E-04 1.002E+00 5.14E-08;
        1.378E-04 1.000E+00 1.12E-07;
        1.771E-04 1.003E+00 2.97E-07;
        2.480E-04 1.002E+00 1.05E-06;
        4.133E-04 1.001E+00 7.51E-06;
        6.199E-04 1.000E+00 3.24E-05;
        7.293E-04 1.000E+00 6.13E-05;
        7.514E-04 1.000E+00 4.62E-05;
        7.749E-04 1.000E+00 2.32E-05;
        7.999E-04 1.000E+00 6.43E-06;
        8.266E-04 1.000E+00 7.14E-06;
        1.240E-03 1.000E+00 3.10E-05;
        1.550E-03 1.000E+00 7.52E-05;
        2.066E-03 9.990E-01 2.27E-04;
        2.480E-03 9.990E-01 4.66E-04;
        3.100E-03 9.980E-01 1.01E-03;
        3.542E-03 9.960E-01 1.54E-03;
        4.133E-03 9.960E-01 2.57E-03;
        4.428E-03 9.960E-01 3.31E-03;
        4.769E-03 9.920E-01 4.26E-03;
        5.166E-03 9.930E-01 5.46E-03;
        5.636E-03 9.930E-01 6.96E-03;
        6.199E-03 9.910E-01 8.76E-03;
        6.525E-03 9.890E-01 9.87E-03;
        6.888E-03 9.890E-01 1.11E-02;
        7.293E-03 9.880E-01 1.37E-02;
        7.749E-03 9.880E-01 1.60E-02;
        8.266E-03 9.890E-01 1.73E-02;
        8.551E-03 9.880E-01 1.81E-02;
        8.856E-03 9.880E-01 1.91E-02;
        8.984E-03 9.880E-01 1.97E-02;
        9.116E-03 9.870E-01 2.04E-02;
        9.253E-03 9.870E-01 2.13E-02;
        9.393E-03 9.860E-01 2.26E-02;
        9.537E-03 9.860E-01 2.44E-02;
        9.686E-03 9.870E-01 2.61E-02;
        9.840E-03 9.880E-01 2.77E-02;
        9.999E-03 9.890E-01 2.88E-02;
        1.016E-02 9.910E-01 2.92E-02;
        1.033E-02 9.920E-01 2.94E-02;
        1.042E-02 9.920E-01 2.97E-02;
        1.051E-02 9.920E-01 3.00E-02;
        1.060E-02 9.930E-01 3.03E-02;
        1.069E-02 9.930E-01 3.01E-02;
        1.078E-02 9.940E-01 3.04E-02;
        1.088E-02 9.940E-01 3.06E-02;
        1.097E-02 9.950E-01 3.07E-02;
        1.107E-02 9.950E-01 3.09E-02;
        1.117E-02 9.950E-01 3.09E-02;
        1.127E-02 9.960E-01 3.07E-02;
        1.137E-02 9.970E-01 3.05E-02;
        1.148E-02 9.970E-01 2.99E-02;
        1.159E-02 9.970E-01 2.93E-02;
        1.170E-02 9.960E-01 2.86E-02;
        1.181E-02 9.960E-01 2.81E-02;
        1.192E-02 9.940E-01 2.81E-02;
        1.204E-02 9.920E-01 2.86E-02;
        1.216E-02 9.920E-01 2.99E-02;
        1.228E-02 9.890E-01 3.24E-02;
        1.240E-02 9.890E-01 3.63E-02;
        1.252E-02 9.910E-01 3.96E-02;
        1.265E-02 9.940E-01 4.27E-02;
        1.278E-02 9.980E-01 4.42E-02;
        1.292E-02 1.001E+00 4.42E-02;
        1.305E-02 1.005E+00 4.39E-02;
        1.319E-02 1.007E+00 4.26E-02;
        1.333E-02 1.009E+00 4.11E-02;
        1.348E-02 1.011E+00 3.98E-02;
        1.362E-02 1.012E+00 3.86E-02;
        1.378E-02 1.013E+00 3.66E-02;
        1.393E-02 1.014E+00 3.51E-02;
        1.409E-02 1.014E+00 3.35E-02;
        1.425E-02 1.014E+00 3.21E-02;
        1.442E-02 1.014E+00 3.15E-02;
        1.459E-02 1.014E+00 3.12E-02;
        1.476E-02 1.016E+00 3.06E-02;
        1.494E-02 1.017E+00 2.89E-02;
        1.512E-02 1.017E+00 2.69E-02;
        1.531E-02 1.016E+00 2.50E-02;
        1.550E-02 1.016E+00 2.43E-02;
        1.569E-02 1.015E+00 2.38E-02;
        1.590E-02 1.014E+00 2.38E-02;
        1.610E-02 1.015E+00 2.45E-02;
        1.631E-02 1.016E+00 2.40E-02;
        1.653E-02 1.018E+00 2.40E-02;
        1.675E-02 1.021E+00 2.36E-02;
        1.687E-02 1.025E+00 2.40E-02;
        1.691E-02 1.028E+00 2.38E-02;
        1.696E-02 1.030E+00 1.96E-02;
        1.701E-02 1.036E+00 1.74E-02;
        1.705E-02 1.031E+00 6.06E-03;
        1.710E-02 1.026E+00 4.79E-03;
        1.722E-02 1.021E+00 4.61E-03;
        1.746E-02 1.015E+00 4.50E-03;
        1.771E-02 1.012E+00 4.45E-03;
        1.797E-02 1.009E+00 4.40E-03;
        1.851E-02 1.005E+00 4.33E-03;
        1.907E-02 1.001E+00 4.27E-03;
        2.066E-02 9.920E-01 4.51E-03;
        2.480E-02 9.740E-01 5.84E-03;
        3.100E-02 9.440E-01 8.64E-03;
        4.133E-02 8.820E-01 1.25E-02;
        6.199E-02 6.730E-01 2.67E-02;
        6.888E-02 5.640E-01 3.55E-02;
        7.749E-02 3.460E-01 6.48E-02;
        7.999E-02 2.520E-01 9.97E-02;
        8.266E-02 1.520E-01 1.87E-01;
        8.551E-02 1.070E-01 3.01E-01;
        8.856E-02 8.400E-02 4.04E-01;
        9.537E-02 6.500E-02 5.93E-01;
        1.033E-01 6.100E-02 7.67E-01;
        1.127E-01 6.100E-02 9.46E-01;
        1.378E-01 7.300E-02 1.35E+00;
        0.15 0.095390828 1.283666394;
        0.155 0.095510386 1.337393822;
        0.16 0.09903925 1.402928641;
        0.165 0.098692838 1.46616516;
        0.17 0.100850207 1.532569987;
        0.175 0.1069563 1.596539899;
        0.18 0.099715746 1.657661977;
        0.185 0.108316112 1.734368006;
        0.19 0.106567769 1.79116071;
        0.195 0.111513266 1.853411775;
        0.2 0.110803374 1.908606137;
        0.205 0.111587326 1.969936987;
        0.21 0.11365555 2.028057589;
        0.215 0.115928445 2.091850713;
        0.22 0.116173424 2.151998203;
        0.225 0.119771906 2.213039835;
        0.23 0.124315543 2.274896559;
        0.235 0.129276519 2.336773934;
        0.24 0.133002743 2.395514502;
        0.245 0.139688588 2.457358928;
        0.25 0.141162655 2.515219055;
        0.255 0.148765766 2.573777623;
        0.26 0.150722638 2.633214255;
        0.265 0.161465056 2.693431077;
        0.27 0.164610587 2.750035753;
        0.275 0.172365826 2.81099804;
        0.28 0.178635303 2.868973784;
        0.285 0.18587018 2.923276881;
        0.29 0.188953314 2.98248275;
        0.295 0.197315218 3.042340102;
        0.3 0.204991638 3.100858199;
        0.305 0.210097266 3.157347124;
        0.31 0.218816718 3.214790935;
        0.315 0.224454799 3.267473023;
        0.32 0.237666454 3.323775766;
        0.325 0.24464928 3.37932734;
        0.33 0.251892832 3.436779015;
        0.335 0.259391824 3.493749292;
        0.34 0.267481852 3.550147397;
        0.345 0.275201252 3.607178361;
        0.35 0.28349792 3.663518946;
        0.355 0.291774119 3.719829584;
        0.36 0.300125667 3.776251503;
        0.365 0.308578012 3.831608068;
        0.37 0.317597538 3.889198914;
        0.375 0.32692637 3.945995546;
        0.38 0.335956002 4.002607455;
        0.385 0.345714203 4.058265431;
        0.39 0.354901676 4.114238107;
        0.395 0.364968364 4.169676475;
        0.4 0.375150842 4.226433266;
        0.405 0.385211589 4.281283449;
        0.41 0.396086448 4.336805792;
        0.415 0.40706511 4.391435073;
        0.42 0.417647849 4.447407079;
        0.425 0.429543735 4.503499508;
        0.43 0.440996226 4.559109307;
        0.435 0.452837879 4.613682059;
        0.44 0.464232752 4.669101846;
        0.445 0.477070026 4.724044433;
        0.45 0.489220122 4.778319404;
        0.455 0.501228231 4.832828338;
        0.46 0.514817248 4.887948117;
        0.465 0.528042125 4.943260916;
        0.47 0.53987657 4.997345967;
        0.475 0.554932886 5.051331558;
        0.48 0.568005038 5.105940631;
        0.485 0.582771419 5.159226931;
        0.49 0.596705366 5.212874187;
        0.495 0.610784373 5.265933683;
        0.5 0.625686295 5.320477736;
        0.505 0.640306464 5.374848125;
        0.51 0.655709839 5.428163169;
        0.515 0.672565753 5.481545831;
        0.52 0.688336416 5.535955836;
        0.525 0.704045108 5.58799423;
        0.53 0.720793584 5.641767754;
        0.535 0.737603948 5.693773756;
        0.54 0.75446839 5.746496546;
        0.545 0.772330366 5.799008035;
        0.55 0.789405353 5.851936501;
        0.555 0.808351698 5.905288517;
        0.56 0.829205097 5.958236408;
        0.565 0.848630853 6.009090973;
        0.57 0.867376853 6.060226283;
        0.575 0.887661988 6.112687118;
        0.58 0.908569739 6.165184423;
        0.585 0.928308533 6.215605902;
        0.59 0.948955518 6.267417058;
        0.595 0.9714896 6.317103137;
        0.6 0.992465612 6.368986418;
        0.605 1.016073317 6.417351867;
        0.61 1.038145667 6.46681888;
        0.615 1.062059906 6.517052343;
        0.62 1.088160063 6.566589093;
        0.625 1.112663572 6.61475776;
        0.63 1.136328574 6.66312559;
        0.635 1.165731637 6.710806975;
        0.64 1.190265203 6.759417967;
        0.645 1.218505245 6.807110623;
        0.65 1.246364405 6.852329839;
        0.655 1.275302761 6.898361917;
        0.66 1.304382818 6.942363757;
        0.665 1.333854457 6.987643931;
        0.67 1.365410391 7.031951759;
        0.675 1.395892303 7.071104251;
        0.68 1.426024482 7.116287362;
        0.685 1.457391234 7.157652493;
        0.69 1.493230683 7.197481356;
        0.695 1.5274935 7.23603734;
        0.7 1.559751729 7.27391404;
        0.705 1.596346402 7.311383032;
        0.71 1.631621525 7.34575992;
        0.715 1.669698976 7.380372258;
        0.72 1.706780893 7.414249862;
        0.725 1.745260386 7.446739811;
        0.73 1.786398844 7.475374294;
        0.735 1.827772224 7.502849546;
        0.74 1.872189153 7.528352557;
        0.745 1.916802427 7.55216072;
        0.75 1.958355454 7.571403838;
        0.755 2.005349601 7.585211073;
        0.76 2.054224115 7.599462923;
        0.765 2.09984451 7.605086914;
        0.77 2.144528834 7.61216158;
        0.775 2.189915668 7.610567038;
        0.78 2.231500036 7.603853787;
        0.785 2.275479945 7.589387919;
        0.79 2.314643646 7.572223927;
        0.795 2.350258044 7.547894599;
        0.8 2.373653298 7.522581337;
        0.805 2.399438967 7.486541749;
        0.81 2.409296598 7.449477169;
        0.815 2.414395973 7.410756231;
        0.82 2.410049326 7.369351244;
        0.825 2.399462094 7.325600203;
        0.83 2.375277499 7.288411569;
        0.835 2.342102357 7.253098886;
        0.84 2.301774531 7.225729444;
        0.845 2.257656113 7.204374444;
        0.85 2.204898553 7.188085811;
        0.855 2.148111401 7.177700369;
        0.86 2.092718783 7.181328004;
        0.865 2.03127715 7.184642431;
        0.87 1.973568962 7.198887028;
        0.875 1.915492386 7.218248477;
        0.88 1.854996083 7.243411521;
        0.885 1.802564727 7.274079983;
        0.89 1.750120449 7.309782595;
        0.895 1.694671727 7.350096072;
        0.9 1.64715608 7.393045594;
        0.905 1.602399525 7.439468008;
        0.91 1.557188172 7.48978204;
        0.915 1.516837858 7.542665118;
        0.92 1.483545494 7.592621895;
        0.925 1.447230511 7.647710517;
        0.93 1.415213787 7.70342976;
        0.935 1.385433286 7.756057335;
        0.94 1.355589911 7.812465936;
        0.945 1.325884954 7.872100794;
        0.95 1.302164463 7.928721983;
        0.955 1.280288795 7.985307721;
        0.96 1.25976713 8.04593147;
        0.965 1.23821822 8.1051949;
        0.97 1.218188621 8.163760378;
        0.975 1.201067407 8.22268587;
        0.98 1.182294465 8.279240316;
        0.985 1.167631515 8.338247019;
        0.99 1.153696234 8.399524689;
        0.995 1.140905235 8.457757667;
        1 1.126639087 8.511598888;
        1.005 1.114680776 8.571960874;
        1.01 1.102964941 8.629920722;
        1.015 1.09358083 8.687515124;
        1.02 1.08727974 8.74843856;
        1.025 1.081445553 8.809967337;
        1.03 1.071397488 8.867857271;
        1.035 1.061537336 8.920081923;
        1.04 1.055494597 8.978255334;
        1.045 1.051806876 9.036830069;
        1.05 1.043454926 9.083765639;
        1.055 1.036482441 9.148731926;
        1.06 1.033133204 9.204040646;
        1.065 1.02747793 9.264432573;
        1.07 1.025591906 9.318521275;
        1.075 1.024418445 9.373122913;
        1.08 1.020244607 9.42862127;
        1.085 1.015304156 9.479917855;
        1.09 1.014378517 9.5368739;
        1.095 1.011890304 9.593431189;
        1.1 1.009905005 9.649917519;
        1.105 1.008451539 9.704482186;
        1.11 1.005344855 9.750882949;
        1.115 1.00382405 9.806001363;
        1.12 1.004294808 9.860152537;
        1.125 1.003882164 9.917498646;
        1.13 1.003972742 9.968398129;
        1.135 1.00500573 10.02382345;
        1.14 1.002828586 10.07500199;
        1.145 1.003941179 10.12708734;
        1.15 1.005317282 10.18335224;
        1.155 1.003754384 10.23706613;
        1.16 1.008818757 10.28628773;
        1.165 1.007932489 10.34146643;
        1.17 1.008357257 10.39215013;
        1.175 1.012688334 10.44200832;
        1.18 1.015213183 10.49336256;
        1.185 1.014759181 10.5517646;
        1.19 1.017683803 10.60545522;
        1.195 1.01813905 10.65817091;
        1.2 1.021736807 10.70579031;
        1.205 1.022859293 10.75807795;
        1.21 1.024257787 10.80782605;
        1.215 1.026450838 10.85926339;
        1.22 1.029738722 10.9110202;
        1.225 1.035047848 10.96181208;
        1.23 1.038667581 11.01266681;
        1.235 1.039978895 11.06849385;
        1.24 1.043389606 11.11569439;
        1.245 1.046657092 11.16984741;
        1.25 1.049913391 11.21616325;
        1.255 1.051297547 11.26655345;
        1.26 1.056495057 11.31619114;
        1.265 1.057188404 11.36079431;
        1.27 1.061492702 11.41081797;
        1.275 1.066383574 11.47288865;
        1.28 1.070500901 11.51937378;
        1.285 1.073080696 11.56250415;
        1.29 1.078091053 11.61775711;
        1.295 1.081263285 11.66829595;
        1.3 1.08667656 11.71967858;
        1.305 1.090905752 11.7737027;
        1.31 1.094227716 11.82020872;
        1.315 1.097106966 11.86438552;
        1.32 1.104427092 11.91930196;
        1.325 1.108138114 11.96737106;
        1.33 1.109282178 12.01542787;
        1.335 1.112492808 12.06389822;
        1.34 1.123987022 12.1162431;
        1.345 1.129233973 12.16532652;
        1.35 1.132981867 12.21775953;
        1.355 1.135566133 12.26700903;
        1.43 1.20897633 12.99198306;
        1.435 1.215774284 13.0402495;
        1.44 1.222359347 13.08984959;
        1.445 1.225637477 13.13438949;
        1.45 1.232348412 13.19350911;
        1.455 1.235053048 13.23424936;
        1.46 1.23989324 13.2825952;
        1.465 1.245674136 13.32973008;
        1.47 1.25590184 13.37524914;
        1.475 1.252421713 13.42119816;
        1.48 1.265153644 13.46516297;
        1.485 1.264792932 13.51960433;
        1.49 1.275023692 13.56523813;
        1.495 1.277693673 13.61515704;
        1.5 1.285699173 13.66027168;
        1.505 1.290149697 13.70964938;
        1.51 1.296809948 13.75760575;
        1.515 1.304059962 13.81378197;
        1.52 1.313924676 13.8577198;
        1.525 1.31627889 13.90548777;
        1.53 1.325289614 13.94978109;
        1.535 1.332565366 13.9955611;
        1.54 1.338103685 14.0463704;
        1.545 1.342317392 14.08906725;
        1.55 1.347399401 14.13278052;
        1.555 1.354287918 14.18605286;
        1.56 1.36162657 14.2300396;
        1.565 1.367953829 14.27204602;
        1.57 1.374881699 14.33214219;
        1.575 1.374949427 14.37325593;
        1.58 1.385662006 14.41873986;
        1.585 1.390425967 14.47146449;
        1.59 1.400170689 14.51965833;
        1.595 1.39786136 14.56653756;
        1.6 1.400984694 14.60865354;
        1.605 1.402928097 14.65292487;
        1.61 1.421414675 14.70647543;
        1.615 1.424742484 14.75635087;
        1.62 1.436002763 14.82100212;
        1.625 1.439358574 14.84862799;
        1.63 1.44609075 14.8956765;
        1.635 1.447904247 14.94712102;
        1.64 1.45643876 14.99170484;
        1.645 1.458873135 15.02392461;
        1.65 1.469917493 15.08179888;
        1.655 1.478177828 15.1223348;
        1.66 1.474977808 15.17582154;
        1.665 1.477843053 15.23233469;
        1.67 1.462292742 15.2773787;
        1.675 1.425778483 15.34740513;
        1.68 1.427335087 15.38951869;
        1.685 1.454908413 15.4724516;
        1.69 1.514176047 15.49621661;
        1.695 1.555668449 15.53351552;
        1.7 1.584018511 15.55632073;
        2.00 2.1962737 20.969371;
        2.11 2.4038118 22.062867;
        2.22 2.6179915 23.257800;
        2.35 2.8842412 24.591905;
        2.50 3.1815630 26.073820;
        2.67 3.5693590 27.735962;
        2.86 4.0097886 29.575977;
        3.08 4.5104289 31.653186;
        3.33 5.0466806 34.058286;
        3.64 5.7888000 36.878956;
        4.00 6.5579910 40.184716;
        4.44 7.6669569 44.210650;
        5.00 8.8866481 49.131664;
        5.71 10.625810 55.544522;
        6.67 13.385937 63.967532;
        8.00 17.708575 75.320562;
        10.0 25.832564 90.720430;
        11.1 30.167884 98.172167;
        12.5 35.492290 107.26425;
        14.3 42.226361 117.86671;
        16.7 50.744867 131.07397;
        20.0 61.871514 148.17314;
        22.2 68.262241 159.23227;
        25.0 76.681189 174.09538;
        28.6 88.918176 192.23917;
        33.3 108.95823 213.37170;
        40.0 134.05166 235.70537;
        44.4 149.32871 248.55977;
        50.0 167.07810 263.04366;
        57.1 188.04844 279.77106;
        66.7 213.38844 299.72607;
        80.0 245.00980 324.59867;
        100 286.08777 357.47647;
        125 329.97183 393.77470;
        154 374.73526 431.44791;
        200 436.98909 485.19932
        ]
    )
    dat[material]
end


text:
logscale(params::AbstractVector{T}) where T <: Real = /(extrema(params)...) < 2e-2 && minimum(params) > 0

logscale(params) = false

@recipe function plot(ho::Hyperoptimizer)
    N = length(ho.params)
    layout --> N
    for i = 1:N
        params = getindex.(ho.history, i)
        if eltype(params) <: Real
            perm = sortperm(params)
            yguide --> "Function value"
            seriestype --> :scatter
            @series begin
                xguide --> ho.params[i]
                subplot --> i
                label --> "Sampled points"
                legend --> false
                xscale --> (logscale(params) ? :log10 : :identity)
                yscale --> (logscale(ho.results) ? :log10 : :identity)
                params[perm], ho.results[perm]
            end
        else
            params = Symbol.(params)
            uniqueparams = sort(unique(params))
            paramidxs = map(x -> findfirst(==(x), uniqueparams), params)
            yguide --> "Function value"
            seriestype --> :scatter
            @series begin
                xguide --> ho.params[i]
                subplot --> i
                label --> "Sampled points"
                legend --> false
                xticks --> (1:length(uniqueparams), uniqueparams)
                xscale --> :identity
                yscale --> (logscale(ho.results) ? :log10 : :identity)
                paramidxs, ho.results
            end
        end
    end
end


text:
using MAT
using Plots
using Statistics
using StatsBase
include("../src/Graph_Louvain.jl")
num = 9
numtimes = 5
Mstart = 4
Mend = 4
mat = matread("Output/Amazon_$(num)_($numtimes)_$(Mstart)_$(Mend).mat")

lin_ari = mat["lin_ari"]
wce_ari = mat["wce_ari"]
sce_ari = mat["sce_ari"]
lin_t = mat["lin_t"]
SC = mat["SC"]
WC = mat["WC"]
LinC = mat["LinC"]
numpoints = length(lin_t)

##
M = matread("../data/Amazon_$num.mat")
names = M["names"]
labels = M["labels"]

# Nine smallest Amazon categories, see 5-core dataset here:
# https://nijianmo.github.io/amazon/index.html#subsets
# for names of categories. This order lists them in increasing size
uLabels = [1 3 2 12 18 25 17 24 15]
uLs = uLabels[1:num]

H = M["H"]
order = vec(sum(H,dims = 2))
m,n = size(H)
vol = sum(H)

## Confirm the hypergraph is one connecte component
A = H'*H
using MatrixNetworks
Acc, p = largest_component(A)

## Compute what percentage of nodes are completely contained in a cluster
mu = 0.75
High = 0
All = 0
for e = 1:m
    global mu, High, All
    edge = findall(x->x>0, H[e,:])
    clusters = labels[edge]
    mm = mode(clusters)
    prcnt = length(findall(x->x==mm,clusters))/length(clusters)
    if prcnt == 1
        All += 1
    elseif prcnt > mu
        High += 1
    end
end

@show All/m, High/m

## Plot
using Plots
using LaTeXStrings
lw = 3
msw = 0
s1 = 300
s2 = 250
ms = 5
lw = 2

lin = mat["lin_ari"]
wce = mat["wce_ari"]
Lams = range(1/(Mstart*vol),Mend/vol,length=numpoints)

## Begin by plotting ARI scores for detecting different category-clusters
xax = L"\lambda * vol(G)"
p = plot(vol*Lams, lin, color=:purple,grid = false,linewidth = lw,xaxis = xax,
legend = :bottom, yaxis= "ARI Scores", label = "HyperLam (Star) ",
markerstrokewidth = msw, markersize = ms, markershape = :circle)

plot!(vol*Lams, wce ,color=:green,linewidth= lw, label = "HyperLam (Clique) ",
markerstrokewidth = msw, markersize = ms, markershape = :circle, size = (s1,s2))

savefig(p,"Plots/Amazon_aris.pdf")

## Now track individual clusters
Fl = zeros(num,numpoints)
Rl = zeros(num,numpoints)
for i = 1:numpoints
    c = round.(Int64,LinC[:,i])
    F1s, PRs, REs =  ClusterTracker(c,labels, uLs)
    Fl[:,i] = F1s
    Rl[:,i] = REs
end
F = zeros(num,numpoints)
R = zeros(num,numpoints)
for i = 1:numpoints
    c = round.(Int64,WC[:,i])
    F1s, PRs, REs =  ClusterTracker(c,labels, uLs)
    F[:,i] = F1s
    R[:,i] = REs
end

## Plot additional figures

plot()
#"Software",
lbs = ["Pantry", "I&S"]
clus = [8 9]
colors = [:red, :black, :blue]
for i = 2:-1:1
    l = clus[i]
    plot!(Lams*vol, Fl[l,:], grid = false, linewidth = lw, markershape = :circle,
    label = lbs[i],legend = :topleft,markerstrokewidth = 0,  size = (s1,s2),
    markersize = ms, color = colors[i])
end
for i = 2:-1:1
    l = clus[i]
    plot!(Lams*vol, F[l,:], grid = false, linewidth = lw, markershape = :circle,
    label = "", legend = :bottomleft,linestyle = :dash, markerstrokewidth = 0,  size = (s1,s2),
    markersize = ms, color = colors[i],xaxis = xax, yaxis= "F1 Scores")
end
savefig("Plots/TwoClusters_F1_final.pdf")

##

plot()

lbs = ["Soft.", "Pantry", "I&S"]
clus = [6 8 9]
colors = [:red, :black, :blue]
for i = 1:3
    l = clus[i]
    plot!(Lams*vol, Fl[l,:], grid = false, linewidth = lw, markershape = :circle,
    label = lbs[i],legend = :topleft,markerstrokewidth = 0,  size = (s1,s2),
    markersize = ms, color = colors[i])
end
for i = 1:3
    l = clus[i]
    plot!(Lams*vol, F[l,:], grid = false, linewidth = lw, markershape = :circle,
    label = "", legend = :topleft,linestyle = :dash, markerstrokewidth = 0,  size = (s1,s2),
    markersize = ms, color = colors[i],xaxis = xax, yaxis= "F1 Scores")
end
savefig("Plots/ThreeClusters_F1_final.pdf")

max_stars_repo_name: nveldt/ParamCC

max_issues_repo_name: nveldt/ParamCC

max_forks_repo_name: nveldt/ParamCC


text:
import ..Properties: @property
import ..Vortex

@property begin
    signature = induce_acc(targ::Target, targvel::Target, src::Source, srcvel::Source)
    preallocator = allocate_acc
    stype = ComplexF64
end

function induce_acc(z::ComplexF64, ż::ComplexF64,
                    p::Vortex.Point, point_vel::ComplexF64)
    -circulation(p)*Points.cauchy_kernel((z - p.z)^2)*conj(ż - point_vel)
end

function induce_acc(z::ComplexF64, ż::ComplexF64,
                           b::Vortex.Blob, blob_vel::ComplexF64)
    induce_acc(z, ż,
               Vortex.Point(Elements.position(b), circulation(b)),
               blob_vel)
end

function surface_pressure_inst(p::Plate, ṗ, ambient_sys, z_new, t, Δt, lesp, tesp)
    @get p (L, C, ss, α, dchebt!)
    @get ṗ (ċ, α̇)

    # Get Ċ from movement of existing vortex blobs (without vortex shedding)
    enforce_no_flow_through!(p, ṗ, ambient_sys, t)

    srcvel = self_induce_velocity(ambient_sys, t)
    induce_velocity!(srcvel, ambient_sys, p, t)

    targvel = fill(ċ, length(p))
    Ċ = zero(targvel)

    induce_acc!(Ċ, p.zs, targvel, ambient_sys, srcvel)

    z₊, z₋ = z_new
    point₊ = Vortex.Point(z₊, 1.0)
    point₋ = Vortex.Point(z₋, 1.0)

    Γ₊, Γ₋, ∂C₊, ∂C₋ = vorticity_flux!(p, point₊, point₋, t, lesp, tesp)

    n̂ = exp(-im*α)
    rmul!(Ċ, n̂)
    dchebt! * Ċ

    @. Ċ += (∂C₊ + ∂C₋)/Δt - im*α̇*C

    Ȧ = MappedVector(imag, Ċ, 1)

    # Plate is moving at a fixed velocity and angle of attack
    Ḃ₀ = Ḃ₁ = 0.0

    Γ̇ = [Γ₋/Δt + _bound_circulation(Ȧ, Ḃ₀, Ḃ₁, L, -(Γ₊ + Γ₋)/Δt, s) for s in p.ss]

    strength(p) .* (Chebyshev.firstkind(real.(C), ss) .- tangent(ċ, α)) .+ Γ̇
end


text:
function play_turn(state::DynamicState, static_state::StaticState,
    decision::Tuple{UInt8, UInt8})
    next_state = state

    fm_pending = get_fast_moves_pending(state)
    active1, active2 = get_active(next_state)
    cmp = get_cmp(state)

    if !iszero(cmp)
        agent, o_agent = isodd(cmp) ? (0x01, 0x02) : (0x02, 0x01)
        next_state = evaluate_charged_move(next_state, static_state, cmp,
            decision[agent] - 0x06, 0x64, decision[o_agent] == 0x01)
        data = next_state.data - UInt32(fm_pending[agent]) *
            UInt32(agent == 0x01 ? 16 : 112)
        if !iszero(fm_pending[o_agent])
            data -= UInt32(fm_pending[o_agent] - 1) *
                UInt32(o_agent == 0x01 ? 16 : 112)
        end
        next_state = DynamicState(next_state[0x01], next_state[0x02], data)
    else
        if fm_pending[1] == 0x01 || fm_pending[2] == 0x01
            next_state = evaluate_fast_moves(next_state, static_state,
                (fm_pending[1] == 0x01 &&
                    get_hp(next_state[0x01][active1]) != 0x0000,
                fm_pending[2] == 0x01 &&
                    get_hp(next_state[0x02][active2]) != 0x0000))
        end

        next_state = step_timers(next_state,
            decision[1] == 0x03 ?
                static_state[0x01][active1].fastMove.cooldown : Int8(0),
            decision[2] == 0x03 ?
                static_state[0x02][active2].fastMove.cooldown : Int8(0))
        if decision[1] == 0x05 || decision[1] == 0x06
            next_state = evaluate_switch(next_state,
                static_state, 0x01, decision[1] - 0x04,
                iszero(get_hp(state[0x01][active1])) ?
                0x18 : 0x00)
        end
        if decision[2] == 0x05 || decision[2] == 0x06
            next_state = evaluate_switch(next_state,
                static_state, 0x02, decision[2] - 0x04,
                iszero(get_hp(state[0x02][active2])) ?
                0x18 : 0x00)
        end
        active1, active2 = get_active(next_state)

        if get_hp(next_state[0x01][active1]) != 0x0000 &&
            get_hp(next_state[0x02][active2]) != 0x0000
            if decision[1] == 0x04
                if decision[2] == 0x04
                    atk_cmp = Base.cmp(
                        static_state[0x01][active1].stats.attack,
                        static_state[0x02][active2].stats.attack
                    )
                    return DynamicState(next_state[0x01],
                        next_state[0x02], next_state.data +
                        (atk_cmp == 1 ? UInt32(2352) :
                        atk_cmp == -1 ? UInt32(3136) : UInt32(19600)))
                else
                    return DynamicState(next_state[0x01],
                        next_state[0x02], next_state.data + UInt32(784))
                end
            elseif decision[2] == 0x04
                return DynamicState(next_state[0x01],
                    next_state[0x02], next_state.data + UInt32(2*784))
            end
        end
    end

    return next_state
end

function get_chance_state_1(state::DynamicState, static_state::StaticState,
    chance::UInt8)
    if chance == 0x05
        return DynamicState(state[0x01], state[0x02], state.data - 0x4360)
    else
        active1, active2 = get_active(state)
        agent, o_agent = chance < 0x03 ? (0x01, 0x02) : (0x02, 0x01)
        move = isodd(chance) ?
            static_state[agent][active1].charged_move_1 :
            static_state[agent][active2].charged_move_2
        a_data = state[agent].data
        d_data = state[o_agent].data
        a_data, d_data = apply_buff(a_data, d_data, move)
        return update_fm_damage(DynamicState(
            DynamicTeam(state[0x01][0x01], state[0x01][0x02],
                state[0x01][0x03], state[0x01].switchCooldown,
                agent == 0x01 ? a_data : d_data),
            DynamicTeam(state[0x02][0x01], state[0x02][0x02],
                state[0x02][0x03], state[0x02].switchCooldown,
                agent == 0x02 ? a_data : d_data),
            state.data - UInt32(chance) * UInt32(3920)), static_state)
    end
end

function get_chance_state_2(state::DynamicState, static_state::StaticState,
    chance::UInt8)
    if chance == 0x05
        return DynamicState(state[0x01], state[0x02], state.data - 0x4050)
    else
        return DynamicState(state[0x01], state[0x02],
            state.data - UInt32(chance) * UInt32(3920))
    end
end

function resolve_chance(state::DynamicState, static_state::StaticState)
    chance = get_chance(state)
    if chance == 0x00
        return state
    elseif chance == 0x05
        return rand() < 0.5 ?
            # subtract chance, add cmp
            get_chance_state_1(state, static_state, chance) :
            get_chance_state_2(state, static_state, chance)
    else
        active1, active2 = get_active(state)
        agent = chance < 0x03 ? 0x01 : 0x02
        buff_chance = isodd(chance) ?
            static_state[agent][active1].charged_move_1.buffChance :
            static_state[agent][active2].charged_move_2.buffChance
        return rand(Int8(0):Int8(99)) < buff_chance ?
            get_chance_state_1(state, static_state, chance) :
            get_chance_state_2(state, static_state, chance)
    end
end

"""
    play_battle(state, static_state)

Play through one battle, starting from the inputted state with random,
equally weighted decisions.
"""
function play_battle(state::DynamicState, static_state::StaticState;
  allow_nothing::Bool = false, allow_overfarming::Bool = false)
    active1, active2 = get_active(state)
    while true
        state = resolve_chance(state, static_state)
        d1, d2 = get_possible_decisions(state, static_state,
            allow_nothing = allow_nothing,
            allow_overfarming = allow_overfarming)
        (iszero(d1) || iszero(d2)) &&
            return battle_score(state, static_state)
        state = play_turn(
            state, static_state, select_random_decision(d1, d2))
    end
end

"""
    battle_scores(state, static_state, N)

Play through N battles, starting from the inputted state with random,
equally weighted decisions.
"""
function battle_scores(starting_state::DynamicState,
  static_state::StaticState, N::Int64)
    return map(x -> play_battle(starting_state, static_state), 1:N)
end


text:
export proxL21!, normL21

"""
    proxL21!(x::Vector{T},λ::Float64; sparseTrafo::Trafo=nothing, slices::Int64=1, kargs...)

group-soft-thresholding for l1/l2-regularization.

# Arguments:
* `x::Array{T}`                 - Vector to apply proximal map to
* `λ::Float64`                  - regularization paramter
* `sparseTrafo::Trafo=nothing`  - sparsifying transform to apply
* `slices::Int64=1`             - number of elements per group
"""
function proxL21!(x::Vector{T},λ::Float64; sparseTrafo::Trafo=nothing, slices::Int64=1, kargs...) where T
  if sparseTrafo != nothing
    z = sparseTrafo*x
  else
    z = x
  end
  if λ != 0
    proxL21!(z, λ, slices)
  end
  if sparseTrafo != nothing
    x[:] = adjoint(sparseTrafo)*z
  else
    x[:] = z
  end
  return x
end

function proxL21!(x::Vector{T}, λ::Float64, slices::Int64) where T
  sliceLength = div(length(x),slices)
  groupNorm = [norm(x[i:sliceLength:end]) for i=1:sliceLength]
  x[:] = [ x[i]*max( (groupNorm[mod1(i,sliceLength)]-λ)/groupNorm[mod1(i,sliceLength)],0 ) for i=1:length(x)]
end

"""
    normL21(x::Vector{T}, λ::Float64; sparseTrafo::Trafo=nothing, slices::Int64=1, kargs...) where T

return the value of the L21-regularization term.
Arguments are the same as in `proxL21!`
"""
function normL21(x::Vector{T}, λ::Float64; sparseTrafo::Trafo=nothing, slices::Int64=1, kargs...) where T
  if sparseTrafo != nothing
    z = sparseTrafo*x
  else
    z = x
  end
  sliceLength = div(length(z),slices)
  groupNorm = [norm(z[i:sliceLength:end]) for i=1:sliceLength]
  return λ*norm(groupNorm,1)
end


text:
#CORRELATION AND SIMULATION TOOLS FOR MC HAMMER
#by Eric Torkia, April 2019

# mch_simtools contains simulation tools and functions to obtain outputs from simulation arrays generated by MCHammer

"""
    cormat(ArrayName, RankOrder=1)

Cormat calculates a symetric correlation matrix using both PPMC and Rank Order. Rank Order is default because this is what it used in the Iman-Conover method for correlating of simulated variables.

**RankOrder = 1** calculates the Spearman rank order correlation used in MCHammer (this argument is optional and defaults to Spearman)

**RankOrder = 0** calculates the Pearson Product Moment Correlation
"""
function cormat(ArrayName, RankOrder=1)
cor_mat = []
M_Size = size(ArrayName,2)

#Converting to a data frame is necessary (and simpler) to perform rankings
if typeof(ArrayName) == Array{Float64,2}
      ArrayName = DataFrame(ArrayName)
end

if RankOrder == 1
      for i=1:M_Size
      cor_vector = []
            for i2=1:M_Size
             cor_i = cor(tiedrank(ArrayName[!, i]),tiedrank(ArrayName[!, i2]));
             push!(cor_vector,cor_i)
             #print(cov_i)
            end
       push!(cor_mat,cor_vector)
      end
else
      for i=1:M_Size
      cor_vector = []
            for i2=1:M_Size
             cor_i = cor(ArrayName[!,i], ArrayName[!,i2]);
             push!(cor_vector,cor_i)
             #print(cov_i)
      end
      push!(cor_mat,cor_vector)
      end
end

return hcat(cor_mat...)
end

"""
    covmat(ArrayName)

Calculates the covariance matrix.
"""
function covmat(ArrayName)
M_Size = size(ArrayName,2)
cov_mat = []

if typeof(ArrayName) == Array{Float64,2}
      ArrayName = DataFrame(ArrayName)

end
      for i=1:M_Size
            for i2=1:M_Size
             cov_i = cov(ArrayName[!,i],ArrayName[!,i2]);
             push!(cov_mat,cov_i)
             #print(cov_i)
            end
      end
return reshape(cov_mat, M_Size, M_Size)
end

#Correlation using Iman-Conover
"""
    corvar(ArrayName, n_trials, correl_matrix)

The corvar function correlates simulation inputs unsing the Iman Conover Method. Your array must contain >2 simulated inputs. **Remember to hcat() your inputs into tables reflecting your input correlation matrices.**

**n_trials**: is the number of trials in the simulation. This must be consistent.

**correl_matrix**: must be defined as a Square Positive Definite correlation matrix. This can be calculated from histroical data using `cormat()` function.
"""
function corvar(ArrayName, n_trials, correl_matrix)

if typeof(ArrayName) == Array{Float64,2}
      ArrayName = DataFrame(ArrayName)

end

#Define how many columns of ISNs are required
array_dims = size(ArrayName,2)

#calc cholesky transform to transfer correlation to ISNs
P = cholesky(correl_matrix)

#Normal() Returns Standard Normals (ISNs)
R = rand(Normal(),n_trials,array_dims)
ISN_Matrix = R*P.U
ISN_Matrix_DF = DataFrame(ISN_Matrix)

#apply ranks to create independant correlation rankings matrix
ISN_Ranked = []
for i = 1:array_dims
      temp_ranks = ordinalrank(ISN_Matrix_DF[!,i])
      push!(ISN_Ranked, temp_ranks)
end
ISN_Ranked_DF = DataFrame(ISN_Ranked)
#Results_Ranked = []
#println(ISN_Ranked_DF)

#Reindex the array of samples using the ISN_Ranks. Sort(Array)[OrderingVector]
final_array=[]
for i = 1:array_dims
      sorted_array = sort(ArrayName[!,i])[ISN_Ranked[i]]
      push!(final_array, sorted_array)
end
Final_DF=DataFrame(final_array)
end

# GetCertainty
"""
    GetCertainty(ArrayName, x, AboveBelow=0)

This function returns the percentage of trials Above (1) or Below(0) a target value of x.
"""
function GetCertainty(ArrayName, x, AboveBelow=0)
      if AboveBelow == 1
            certainty = count(i ->(i>=x), ArrayName)/size(ArrayName,1)
      else
            certainty = count(i ->(i<=x), ArrayName)/size(ArrayName,1)

      end
      return certainty
end

"""
    cmd(x)

Shell /Dos Command wrapper to run batch and shell commands in script. This is used to process SQL from the command line or perform system level operation in a script using a command prompt.
"""
cmd(x) = run(`cmd /C $x`)

"""
    function  truncate_digit(num, digits=2)
Truncation algorithim to remove decimals (ported by anonymous author from Maple) e.g.

      0.066 = 0.06
      0.063 = 0.06

"""
function  truncate_digit(num, digits=2)
    if num == 0.0 then
        return num
    else
        e = ceil(log10(abs(num)))
        scale = 10^(digits - e)
        return trunc(num * scale) / scale
    end
end

#fractiles function calculates percentiles at equal increments.
"""
    fractiles(ArrayName, Increment=0.1)

The fractiles function calculates percentiles at equal increments. The default optional argument for Increments is 0.1 for deciles but can be set to anything such as 0.05 for quintiles or 0.01 for percentiles.
"""
function fractiles(ArrayName, Increment=0.1)
names = map(i -> "P" * string(i), collect(0:Increment:1)*100)
P_Vals = collect(0:Increment:1)
#convert(Array{Float64}, ArrayName)
fractiles = quantile(collect(Float64, ArrayName), P_Vals)
return hcat(names, fractiles)
end

#convert(Array{Float64}, ArrayName)
# convert(Array{Int64}, collect(0:Increment:1)*100)

#Testing Script
#--------------------------
#Must setup as vectors using commas in the []

#Final_Ranks = []
# for i = 1:array_dims
#       temp_ranks = ordinalrank(Final_DF[i])
#       push!(Final_Ranks, temp_ranks)
# end
# Final_Ranks_DF = DataFrame(Final_Ranks)

#Useful Variables
#Phi = (sqrt(5)+1)/2
# n_trials = 1000
# correl_matrix = [1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0;0 0 0 1 0.75 -0.7; 0 0 0 0.75 1 -0.95; 0 0 0 -0.7 -0.95 1 ]
#
# #Test Data
# test_logNormal = [rand(LogNormal(0, 0.5),n_trials) rand(Normal(3,2),n_trials) rand(Gamma(1, 0.5),n_trials) rand(LogNormal(0, 0.5),n_trials) rand(Normal(3,2),n_trials) rand(Gamma(1, 0.5),n_trials)]


text:
# random functions that dispatch either to CURAND or GPUArrays' generic RNG

using Random

export rand_logn!, rand_poisson!

# native RNG

"""
    CUDA.RNG()

A random number generator using `rand()` in a device kernel.

See also: [CUDA.Philox2x32](@ref)
"""
mutable struct RNG <: AbstractRNG
    seed::UInt32
    counter::UInt32

    function RNG(seed::Integer)
        new(seed%UInt32, 0)
    end
end

RNG() = RNG(Random.rand(UInt32))

function Random.seed!(rng::RNG, seed::Integer)
    rng.seed = seed % UInt32
    rng.counter = 0
end

Random.seed!(rng::RNG) = seed!(rng, Random.rand(UInt32))

function Random.rand!(rng::RNG, A::AnyCuArray)
    function kernel(A::AbstractArray{T}, seed::UInt32, counter::UInt32) where {T}
        device_rng = Random.default_rng()

        # initialize the state
        @inbounds Random.seed!(device_rng, seed, counter)

        # grid-stride loop
        threadId = threadIdx().x
        offset = (blockIdx().x - 1) * blockDim().x
        while offset < length(A)
            i = threadId + offset
            if i <= length(A)
                @inbounds A[i] = Random.rand(device_rng, T)
            end

            offset += (blockDim().x - 1) * gridDim().x
        end

        return
    end

    kernel = @cuda launch=false name="rand!" kernel(A, rng.seed, rng.counter)
    config = launch_configuration(kernel.fun; max_threads=64)
    threads = max(32, min(config.threads, length(A)))
    blocks = min(config.blocks, cld(length(A), threads))
    kernel(A, rng.seed, rng.counter; threads=threads, blocks=blocks)

    new_counter = Int64(rng.counter) + length(A)
    overflow, remainder = fldmod(new_counter, typemax(UInt32))
    rng.seed += overflow     # XXX: is this OK?
    rng.counter = remainder

    A
end

# TODO: `randn!`; cannot reuse from Base or RandomNumbers, as those do scalar indexing

# generic functionality

function Random.rand!(rng::Union{RNG,CURAND.RNG,GPUArrays.RNG}, A::AbstractArray{T}) where {T}
    B = CuArray{T}(undef, size(A))
    Random.rand!(rng, B)
    copyto!(A, B)
end

function Random.rand(rng::Union{RNG,CURAND.RNG,GPUArrays.RNG}, T::Type)
    assertscalar("scalar rand")
    A = CuArray{T}(undef, 1)
    Random.rand!(rng, A)
    A[]
end

# RNG-less interface

# the interface is split in two levels:
# - functions that extend the Random standard library, and take an RNG as first argument,
#   will only ever dispatch to CURAND and as a result are limited in the types they support.
# - functions that take an array will dispatch to either CURAND or GPUArrays
# - non-exported functions are provided for constructing GPU arrays from only an eltype

const curand_rng = CURAND.default_rng
gpuarrays_rng() = GPUArrays.default_rng(CuArray)

function seed!(seed=Base.rand(UInt64))
    Random.seed!(curand_rng(), seed)
    Random.seed!(gpuarrays_rng(), seed)
end

# CURAND in-place
Random.rand!(A::CURAND.UniformArray) = Random.rand!(curand_rng(), A)
Random.randn!(A::CURAND.NormalArray; kwargs...) = Random.randn!(curand_rng(), A; kwargs...)
rand_logn!(A::CURAND.LognormalArray; kwargs...) = CURAND.rand_logn!(curand_rng(), A; kwargs...)
rand_poisson!(A::CURAND.PoissonArray; kwargs...) = CURAND.rand_poisson!(curand_rng(), A; kwargs...)

# CURAND out-of-place
rand(T::CURAND.UniformType, dims::Dims) = Random.rand(curand_rng(), T, dims)
randn(T::CURAND.NormalType, dims::Dims; kwargs...) = Random.randn(curand_rng(), T, dims; kwargs...)
rand_logn(T::CURAND.LognormalType, dims::Dims; kwargs...) = CURAND.rand_logn(curand_rng(), T, dims; kwargs...)
rand_poisson(T::CURAND.PoissonType, dims::Dims; kwargs...) = CURAND.rand_poisson(curand_rng(), T, dims; kwargs...)

# support all dimension specifications
rand(T::CURAND.UniformType, dim1::Integer, dims::Integer...) =
    Random.rand(curand_rng(), T, Dims((dim1, dims...)))
randn(T::CURAND.NormalType, dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn(curand_rng(), T, Dims((dim1, dims...)); kwargs...)
rand_logn(T::CURAND.LognormalType, dim1::Integer, dims::Integer...; kwargs...) =
    CURAND.rand_logn(curand_rng(), T, Dims((dim1, dims...)); kwargs...)
rand_poisson(T::CURAND.PoissonType, dim1::Integer, dims::Integer...; kwargs...) =
    CURAND.rand_poisson(curand_rng(), T, Dims((dim1, dims...)); kwargs...)

# GPUArrays in-place
Random.rand!(A::AnyCuArray) = Random.rand!(gpuarrays_rng(), A)
Random.randn!(A::AnyCuArray; kwargs...) =
    error("CUDA.jl does not support generating normally-distributed random numbers of type $(eltype(A))")
# FIXME: GPUArrays.jl has a randn! nowadays, but it doesn't work with e.g. Cuint
rand_logn!(A::AnyCuArray; kwargs...) =
    error("CUDA.jl does not support generating lognormally-distributed random numbers of type $(eltype(A))")
rand_poisson!(A::AnyCuArray; kwargs...) =
    error("CUDA.jl does not support generating Poisson-distributed random numbers of type $(eltype(A))")

# GPUArrays out-of-place
rand(T::Type, dims::Dims) = Random.rand!(CuArray{T}(undef, dims...))
randn(T::Type, dims::Dims; kwargs...) = Random.randn!(CuArray{T}(undef, dims...); kwargs...)
rand_logn(T::Type, dims::Dims; kwargs...) = rand_logn!(CuArray{T}(undef, dims...); kwargs...)
rand_poisson(T::Type, dims::Dims; kwargs...) = rand_poisson!(CuArray{T}(undef, dims...); kwargs...)

# support all dimension specifications
rand(T::Type, dim1::Integer, dims::Integer...) =
    Random.rand!(CuArray{T}(undef, dim1, dims...))
randn(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn!(CuArray{T}(undef, dim1, dims...); kwargs...)
rand_logn(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    rand_logn!(CuArray{T}(undef, dim1, dims...); kwargs...)
rand_poisson(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    rand_poisson!(CuArray{T}(undef, dim1, dims...); kwargs...)

# untyped out-of-place
rand(dim1::Integer, dims::Integer...) =
    Random.rand(curand_rng(), Dims((dim1, dims...)))
randn(dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn(curand_rng(), Dims((dim1, dims...)); kwargs...)
rand_logn(dim1::Integer, dims::Integer...; kwargs...) =
    CURAND.rand_logn(curand_rng(), Dims((dim1, dims...)); kwargs...)
rand_poisson(dim1::Integer, dims::Integer...; kwargs...) =
    CURAND.rand_poisson(curand_rng(), Dims((dim1, dims...)); kwargs...)


text:
## these are differently named than SymPy or missing or ...

Base.abs(x::SymbolicObject) = sympy.Abs(x)
Base.abs2(x::SymbolicObject) = x * conj(x)
Base.max(x::Sym, a) = sympy.Max(x, a)
Base.min(x::Sym, a) = sympy.Min(x, a)

Base.cbrt(x::Sym) = x^(1//3)
Base.ceil(x::Sym) = sympy.ceiling(x)

## Trig
Base.asech(z::Sym) = log(sqrt(1/z-1)*sqrt(1/z+1) + 1/z)
Base.acsch(z::Sym) = log(sqrt(1+1/z^2) + 1/z) ## http://mathworld.wolfram.com/InverseHyperbolicCosecant.html
Base.atan(y::Sym, x) = sympy.atan2(y,x)

Base.sinc(x::Sym) = iszero(x) ? one(x) : sin(PI*x)/(PI*x)
cosc(x::Sym) = diff(sinc(x))

Base.sincos(x::Sym) = (sin(x), cos(x))
Base.sinpi(x::Sym) = sympy.sin(x*PI)
Base.cospi(x::Sym) = sympy.cos(x*PI)
degree_variants = (:sind, :cosd, :tand, :cotd, :secd, :cscd,
                   :asind, :acosd, :atand, :acotd, :asecd, :acscd)

for  methvar in degree_variants
    meth = Symbol(String(methvar)[1:end-1])
    @eval begin
        (Base.$methvar)(ex::SymbolicObject) = ($meth)((PI/180)*ex)
    end
end

Base.rad2deg(x::Sym) = (x * 180) / PI
Base.deg2rad(x::Sym) = (x * PI) / 180

Base.hypot(x::Sym, y::Number) = hypot(promote(x,y)...)
Base.hypot(xs::Sym...) = sqrt(sum(abs(xᵢ)^2 for xᵢ ∈ xs))

## exponential
Base.log1p(x::Sym) = sympy.log(1 + x)
Base.log(x::Sym) = sympy.log(x)
Base.log(b::Number, x::Sym) = sympy.log(x, b)
Base.log2(x::SymbolicObject) = log(2,x)
Base.log10(x::SymbolicObject) = log(10,x)

## calculus.
## use a pair for limit x=>0
limit(x::SymbolicObject, xc::Pair, args...;kwargs...) = limit(x, xc[1], xc[2], args...;kwargs...)
## allow a function
limit(f::Function, x::Sym, c;kwargs...) = limit(Sym(f(x)), x, c; kwargs...)
function  limit(f::Function, c;kwargs...)
    @vars  x
    limit(f, x, c; kwargs...)
end

## This is  type piracy  and a bad   idea
function Base.diff(f::Function, n::Integer=1)
    @vars x
    sympy.diff(f(x), x, n)
end

## integrate(ex,a,b)
function integrate(ex::SymbolicObject, a::Number, b::Number)
    fs = free_symbols(ex)
    if length(fs) !== 1
        @warn "Need exactly on free symbol. Use `integrate(ex, (x, a, b))` instead"
        return
    end
    integrate(ex, (fs[1], a, b))
end
function integrate(f::Function, a::Number, b::Number)
    @vars x
    sympy.integrate(f(x), (x, a, b))
end
function integrate(f::Function)
    @syms  x
    sympy.integrate(f(x), x)
end

## Add interfaces for solve, nonlinsolve when vector of equations passed in

## An alternative to Eq(lhs, rhs) following Symbolics.jl
"""
    lhs ~ rhs

Specify an equation.

Alternative syntax to `Eq(lhs, rhs)` or `lhs ⩵ rhs` (`\\Equal[tab]`) following `Symbolics.jl`.
"""
Base.:~(lhs::Number, rhs::SymbolicObject) = Eq(lhs, rhs)
Base.:~(lhs::SymbolicObject, rhs::Number) = Eq(lhs, rhs)
Base.:~(lhs::SymbolicObject, rhs::SymbolicObject) = Eq(lhs, rhs)

"""
    solve

Use `solve` to solve algebraic equations.

Examples:

```julia
julia> using SymPy

julia> @syms x y a b c d
(x, y, a, b, c, d)

julia> solve(x^2 + 2x + 1, x) # [-1]
1-element Vector{Sym}:
 -1

julia> solve(x^2 + 2a*x + a^2, x) # [-a]
1-element Vector{Sym}:
 -a

julia> solve([a*x + b*y-3, c*x + b*y - 1], [x,y]) # Dict(y => (a - 3*c)/(b*(a - c)),x => 2/(a - c))
Dict{Any, Any} with 2 entries:
  y => (a - 3*c)/(a*b - b*c)
  x => 2/(a - c)

```

!!! note
    A very nice example using `solve` is a [blog](https://newptcai.github.io/euclidean-plane-geometry-with-julia.html) entry on [Napolean's theorem](https://en.wikipedia.org/wiki/Napoleon%27s_theorem) by Xing Shi Cai.
"""
solve() = ()
solve(V::Vector{T}, args...; kwargs...) where {T <: SymbolicObject} =
    sympy.solve(V, args...; kwargs...)

"""
    nonlinsolve

Note: if passing variables in use a tuple (e.g., `(x,y)`) and *not* a vector (e.g., `[x,y]`).
"""
nonlinsolve(V::AbstractArray{T,N}, args...; kwargs...) where {T <: SymbolicObject, N} =
    sympy.nonlinsolve(V, args...; kwargs...)

linsolve(V::AbstractArray{T,N}, args...; kwargs...) where {T <: SymbolicObject, N} =
    sympy.linsolve(V, args...; kwargs...)
linsolve(Ts::Tuple, args...; kwargs...) where {T <: SymbolicObject} =
    sympy.linsolve(Ts, args...; kwargs...)

## dsolve allowing initial condiation to be specified

"""
   dsolve(eqn, var, args..,; ics=nothing, kwargs...)

Call `sympy.dsolve`.

The initial conditions are specified with a dictionary.

Example:

```jldoctest dsolve
julia> using SymPy

julia> @syms α, x, f(), g()
(α, x, f, g)

julia> ∂ = Differential(x)
Differential(x)

julia> eqn = ∂(f(x)) ~ α * x
d
──(f(x)) = x⋅α
dx
```

```julia
julia> dsolve(eqn)
             2
            x ⋅α
f(x) = C₁ + ────
             2
```

```jldoctest dsolve
julia> dsolve(eqn(α=>2); ics=Dict(f(0)=>1)) |> print # fill in parameter, initial condition
Eq(f(x), x^2 + 1)

julia> eqn = ∂(∂(f(x))) ~ -f(x); print(eqn)
Eq(Derivative(f(x), (x, 2)), -f(x))

julia> dsolve(eqn)
f(x) = C₁⋅sin(x) + C₂⋅cos(x)

julia> dsolve(eqn; ics = Dict(f(0)=>1, ∂(f)(0) => -1))
f(x) = -sin(x) + cos(x)

julia> eqn = ∂(∂(f(x))) - f(x) - exp(x);

julia> dsolve(eqn, ics=Dict(f(0) => 1, f(1) => Sym(1//2))) |> print # not just 1//2
Eq(f(x), (x/2 + (-exp(2) - 2 + E)/(-2 + 2*exp(2)))*exp(x) + (-E + 3*exp(2))*exp(-x)/(-2 + 2*exp(2)))
```

Systems

```jldoctest dsolve
julia> @syms x() y() t g
(x, y, t, g)

julia> ∂ = Differential(t)
Differential(t)

julia> eqns = [∂(x(t)) ~ y(t), ∂(y(t)) ~ x(t)]
2-element Vector{Sym}:
 Eq(Derivative(x(t), t), y(t))
 Eq(Derivative(y(t), t), x(t))

julia> dsolve(eqns)
2-element Vector{Sym}:
 Eq(x(t), -C1*exp(-t) + C2*exp(t))
  Eq(y(t), C1*exp(-t) + C2*exp(t))

julia> dsolve(eqns, ics = Dict(x(0) => 1, y(0) => 2))
2-element Vector{Sym}:
 Eq(x(t), 3*exp(t)/2 - exp(-t)/2)
 Eq(y(t), 3*exp(t)/2 + exp(-t)/2)

julia> eqns = [∂(∂(x(t))) ~ 0, ∂(∂(y(t))) ~ -g]
2-element Vector{Sym}:
  Eq(Derivative(x(t), (t, 2)), 0)
 Eq(Derivative(y(t), (t, 2)), -g)

julia> dsolve(eqns)  # can't solve for initial conditions though! (NotAlgebraic)
2-element Vector{Sym}:
           x(t) = C₁ + C₂⋅t
 Eq(y(t), C3 + C4*t - g*t^2/2)

julia> @syms t x() y()
(t, x, y)

julia> eq = (∂(x)(t) ~ x(t)*y(t)*sin(t), ∂(y)(t) ~ y(t)^2 * sin(t))
(Eq(Derivative(x(t), t), x(t)*y(t)*sin(t)), Eq(Derivative(y(t), t), y(t)^2*sin(t)))
```

```julia
julia> dsolve(eq)  # returns a set to be `collect`ed:
PyObject {Eq(x(t), -exp(C1)/(C2*exp(C1) - cos(t))), Eq(y(t), -1/(C1 - cos(t)))}
```

```julia
julia> dsolve(eq) |> collect
2-element Vector{Any}:
 Eq(x(t), -exp(C1)/(C2*exp(C1) - cos(t)))
               Eq(y(t), -1/(C1 - cos(t)))
```

"""
function dsolve(eqn, args...;
                ics::Union{Nothing, AbstractDict, Tuple}=nothing,
                kwargs...)
    if isa(ics, Tuple) # legacy
        _dsolve(eqn, args...; ics=ics, kwargs...)
    else
        sympy.dsolve(eqn, args...; ics=ics, kwargs...)
    end
end

rhs(x::SymbolicObject) = pycall_hasproperty(x, :rhs) ? x.rhs : x
lhs(x::SymbolicObject) = pycall_hasproperty(x, :lhs) ? x.lhs : x

export dsolve, rhs, lhs

## ---- deprecate ----
## used with ics=(u,0,1) style
function _dsolve(eqn::Sym, args...; ics=nothing, kwargs...)

    Base.depwarn("Use of tuple(s), `(u, x₀, u₀)`, to specify initial conditions is deprecated. Use a dictionary: `ics=Dict(u(x₀) => u₀)`.", :_dsolve)

    if isempty(args)
        var = first(free_symbols(eqn))
    else
        var = first(args)
    end
    # var might be f(x) or x, we want `x`
    if Introspection.classname(var) != "Symbol"
        var = first(var.args)
    end
    ## if we have one initial condition, can be passed in a (u,x0,y0) *or* ((u,x0,y0),)
    ## if more than oneq a tuple of tuples
    if eltype(ics) <: Tuple
        __dsolve(eqn, var, ics; kwargs...)
    else
        __dsolve(eqn, var, (ics,); kwargs...)
        end
end

function __dsolve(eqn::Sym, var::Sym, ics; kwargs...)
    if length(ics) == 0
        throw(ArgumentError("""Some initial value specification is needed.
Specifying the function, as in `dsolve(ex, f(x))`, is deprecated.
Use `sympy.dsolve(ex, f(x); kwargs...)` directly for that underlying interface.
"""))
    end

    out = sympy.dsolve(eqn; kwargs...)
    ord = sympy.ode_order(eqn, var)

    ## `out` may be an array of solutions. If so we do each one.
    ## we want to use an array for output only if needed
    if !isa(out, Array)
        return _solve_ivp(out, var, ics,ord)
    else
        output = Sym[]
        for o in out
            a = _solve_ivp(o, var, ics,ord)
            a != nothing && push!(output, a)
        end
        return length(output) == 1 ? output[1] : output
    end
end

## Helper.
## out is an equation in var with constants. Args are intial conditions
## Return `nothing` if initial condition is not satisfied (found by `solve`)
function _solve_ivp(out, var, args, o)

    eqns = Sym[(diff(out.rhs(), var, f.n))(var=>x0) - y0 for (f, x0, y0) in args]
    sols = solve(eqns, Sym["C$i" for i in 1:o], dict=true)
    if length(sols) == 0
       return nothing
    end

    ## massage output
    ## Might have more than one solution, though unlikely. But if we substitute a variable
    ## for y0 we will get an array back from solve which may have length 1.
    if isa(sols, Array)
        if length(sols) == 1
            sols = sols[1]
        else
            return [out([Pair(k,v) for (k,v) in sol]...) for sol in sols]
        end
    end

    out([Pair(k,v) for (k,v) in sols]...)
end

## For System Of Ordinary Differential Equations
## may need to collect return values
# dsolve(eqs::Union{Array, Tuple}, args...; kwargs...) = sympy.dsolve(eqs, args...; kwargs...)


text:
import Base.(>)
#=
== : get a new Assoc where all of the elements of input Assoc mataches the given Element.
=#
function >(A::Assoc, E::Union{AbstractString,Number})
    tarIndex = searchsortedfirst(Val(A),E)
    if !(E == Val(A)[tarIndex])
        tarIndex = tarIndex-1

    end
    if (isa(E,Number) & (A.val ==[1.0])  )
        print("Hello!")
        tarIndex = E
    end
    ~, ~, valkey = findnz(A.A)
    mapping = find( x-> x > tarIndex, valkey)
    rows, cols, vals = find(A)
    return Assoc(rows[mapping],cols[mapping],vals[mapping])
end

>(E::Union{AbstractString,Number},A::Assoc) = (A < E)

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


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

@doc raw"""CustomResourceDefinitionStatus indicates the state of the CustomResourceDefinition

    IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionStatus(;
        acceptedNames=nothing,
        conditions=nothing,
        storedVersions=nothing,
    )

    - acceptedNames::IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionNames : acceptedNames are the names that are actually being used to serve discovery. They may be different than the names in spec.
    - conditions::Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionCondition} : conditions indicate state for particular aspects of a CustomResourceDefinition
    - storedVersions::Vector{String} : storedVersions lists all versions of CustomResources that were ever persisted. Tracking these versions allows a migration path for stored versions in etcd. The field is mutable so a migration controller can finish a migration to another version (ensuring no old objects are left in storage), and then remove the rest of the versions from this list. Versions may not be removed from `spec.versions` while they exist in this list.
"""
mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionStatus <: SwaggerModel
    acceptedNames::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionNames } # spec name: acceptedNames
    conditions::Any # spec type: Union{ Nothing, Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionCondition} } # spec name: conditions
    storedVersions::Any # spec type: Union{ Nothing, Vector{String} } # spec name: storedVersions

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionStatus(;acceptedNames=nothing, conditions=nothing, storedVersions=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionStatus, Symbol("acceptedNames"), acceptedNames)
        setfield!(o, Symbol("acceptedNames"), acceptedNames)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionStatus, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionStatus, Symbol("storedVersions"), storedVersions)
        setfield!(o, Symbol("storedVersions"), storedVersions)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionStatus

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionStatus = Dict{Symbol,Symbol}(Symbol("acceptedNames")=>Symbol("acceptedNames"), Symbol("conditions")=>Symbol("conditions"), Symbol("storedVersions")=>Symbol("storedVersions"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionStatus = Dict{Symbol,String}(Symbol("acceptedNames")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionNames", Symbol("conditions")=>"Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionCondition}", Symbol("storedVersions")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionStatus }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionStatus))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionStatus }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionStatus[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionStatus)
    (getproperty(o, Symbol("acceptedNames")) === nothing) && (return false)
    (getproperty(o, Symbol("storedVersions")) === nothing) && (return false)
    true
end

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


text:
# ---
# title: 380. Insert Delete GetRandom O(1)
# id: problem380
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Hash Table, Design
# link: <https://leetcode.com/problems/insert-delete-getrandom-o1/description/>
# hidden: true
# ---
#
# Implement the `RandomizedSet` class:
#
#   * `bool insert(int val)` Inserts an item `val` into the set if not present. Returns `true` if the item was not present, `false` otherwise.
#   * `bool remove(int val)` Removes an item `val` from the set if present. Returns `true` if the item was present, `false` otherwise.
#   * `int getRandom()` Returns a random element from the current set of elements (it's guaranteed that at least one element exists when this method is called). Each element must have the **same probability** of being returned.
#
# **Follow up:** Could you implement the functions of the class with each
# function works in **average** `O(1)` time?
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["RandomizedSet", "insert", "remove", "insert", "getRandom", "remove", "insert", "getRandom"]
#     [[], [1], [2], [2], [], [1], [2], []]
#     **Output**
#     [null, true, false, true, 2, true, false, 2]
#
#     **Explanation**
#     RandomizedSet randomizedSet = new RandomizedSet();
#     randomizedSet.insert(1); // Inserts 1 to the set. Returns true as 1 was inserted successfully.
#     randomizedSet.remove(2); // Returns false as 2 does not exist in the set.
#     randomizedSet.insert(2); // Inserts 2 to the set, returns true. Set now contains [1,2].
#     randomizedSet.getRandom(); // getRandom() should return either 1 or 2 randomly.
#     randomizedSet.remove(1); // Removes 1 from the set, returns true. Set now contains [2].
#     randomizedSet.insert(2); // 2 was already in the set, so return false.
#     randomizedSet.getRandom(); // Since 2 is the only number in the set, getRandom() will always return 2.
#
#
#
#
# **Constraints:**
#
#   * `-231 <= val <= 231 - 1`
#   * At most `105` calls will be made to `insert`, `remove`, and `getRandom`.
#   * There will be **at least one** element in the data structure when `getRandom` is called.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# perfctr_convenience.jl
using LIKWID
using LinearAlgebra

Marker.init()

A = rand(128, 64)
B = rand(64, 128)
C = zeros(128, 128)

for _ in 1:100
    @region "scalar_add" 1+2
end

@region "matmul" begin
    for _ in 1:100
        mul!(C, A, B)
    end
end

Marker.close()


text:
include("Map/Map.jl")
include("Map/ResidueRingPoly.jl")
include("Map/NfOrd.jl")
include("Map/NumberField.jl")
include("Map/AbGrp.jl")
include("Map/NfRel.jl")
include("Map/FqPolyRing.jl")