text:
using Test
using WignerSymbols
using LinearAlgebra
using Random

Random.seed!(1234)

smalljlist = 0:1//2:10
largejlist = 0:1//2:1000

@testset "triangle coefficient" begin
    for j1 in smalljlist, j2 in smalljlist
        for j3 = abs(j1-j2):(j1+j2)
            @test Δ(j1,j2,j3) ≈ sqrt(factorial(big(Int(j1+j2-j3)))*
                                factorial(big(Int(j1-j2+j3)))*factorial(big(Int(j2+j3-j1)))/
                                factorial(big(Int(j1+j2+j3+1))))
        end
    end
end

# test 3j:
#--------
@testset "clebschgordan: test orthogonality" begin
    for j1 in smalljlist, j2 in smalljlist
        d1::Int = 2*j1+1
        d2::Int = 2*j2+1
        M = zeros(Float64, (d1*d2, d1*d2))
        ind1 = 1
        for m1 in -j1:j1, m2 in -j2:j2
            ind2 = 1
            @inbounds for j3 in abs(j1-j2):(j1+j2), m3 in -j3:j3
                M[ind1,ind2] = clebschgordan(j1,m1,j2,m2,j3,m3)
                ind2 += 1
            end
            ind1 += 1
        end
        @test M'*M ≈ one(M)
    end
end

# test recurrence relations: Phys Rev E 57, 7274 (1998)
@testset "wigner3j: test recurrence relations" begin
    for k = 1:10
        j2 = convert(BigFloat, rand(0:1//2:1000))
        j3 = convert(BigFloat, rand(0:1//2:1000))
        m2 = convert(BigFloat, rand(-j2:j2))
        m3 = convert(BigFloat, rand(-j3:j3))

        for j in max(abs(j2-j3),abs(m2+m3))+1:(j2+j3)-1
            X = j*sqrt(((j+1)^2-(j2-j3)^2)*((j2+j3+1)^2-(j+1)^2)*((j+1)^2-(m2+m3)^2))
            Y = (2*j+1)*((m2+m3)*(j2*(j2+1)-j3*(j3+1)) - (m2-m3)*j*(j+1))
            Z = (j+1)*sqrt((j^2-(j2-j3)^2)*((j2+j3+1)^2-j^2)*(j^2-(m2+m3)^2))
            tol = 10*max(abs(X),abs(Y),abs(Z))*eps(BigFloat)
            @test (X*wigner3j(BigFloat,j+1,j2,j3,-m2-m3,m2,m3) +
                        Z*wigner3j(BigFloat,j-1,j2,j3,-m2-m3,m2,m3)) ≈
                            (-Y*wigner3j(BigFloat,j,j2,j3,-m2-m3,m2,m3)) atol=tol
        end
    end
end

@testset "wigner3j: test orthogonality relations" begin
    # equivalent to Clebsch-Gordan orthogonality, now test using Float32
    for j1 in smalljlist, j2 in smalljlist
        d1::Int = 2*j1+1
        d2::Int = 2*j2+1
        M = zeros(Float32, (d1*d2, d1*d2))
        ind2 = 1
        for m1 in -j1:j1, m2 in -j2:j2
            ind1 = 1
            @inbounds for j3 in abs(j1-j2):(j1+j2), m3 in -j3:j3
                d3::Int = 2*j3+1
                M[ind1,ind2] += sqrt(d3) * wigner3j(Float32, j1, j2, j3, m1, m2, m3)
                ind1 += 1
            end
            ind2 += 1
        end
        @test M'*M ≈ one(M) # orthogonality relation type 1
        @test M*M' ≈ one(M) # orthogonality relation type 2
    end
end

# test 6j
#----------
@testset "wigner6j: test orthogonality" begin
    for j1 in smalljlist, j2 in smalljlist, j4 in smalljlist
        for j5 in max(abs(j1-j2-j4),abs(j1-j2+j4),abs(j1+j2-j4)):(j1+j2+j4)
            j6range = max(abs(j2-j4),abs(j1-j5)):min((j2+j4),(j1+j5))
            j3range = max(abs(j1-j2),abs(j4-j5)):min((j1+j2),(j4+j5))
            @test length(j6range) == length(j3range)
            M = zeros(Float64, (length(j3range), length(j6range)))
            for (k2,j6) in enumerate(j6range)
                for (k1,j3) in enumerate(j3range)
                    M[k1,k2] = sqrt(2*j3+1)*sqrt(2*j6+1)*wigner6j(j1,j2,j3,j4,j5,j6)
                end
            end
            @test M'*M ≈ one(M)
        end
    end
end

@testset "wigner6j: test special cases" begin
    for j1 in smalljlist, j2 in smalljlist
        j6 = 0
        j4 = j2
        j5 = j1
        for j3 in abs(j1-j2):(j1+j2)
            @test wigner6j(j1,j2,j3,j4,j5,j6) ≈ (-1)^(j1+j2+j3)/sqrt((2*j1+1)*(2*j2+1))
        end
    end
end

@testset "wigner6j: test recurrence relation" begin
    for k = 1:10
        j2 = convert(BigFloat,rand(largejlist))
        j3 = convert(BigFloat,rand(largejlist))
        l1 = convert(BigFloat,rand(largejlist))
        l2 = convert(BigFloat,rand(abs(l1-j3):(l1+j3)))
        l3 = convert(BigFloat,rand(abs(l1-j2):min(l1+j2)))

        for j in intersect(abs(j2-j3):(j2+j3), abs(l2-l3):(l2+l3))
            X = j * sqrt( ((j+1)^2-(j2-j3)^2) * ((j2+j3+1)^2-(j+1)^2) *
                            ((j+1)^2-(l2-l3)^2) * ((l2+l3+1)^2 - (j+1)^2) )
            Y = (2*j+1) * ( j*(j+1)*( -j*(j+1) + j2*(j2+1) + j3*(j3+1) - 2*l1*(l1+1)) +
                            l2*(l2+1)*( j*(j+1) + j2*(j2+1) - j3*(j3+1) ) +
                            l3*(l3+1)*( j*(j+1) - j2*(j2+1) + j3*(j3+1) ) )
            Z = (j+1) * sqrt( (j^2-(j2-j3)^2) * ((j2+j3+1)^2-j^2) *
                                (j^2-(l2-l3)^2) * ((l2+l3+1)^2-j^2) )
            tol = 10 * max(abs(X), abs(Y), abs(Z)) * eps(BigFloat)
            @test (X*wigner6j(BigFloat, j+1, j2, j3, l1, l2, l3) +
                    Z*wigner6j(BigFloat, j-1, j2, j3, l1, l2, l3)) ≈
                        (-Y*wigner6j(BigFloat, j, j2, j3, l1, l2, l3)) atol=tol
        end
    end
end

@testset "test recoupling relation between 3j/clebschgordan and 6j/racahW symbols" begin
    smallerjlist = 0:1//2:5
    for j1 in smallerjlist, j2 in smallerjlist, j3 in smallerjlist
        m1range = -j1:j1
        m2range = -j2:j2
        m3range = -j3:j3
        V1 = Array{Float64}(undef, length(m1range),length(m2range),length(m3range))
        V2 = Array{Float64}(undef, length(m1range),length(m2range),length(m3range))
        for J in max(abs(j1-j2-j3),abs(j1-j2+j3),abs(j1+j2-j3)):(j1+j2+j3)
            J12range = max(abs(j1-j2),abs(J-j3)):min((j1+j2),(J+j3))
            J23range = max(abs(j2-j3),abs(j1-J)):min((j2+j3),(j1+J))
            for J12 in J12range, J23 in J23range
                M = rand(-J:J) # only test for one instance of M in -J:J, should be independent of M anyway
                fill!(V1,0)
                fill!(V2,0)
                for (k1,m1) in enumerate(m1range), (k2,m2) in enumerate(m2range)
                    abs(m1+m2)<=J12 || continue
                    for (k3,m3) in enumerate(m3range)
                        abs(m2+m3)<=J23 || continue
                        m1+m2+m3==M || continue
                        V1[k1,k2,k3] = clebschgordan(j1,m1,j2,m2,J12) *
                                        clebschgordan(J12,m1+m2,j3,m3,J)
                        V2[k1,k2,k3] = clebschgordan(j2,m2,j3,m3,J23) *
                                        clebschgordan(j1,m1,J23,m2+m3,J)
                    end
                end
                @test racahW(j1,j2,J,j3,J12,J23) ≈
                        dot(V2,V1)/sqrt((2*J12+1)*(2*J23+1)) atol=10*eps(Float64)
            end
        end
    end
end


text:
# Generates samples of dense random matrices that are distributed according
# to the classical Gaussian random matrix ensembles
#
# The notation follows closely that in:
#
# Ioana Dumitriu and Alan Edelman, " Models for Beta s"
# Journal of Mathematical Physics, vol. 43 no. 11 (2002), pp. 5830--5547
# doi: 10.1063/1.1507823
# arXiv: math-ph/0206043
#
# References
#
# Alan Edelman, Per-Olof Persson and Brian D Sutton, "The fourfold way"
# http://www-math.mit.edu/~edelman/homepage/papers/ffw.pdf

export GaussianHermite, GaussianLaguerre, GaussianJacobi,
       Wigner, Wishart, MANOVA #Synonyms for Hermite, Laguerre and Jacobi

#A convenience function to define a chi scalar random variable
chi(df::Real) = sqrt(rand(Chisq(df)))

#####################
# Hermite ensemble #
#####################

type GaussianHermite <: ContinuousMatrixDistribution
  beta::Real
end
typealias Wigner  GaussianHermite

# Generates a NxN symmetric Wigner matrix
function rand(d::GaussianHermite, n::Integer)
  if d.beta == 1 #real
    A = randn(n, n)
  elseif d.beta == 2 #complex
    A = randn(n, n) + im*randn(n, n)
  elseif d.beta == 4 #quaternion
    #Employs 2x2 matrix representation of quaternions
    X = randn(n, n) + im*randn(n, n)
    Y = randn(n, n) + im*randn(n, n)
    A = [X Y; -conj(Y) conj(X)]
  else
    error(@sprintf("beta = %d is not implemented", d.beta))
  end
  normalization = sqrt(2*d.beta*n)
  return (A + A') / normalization
end

rand(d::GaussianHermite, dims::Dim2) = dims[1]==dims[2] ? rand(d, dims[1]) : error("Can only generate square matrices")

#Generates a NxN tridiagonal Wigner matrix
#The beta=infinity case is defined in Edelman, Persson and Sutton, 2012
function tridrand(d::GaussianHermite, n::Integer)
  if d.beta<=0 error("beta must be positive") end
  if d.beta==Inf return SymTridiagonal(zeros(n), [sqrt(x/2) for x=n-1:-1:1]) end
  Hdiag = randn(n)/sqrt(n)
  Hsup = [chi(d.beta*i)/sqrt(2*n) for i=n-1:-1:1]
  return SymTridiagonal(Hdiag, Hsup)
end

tridrand(d::GaussianHermite, dims::Dim2) = dims[1]==dims[2] ? tridrand(d, dims[1]) : error("Can only generate square matrices")

#Return n eigenvalues distributed according to the Hermite ensemble
eigvalrand(d::GaussianHermite, n::Integer) = eigvals(tridrand(d, n))

#Calculate Vandermonde determinant term
function VandermondeDeterminant{Eigenvalue<:Number}(lambda::Vector{Eigenvalue}, beta::Real)
  n = length(lambda)
  Vandermonde = 1.0
  for j=1:n
    for i=1:j-1
      Vandermonde *= abs(lambda[i] - lambda[j])^beta
    end
  end
  Vandermonde
end

function eigvaljpdf{Eigenvalue<:Number}(d::GaussianHermite, lambda::Vector{Eigenvalue})
  n = length(lambda)
  #Calculate normalization constant
  c = (2pi)^(-n/2)
  for j=1:n
    c *= gamma(1 + beta/2)/gamma(1 + beta*j/2)
  end
  Energy = sum(lambda.^2/2) #Calculate argument of exponential
  VandermondeDeterminant(lambda, beta) * exp(-Energy)
end

#####################
# Laguerre ensemble #
#####################

type GaussianLaguerre <: ContinuousMatrixDistribution
  beta::Real
  a::Real
end
typealias Wishart GaussianLaguerre

#Generates a NxN Hermitian Wishart matrix
#TODO Check - the eigenvalue distribution looks funky
#TODO The appropriate matrix size should be calculated from a and one matrix dimension
function rand(d::GaussianLaguerre, dims::Dim2)
  n = 2.0*a/d.beta
  if d.beta == 1 #real
    A = randn(dims)
  elseif d.beta == 2 #complex
    A = randn(dims) + im*randn(dims)
  elseif d.beta == 4 #quaternion
    #Employs 2x2 matrix representation of quaternions
    X = randn(dims) + im*randn(dims)
    Y = randn(dims) + im*randn(dims)
    A = [X Y; -conj(Y) conj(X)]
    error(@sprintf("beta = %d is not implemented", d.beta))
  end
  return (A * A') / dims[1]
end

#Generates a NxN bidiagonal Wishart matrix
function bidrand(d::GaussianLaguerre, m::Integer)
  if d.a <= d.beta*(m-1)/2.0
    error(@sprintf("Given your choice of m and beta, a must be at least %f (You said a = %f)", d.beta*(m-1)/2.0, d.a))
  end
  Bidiagonal([chi(2*d.a-i*d.beta) for i=0:m-1], [chi(d.beta*i) for i=m-1:-1:1], true)
end

#Generates a NxN tridiagonal Wishart matrix
function tridrand(d::GaussianLaguerre, m::Integer)
  B = bidrand(d, m)
  L = B * B'
  SymTridiagonal(diag(L)/sqrt(m), diag(L,1)/sqrt(m))
end

#Return n eigenvalues distributed according to the Laguerre ensemble
eigvalrand(d::GaussianLaguerre, m::Integer) = eigvals(tridrand(d, m))

#TODO Check m and ns
function eigvaljpdf{Eigenvalue<:Number}(d::GaussianLaguerre, lambda::Vector{Eigenvalue})
  m = length(lambda)
  #Laguerre parameters
  p = 0.5*d.beta*(m-1) + 1.0
  #Calculate normalization constant
  c = 2.0^-(m*d.a)
  z = (d.a - d.beta*(m)*0.5)
  for j=1:m
     z += 0.5*d.beta
     if z < 0 && (int(z) - z) < eps()
         #Pole of gamma function, there is no density here no matter what
         return 0.0
     end
     c *= gamma(1 + beta/2)/(gamma(1 + beta*j/2)*gamma(z))
  end

  Prod = prod(lambda.^(a-p)) #Calculate Laguerre product term
  Energy = sum(lambda)/2 #Calculate argument of exponential
  return c * VandermondeDeterminant(lambda, beta) * Prod * exp(-Energy)
end

###################
# Jacobi ensemble #
###################

#Generates a NxN self-dual MANOVA matrix
type GaussianJacobi <: ContinuousMatrixDistribution
  beta::Real
  a::Real
  b::Real
end
typealias MANOVA  GaussianJacobi

function rand(d::GaussianJacobi, m::Integer)
  w1 = Wishart(m, int(2.0*d.a/d.beta), d.beta)
  w2 = Wishart(m, int(2.0*d.b/d.beta), d.beta)
  return (w1 + w2) \ w1
end

# A helper function for Jacobi samples
function SampleCSValues(n::Integer, a::Real, b::Real, beta::Real)
  if beta == Inf
    c =sqrt((a+[1:n])./(a+b+2*[1:n]))
    s =sqrt((b+[1:n])./(a+b+2*[1:n]))
    cp=sqrt([1:n-1]./(a+b+1+2*[1:n-1]))
    sp=sqrt((a+b+1+[1:n-1])./(a+b+1+2*[1:n-1]))
  else
    #Generate cosine-squared values
    csq  = [rand(Beta(beta*(a+i)/2,beta*(b+i)/2)) for i=1:n]
    cpsq = [rand(Beta(beta*i/2,beta*(a+b+1+i)/2)) for i=1:n]
    #Cosine-sine pairs
    c , s  = sqrt(csq) , sqrt(1-csq)
    cp, sp = sqrt(cpsq), sqrt(1-cpsq)
  end
  return c, s, cp, sp
end

#Generates a 2Mx2M sparse MANOVA matrix
#Jacobi ensemble
#
# Reference:
#   A. Edelman and B. D. Sutton, "The beta-Jacobi matrix model, the CS decomposition,
#   and generalized singular value problems", Foundations of Computational Mathematics,
#   vol. 8 iss. 2 (2008), pp 259-285.
#TODO check normalization
function sprand(d::GaussianJacobi, n::Integer, a::Real, b::Real)
  CoordI = zeros(8n-4)
  CoordJ = zeros(8n-4)
  Values = zeros(8n-4)

  c, s, cp, sp = SampleCSValues(n, a, b, d.beta)

  #Diagonals of each block
  for i=1:n
    CoordI[i], CoordJ[i] = i, i
    Values[i] = i==1 ? c[n] : c[n+1-i] * sp[n+1-i]
  end
  for i=1:n
    CoordI[n+i], CoordJ[n+i] = i, n+i
    Values[n+i] = i==n ? s[1] : s[n+1-i] * sp[n-i]
  end
  for i=1:n
    CoordI[2n+i], CoordJ[2n+i] = n+i, i
    Values[2n+i] = i==1 ? -s[n] : -s[n+1-i] * sp[n+1-i]
  end
  for i=1:n
    CoordI[3n+i], CoordJ[3n+i] = n+i, n+i
    Values[3n+i] = i==n ? c[1] : c[n+1-i] * sp[n-i]
  end
  #Off-diagonals of each block
  for i=1:n+1
    CoordI[4n+i], CoordJ[4n+i] = i,i+1
    Values[4n+i] = -s[n+1-i]*cp[n-i]
  end
  for i=1:n+1
    CoordI[5n-1+i], CoordJ[5n-1+i] = i+1,n+i
    Values[5n-1+i] = c[n-i]*cp[n-i]
  end
  for i=1:n+1
    CoordI[6n-2+i], CoordJ[6n-2+i] = n+i,i+1
    Values[6n-2+i] = -c[n+1-i]*cp[n-i]
  end
  for i=1:n+1
    CoordI[7n-3+i], CoordJ[7n-3+i] = n+i,i+1
    Values[7n-3+i] = -s[n-i]*cp[n-i]
  end

  return sparse(CoordI, CoordJ, Values)
end

#Return n eigenvalues distributed according to the Jacobi ensemble
function eigvalrand(d::GaussianJacobi, n::Integer)
  #Generate just the upper left quadrant of the matrix
  c, s, cp, sp = SampleCSValues(n, d.a, d.b, d.beta)
  dv = [i==1 ? c[n] : c[n+1-i] * sp[n+1-i] for i=1:n]
  ev = [-s[n+1-i]*cp[n-i] for i=1:n-1]

  ##TODO: understand why dv and ev are returned as Array{Any,1}
  M = Bidiagonal(convert(Array{Float64,1},dv), convert(Array{Float64,1},ev), false)
  return svdvals(M)
end

#TODO Check m and ns
function eigvaljpdf{Eigenvalue<:Number}(d::GaussianJacobi, lambda::Vector{Eigenvalue})
  m = length(lambda)
  #Jacobi parameters
  a1, a2 = d.a, d.b
  p = 1.0 + d.beta*(m-1)/2.0
  #Calculate normalization constant
  c = 1.0
  for j=1:m
    z1 = (a1 - beta*(m-j)/2.0)
    if z1 < 0 && (int(z1) - z1) < eps()
      return 0.0 #Pole of gamma function, there is no density here no matter what
    end
    z2 = (a2 - beta*(m-j)/2.0)
    if z2 < 0 && (int(z2) - z2) < eps()
      return 0.0 #Pole of gamma function, there is no density here no matter what
    end
    c *= gamma(1 + beta/2)*gamma(a1+a2-beta*(m-j)/2)
    c /= gamma(1 + beta*j/2)*gamma(z1)*gamma(z2)
  end

  Prod = prod(lambda.^(a1-p))*prod((1-lambda).^(a2-p)) #Calculate Laguerre product term
  Energy = sum(lambda/2) #Calculate argument of exponential

  return c * VandermondeDeterminant(lambda, beta) * Prod * exp(-Energy)
end

ext: jl


text:
module Day12

# Part 1
instructions = readlines("day12_input.txt")

mutable struct Ship
    pos
    angle
end

function move!(ship::Ship, instruction)
    rules = Dict(
        'N' => (s::Ship, v) -> s.pos += [0, v],
        'S' => (s::Ship, v) -> s.pos -= [0, v],
        'E' => (s::Ship, v) -> s.pos += [v, 0],
        'W' => (s::Ship, v) -> s.pos -= [v, 0],
        'L' => (s::Ship, v) -> s.angle += v,
        'R' => (s::Ship, v) -> s.angle -= v,
        'F' => (s::Ship, v) -> s.pos += v * [round(Int, cosd(s.angle)), round(Int, sind(s.angle))]
    )
    rules[first(instruction)](ship, parse(Int32, instruction[2:end]))
end

ship = Ship([0, 0], 0)
map(i -> move!(ship, i), instructions)
println("Part 1: Manhattan dist = $(sum(abs.(ship.pos)))")

# Part 2
mutable struct Waypoint
    pos
end

function rotationmatrix(angle)
    c = round(Int, cosd(angle))
    s = round(Int, sind(angle))
    [c -s
     s c]
end

function move!(ship::Ship, waypoint::Waypoint, instruction)
    rules = Dict(
        'N' => (s::Ship, w::Waypoint, v) -> w.pos += [0, v],
        'S' => (s::Ship, w::Waypoint, v) -> w.pos -= [0, v],
        'E' => (s::Ship, w::Waypoint, v) -> w.pos += [v, 0],
        'W' => (s::Ship, w::Waypoint, v) -> w.pos -= [v, 0],
        'L' => (s::Ship, w::Waypoint, v) -> w.pos = rotationmatrix(v) * w.pos,
        'R' => (s::Ship, w::Waypoint, v) -> w.pos = rotationmatrix(-v) * w.pos,
        'F' => (s::Ship, w::Waypoint, v) -> s.pos += v * w.pos
    )
    rules[first(instruction)](ship, waypoint, parse(Int32, instruction[2:end]))
end

ship = Ship([0, 0], 0)
waypoint = Waypoint([10, 1])
map(i -> move!(ship, waypoint, i), instructions)
println("Part 2: Manhattan dist = $(sum(abs.(ship.pos)))")

end

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
#!/usr/bin/env julia

include("rsa.jl")

#=
Transform a string into a BigInt, add in 0xAA to avoid unpleasantness. Proper PKCS
padding will have to wait for now. We treat it as a base-256 integer and just pull
off the digits.
=#

function str2bigint(s)
    sum::BigInt = 0
    pow::BigInt = 1
    for c in s
        sum += pow * (0xAA ⊻ BigInt(c))
        pow *= 256
    end
    sum
end

#=
Transform a BigInt back into a string, subtracting off the 0xAA. We treat it as a base-256
integer and just pull off the digits.
=#

function bigint2str(n)
    s = ""
    while n > 0
        s = s * Char(0xAA ⊻ (n % 256))
        n ÷= 256
    end
    s
end

print("How many bits? ")

bits = parse(Int64, readline())

(e, d, n) = makeKey(bits);

println("e = $e")
println("d = $d")
println("n = $n")

print(">> ")
for m in eachline()
    c = encrypt(str2bigint(m), e, n); println("En[$m] = $c")
    t = bigint2str(decrypt(c, d, n)); println("De[$c] = $t")
    print(">> ")
end

max_stars_repo_name: darrelllong/RSA-in-Julia

max_issues_repo_name: darrelllong/RSA-in-Julia

max_forks_repo_name: darrelllong/RSA-in-Julia


text:
function getTreasuries()
    (intRates, metaData) = EconDatasets.readUSTreasuries()

    ## save on disk
    fileName = joinpath(Pkg.dir("EconDatasets"), "data", "Treasuries.csv")
    writetimearray(intRates, fileName)
end


text:
function check_asy_present()
    global _HAVE_ASY
    if !_HAVE_ASY
        error("3D Output not available with Cairo backend")
    end
end

extension(f) = split(f,'.')[end]
sans_extension(f) = join(split(f,'.')[1:end-1],'.')

"""
    save(filename::AbstractString,
         P::Union{Plot2D,Plot3D};
         runasy=true,
         forcepdf=false)

Save Asymptote figure. If `filename` has extension `.asy`, then
an asy file is saved together with any auxiliary data files.

If `filename` has extension `.pdf`, `.svg` or `.png`, then
only the resulting image file is saved to the location `filename`
"""
function save(filename::AbstractString,
              P::Union{Plot2D,Plot3D};
              runasy=true,
              forcepdf=false,
              bbox=false)
    global _BACKEND
    if !any(endswith(filename,ext) for ext in (".asy",".pdf",".svg",".png"))
        error("Extension should be asy, pdf, svg, or png")
    end
    if _BACKEND == "asy"
        if endswith(filename,".asy")
            check_asy_present()
            if forcepdf
                Q = deepcopy(P)
                filter!(p -> first(p) ≠ :pdf, Q.options)
                push!(Q.options,(:pdf,true))
            else
                Q = P
            end
            if bbox ≠ false
                Q += Plot(box(bbox.xmin,bb.max,bb.ymin,bb.ymax;other=invisible))
            end
            S = AsyString(Q)
            write(filename,S.str)
            directory = dirname(filename)
            if directory == ""
                directory = "."
            end
            for (k,v) in S.data
                DelimitedFiles.writedlm("$directory/$k",v,',')
            end
            if runasy
                originaldir = pwd()
                try
                    cd(directory)
                    run(`asy $(basename(filename))`)
                finally
                    cd(originaldir)
                end
            end
        else
            tempdir = mktempdir()
            asyfile = "$tempdir/myplot.asy"
            pdffile = "$tempdir/myplot.pdf"
            runasy = !endswith(filename,".asy")
            save(asyfile,P;runasy=runasy,forcepdf=true)
            if endswith(filename,".asy")
                mv(asyfile,filename,force=true)
            elseif endswith(filename,".pdf")
                mv(pdffile,filename,force=true)
            elseif endswith(filename,".png")
                global _DEFAULT_WIDTH
                D = Dict{Symbol,Any}(P.options)
                width = :width in keys(D) ? D[:width] : _DEFAULT_WIDTH
                base = sans_extension(filename)
                run(`pdftoppm -png -scale-to-x $width -scale-to-y
                                        -1 -singlefile $pdffile $base`)
            elseif endswith(filename,".svg")
                run(`pdf2svg $tempdir/myplot.pdf $filename`)
            end
        end
    else
        if isa(P,Plot3D)
            error("3D output not supported with Cairo backend")
        end
        D = Dict{Symbol,Any}(P.options)
        border = (:border in keys(D) ? D[:border]
                        : _DEFAULT_PLOT2D_KWARGS[:border])
        bufferdata = bytes(P;format=Symbol(extension(filename)),
                             bbox=bbox,
                             border=border)
        write(filename,bufferdata)
    end
end

save(S::GraphicElement;kwargs...) = save(Plot([S],kwargs))

function save(P::Union{Plot2D,Plot3D},
              filename::AbstractString; kwargs...)
    @warn "file name should be the first argument to save"
    save(filename, P; kwargs...)
end

_SHOW_PLOTS = true
function showplots(b::Bool)
    global _SHOW_PLOTS
    _SHOW_PLOTS = b
    """Plot display turned $(b ? "on" : "off")"""
end

Base.showable(::MIME"image/svg+xml",P::Plot2D) = true
Base.showable(::MIME"image/svg+xml",P::Plot3D) = false
Base.showable(::MIME"image/png",P::Plot2D) = true

function png()
    @eval Base.showable(::MIME"image/svg+xml",P::Plot2D) = false
    @eval Base.showable(::MIME"image/svg+xml",P::Plot3D) = false
    "Turning svg display off"
end

function svg()
    @eval Base.showable(::MIME"image/svg+xml",P::Plot2D) = true
    @eval Base.showable(::MIME"image/svg+xml",P::Plot3D) = true
    "Turning svg display on"
end

function Base.show(io::IO, ::MIME"image/svg+xml", P::Plot2D)
    if length(P.elements) == 0
        return nothing
    end
    global _BACKEND
    if _BACKEND == "asy"
        directory = mktempdir()
        filename = "$directory/myplot.svg"
        save(filename, P)
        write(io, read("$directory/myplot.svg"))
    else
        D = Dict{Symbol,Any}(P.options)
        border = (:border in keys(D) ? D[:border]
                        : _DEFAULT_PLOT2D_KWARGS[:border])
        write(io,bytes(P;format=:svg,border=border))
    end
end

function Base.show(io::IO, ::MIME"image/png", P::Plot2D)
    if length(P.elements) == 0
        return nothing
    end
    global _BACKEND
    if _BACKEND == "asy"
        directory = mktempdir()
        filename = "$directory/myplot.png"
        save(filename, P)
        write(io, read("$directory/myplot.png"))
    elseif _BACKEND == "cairo"
        D = Dict{Symbol,Any}(P.options)
        border = (:border in keys(D) ? D[:border]
                        : _DEFAULT_PLOT2D_KWARGS[:border])
        write(io,bytes(P;format=:png,border=border))
    end
end

function Base.show(io::IO, ::MIME"image/svg+xml", P::Plot3D)
    if length(P.elements) == 0
        return
    end
    global _HAVE_ASY
    if !_HAVE_ASY
        l = length(P.elements)
        show(io,"""Plot3D(<$(l)element$(l == 1 ? "" : "s")>)""")
        return
    end
    tempdir = mktempdir()
    filename = "$tempdir/myplot.svg"
    save(filename, P)
    write(io,read(filename))
    write("/Users/sswatson/Desktop/otherstuff.svg",read(filename))
end

function Base.show(io::IO, ::MIME"image/png", P::Plot3D)
    if length(P.elements) == 0
        return nothing
    end
    global _HAVE_ASY
    if !_HAVE_ASY
        l = length(P.elements)
        show(io,"""Plot3D(<$(l)element$(l == 1 ? "" : "s")>)""")
        return
    end
    tempdir = mktempdir()
    filename = "$tempdir/myplot.png"
    save(filename, P)
    write(io, read("$tempdir/myplot.png"))
end

_is_ijulia() = isdefined(Main, :IJulia) && Main.IJulia.inited
_is_juno()   = isdefined(Main, :Juno) && Main.Juno.isactive()

function Base.show(io::IO,::MIME"text/plain",P::Plot)
    if length(P.elements) == 0
        return nothing
    end
    global _SHOW_PLOTS
    if _SHOW_PLOTS && !_is_ijulia() && !_is_juno() && isdefined(Base, :active_repl)
        tempdir = mktempdir()
        filename = "$tempdir/myplot.pdf"
        save(filename,P)
        try
            if Sys.isapple()
                run(`open $filename`)
            elseif Sys.islinux() || is_bsd()
                run(`xdg-open $filename`)
            elseif Sys.iswindows()
                run(`start $filename`)
            end
        catch e
            error(string(
            "Failed to show the generated pdf.",
            "Run `showplots(false)` to stop trying to show pdfs.\n",
                "Error: "), sprint(Base.showerror, e))
        end
    else
        n = length(P.elements)
        s = n == 1 ? "" : "s"
        print(io,"$(split(string(typeof(P)),".")[end])(<$n elements>)")
    end
end


text:
for OP in (:(Base.eigvals),:(Base.eigs))
    @eval $OP(A::Operator,n::Integer;tolerance::Float64=100eps())=$OP([A],n;tolerance=tolerance)
end

Base.eigvals{T<:Operator}(A::Vector{T},n::Integer;tolerance::Float64=100eps())=eigs(A,n;tolerance=tolerance)[1]

function Base.eigs{T<:Operator}(A::Vector{T},n::Integer;tolerance::Float64=100eps())
    nf = length(A)-1
    for k=1:nf
        @assert isafunctional(A[k])
    end
    A = promotedomainspace(A,choosedomainspace(A))
    typ = mapreduce(eltype,promote_type,A)

    ds=domainspace(A[end])
    C = Conversion(ds,rangespace(A[end]))

    A1,C1 = zeros(typ,n,n),zeros(typ,n,n)
    for k=1:nf
        A1[k,1:n] = A[k][1:n]
    end
    A1[1+nf:end,1:n] = A[end][1:n-nf,1:n]
    C1[1+nf:end,1:n] = C[1:n-nf,1:n]

    λ,V=eig(A1,C1)

    pruneeigs(λ,V,ds,tolerance)
end

function pruneeigs(λ,V,ds,tolerance)
    retλ=eltype(λ)[]
    retV=Fun{typeof(ds),eltype(V)}[]
    n=length(λ)
    for k=1:n
        if slnorm(V,n-3:n,k)≤tolerance
            push!(retλ,λ[k])
            push!(retV,Fun(V[:,k],ds))
        end
    end
    retλ,retV
end


text:
using Documenter, Inpaintings

makedocs(
    sitename="Inpaintings Documentation",
    # options
    modules = [Inpaintings]
)

deploydocs(
    repo = "github.com/briochemc/Inpaintings.jl.git",
)


text:
mutable struct Die
    face::Vector{String}
end

function roll!(die, dir)
    face = die.face
    if dir == 'E'
        face[1], face[3], face[4], face[6] = face[4], face[1], face[6], face[3]
    elseif dir == 'W'
        face[1], face[3], face[4], face[6] = face[3], face[6], face[1], face[4]
    elseif dir == 'N'
        face[1], face[2], face[5], face[6] = face[2], face[6], face[1], face[5]
    else
        face[1], face[2], face[5], face[6] = face[5], face[1], face[6], face[2]
    end
end

die = Die(split(readline()))
dirs = readline()
for c in dirs
    roll!(die, c)
end
println(die.face[1])


text:
using LinearAlgebra, Polynomials
@testset "Polynomial evaluation" begin
    A = Polynomial("x")

    coeff = collect(1:21.0)

    poly_gens = Dict(
        "Monomial" => (:graph_monomial, :graph_monomial_degopt),
        "Horner"   => (:graph_horner, :graph_horner_degopt),
        "PS"       => (:graph_ps, :graph_ps_degopt),
    )

    for key in keys(poly_gens)
        @testset "$key" begin
            for n = 2:length(coeff)
                (graph_1, cref_1) = eval(poly_gens[key][1])(coeff[1:n])
                (graph_2, cref_2) = eval(poly_gens[key][2])(coeff[1:n])
                @test eval_graph(Compgraph(Any, graph_2), A) ==
                      eval_graph(Compgraph(Any, graph_1), A)
                @test sum(values(graph_2.operations) .== :mult) ==
                      sum(values(graph_1.operations) .== :mult)
            end

            if key == "Horner" # Add test here if scaling is available
                s = -0.89
                (graph_1, cref_1) =
                    eval(poly_gens[key][2])(coeff .* s .^ (0:length(coeff)-1))
                (graph_2, cref_2) = eval(poly_gens[key][2])(coeff, scaling = s)
                @test eval_graph(graph_1, A) ≈ eval_graph(graph_2, A)
            end
        end
    end
end


text:
__precompile__()
module eps_InAsSbP
export G_E1_E2,func_Q,func_Q2,E_T,func_f,InAsSbP_struct, eps_InAsSbP_xy,eps_InAsSbP_imag_xy
#Used this version in the end, accurate imaginary part from the range of 0.1 to 1.5, above 1.5, imaginary part becomes negative"""
#import numpy as np
#import matplotlib.pyplot as plt

e= 1.602*10^(-19)
eV = 1.60218*10^(-19)
eps_0 = 8.8542*10^(-12)
c = 3.0*10^8
hbar = 1.05457*10^(-34) #m2kg/s
hbar_eV = 6.5821*10^(-16) #eV s
m_e =9.109*10^(-31)#kg
m_0 = 0.9109*10^(-30)
m_0_eV = 0.51099906*10^6 #eV
T = 300.0

#list of constants: list of constants for InAsSbP, AD=InAs, BD=InSb, CD=InP ABC=AsSbP  #Cuevas Paper Table1
#can't find bowing constants for AsSbP
#Bowing Constants for GaIn, AsSb from ZhangThermoConversion
A_InAs = 19.950   #A (eV^1.5) for eps_1, Cuevas
A_InSb = 16.830   #A (eV^1.5) for eps_1, Cueva #A_InP =  6.57     #Need fitting Adachi Table1
C_InAs = 1.750  #C(Dimension-less) for eps_3 Cuevas
C_InSb = 2.64#C_InP = 1.49 #Need fitting Adachi Table1
D_InAs =21.511 #D(Dimension-less) for eps_4 Cuevas
D_InSb = 38.830 #D_InP = 60.4 #Need fitting Adachi Table1
G_E0_InAs = 8.660 #in eV Cuevas
G_E0_InSb = 4.20

@inline function  G_E1_E2(gamma_L,g,T)
    return gamma_L + g*T
end
precompile(G_E1_E2,(Float64,Float64,Float64))
G_E1_InAs = G_E1_E2(0.25,0.3,T)*0.001 #gamma in meV, change to eV #Cuevas Table 3
G_E1_InSb = G_E1_E2(0.15,0.28,T)*0.001
G_E1_InP  = 0.094 # Adachi Table1

G_E2_InAs =  G_E1_E2(403,0.3,T)*0.001 #gamma in meV, change to eV #Cuevas Table 3
G_E2_InSb =  G_E1_E2(633.9,0.49,T)*0.001
G_E2_InP =  0.10 # Adachi Table1

G_Eg_InAs =  0.350#(eV) Cuevas
G_Eg_InSb = 0.56
#Cuevas Table2, The energies are in eV, NO bowing constants for E0??Only temperature depence?
E0_InAs = 0.42  #at 300K, E0_T = 1.5,at 0K 1.4
dt0_InAs = 0.28*0.001
b0_InAs = 93.0
E0_InSb = 0.24
dt0_InSb = 0.32*0.001
b0_InSb = 170.0 #E0_InP  = 1.35 #Need fitting Adachi Table1

E0Delta0_InAs = 0.79#-0.42
dtE0Delta0_InAs = 0.34*0.001
bE0Delta0_InAs  = 248.0
E0Delta0_InSb = 1.20#-0.24
dtE0Delta0_InSb = 0.32*0.001
bE0Delta0_InSb = 170.0 #E0Delta0_GaAs = 1.85#-0.152 dtE0Delta0_GaAs = 0.35*0.001 bE0Delta0_GaAs = 225.0 E0Delta0_GaSb = 1.57#-0.81 dtE0Delta0_GaSb = 0.42*0.001
#E0Delta0_InP = 1.35 + 0.10 # Adachi Table1

E1_InAs = 2.5 #Cuevas Table2
dt1_InAs = 0.5*0.001
b1_InAs = 0.0
E1_InSb = 2.0
dt1_InSb = 0.68*0.001
b1_InSb = 132.0
#E1_InP = 3.10 # Adachi Table1

E1Delta1_InAs =2.88#-2.61 #Cuevas Table 2
dtE1Delta1_InAs = 0.5*0.001
bE1Delta1_InAs = 0.0
E1Delta1_InSb = 2.49#-2.00
dtE1Delta1_InSb = 0.65*100
bE1Delta1_InSb = 170.0 #E1Delta1_InP = 3.10+0.15 #Adachi  Table1

Delta1_InAs = 2.88-2.61 #Cuevas
Delta1_InSb = 2.49-2.00 #Delta1_InP = 0.15 #Adachi  Table1

E2_InAs = 4.70 #from Zhang Thermo Conversion
dt2_InAs = 0.56*0.001
b2_InAs = 0.0
E2_InSb = 3.9
dt2_InSb = 0.54*0.001
b2_InSb = 0.0
E2_InP = 4.7 #Adachi  Table1

Eg_InAs = 1.13 #Cuevas Table 2
dtg_InAs = 0.28*0.001
bg_InAs = 93.0
Eg_InSb = 0.93
dtg_InSb = 0.32*0.001
bg_InSb = 170.0
Eg_InP  = 2.05 #Adachi Table1

#Cuevas paper Table IV in angstrong
alpha_lc_InAs = 6.06
alpha_lc_InSb = 6.48
alpha_lc_InP =  5.87
#alpha_lc_GaAs = 5.65 alpha_lc_GaSb = 6.10 alpha_lc_GaIn = 0.055  #Cuevas no value #zhang 0.055 alpha_lc_AsSb = -0.023 #Cuevas no value #zhang -0.023

#Fitted Parameters
A_InP =  31.199833249457633
C_InP=  0.5600889435237922
D_InP=  -147.6161185407525
E0_InP = 2.106123911548122
E0Delta0_InP = 1.422040863793526
E1_InP = 3.030856569453735
E1Delta1_InP = 2.88535017475427
Delta1_InP =  -0.3610122449508903
G_E0_InP= 0.12337661235349834
G_Eg_InP=  0.8902957174564856

#Q Parameter in Notes
@inline function  func_Q(x,y,B_AD,B_BD,B_AC,B_CD,C_ABC)
    return x*B_AD + y*B_BD + (1-x-y)*B_CD  + C_ABC*x*y*(1-x-y)
end
precompile(func_Q,(Float64,Float64,Float64,Float64,Float64,Float64,Float64))

 #Q Parameter in Notes, No bowing parameters,  AD=InAs, BD=InSb, CD=InP ABC=AsSbP
@inline function  func_Q2(x,y,B_AD,B_BD,B_CD)
    return x*B_AD + y*B_BD + (1-x-y)*B_CD
end
precompile(func_Q2,(Float64,Float64,Float64,Float64,Float64))

@inline function  E_T(E0,dt,beta)
    return E0-dt*300.0^2/(300.0+beta)
end
precompile(E_T,(Float64,Float64,Float64))

@inline function  func_f(z)
    return z^(-2)*(2.0-(1.0+z)^0.5-(1.0-z)^0.5)
end
precompile(func_f,(Float64,))

struct InAsSbPDsc
    x::Float64
    y::Float64
    A_InAsSbP::Float64
    C_InAsSbP::Float64
    D_InAsSbP::Float64
    G_E0_InAsSbP::Float64 #No value found for Gamma_E1
    G_E1_InAsSbP::Float64
    G_Eg_InAsSbP::Float64
    G_E2_InAsSbP::Float64
    E0_InAsSbP::Float64
    E0Delta0_InAsSbP::Float64
    E1_InAsSbP::Float64
    E1Delta1_InAsSbP::Float64
    Delta1_InAsSbP::Float64
    E2_InAsSbP::Float64
    Eg_InAsSbP::Float64
    alpha_lc_InAsSbP::Float64
end

@inline function InAsSbP_struct(x,y)
     #assume equal proportion of As,Sb and P
    #list of constants for InAsSbP, AD=InAs, BD=InSb, CD=InP ABC=AsSbP  #Cuevas Paper Table1
    #Bowing Constants for GaIn, AsSb from ZhangThermoConversion
	E0_T_InAs = E_T(E0_InAs,dt0_InAs,b0_InAs)
	E0_T_InSb = E_T(E0_InSb,dt0_InSb,b0_InSb)
	E0Delta0_T_InAs = E_T(E0Delta0_InAs,dtE0Delta0_InAs,bE0Delta0_InAs)
	E0Delta0_T_InSb = E_T(E0Delta0_InSb,dtE0Delta0_InSb,bE0Delta0_InSb)
	E1_T_InAs = E_T(E1_InAs,dt1_InAs,b1_InAs)
	E1_T_InSb = E_T(E1_InSb,dt1_InSb,b1_InSb)
	E1Delta1_T_InAs = E_T(E1Delta1_InAs,dtE1Delta1_InAs,bE1Delta1_InAs)
	E1Delta1_T_InSb =E_T(E1Delta1_InSb,dtE1Delta1_InSb,bE1Delta1_InSb)
	E2_T_InAs = E_T(E2_InAs,dt2_InAs,b2_InAs)
	E2_T_InSb = E_T(E2_InSb,dt2_InSb,b2_InSb)
	Eg_T_InAs = E_T(Eg_InAs,dtg_InAs,bg_InAs)
	Eg_T_InSb = E_T(Eg_InSb,dtg_InSb,bg_InSb)

	A = func_Q2(x,y,A_InAs,A_InSb,A_InP)
    C= func_Q2(x,y,C_InAs,C_InSb,C_InP)
    D = func_Q2(x,y,D_InAs,D_InSb,D_InP)
    G_E0 = func_Q2(x,y,G_E0_InAs,G_E0_InSb,G_E0_InP) #No value found for Gamma_E1
    G_E1 = func_Q2(x,y,G_E1_InAs,G_E1_InSb,G_E1_InP)
    G_Eg = func_Q2(x,y,G_Eg_InAs,G_Eg_InSb,G_Eg_InP)
    G_E2 = func_Q2(x,y,G_E2_InAs,G_E2_InSb,G_E2_InP)
    E0 = func_Q2(x,y,E0_T_InAs,E0_T_InSb,E0_InP)    #No dt0_InP
    E0Delta0 = func_Q2(x,y,E0Delta0_T_InAs, E0Delta0_T_InSb,E0Delta0_InP) #no E0Delta0_T_InP
    E1 = func_Q2(x,y,E1_T_InAs,E1_T_InSb,E1_InP) #no E1_T_InP
    E1Delta1 = func_Q2(x,y,E1Delta1_T_InAs,E1Delta1_T_InSb,E1Delta1_InP)
    Delta1= func_Q2(x,y,Delta1_InAs,Delta1_InSb,Delta1_InP)
    E2= func_Q2(x,y,E2_T_InAs,E2_T_InSb,E2_InP)
    Eg= func_Q2(x,y,Eg_T_InAs,Eg_T_InSb,Eg_InP)
    alpha_lc = func_Q2(x,y,alpha_lc_InAs,alpha_lc_InSb,alpha_lc_InP)
    return InAsSbPDsc(x,y,A,C,D,G_E0,G_E1,G_Eg,G_E2,E0,E0Delta0,E1,E1Delta1, Delta1, E2,Eg ,alpha_lc)
end
precompile(InAsSbP_struct,(Float64,Float64))

@inline function  eps_InAsSbP_xy(E_photon,InAsSbPstruct)
    #GaInAsSb in the energy range 0.5-6eV omega  (hbar*2*math.pi)*c/(0.5*eV) = 2.5um, frequency
    # m_GaInAsSb = 0.022+0.032*x -0.012*x^2  #Equation 4a) of Zhang thermal Conversion, close to m_0
    m_InAsSbP = 0.028     #checked online, cannot find
    eps_1 =InAsSbPstruct.A_InAsSbP*InAsSbPstruct.E0_InAsSbP^(-1.5)*(func_f((E_photon + im*InAsSbPstruct.G_E0_InAsSbP)/InAsSbPstruct.E0_InAsSbP) + 0.5*(InAsSbPstruct.E0_InAsSbP/(InAsSbPstruct.E0Delta0_InAsSbP))^1.5*func_f((E_photon + im*InAsSbPstruct.G_E0_InAsSbP)/InAsSbPstruct.E0Delta0_InAsSbP))
    B1=44*(InAsSbPstruct.E1_InAsSbP+InAsSbPstruct.Delta1_InAsSbP/3)/(InAsSbPstruct.alpha_lc_InAsSbP*InAsSbPstruct.E1_InAsSbP^2)
    B2=44*(InAsSbPstruct.E1_InAsSbP+2*InAsSbPstruct.Delta1_InAsSbP/3)/(InAsSbPstruct.alpha_lc_InAsSbP*(InAsSbPstruct.E1_InAsSbP+InAsSbPstruct.Delta1_InAsSbP)^2)
    eps_2term1 = -B1*((E_photon + im*InAsSbPstruct.G_E1_InAsSbP)/InAsSbPstruct.E1_InAsSbP)^(-2)*log(1-((E_photon + im*InAsSbPstruct.G_E1_InAsSbP)/InAsSbPstruct.E1_InAsSbP)^2)
    eps_2term2 = -B2*((E_photon + im*InAsSbPstruct.G_E1_InAsSbP)/(InAsSbPstruct.E1Delta1_InAsSbP))^(-2)*log(1-((E_photon + im*InAsSbPstruct.G_E1_InAsSbP)/(InAsSbPstruct.E1Delta1_InAsSbP))^2)
    eps_2 = eps_2term1 + eps_2term2
    eps_3 = InAsSbPstruct.C_InAsSbP*InAsSbPstruct.E2_InAsSbP^2/(InAsSbPstruct.E2_InAsSbP^2-E_photon^2 + im*(-E_photon*InAsSbPstruct.G_E2_InAsSbP))
    #changed from G_E0_GaInAsSb to G_Eg_GaInAsSb
    E_complex = (E_photon + im*InAsSbPstruct.G_Eg_InAsSbP)
    eps_4_1st =-log(InAsSbPstruct.E1_InAsSbP/InAsSbPstruct.Eg_InAsSbP)*InAsSbPstruct.Eg_InAsSbP^2/E_complex^2
    eps_4_2nd = (1+InAsSbPstruct.Eg_InAsSbP/E_complex)^2*log((E_complex+InAsSbPstruct.E1_InAsSbP)/(E_complex+InAsSbPstruct.Eg_InAsSbP))/2
    eps_4_3rd = (1-InAsSbPstruct.Eg_InAsSbP/E_complex)^2*log((E_complex-InAsSbPstruct.E1_InAsSbP)/(E_complex-InAsSbPstruct.Eg_InAsSbP))/2
    eps_4 = (2*InAsSbPstruct.D_InAsSbP/pi)*(eps_4_1st+eps_4_2nd+eps_4_3rd)
    eps = eps_1 + eps_2 + eps_3 + eps_4
    return eps
end
precompile(eps_InAsSbP_xy,(Float64,InAsSbPDsc))

@inline function eps_InAsSbP_imag_xy(enr,InAsSbPstruct)
    #In the code InAs_xSb_yP_1-x-y,  from III-V Ternary and Quaternary Compounds Table 30.14 uses InP_xAs_ySb_1-x-y, and E0 = 0.512+0.030*y-0.183*y^2
    E_bandgap = 0.512 + 0.030*InAsSbPstruct.x-0.183*InAsSbPstruct.x^2
    if(enr>E_bandgap)
        return imag(eps_InAsSbP_xy(enr,InAsSbPstruct))
    else
        return 0
    end
end
precompile(eps_InAsSbP_imag_xy,(Float64,InAsSbPDsc))
end


text:
"""
	clean_proba_(probability:Number, atol=ATOL)

Checks whether a (complex) number is close enough to a valid probability with tolerance
`ATOL`. If so, convert it to a positive real number.
"""
function clean_proba(probability::Number, atol=ATOL)

    """checks if a (complex) number is a close enough probability
    converts it to a positive real number if complex or just positive if real"""

    not_a_proba() = error(probability , " is not a probability")

    if real(probability) >= -ATOL && real(probability) <= 1 + ATOL
        if isa(probability,Complex)
            if abs(imag(probability)) <= ATOL
                return abs(probability)
            else
                not_a_proba()
            end
        elseif isa(probability,Real)
            return abs(probability)
        else
            error(typeof(probability), " probability type undefined")
        end
    else
        not_a_proba()
    end
end

"""
	clean_pdf(A::Array, atol=ATOL)

Checks wether an array is an acceptable discrete probability distribution with
tolerance `ATOL`. If so, converts its elements to normalized positive real numbers.
"""
function clean_pdf(A::Array, atol = ATOL)

	"""checks if an array has all elements as acceptable probabilities within atol
	and converts them to that and summing to one within length(A) * atol
	and renormalizes"""

	A .= clean_proba.(A)
	normalization = sum(A)
	if isapprox(normalization, 1, atol = length(A) * atol)
		A .= 1/normalization * A
		return convert(Vector{Real}, A)
	else
		error("A not normalized")
	end
end

"""
	isa_pdf(pdf)

Asserts if `pdf`	is a valid probability distribution.
"""
function isa_pdf(pdf)

	"""asserts if pdf is a probability distribution"""
	clean_pdf(pdf)

end

"""
	tvd(a,b)

Computes the total variation distance between two probability distributions.
"""
function tvd(a,b)
	"""total variation distance"""
	sum(abs.(a-b))
end

"""
	sqr(a,b)

Computes the euclidian distance between two probability distributions.
"""
function sqr(a,b)
	"""euclidian distance"""
	sqrt(sum((a-b).^2))
end


text:
using PlanktonIndividuals

grid = RectilinearGrid(size = (16, 16, 16), x = (0,32), y = (0,32), z = (0,-32))

model = PlanktonModel(CPU(), grid; mode = CarbonMode())

function tot_mass(nut, g)
    mass = zeros(g.Nx, g.Ny, g.Nz)
    for i in 1:g.Nx
        for j in 1:g.Ny
            for k in 1:g.Nz
                mass[i,j,k] = nut[i+g.Hx, j+g.Hy, k+g.Hz] * PlanktonIndividuals.Grids.volume(i+g.Hx, j+g.Hy, k+g.Hz, g)
            end
        end
    end
    return sum(mass)
end

TC = tot_mass(model.nutrients.DIC.data, grid) +
     tot_mass(model.nutrients.DOC.data, grid) +
     tot_mass(model.nutrients.POC.data, grid)
TC = TC + sum(model.individuals.phytos.sp1.data.Bm)

uvel = zeros(16,16,16,11)
vvel = zeros(16,16,16,11)
wvel = zeros(16,16,17,11)

for i in 1:11
    uvel[:,:,:,i] .= randn(16,16,16) .* 1e-4
    vvel[:,:,:,i] .= randn(16,16,16) .* 1e-4
    wvel[:,:,:,i] .= randn(16,16,17) .* 1e-4
end

# add boundary conditions for DOC
model.nutrients.DON.bc.west  = 1.0e-3 # west boundary condition of 0.1 mmol/m^2/second
model.nutrients.DON.bc.east  = randn(16,16) .* 1e-3 # east boundary condition of -0.1 mmol/m^2/second
model.nutrients.DON.bc.south = randn(16,16,10) .* 1e-3 # south boundary condition of 0.1 mmol/m^2/second
model.nutrients.DON.bc.north = randn(16,16,10) .* 1e-3 # north boundary condition of -0.1 mmol/m^2/second

sim = PlanktonSimulation(model, ΔT = 60, iterations = 10, vels=(u=uvel, v=vvel, w=wvel))

update!(sim)

TCt = tot_mass(model.nutrients.DIC.data, grid) +
      tot_mass(model.nutrients.DOC.data, grid) +
      tot_mass(model.nutrients.POC.data, grid)
TCt = TCt + sum(model.individuals.phytos.sp1.data.Bm)

@testset "PlanktonIndividuals 3D tests:" begin
    @test isapprox(TC,TCt; atol=1e2)
end


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

"""
    Kriging(var₁=>param₁, var₂=>param₂, ...)

A polyalgorithm Kriging estimation solver.

Each pair `var=>param` specifies the `KrigingParam` `param`
for the Kriging variable `var`. In order to avoid boilerplate
code, the constructor expects pairs of `Symbol` and `NamedTuple`
instead.

## Parameters

* `variogram` - Variogram model (default to `GaussianVariogram()`)
* `mean`      - Simple Kriging mean
* `degree`    - Universal Kriging degree
* `drifts`    - External Drift Kriging drift functions

Latter options override former options. For example, by specifying
`drifts`, the user is telling the algorithm to ignore `degree` and
`mean`. If no option is specified, Ordinary Kriging is used by
default with the `variogram` only.

* `maxneighbors` - Maximum number of neighbors (default to `nothing`)
* `neighborhood` - Search neighborhood (default to `nothing`)
* `distance`     - Distance used to find nearest neighbors (default to `Euclidean()`)

The `maxneighbors` option can be used to perform approximate Kriging
with a subset of data points per estimation location. Two neighborhood
search methods are available depending on the value of `neighborhood`:

* If a `neighborhood` is provided, local Kriging is performed by sliding
  the `neighborhood` in the domain.

* If `neighborhood` is not provided, the Kriging system is built using
  `maxneighbors` nearest neighbors according to a `distance`.

## Examples

Solve the variable `:var₁` with Simple Kriging by specifying
the `mean`, and the variable `:var₂` with Universal Kriging
by specifying the `degree` and the `variogram` model.

```julia
julia> Kriging(
  :var₁ => (mean=1.,),
  :var₂ => (degree=1, variogram=SphericalVariogram(range=20.))
)
```

Solve all variables of the problem with the default parameters
(i.e. Ordinary Kriging with unit Gaussian variogram):

```julia
julia> Kriging()
```
"""
@estimsolver Kriging begin
  @param variogram = GaussianVariogram()
  @param mean = nothing
  @param degree = nothing
  @param drifts = nothing
  @param minneighbors = 1
  @param maxneighbors = nothing
  @param neighborhood = nothing
  @param distance = Euclidean()
end

function preprocess(problem::EstimationProblem, solver::Kriging)
  # retrieve problem info
  pdomain = domain(problem)

  # result of preprocessing
  preproc = Dict{Symbol,NamedTuple}()

  for (var, V) in variables(problem)
    # get user parameters
    if var ∈ keys(solver.params)
      varparams = solver.params[var]
    else
      varparams = KrigingParam()
    end

    # determine which Kriging variant to use
    if varparams.drifts ≠ nothing
      estimator = ExternalDriftKriging(varparams.variogram, varparams.drifts)
    elseif varparams.degree ≠ nothing
      estimator = UniversalKriging(varparams.variogram, varparams.degree, ndims(pdomain))
    elseif varparams.mean ≠ nothing
      estimator = SimpleKriging(varparams.variogram, varparams.mean)
    else
      estimator = OrdinaryKriging(varparams.variogram)
    end

    # determine minimum/maximum number of neighbors
    minneighbors = varparams.minneighbors
    maxneighbors = varparams.maxneighbors

    # determine path and neighborhood search method
    if varparams.maxneighbors ≠ nothing
      # locations with data for variable
      varlocs = collect(keys(datamap(problem, var)))

      if varparams.neighborhood ≠ nothing
        # local search with a neighborhood
        neigh = varparams.neighborhood

        # create a path from the data and outwards
        # use at most 10^2 points to generate path
        N = length(varlocs); M = ceil(Int, N/10^2)
        path = SourcePath(pdomain, view(varlocs,1:M:N))

        searcher  = NeighborhoodSearcher(pdomain, neigh)
        bsearcher = BoundedSearcher(searcher, maxneighbors)
      else
        # nearest neighbor search with a distance
        distance = varparams.distance
        path = LinearPath(pdomain)
        bsearcher = NearestNeighborSearcher(pdomain, maxneighbors,
                                            locations=varlocs, metric=distance)
      end
    else
      # use all data points as neighbors
      path = LinearPath(pdomain)
      bsearcher = nothing
    end

    # save preprocessed input
    preproc[var] = (estimator=estimator, path=path,
                    minneighbors=minneighbors,
                    maxneighbors=maxneighbors,
                    bsearcher=bsearcher)
  end

  preproc
end

function solve(problem::EstimationProblem, solver::Kriging)
  # preprocess user input
  preproc = preprocess(problem, solver)

  # results for each variable
  μs = []; σs = []
  for (var, V) in variables(problem)
    if preproc[var].maxneighbors ≠ nothing
      # perform Kriging with reduced number of neighbors
      varμ, varσ = solve_approx(problem, var, preproc)
    else
      # perform Kriging with all data points as neighbors
      varμ, varσ = solve_exact(problem, var, preproc)
    end

    push!(μs, var => varμ)
    push!(σs, var => varσ)
  end

  EstimationSolution(domain(problem), Dict(μs), Dict(σs))
end

function solve_approx(problem::EstimationProblem, var::Symbol, preproc)
    # retrieve problem info
    pdata = data(problem)
    pdomain = domain(problem)

    # unpack preprocessed parameters
    estimator, path, minneighbors, maxneighbors, bsearcher = preproc[var]

    # determine value type
    V = variables(problem)[var]

    # pre-allocate memory for result
    varμ = Vector{V}(undef, npoints(pdomain))
    varσ = Vector{V}(undef, npoints(pdomain))

    # pre-allocate memory for coordinates
    xₒ = MVector{ndims(pdomain),coordtype(pdomain)}(undef)

    # pre-allocate memory for neighbors
    neighbors = Vector{Int}(undef, maxneighbors)
    X = Matrix{coordtype(pdomain)}(undef, ndims(pdomain), maxneighbors)

    # keep track of estimated locations
    estimated = falses(npoints(pdomain))
    for (loc, datloc) in datamap(problem, var)
      varμ[loc] = pdata[datloc,var]
      varσ[loc] = zero(V)
      estimated[loc] = true
    end

    # estimation loop
    for location in path
      if !estimated[location]
        # coordinates of neighborhood center
        coordinates!(xₒ, pdomain, location)

        # find neighbors with previously estimated values
        nneigh = search!(neighbors, xₒ, bsearcher, mask=estimated)

        # skip location in there are too few neighbors
        if nneigh < minneighbors
          varμ[location] = NaN
          varσ[location] = NaN
          estimated[location] = true
          continue
        end

        # final set of neighbors
        nview = view(neighbors, 1:nneigh)

        # update neighbors coordinates
        coordinates!(X, pdomain, nview)

        Xview = view(X,:,1:nneigh)
        zview = view(varμ, nview)

        # fit estimator to data
        krig = fit(estimator, Xview, zview)

        # save mean and variance
        μ, σ² = predict(krig, xₒ)
        varμ[location] = μ
        varσ[location] = σ²

        # mark location as estimated and continue
        estimated[location] = true
      end
    end

    varμ, varσ
end

function solve_exact(problem::EstimationProblem, var::Symbol, preproc)
    # retrieve problem info
    pdata = data(problem)
    pdomain = domain(problem)

    # unpack preprocessed parameters
    estimator, path, minneighbors, maxneighbors, bsearcher = preproc[var]

    # determine value type
    V = variables(problem)[var]

    # pre-allocate memory for result
    varμ = Vector{V}(undef, npoints(pdomain))
    varσ = Vector{V}(undef, npoints(pdomain))

    # pre-allocate memory for coordinates
    xₒ = MVector{ndims(pdomain),coordtype(pdomain)}(undef)

    # fit estimator to data
    X, z = valid(pdata, var)
    krig = fit(estimator, X, z)

    for location in path
      coordinates!(xₒ, pdomain, location)

      μ, σ² = predict(krig, xₒ)

      varμ[location] = μ
      varσ[location] = σ²
    end

    varμ, varσ
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://discourse.julialang.org/t/fast-performance-of-array-comprehension-without-allocations/65352
#
# https://github.com/JuliaSIMD/LoopVectorization.jl/issues/307

# %%
println("Julia v$VERSION")
using Pkg
Pkg.status("SLEEFPirates")
Pkg.status("LoopVectorization")

# %%
using BenchmarkTools
using LoopVectorization

function assign!(u, x, y, z)
    @turbo for k in 1:length(z)
        for j in 1:length(y)
            for i in 1:length(x)
                u[i, j, k] = sin(x[i]) + sin(y[j]) + sin(z[k])
            end
        end
    end
end

function assign_broadcast!(uv, X, Y, Z, xx, yy, zz)
    @. X = sin(xx)
    @. Y = sin(yy)
    @. Z = sin(zz)
    @. uv = X + Y + Z
end

itot = 384
dx = 1. / itot
x = dx*collect(0:itot-1); y = dx*collect(0:itot-1); z = dx*collect(0:itot-1)
u = zeros(itot+8, itot+8, itot+8)

uv = @view u[5:5+itot-1, 5:5+itot-1, 5:5+itot-1]
xx, yy, zz = reshape(x, (:, 1, 1)), reshape(y, (1, :, 1)), reshape(z, (1, 1, :))
X, Y, Z = similar(xx), similar(yy), similar(zz)

assign!(uv, x, y, z)
a = deepcopy(uv)
uv[:, :, :] = [ sin(x) + sin(y) + sin(z) for x=x, y=y, z=z ]
b = deepcopy(uv)
@. uv = sin(xx) + sin(yy) + sin(zz)
c = deepcopy(uv)
assign_broadcast!(uv, X, Y, Z, xx, yy, zz)
d = deepcopy(uv)
@show a ≈ b ≈ c ≈ d

print("LoopVectorization.@turbo:")
@btime assign!($uv, $x, $y, $z)
print("comprehension:           ")
@btime $uv[:, :, :] = [ sin(x) + sin(y) + sin(z) for x=$x, y=$y, z=$z ]
print("simple broadcast:        ")
@btime @. $uv = sin($xx) + sin($yy) + sin($zz)
print("assign_broadcast!:       ")
@btime assign_broadcast!($uv, $X, $Y, $Z, $xx, $yy, $zz);

# %%

max_stars_repo_path: 0013/LoopVectorization example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0013/LoopVectorization example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0013/LoopVectorization example.jl

max_forks_repo_name: genkuroki/public


text:
##

using ACE
using StaticArrays, Random, Printf, Test, LinearAlgebra, ACE.Testing
using ACE: evaluate, evaluate_d, SymmetricBasis, PIBasis, O3, rand_vec3
using ACE.Random: rand_rot, rand_refl
using ACEbase.Testing: fdtest, println_slim
using ACE.Testing: __TestSVec

# using Profile, ProfileView

# Extra using Wigner for computing Wigner Matrix
using ACE.Wigner
using ACE.Wigner: get_orbsym

##
# construct the 1p-basis
maxdeg = 6
ord = 3
Bsel = SimpleSparseBasis(ord, maxdeg)

B1p = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg)

# generate a configuration
nX = 10

_randX() = State(rr = rand_vec3(B1p["Rn"]))
Xs = [ _randX() for _=1:nX ]
cfg = ACEConfig(Xs)

## Keep for futher profiling

maxdeg = 10
ord = 2
Bsel = SimpleSparseBasis(ord, maxdeg)

φ = ACE.SphericalMatrix(2,2; T = ComplexF64)
B1p = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg)
basis = SymmetricBasis(φ, B1p, Bsel)
@show length(basis)

Xs = [ _randX() for _=1:200 ]

##

using BenchmarkTools

# @btime evaluate($basis, $cfg)

# @profview begin
#    for _ = 1:10_000
#       evaluate(basis, cfg)
#    end
# end

##

c = randn(length(basis)) ./ (1:length(basis)).^2
m = ACE.LinearACEModel(basis, c)

# @btime evaluate($m, $cfg)

# @profview begin
#    for _ = 1:40_000
#       evaluate(m, cfg)
#    end
# end

##

AA = evaluate(basis.pibasis, cfg)
C̃ = m.evaluator.coeffs
C = [ c.val for c in C̃ ]
C1 = randn(length(AA))

using LoopVectorization

##

_contract1(AA, C) = sum( AA[i] * C[i] for i = 1:length(AA) )

function _contract2(AA, C)
   @assert length(AA) == length(C)
   val = AA[1] * C[1]
   @inbounds for i = 2:length(AA)
      val += AA[i] * C[i]
   end
   return val
end

function _contract3(AA, C)
   @assert length(AA) == length(C)
   val = AA[1] * C[1]
   for i = 2:length(AA)
      @fastmath @inbounds val += AA[i] * C[i]
   end
   return val
end

##

@info("With Matrix")
@btime _contract1($AA, $C)
@btime _contract2($AA, $C)
@btime _contract3($AA, $C)

@info("With Scalar")
@btime _contract1($AA, $C1)
@btime _contract2($AA, $C1)
@btime _contract3($AA, $C1)

##

C = randn(296, length(AA))
C * AA

@btime $C * $AA

AAA = zeros(eltype(AA), length(AA), 512)
for i = 1:512
   AAA[:, i] = AA
end

@btime $C * $AAA

AAA1 = collect(AAA')
BB = AAA1 * C1
@btime mul!($BB, $AAA1, $C1)
@btime dot($AA, $C1)

C10 = rand(4089, 16)
BB10 = AAA1 * C10
@btime mul!($BB10, $AAA1, $C10)


text:
#
# Grid based A* shorest path planning
#
# author: Atsushi Sakai (@Atsushi_twi)
#

module grid_a_star

using PyPlot
using NearestNeighbors
using DataStructures

type Node
    x::Int64 #x index
    y::Int64 #y index
    cost::Float64 # cost
    pind::Int64 # parent index
end

function calc_dist_policy(gx::Float64, gy::Float64,
                     ox::Array{Float64}, oy::Array{Float64},
                     reso::Float64, vr::Float64)
    """
    gx: goal x position [m]
    gx: goal x position [m]
    ox: x position list of Obstacles [m]
    oy: y position list of Obstacles [m]
    reso: grid resolution [m]
    vr: vehicle radius[m]
    """

    ngoal = Node(round(Int64, gx/reso),round(Int64, gy/reso),0.0, -1)

    ox = [iox/reso for iox in ox]
    oy = [ioy/reso for ioy in oy]

    obmap, minx, miny, maxx, maxy, xw, yw = calc_obstacle_map(ox, oy, reso, vr)

    #open, closed set
    open, closed = Dict{Int64, Node}(), Dict{Int64, Node}()
    open[calc_index(ngoal, xw, minx, miny)] = ngoal

    motion = get_motion_model()
    nmotion = length(motion[:,1])
    pq = PriorityQueue()
    enqueue!(pq, calc_index(ngoal, xw, minx, miny), ngoal.cost)

    while true
        if length(open) == 0;println("Finish Search");break;end

        c_id = dequeue!(pq)
        current = open[c_id]

        delete!(open, c_id)
        closed[c_id] = current

        for i in 1:nmotion # expand search grid based on motion model
            node = Node(current.x+motion[i,1], current.y+motion[i,2], current.cost+motion[i,3], c_id)

            if !verify_node(node, minx, miny, xw, yw, obmap)
                continue
            end

            node_ind = calc_index(node, xw, minx, miny)

            # If it is already in the closed set, skip it
            if haskey(closed,node_ind)  continue end

            if haskey(open, node_ind)
                if open[node_ind].cost > node.cost
                    # If so, update the node to have a new parent
                    open[node_ind].cost = node.cost
                    open[node_ind].pind = c_id
                end
            else # add to open set
                open[node_ind] = node
                enqueue!(pq, calc_index(node, xw, minx, miny), node.cost)
            end
        end
    end

    pmap = calc_policy_map(closed, xw, yw, minx, miny)

    return pmap
end

function calc_policy_map(closed, xw, yw, minx, miny)

    pmap = fill(Inf, (xw,yw))

    for n in values(closed)
        pmap[n.x-minx, n.y-miny] = n.cost
    end
    # println(pmap)

    return pmap
end

function calc_astar_path(sx::Float64, sy::Float64, gx::Float64, gy::Float64,
                         ox::Array{Float64}, oy::Array{Float64}, reso::Float64, vr::Float64)
    """
    sx: start x position [m]
    sy: start y position [m]
    gx: goal x position [m]
    gx: goal x position [m]
    ox: x position list of Obstacles [m]
    oy: y position list of Obstacles [m]
    reso: grid resolution [m]
    """

    nstart = Node(round(Int64,sx/reso),round(Int64, sy/reso),0.0, -1)
    ngoal = Node(round(Int64, gx/reso),round(Int64, gy/reso),0.0, -1)

    ox = [iox/reso for iox in ox]
    oy = [ioy/reso for ioy in oy]

    obmap, minx, miny, maxx, maxy, xw, yw = calc_obstacle_map(ox, oy, reso, vr)

    #open, closed set
    open, closed = Dict{Int64, Node}(), Dict{Int64, Node}()
    open[calc_index(nstart, xw, minx, miny)] = nstart

    motion = get_motion_model()
    nmotion = length(motion[:,1])
    pq = PriorityQueue()
    enqueue!(pq, calc_index(nstart, xw, minx, miny), calc_cost(nstart, ngoal))

    while true
        if length(open) == 0;println("Error: No open set");break;end

        c_id = dequeue!(pq)
        current = open[c_id]

        if current.x == ngoal.x && current.y == ngoal.y # check goal
            println("Goal!!")
            closed[c_id] = current
            break
        end

        delete!(open, c_id)
        closed[c_id] = current

        for i in 1:nmotion # expand search grid based on motion model
            node = Node(current.x+motion[i,1], current.y+motion[i,2], current.cost+motion[i,3], c_id)

            if !verify_node(node, minx, miny, xw, yw, obmap)
                continue
            end

            node_ind = calc_index(node, xw, minx, miny)

            # If it is already in the closed set, skip it
            if haskey(closed,node_ind)  continue end

            if haskey(open, node_ind)
                if open[node_ind].cost > node.cost
                    # If so, update the node to have a new parent
                    open[node_ind].cost = node.cost
                    open[node_ind].pind = c_id
                end
            else # add to open set
                open[node_ind] = node
                enqueue!(pq, calc_index(node, xw, minx, miny), calc_cost(node, ngoal))
            end
        end
    end

    rx, ry = get_final_path(closed, ngoal, nstart, xw, minx, miny, reso)

    return rx, ry
end

function verify_node(node::Node, minx::Int64, miny::Int64, xw::Int64, yw::Int64, obmap::Array{Bool,2})

    if (node.x - minx) >= xw
        return false
    elseif (node.x - minx) <= 0
        return false
    end
    if (node.y - miny) >= yw
        return false
    elseif (node.y - miny) <= 0
        return false
    end

    #collision check
    if obmap[node.x-minx, node.y-miny]
        return false
    end

    return true
end

function calc_cost(n::Node, ngoal::Node)
    return (n.cost + h(n.x - ngoal.x, n.y - ngoal.y))
end

function get_motion_model()
    # dx, dy, cost
    motion=[1 0 1;
          0 1 1;
          -1 0 1;
          0 -1 1;
          -1 -1 sqrt(2);
          -1 1 sqrt(2);
          1 -1 sqrt(2);
          1 1 sqrt(2);]

    return motion
end

function calc_index(node::Node, xwidth::Int64, xmin::Int64, ymin::Int64)
    return (node.y - ymin)*xwidth + (node.x - xmin)
end

function calc_obstacle_map(ox::Array{Float64}, oy::Array{Float64}, reso::Float64, vr::Float64)

    minx = round(Int64, minimum(ox))
    miny = round(Int64, minimum(oy))
    maxx = round(Int64, maximum(ox))
    maxy = round(Int64, maximum(oy))
    # println("minx:", minx)
    # println("miny:", miny)
    # println("maxx:", maxx)
    # println("maxy:", maxy)

    xwidth = round(Int64, maxx - minx)
    ywidth = round(Int64, maxy - miny)
    # println("xwidth:", xwidth)
    # println("ywidth:", ywidth)

    obmap = fill(false, (xwidth,ywidth))

    kdtree = KDTree(hcat(ox, oy)')
    for ix in 1:xwidth
        x = ix + minx
        for iy in 1:ywidth
            y = iy + miny
            idxs, onedist = knn(kdtree, [x, y] , 1)
            if onedist[1] <= vr/reso
                obmap[ix,iy] = true
            end
        end
    end

    # println(" calc_obstacle_map done")

    return obmap, minx, miny, maxx, maxy, xwidth, ywidth
end

function get_final_path(closed::Dict{Int64, Node},
                        ngoal::Node,
                        nstart::Node,
                        xw::Int64,
                        minx::Int64,
                        miny::Int64,
                        reso::Float64)

    rx, ry = [ngoal.x],[ngoal.y]
    nid = calc_index(ngoal, xw, minx, miny)
    while true
        n = closed[nid]
        push!(rx, n.x)
        push!(ry, n.y)
        nid = n.pind

        if rx[end] == nstart.x && ry[end] == nstart.y
            println("done")
            break
        end
    end

    rx = reverse(rx) .* reso
    ry = reverse(ry) .* reso

    return rx, ry
end

function search_min_cost_node(open::Dict{Int64, Node}, ngoal::Node)
    mnode = nothing
    mcost = Inf
    for n in values(open)
        # println(n)
        cost = n.cost + h(n.x - ngoal.x, n.y - ngoal.y)
        if mcost > cost
            mnode = n
            mcost = cost
        end
    end
    # println("minnode:", mnode)

    return mnode
end

function h(x::Int64, y::Int64)
    """
    Heuristic cost function
    """
    return sqrt(x^2+y^2);
end

function main()
    println(PROGRAM_FILE," start!!")

    sx = 10.0  # [m]
    sy = 10.0  # [m]
    gx = 50.0  # [m]
    gy = 50.0  # [m]

    ox = Float64[]
    oy = Float64[]

    for i in 0:60
        push!(ox, Float64(i))
        push!(oy, 0.0)
    end
    for i in 0:60
        push!(ox, 60.0)
        push!(oy, Float64(i))
    end
    for i in 0:60
        push!(ox, Float64(i))
        push!(oy, 60.0)
    end
    for i in 0:60
        push!(ox, 0.0)
        push!(oy, Float64(i))
    end
    for i in 0:40
        push!(ox, 20.0)
        push!(oy, Float64(i))
    end
    for i in 0:40
        push!(ox, 40.0)
        push!(oy, 60.0-Float64(i))
    end

    const VEHICLE_RADIUS = 5.0 #[m]
    const GRID_RESOLUTION = 1.0 #[m]

    @time rx, ry = calc_astar_path(sx, sy, gx, gy, ox, oy, GRID_RESOLUTION, VEHICLE_RADIUS)

    plot(ox, oy, ".k",label="obstacles")
    plot(sx, sy, "xr",label="start")
    plot(gx, gy, "xb",label="goal")
    plot(rx, ry, "-r",label="A* path")
    legend()
    grid(true)
    axis("equal")
    show()

    println(PROGRAM_FILE," Done!!")
end

if length(PROGRAM_FILE)!=0 &&
    contains(@__FILE__, PROGRAM_FILE)

    main()
end

end #module

max_stars_repo_name: Oasis-Guo/HybridAStarTrailer

max_issues_repo_name: Oasis-Guo/HybridAStarTrailer

max_forks_repo_name: Oasis-Guo/HybridAStarTrailer


text:
# ---
# title: 7. Reverse Integer
# id: problem7
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/reverse-integer/description/>
# hidden: true
# ---
#
# Given a 32-bit signed integer, reverse digits of an integer.
#
# **Note:**
# Assume we are dealing with an environment that could only store integers
# within the 32-bit signed integer range: [−231,  231 − 1]. For the purpose of
# this problem, assume that your function returns 0 when the reversed integer
# overflows.
#
#
#
# **Example 1:**
#
#
#
#     Input: x = 123
#     Output: 321
#
#
# **Example 2:**
#
#
#
#     Input: x = -123
#     Output: -321
#
#
# **Example 3:**
#
#
#
#     Input: x = 120
#     Output: 21
#
#
# **Example 4:**
#
#
#
#     Input: x = 0
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `-231 <= x <= 231 - 1`
#
#
## @lc code=start
using LeetCode

function reverse_integer(x::Int32)::Int
    s = sign(x)
    r = parse(Int, reverse(string(s * x)))
    return r <= typemax(Int32) ? s * r : 0
end
## @lc code=end


text:
"""
    calculate_conjugacy_classes(group::Bijection{String, T}; inverse=inv, prnt=true) -> conjugacy_classes::Vector{Set{String}}

Calculate all conjugacy classes {Cl(g) = {h * g * h^{-1} | h ∈ 𝒢}} for a given finite, closed group, given the `inverse` function for elements of type `T` in the group `group`.
The `inverse` function defaults to the Julia `inv` function.
`prnt` enables/disables output of the class-sizes and remaining elements while computing the classes.
"""
function calculate_conjugacy_classes(group::Bijection{String, T}; inverse=inv, prnt=true)::Vector{Set{String}} where T
    # store the labels of the conjugacy classes in an Array/Vector of Sets
    conj = Vector{Set{String}}([Set([])])

    # for performance reasons, we don't directly iterate over the keys
    all_keys = sort(collect(keys(group)))
    max_key_idx = length(all_keys)

    #=
    Direct iteration over the keys makes the function slow down at the last few conjugacy classes,
    because for large groups >95% of the elements where already calculated and are unnecessary to check again.
    Instead, we remove all keys that are already found from the all_keys array and check only the remaining ones.
    =#
    while max_key_idx > 0
        k = all_keys[1]

        g = group[k]
        g_conj = Set([k])

        # calculate all conjugate elements to g
        for λ in keys(group)
            h = group[λ]
            el = h * g * inverse(h)
            push!(g_conj, group(el))
        end
        if prnt print("Cl($k): \t|", length(g_conj), "|") end
        # store g_conj
        push!(conj, g_conj)

        # delete all keys in g_conj from all_keys and reset the maximum index
        setdiff!(all_keys, g_conj)
        max_key_idx = length(all_keys)
        if prnt println("\t remaining elements: $max_key_idx") end
    end

    # remove empty set from first position
    popfirst!(conj)
    # return all conjugacy classes
    return conj
end


text:
export Conn
import Base.show

abstract type Obj end

mutable struct Conn <: Obj
    input::Float64
    output::Float64
    state::Array{Float64}
    params::Array{Float64}
    source::Ref
    dest::Ref
end

function Conn(cfg::Dict, source::Obj, dest::Obj, params::Array{Float64})
    Conn(0.0, 0.0, zeros(cfg["n_conn_state"]), params, Ref(source), Ref(dest))
end

function get_params(c::Conn, n::Int64)
    vcat(c.params, zeros(n - length(c.params)))
end

function show(io::IO, x::Conn)
    print(io, string("Conn(", x.input, ", ", x.output, ", ", x.state, ", ",
                     x.params, ", ", x.source[].params, ", ", x.dest[].params, ")"));
end


text:
using StaticArrays

const Vec3{T} = SVector{3, T} where T


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

struct RtmApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_rtmConnect(_api::RtmApi; batch_presence_aware=nothing, presence_sub=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/rtm.connect", ["slackAuth"])
    Swagger.set_param(_ctx.query, "batch_presence_aware", batch_presence_aware)  # type Bool
    Swagger.set_param(_ctx.query, "presence_sub", presence_sub)  # type Bool
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Starts a Real Time Messaging session.
Param: batch_presence_aware::Bool
Param: presence_sub::Bool
Return: Dict{String, Any}
"""
function rtmConnect(_api::RtmApi; batch_presence_aware=nothing, presence_sub=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_rtmConnect(_api; batch_presence_aware=batch_presence_aware, presence_sub=presence_sub, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function rtmConnect(_api::RtmApi, response_stream::Channel; batch_presence_aware=nothing, presence_sub=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_rtmConnect(_api; batch_presence_aware=batch_presence_aware, presence_sub=presence_sub, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export rtmConnect


text:
#This function computes the pressure.

include(mapping_functions.jl)
include(pressure_functions.jl)

function GasLatticePressure(steps=20,startarr=0,randomize=true,dim=(100,100))
    randomize ? Grid=rand([1 2 3 5 7],dim) : Grid=startarr

    filtertop  =   [0 1 0;0 0 0;0 0 0]
    filterbot  =   [0 0 0;0 0 0;0 1 0]
    filterleft =   [0 0 0;1 0 0;0 0 0]
    filterright=   [0 0 0;0 0 1;0 0 0]
    rang=2:dim[1]+1

    p1=Progress(400)
    collisionscount=[]   #Array for storing the number of collision counts per n. The returned array.
    for number=100:600
        pressurearray=[]     #Array for storing the number of collision counts per iterations for fixed n.
        for iterations=1:10
            iteration_pressure_array=[]  #Array for storing the number of collisions per step for an iteration.
            Grid=arraygenerator(number,dim)
            for k=1:steps
                    Grid_top=map(hasbot,conv2(Grid,filtertop)[rang,rang])

                    Grid_left=map(hasright,conv2(Grid,filterleft)[rang,rang])
                    Grid_bot=map(hastop,conv2(Grid,filterbot)[rang,rang])
                    Grid_right=map(hasleft,conv2(Grid,filterright)[rang,rang])

                    Grid_post_filter=Grid_top .* Grid_left .* Grid_right .* Grid_bot
                    Grid_pre_collision=mapzeroes(Grid,Grid_post_filter)
                    Grid_pre_wall=map(collision,Grid_pre_collision)

                    #This is the function for calling the collisions per frame/step. It has to be called before
                    #the border collisions function.
                    iteration_pressure_array=push!(iteration_pressure_array,countcollisions(Grid_pre_wall))
                    Grid=map(bordercollisions,Grid_pre_wall)

            end
            #The average of the collisions per step is taken.
            pressurearray=push!(pressurearray,average_of_array(iteration_pressure_array))
        end
        #The average of the collisions per iteration is taken.
        collisionscount=push!(collisionscount,average_of_array(pressurearray))
        next!(p1)
    end

    return collisionscount
end


text:
struct PlsdaAgg
    fm
    nlv
    w_mod
    lev::AbstractVector
    ni::AbstractVector
end

"""
    plsrda_agg(X, y, weights = ones(size(X, 1)); nlv)
Aggregation of PLSR-DA models with different numbers of LVs.
* `X` : X-data.
* `y` : y-data (class membership).
* weights : Weights of the observations.
* `nlv` : A character string such as "5:20" defining the range of the numbers of LVs
    to consider ("5:20": the predictions of models with nb LVS = 5, 6, ..., 20
    are averaged). Syntax such as "10" is also allowed ("10": correponds to
    the single model with 10 LVs).

Ensemblist method where the predictions are calculated by "averaging"
the predictions of a set of models built with different numbers of
latent variables (LVs).

For instance, if argument `nlv` is set to `nlv = "5:10"`, the prediction for
a new observation is the most occurent class within the predictions
returned by the models with 5 LVS, 6 LVs, ... 10 LVs, respectively.
"""
function plsrda_agg(X, y, weights = ones(size(X, 1)); nlv)
    n = size(X, 1)
    p = size(X, 2)
    nlv = eval(Meta.parse(nlv))
    nlvmax = maximum(nlv)
    nlv = (max(minimum(nlv), 0):min(nlvmax, n, p))
    w = ones(nlvmax + 1)
    w_mod = mweights(w[collect(nlv) .+ 1])   # uniform weights
    fm = plsrda(X, y, weights; nlv = nlvmax)
    PlsdaAgg(fm, nlv, w_mod, fm.lev, fm.ni)
end

"""
    plslda_agg(X, y, weights = ones(size(X, 1)); nlv)
Aggregation of PLSR models with different numbers of LVs.
* `X` : X-data.
* `y` : y-data (class membership).
* weights : Weights of the observations.
* `nlv` : A character string such as "5:20" defining the range of the numbers of LVs
    to consider ("5:20": the predictions of models with nb LVS = 5, 6, ..., 20
    are averaged). Syntax such as "10" is also allowed ("10": correponds to
    the single model with 10 LVs).

Ensemblist method where the predictions are calculated by "averaging"
the predictions of a set of models built with different numbers of
latent variables (LVs).

For instance, if argument `nlv` is set to `nlv = "5:10"`, the prediction for
a new observation is the most occurent class within the predictions
returned by the models with 5 LVS, 6 LVs, ... 10 LVs, respectively.
"""
function plslda_agg(X, y, weights = ones(size(X, 1)); nlv)
    n = size(X, 1)
    p = size(X, 2)
    nlv = eval(Meta.parse(nlv))
    nlvmax = maximum(nlv)
    nlv = (max(minimum(nlv), 0):min(nlvmax, n, p))
    w = ones(nlvmax + 1)
    w_mod = mweights(w[collect(nlv) .+ 1])   # uniform weights
    fm = plslda(X, y, weights; nlv = nlvmax)
    PlsdaAgg(fm, nlv, w_mod, fm.lev, fm.ni)
end

"""
    plsqda_agg(X, y, weights = ones(size(X, 1)); nlv)
Aggregation of PLSR models with different numbers of LVs.
* `X` : X-data.
* `y` : y-data (class membership).
* weights : Weights of the observations.
* `nlv` : A character string such as "5:20" defining the range of the numbers of LVs
    to consider ("5:20": the predictions of models with nb LVS = 5, 6, ..., 20
    are averaged). Syntax such as "10" is also allowed ("10": correponds to
    the single model with 10 LVs).

Ensemblist method where the predictions are calculated by "averaging"
the predictions of a set of models built with different numbers of
latent variables (LVs).

For instance, if argument `nlv` is set to `nlv = "5:10"`, the prediction for
a new observation is the most occurent class within the predictions
returned by the models with 5 LVS, 6 LVs, ... 10 LVs, respectively.
"""
function plsqda_agg(X, y, weights = ones(size(X, 1)); nlv)
    n = size(X, 1)
    p = size(X, 2)
    nlv = eval(Meta.parse(nlv))
    nlvmax = maximum(nlv)
    nlv = (max(minimum(nlv), 0):min(nlvmax, n, p))
    w = ones(nlvmax + 1)
    w_mod = mweights(w[collect(nlv) .+ 1])   # uniform weights
    fm = plsqda(X, y, weights; nlv = nlvmax)
    PlsdaAgg(fm, nlv, w_mod, fm.lev, fm.ni)
end

"""
    predict(object::PlsrAgg, X)
Compute y-predictions from a fitted model.
* `object` : The fitted model.
* `X` : X-data for which predictions are computed.
"""
function predict(object::PlsdaAgg, X)
    X = ensure_mat(X)
    m = size(X, 1)
    nlv = object.nlv
    le_nlv = length(nlv)
    zpred = predict(object.fm, X; nlv = nlv).pred
    if(le_nlv == 1)
        pred = zpred
    else
        z = reduce(hcat, zpred)
        pred = similar(object.fm.lev, m, 1)
        @inbounds for i = 1:m
            pred[i, :] .= findmax_cla(z[i, :], object.w_mod)
        end
    end
    (pred = pred, predlv = zpred)
end


text:
export factorialseries, stirlingseries

# evaluate p[1]/(z)₁ + p[2]/(z)₂ + p[3]/(z)₃ + ⋯, i.e. a rising factorial series via downward recurrence
macro factorialseries(z, p...)
    ex = esc(p[end])
    for i = length(p)-1:-1:1
        ex = :(muladd(inv(x+$i), $ex, $(esc(p[i]))))
    end
    ex = :($ex/x)
    Expr(:block, :(x = $(esc(z))), ex)
end

# stirlingseries is a factorially resummed (convergent) version of Stirling's series for the gamma function.

stirlingseries(z::Float64,n::Val{1}) = @factorialseries(z+1,1/12)
stirlingseries(z::Float64,n::Val{2}) = @factorialseries(z+1,1/12,1/12)
stirlingseries(z::Float64,n::Val{3}) = @factorialseries(z+1,1/12,1/12,59/360)
stirlingseries(z::Float64,n::Val{4}) = @factorialseries(z+1,1/12,1/12,59/360,29/60)
stirlingseries(z::Float64,n::Val{5}) = @factorialseries(z+1,1/12,1/12,59/360,29/60,533/280)
stirlingseries(z::Float64,n::Val{6}) = @factorialseries(z+1,1/12,1/12,59/360,29/60,533/280,1577/168)
stirlingseries(z::Float64,n::Val{7}) = @factorialseries(z+1,1/12,1/12,59/360,29/60,533/280,1577/168,280361/5040)
stirlingseries(z::Float64,n::Val{8}) = @factorialseries(z+1,1/12,1/12,59/360,29/60,533/280,1577/168,280361/5040,69311/180)
stirlingseries(z::Float64,n::Val{9}) = @factorialseries(z+1,1/12,1/12,59/360,29/60,533/280,1577/168,280361/5040,69311/180,36226519/11880)
stirlingseries(z::Float64,n::Val{10}) = @factorialseries(z+1,1/12,1/12,59/360,29/60,533/280,1577/168,280361/5040,69311/180,36226519/11880,7178335/264)
stirlingseries(z::Float64,n::Val{11}) = @factorialseries(z+1,1/12,1/12,59/360,29/60,533/280,1577/168,280361/5040,69311/180,36226519/11880,7178335/264,64766889203/240240)
stirlingseries(z::Float64,n::Val{12}) = @factorialseries(z+1,1/12,1/12,59/360,29/60,533/280,1577/168,280361/5040,69311/180,36226519/11880,7178335/264,64766889203/240240,32128227179/10920)
stirlingseries(z::Float64,n::Val{13}) = @factorialseries(z+1,1/12,1/12,59/360,29/60,533/280,1577/168,280361/5040,69311/180,36226519/11880,7178335/264,64766889203/240240,32128227179/10920,459253205417/13104)
stirlingseries(z::Float64,n::Val{14}) = @factorialseries(z+1,1/12,1/12,59/360,29/60,533/280,1577/168,280361/5040,69311/180,36226519/11880,7178335/264,64766889203/240240,32128227179/10920,459253205417/13104,325788932161/720)
stirlingseries(z::Float64,n::Val{15}) = @factorialseries(z+1,1/12,1/12,59/360,29/60,533/280,1577/168,280361/5040,69311/180,36226519/11880,7178335/264,64766889203/240240,32128227179/10920,459253205417/13104,325788932161/720,2311165698322609/367200)

stirlingseries(z::AbstractVector,n::Val{N}) where N = [ stirlingseries(z[i],n) for i=1:length(z) ]

function stirlingnumbers(n::Int)
    s = zeros(BigInt,n)
    s[1] = 1
    for nn = 2:n
        for k = nn:-1:2
            s[k] = (nn-1)*s[k] + s[k-1]
        end
        s[1] = (nn-1)*s[1]
    end
    s
end

function stirlingcoefficient(n::Int)
    s = stirlingnumbers(n)
    cₙ = s[1]//2//3
    for k = 2:n
        cₙ += k*s[k]//(k+1)//(k+2)
    end
    cₙ//2n
end

# 1/12,1/12,59/360,29/60,533/280,1577/168,280361/5040,69311/180,36226519/11880,7178335/264,64766889203/240240,32128227179/10920,459253205417/13104,325788932161/720,2311165698322609/367200,287144996287039/3060,1215091897184850539/813960,402833263943353393/15960,476099430416027805187/1053360,236881416523193720213/27720,650730651653461090091101/3825360,58892403426105875913001/16560,391322630274754041620578051/5023200,292283822581742493926985629/163800,41922362896177303908488106751/982800,3212518869076718326256692951/3024,72400537401483369096899447543287/2630880,5153156706882006699283666106197/6960,686476646948153636599886190355425587/33227040,136841160970757197040629784385718445/229152,2086964833758413423189466175604216292877/116867520,130030515625408962047441059505505281107/235620,4978142183290000320492136142405979085829/282744,1459881140954415563253751607947838024953/2520,31663466129700972962247259201794272834399495731/1612119600
# 1//12,1//12,59//360,29//60,533//280,1577//168,280361//5040,69311//180,36226519//11880,7178335//264,64766889203//240240,32128227179//10920,459253205417//13104,325788932161//720,2311165698322609//367200,287144996287039//3060,1215091897184850539//813960,402833263943353393//15960,476099430416027805187//1053360,236881416523193720213//27720,650730651653461090091101//3825360,58892403426105875913001//16560,391322630274754041620578051//5023200,292283822581742493926985629//163800,41922362896177303908488106751//982800,3212518869076718326256692951//3024,72400537401483369096899447543287//2630880,5153156706882006699283666106197//6960,686476646948153636599886190355425587//33227040,136841160970757197040629784385718445//229152,2086964833758413423189466175604216292877//116867520,130030515625408962047441059505505281107//235620,4978142183290000320492136142405979085829//282744,1459881140954415563253751607947838024953//2520,31663466129700972962247259201794272834399495731//1612119600

max_stars_repo_name: MikaelSlevinsky/SequenceTransformations

max_issues_repo_name: MikaelSlevinsky/SequenceTransformations

max_forks_repo_name: MikaelSlevinsky/SequenceTransformations


text:
using Grader
using Test

@testset "Grader.jl" begin
    # Test whether evalasmodule is working correctly for a basic problem.
    @testset "evalasmodule" begin
        golden_code = """x=2
        y = x * 3
        """

        student_code = """x=2
        y = x + x + x
        z = 4
        """

        golden_result = Grader.evalasmodule(golden_code)
        student_result = Grader.evalasmodule(student_code)

        @testset "correct x" begin
            @test golden_result.x == student_result.x
        end
        @testset "correct y" begin
            @test golden_result.y ≈ student_result.y
        end
        @testset "correct z only in student code" begin
            @test_throws UndefVarError golden_result.z
            @test student_result.z == 4
        end

        @testset "Throws ParseError for syntax error" begin
            syntaxerr = "x= @@@"
            @test_throws Base.Meta.ParseError Grader.evalasmodule(syntaxerr)
        end
    end

    @testset "grade problem" begin
        @testset "correct answer" begin
            goldencode = """x=2
            y = x * 3
            """
            studentcode = """x=2
            y = x + x + x
            """

            p = Problem()
            golden = rungolden!(p, goldencode)
            student = runstudent!(p, studentcode)

            grade!(p, "y", "check y", 2, :($student.y ≈ $golden.y), "y is incorrect")

            @test length(p.tests) == 1
            @test p.tests[1].max_points ≈ 2
            @test p.tests[1].points ≈ 2
            @test p.score ≈ 1
        end

        @testset "partially correct answer" begin
            goldencode = """x=2
            y = x * 3
            z = y + 2
            """
            studentcode = """x=2
            y = x + x + x
            z = y + 3
            """

            p = Problem()
            golden = rungolden!(p, goldencode)
            student = runstudent!(p, studentcode)

            grade!(p, "y", "check y", 2, :($student.y ≈ $golden.y), "y is incorrect")
            grade!(p, "z", "check z", 8, :($student.z ≈ $golden.z), "z is incorrect")

            @test length(p.tests) == 2
            @test p.tests[1].max_points ≈ 2
            @test p.tests[2].max_points ≈ 8
            @test p.tests[1].points ≈ 2
            @test p.tests[2].points ≈ 0
            @test p.tests[2].message == "z is incorrect"
            @test p.score ≈ 0.2
        end

        @testset "error in student code" begin
            goldencode = """x=2
            y = x * 3
            """
            studentcode = """x=2
            y = x + x + xx
            """

            p = Problem()
            golden = rungolden!(p, goldencode)
            student = runstudent!(p, studentcode)

            grade!(p, "y", "check y", 2, :($student.y ≈ $golden.y), "y is incorrect")

            @test !p.gradable
            @test p.message == "There was an error running your code, please see information below."
        end

        @testset "missing value in student code" begin
            goldencode = """x=2
            y = x * 3
            """
            studentcode = """x=2
            y = missing
            """

            p = Problem()
            golden = rungolden!(p, goldencode)
            student = runstudent!(p, studentcode)

            grade!(p, "y", "check y", 2, :($student.y ≈ $golden.y), "y is incorrect")

            @test p.gradable
            @test p.tests[1].message == "your answer contains a 'missing' value"
        end

        @testset "error in golden code" begin
            goldencode = """x=2
            y = x * 3y
            """
            studentcode = """x=2
            y = x + x + x
            """

            p = Problem()
            golden = rungolden!(p, goldencode)
            student = runstudent!(p, studentcode)

            @test !p.gradable
            @test p.message == "Internal grading error, please notify instructor."
        end

        @testset "student missing variable" begin
            goldencode = """x=2
            y = x * 3
            """
            studentcode = """x=2
            """

            p = Problem()
            golden = rungolden!(p, goldencode)
            student = runstudent!(p, studentcode)

            grade!(p, "y", "check y", 2, :($student.y ≈ $golden.y), "y is incorrect")

            @test p.tests[1].output == "\nUndefVarError: y not defined"
            @test p.tests[1].message == "y is incorrect"
        end

        @testset "type conversion error" begin
            studentcode = """
            function fib(n::Int)::Int
                nothing
            end
            """
            goldencode = """
            function fib(n::Int)::Int
                if n <= 1
                    return n
                else
                    return fib(n - 1) + fib(n - 2)
                end
            end
            """

            p = Problem()
            golden = rungolden!(p, goldencode)
            student = runstudent!(p, studentcode)

            Grader.grade!(p, "fib(1)", "Check fib(1)", 1, :($student.fib(1) ≈ $golden.fib(1)), "fib(1) is incorrect")
            Grader.grade!(p, "fib(7)", "Check fib(7)", 2, :($student.fib(7) ≈ $golden.fib(7)), "fib(7) is incorrect")
            Grader.grade!(p, "fib(7)", "Check fib(7)", 2, :($golden.fib(7) ≈ $golden.fib(7)), "fib(7) is incorrect")

            @test p.score ≈ 2.0/5.0
            @test occursin( "Cannot `convert`", p.tests[1].output)
            @test p.tests[1].message == "fib(1) is incorrect"
        end
    end

    @testset "pl_JSON" begin
        goldencode = """x=2
        y = x * 3
        z = y + 2
        """
        studentcode = """x=2
        y = x + x + x
        z = y + 3
        """

        p = Problem()
        golden = rungolden!(p, goldencode)
        student = runstudent!(p, studentcode)

        grade!(p, "y", "check y", 2, :($student.y ≈ $golden.y), "y is incorrect")
        grade!(p, "z", "check z", 8, :($student.z ≈ $golden.z), "z is incorrect")

        s = IOBuffer()
        pl_JSON(s, p)
        jsondata = String(take!(s))

        @test jsondata == """{"gradable":true,"score":0.2,"message":"","output":"","images":[],"tests":[{"name":"y","description":"check y","points":2.0,"max_points":2.0,"message":"","output":"","images":[]},{"name":"z","description":"check z","points":0.0,"max_points":8.0,"message":"z is incorrect","output":"","images":[]}]}"""

        @test length(p.tests) == 2
        @test p.tests[1].max_points ≈ 2
        @test p.tests[2].max_points ≈ 8
        @test p.tests[1].points ≈ 2
        @test p.tests[2].points ≈ 0
        @test p.tests[2].message == "z is incorrect"
        @test p.score ≈ 0.2
    end
end