text:
# ------------------------------------------------------------------------------------------
# # Getting started
#
# Topics:
# 1.  How to print
# 2. How to assign variables
# 3. How to comment
# 4. Syntax for basic math
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## How to print
#
# In Julia we usually use `println()` to print
# ------------------------------------------------------------------------------------------

println("I'm excited to learn Julia!")

# ------------------------------------------------------------------------------------------
# ## How to assign variables
#
# All we need is a variable name, value, and an equal's sign!<br>
# Julia will figure out types for us.
# ------------------------------------------------------------------------------------------

my_answer = 42
typeof(my_answer)

my_pi = 3.14159
typeof(my_pi)

😺 = "smiley cat!"
typeof(😺)

# ------------------------------------------------------------------------------------------
# To type a smiley cat, use tab completion to select the emoji name and then tab again
# ------------------------------------------------------------------------------------------

# \:smi + <tab> --> select with down arrow + <enter> ---> <tab> + <enter> to complete

# ------------------------------------------------------------------------------------------
# After assigning a value to a variable, we can reassign a value of a different type to that
# variable without any issue.
# ------------------------------------------------------------------------------------------

😺 = 1

typeof(😺)

# ------------------------------------------------------------------------------------------
# Note: Julia allows us to write super generic code, and 😺 is an example of this.
#
# This allows us to write code like
# ------------------------------------------------------------------------------------------

😀 = 0
😞 = -1

😺 + 😞 == 😀

# ------------------------------------------------------------------------------------------
# ## How to comment
# ------------------------------------------------------------------------------------------

# You can leave comments on a single line using the pound/hash key

#=

For multi-line comments,
use the '#= =#' sequence.

=#

# ------------------------------------------------------------------------------------------
# ## Syntax for basic math
# ------------------------------------------------------------------------------------------

sum = 3 + 7

difference = 10 - 3

product = 20 * 5

quotient = 100 / 10

power = 10 ^ 2

modulus = 101 % 2

# ------------------------------------------------------------------------------------------
# ### Exercises
#
# #### 1.1
# Look up docs for the `convert` function.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### 1.2
# Assign `365` to a variable named `days`. Convert `days` to a float and assign it to
# variable `days_float`
# ------------------------------------------------------------------------------------------

@assert days == 365
@assert days_float == 365.0

# ------------------------------------------------------------------------------------------
# #### 1.3
# See what happens when you execute
#
# ```julia
# convert(Int64, "1")
# ```
# and
#
# ```julia
# parse(Int64, "1")
# ```
# ------------------------------------------------------------------------------------------

Please click on `Validate` on the top, once you are done with the exercises.

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using LinearAlgebra
using Statistics
using StatsBase
using Distributions
using ProgressMeter

"""
	metacommunity(X, nsim=200)

	Compute coherence, turnover, clumping of a metacommunity. Equivilant to `metacommunity(X,ord=1,model=5,Comm=0,iter=nsim,ax=1)` in EMS by CHRIS L. HIGGINS (Presley et al. 2010)

	`X` is assumed to be site-by-species incidence matrix. Site as row, species as col.
	`nsim` specifies numbers of simulation (null model generation).
	return values are identically named and ordered to EMS.
"""
function metacommunity(X, nsim=200)
	Abs, Apr, MA, SA = coherence(X, nsim)
	Re, Rpr, MR, SR = turnover(X, nsim)
	ret = [Abs, Apr, MA, SA, Re, Rpr, MR, SR, M, Mpr]
	return ret
end

function summarize(ret)
	Abs, Apr, MA, SA, Re, Rpr, MR, SR, M, Mpr = ret
	coh = MA - Abs
	clu = M - 1
	tur = Re - MR
	if coh > 0; print("POSITIVE coherence");else;print("NEGATIVE coherence");end
	print(", p = ");println(Apr);
	print("\tEmbedded absences: "); println(Abs)
	print("\tNull model expected: "); println(MA)
	if tur > 0; print("POSITIVE turnover");else;print("NEGATIVE turnover");end
	print(", p = ");println(Rpr);
	print("\tSpecies replacement: "); println(Re)
	print("\tNull model expected: "); println(MR)
	if clu > 0; print("POSITIVE clumping");else;print("NEGATIVE clumping");end
	print(", p = ");println(Mpr);
	print("\tMorista's I: "); print(M)
end

function coherence(X, nsim=200)
	X = makebinary(X)
	X = allzeros(X)
	Y = ordinate(X)
	Abs = embedabs(Y)
	null = randomnull_model5(X, nsim)
	EAbs = zeros(Int, nsim)
	EX = zeros(Int64, size(X))
	@showprogress 1 "examing coherence of null..." for i = 1:nsim
		EX = allzeros(@view null[:,:,i])
		EY = ordinate(EX)
		EAbs[i] = embedabs(EY)
	end
	MA = mean(EAbs)
	SA = std(EAbs)
	Apr = ztest(Abs, MA, SA, nsim)
	return Abs, Apr, MA, SA
end

function turnover(X, nsim=200)
	X = makebinary(X)
	X = allzeros(X)
	Y = ordinate(X)
	Y = boundfill_comm0(Y)
	Re = checkstat(Y)
	null = replacenull_comm0(Y, nsim)
	ERe = zeros(Int, nsim)
	@showprogress 1 "examing turnover of null... " for i = 1:nsim
		ERe[i] = checkstat(@view null[:,:,i])
	end
	MR = mean(ERe)
	SR = std(ERe)
	Rpr = ztest(Re, MR, SR, nsim)
	return Re, Rpr, MR, SR
end

function clumping(X)
	X = transpose(X)
	X = makebinary(X)
	X = allzeros(X)
	Y = ordinate(X)
	r, c = size(Y)
	ComBnd = zeros(c, 1)
	M = 0
	ComBndChi = 0
	@inbounds @simd for i = 1:r
		ind = findall(Y[i,:] .== 1)
		First = ind[1];
		Last = ind[length(ind)]
		for j = 1:c
			if First == j
				ComBnd[j] += 1
			end
			if Last == j
				ComBnd[j] += 1
			end
		end
	end
	TotComBnd = r * 2 - ComBnd[1] - ComBnd[c]
	ExpComBnd = TotComBnd / (c - 2)
	df = -1
	@inbounds @simd for i = 2:(c-1)
		M += ((ComBnd[i]/TotComBnd)*((ComBnd[i]-1)/(TotComBnd-1)))
		ComBndChi += ((ComBnd[i] - ExpComBnd)^2/ExpComBnd)
		df += 1
	end
	M *= c-2
	if M < 1
		Mpr = chi2test(ComBndChi, df, "left")
	else
		Mpr = chi2test(ComBndChi, df, "right")
	end
	return M, Mpr
end

# two-tailed z-test
function ztest(mu, xbar, stddev, n)
    z = (xbar-mu)/(stddev/sqrt(n))
    p = 2 * cdf(Normal(), -abs(z))
    return p
end

function chi2test(x, v, tail="left")
	if tail == "left"
		p = cdf(Chisq(v), x)
	elseif tail == "right"
		p = 1 - cdf(Chisq(v), x)
	end
	return p
end

function boundfill_comm0(X)
	sites, species = size(X)
	Y = copy(X)
	for i = 1:species
		ind = findall(X[:,i] .== 1)
		First = ind[1]
		Last = ind[length(ind)]
		Y[First:Last,i] .= 1
	end
	return Y
end

function checkstat(X)
	r, c = size(X)
	C = 0
	Xij = zeros(Int64, r)
	@inbounds @simd for i = 1:(c-1)
		Xi = @view X[:,i]
		R1 = sum(Xi)
		for j = (i+1):c
			Xj = @view X[:,j]
			S = sum(Xi .& Xj)
			R2 = sum(Xj)
			C += (R1-S)*(R2-S)
		end
	end
	return C
end

function embedabs(X)
	sites, species = size(X)
	FirstSpecies, LastSpecies, FirstSite, LastSite = boundfind(X)
	Abs = 0
	for i = 1:sites
		for j = 1:species
			if X[i,j] == 0
				flag = 0
				if j > FirstSpecies[i] && j < LastSpecies[i]
					Abs += 1
					flag = 1
				end
				if flag == 0 && i > FirstSite[j] && i < LastSite[j]
					Abs += 1
				end
			end
		end
	end
	return Abs
end

function boundfind(X)
	r,c = size(X)
	FirstSpecies = zeros(r,1)
    LastSpecies = zeros(r,1)
    FirstSite = zeros(1,c)
	LastSite = zeros(1,c)
	for i = 1:r
        ind = findall(X[i,:] .== 1)
        FirstSpecies[i] = ind[1];
        LastSpecies[i] = ind[length(ind)];
	end
	for i = 1:c
        ind = findall(X[:,i] .== 1);
        FirstSite[i] = ind[1];
        LastSite[i] = ind[length(ind)];
	end
	return FirstSpecies, LastSpecies, FirstSite, LastSite
end

# EMS randomnull model=5
function randomnull_model5(X, nsim)
    nrow, ncol = size(X)
    Incidence = sum(X, dims=1)
    Richness = sum(X, dims=2)
    FixRichness = ones(nrow, 1)
    PropIncidence = Incidence / sum(Incidence)
    P = FixRichness * PropIncidence
    null = zeros(UInt8, nrow, ncol, nsim)
    stop = 0
	wv = FrequencyWeights(@view P[1,:])
	pbar = Progress(nsim, 1, "generating random null...   ")
    @inbounds @simd for k = 1:nsim
        for i = 1:nrow
            wv = FrequencyWeights(@view P[i,:])
            for j = 1:Richness[i]
                stop = 0
                while stop == 0
                    ind = sample(1:ncol, wv)
                    if null[i, ind, k] != 1
                        null[i, ind, k] = 1
                        stop = 1
                    end
				end
            end
		end
		next!(pbar)
    end
    return null
end

function replacenull_comm0(X, nsim)
	sites, species = size(X)
	null = zeros(Int, sites, species, nsim)
	pbar = Progress(nsim, 1, "generating replace null...  ")
	@inbounds @simd for k = 1:nsim
		for i = 1:species
			len = sum(@view X[:,i])
			rd = Int(ceil(rand()*(sites-len+1)))
			null[rd:(rd+len-1),i,k] .= 1
		end
		next!(pbar)
	end
	return null
end

function makebinary(X)
	return Int.(X .> 0)
end

function allzeros(X)
	Rsum = sum(X, dims=2)
	Rind = findall(Rsum[:] .> 0)
	Y = @view X[Rind,:]
	Csum = sum(X, dims=1)
	Cind = findall(Csum[:] .> 0)
	Y = @view Y[:,Cind]
	return Y
end

function ordinate(X, ax=1)
	V, W = ca(X)
	ir = sortperm(@view V[:,ax])
	ic = sortperm(@view W[:,ax])
	Y = @view X[ir,ic]
	#Y = @view Y[size(Y,1):-1:1,:]
	return Y
end

function ca(X)
	tot = sum(X)
	Y = X / tot
	rw = sum(Y, dims=2)
	cw = sum(Y, dims=1)
	rc = rw * cw
	Y = (Y-rc) ./ sqrt.(rc)
	F = svd(Y, full=true)
	V = F.U ./ sqrt.(rw)
	# W = F.V ./ sqrt.(cw)
	W = transpose(F.Vt ./ sqrt.(cw))
	return V, W
end

function eig(X)
	vals, vecs = eigen(X)
	if !isreal(vals)
		mask = isreal.(vals)
		vals = Float64.(mask .* vals + abs.(.!mask .* vals))
	end
	if !isreal(vecs)
		mask = isreal.(vecs)
		vecs = Float64.(mask .* vecs + abs.(.!mask .* vecs))
	end
	ind = sortperm(-vals)
	vecs = @view vecs[:,ind]
	for j = 1:size(vecs, 2)
		if sum(@view vecs[:,j]) < 0
			vecs[:,j] = - @view vecs[:,j]
		end
	end
	return vecs
end

println()
println("########### fastEMS.jl ###########")
println("#       AUTHOR: YUANFEI PAN      #")
println("#         DATE: 2020-02-27       #")
println("########### fastEMS.jl ###########")
println()
println("INFO: Precompiling fastEMS library.")
println("INFO: This may take a few seconds.")
metacommunity(rand(0:1, 10, 10))
println("INFO: fastEMS library loaded.")
println("INFO: See ?metacommunity for help.")

ext: jl

max_stars_repo_name: Augustpan/fastEMS

max_issues_repo_name: Augustpan/fastEMS

max_forks_repo_name: Augustpan/fastEMS


text:
#!/usr/bin/env julia

using CSVFiles
using DataFrames
using MAT

function load_echoview_matrix(filename)

    df = CSVFiles.load(filename, header_exists=false,
                       skiplines_begin=1) |> DataFrame

    info("Transposing ..")
    A = transpose(convert(Array,df[:, 14:end]))

    info("Converting ..")
    A= Float64.(A)
end

function main(args)
    filenames = args

    if length(filenames) == 1
        dir = filenames[1]
        if isdir(dir)
            filenames =  filter(x->endswith(x,".csv"), readdir(dir))
            filenames = ["$(dir)$(x)" for x in filenames]
        end
    end

    for filename in filenames
        out ="$(basename(filename)).mat"

        if isfile(out)
            info("Skipping $out")
            continue
        end

        try
            info("Loading $filename ...")
            a = load_echoview_matrix(filename)

            dict = Dict()
            dict["data"] = a
            dict["DESCR"] = "Converted from $(basename(filename)) by EchoJulia."

            info("Writing $out ...")
            matwrite(out, dict)
        catch ex
            info("$ex")
        end
    end

end

main(ARGS)


text:
include("gammq.jl")

function smear1_ttm3(drsq::Float64, α_12::Float64, a::Float64)
    dd::Float64    = sqrt(drsq)
    dri::Float64   = 1.0/dd
    drsqi::Float64 = dri^2

    if (α_12 > eps(Float64))
        g23::Float64      = exp(gammln(2.0/3.0))
        AA::Float64       = α_12^(1.0/6.0)
        rA::Float64       = dd/AA
        rA3::Float64      = rA^3
        exp1::Float64     = exp(-a*rA3)
        a_cubert::Float64 = a^(1.0/3.0)

        ts0::Float64 = (1.0 - exp1 + a_cubert*rA*g23*gammq(2.0/3.0, a*rA3))*dri
        ts1::Float64 = (1.0 - exp1)*dri*drsqi
        return ts0, ts1
    end

    ts0 = dri
    ts1 = dri*drsqi
    return ts0, ts1
end

function smear2_ttm3(drsq::Float64, α_12::Float64, a::Float64)
    dd::Float64 = sqrt(drsq)

    if (α_12 > eps(Float64))
        exp1::Float64 = exp(-a*dd^3/sqrt(α_12))
        ts1::Float64 = (1.0 - exp1)/dd^3
        ts2::Float64 = (ts1 - exp1*a/(sqrt(α_12)))/drsq
        return ts1, ts2
    end
    ts1 = dd^-3
    ts2 = ts1/drsq
    return ts1, ts2
end

function smear3_ttm3(drsq::Float64, α_12::Float64, a::Float64)
    dd::Float64    = sqrt(drsq)

    if (α_12 > eps(Float64))
        exp1::Float64 = exp(-a*dd^3/sqrt(α_12))

        ts1::Float64 = (1.0 - exp1)/dd^3
        ts2::Float64 = (ts1 - exp1*a/(sqrt(α_12)))/drsq
        ts3::Float64 = (ts2 - 0.6*exp1*dd*a*a/α_12)/drsq
        return ts1, ts2, ts3
    end

    ts1 = dd^-3
    ts2 = ts1 / drsq
    ts3 = ts2 / drsq
    return ts1, ts2, ts3
end

function smear1_ttm4(drsq::Float64, α_12::Float64, a::Float64)
    dd::Float64    = sqrt(drsq)
    dri::Float64   = 1.0/dd
    drsqi::Float64 = dri^2

    if (α_12 > eps(Float64))
        g34::Float64   = exp(gammln(3.0/4.0))
        AA::Float64    = α_12^(1.0/6.0)
        rA::Float64    = dd/AA
        rA4::Float64   = rA^4
        exp1::Float64  = exp(-a*rA4)
        a_mrt::Float64 = a^(1.0/4.0)

        ts0 = (1.0 - exp1 + a_mrt*rA*g34*ttm::gammq(3.0/4.0, a*rA4))*dri
        ts1 = (1.0 - exp1)*dri*drsqi
        return ts0, ts1
    end
    ts0 = dri
    ts1 = dri*drsqi
    return ts0, ts1
end

function smear2_ttm4(drsq::Float64, α_12::Float64, a::Float64)
    dd::Float64    = sqrt(drsq)
    dri::Float64   = 1.0/dd
    drsqi::Float64 = dri^2

    if (α_12 > eps(Float64))
        AA::Float64   = α_12^(1.0/6.0)
        rA::Float64   = dd/AA
        rA4::Float64  = rA^4
        exp1::Float64 = exp(-a*rA4)

        ts1 = (1.0 - exp1)*dri*drsqi
        ts2 = (ts1 - (4.0/3.0)*a*exp1*rA4*dri*drsqi)*drsqi
        return ts1, ts2
    end
    ts1 = dri*drsqi
    ts2 = ts1*drsqi
    return ts1, ts2
end

function smear3_ttm4(drsq::Float64, α_12::Float64, a::Float64)
    dd::Float64    = sqrt(drsq)
    dri::Float64   = 1.0/dd
    drsqi::Float64 = dri^2

    if (α_12 > eps(Float64))
        AA::Float64   = α_12^(1.0/6.0)
        AA4::Float64  = AA^4
        rA::Float64   = dd/AA
        rA4::Float64  = rA^4
        exp1::Float64 = exp(-a*rA4)

        ts1 = (1.0 - exp1)*dri*drsqi
        ts2 = (ts1 - (4.0/3.0)*a*exp1*rA4*dri*drsqi)*drsqi
        ts3 = (ts2 - (4.0/15.0)*dri*a*(4*a*rA4 - 1.0)*exp1/AA4)*drsqi
        return ts1, ts2, ts3
    end
    ts1 = dri*drsqi
    ts2 = ts1*drsqi
    ts3 = ts2*drsqi
    return ts1, ts2, ts3
end

max_stars_repo_name: heindelj/MolecularPotentials

max_issues_repo_name: heindelj/MolecularPotentials

max_forks_repo_name: heindelj/MolecularPotentials


text:
path = "docs/src/assets/"

using ISA, Plots
ψ₀ = AMFMcomp(t->exp(-t^2),t->200.0,0.0)
ψ₁ = AMFMcomp(t->1.0,t->100*t+50,0.1)
ψ₂ = AMFMcomp(t->1.0+0.8*cos(2t),t->100 + 7.5*sin(t),π)
fs = 16_000
t = 0:1/fs:3
𝚿₀ = numComp( ψ₀(t), fs )
𝚿₁ = numComp( ψ₁(t), fs )
𝚿₂ = numComp( ψ₂(t), fs )
𝐳 = numModel([𝚿₀,𝚿₁,𝚿₂])
plot(𝐳)
png(path*"IS_num_exAMFMmodel1.png")


text:
abstract type EoS end
abstract type SAFT <: EoS end

abstract type PCSAFTFamily <: SAFT end
abstract type sPCSAFTFamily <: SAFT end
abstract type ogSAFTFamily <: SAFT end
abstract type SAFTVRMieFamily <: SAFT end

struct SAFTVRMie <: SAFTVRMieFamily;
    components
    params::SAFTVRMieParams
end

function SAFTVRMie(components::T1) where T1 <: AbstractVector{T2} where T2 <: AbstractString
    set_components = [Set([components[i]]) for i in 1:length(components)]
    raw_params = retrieveparams(components, "SAFTVRMie"; kwargs...)
    model = SAFTVRMie(set_components, create_SAFTVRMieParams(raw_params))
end

struct PCSAFT <: PCSAFTFamily
    components
    params::PCSAFTParams
end

function PCSAFT(components::T1) where T1 <: AbstractVector{T2} where T2 <: AbstractString
    set_components = [Set([components[i]]) for i in 1:length(components)]
    raw_params = retrieveparams(components, "PCSAFT"; kwargs...)
    model = PCSAFT(set_components, create_PCSAFTParams(raw_params))
end

struct sPCSAFT <: sPCSAFTFamily
    components
    params::sPCSAFTParams
end

function sPCSAFT(components::T1) where T1 <: AbstractVector{T2} where T2 <: AbstractString
    set_components = [Set([components[i]]) for i in 1:length(components)]
    raw_params = retrieveparams(components, "sPCSAFT"; kwargs...)
    model = sPCSAFT(set_components, create_sPCSAFTParams(raw_params))
end

struct ogSAFT <: ogSAFTFamily
    components
    params::ogSAFTParams
end

function ogSAFT(components::T1) where T1 <: AbstractVector{T2} where T2 <: AbstractString
    set_components = [Set([components[i]]) for i in 1:length(components)]
    raw_params = retrieveparams(components, "ogSAFT"; kwargs...)
    model = ogSAFT(set_components, create_ogSAFTParams(raw_params))
end


text:
if !(VERSION.major==0&&VERSION.minor<=6)
	using Test
else
	using Base.Test
end
using FinancialToolbox

print_colored("Starting Complex Number Test\n",:green)

#Test Parameters
testToll=1e-14;
spot=10;K=10;r=0.02;T=2.0;sigma=0.2;d=0.01;
assert_(value,toll)=@test abs(value)<toll

#EuropeanCall Option
PriceCall=blsprice(spot,K,r,T,sigma,d);
DeltaCall=blsdelta(spot,K,r,T,sigma,d);
ThetaCall=blstheta(spot,K,r,T,sigma,d);
LambdaCall=blslambda(spot,K,r,T,sigma,d);
RhoCall=blsrho(spot,K,r,T,sigma,d);

#EuropeanPut Option
PricePut=blsprice(spot,K,r,T,sigma,d,false);
DeltaPut=blsdelta(spot,K,r,T,sigma,d,false);
ThetaPut=blstheta(spot,K,r,T,sigma,d,false);
RhoPut=blsrho(spot,K,r,T,sigma,d,false);
LambdaPut=blslambda(spot,K,r,T,sigma,d,false);

#Equals for both Options
Gamma=blsgamma(spot,K,r,T,sigma,d);
Vega=blsvega(spot,K,r,T,sigma,d);

## Complex Test with Complex Step Approximation for European Call
#Test parameters
DerToll=1e-13;
di=1e-15;
df(f,x)=f(x+1im*di)/di;
#Function definition
Fcall1(spot)=blsprice(spot,K,r,T,sigma,d);
Gcall1(r)=blsprice(spot,K,r,T,sigma,d);
Hcall1(T)=blsprice(spot,K,r,T,sigma,d);
Lcall1(sigma)=blsprice(spot,K,r,T,sigma,d);
Pcall1(spot)=blsprice(spot,K,r,T,sigma,d)*spot.re/blsprice(spot.re,K,r,T,sigma,d);
#TEST
print_colored("--- European Call Sensitivities: Complex Step Approximation\n",:yellow)
print_colored("-----Testing Delta\n",:blue);
assert_(df(Fcall1,spot).im-DeltaCall,DerToll)
print_colored("-----Testing Rho\n",:blue);
assert_(df(Gcall1,r).im-RhoCall,DerToll)
print_colored("-----Testing Theta\n",:blue);
assert_(-df(Hcall1,T).im-ThetaCall,DerToll)
print_colored("-----Testing Lambda\n",:blue);
assert_(df(Pcall1,spot).im-LambdaCall,DerToll)

## Complex Test with Complex Step Approximation for European Put
#Function definition
Fput1(spot)=blsprice(spot,K,r,T,sigma,d,false);
Gput1(r)=blsprice(spot,K,r,T,sigma,d,false);
Hput1(T)=blsprice(spot,K,r,T,sigma,d,false);
Pput1(spot)=blsprice(spot,K,r,T,sigma,d,false)*spot.re/blsprice(spot.re,K,r,T,sigma,d,false);
#TEST
print_colored("--- European Put Sensitivities: Complex Step Approximation\n",:yellow)
print_colored("-----Testing Delta\n",:blue);
assert_(df(Fput1,spot).im-DeltaPut,DerToll)
print_colored("-----Testing Rho\n",:blue);
assert_(df(Gput1,r).im-RhoPut,DerToll)
print_colored("-----Testing Theta\n",:blue);
assert_(-df(Hput1,T).im-ThetaPut,DerToll)
print_colored("-----Testing Lambda\n",:blue);
assert_(df(Pput1,spot).im-LambdaPut,DerToll)

print_colored("-----Testing Vega\n",:blue);
assert_(df(Lcall1,sigma).im-Vega,DerToll)
print_colored("Complex Number Test Passed\n",:green)
println("")

#TEST OF INPUT VALIDATION
print_colored("Starting Input Validation Test Complex\n",:magenta)
print_colored("----Testing Negative  Spot Price \n",:cyan)
@test_throws(ErrorException, blsprice(-spot+1im,K,r,T,sigma,d))
@test_throws(ErrorException, blkprice(-spot+1im,K,r,T,sigma))
@test_throws(ErrorException, blsdelta(-spot+1im,K,r,T,sigma,d))
@test_throws(ErrorException, blsgamma(-spot+1im,K,r,T,sigma,d))
@test_throws(ErrorException, blstheta(-spot+1im,K,r,T,sigma,d))
@test_throws(ErrorException, blsrho(-spot+1im,K,r,T,sigma,d))
@test_throws(ErrorException, blspsi(-spot+1im,K,r,T,sigma,d))
@test_throws(ErrorException, blslambda(-spot+1im,K,r,T,sigma,d))
@test_throws(ErrorException, blsvanna(-spot+1im,K,r,T,sigma,d))
@test_throws(ErrorException, blsvega(-spot+1im,K,r,T,sigma,d))

print_colored("----Testing Negative  Strike Price \n",:cyan)
@test_throws(ErrorException, blsprice(spot,-K+1im,r,T,sigma,d))
@test_throws(ErrorException, blkprice(spot,-K+1im,r,T,sigma))
@test_throws(ErrorException, blsdelta(spot,-K+1im,r,T,sigma,d))
@test_throws(ErrorException, blsgamma(spot,-K+1im,r,T,sigma,d))
@test_throws(ErrorException, blstheta(spot,-K+1im,r,T,sigma,d))
@test_throws(ErrorException, blsrho(spot,-K+1im,r,T,sigma,d))
@test_throws(ErrorException, blspsi(spot,-K+1im,r,T,sigma,d))
@test_throws(ErrorException, blslambda(spot,-K+1im,r,T,sigma,d))
@test_throws(ErrorException, blsvanna(spot,-K+1im,r,T,sigma,d))
@test_throws(ErrorException, blsvega(spot,-K+1im,r,T,sigma,d))

print_colored("----Testing Negative  Time to Maturity \n",:cyan)
@test_throws(ErrorException, blsprice(spot,K,r,-T+1im,sigma,d))
@test_throws(ErrorException, blkprice(spot,K,r,-T+1im,sigma))
@test_throws(ErrorException, blsdelta(spot,K,r,-T+1im,sigma,d))
@test_throws(ErrorException, blsgamma(spot,K,r,-T+1im,sigma,d))
@test_throws(ErrorException, blstheta(spot,K,r,-T+1im,sigma,d))
@test_throws(ErrorException, blsrho(spot,K,r,-T+1im,sigma,d))
@test_throws(ErrorException, blspsi(spot,K,r,-T+1im,sigma,d))
@test_throws(ErrorException, blslambda(spot,K,r,-T+1im,sigma,d))
@test_throws(ErrorException, blsvanna(spot,K,r,-T+1im,sigma,d))
@test_throws(ErrorException, blsvega(spot,K,r,-T+1im,sigma,d))

print_colored("----Testing Negative  Volatility \n",:cyan)
@test_throws(ErrorException, blsprice(spot,K,r,T,-sigma+1im,d))
@test_throws(ErrorException, blkprice(spot,K,r,T,-sigma+1im))
@test_throws(ErrorException, blsdelta(spot,K,r,T,-sigma+1im,d))
@test_throws(ErrorException, blsgamma(spot,K,r,T,-sigma+1im,d))
@test_throws(ErrorException, blstheta(spot,K,r,T,-sigma+1im,d))
@test_throws(ErrorException, blsrho(spot,K,r,T,-sigma+1im,d))
@test_throws(ErrorException, blspsi(spot,K,r,T,-sigma+1im,d))
@test_throws(ErrorException, blslambda(spot,K,r,T,-sigma+1im,d))
@test_throws(ErrorException, blsvanna(spot,K,r,T,-sigma+1im,d))
@test_throws(ErrorException, blsvega(spot,K,r,T,-sigma+1im,d))

print_colored("Complex Input Validation Test Passed\n",:magenta)


text:
using Agents, Test, Random

model1 = ABM(Agent1, Space((3,3)))

agent = add_agent!((1,1), model1)
@test agent.pos == (1, 1)
@test agent.id == 1
pos1 = model1.space.agent_positions[coord2vertex((1,1), model1)]
@test length(pos1) == 1
@test pos1[1] == 1

move_agent!(agent, (2,2), model1)

@test agent.pos == (2,2)
pos1 = model1.space.agent_positions[coord2vertex((1,1), model1)]
@test length(pos1) == 0
pos2 = model1.space.agent_positions[coord2vertex((2,2), model1)]
@test pos2[1] == 1

# %% Scheduler tests
N = 1000

# by_id
model = ABM(Agent0; scheduler = by_id)
for i in 1:N; add_agent!(model); end
@test sort!(collect(keys(model.agents))) == 1:N
@test model.scheduler(model) == 1:N

# fastest
Random.seed!(12)
model = ABM(Agent0; scheduler = fastest)
for i in 1:N; add_agent!(model); end
@test sort!(collect(model.scheduler(model))) == 1:N

# random
Random.seed!(12)
model = ABM(Agent0; scheduler = random_activation)
for i in 1:N; add_agent!(model); end
@test model.scheduler(model)[1:3] == [913, 522, 637] # reproducibility test

# partial
Random.seed!(12)
model = ABM(Agent0; scheduler = partial_activation(0.1))
for i in 1:N; add_agent!(model); end

a = model.scheduler(model)
@test length(a) < N
@test a[1] == 74 # reproducibility test

# by property
model = ABM(Agent2; scheduler = property_activation(:weight))
for i in 1:N; add_agent!(model, rand()/rand()); end

Random.seed!(12)
a = model.scheduler(model)

ids = collect(keys(model.agents))
properties = [model.agents[id].weight for id in ids]

@test ids[sortperm(properties)] == a

@testset "sample!" begin
  model = ABM(Agent2)
  for i in 1:20; add_agent!(model, rand()/rand()); end
  allweights = [i.weight for i in values(model.agents)]
  mean_weights = sum(allweights)/length(allweights)
  sample!(model, 12, :weight)
  @test Agents.nagents(model) == 12
  allweights = [i.weight for i in values(model.agents)]
  mean_weights_new = sum(allweights)/length(allweights)
  @test mean_weights_new > mean_weights

  model.agents[1].weight = 2
  model.agents[2].weight = 3
  @test model.agents[1].weight == 2
  @test model.agents[2].weight == 3

  sample!(model, 40, :weight)
  @test Agents.nagents(model) == 40

  model2 = ABM(Agent5, Space((10, 10)))
  for i in 1:20; add_agent!(Agent5(i, i, rand()/rand()), model2); end
  allweights = [i.weight for i in values(model2.agents)]
  mean_weights = sum(allweights)/length(allweights)
  sample!(model2, 12, :weight)
  @test Agents.nagents(model2) == 12
  allweights = [i.weight for i in values(model2.agents)]
  mean_weights_new = sum(allweights)/length(allweights)
  @test mean_weights_new > mean_weights

  sample!(model2, 40, :weight)
  @test Agents.nagents(model2) == 40
end


text:
using Knet, ArgParse, JLD
using Knet: clear!

@knet function droplstm(x0; fbias=1, drop1=0.2, drop2=0.2, o...)
    x = drop(x0; pdrop=drop1)
    input  = add2(x,h; o..., f=sigm)
    forget = add2(x,h; o..., f=sigm, binit=Constant(fbias))
    output = add2(x,h; o..., f=sigm)
    newmem = add2(x,h; o..., f=tanh)
    ig = mul(input,newmem)
    fc = mul(forget,cell)
    cell = add(ig,fc)
    tc = tanh(cell)
    h  = mul(tc,output)
    h0 = drop(h; pdrop=drop2)
end

function ipa(args=ARGS)
    s = ArgParseSettings()
    @add_arg_table s begin
        ("--datafiles"; nargs='+'; required=true; help="First pair used for training, second pair if any for dev, others test")
        ("--dictfiles"; nargs='+'; required=true; help="Dictionary file pair")
        ("--loadfile"; help="model input file name")
        ("--bestfile"; help="model best file name")
        ("--savefile"; help="model output file name")
        ("--epochs"; arg_type=Int; default=30; help="number of epochs to train")
        ("--batchsize"; arg_type=Int; default=128; help="training minibatch size")
        ("--lossreport"; arg_type=Int; default=100000; help="frequency of loss output")
        ("--hidden"; arg_type=Int; default=512; help="hidden layer size")
        ("--lr"; arg_type=Float64; default=2.0; help="learning rate")
        ("--decay"; arg_type=Float64; default=1.0; help="lr *= decay if an epoch does not improve dev")
        ("--gclip"; arg_type=Float64; default=5.0; help="gradient clipping")
        ("--drop1"; arg_type=Float64; default=0.8; help="dropout probability for input embedding")
        ("--drop2"; arg_type=Float64; default=0.4; help="dropout probability for output embedding")
        ("--fbias"; arg_type=Float64; default=1.0; help="forget gate bias for lstm")
        ("--ftype"; default="Float32"; help="floating point type")
        # ("--winit"; arg_type=Float64; default=0.15)
        ("--usesparse"; action=:store_true; help="use sparse matrices if true")
        ("--fast"; help="skip norm and loss calculations."; action=:store_true)
        ("--gcheck"; arg_type=Int; default=0; help="perform gradient checking on gcheck samples if nonzero")
        ("--seed"; arg_type=Int; default=42; help="random seed")
    end
    isa(args, AbstractString) && (args=split(args))
    opts = parse_args(args, s)
    println(opts)
    for (k,v) in opts; @eval ($(symbol(k))=$v); end
    seed > 0 && setseed(seed)
    @assert length(dictfiles)==2
    global data = Any[]
    for i in 1:2:length(datafiles)
        push!(data, S2SData(datafiles[i], datafiles[i+1]; dict1=dictfiles[1], dict2=dictfiles[2], dense=!usesparse))
    end
    vocab = maxtoken(data[1],2)
    global model = loadfile == nothing ?
    S2S(droplstm; fbias=fbias, hidden=hidden, vocab=vocab, winit=Xavier(), drop1=drop1, drop2=drop2) :
    gpucopy(load(loadfile, "model"))
    lrate = lr; setopt!(model; lr=lrate)

    perp = zeros(length(data))
    (maxnorm,losscnt) = fast ? (nothing,nothing) : (zeros(2),zeros(2))
    lastloss = bestloss = Inf
    t0 = time_ns()
    for epoch=1:epochs
        fast || (fill!(maxnorm,0); fill!(losscnt,0))
        train(model, data[1], softloss; gclip=gclip, maxnorm=maxnorm, losscnt=losscnt, lossreport=lossreport) # use first pair for training
        fast || (perp[1] = exp(losscnt[1]/losscnt[2]))
        for d=2:length(data)
            loss = test(model, data[d], softloss)
            perp[d] = exp(loss)
            if d==2              # use second pair for validation
                bestfile!=nothing && loss < bestloss && (bestloss=loss; save(bestfile, "model", clear!(cpucopy(model)); compress=true))
                loss < lastloss || (@show lrate = decay*lrate; setopt!(model; lr=lrate))
                lastloss = loss
            end
        end
        println("epoch  secs    ptrain  ptest.. wnorm  gnorm")
        myprint(epoch, (time_ns()-t0)/1e9, perp..., (fast ? [] : maxnorm)...)
        gcheck > 0 && gradcheck(model, data[1], softloss; gcheck=gcheck)
    end
    savefile==nothing || save(savefile, "model", clear!(cpucopy(model)); compress=true)
    return (fast ? (perp...) :  (perp..., maxnorm...))
end

myprint(a...)=(for x in a; @printf("%-6g ",x); end; println(); flush(STDOUT))

!isinteractive() && !isdefined(:load_only) && ipa(ARGS)


text:
using LRUCache
export play, stream, stop, setup_sound, current_sound_latency, resume_sounds,
  pause_sounds, tick, sound_is_setup, clear_sound_cache

const weber_sound_version = 4

let
  version_in_file =
    match(r"libweber-sound\.([0-9]+)\.(dylib|dll)",weber_sound).captures[1]
  if parse(Int,version_in_file) != weber_sound_version
    error("Versions for sound driver do not match. Please run ",
          "Pkg.build(\"TimedPlayback\").")
  end
end

const default_sample_rate = 44100Hz

abstract type Hooks end
struct DefaultHooks <: Hooks end

mutable struct SoundSetupState
  samplerate::Freq{Int}
  cache::Bool
  playing::LRU{UInt,Sound}
  state::Ptr{Void}
  num_channels::Int
  queue_size::Int
  stream_unit::Int
  hooks::Hooks
end

const default_stream_unit = 2^12
const sound_setup_state =
    SoundSetupState(0Hz,false,LRU{UInt,Sound}(1),C_NULL,0,0,default_stream_unit,
                    DefaultHooks())
sound_is_setup() = sound_setup_state.samplerate != 0Hz

"""
With no argument samplerate reports the current playback sample rate, as
defined by [`setup_sound`](@ref).
"""
function samplerate(s::SoundSetupState=sound_setup_state)
  if s.samplerate == 0Hz
    default_sample_rate
  else
    s.samplerate
  end
end

# Give some time after the sound stops playing to clean it up.
# This ensures that even when there is some latency
# the sound will not be GC'ed until it is done playing.
const sound_cleanup_wait = 2

# register_sound: ensures that sounds are not GC'ed while they are
# playing. Whenever a new sound is registered it removes sounds that are no
# longer playing. This is called internally by all methods that send requests to
# play sounds to the weber-sound library (implemented in weber_sound.c)
function register_sound(current::Sound,done_at::Float64,wait=sound_cleanup_wait)
  setstate = sound_setup_state
  setstate.playing[object_id(current)] = current
end

show_latency_warnings() = show_latency_warnings(sound_setup_state.hooks)
show_latency_warnings(::Hooks) = false

function ws_if_error(msg)
  if sound_setup_state.state != C_NULL
    str = unsafe_string(ccall((:ws_error_str,weber_sound),Cstring,
                              (Ptr{Void},),sound_setup_state.state))
    if !isempty(str) error(msg*" - "*str) end

    str = unsafe_string(ccall((:ws_warn_str,weber_sound),Cstring,
                              (Ptr{Void},),sound_setup_state.state))
    if !isempty(str) && show_latency_warnings()
      warn(msg*" - "*str)
    end
  end
end

"""
    setup_sound(;[sample_rate=samplerate()],[num_channels=8],[queue_size=8],
                [stream_unit=2^11],[caching=false]
                [hooks=TimedPlayback.DefaultHooks()])

Initialize format and capacity of audio playback.

This function is called automatically (using the default settings) the first
time a `Sound` object is created (e.g. during [`play`](@ref)).  It need not
normally be called explicitly, unless you wish to change one of the default
settings.

# Sample Rate

Sample rate determines the maximum playable frequency (max freq is ≈
sample_rate/2). Changing the sample rate from the default 44100 to a new value
will also change the default sample rate sounds will be created at, to match
this new sample rate.

# Channel Number

The number of channels determines the number of sounds and streams that can be
played concurrently. Note that discrete sounds and streams use a distinct set of
channels.

# Queue Size

Sounds can be queued to play ahead of time (using the `time` parameter of
[`play`](@ref)). When you request that a sound be played it may be queued to
play on a channel where a sound is already playing. The number of sounds that
can be queued to play at once is determined by queue size. The number of
channels times the queue size determines the number of sounds that you can queue
up to play ahead of time.

# Stream Unit

The stream unit determines the number of samples that are streamed at one time.
If this value is too small for your hardware, streams will sound jumpy. However
the latency of streams will increase as the stream unit increases.

# Caching

If you enable caching, by default, loading or converting the same object into a
sound after the first time will result in the same sound object. You can always
disable caching for a specific call to generate a sound (see documentation for
[`sound`](@ref) and [`playable`](@ref)).

# Playback Hooks

There are a number of hooks into the playback mechanism
which you can use to customize playback. To do so, define a
child of the type `TimedPlayback.Hooks` and define as many of the following
methods on this type as you wish. All methods have default fallbacks.

* `show_latency_warnings(hooks)` - true if you wish warnings about slow
playback latency to be displayed. (Default is `false`)
* `on_high_latency(hooks,latency)` - called following a high latency warning;
  specified the latency of playback. (Default does nothing)
* `on_play(hooks)` - called each `play` is called (Default does nothing)
* `tick` - returns the current time in seconds since epochs (Default uses
  `TimedPlayback.precise_time()`)
* `on_no_timing(hooks)` - called when there is no timing specified during
   `play`. (Default does nothing)
* `get_streamers` - returns a dictionary of channel numbers to streamers. See
   documentation for `Streamer`. Currently undocumented.
   (Default manages streamers internally)
"""
function setup_sound(;sample_rate=samplerate(),
                     buffer_size=nothing,queue_size=8,num_channels=8,
                     stream_unit=default_stream_unit,
                     caching=false,
                     hooks=DefaultHooks())
  sound_setup_state.hooks = hooks
  sample_rate_Hz = inHz(Int,sample_rate)
  empty!(sound_cache)

  if sound_is_setup()
    ccall((:ws_close,weber_sound),Void,(Ptr{Void},),sound_setup_state.state)
    ws_if_error("While closing old audio stream during setup")
    ccall((:ws_free,weber_sound),Void,(Ptr{Void},),sound_setup_state.state)
  else
    if !sound_is_initialized[]
      sound_is_initialized[] = true
      atexit() do
        sleep(0.1)
        ccall((:ws_close,weber_sound),Void,
              (Ptr{Void},),sound_setup_state.state)
        ws_if_error("While closing audio stream at exit.")
        ccall((:ws_free,weber_sound),Void,
              (Ptr{Void},),sound_setup_state.state)
      end
    end
  end

  sound_setup_state.samplerate = sample_rate_Hz
  sound_setup_state.cache = caching
  sound_setup_state.state = ccall((:ws_setup,weber_sound),Ptr{Void},
                                  (Cint,Cint,Cint,),ustrip(sample_rate_Hz),
                                  num_channels,queue_size)
  sound_setup_state.playing =
    LRU{UInt,Sound}(2 * (queue_size*num_channels + 2*num_channels))
  sound_setup_state.num_channels = num_channels
  sound_setup_state.queue_size = queue_size
  sound_setup_state.stream_unit = stream_unit
  ws_if_error("While trying to initialize sound")
end
const sound_is_initialized = fill(false)

"""
    TimedPlayback.usecache()

Reports the default state of caching. If true, then sounds will be cached when
created through [`sound`](@ref) or [`playable`](@ref) by default.
"""
usecache() = sound_setup_state.cache

"""
    current_sound_latency()

Reports the current, minimum latency of audio playback.

The current latency depends on your hardware and software drivers. This
estimate does not include the time it takes for a sound to travel from
your sound card to speakers or headphones. This latency estimate is used
internally by [`play`](@ref) to present sounds at accurate times.
"""
function current_sound_latency()
  ccall((:ws_cur_latency,weber_sound),Cdouble,
        (Ptr{Void},),sound_setup_state.state)
end

"""
    play(x;[channel=0],[time=0s])

Plays a sound (created via [`sound`](@ref)).

For convenience, play can also can be called on any object that can be turned
into a sound (via `sound`).

If a time is specified, it indicates the amount of time since epoch that the
sound should start playing (see [`precise_time`](@ref)).

This function returns immediately with the channel the sound is playing on. You
may provide a specific channel that the sound plays on: only one sound can be
played per channel. Normally it is unecessary to specify a channel, because an
appropriate channel is selected for you. However, pausing and resuming of
sounds occurs on a per channel basis, so if you plan to pause a specific
sound, you can do so by specifying its channel.

# Streams

Play can also be used to present a continuous stream of sound.  In this case,
the channel defaults to channel 1 (there is no automatic selection of channels
for streams). Streams are usually created by specifying an infinite length
during sound generation using [`tone`](@ref), [`noise`](@ref),
[`harmonic_complex`](@ref) or [`audible`](@ref).
"""
function play(x;time=0.0s,channel=0)
  if !sound_is_setup()
    setup_sound()
  end

  on_play(sound_setup_state.hooks)
  play_(playable(x),ustrip(inseconds(time,samplerate(x))),channel)
end
on_play(::Hooks) = nothing

immutable WS_Sound
  buffer::Ptr{Void}
  len::Cint
end
WS_Sound{R}(x::Sound{R,Q0f15,2}) = WS_Sound(pointer(x.data),size(x,1))

tick(::Hooks) = precise_time()
tick(x::SoundSetupState) = tick(x.hooks)
tick() = tick(sound_setup_state)

"""
    at(offset)

Returns a precise time since epoch (in seconds) plus some offset.

Use this method to help specify exactly when a given sound should be played. You
can subsently add aditional offsets via
[`Unitful`](https://github.com/ajkeller34/Unitful.jl) values, e.g.

    sound1_time = at(10s) # 10 seconds from now
    sound2_time = sound1_time + 200ms # 10.2 seconds from now

In addition to adding a specificied onset, This differs from [`Base.time`](@ref)
in that it is generally more precise on windows machines (as of Julia v0.6), and
in that it makes use of the `Unitful` seconds type.

"""
at(time) = time + s*tick()

function play_{R}(x::Sound{R,Q0f15,2},time::Float64=0.0,channel::Int=0)
  if R != ustrip(samplerate())
    error("Sample rate of sound ($(R*Hz)) and audio playback ($(samplerate()))",
          " do not match. Please resample this sound by calling `resample` ",
          "or `playable`.")
  end
  if !(1 <= channel <= sound_setup_state.num_channels || channel <= 0)
    error("Channel $channel does not exist. Must fall between 1 and",
          " $(sound_setup_state.num_channels)")
  end

  # first, verify the sound can be played when we want to
  if time > 0.0
    latency = current_sound_latency()
    now = tick()
    if now + latency > time && show_latency_warnings()
      if latency > 0
        warn("Requested timing of sound cannot be achieved. ",
             "With your hardware you cannot request the playback of a sound ",
             "< $(round(1000*latency,2))ms before it begins.")
      else
        warn("Requested timing of sound cannot be achieved. ",
             "Give more time for the sound to be played.")
      end
      on_high_latency(sound_setup_state.hooks,(now + latency) - time)
    end
  else
    on_no_timing(sound_setup_state.hooks)
  end

  # play the sound
  channel = ccall((:ws_play,weber_sound),Cint,
                  (Cdouble,Cdouble,Cint,Ref{WS_Sound},Ptr{Void}),
                  tick(),time,channel-1,WS_Sound(x),
                  sound_setup_state.state) + 1
  ws_if_error("While playing sound")
  register_sound(x,(time > 0.0 ? time : tick()) +
                 ustrip(duration(x)))

  channel
end

on_high_latency(::Hooks,latency) = nothing
on_no_timing(::Hooks) = nothing

"""
    play(fn::Function)

Play the sound that's returned by calling `fn`.
"""
function play(fn::Function;keys...)
  play(fn();keys...)
end

"""
    Streamer

TODO: document streamer implementation
"""
mutable struct Streamer
  next_stream::Float64
  channel::Int
  stream::AbstractStream
  cache::Nullable{Sound}
  done_at::Float64
  start_at::Float64
end

const num_channels = 8
const streamers = Dict{Int,Streamer}()

function setup_streamers()
  empty = EmptyStream{ustrip(samplerate()),Q0f15}()
  streamers[-1] = Streamer(0.0,1,empty,Nullable(),0.0,0.0)
  Timer(t -> map(process,values(streamers)),1/60,1/60)
end

function get_streamers(::Hooks)
  if isempty(streamers)
    setup_streamers()
  end
  streamers
end

function play_{R}(stream::AbstractStream{R},time::Float64=0.0,channel::Int=1)
  channel = channel == 0 ? 1 : channel
  @assert 1 <= channel <= sound_setup_state.num_channels
  if R != ustrip(samplerate())
    error("Sample rate of sound ($(R*Hz)) and audio playback ($(samplerate()))",
          " do not match. Please resample this sound by calling `resample`.")
  end

  cur_streamers = get_streamers(sound_setup_state.hooks)
  unit_s = sound_setup_state.stream_unit / R

  if channel in keys(cur_streamers)
    streamer = cur_streamers[channel]

    if time > 0
      if streamer.done_at < time && show_latency_warnings()
        offset = time - streamer.done_at

        cat_stream = [limit(streamer.stream,offset*s); stream]
        cur_streamers[channel] =
          Streamer(streamer.next_stream,channel,cat_stream,Nullable(),
                   streamer.done_at,-1)
      else
        warn("Requested timing of stream cannot be achieved. ",
             "With the current streaming settings you cannot request playback ",
             "more than $(round(1000*unit_s,2))ms beforehand.")
        on_high_latency(sound_setup_state.hooks,streamer.done_at - time)
        cur_streamers[channel] = Streamer(tick(),channel,stream,Nullable(),
                                          streamer.done_at + unit_s,-1)
      end
    else
      cur_streamers[channel] = Streamer(tick(),channel,stream,Nullable(),
                                        streamer.done_at + unit_s,-1)
    end
  else
    cur_streamers[channel] = Streamer(tick(),channel,stream,Nullable(),
                                      tick()+unit_s,time)
  end
end

"""
    stop(channel)

Stop the stream that is playing on the given channel.
"""
function stop(channel::Int,stream::Bool=false)
  @assert 1 <= channel <= sound_setup_state.num_channels
  if stream
    delete!(get_streamers(sound_setup_state.hooks),channel)
  else
    ccall((:ws_stop,weber_sound),Void,(Ptr{Void},Cint),
          sound_setup_state.state,channel-1)
  end
  nothing
end

function process(streamer::Streamer)
  if streamer.stream isa EmptyStream
    return
  elseif done(streamer.stream) && isnull(streamer.cache)
    stop(streamer.channel)
    return
  end

  toplay = if !isnull(streamer.cache) get(streamer.cache) else
    x = sound(streamer.stream,sound_setup_state.stream_unit)
    result = playable(x)
    streamer.cache = Nullable(result)
    result
  end

  done_at = ccall((:ws_play_next,weber_sound),Cdouble,
                  (Cdouble,Cdouble,Cint,Ref{WS_Sound},Ptr{Void}),
                  tick(),streamer.start_at,streamer.channel-1,WS_Sound(toplay),
                  sound_setup_state.state)

  ws_if_error("While playing sound")
  if done_at < 0
    # sound not ready to be queued for playing, wait a bit and try again
    streamer.next_stream += ustrip(0.05duration(toplay))
  else
    # sound was queued to play, wait until this queued sound actually
    # starts playing to queue the next stream unit
    streamer.cache = Nullable()
    register_sound(toplay,done_at,4sound_setup_state.stream_unit / samplerate())
    streamer.next_stream += ustrip(0.75duration(toplay))
    streamer.done_at = done_at
    streamer.start_at = -1
  end
end

"""
    pause_sounds([channel],[isstream=false])

Pause all sounds (or a stream) playing on a given channel.

If no channel is specified, then all sounds are paused.
"""
function pause_sounds(channel=-1,isstream=false)
  if sound_is_setup()
    @assert 1 <= channel <= sound_setup_state.num_channels || channel <= 0
    ccall((:ws_pause,weber_sound),Void,(Ptr{Void},Cint,Cint,Cint),
          sound_setup_state.state,channel-1,isstream,true)
    ws_if_error("While pausing sounds")
  end
end

"""
    resume_sounds([channel],[isstream=false])

Resume all sounds (or a stream) playing on a given channel.

If no channel is specified, then all sounds are resumed.
"""
function resume_sounds(channel=-1,isstream=false)
  if sound_is_setup()
    @assert 1 <= channel <= sound_setup_state.num_channels || channel <= 0
    ccall((:ws_pause,weber_sound),Void,(Ptr{Void},Cint,Cint,Cint),
        sound_setup_state.state,channel-1,isstream,false)
    ws_if_error("While resuming audio playback")
  end
end


text:
using Symbolics
using SymbolicUtils
using RuntimeGeneratedFunctions
# using ModelingToolkit
using Unitful, UnitfulEquivalences
using LinearAlgebra
using StaticArrays
using SparseArrays
using Tullio
using Symbolics: get_variables
using ChainRules
using Zygote
using ForwardDiff
using Symbolics: destructure_arg, DestructuredArgs, make_array, _make_array, SerialForm, Func, toexpr, _build_and_inject_function, @__MODULE__, MultithreadedForm, tosymbol, Sym
using SymbolicUtils.Code: MakeArray
using Rotations
import Symbolics: substitute, simplify
Symbolics.substitute(A::AbstractArray{Num},d::Dict) = Symbolics.substitute.(A,(d,))
Symbolics.simplify(A::AbstractArray{Num}) = Symbolics.simplify.(A)

# short term manual fix for bug in SymbolicUtils.Code.create_array(::Type{<:Matrix},...) methods
# see https://github.com/JuliaSymbolics/SymbolicUtils.jl/issues/270
# import SymbolicUtils.Code: create_array
# @inline function SymbolicUtils.Code.create_array(::Type{<:Matrix}, ::Nothing, ::Val{dims}, elems...) where dims
# 	rows = Tuple(repeat([dims[2],],dims[1]))
# 	Base.hvcat(rows, elems...)
# end
# @inline function SymbolicUtils.Code.create_array(::Type{<:Matrix}, T, ::Val{dims}, elems...) where dims
# 	rows = Tuple(repeat([dims[2],],dims[1]))
#     Base.typed_hvcat(T, rows, elems...)
# end

# adjoint/rrule for SymbolicUtils.Code.create_array
# https://github.com/JuliaSymbolics/SymbolicUtils.jl/pull/278/files
function ChainRulesCore.rrule(::typeof(SymbolicUtils.Code.create_array), A::Type{<:AbstractArray}, T, u::Val{j}, d::Val{dims}, elems...) where {dims, j}
  y = SymbolicUtils.Code.create_array(A, T, u, d, elems...)
  function create_array_pullback(Δ)
    dx = Δ
    (NO_FIELDS, DoesNotExist(), DoesNotExist(), DoesNotExist(), DoesNotExist(), dx..., ntuple(_ -> DoesNotExist(), length(elems) - prod(dims) + j)...)
  end
  y, create_array_pullback
end

# add Symbolics.get_variables for arrays of `Num`s
import Symbolics.get_variables
function Symbolics.get_variables(A::AbstractArray{Num})
	unique(vcat(get_variables.(A)...))
end

RuntimeGeneratedFunctions.init(@__MODULE__)

# add minimal Unitful+Symbolics interoperability
import Base:*
*(x::Unitful.AbstractQuantity,y::Num) =  Quantity(x.val*y, unit(x))
*(y::Num,x::Unitful.AbstractQuantity) = x*y

# rrules
ChainRulesCore.rrule(T::Type{<:SArray}, xs::Number...) = ( T(xs...), dv -> (NO_FIELDS, dv...) )
ChainRulesCore.rrule(T::Type{<:SArray}, x::AbstractArray) = ( T(x), dv -> (NO_FIELDS, dv) )
ChainRulesCore.rrule(T::Type{<:SMatrix}, xs::Number...) = ( T(xs...), dv -> (NO_FIELDS, dv...) )
ChainRulesCore.rrule(T::Type{<:SMatrix}, x::AbstractMatrix) = ( T(x), dv -> (NO_FIELDS, dv) )
ChainRulesCore.rrule(T::Type{<:SVector}, xs::Number...) = ( T(xs...), dv -> (NO_FIELDS, dv...) )
ChainRulesCore.rrule(T::Type{<:SVector}, x::AbstractVector) = ( T(x), dv -> (NO_FIELDS, dv) )

function ChainRulesCore.rrule(::typeof(SparseArrays.SparseMatrixCSC),
    m::Integer, n::Integer, pp::Vector, ii::Vector, Av::Vector)
    A = SparseMatrixCSC(m,n,pp,ii,Av)
    function SparseMatrixCSC_pullback(dA::AbstractMatrix)
        dAv = Vector{eltype(dA)}(undef, length(Av))
        for j = 1:n, p = pp[j]:pp[j+1]-1
            dAv[p] = dA[ii[p],j]
        end
        return (NO_FIELDS, DoesNotExist(), DoesNotExist(), DoesNotExist(), DoesNotExist(), dAv)
    end
    function SparseMatrixCSC_pullback(dA::SparseMatrixCSC)
        @assert getproperty.(Ref(A), (:m,:n,:colptr,:rowval)) == getproperty.(Ref(dA), (:m,:n,:colptr,:rowval))
        return (NO_FIELDS, DoesNotExist(), DoesNotExist(), DoesNotExist(), DoesNotExist(), dA.nzval)
    end

    return A, SparseMatrixCSC_pullback
end

function mult(χ::AbstractArray{T,3},v₁::AbstractVector,v₂::AbstractVector) where T<:Real
	@tullio v₃[i] := χ[i,j,k] * v₁[j] * v₂[k]
end

function mult(χ::AbstractArray{T,4},v₁::AbstractVector,v₂::AbstractVector,v₃::AbstractVector) where T<:Real
	@tullio v₄[i] := χ[i,j,k,l] * v₁[j] * v₂[k] * v₃[l]
end

function Δₘ_factors(λs,ε_sym)
	λ = Num(first(get_variables(sum(ε_sym))))
	diagε_m1 = Vector(diag(ε_sym)) .- 1
	mapreduce(lm->substitute.( diagε_m1, ([λ=>lm],)), .*, λs)
end

function Δₘ(λs::AbstractVector,ε_sym, λᵣs::AbstractVector, χᵣ::AbstractArray{T,3}) where T
	dm = Δₘ_factors(λs,ε_sym) ./ Δₘ_factors(λᵣs,ε_sym)
	@tullio χ[i,j,k] := χᵣ[i,j,k] * dm[i] * dm[j] * dm[k] fastmath=true
end

function generate_array_fn(args::Vector{Num},A::AbstractMatrix; expr_module=@__MODULE__(), parallel=SerialForm())
	return fn = _build_and_inject_function(expr_module,toexpr(Func(args,[],make_array(parallel,args,A,Matrix))))
end

function generate_array_fn(arg::Num,A::AbstractMatrix; expr_module=@__MODULE__(), parallel=SerialForm())
	return fn = generate_array_fn([arg,], A; expr_module, parallel)
end

function generate_array_fn(args::Vector{Num},A::SArray; expr_module=@__MODULE__(), parallel=SerialForm())
	return fn = _build_and_inject_function(expr_module,toexpr(Func(args,[],make_array(parallel,args,A,SArray))))
end

function generate_array_fn(arg::Num,A::SArray; expr_module=@__MODULE__(), parallel=SerialForm())
	return fn = generate_array_fn([arg,], A; expr_module, parallel)
end

function generate_array_fn(args::Vector{Num},A::TA; expr_module=@__MODULE__(), parallel=SerialForm()) where TA<:AbstractArray
	return fn = _build_and_inject_function(expr_module,toexpr(Func(args,[],make_array(parallel,args,A,TA))))
end

function generate_array_fn(arg::Num,A::TA; expr_module=@__MODULE__(), parallel=SerialForm()) where TA<:AbstractArray
	return fn = generate_array_fn([arg,], A; expr_module, parallel)
end

"""
################################################################################
#																			   #
#							   	   Materials							   	   #
#																			   #
################################################################################
"""

abstract type AbstractMaterial end

struct Material <: AbstractMaterial
	models::Dict
	defaults::Dict
end

function get_model(mat::AbstractMaterial,model_name::Symbol,args...)
	model = mat.models[model_name]
	missing_var_defaults = filter(x->!in(first(x),tosymbol.(args)),mat.defaults)
	subs =  Dict([(Sym{Real}(k),v) for (k,v) in missing_var_defaults])
	if typeof(model)<:AbstractArray
		model_subs = substitute.(model, (subs,))
	else
		model_subs = substitute(model, subs)
	end
	return model_subs
end

function generate_fn(mat::AbstractMaterial,model_name::Symbol,args...)
	model = get_model(mat,model_name,args...)
	if typeof(model)<:AbstractArray
		fn = generate_array_fn([Num(Sym{Real}(arg)) for arg in args],model)
	else
		fn = build_function(model,args...;expression=Val{false})
	end
	return fn
end

"""
################################################################################
#																			   #
#							   	  Rotations							   		   #
#																			   #
################################################################################
"""

struct RotatedMaterial{TM,TR} <: AbstractMaterial
	parent::TM
	rotation::TR
	rotation_defaults::Dict
end

function rotate(χ::AbstractMatrix,𝓡::AbstractMatrix)
	@tullio χᵣ[i,j] := 𝓡[a,i] * 𝓡[b,j] * χ[a,b]  fastmath=true
end

function rotate(χ::AbstractArray{T,3},𝓡::AbstractMatrix) where {T<:Real}
	@tullio χᵣ[i,j,k] := 𝓡[a,i] * 𝓡[b,j] * 𝓡[c,k] * χ[a,b,c]  fastmath=true
end

function rotate(χ::AbstractArray{T,4},𝓡::TR) where {T<:Real, TR<:StaticMatrix{3,3}}
	@tullio χᵣ[i,j,k,l] := 𝓡[a,i] * 𝓡[b,j] * 𝓡[c,k] * 𝓡[d,l] * χ[a,b,c,d]  fastmath=true
end

rotate(χ::Real,𝓡::StaticMatrix{3,3}) = χ

function rotate(mat::TM,𝓡::TR) where {TM<:AbstractMaterial,TR<:AbstractMatrix}
	if eltype(𝓡)<:Num
		vars = get_variables(𝓡)
		defs = Dict{Symbol,Real}([ tosymbol(var) => 0.0 for var in vars])
	else
		defs = Dict{Symbol,Real}([])
	end
	RotatedMaterial{TM,TR}(mat,𝓡,defs)
end

function rotate(mat::TM,𝓡::TR,defs::Dict) where {TM<:AbstractMaterial,TR<:AbstractMatrix}
	RotatedMaterial{TM,TR}(mat,𝓡,defs)
end

function get_model(mat::RotatedMaterial,model_name::Symbol,args...)
	model = rotate(mat.parent.models[model_name],mat.rotation)
	defs = merge(mat.parent.defaults,mat.rotation_defaults)
	missing_var_defaults = filter(x->!in(first(x),tosymbol.(args)),defs)
	subs =  Dict([(Sym{Real}(k),v) for (k,v) in missing_var_defaults])
	if typeof(model)<:AbstractArray
		model_subs = substitute.(model, (subs,))
	else
		model_subs = substitute(model, subs)
	end
	return model_subs
end

"""
################################################################################
#																			   #
#					  Dispersion (group index, GVD) models					   #
#																			   #
################################################################################
"""

function ng_model(n_model::Num, λ::Num)
	Dλ = Differential(λ)
	return n_model - ( λ * expand_derivatives(Dλ(n_model)) )
end

function gvd_model(n_model::Num, λ::Num)
	Dλ = Differential(λ)
	return λ^3 * expand_derivatives(Dλ(Dλ(n_model)))
end

ng_model(n_model::AbstractArray{Num}, λ::Num) = ng_model.(n_model,(λ,))
gvd_model(n_model::AbstractArray{Num}, λ::Num) = gvd_model.(n_model,(λ,))

function ng_model(mat::AbstractMaterial; symbol=:λ)
	λ = Num(Sym{Real}(symbol))
	n_model = sqrt.(get_model(mat,:ε,symbol))
	return ng_model(n_model,λ)
end

function gvd_model(mat::AbstractMaterial; symbol=:λ)
	λ = Num(Sym{Real}(symbol))
	n_model = sqrt.(get_model(mat,:ε,symbol))
	return gvd_model(n_model,λ)
end

##

function n²_MgO_LiNbO₃_sym(λ, T; a₁, a₂, a₃, a₄, a₅, a₆, b₁, b₂, b₃, b₄, T₀)
    f = (T - T₀) * (T + T₀ + 2*273.16)  # so-called 'temperature dependent parameter'
    λ² = λ^2
    a₁ + b₁*f + (a₂ + b₂*f) / (λ² - (a₃ + b₃*f)^2) + (a₄ + b₄*f) / (λ² - a₅^2) - a₆*λ²
end

pₑ = (
    a₁ = 5.756,
    a₂ = 0.0983,
    a₃ = 0.202,
    a₄ = 189.32,
    a₅ = 12.52,
    a₆ = 1.32e-2,
    b₁ = 2.86e-6,
    b₂ = 4.7e-8,
    b₃ = 6.113e-8,
    b₄ = 1.516e-4,
    T₀ = 24.5,      # reference temperature in [Deg C]
)
pₒ = (
    a₁ = 5.653,
    a₂ = 0.1185,
    a₃ = 0.2091,
    a₄ = 89.61,
    a₅ = 10.85,
    a₆ = 1.97e-2,
    b₁ = 7.941e-7,
    b₂ = 3.134e-8,
    b₃ = -4.641e-9,
    b₄ = -2.188e-6,
    T₀ = 24.5,      # reference temperature in [Deg C]
)

pᵪ₂ = (
	d₃₃ =   20.3,    #   pm/V
	d₃₁ =   -4.1,    #   pm/V
	d₂₂ =   2.1,     #   pm/V
	λs  =  [1.313, 1.313, 1.313/2.0]
)

function make_LN(;pₒ=pₒ,pₑ=pₑ,pᵪ₂=pᵪ₂)
	@variables λ, T, λs[1:3]
	nₒ² = n²_MgO_LiNbO₃_sym(λ, T; pₒ...)
	nₑ² = n²_MgO_LiNbO₃_sym(λ, T; pₑ...)
	ε 	= diagm([nₒ², nₒ², nₑ²])
	d₃₃, d₃₁, d₂₂, λᵣs = pᵪ₂
	χ⁽²⁾ᵣ = cat(
		[ 	0.0	 	-d₂₂ 	d₃₁			#	xxx, xxy and xxz
		 	-d₂₂	0.0 	0.0			#	xyx, xyy and xyz
			d₃₁	 	0.0		0.0		],	#	xzx, xzy and xzz
		[ 	-d₂₂	0.0 	0.0			#	yxx, yxy and yxz
			0.0	 	d₂₂ 	d₃₁			#	yyx, yyy and yyz
			0.0	 	d₃₁		0.0		],	#	yzx, yzy and yzz
		[ 	d₃₁	 	0.0 	0.0			#	zxx, zxy and zxz
			0.0	 	d₃₁ 	0.0			#	zyx, zyy and zyz
			0.0	 	0.0 	d₃₃		],	#	zzx, zzy and zzz
		 dims = 3
	)
	nₒ = sqrt(nₒ²)
	ngₒ = ng_model(nₒ,λ)
	gvdₒ = gvd_model(nₒ,λ)
	nₑ = sqrt(nₑ²)
	ngₑ = ng_model(nₑ,λ)
	gvdₑ = gvd_model(nₑ,λ)
	models = Dict([
		:nₒ		=>	nₒ,
		:ngₒ	=>	ngₒ,
		:gvdₒ	=>	gvdₒ,
		:nₑ		=>	nₑ,
		:ngₑ	=>	ngₑ,
		:gvdₑ	=>	gvdₑ,
		:ng		=>	diagm([ngₒ, ngₒ, ngₑ]),
		:gvd	=>	diagm([gvdₒ, gvdₒ, gvdₑ]),
		:ε 		=> 	ε,
		:χ⁽²⁾	=>	SArray{Tuple{3,3,3}}(Δₘ(λs,ε, λᵣs, χ⁽²⁾ᵣ)),
	])
	defaults =	Dict([
		:λ		=>		0.8,	# μm
		:T		=>		24.5,	# °C
		:λs₁	=>		1.064,	# μm
		:λs₂	=>		1.064,	# μm
		:λs₃	=>		0.532,	# μm

	])
	Material(models, defaults)
end

LN = make_LN()

LN.models

LN.models[:ε]
get_model(LN,:ε, :λ)
feps1 = generate_fn(LN,:ε, :λ)
feps1(0.8)
feps1(0.7)
Zygote.gradient(x->sum(feps1(x)), 0.7)
ForwardDiff.derivative(x->sum(feps1(x)), 0.7)
Zygote.hessian(x->sum(feps1(x)), 0.7)
ForwardDiff.derivative(y->ForwardDiff.derivative(x->sum(feps1(x)),y), 0.7)
ForwardDiff.derivative(feps1,0.7)
ForwardDiff.derivative(x->ForwardDiff.derivative(feps1,x),0.7)

get_model(LN,:χ⁽²⁾, :λs₁, :λs₂, :λs₃, :T)
get_model(LN,:χ⁽²⁾)
get_model(LN,:χ⁽²⁾, :λs₁, :λs₂, :λs₃)
fchi21 = generate_fn(LN, :χ⁽²⁾, :λs₁, :λs₂, :λs₃)
Zygote.gradient(x->sum(fchi21(2x,2x,x)), 0.7)
ForwardDiff.derivative(x->sum(fchi21(2x,2x,x)), 0.7)
Zygote.hessian(x->sum(fchi21(2x,2x,x)), 0.7)
ForwardDiff.derivative(y->ForwardDiff.derivative(x->sum(fchi21(2x,2x,x)),y), 0.7)

LN.models[:nₒ]
fno = generate_fn(LN,:nₒ, :λ)
fno(0.8)
fno''(0.8)
Zygote.hessian(fno,0.8)

fno'.([0.7,0.8,0.9])
##
𝓡 = MRP(RotY(π/2)) |> Matrix
LNr = rotate(LN,𝓡)
𝓡2 = MRP(RotY(π/4)) |> Matrix
LNr2 = rotate(LN,𝓡2)
epsmr1 = get_model(LNr,:ε)
epsmr2 = get_model(LNr2,:ε)

epsmr1 = get_model(LNr,:ε,:λ)
epsmr2 = get_model(LNr2,:ε,:λ)
epsmr2 = get_model(LNr2,:ε,:λ,:T)

@variables θ
𝓡3 = RotY(θ) |> Matrix
LNr3 = rotate(LN,𝓡3)
LNr3.rotation

feps1 = generate_fn(LNr,:ε, :λ)
feps1(0.8)

chi2r1 = get_model(LNr,:χ⁽²⁾, :λs₁, :λs₂, :λs₃)
get_model(LNr,:χ⁽²⁾)[1,1,1]
fchi2r1 = generate_fn(LNr,:χ⁽²⁾, :λs₁, :λs₂, :λs₃)

fchi2r1(1.6,1.6,0.8)[1,1,1]
fchi2r1(1.064,1.064,0.532)[1,1,1]

Zygote.gradient(x->fchi2r1(x,x,x/2)[1,1,1],1.064)
ForwardDiff.derivative(x->fchi2r1(x,x,x/2)[1,1,1],1.064)
Zygote.hessian(x->fchi2r1(x,x,x/2)[1,1,1],1.064)
ForwardDiff.derivative(y->ForwardDiff.derivative(x->fchi2r1(x,x,x/2)[1,1,1],y),1.064)

feps3 = generate_fn(LNr3,:ε, :λ, :θ)

feps3(0.8,0.0)

feps3(0.8,π/2.0)

feps3

##

ng_fn(mat::AbstractMaterial; symbol=:λ)

ng_model(LNr)
@variables λ
Dλ = Differential(λ)
Dλ |> typeof
# ng(n_sym::Num) = n_sym - λ * expand_derivatives(Dλ(n_sym))
# gvd(n_sym::Num) = λ^3 * expand_derivatives(Dλ(Dλ(n_sym)))

nom = get_model(LN,:nₒ,:λ)
nem = get_model(LN,:nₑ,:λ)
typeof(nom)
Dλ(nom)|>expand_derivatives

ns = [nom, nom, nem]
Dλ.(ns)|>expand_derivatives
##

𝓡 = MRP(RotY(π/2))
𝓡
typeof(𝓡)<:AbstractMatrix
LN.models[:ε]

@tullio χᵣ[i,j] := 𝓡[a,i] * 𝓡[b,j] * (LN.models[:ε])[a,b]  fastmath=true

rotate((LN.models[:ε]),Matrix{3,3}(𝓡))
rotate((get_model(LN,:χ⁽²⁾)),Matrix{3,3}(𝓡))

##
ε_syms = tosymbol.(ε_vars)
ε_syms2 = tosymbol.(get_variables(LN.models[:ε]))
ε_syms[1] |> tosymbol

getindex.((LN.defaults,),ε_syms)

@variables x,y
z = 2x - y^2

substitute(z,Dict([Sym{Real}(:y)=>0.3]))

Num(:x)

##
function make_LN(;pₒ=pₒ,pₑ=pₑ,pᵪ₂=pᵪ₂)
	@variables λ, T
	nₒ² = n²_MgO_LiNbO₃_sym(λ, T; pₒ...)
	nₑ² = n²_MgO_LiNbO₃_sym(λ, T; pₑ...)
	ε 	= diagm([nₒ², nₒ², nₑ²])
	Material(
		ε,							# dielectric tensor
		[λ, T],						# variables
		[pₒ,pₑ,pᵪ₂],				# parameters
		[(0.0,8.0),(-20.0,300.0)],	# domains
		Dict([λ=>0.8,T=>pₒ.T₀]),	# Defaults
	)
end

LN = make_LN()
feps_LN = generate_ε(LN)
feps_LN(0.8)
feps_LN.([0.8,0.85,0.9])
ForwardDiff.derivative(feps_LN,0.7)
ForwardDiff.derivative(x->ForwardDiff.derivative(feps_LN,x),0.7)
Zygote.gradient(x->sum(feps_LN(x)),0.7)
ForwardDiff.derivative(x->sum(feps_LN(x)),0.7)

fchi2_LN = generate_χ⁽²⁾(LN)

fchi2_LN(1.0,1.0,0.5)
fchi2_LN(1.0,1.0,0.5)[3,3,3]

Zygote.gradient(x->sum(fchi2_LN(x,x,x/2)),1.0)
Zygote.hessian(x->sum(fchi2_LN(x,x,x/2)),1.0)
ForwardDiff.derivative(x->sum(fchi2_LN(x,x,x/2)),1.0)
ForwardDiff.derivative(x->fchi2_LN(x,x,x/2),1.0)
ForwardDiff.derivative(y->ForwardDiff.derivative(x->sum(fchi2_LN(x,x,x/2)),y),1.0)
ForwardDiff.derivative(y->ForwardDiff.derivative(x->fchi2_LN(x,x,x/2),y),1.0)
##

@variables λ, T, εᵣ[1:3,1:3], χ⁽²⁾

name 			=	:MgO_LiNbO₃
props 			=	[	:εᵣ, 		:χ⁽²⁾ 	]
prop_units 		= 	[	NoUnits, 	pm/V	]
params 			= 	[	:λ, 		:T    	]
param_units 	= 	[	μm, 		°C		]
param_defaults 	= 	[	0.8μm,		20°C 	]

prop_models		=	Dict([
						:εᵣ 	=> 	[ 	2.0*λ	0.0		0.0
						 				0.0		2.1*λ	0.0
										0.0		0.0		2.2*λ	],

						:χ⁽²⁾ => 	cat( [ 	4.0*λ	0.0		0.0
								 			0.0		4.1*λ	0.0
											0.0		0.0		2.2*λ	],
										[ 	4.0*λ+T	0.0		0.0
								 			0.0		4.1*λ	0.0
											0.0		0.0		2.2*λ	],
										dims=3,
										),
])

using Symbolics: get_variables, tosymbol

get_variables.(prop_models)

Dict( zip(keys(prop_models), tosymbol.(vars) for vars in (get_variables.(values(prop_models))) ) )

import Symbolics.get_variables
function Symbolics.get_variables(A::AbstractArray{Num})
	unique(vcat(get_variables.(A)...))
end

get_variables.(values(prop_models))

rhss(value.(εᵣ))

using Symbolics: rhss, lhss, value

# LN = Material(
# 	[	εᵣ	~	[	]
#
# 	]
#
# )

##
## try out Symbolics+Unitful (+FieldMetadata?)
@variables λ

λ

lm = 1.0u"μm"

uconvert(u"cm", lm)

ustrip(u"nm",lm)

lm - λ*1.0u"μm"

λ.val

using SymbolicUtils: symtype, @syms, getmetadata
symtype(λ.val)

typeof(λ.val)
tosymbol(λ)
getmetadata(λ.val)

u"cm" |> typeof <: Unitful.Units

lm
u"cm"(lm)
u"Hz"(lm,Spectral())

uconvert(u"Hz",lm,Spectral())

using Unitful.DefaultSymbols

using FieldMetadata
@metadata describe ""

@describe mutable struct Described
   a::Int     | "an Int with a description"
   b::Float64 | "a Float with a description"
end

d = Described(1, 1.0)

describe(d,:a)

using Parameters
@metadata describe "" String
@metadata bounds (0, 1) Tuple

@bounds @describe @with_kw struct WithKeyword{T}
    a::T = 3 | (0, 100) | "a field with a range, description and default"
    b::T = 5 | (2, 9)   | "another field with a range, description and default"
end

k = WithKeyword()

describe(k, :b)

bounds(k,:b)

##

d₃₃, d₃₁, d₂₂, λᵣs = LN.parameters[3]
chi2 = cat(
	[ 	0.0	 	-d₂₂ 	d₃₁			#	xxx, xxy and xxz
		-d₂₂	0.0 	0.0			#	xyx, xyy and xyz
		d₃₁	 	0.0		0.0		],	#	xzx, xzy and xzz
	[ 	-d₂₂	0.0 	0.0			#	yxx, yxy and yxz
		0.0	 	d₂₂ 	d₃₁			#	yyx, yyy and yyz
		0.0	 	d₃₁		0.0		],	#	yzx, yzy and yzz
	[ 	d₃₁	 	0.0 	0.0			#	zxx, zxy and zxz
		0.0	 	d₃₁ 	0.0			#	zyx, zyy and zyz
		0.0	 	0.0 	d₃₃		],	#	zzx, zzy and zzz
	 dims = 3
)

SArray{Tuple{3,3,3}}(chi2)

function gen_A()
	@variables a
	A	= 	[	a  			0.2
		 		-1.0*a   	2.0*a	]
	fA1 = build_function(A1,λ₁)[1]|>eval
end

##

@variables a
eps_LN = calculate_ε(a,LN)
f1 = generate_array_fn([a,],eps_LN)
f1(0.9)

f2 = _build_and_inject_function(Main,toexpr(Func([a,],[],make_array(SerialForm(),[a,],eps_LN,Matrix))))
f2(0.9)

f3 = _build_and_inject_function(Main,toexpr(Func([a,],[],MakeArray(eps_LN::Matrix,Matrix))))
f3(0.9)

using SymbolicUtils.Code: create_array

f4 = _build_and_inject_function(Main,toexpr(Func([a,],[],create_array(Matrix,nothing, Val{(3,3,3)}(),A...))))
create_array(Matrix, nothing, Val{(3, 3, 3)}(), A...)
f4(0.9)

sum()
Zygote.gradient(x->sum(f4(x)),0.9)

make_array(SerialForm(),[a,],eps_LN,Matrix)

using SymbolicUtils.Code: MakeArray

toexpr(Func([a,],[],make_array(SerialForm(),[a,],eps_LN,Matrix)))

make_array(SerialForm(),[a,],eps_LN,Matrix)))

MakeArray(eps_LN,Matrix)

## MBE
using Symbolics: @variables
using SymbolicUtils.Code: MakeArray, toexpr, create_array
@variables a
A	= 	[	a  		0		0
			0   	2*a		0
			0		0		3*a		]

A_expr = toexpr(MakeArray(A,Matrix))
eval(A_expr)

# ERROR: ArgumentError: argument count 9 does not match specified shape (2, 3)
# Stacktrace:
# {Number})
#    @ Base ./abstractarray.jl:1938
#  [2] hvcat(::Tuple{Int64, Int64}, ::Symbolics.Num, ::Vararg{Number})
#    @ Base ./abstractarray.jl:1925
#  [3] create_array(::Type{Matrix{T} where T}, ::Nothing, ::Val{(3, 3)}, ::Symbolics.Num, ::Vararg{Any})
#    @ SymbolicUtils.Code ~/.julia/packages/SymbolicUtils/fikCA/src/code.jl:411
#  [4] top-level scope
#    @ ~/.julia/packages/SymbolicUtils/fikCA/src/code.jl:374
#  [5] eval
# ...

##
using StaticArrays
using ChainRulesCore
using Zygote
using Symbolics: Num, make_array, SerialForm, Func, toexpr, _build_and_inject_function, @variables
using RuntimeGeneratedFunctions: @__MODULE__, init
init(@__MODULE__)
# generic rrules for StaticArray construction, might also want these for SVector and SMatrix
ChainRulesCore.rrule(T::Type{<:SArray}, xs::Number...) = ( T(xs...), dv -> (NO_FIELDS, dv...) )
ChainRulesCore.rrule(T::Type{<:SArray}, x::AbstractArray) = ( T(x), dv -> (NO_FIELDS, dv) )

# hacky methods bypassing Symbolics.build_function to explicitly specify `SArray` in `make_array`
function generate_array_fn(args::Vector{Num},A::SArray; expr_module=@__MODULE__(), parallel=SerialForm())
	return fn = _build_and_inject_function(expr_module,toexpr(Func(args,[],make_array(parallel,args,A,SArray))))
end

function generate_array_fn(arg::Num,A::SArray; expr_module=@__MODULE__(), parallel=SerialForm())
	return fn = generate_array_fn([arg,], A; expr_module, parallel)
end

function generate_A_fn()
	@variables a
	A	= 		cat(	[	a  		0		0
							0   	2*a		0
							0		0		3*a		],

						[	a^2  	0		0
							0   	sin(a)	0
							0		0		3.3		],
						dims=3,
			)
	return f_A = generate_array_fn(a,SArray{Tuple{3,3,2}}(A))
end
f_A = generate_A_fn()
f_A(0.3)

Zygote.gradient(x->sum(f_A(x)),0.3)

##
# M1 = [ 1.0a 2.0
# 	   3.0	4.0	]

fA1 = build_function(As,a)[1]|>eval
# fA2 = build_function(A,[M1,a])[1]|>eval
fA1(0.5)
Zygote.gradient(x->sum(fA1(x)),0.5)

args = [a,]
dargs = map(arg -> destructure_arg(arg, !checkbounds), [args...])
i = findfirst(x->x isa DestructuredArgs, dargs) # nothing

make_array(SerialForm(),[a,],A,Array)

make_array(SerialForm(),[a,],A,Matrix)

make_array(MultithreadedForm(),[a,],A,Matrix)

fA4 = _build_and_inject_function(Main,toexpr(Func([a,],[],make_array(SerialForm(),[a,],A,SMatrix)))) # == oop_expr2

fA5 = _build_and_inject_function(Main,toexpr(Func([a,],[],make_array(SerialForm(),[a,],A,SArray))))

fA5(0.7)
Zygote.refresh()
Zygote.gradient(x->sum(fA4(x)),0.7)
ForwardDiff.derivative(x->sum(fA5(x)),0.7)
ForwardDiff.derivative(fA4,0.7)

Zygote.gradient(x->sum(fA5(x)),0.7)
ForwardDiff.derivative(x->sum(fA5(x)),0.7)
ForwardDiff.derivative(fA5,0.7)

oop_expr1 = Func(dargs,[],make_array(SerialForm(),[a,],A,Array))

oop_expr2 = Func(dargs,[],make_array(SerialForm(),[a,],A,Matrix))

oop_expr3 = Func(dargs,[],make_array(SerialForm(),[a,],A,Matrix))

expr1 = oop_expr1 |> toexpr

expr2 = oop_expr2 |> toexpr

expr3 = oop_expr2 |> toexpr

MA1 = make_array(SerialForm(),[a,],A,Matrix)

function toexpr2(a::MakeArray, st)
    similarto = toexpr(a.similarto, st)
    T = similarto isa Type ? similarto : :(typeof($similarto))
    elT = a.output_eltype
    quote
        $create_array($T,
                     $elT,
                     Val{$(size(a.elems))}(),
                     $(map(x->toexpr(x, st), a.elems)...),)
    end
end

toexpr2(MA1,Matrix{Float64}())

expr2.head

expr2.args[1].head

expr_mod = @__MODULE__()

# using RuntimeGeneratedFunctions: _tagname,
# using RuntimeGeneratedFunctions

# expr1
(SymbolicUtils.Code.create_array)(Array, nothing, Val{(3, 3)}(), a, (*)(-1.0, a), 0.0, 0.2, (*)(2.0, a), 1.5, (*)(0.5, a), (*)(-2.3, a), (^)(a, 2))
(SymbolicUtils.Code.create_array)(Matrix{T} where T, nothing, Val{(3, 3, 3)}(), a, (*)(-1.0, a), 0.0, 0.2, (*)(2.0, a), 1.5, (*)(0.5, a), (*)(-2.3, a), (^)(a, 2))
#expr2 ^

create_array(Matrix, nothing, Val{(3, 3, 3)}(), a, (*)(-1.0, a), 0.0, 0.2, (*)(2.0, a), 1.5, (*)(0.5, a), (*)(-2.3, a), (^)(a, 2))
create_array(Matrix, nothing, Val{(3, 3, 3)}(), A...)

# expr1
(create_array2)(Array, nothing, Val{(3, 3)}(), a, (*)(-1.0, a), 0.0, 0.2, (*)(2.0, a), 1.5, (*)(0.5, a), (*)(-2.3, a), (^)(a, 2))
(create_array2)(Matrix, nothing, Val{(3, 3)}(), a, (*)(-1.0, a), 0.0, 0.2, (*)(2.0, a), 1.5, (*)(0.5, a), (*)(-2.3, a), (^)(a, 2))
#expr2 ^

els1 = (a, (*)(-1.0, a), 0.0, 0.2, (*)(2.0, a), 1.5, (*)(0.5, a), (*)(-2.3, a), (^)(a, 2))
hvcat((3,3,3),els1...)

@inline function create_array2(::Type{<:Matrix}, ::Nothing, ::Val{dims}, elems...) where dims
	println("")
	println("dims:")
	@show dims
	Base.hvcat(dims, elems...)
end

(SymbolicUtils.Code.create_array)(Matrix{T} where T, nothing, Val{(3, 3)}(), a, (*)(-1.0, a), 0.0, 0.2, (*)(2.0, a), 1.5, (*)(0.5, a), (*)(-2.3, a), (^)(a, 2))

module_tag = getproperty(@__MODULE__(), _tagname)

fA2 = _build_and_inject_function(Main,expr1)
fA3 = _build_and_inject_function(Main,expr2)

fA2(0.8)
fA3(0.8)

Zygote.gradient(x->sum(fA2(x)),0.8)
ForwardDiff.derivative(x->sum(fA2(x)),0.8)

##
eps_sym1 = calculate_ε(λ₁,LN)

using StaticArrays

#substitute.(A1, λ₁=>0.21)

fA1 = build_function(A1,λ₁)[1]|>eval
Zygote.gradient(x->sum(fA1(x)),0.5)

fA2 = build_function(A1::Matrix{Num},λ₁)[1]|>eval
Zygote.gradient(x->sum(fA2(x)),0.5)

A1s = SMatrix{2,2}(A1)
fA3 = build_function(A1s,λ₁)[1]|>eval
Zygote.gradient(x->sum(fA3(x)),0.5)

(SymbolicUtils.Code.create_array)(Array, nothing, Val{(2, 2)}(), (*)(3.2, λ₁), (*)(-1.2, λ₁), 0.2, λ₁)

(SymbolicUtils.Code.create_array)(Matrix, nothing, Val{(2, 2)}(), (*)(3.2, λ₁), (*)(-1.2, λ₁), 0.2, λ₁)

foo2(x) = (SymbolicUtils.Code.create_array)(Matrix, nothing, Val{(2, 2)}(), (*)(3.2, x), (*)(-1.2, x), 0.2, x)
foo2(3.3)
Zygote.gradient(x->sum(foo2(x)),3.3)
foo3(x) = (SymbolicUtils.Code.create_array)(Array, nothing, Val{(2, 2)}(), (*)(3.2, x), (*)(-1.2, x), 0.2, x)
foo3(3.3)
Zygote.gradient(x->sum(foo3(x)),3.3)

using SymbolicUtils
using SymbolicUtils.Code
toexpr.(eps_sym1)

eps_code1 = build_function(eps_sym1,λ₁)[1]

##
@variables λ, T
c = Unitful.c0      # Unitful.jl speed of light
Dλ = Differential(λ)
DT = Differential(T)
ng(n_sym::Num) = n_sym - λ * expand_derivatives(Dλ(n_sym))
gvd(n_sym::Num) = λ^3 * expand_derivatives(Dλ(Dλ(n_sym))) # gvd = uconvert( ( 1 / ( 2π * c^2) ) * _gvd(lm_um,T_C)u"μm", u"fs^2 / mm" )

##

s1 = λ + 3T

typeof(s1)

Symbolics.toexpr(s1)

nₒ_MgO_LiNbO₃_sym

simplify(nₒ_MgO_LiNbO₃_sym)

ngo = ng(nₒ_MgO_LiNbO₃_sym)

simplify(ngo)

gvdo = gvd(nₒ_MgO_LiNbO₃_sym)

epsLN = sparse(Diagonal([nₒ²_MgO_LiNbO₃_sym, nₒ²_MgO_LiNbO₃_sym, nₑ²_MgO_LiNbO₃_sym]))
χ⁽²⁾_MgO_LiNbO₃_sym(λ)
χ⁽²⁾_MgO_LiNbO₃_sym(0.8)

feps_oop = build_function(epsLN,λ, expression=Val{false})[1]
feps_ip = build_function(epsLN,λ, expression=Val{false})[2]

using ChainRules, Zygote, ForwardDiff

feps_oop(0.8)[1,1]
Zygote.refresh()
Zygote.gradient(x->feps_oop(x)[1,1],0.8)[1]
Zygote.hessian(x->feps_oop(x)[1,1],0.8)[1]
ForwardDiff.derivative(feps_oop,0.8)
ForwardDiff.derivative(x->ForwardDiff.derivative(feps_oop,x),0.8)
ForwardDiff.derivative(y->ForwardDiff.derivative(x->ForwardDiff.derivative(feps_oop,x),y),0.8)
##

LN2 = make_LN()
LN3 = make_LN()

LN4 = make_LN()

LN_defs = [var=>LN.defaults[var] for var in LN.variables[2:end]]
substitute.(LN.ε,LN_defs)

f_χ⁽²⁾_LN = eval(generate_χ⁽²⁾(LN4))
expression = Val{true}
f_chiLN = generate_χ⁽²⁾(LN)
f_chiLN2 = generate_χ⁽²⁾(LN2)
f_chiLN4 = generate_χ⁽²⁾(LN3)
f_chiLN5 = generate_χ⁽²⁾(LN2)
f_chiLN(1.064,1.064,0.532)
f_chiLN(1.064,1.064,0.532)[3,3,3]
Zygote.gradient(x->f_χ⁽²⁾_LN(x,x,x/2)[3,3,3],1.1)

Zygote.gradient(x->f_chiLN2(x,x,x/2)[3,3,3],1.1)
Zygote.gradient(x->f_chiLN5(x,x,x/2)[3,3,3],1.1)

@variables λ₁, λ₂, λ₃
χ⁽²⁾_sym = calculate_χ⁽²⁾([λ₁, λ₂, λ₃],LN)

get_variables(sum(χ⁽²⁾_sym))[1]

get_variables(sum(χ⁽²⁾_sym))[1] |> typeof
get_variables(sum(χ⁽²⁾_sym))[1] |> Symbolics.arguments
Num(get_variables(sum(χ⁽²⁾_sym))[1]) === λ₁
λ₁ |> typeof
λ₁ |> Symbolics.istree
χ⁽²⁾_sym |> Symbolics.istree

χ⁽²⁾_sym = calculate_χ⁽²⁾([λ₁, λ₂, λ₃],LN)
χ⁽²⁾_sym = calculate_χ⁽²⁾([1.064,1.064,0.532],LN)[3,3,3]

f_χ⁽²⁾_LN(1.064,1.064,0.532)[3,3,3]

substitute.(χ⁽²⁾_sym,[λ₁=>1.064, λ₂=>1.064, λ₃=>1.064/2])

Dict(zip(get_variables(sum(χ⁽²⁾_sym)),[1.064,1.064,0.532]))

LN.parameters[3][:d₃₃]
d₃₃3, d₃₁3, d₂₂3, λsᵣ3 = LN.parameters[3]
d₃₃3
d₃₁3
Symbolics.get_variables(sum(LN.ε))
Symbolics.get_varnumber(sum(LN.ε))
Symbolics.tosymbol.(Symbolics.get_variables(sum(LN.ε)))

##
pₑ_MgO_LiNbO₃ = (
    a₁ = 5.756,
    a₂ = 0.0983,
    a₃ = 0.202,
    a₄ = 189.32,
    a₅ = 12.52,
    a₆ = 1.32e-2,
    b₁ = 2.86e-6,
    b₂ = 4.7e-8,
    b₃ = 6.113e-8,
    b₄ = 1.516e-4,
    T₀ = 24.5,      # reference temperature in [Deg C]
)
pₒ_MgO_LiNbO₃ = (
    a₁ = 5.653,
    a₂ = 0.1185,
    a₃ = 0.2091,
    a₄ = 89.61,
    a₅ = 10.85,
    a₆ = 1.97e-2,
    b₁ = 7.941e-7,
    b₂ = 3.134e-8,
    b₃ = -4.641e-9,
    b₄ = -2.188e-6,
    T₀ = 24.5,      # reference temperature in [Deg C]
)

@variables λ, T

nₒ²_MgO_LiNbO₃_λT_sym = n²_MgO_LiNbO₃_sym(λ, T; pₒ_MgO_LiNbO₃...)
nₒ_MgO_LiNbO₃_λT_sym = sqrt(nₒ²_MgO_LiNbO₃_λT_sym)
nₒ²_MgO_LiNbO₃_sym = substitute(nₒ²_MgO_LiNbO₃_λT_sym,[T=>pₒ_MgO_LiNbO₃.T₀])
nₒ_MgO_LiNbO₃_sym = sqrt(nₒ²_MgO_LiNbO₃_sym)

nₑ²_MgO_LiNbO₃_λT_sym = n²_MgO_LiNbO₃_sym(λ, T; pₑ_MgO_LiNbO₃...)
nₑ_MgO_LiNbO₃_λT_sym = sqrt(nₑ²_MgO_LiNbO₃_λT_sym)
nₑ²_MgO_LiNbO₃_sym = substitute(nₑ²_MgO_LiNbO₃_λT_sym,[T=>pₑ_MgO_LiNbO₃.T₀])
nₑ_MgO_LiNbO₃_sym = sqrt(nₑ²_MgO_LiNbO₃_sym)

# 	Values to use for green-pumped processes
#    d₃₃ = 27

d₃₃ =   20.3    #   pm/V
d₃₁ =   -4.1    #   pm/V
d₂₂ =   2.1     #   pm/V

χ⁽²⁾ᵣ_MgO_LiNbO₃ = cat(
	[ 	0.0	 	-d₂₂ 	d₃₁			#	xxx, xxy and xxz
	 	-d₂₂	0.0 	0.0			#	xyx, xyy and xyz
		d₃₁	 	0.0		0.0		],	#	xzx, xzy and xzz
	[ 	-d₂₂	0.0 	0.0			#	yxx, yxy and yxz
		0.0	 	d₂₂ 	d₃₁			#	yyx, yyy and yyz
		0.0	 	d₃₁		0.0		],	#	yzx, yzy and yzz
	[ 	d₃₁	 	0.0 	0.0			#	zxx, zxy and zxz
		0.0	 	d₃₁ 	0.0			#	zyx, zyy and zyz
		0.0	 	0.0 	d₃₃		],	#	zzx, zzy and zzz
	 dims = 3
)
λs_χ⁽²⁾ᵣ_MgO_LiNbO₃ = [1.313,1.313,1.313/2.0]

nₒ²ᵣs = [ substitute(nₒ²_MgO_LiNbO₃_sym,[λ=>lm]).val for lm in λs_χ⁽²⁾ᵣ_MgO_LiNbO₃ ]
nₑ²ᵣs = [ substitute(nₑ²_MgO_LiNbO₃_sym,[λ=>lm]).val for lm in λs_χ⁽²⁾ᵣ_MgO_LiNbO₃ ]
εᵣs = [sparse(Diagonal([nosq, nosq, nesq])) for (nosq, nesq) in zip(nₒ²ᵣs,nₑ²ᵣs)]

ε_LN = sparse(Diagonal([nₒ²_MgO_LiNbO₃_sym, nₒ²_MgO_LiNbO₃_sym, nₑ²_MgO_LiNbO₃_sym]))
[eps->(eps-1) for eps in diag(ε_LN)]

@variables λ₁, λ₂
λ₃ = λ₁ + λ₂

[ substitute.( diag(ε_LN - 1I), ([λ=>lm],)) for lm in (λ₁,λ₂,λ₃) ]

reduce(.*,[ substitute.( diag(ε_LN - 1I), ([λ=>lm],)) for lm in (λ₁,λ₂,λ₃) ] )
diagε_m1 = Vector(diag(ε_LN)) .- 1
map(lm->substitute.( diagε_m1, ([λ=>lm],)), (λ₁,λ₂,λ₃) )
mapreduce(lm->substitute.( diagε_m1, ([λ=>lm],)), .*, (λ₁,λ₂,λ₃) )

Δₘ_factors([λ₁,λ₂,λ₃],ε_LN)
Δₘ_factors([1.313,1.313,1.313/2.0],ε_LN)

χᵣ = χ⁽²⁾ᵣ_MgO_LiNbO₃
dm = Δₘ_factors([λ₁,λ₂,λ₃],ε_LN) ./ Δₘ_factors([1.313,1.313,1.313/2.0],ε_LN)
@tullio χ[i,j,k] := χᵣ[i,j,k] * dm[i] * dm[j] * dm[k] fastmath=true

dm = Δₘ_factors([1.0,1.0,0.5],ε_LN) ./ Δₘ_factors([1.313,1.313,1.313/2.0],ε_LN)
@tullio χ[i,j,k] := χᵣ[i,j,k] * dm[i] * dm[j] * dm[k]
χ⁽²⁾ᵣ_MgO_LiNbO₃
χ⁽²⁾_MgO_LiNbO = Δₘ([λ₁,λ₂,λ₃],ε_LN,[1.313,1.313,1.313/2.0],χ⁽²⁾ᵣ_MgO_LiNbO₃)

Δₘ([1.0,1.0,0.5],ε_LN,[1.313,1.313,1.313/2.0],χ⁽²⁾ᵣ_MgO_LiNbO₃)
Δₘ_factors(,ε_LN)

substitute( diag(ε_LN - 1I), [λ=>λ₁])

Δₘ(λs, εs, λs_χ⁽²⁾ᵣ_MgO_LiNbO₃, εᵣs,χ⁽²⁾ᵣ_MgO_LiNbO₃)

χ⁽²⁾_MgO_LiNbO₃_sym(λ::Real) =  χ⁽²⁾_MgO_LiNbO₃_sym([λ,λ,λ/2])

using Tullio
function mult(χ::AbstractArray{T,3},v₁::AbstractVector,v₂::AbstractVector) where T<:Real
	@tullio v₃[i] := χ[i,j,k] * v₁[j] * v₂[k]
end

function mult(χ::AbstractArray{T,4},v₁::AbstractVector,v₂::AbstractVector,v₃::AbstractVector) where T<:Real
	@tullio v₄[i] := χ[i,j,k,l] * v₁[j] * v₂[k] * v₃[l]
end

function Δₘ(λs::AbstractVector, fε::Function, λᵣs::AbstractVector, χᵣ::AbstractArray{T,3}) where T
	dm = flat(map( (lm,lmr) -> (diag(fε(lm)).-1.) ./ (diag(fε(lmr)).-1.), λs, λᵣs ))
	@tullio χ[i,j,k] := χᵣ[i,j,k] * dm[i,1] * dm[j,2] * dm[k,3] fastmath=true
end

function Δₘ(λs::AbstractVector,ε_sym, λᵣs::AbstractVector, χᵣ::AbstractArray{T,3}) where T
	# dm = Symbolics.value.(Δₘ_factors(λs,ε_sym) ./ Δₘ_factors(λᵣs,ε_sym))
	dm = Δₘ_factors(λs,ε_sym) ./ Δₘ_factors(λᵣs,ε_sym)
	@tullio χ[i,j,k] := χᵣ[i,j,k] * dm[i] * dm[j] * dm[k] fastmath=true
end

##
struct Material{T}
	ε::T
	fε::Function
	fng::Function
	fgvd::Function
	# fχ⁽²⁾::Function
end

##

ngLN_sym = sparse(Diagonal(ng.([nₒ_MgO_LiNbO₃_sym,nₒ_MgO_LiNbO₃_sym,nₑ_MgO_LiNbO₃_sym])))

Num |> supertypes

ngLN_expr = build_function(ng(nₒ_MgO_LiNbO₃_sym), λ) #expresssion = false, parallel=Symbolics.MultithreadedForm())

ngLN = eval(build_function(ngLN_sym,
	λ;
	expresssion = Val{false}, )[1])
	# parallel=Symbolics.MultithreadedForm())[1])
ngLNfast = build_function(ngLN_sym,
	λ;
	expresssion = Val{false})[2] |> eval

ngLN

ngLN(0.9)

using ChainRules, Zygote, ForwardDiff

ForwardDiff.derivative(x->ngLN(x),0.9)
ForwardDiff.derivative(x->sum(ngLN(x)),0.9)
Zygote.gradient(x->sum(ngLN(x)),0.9)

struct Material{T}
	ε::T
	fε::Function
	fng::Function
	fgvd::Function
end

substitute(nₒ²_MgO_LiNbO₃_sym,[λ=>0.8])

@variables θ
using Rotations, StaticArrays, SparseArrays

R1 = RotY(0.1)

Matrix(RotY(θ))

SMatrix{3,3}(R1)

##
#
# # χ⁽²⁾ Nonlinear dielectric susceptibility tensor debug
#
# d₃₃ =   25.0    #   pm/V
# d₃₁ =   -4.1    #   pm/V
# d₂₂ =   2.1     #   pm/V
#
# #          xx      yy       zz      zy      zx      xy
# deff = [    0.      0.      0.      0.      d₃₁     -d₂₂     #   x
#             -d₂₂    d₂₂     0.      d₃₁     0.      0.       #   y
#             d₃₁     d₃₁     d₃₃     0.      0.      0.   ]   #   z
#
# # χ⁽²⁾ LiNbO₃ (3m) with extraordinary axis along ẑ (3)
# χ⁽²⁾₃ = cat(
# 	[ 	0.0	 	-d₂₂ 	d₃₁			#	xxx, xxy and xxz
# 	 	-d₂₂	0.0 	0.0			#	xyx, xyy and xyz
# 		d₃₁	 	0.0		0.0		],	#	xzx, xzy and xzz
# 	[ 	-d₂₂	0.0 	0.0			#	yxx, yxy and yxz
# 		0.0	 	d₂₂ 	d₃₁			#	yyx, yyy and yyz
# 		0.0	 	d₃₁		0.0		],	#	yzx, yzy and yzz
# 	[ 	d₃₁	 	0.0 	0.0			#	zxx, zxy and zxz
# 		0.0	 	d₃₁ 	0.0			#	zyx, zyy and zyz
# 		0.0	 	0.0 	d₃₃		],	#	zzx, zzy and zzz
# 	 dims = 3
# )
#
# function Δₘ(fε::Function,λs,λᵣs)
# 	# λ₁,λ₂,λ₃ 	= 	λs
# 	# λ₁ᵣ,λ₂ᵣ,λ₃ᵣ	= 	λᵣs
# 	#
# 	*.(diag.(fε.(λs)) .- 1.0)...) / *.((diag.(fε.(λᵣs)) .- 1.0)...)
# end
#
# ##
# function fε(λ::T) where T<:Real
#     nₑ² = nₑ²_MgO_LiNbO₃(λ)
#     nₒ² = nₒ²_MgO_LiNbO₃(λ)
#     # Diagonal( [ nₑ², nₒ², nₒ² ] )
#     SMatrix{3,3,T,9}( nₒ²,    0.,     0.,
#                       0.,     nₒ²,    0.,
#                       0.,     0.,     nₑ², )
# end
#
# fε(0.8)
#
# χ⁽²⁾ᵣ	= 	20.3 * (2.0/π)		#	[pm/V]
# λ₁ᵣ		= 	1.313				#	[μm]
# λ₂ᵣ		= 	λ₁ᵣ					#	[μm]
# λ₃ᵣ		= 	λ₁ᵣ/2				#	[μm]
#
# λ₁		= 	1.56				#	[μm]
# λ₂		= 	λ₁					#	[μm]
# λ₃		= 	λ₁/2				#	[μm]
#
# λs		=	[ 	λ₁,		λ₂,		λ₃	]
# λᵣs		=	[	λ₁ᵣ,	λ₂ᵣ,	λ₃ᵣ	]
#
# using Tullio
# dm0 = flat(map( (lm,lmr) -> (diag(fε(lm)).-1.) ./ (diag(fε(lmr)).-1.), λs, λᵣs ))
# @tullio χ⁽²⁾[i,j,k] := χ⁽²⁾₃[i,j,k] * dm0[i,1] * dm0[j,2] * dm0[k,3]
#
# function Δₘ(λs::AbstractVector, χᵣ::AbstractArray{T,3}, λᵣs::AbstractVector) where T
# 	dm = flat(map( (lm,lmr) -> (diag(fε(lm)).-1.) ./ (diag(fε(lmr)).-1.), λs, λᵣs ))
# 	@tullio χ[i,j,k] := χᵣ[i,j,k] * dm[i,1] * dm[j,2] * dm[k,3] verbose=true
# end
#
# χ⁽²⁾_LN(λs::AbstractVector) =  Δₘ(λs, χ⁽²⁾₃, [1.313,1.313,1.313/2.0])
# χ⁽²⁾_LN(λ::Real) =  Δₘ([λ,λ,λ/2.0], χ⁽²⁾₃, [1.313,1.313,1.313/2.0])
#
# χ⁽²⁾_MgO_LiNbO₃([1.064,1.064,0.532])[3,3,3]
# χ⁽²⁾_MgO_LiNbO₃(1.064) #[3,3,3]
# χ⁽²⁾_MgO_LiNbO₃([2.1,2.1,1.05])[3,3,3]
#
# Zygote.gradient(x->χ⁽²⁾_LN([x,x,x/2.0])[3,3,3],2.1)
# Zygote.gradient(x->χ⁽²⁾_LN(x)[3,3,3],2.1)[1]
# ForwardDiff.derivative(x->χ⁽²⁾_LN(x)[3,3,3],2.1)[1]
# Zygote.hessian(x->χ⁽²⁾_LN(x)[3,3,3],2.1)[1]
#
# using Rotations
# ForwardDiff.derivative(t->rotate(χ⁽²⁾_LN(1.064),RotY(t))[3,3,3],0.1)
# ForwardDiff.gradient(xx->rotate(χ⁽²⁾_LN(xx[1]),RotY(xx[2]))[3,3,3],[1.064,π/6.0])
# Zygote.gradient(xx->rotate(χ⁽²⁾_LN(xx[1]),Zygote.forwarddiff(RotY,xx[2]))[3,3,3],[1.064,π/6.0])
# Zygote.gradient((xx,tt)->rotate(χ⁽²⁾_LN(xx),Zygote.forwarddiff(RotY,tt))[3,3,3],1.064,π/6.0)
# Zygote.hessian((xx,tt)->rotate(χ⁽²⁾_LN(xx),Zygote.forwarddiff(RotY,tt))[3,3,3],1.064,π/6.0)
# Zygote.hessian(xx->rotate(χ⁽²⁾_LN(xx[1]),Zygote.forwarddiff(RotY,xx[2]))[3,3,3],[1.064,π/6.0])
#
# χ⁽²⁾_LN([1.05,2.1,2.1])[1,1,2]
# χ⁽²⁾_LN([2.1,2.1,1.05])[1,1,2]
#
# lms = collect(range(1.9,2.2,length=20))
# chi2LN333 = getindex.(χ⁽²⁾_LN.(lms),(3,),(3,),(3,))
#
# lineplot(lms,chi2LN333)
#
# chi2[3,3,3] ≈  χ⁽²⁾[3,3,3]/χ⁽²⁾₃[3,3,3]
#
# function mult(χ::AbstractArray{T,3},v₁::AbstractVector,v₂::AbstractVector) where T<:Real
# 	@tullio v₃[i] := χ[i,j,k] * v₁[j] * v₂[k]
# end
#
# function mult(χ::AbstractArray{T,4},v₁::AbstractVector,v₂::AbstractVector,v₃::AbstractVector) where T<:Real
# 	@tullio v₄[i] := χ[i,j,k,l] * v₁[j] * v₂[k] * v₃[l]
# end
#
#
#
#
#
#
#
#
#
#
#
# dm_check1 = ((nₑ²_MgO_LiNbO₃(1.56) - 1)^2 * (nₑ²_MgO_LiNbO₃(1.56/2.0) - 1)) / ((nₑ²_MgO_LiNbO₃(1.313) - 1)^2 * (nₑ²_MgO_LiNbO₃(1.313/2.0) - 1))
# dm_check1
#
# χ⁽²⁾[3,3,3] / χ⁽²⁾₃[3,3,3]
#
# dm_check2 = ((nₒ²_MgO_LiNbO₃(1.56) - 1)^2 * (nₒ²_MgO_LiNbO₃(1.56/2.0) - 1)) / ((nₒ²_MgO_LiNbO₃(1.313) - 1)^2 * (nₒ²_MgO_LiNbO₃(1.313/2.0) - 1))
# dm_check2
# χ⁽²⁾[1,2,1] / χ⁽²⁾₃[1,2,1]
#
#
# .*(dm0)
#
# (diag.(fε.(λs).-1.0I)) ./ (diag.(fε.(λᵣs).-1.0I))
#
# ##
#
# # χ⁽²⁾ LiNbO₃ (3m) with extraordinary axis along x̂ (1)
# χ⁽²⁾₁ = cat(
# 	[ 	d₃₃	 	0.0 	0.0			#	xxx, xxy and xxz
# 	 	0.0	 	d₃₁ 	0.0			#	xyx, xyy and xyz
# 		0.0	 	0.0		d₃₁		],	#	xzx, xzy and xzz
# 	[ 	0.0	 	d₃₁ 	0.0			#	yxx, yxy and yxz
# 		d₃₁	 	d₂₂ 	0.0			#	yyx, yyy and yyz
# 		0.0	 	0.0		-d₂₂	],	#	yzx, yzy and yzz
# 	[ 	0.0	 	0.0 	d₃₁			#	zxx, zxy and zxz
# 		0.0	 	0.0 	-d₂₂		#	zyx, zyy and zyz
# 		d₃₁	 	-d₂₂ 	0.0		],	#	zzx, zzy and zzz
# 	 dims = 3
# )
#
# E1 = randn(3)
# E2 = randn(3)
# E12 = [
# 	E1[1]*E2[1],
# 	E1[2]*E2[2],
# 	E1[3]*E2[3],
# 	E1[2]*E2[3] + E1[3]*E2[2],
# 	E1[1]*E2[3] + E1[3]*E2[1],
# 	E1[1]*E2[2] + E1[2]*E2[1],
# ]
# using Tullio
#
# @tullio P3[k] := χ⁽²⁾₃[i,j,k] * E1[i] * E2[j] verbose=1
# P3d = deff * E12
# @assert P3d ≈ P3
#
# using Rotations
# R1 = MRP(RotY(π/2))		# 90° rotation around ŷ
# R1i = transpose(R1)
# nₑ = 2.1
# nₒ = 2.2
# eps1 = [	nₑ^2	0.0		0.0
# 			0.0		nₒ^2	0.0
# 			0.0		0.0		nₒ^2	]
# eps2 = [	nₒ^2	0.0		0.0
# 			0.0		nₒ^2	0.0
# 			0.0		0.0		nₑ^2	]
# eps3 = diagm( [nₒ^2,nₒ^2,nₑ^2] )
#
# @assert R1 * eps1 / R1 ≈ eps2
# @assert R1 * eps3 / R1 ≈ eps1
#
# R2 = RotY(π/6)
# R2 * eps3 / R2
#
# @assert rotate(eps1,R1) ≈ MArray{Tuple{3,3}}(eps2)
#
#
# @tullio P11[i] := eps1[i,j] * E1[j] verbose=true
# @tullio P21[i] := eps2[i,j] * E1[j] verbose=true
# @tullio P212[i] := R1[i,j] * eps1[j,k] * R1i[k,h] * E1[h] verbose=true
# @assert P212 ≈ P21
#
# @tullio P3[i] := χ⁽²⁾₃[i,j,k] * E1[j] * E2[k] verbose=1
# P3d = deff * E12
# @assert P3d ≈ P3
# @tullio P32[i] := R1i[i,j] * χ⁽²⁾₁[j,k,h] * R1[k,a] * E1[a] * R1[h,b] * E2[b] verbose=1
# @tullio χ⁽²⁾₁2[i,j,k] := R1i[i,a] * χ⁽²⁾₁[a,b,c] * R1[b,j] * R1[c,k] verbose=1
# @tullio P33[i] := χ⁽²⁾₁2[i,j,k] * E1[j] * E2[k] verbose=1
# # @tullio P33[i] := R1[i,j] * χ⁽²⁾₁[j,k,h] * R1i[k,a] * E1[a] * R1i[h,b] * E2[b] verbose=1
# @assert Vector(P32) ≈ P3
# @assert Vector(P33) ≈ P3
# @assert χ⁽²⁾₁2 ≈ χ⁽²⁾₃
#
# import Rotations: Rotation
#
# function rotate(χ::AbstractArray{T,2},𝓡::TR) where {T<:Real, TR<:StaticMatrix{3,3}}
# 	@tullio χᵣ[i,j] := 𝓡[a,i] * 𝓡[b,j] * χ[a,b]  fastmath=true
# end
#
# function rotate(χ::AbstractArray{T,3},𝓡::TR) where {T<:Real, TR<:StaticMatrix{3,3}}
# 	@tullio χᵣ[i,j,k] := 𝓡[a,i] * 𝓡[b,j] * 𝓡[c,k] * χ[a,b,c]  fastmath=true
# end
#
# function rotate(χ::AbstractArray{T,4},𝓡::TR) where {T<:Real, TR<:StaticMatrix{3,3}}
# 	@tullio χᵣ[i,j,k,l] := 𝓡[a,i] * 𝓡[b,j] * 𝓡[c,k] * 𝓡[d,l] * χ[a,b,c,d]  fastmath=true
# end
# Zygote.gradient(x->dot(E1,rotate(eps1,Zygote.forwarddiff(y->MRP(RotY(y)),x)),E1),0.1)[1]
# ForwardDiff.derivative(x->dot(E1,rotate(eps1,MRP(RotY(x))),E1),0.1)
# Zygote.hessian(x->dot(E1,rotate(eps1,Zygote.forwarddiff(y->MRP(RotY(y)),x)),E1),0.1)[1]
#
# rotate(χ⁽²⁾₃,R1)
# @assert rotate(χ⁽²⁾₁,R1) ≈ MArray{Tuple{3,3,3}}(χ⁽²⁾₃)
# @assert rotate(χ⁽²⁾₃,R1') ≈ MArray{Tuple{3,3,3}}(χ⁽²⁾₁)
# @assert rotate(rotate(χ⁽²⁾₁,R1),R1i) ≈ MArray{Tuple{3,3,3}}(χ⁽²⁾₁)
# @assert rotate(rotate(χ⁽²⁾₁,R1),inv(R1)) ≈ MArray{Tuple{3,3,3}}(χ⁽²⁾₁)
# @assert rotate(χ⁽²⁾₃,R3) ≈ MArray{Tuple{3,3,3}}(χ⁽²⁾₁)
#
# round.(rotate(rotate(χ⁽²⁾₁,R1),inv(R1)),digits=3) ≈ MArray{Tuple{3,3,3}}(χ⁽²⁾₁)
#
# round.(rotate(rotate(χ⁽²⁾₁,R1),inv(R1)),digits=3) - MArray{Tuple{3,3,3}}(χ⁽²⁾₁)
#
#
#
#
#
# rotate(χ⁽²⁾₁,R3) - MArray{Tuple{3,3,3}}(χ⁽²⁾₃)
#
#
# rotate(rotate(χ⁽²⁾₁,R1),R1i) - MArray{Tuple{3,3,3}}(χ⁽²⁾₁)
#
# χ⁽²⁾₁
#
# AngleAxis(R1)
#
#
#
#
# transpose(MRP(R1))
#
#
#
#
# R1
#
#
# R3 = RotY(-π/2)
# R3 ≈ inv(R1)
# R1 * χ⁽²⁾₃
# R1
#
#
#
#
# inv(R1)
#
#
#
#
# R1 ≈ inv(inv(R1))
#
# R1
#
# inv(R1)


text:
function quicsearch(A, b)
    first_index = firstindex(A)
    last_index = lastindex(A)
    current_index = -1
    flag = false
    while (first_index + 1 < last_index) &&(current_index == -1)
        mid = div((first_index + last_index), 2)
        if A[mid] == b
            current_index = mid
            flag = true
        else
            if b < A[mid]
                last_index = mid -1
            else
                first_index = mid +1
            end
        end
    end
    if A[first_index] > b
        current_index = first_index
        flag = false
    elseif A[last_index] > b
        current_index = last_index
        flag = false
    elseif A[last_index] < b
        current_index = last_index + 1
        flag = false
    end
    return flag, current_index
end

function insertsort!(A)
    for i in 2:length(A)
        key = A[i]
        j = i
        c,d = quicsearch(A[1:i-1],key)
        while d != j
            A[j],A[j-1] = A[j-1],A[j]
            j -= 1
        end
    end
    return A
end