text:
module JlExprPath
# preparing a special type represents Expr because leaf nodes with non-boxing value can't be referenced
ExprRef = Tuple{Expr, Int64} # (expr, idx) if idx == 0 represents expr itself, otherwise expr.args[idx]
getValue(x::ExprRef) = x[2] == 0? x[1]: x[1].args[x[2]]
setValue!(x::ExprRef, v) = if x[2] == 0; x = (v, 0) else; x[1].args[x[2]] = v; end
#Base.show(io, x::Tuple{Expr, Int64}) = print(io, get(x))

ExprVec = Vector{ExprRef}
getValue(x::ExprVec) = x[end][2] == 0? x[end][1]: x[end][1].args[x[end][2]]
setValue!(x::ExprVec, v) = if x[end][2] == 0; x[end] = (v, 0) else; x[end][1].args[x[end][2]] = v; end

ExprVecVec = Vector{ExprVec}
getValue(x::ExprVecVec) = map(getValue, x)
# make(x::Vector{Expr})
# Base.show(io, x::Vector{Tuple{Expr, Int64}}) = print(io, get(x))
Pattern = Tuple{AbstractString, AbstractString, AbstractString}
PatternVec = Vector{Pattern}
export getValue, setValue!

function evalVerbFunc(name::AbstractString, nodes::ExprVecVec, verb::AbstractString)::String
	if name == ""
		if isempty(nodes)
			return "\"\""
		else
			if nodes|>length != 1
				warn("Got 2 or more nodes from query '$verb'. the query at here should return only one node.")
			end
			v = nodes[1]|>getValue
			return  v isa Expr? "\"$(v.head)\"": "\"$v\""
		end
	elseif name == "length"
		return string(length(nodes))
	elseif name == "pos"
		if isempty(nodes)
			error("function '$name' requires just one node but query '$verb' returned no nodes")
		else
			if nodes|>length != 1
				warn("Got 2 or more nodes from query '$verb' in function '$name'. the query at here should return only one node.")
			end
			v = nodes[1]
			if length(v) < 2 # root node
				return "0"
			end
			return string(v[end][2] == 0? v[end - 1][2]: v[end][2])
		end
	end
end

function evalVerb(verb::AbstractString, curnode::ExprVec)::Bool
	if verb == nothing || verb == ""
		return true
	end
	verb = replace(verb, r"@(\d+)", s"@(*[\1])")
	r = r"@(?<name>[a-z]*)
		\(
			(?<body>(?:
				(?P<v>
					\(((?>[^\(\)]+)|(?P>v))*\)
				)
				|[^\(]
			)+?)
		\)"x
	verb = replace(verb, r, s->begin
		m = match(r, s)
		return evalVerbFunc(m[:name], findExpr(m[:body], curnode), verb)
	end)

	local ret
	try
		ret = verb|>parse|>eval
	catch
		error("$verb is illegal expression")
	end
	if ret isa Bool
		return ret
	else
		error("$verb is not boolean expression")
	end
end

function testNode(ptn::Tuple{AbstractString, AbstractString, AbstractString}, er::ExprVec)::Bool
	if ptn |> isempty; return false; end
	axis, name, verb = ptn
	m = match(r"^\s*(\d+)\s*$", verb)
	if m != nothing
		nodenum = (length(er) > 1 && er[end][2] == 0? er[end - 1][2]: er[end][2])
		if parse(Int64, m.captures[1]) != nodenum
			return false
		end
		verb = ""
	end
	m = match(r"^\s*end\s*$", verb)
	if m != nothing
		nodenum = (length(er) > 1 && er[end][2] == 0? er[end - 1][2]: er[end][2])
		nodelen = (length(er) > 1 && er[end][2] == 0? er[end - 1][1].args|>length: er[end][1].args|>length)
		if nodenum != nodelen
			return false
		end
		verb = ""
	end
	if name == "*"
		return evalVerb(verb, er)
	end

    e = getValue(er)
    if e isa Expr
        if "$(e.head)" != name; return false; end
    else
        if "$e" != name; return false; end
	end
	return evalVerb(verb, er)
end

function childAxis(ptn::Pattern, src::ExprVec)::ExprVecVec
	ret = []
	s = getValue(src)
	if !(s isa Expr)
		return []
	end
	for i in 1:length(s.args)
		cv = [src..., (s, i)]
		if testNode(ptn, cv)
			push!(ret, cv)
		end
	end
	ret
end

function ancestorAxis(ptn::Pattern, src::ExprVec)::ExprVecVec
	if src|>length < 2
		throw("no ancestors")
		return []
	end
	ret::ExprVecVec = []
	for i = 1:(src|>length - 1)
		if testNode(ptn, src[1:i])
			push!(ret,src[1:i])
		end
	end
	ret
end

function ancestorOrSelfAxis(ptn::Pattern, src::ExprVec)::ExprVecVec
	if src|>length < 1
		throw("source is empty")
		return []
	end
	ret::ExprVecVec = []
	for i = 1:length(src)
		if testNode(ptn, src[1:i])
			push!(ret,src[1:i])
		end
	end
	ret
end

# filter src and its all descendants by ptn and return vector of them
function descendantOrSelfAxis(ptn::Pattern, src::ExprVec)::ExprVecVec
	ret = []
	s = getValue(src)
	if testNode(ptn, src)
		push!(ret, src)
	end
	if !(s isa Expr)
		return ret
	end
	for i in 1:length(s.args)
		append!(ret, descendantOrSelfAxis(ptn, [src..., (s, i)]))
	end
	ret
end

function descendantAxis(ptn::Pattern, src::ExprVec)::ExprVecVec
	if !(s isa Expr)
		return []
	end
	ret = []
	s = getValue(src)
	for i in 1:length(s.args)
		append!(ret, descendantOrSelfAxis(ptn, [src..., (s, i)]))
	end
	ret
end

function findExprImp(ptn::Pattern, srcs::ExprVecVec)::ExprVecVec
    if ptn|>isempty; return srcs; end
    axis, name, verb = ptn
	ret = []
	for src in srcs
		s = getValue(src)
		if axis == "self"
			if testNode(ptn, src)
				push!(ret, src)
			end
		elseif axis == "parent"
			if src|>length > 1 && testNode(ptn, src[1:(end - 1)])
				push!(ret, src[1:(end - 1)])
			end
		elseif axis == "child"
			append!(ret, childAxis(ptn, src))
		elseif axis == "ancestor"
			append!(ret, ancestorAxis(ptn, src))
		elseif axis == "ancestor-or-self"
			append!(ret, ancestorOrSelfAxis(ptn, src))
		elseif axis == "descendant"
			append!(ret, descendantAxis(ptn, src))
		elseif axis == "descendant-or-self"
			append!(ret, descendantOrSelfAxis(ptn, src))
		elseif axis == "root"
			ret = testNode(ptn, src[1:1])? [src[1:1]]: []
			break
		elseif axis == "func"
			error("func axis is available only in verb clause(i.e. in [...])")
		else
			throw("unknown axis $axis")
		end
	end
    ret
end

function parsePathString(str::AbstractString)::PatternVec
	ret = []
	if isempty(str)
		throw("invalid path")
	end

	# split path by `/` considering the `/` is out of `[...]`
	for (i, m) in enumerate(eachmatch(r"(?<=^|/)(?:(?P<v>\[((?>[^\[\]]+)|(?P>v))*\])|[^\[/])*", str))
		s = m.match
		if s == ""
			push!(ret, ((i == 1? "root": "descendant-or-self"), "*", ""))
			continue
		end
		rm = match(r"^(?:([-a-z]+?)\s*::\s*)?([^\s\[]+?)\s*(?:\[(.+)\])?\s*$", s)
		if rm isa Void
			throw("syntax error in $(s)")
		end
		axis, name, verb = rm.captures[1:3]
		if name == "."
			axis = "self"
			name = "*"
		end
		if name == ".."
			axis = "parent"
			name = "*"
		end
		axis = axis isa Void? "child": axis
		verb = verb isa Void? "": verb
		push!(ret, (axis, name, verb))
	end
	ret
end

function findExpr(ptn::AbstractString, src::ExprVecVec)::ExprVecVec
	v = copy(src)
	for p in parsePathString(ptn)
		v = findExprImp(p, v)
		if isempty(v)
			return []
		end
	end
	v|>unique
end
findExpr(ptn::AbstractString, src::ExprVec)::ExprVecVec = findExpr(ptn, [src])
findExpr(ptn::AbstractString, src::Expr)::ExprVecVec = findExpr(ptn, [[(src, 0)]])
export findExpr
end


text:
using OmegaCore
using Test

(x::StdUnif)(ω) = ω.data[x]

function test()
  x = 1 ~ Normal(0, 1)
  y = 2 ~ Normal(0, 1)
  z = x +ₚ y

  function propose(::typeof(z), z_, ω)
    if x in ω && y ∈ ω
      return nothing
    elseif x ∈ ω

      c(x) = cdf(Normal(0, 1), x)
      Dict{Any, Any}(StdUnif(2) => c(x_ - ω.data[μ]))
    elseif StdUnif(2) in keys(ω.data)
      nothing
    end
      ## Provide a value for μ
      # @assert false
  end

  x_ = 3.5
  ω = Ω(Dict{Any, Any}(x => x_))
  propose(μ, ω)
end

function test()
  μ = 1 ~ Normal(0, 1)
  x(ω) = (2 ~ Normal(x(ω), 1))(ω)
  joint = OmegaCore.@joint μ x

  function subpropose1(qω, f, ω)
    # 1. Return new values, which are just merged,
    # with null case being empty set
    # 2. Make it mutative, modify ω
    # 3.
  end
  qω = OmegaCore.defω()
  ω = OmegaCore.defω()
  cp = CompProposal((stdpropose,))
  ω_ = propose(qω, joint, ω, cp)
end

ext: jl

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
module IsingMCMC

function genFunctor( Jd, Hd)
    Jd = Jd
    Hd = Hd

    # E(S) = - 1/2\sum_{i,j, i not j} J_{i,j} si sj - \sum_i H_i si
    function E(S, N)
        local Et = 0
        for i in 1:N
            for (j,Jv) in Jd[i]
                @inbounds Et += -S[i]*Jv*S[j]/2
            end
            @inbounds Et += -S[i]*Hd[i]
        end
        Et
    end

    # differential energy
    function dE(S, k, N)
        local dEt = 0
        for (j,Jv) in Jd[k]
            @inbounds dEt += S[k]*Jv*S[j]*2
        end
        @inbounds dEt += S[k]*Hd[k]*2
        dEt
    end

    # probability for dE
    PrE(dE, T) = exp(-dE/T)

    # flip spin at x
    function flipx!(list, x)
        list[x] *= -1
    end

    # initialize
    function initS(S, N)
        for i in 1:N
            if rand() < 0.5
                flipx!(S,i) # random flip at first
            end
        end
    end

    # MCMC
    function MCMC( S, T, N, trial)
        # initialize
        #simE = zeros(Int, div(trial,1000)+1)
        simE = zeros(Float32, div(trial,1000)+1)

        # MCMC trial
        Ec = E(S,N) # current energy
        simE[1] = Ec
        local k, de
        @inbounds for t in 1:trial
            k  = rand(1:N) # Gibbs sampling position
            de = dE(S, k, N)
            if rand() < PrE(de, T) # MH criteria
                flipx!(S, k) # change k
                Ec += de     # change E
            end
            (t % 1000 == 0) && (simE[div(t,1000)+1] = Ec)
        end

        (simE)
    end

    return (E, initS, MCMC)
end # of functor

end # of module


text:
# script for invoking the solver using the input file specified on the
# command line

using PDESolver
using EllipticEquationMod

mesh, sbp, eqn, opts = run_solver(ARGS[1])


text:
module Plotting

using Plots

using ..Structs

export plot_rectangles

rectangle(x, y, w, h) = Shape(x .+ [0, w, w, 0], y .+ [0, 0, h, h])

function plot_rectangles(rect_sizes, positions)
    p = Plots.plot(legend = false)

    for (dims, pos) in zip(rect_sizes, positions)
        Plots.plot!(
            p,
            rectangle(pos..., dims...),
            #fillcolor=nothing,
            # linewidth=0
        )
    end

    return p
end

end # module


text:
using PkgSkeleton, Test, Dates, UUIDs

# import internals for testing
using PkgSkeleton: get_replacement_values, resolve_template_dir, pkg_name_from_path,
    replace_multiple

####
#### Command line git should be installed for tests (so that they don't depend in LibGit2).
####

if !success(`git --help`)
    @info "Command line git should be installed."
    exit(1)
end

####
#### For CI, set up environment, otherwise use local settings (and assume they are defined).
####

const CI = parse(Bool, lowercase(get(ENV, "CI", "false")))

function getgitopt(opt)
    try
        chomp(read(`git config --get $(opt)`, String))
    catch
        error("couldn't get git option $(opt)")
    end
end

setgitopt(name, value) = run(`git config --global --add $(name) $(value)`)

if CI
    USERNAME = "Joe H. User"
    USEREMAIL = "test@email.domain"
    GHUSER = "somethingclever"
    setgitopt("user.name", USERNAME)
    setgitopt("user.email", USEREMAIL)
    setgitopt("github.user", GHUSER)
else
    USERNAME = getgitopt("user.name")
    USEREMAIL = getgitopt("user.email")
    GHUSER = getgitopt("github.user")
end

####
#### Mock _confirm_default since no user should be assumed to be in front of the keyboard.
####

function PkgSkeleton._confirm_default(prompt, default)
    println("$prompt ($default)>")
    default
end

####
#### test components
####

@testset "replacement values" begin
    d = Dict(get_replacement_values(; pkg_name = "FOO"))
    @test d["{UUID}"] isa UUID
    @test d["{GHUSER}"] == GHUSER
    @test d["{USERNAME}"] == USERNAME
    @test d["{USEREMAIL}"] == USEREMAIL
    @test d["{YEAR}"] == year(now())
end

@testset "template directories" begin
    default_template = abspath(joinpath(@__DIR__, "..", "templates", "default"))
    @test resolve_template_dir(:default) == default_template
    @test_throws ArgumentError resolve_template_dir(:nonexistent_builtin)
    @test resolve_template_dir(default_template) == default_template
    @test_throws ArgumentError resolve_template_dir(tempname()) # nonexistent
end

@testset "package name from path" begin
    @test pkg_name_from_path("/tmp/FooBar") == "FooBar"
    @test pkg_name_from_path("/tmp/FooBar.jl") == "FooBar"
    @test pkg_name_from_path("/tmp/FooBar/") == "FooBar"
    @test pkg_name_from_path("/tmp/FooBar.jl/") == "FooBar"
    @test_throws ArgumentError pkg_name_from_path("/tmp/FooBar.bin/")
end

@testset "multiple replaces" begin
    @test replace_multiple("{COLOR} {DISH}", ["{COLOR}" => "green", "{DISH}" => "curry"]) ==
        "green curry"
end

####
#### test generation
####

function check_dest_dir(pkg_name, dest_dir;)
    # NOTE: very rudimentary check and should be improved, checking all valid substitutions
    readme = joinpath(dest_dir, "README.md")
    @test isfile(readme)
    @test occursin(pkg_name, read(readme, String))
    mainsrc = joinpath(dest_dir, "src", pkg_name * ".jl")
    @test isfile(mainsrc)
end

@testset "package generation and checks" begin
    mktempdir() do tempdir
        dest_dir = joinpath(tempdir, "Foo")

        # test generated structure
        @test PkgSkeleton.generate(dest_dir) == true
        check_dest_dir("Foo", dest_dir)

        # run various sanity checks (mostly test contents of the template, CI will error)
        cd(dest_dir) do
            @info "test documentation (instantiation)"
            run(`julia --project=docs -e 'using Pkg; Pkg.instantiate(); Pkg.develop(PackageSpec(path=pwd()))'`)
            @info "test documentation (generation)"
            run(`julia --project=docs --color=yes docs/make.jl`)
            @info "test coverage (only instantiation)"
            run(`julia --project=test/coverage -e 'using Pkg; Pkg.instantiate()'`)
        end

        @test PkgSkeleton.generate(dest_dir) == false # will not overwrite
        @test PkgSkeleton.generate(dest_dir; skip_existing_dir = false) # will overwrite
    end
end


text:
# Michelle Medina
# August 11, 2020
# Testing Convergence of Derivative

# USING L2 NORM
H = []
E = []

for N in 50:200
    grid = Params(N,4)

	xVal = collect(range(0.0, stop=1.0, length = 129))
	y = calcNodal(x -> sin(4*pi*x), xVal)
	cosx = calcNodal(x -> 4*pi*cos(4*pi*x), xVal)
	modal = Nodal_2_H(y)

	dy = zeros(length(xVal))
	for index in CartesianIndices(xVal)
		dy[index] = derivEvaluate(modal, xVal[index])
	end

    # u = functionOnGrid(x->sin(x),grid)
    # dUdx = functionOnGrid(x->cos(x),grid) #known derivative of sin(x)
    # derivative = Deriv1D(u,grid)

    #L 2 norm: sum of all the differences between the numerical derivative and
    #the analytical one squared and then take the square root
    error = sqrt(sum((dy[2:end-1] - cosx[2:end-1]).^2))

    h = dX(grid)
    #append h to the array H and error to the array E
    append!(H,h)
    append!(E,error)

end

#plotting error vs h
loglog(H,E,"r-")

#first order convergence
loglog(H,H,"g--",label="1")

#second order convergence
loglog(H,H.^2,"b--",label="2")
xlabel("h")
ylabel("L 2 Norm")
legend(loc="best",frameon=false)
show()


text:
using PyPlot
# Set up a time grid for the plot:
ntime = 10000
time = collect(linspace(0,4pi,ntime))

# exp-sine-squared:
x1 = [0.24998970200450718,0.7528029456942368,0.9604772168743339,0.6433934213976025]
x2 = [0.4508506501130407,-0.44533765457178637,1.0810460905919905,1.1810193449724977]
x3 = [0.20659251082649457,-0.20736062231716337,1.659019056376522,0.21064077694469815]
x4 = [0.07284486969656653,-0.0729380659597069,2.431485661375288,0.047122405324731394]
x5 = [0.019204025483424515,-0.01921291003735305,3.2797651575680282,0.011419793258824862]

gamma = 1.0
model_expsin2_g1=zeros(ntime)
model_expsin2_g1 += x1[1]+x1[2].*exp(-gamma.^x1[3].*x1[4])
model_expsin2_g1 += (x2[1]+x2[2].*exp(-gamma.^x2[3].*x2[4])).*cos(2.*time)
model_expsin2_g1 += (x3[1]+x3[2].*exp(-gamma.^x3[3].*x3[4])).*cos(4.*time)
model_expsin2_g1 += (x4[1]+x4[2].*exp(-gamma.^x4[3].*x4[4])).*cos(6.*time)
model_expsin2_g1 += (x5[1]+x5[2].*exp(-gamma.^x5[3].*x5[4])).*cos(8.*time)
expsin2_g1 = exp.(-gamma.*sin(time).^2)

gamma = 3.0
model_expsin2_g3=zeros(ntime)
model_expsin2_g3 += x1[1]+x1[2].*exp(-gamma.^x1[3].*x1[4])
model_expsin2_g3 += (x2[1]+x2[2].*exp(-gamma.^x2[3].*x2[4])).*cos(2.*time)
model_expsin2_g3 += (x3[1]+x3[2].*exp(-gamma.^x3[3].*x3[4])).*cos(4.*time)
model_expsin2_g3 += (x4[1]+x4[2].*exp(-gamma.^x4[3].*x4[4])).*cos(6.*time)
model_expsin2_g3 += (x5[1]+x5[2].*exp(-gamma.^x5[3].*x5[4])).*cos(8.*time)
expsin2_g3 = exp.(-gamma.*sin(time).^2)

# Gaussian (squared-exponential):
#abest=[2.7518074963750734,1.5312499407621147,-0.34148888178746056,-0.3546256365636662,4.753504292338284,2.042571639590954,-1.3968344706386124,1.7967431129812566,-1.7641684073510928]
#abest=[2.670143,1.500258, 0.346731,-0.458705,4.125588, 1.390912, -1.211438,1.744602, 1.798436]
abest = [-0.519196,0.285478,1.383928,1.83396, 63.119512,0.0,1.479579,0.0, -32.173615,0.0,1.593691,0.0, -29.425156,0.0,1.388929,0.0]

gaussian=exp.(-time.^2/2.)
nexp = 4

model_gaussian = zeros(ntime)
for k=1:nexp
  model_gaussian += abest[1+(k-1)*4].*exp.(-abest[3+(k-1)*4].*time).*cos(abest[4+(k-1)*4].*time)
  model_gaussian += abest[2+(k-1)*4].*exp.(-abest[3+(k-1)*4].*time).*sin(abest[4+(k-1)*4].*time)
end

# Matern 5/2:
matern52=(1.0+sqrt(5.).*time+5.*time.^2./3.).*exp.(-sqrt(5.).*time)
nexp = 4
atrial = [-147.6965573,1.0097610,85.8875034,1.0621734,-4.6540017,0.8353745,67.4630744,0.9160394]
model_52 = zeros(ntime)
for k=1:nexp
  model_52 += atrial[1+(k-1)*2].*exp.(-atrial[2+(k-1)*2].*sqrt(5.)*time)
end

colors = [ "#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd", "#8c564b",
               "#e377c2", "#7f7f7f", "#bcbd22", "#17becf", ]
linewidth=3
fontsize=10

#fix,axes = subplots(2,2,figsize=[10,8])
fix,axes = subplots(1,4,figsize=[20,4])
#fix,ax = subplots()

ax = axes[1]
ax[:set_title]("(a). Simple, exact kernels",fontsize=2*fontsize)
ax[:plot](time,1.+0.*time,alpha=0.5,label=L"$1$",linewidth=linewidth,color=colors[1])
ax[:plot](time,exp.(-time/3),alpha=0.5,label=L"$e^{-t/3}$",linewidth=linewidth,color=colors[2])
ax[:plot](time,cos(time),alpha=1.0,alpha=0.5,label=L"$\cos{(\tau)}$",color=colors[3],linewidth=linewidth)
ax[:plot](time,exp.(-time./3).*cos(time),alpha=1.0,alpha=0.5,label=L"$e^{-t/3}\cos{(\tau)}$",color=colors[4],linewidth=linewidth)
ax[:legend](loc="lower center",fontsize=fontsize)
ax[:set_xlabel](L"$\tau/\tau_0$",fontsize=2*fontsize)
ax[:set_ylabel]("Dimensional ACF",fontsize=2*fontsize)
ax[:axis]([0,4pi,-1.1,1.1])

ax = axes[3]
ax[:set_title]("(b). Gaussian",fontsize=2*fontsize)
ax[:plot](time,gaussian,alpha=0.25,label=L"$G(z)$, exact",linewidth=linewidth,color=colors[1])
ax[:plot](time,model_gaussian,alpha=1.0,ls="dashed",alpha=0.5,label=L"$G(z)$, approx",color=colors[1],linewidth=linewidth)
print("ExpSq error: ",std(gaussian-model_gaussian)," ",maximum(abs(gaussian-model_gaussian)))
ax[:legend](loc="upper right",fontsize=fontsize)
ax[:set_xlabel](L"$\tau/\tau_0$",fontsize=2*fontsize)
#ax[:set_ylabel]("Dimensional ACF")
ax[:axis]([0,pi,0,1.1])

ax = axes[2]
ax[:plot](time,expsin2_g1,alpha=0.25,label=L"$E(w,1)$, exact",linewidth=linewidth,color=colors[1])
ax[:plot](time,model_expsin2_g1,alpha=1.0,label=L"$E(w,1)$, approx",ls="dashed",alpha=0.5,color=colors[1],linewidth=linewidth)
ax[:plot](time,expsin2_g3,alpha=0.25,label=L"$E(w,3)$, exact",linewidth=linewidth,color=colors[2])
ax[:plot](time,model_expsin2_g3,alpha=1.0,label=L"$E(w,3)$, approx",ls="dashed",alpha=0.5,color=colors[2],linewidth=linewidth)
ax[:set_title]("(c). Exp-sine-squared",fontsize=2*fontsize)
ax[:legend](loc="upper center",fontsize=fontsize)
ax[:set_xlabel](L"$\tau/\tau_0$",fontsize=2*fontsize)
#ax[:set_ylabel]("Dimensional ACF")
ax[:axis]([0,pi,0,1.1])

ax = axes[4]
matern32 = (1.+sqrt(3.).*time).*exp.(-sqrt(3.).*time)
b = 100.
model_32 = (1-b)*exp.(-sqrt(3.).*time)+b*exp.(-(b-1)/b.*sqrt(3.).*time)
ax[:plot](time,matern32,alpha=0.25,label=L"$K_{3/2}(\tau)$, exact",linewidth=linewidth,color=colors[1])
ax[:plot](time,model_32,alpha=1.0, label=L"$K_{3/2}(\tau)$, approx",ls="dashed",alpha=0.5,color=colors[1],linewidth=linewidth)
ax[:plot](time,matern52,alpha=0.25,label=L"$K_{5/2}(\tau)$, exact",linewidth=linewidth,color=colors[2])
ax[:plot](time,model_52,alpha=1.0, label=L"$K_{5/2}(\tau)$, approx",ls="dashed",alpha=0.5,color=colors[2],linewidth=linewidth)
ax[:set_title]("(d). Matern",fontsize=2*fontsize)
ax[:legend](loc="upper right",fontsize=fontsize)
ax[:set_xlabel](L"$\tau/\tau_0$",fontsize=2*fontsize)
#ax[:set_ylabel]("Dimensional ACF",fontsize=2*fontsize)
ax[:axis]([0,pi,0,1.1])

tight_layout()
savefig("kernel_approx.pdf", bbox_inches="tight")

max_stars_repo_name: ericagol/celerite

max_issues_repo_name: ericagol/celerite

max_forks_repo_name: ericagol/celerite


text:
using HSA

function check(message)
    println("✓ $message")
end

# Use the HSAIL Code Generator to compile a
# kernel function to BRIG

@hsa_kernel function vector_add_kernel(a::Ptr{Int64},b::Ptr{Int64})
    idx = get_global_id(Int32(0)) + 1

    x = Base.unsafe_load(b, idx)
    x += Base.unsafe_load(a, idx)
    Base.unsafe_store!(b, x, idx)
    return nothing
end
check("defining the kernel")

N =  1024*1024

a_in = Array(Int, N)
rand!(a_in)

b_out = Array(Int, N)
rand!(b_out)

expected = a_in + b_out

@hsa (N) vector_add_kernel(a_in, b_out)
check("dispatching the kernel call via @hsa")

failed = 0
for (i in 1:N)
    if (b_out[i]!=expected[i])
        failed = failed + 1

        if failed > 10
            println("...")
            break
        end

        println("Bad index: $i\n$(a_in[i]) --> $(expected[i])")
    end
end

if(failed == 0)
    println("Passed validation.")
else
    println("VALIDATION FAILED!\n")
end


text:
const MAX_ITERATION = Ref{Int}(1000)

"""
    replace!(r, zxd)
Match and replace with the rule `r`.
"""
function Base.replace!(r::AbstractRule, zxd::AbstractZXDiagram)
    matches = match(r, zxd)
    rewrite!(r, zxd, matches)
    return zxd
end

"""
    simplify!(r, zxd)
Simplify `zxd` with the rule `r`.
"""
function simplify!(r::AbstractRule, zxd::AbstractZXDiagram)
    i = 1
    matches = match(r, zxd)
    while length(matches) > 0
        rewrite!(r, zxd, matches)
        matches = match(r, zxd)
        i += 1
        if i > MAX_ITERATION.x && r in (Rule{:p2}(), Rule{:p3}(), Rule{:pab}())
            @warn "Try to simplify this ZX-diagram with rule $r more than $(MAX_ITERATION.x) iterarions"
            break
        end
    end
    return zxd
end

"""
    clifford_simplification(zxd)
Simplify `zxd` with the algorithms in [arXiv:1902.03178](https://arxiv.org/abs/1902.03178).
"""
function clifford_simplification(circ::ZXDiagram)
    zxg = ZXGraph(circ)
    zxg = clifford_simplification(zxg)
    return zxg
end

function clifford_simplification(zxg::ZXGraph)
    simplify!(Rule{:lc}(), zxg)
    simplify!(Rule{:p1}(), zxg)
    match_id = match(Rule{:id}(), zxg)
    while length(match_id) > 0
        rewrite!(Rule{:id}(), zxg, match_id)
        simplify!(Rule{:lc}(), zxg)
        simplify!(Rule{:p1}(), zxg)
        match_id = match(Rule{:id}(), zxg)
    end
    replace!(Rule{:pab}(), zxg)

    return zxg
end

function clifford_simplification(bir::BlockIR)
    zxd = convert_to_zxd(bir)
    zxg = clifford_simplification(zxd)
    chain = circuit_extraction(zxg)
    return BlockIR(bir.parent, bir.nqubits, chain)
end

function full_reduction(cir::ZXDiagram)
    zxg = ZXGraph(cir)
    zxg = full_reduction(zxg)
    return zxg
end

function full_reduction(zxg::ZXGraph)
    simplify!(Rule{:lc}(), zxg)
    simplify!(Rule{:p1}(), zxg)
    simplify!(Rule{:p2}(), zxg)
    simplify!(Rule{:p3}(), zxg)
    replace!(Rule(:pab), zxg)
    simplify!(Rule{:p1}(), zxg)
    match_id = match(Rule{:id}(), zxg)
    match_gf = match(Rule{:gf}(), zxg)
    while length(match_id) + length(match_gf) > 0
        rewrite!(Rule{:id}(), zxg, match_id)
        rewrite!(Rule{:gf}(), zxg, match_gf)
        simplify!(Rule{:lc}(), zxg)
        simplify!(Rule{:p1}(), zxg)
        simplify!(Rule{:p2}(), zxg)
        simplify!(Rule{:p3}(), zxg)
        replace!(Rule(:pab), zxg)
        simplify!(Rule{:p1}(), zxg)
        match_id = match(Rule{:id}(), zxg)
        match_gf = match(Rule{:gf}(), zxg)
    end

    return zxg
end

function full_reduction(bir::BlockIR)
    zxd = convert_to_zxd(bir)
    zxg = full_reduction(zxd)
    chain = circuit_extraction(zxg)
    return BlockIR(bir.parent, bir.nqubits, chain)
end

function compose_permutation(p1::Dict{Int, Int}, p2::Dict{Int, Int})
    p2 = copy(p2)
    p = Dict{Int, Int}()
    for (k1, v1) in p1
        if haskey(p2, v1)
            p[k1] = p2[v1]
            delete!(p2, v1)
        else
            p[k1] = v1
        end
    end
    for (k2, v2) in p2
        p[k2] = v2
    end
    return p
end

map_locations(d::Dict{T, T}, g::Gate) where {T <: Integer} =
    Gate(g.operation, map_locations(d, g.locations))
map_locations(d::Dict{T, T}, cg::Ctrl) where {T <: Integer} =
    Ctrl(map_locations(d, cg.gate), map_locations(d, cg.ctrl))
map_locations(d::Dict{T, T}, locs::Locations) where {T <: Integer} =
    Locations(Tuple(haskey(d, i) ? d[i] : i for i in locs.storage))
map_locations(d::Dict{T, T}, locs::CtrlLocations) where {T <: Integer} =
    CtrlLocations(map_locations(d, locs.storage))

function simplify_swap!(qc::Chain; replace_swap::Bool = true)
    chain_after_swap = Chain()
    loc_map = Dict{Int, Int}()
    while length(qc.args) > 0
        g = pop!(qc.args)
        if g isa Gate
            if g.operation === SWAP
                loc1, loc2 = plain(g.locations)[1:2]
                loc_map = compose_permutation(Dict{Int, Int}(loc1 => loc2, loc2 => loc1), loc_map)
                continue
            end
        end
        pushfirst!(chain_after_swap.args, map_locations(loc_map, g))
    end
    chain_swap = generate_swap(loc_map)
    replace_swap && replace_swap!(chain_swap, chain_after_swap)
    for g in chain_swap.args
        push!(qc.args, g)
    end
    for g in chain_after_swap.args
        push!(qc.args, g)
    end
    return qc
end

function generate_swap(loc_map::Dict{T, T}) where {T <: Integer}
    perm = copy(loc_map)
    subperms = []
    while length(perm) > 0
        k, v = pop!(perm)
        k == v && continue
        sp = [k]
        while !(v in sp)
            push!(sp, v)
            v = pop!(perm, v)
        end
        push!(subperms, sp)
    end

    qc_swap = []
    for subperm in subperms
        for k in 2:length(subperm)
            push!(qc_swap, Gate(SWAP, Locations((subperm[1], subperm[k]))))
        end
    end
    return Chain(qc_swap...)
end

function replace_swap!(chain_swap::Chain, chain_after_swap::Chain)
    qc_swap = chain_swap.args
    qc_after_swap = chain_after_swap.args
    while length(qc_swap) > 0
        g_swap = pop!(qc_swap)
        loc1, loc2 = g_swap.locations.storage[1:2]
        qmap = Dict(loc1 => loc2, loc2 => loc1)
        j = 1
        while j <= length(qc_after_swap)
            g = qc_after_swap[j]
            j += 1
            if g isa Ctrl
                if g.gate === X && length(g.gate.locations) == 1 && length(g.ctrl) == 1
                    loc = plain(g.gate.locations)[]
                    ctrl = plain(g.ctrl.storage)[]
                    if haskey(qmap, loc) && haskey(qmap, ctrl)
                        insert!(qc, j, convert_to_gate(Val(:CNOT), ctrl, loc))
                        break
                    end
                end
            end
            qc_after_swap[j-1] = map_locations(qmap, qc_after_swap[j-1])
        end
        if j > length(qc_after_swap)
            push!(qc_after_swap, convert_to_gate(Val(:CNOT), loc1, loc2))
            push!(qc_after_swap, convert_to_gate(Val(:CNOT), loc2, loc1))
            push!(qc_after_swap, convert_to_gate(Val(:CNOT), loc1, loc2))
        end
    end

    return chain_swap, chain_after_swap
end


text:
using Gadfly, DataArrays, Datetime, RDatasets

economics = dataset("ggplot2", "economics")
dates = Date[date(d) for d in economics[:Date]]

try
    economics[:Date] = dates
catch
    economics["Date"] = dates
end

p = plot(economics, x=:Date, y=:Unemploy, Geom.line)


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((e*x)^m*(a+b*x^2)^3*(A+B*x^2)*(c+d*x^2), x) == :(a*e^-5*(e*x)^(5+m)*(5+m)^-1*(B*a*(a*d+3*b*c)+3*A*b*(a*d+b*c))+b*e^-7*(e*x)^(7+m)*(7+m)^-1*(A*b*(b*c+3*a*d)+3*B*a*(a*d+b*c))+a^2*e^-3*(e*x)^(3+m)*(3+m)^-1*(A*a*d+B*a*c+3*A*b*c)+b^2*e^-9*(e*x)^(9+m)*(9+m)^-1*(A*b*d+B*b*c+3*B*a*d)+A*c*a^3*e^-1*(e*x)^(1+m)*(1+m)^-1+B*d*b^3*e^-11*(e*x)^(11+m)*(11+m)^-1)
@test integrate((e*x)^m*(a+b*x^2)^2*(A+B*x^2)*(c+d*x^2), x) == :(e^-5*(e*x)^(5+m)*(5+m)^-1*(A*b*(b*c+2*a*d)+B*a*(a*d+2*b*c))+a*e^-3*(e*x)^(3+m)*(3+m)^-1*(A*a*d+B*a*c+2*A*b*c)+b*e^-7*(e*x)^(7+m)*(7+m)^-1*(A*b*d+B*b*c+2*B*a*d)+A*c*a^2*e^-1*(e*x)^(1+m)*(1+m)^-1+B*d*b^2*e^-9*(e*x)^(9+m)*(9+m)^-1)
@test integrate((e*x)^m*(A+B*x^2)*(a+b*x^2)*(c+d*x^2), x) == :(e^-5*(e*x)^(5+m)*(5+m)^-1*(A*b*d+B*a*d+B*b*c)+e^-3*(e*x)^(3+m)*(3+m)^-1*(A*a*d+A*b*c+B*a*c)+A*a*c*e^-1*(e*x)^(1+m)*(1+m)^-1+B*b*d*e^-7*(e*x)^(7+m)*(7+m)^-1)
@test integrate((e*x)^m*(A+B*x^2)*(c+d*x^2), x) == :(e^-3*(e*x)^(3+m)*(3+m)^-1*(A*d+B*c)+A*c*e^-1*(e*x)^(1+m)*(1+m)^-1+B*d*e^-5*(e*x)^(5+m)*(5+m)^-1)
@test integrate((e*x)^m*(a+b*x^2)^3*(c+d*x^2)^2*(A+B*x^2), x) == :(e^-7*(e*x)^(7+m)*(7+m)^-1*(A*b*(b^2*c^2+3*a^2*d^2+6*a*b*c*d)+B*a*(a^2*d^2+3*b^2*c^2+6*a*b*c*d))+a*e^-5*(e*x)^(5+m)*(5+m)^-1*(A*(a^2*d^2+3*b^2*c^2+6*a*b*c*d)+B*a*c*(2*a*d+3*b*c))+b*e^-9*(e*x)^(9+m)*(9+m)^-1*(c*b^2*(B*c+2*A*d)+3*B*a^2*d^2+3*a*b*d*(A*d+2*B*c))+A*a^3*c^2*e^-1*(e*x)^(1+m)*(1+m)^-1+B*b^3*d^2*e^-13*(e*x)^(13+m)*(13+m)^-1+c*a^2*e^-3*(e*x)^(3+m)*(3+m)^-1*(B*a*c+2*A*a*d+3*A*b*c)+d*b^2*e^-11*(e*x)^(11+m)*(11+m)^-1*(A*b*d+2*B*b*c+3*B*a*d))
@test integrate((e*x)^m*(a+b*x^2)^2*(c+d*x^2)^2*(A+B*x^2), x) == :(e^-7*(e*x)^(7+m)*(7+m)^-1*(B*a^2*d^2+c*b^2*(B*c+2*A*d)+2*a*b*d*(A*d+2*B*c))+e^-5*(e*x)^(5+m)*(5+m)^-1*(A*(a^2*d^2+b^2*c^2+4*a*b*c*d)+2*B*a*c*(a*d+b*c))+A*a^2*c^2*e^-1*(e*x)^(1+m)*(1+m)^-1+B*b^2*d^2*e^-11*(e*x)^(11+m)*(11+m)^-1+a*c*e^-3*(e*x)^(3+m)*(3+m)^-1*(2*A*(a*d+b*c)+B*a*c)+b*d*e^-9*(e*x)^(9+m)*(9+m)^-1*(A*b*d+2*B*a*d+2*B*b*c))
@test integrate((e*x)^m*(c+d*x^2)^2*(A+B*x^2)*(a+b*x^2), x) == :(e^-5*(e*x)^(5+m)*(5+m)^-1*(a*d*(A*d+2*B*c)+b*c*(B*c+2*A*d))+c*e^-3*(e*x)^(3+m)*(3+m)^-1*(A*b*c+B*a*c+2*A*a*d)+d*e^-7*(e*x)^(7+m)*(7+m)^-1*(A*b*d+B*a*d+2*B*b*c)+A*a*c^2*e^-1*(e*x)^(1+m)*(1+m)^-1+B*b*d^2*e^-9*(e*x)^(9+m)*(9+m)^-1)
@test integrate((e*x)^m*(c+d*x^2)^2*(A+B*x^2), x) == :(A*c^2*e^-1*(e*x)^(1+m)*(1+m)^-1+B*d^2*e^-7*(e*x)^(7+m)*(7+m)^-1+c*e^-3*(e*x)^(3+m)*(3+m)^-1*(B*c+2*A*d)+d*e^-5*(e*x)^(5+m)*(5+m)^-1*(A*d+2*B*c))
@test integrate((e*x)^m*(a+b*x^2)^3*(c+d*x^2)^3*(A+B*x^2), x) == :(e^-9*(e*x)^(9+m)*(9+m)^-1*(B*a^3*d^3+b^3*c^2*(B*c+3*A*d)+3*b*a^2*d^2*(A*d+3*B*c)+9*a*c*d*b^2*(A*d+B*c))+e^-7*(e*x)^(7+m)*(7+m)^-1*(A*(a^3*d^3+b^3*c^3+9*a*d*b^2*c^2+9*b*c*a^2*d^2)+3*B*a*c*(a^2*d^2+b^2*c^2+3*a*b*c*d))+A*a^3*c^3*e^-1*(e*x)^(1+m)*(1+m)^-1+B*b^3*d^3*e^-15*(e*x)^(15+m)*(15+m)^-1+a^2*c^2*e^-3*(e*x)^(3+m)*(3+m)^-1*(3*A*(a*d+b*c)+B*a*c)+b^2*d^2*e^-13*(e*x)^(13+m)*(13+m)^-1*(A*b*d+3*B*a*d+3*B*b*c)+3*a*c*e^-5*(e*x)^(5+m)*(5+m)^-1*(A*(a^2*d^2+b^2*c^2+3*a*b*c*d)+B*a*c*(a*d+b*c))+3*b*d*e^-11*(e*x)^(11+m)*(11+m)^-1*(B*a^2*d^2+c*b^2*(A*d+B*c)+a*b*d*(A*d+3*B*c)))
@test integrate((e*x)^m*(a+b*x^2)^2*(c+d*x^2)^3*(A+B*x^2), x) == :(e^-7*(e*x)^(7+m)*(7+m)^-1*(a^2*d^2*(A*d+3*B*c)+b^2*c^2*(B*c+3*A*d)+6*a*b*c*d*(A*d+B*c))+c*e^-5*(e*x)^(5+m)*(5+m)^-1*(A*(b^2*c^2+3*a^2*d^2+6*a*b*c*d)+B*a*c*(2*b*c+3*a*d))+d*e^-9*(e*x)^(9+m)*(9+m)^-1*(B*a^2*d^2+3*c*b^2*(A*d+B*c)+2*a*b*d*(A*d+3*B*c))+A*a^2*c^3*e^-1*(e*x)^(1+m)*(1+m)^-1+B*b^2*d^3*e^-13*(e*x)^(13+m)*(13+m)^-1+a*c^2*e^-3*(e*x)^(3+m)*(3+m)^-1*(B*a*c+2*A*b*c+3*A*a*d)+b*d^2*e^-11*(e*x)^(11+m)*(11+m)^-1*(A*b*d+2*B*a*d+3*B*b*c))
@test integrate((e*x)^m*(c+d*x^2)^3*(A+B*x^2)*(a+b*x^2), x) == :(c*e^-5*(e*x)^(5+m)*(5+m)^-1*(b*c*(B*c+3*A*d)+3*a*d*(A*d+B*c))+d*e^-7*(e*x)^(7+m)*(7+m)^-1*(a*d*(A*d+3*B*c)+3*b*c*(A*d+B*c))+c^2*e^-3*(e*x)^(3+m)*(3+m)^-1*(A*b*c+B*a*c+3*A*a*d)+d^2*e^-9*(e*x)^(9+m)*(9+m)^-1*(A*b*d+B*a*d+3*B*b*c)+A*a*c^3*e^-1*(e*x)^(1+m)*(1+m)^-1+B*b*d^3*e^-11*(e*x)^(11+m)*(11+m)^-1)
@test integrate((e*x)^m*(c+d*x^2)^3*(A+B*x^2), x) == :(A*c^3*e^-1*(e*x)^(1+m)*(1+m)^-1+B*d^3*e^-9*(e*x)^(9+m)*(9+m)^-1+c^2*e^-3*(e*x)^(3+m)*(3+m)^-1*(B*c+3*A*d)+d^2*e^-7*(e*x)^(7+m)*(7+m)^-1*(A*d+3*B*c)+3*c*d*e^-5*(e*x)^(5+m)*(5+m)^-1*(A*d+B*c))
@test integrate(x^-1*(a+b*x^2)^(1/2)*(A+B*x^2)*(c+d*x^2), x) == :(A*c*(a+b*x^2)^(1/2)+-1//15*b^-2*(a+b*x^2)^3//2*(-5*b*(A*d+B*c)+2*B*a*d+-3*B*b*d*x^2)+-1*A*c*a^(1/2)*arctanh(a^-1//2*(a+b*x^2)^(1/2)))
@test integrate(x^-1*(c+d*x^2)^(1/2)*(A+B*x^2)*(a+b*x^2), x) == :(A*a*(c+d*x^2)^(1/2)+-1//15*d^-2*(c+d*x^2)^3//2*(-1*d*(5*A*b+5*B*a)+2*B*b*c+-3*B*b*d*x^2)+-1*A*a*c^(1/2)*arctanh(c^-1//2*(c+d*x^2)^(1/2)))

max_stars_repo_path: test/1 Algebraic functions/1.1 Binomial products/1.1.2 Quadratic/1.1.2.6 (g x)^m (a+b x^2)^p (c+d x^2)^q (e+f x^2)^r.jl


text:
# Enter your code here
n = parse(Int, readline())
arr = parse.(Int,split(readline()))
numswaps = 0
for i = 1:n
    for j = 1:n-1
        if arr[j]>arr[j+1]
            dummy = arr[j]
            arr[j] = arr[j+1]
            arr[j+1] = dummy
            numswaps += 1
        end
    end
end
print("Array is sorted in $(numswaps) swaps.",'\n')
print("First Element: $(arr[1])",'\n')
print("Last Element: $(arr[end])",'\n')

max_stars_repo_path: Hackerrank/30 Days of Code/Julia/day 20.jl

max_stars_repo_name: Next-Gen-UI/Code-Dynamics

max_issues_repo_name: Next-Gen-UI/Code-Dynamics

max_forks_repo_name: Next-Gen-UI/Code-Dynamics


text:
using LassoPath
using Test
import UnicodePlots as UnicodePlots
using Plots
using LinearAlgebra

function Test1()
    """
        Test 1, some basic test on a random matrix, copied columns of matrix.
    """
    deg = 3
    N = 2000
    B = rand(N, deg)
    A = fill(0.0, (N, 2deg))
    A[:, 1:deg] = B
    A[:, deg + 1:end] = 2B.^2 + B + 0.1.*randn(size(B))

    b = A*ones(2deg, 1)
    instance = LassoSCOP(A, b)
    results, λs = GetLassoPath(instance)
    display(results)
    display(λs)

    # Let's plot this, and a global variable will be defined.
    global Test1Results = results

    ylim = (minimum(Test1Results), maximum(Test1Results))
    ThePlot = UnicodePlots.lineplot(log2.(λs), Test1Results[1, :], ylim=ylim, title="Lasso Path")

    for II ∈ 2:size(Test1Results, 1)
        UnicodePlots.lineplot!(ThePlot, log2.(λs), Test1Results[II, :])
    end
    display(ThePlot)
    ThePlot = UnicodePlots.heatmap(results[end:-1:begin, :])
    display(ThePlot)

    # Visualize the path.
    VisualizeLassoPath(instance)

    # Finding the important weights for it.
    PrintTitle("These are the important weights: ")
    WeightsIndices = CaptureImportantWeights(instance)
    println(WeightsIndices)
    return true

end

function Test2()
    """
        Testing the iterator for the lasso path object.
    """
    PrintTitle("Testing the Iterator for the Lasso Path type. ")
    deg = 3
    N = 2000
    B = rand(N, deg)
    A = fill(0.0, (N, 2deg))
    A[:, 1:deg] = B
    A[:, deg + 1:end] = 2B.^2 + B + 0.1.*randn(size(B))

    # Use this to test the lasso iterate.
    let A = A, b = A*ones(2deg, 1)
        l = LassoSCOP(A, b )
        for (λ, x) ∈ l
            println(λ)
            println(x)
        end
    end

    return true
end

function Test3()
    """
        Test proximal Gradient, internal methods
    """
    PrintTitle("Testing the proximal operator with simple regression. ")
    N, M = 10, 100
    A = randn(N, M)
    b = randn(N, 1)
    instance = LassoPath.BuildPG2NormLasso(A, b, 1/M)
    instance.tol= 1e-4

    sol = LassoPath.OptimizeProximalGradient(instance, A\b)
    error = norm(A*sol - b)
    println("The error of the solution is: $(error)")
    println("Non-zero entries: $(sum(sol .!= 0))")

    return error != Inf

end

function Test4()
    """
        Testing the proximal lasso path method.

    """
    """
    Test 1, some basic test on a random matrix, copied columns of matrix.
"""
    deg = 3
    N = 2000
    B = rand(N, deg)
    A = fill(0.0, (N, 2deg))
    A[:, 1:deg] = B
    A[:, deg + 1:end] = 2B.^2 + B + 0.1.*randn(size(B))

    b = A*ones(2deg, 1)
    instance = Lasso2NormProx(A, b)

    results, λs = GetLassoPath(instance)
    display(results)
    display(λs)

    # Let's plot this, and a global variable will be defined.
    global Test1Results = results

    ylim = (minimum(Test1Results), maximum(Test1Results))
    ThePlot = UnicodePlots.lineplot(log2.(λs), Test1Results[1, :], ylim=ylim, title="Lasso Path")

    for II ∈ 2:size(Test1Results, 1)
        UnicodePlots.lineplot!(ThePlot, log2.(λs), Test1Results[II, :])
    end
    display(ThePlot)
    ThePlot = UnicodePlots.heatmap(results[end:-1:begin, :])
    display(ThePlot)

    # Visualize the path.
    VisualizeLassoPath(instance)

    # Finding the important weights for it.
    PrintTitle("These are the important weights: ")
    WeightsIndices = CaptureImportantWeights(instance)
    println(WeightsIndices)

    return true
end

function Test5()
    """
        Test High dimensional Array for the proximal gradient lasso method.
    """

end

# A set of basic tests for Lasso.
@testset "LassoPath.jl" begin
    PrintTitle("Test Set 1")
    @test Test1()
    @test Test2()
    @test Test3()
    @test Test4()
end


text:
using BitBasis, LegibleLambdas

export AbstractRegister, @λ, @lambda

"""
    AbstractRegister{B}

Abstract type for quantum registers. `B` is the batch size.
"""
abstract type AbstractRegister{B} end

# properties
"""
    nactive(register) -> Int

Returns the number of active qubits.

!!! note

    Operators always apply on active qubits.
"""
@interface nactive(::AbstractRegister)

"""
    nqubits(register) -> Int

Returns the (total) number of qubits. See [`nactive`](@ref), [`nremain`](@ref)
for more details.
"""
@interface nqubits(::AbstractRegister)

"""
    nremain(register) -> Int

Returns the number of non-active qubits.
"""
@interface nremain(r::AbstractRegister) = nqubits(r) - nactive(r)

"""
    nbatch(register) -> Int

Returns the number of batches.
"""
@interface nbatch(r::AbstractRegister{B}) where B = B

# same with nbatch
Base.length(r::AbstractRegister{B}) where B = B

"""
    addbits!(register, n::Int) -> register
    addbits!(n::Int) -> λ(register)

Add `n` qubits to given register in state |0>.
i.e. |psi> -> |000> ⊗ |psi>, increased bits have higher indices.

If only an integer is provided, then returns a lambda function.
"""
@interface addbits!(::AbstractRegister, n::Int)

addbits!(n::Int) = @λ(register -> addbits!(register, n))

"""
    insert_qubits!(register, loc::Int; nqubits::Int=1) -> register
    insert_qubits!(loc::Int; nqubits::Int=1) -> λ(register)

Insert `n` qubits to given register in state |0>.
i.e. |psi> -> |psi> ⊗ |000> ⊗ |psi>, increased bits have higher indices.

If only an integer is provided, then returns a lambda function.
"""
@interface insert_qubits!(::AbstractRegister, loc::Int; nqubits::Int=1)
insert_qubits!(loc::Int; nqubits::Int=1) = @λ(register -> insert_qubits!(register, loc; nqubits=n))

"""
    focus!(register, locs) -> register

Focus the wires on specified location.

# Example

```julia
julia> focus!(r, (1, 2, 4))

```
"""
@interface focus!(r::AbstractRegister, locs)

"""
    focus!(locs...) -> f(register) -> register

Lazy version of [`focus!`](@ref), this returns a lambda which requires a register.
"""
focus!(locs::Int...) = focus!(locs)
focus!(locs::NTuple{N, Int}) where N = @λ(register -> focus!(register, locs))
focus!(locs::UnitRange) = @λ(register -> focus!(register, locs))
"""
    focus(f, register, locs...)

Call a callable `f` under the context of `focus`. See also [`focus!`](@ref).

# Example

print the focused register

```julia
julia> r = ArrayReg(bit"101100")
ArrayReg{1,Complex{Float64},Array...}
    active qubits: 6/6

julia> focus(x->(println(x);x), r, 1, 2);
ArrayReg{1,Complex{Float64},Array...}
    active qubits: 2/6
```
"""
@interface focus!(f::Base.Callable, r::AbstractRegister, locs::Int...) = focus(f, r, locs)

focus!(f::Base.Callable, r::AbstractRegister, loc::Int) = focus(f, r, (loc, ))
focus!(f::Base.Callable, r::AbstractRegister, locs) =
    relax!(f(focus!(r, locs)), locs; to_nactive=nqubits(r))

"""
    relax!(register[, locs]; to_nactive=nqubits(register)) -> register

Inverse transformation of [`focus!`](@ref), where `to_nactive` is the number
 of active bits for target register.
"""
@interface relax!(r::AbstractRegister, locs; to_nactive::Int=nqubits(r))
relax!(r::AbstractRegister; to_nactive::Int=nqubits(r)) = relax!(r, (); to_nactive=to_nactive)

"""
    relax!(locs::Int...; to_nactive=nqubits(register)) -> f(register) -> register

Lazy version of [`relax!`](@ref), it will be evaluated once you feed a register
to its output lambda.
"""
relax!(locs::Int...; to_nactive::Union{Nothing, Int}=nothing) =
    relax!(locs; to_nactive=to_nactive)

function relax!(locs::NTuple{N, Int}; to_nactive::Union{Nothing, Int}=nothing) where N
    lambda = function (r::AbstractRegister)
        if to_nactive === nothing
            return relax!(r, locs; to_nactive=nqubits(r))
        else
            return relax!(r, locs; to_nactive=to_nactive)
        end
    end

    @static if VERSION < v"1.1.0"
        return LegibleLambda(
            "(register->relax!(register, locs...; to_nactive))",
            lambda
            )
    else
        return LegibleLambda(
                lambda,
                :(register->relax!(register, locs...; to_nactive)),
                Dict(:locs=>locs, :to_nactive=>to_nactive)
        )
    end
end

## Measurement

export ComputationalBasis
struct ComputationalBasis end

export AllLocs
struct AllLocs end

export measure, measure!, measure_remove!, measure_collapseto!

"""
    measure([rng,] register[, operator][, locs]; nshots=1) -> Vector{Int}

Return measurement results of current active qubits (regarding to active qubits,
see [`focus!`](@ref) and [`relax!`](@ref)).
"""
function measure end

"""
    measure!([rng,] [operator, ]register[, locs])

Measure current active qubits or qubits at `locs` and collapse to result state.
"""
function measure! end

"""
    measure_remove!([rng,] [operator, ]reg::AbstractRegister[, locs])

Measure current active qubits or qubits at `locs` and remove them.
"""
function measure_remove! end

"""
    measure_collapseto!([rng,] [operator, ]reg::AbstractRegister[, locs]; config) -> Int

Measure current active qubits or qubits at `locs` and set the register to specific value.
"""
function measure_collapseto! end

# focus context
for FUNC in [:measure!, :measure_collapseto!, :measure_remove!, :measure]
    rotback = FUNC == :measure! ? :(reg.state = V*reg.state) : :()
    @eval function $FUNC(rng::AbstractRNG, op::Eigen, reg::AbstractRegister, locs::AllLocs; kwargs...)
        E, V = op
        reg.state = V'*reg.state
        res = $FUNC(rng, ComputationalBasis(), reg, locs; kwargs...)
        $rotback
        E[Int64.(res) .+ 1]
    end
    @eval $FUNC(rng::AbstractRNG, op, reg::AbstractRegister; kwargs...) = $FUNC(rng, op, reg, AllLocs(); kwargs...)
    @eval $FUNC(rng::AbstractRNG, reg::AbstractRegister, locs; kwargs...) = $FUNC(rng, ComputationalBasis(), reg, locs; kwargs...)
    @eval $FUNC(rng::AbstractRNG, reg::AbstractRegister; kwargs...) = $FUNC(rng, ComputationalBasis(), reg, AllLocs(); kwargs...)
    @eval $FUNC(args...; kwargs...) = $FUNC(Random.GLOBAL_RNG, args...; kwargs...)
end

for FUNC in [:measure_collapseto!, :measure!, :measure]
    @eval function $FUNC(rng::AbstractRNG, op, reg::AbstractRegister, locs::Union{Tuple, Vector, Integer, UnitRange}; kwargs...)
        nbit = nactive(reg)
        focus!(reg, locs)
        res = $FUNC(rng, op, reg, AllLocs(); kwargs...)
        relax!(reg, locs; to_nactive=nbit)
        res
    end
end

function measure_remove!(rng::AbstractRNG, op, reg::AbstractRegister, locs)
    nbit = nactive(reg)
    focus!(reg, locs)
    res = measure_remove!(rng, op, reg, AllLocs())
    relax!(reg; to_nactive=nbit-length(locs))
    res
end

"""
    select!(dest::AbstractRegister, src::AbstractRegister, bits::Integer...) -> AbstractRegister
    select!(register::AbstractRegister, bits::Integer...) -> register

select a subspace of given quantum state based on input eigen state `bits`.
See also [`select`](@ref).

## Example

`select!(reg, 0b110)` will select the subspace with (focused) configuration `110`.
After selection, the focused qubit space is 0, so you may want call `relax!` manually.

!!! tip

    Developers should overload `select!(r::RegisterType, bits::NTuple{N, <:Integer})` and
    do not assume `bits` has specific number of bits (e.g `Int64`), or it will restrict the
    its maximum available number of qubits.
"""
@interface select!(r::AbstractRegister, bits)

"""
    select!(b::Integer) -> f(register)

Lazy version of [`select!`](@ref). See also [`select`](@ref).
"""
select!(bits...) = @λ(register->select!(register, bits...))

"""
    partial_tr(register, locs)

Return a register which is the partial traced on `locs`.
"""
@interface partial_tr(r::AbstractRegister, locs)

"""
    partial_tr(locs) -> f(register)

Curried version of `partial_tr(register, locs)`.
"""
partial_tr(locs) = @λ(register->partial_tr(register, locs))

"""
    select(register, bits) -> AbstractRegister

Non-inplace version of [`select!`](@ref).
"""
@interface select(register::AbstractRegister, bits)

"""
    join(::AbstractRegister...) -> register

Merge several registers as one register via tensor product.
"""
@interface Base.join(::AbstractRegister...)

"""
    repeat(r::AbstractRegister, n::Int) -> register

Repeat register `r` for `n` times on batch dimension.
"""
@interface Base.repeat(::AbstractRegister, n::Int)

"""
    basis(register) -> UnitRange

Returns an `UnitRange` of the all the bits in the Hilbert space of given register.
"""
@interface BitBasis.basis(r::AbstractRegister) = basis(nqubits(r))

"""
    probs(register)

Returns the probability distribution of computation basis, aka ``|<x|ψ>|^2``.
"""
@interface probs(r::AbstractRegister)

"""
    reorder!(reigster, orders)

Reorder the locations of register by input orders.
"""
@interface reorder!(r::AbstractRegister, orders)

"""
    invorder(register)

Inverse the locations of register.
"""
@interface invorder!(r::AbstractRegister) = reorder!(r, Tuple(nactive(r):-1:1))

"""
    collapseto!(register, config)

Set the `register` to bit string literal `bit_str` (or an equivalent integer). About bit string literal,
see more in [`@bit_str`](@ref).
"""
collapseto!(r::AbstractRegister, config::Integer=0) = collapseto!(r, Int64(config))
@interface collapseto!(r::AbstractRegister, config::BitStr) = collapseto!(r, config)

"""
    fidelity(register1, register2)

Return the fidelity between two states.

# Definition
The fidelity of two quantum state for qubits is defined as:

```math
F(ρ, σ) = tr(\\sqrt{\\sqrt{ρ}σ\\sqrt{ρ}})
```

Or its equivalent form (which we use in numerical calculation):

```math
F(ρ, σ) = sqrt(tr(ρσ) + 2 \\sqrt{det(ρ)det(σ)})
```

# Reference

- Jozsa R. Fidelity for mixed quantum states[J]. Journal of modern optics, 1994, 41(12): 2315-2323.
- Nielsen M A, Chuang I. Quantum computation and quantum information[J]. 2002.

!!! note

    The original definition of fidelity ``F`` was from "transition probability",
    defined by Jozsa in 1994, it is the square of what we use here.
"""
@interface fidelity(r1::AbstractRegister, r2::AbstractRegister)

"""
    tracedist(register1, register2)

Return the trace distance of `register1` and `register2`.

# Definition
Trace distance is defined as following:

```math
\\frac{1}{2} || A - B ||_{tr}
```

# Reference

- https://en.wikipedia.org/wiki/Trace_distance
"""
@interface tracedist(r1::AbstractRegister, r2::AbstractRegister)

"""
    density_matrix(register)

Returns the density matrix of current active qubits.
"""
@interface density_matrix(::AbstractRegister)

"""
    ρ(register)

Returns the density matrix of current active qubits. This is the same as
[`density_matrix`](@ref).
"""
@interface ρ(x) = density_matrix(x)

"""
    viewbatch(register, i::Int) -> AbstractRegister{1}

Returns a view of the i-th slice on batch dimension.
"""
@interface viewbatch(::AbstractRegister, ::Int)

function Base.iterate(it::AbstractRegister{B}, state=1) where B
    if state > B
        return nothing
    else
        return viewbatch(it, state), state + 1
    end
end

# fallback printing
function Base.show(io::IO, reg::AbstractRegister)
    summary(io, reg)
    print(io, "\n    active qubits: ", nactive(reg), "/", nqubits(reg))
end


text:
using ConcreteStructs
using DynamicIterators
import DynamicIterators: dub, dyniterate, evolve
using Base.Iterators: SizeUnknown, IsInfinite

import MeasureBase: For

export Chain

@concrete terse struct Chain{K,M} <: AbstractMeasure
    κ::K
    μ::M
end

Pretty.quoteof(c::Chain) = :(Chain($(Pretty.quoteof(c.κ)), $(Pretty.quoteof(c.μ))))

Base.length(::Chain) = ∞

function basemeasure(mc::Chain)
    Chain(basemeasure ∘ mc.κ, basemeasure(mc.μ))
end

Base.IteratorEltype(mc::Chain) = Base.HasEltype()

Base.eltype(::Type{C}) where {K,M,C<:Chain{K,M}} = eltype(M)

function logdensity(mc::Chain, x)
    μ = mc.μ
    ℓ = 0.0
    for xj in x
        ℓ += logdensity(μ, xj)
        μ = mc.κ(xj)
    end
    return ℓ
end

DynamicIterators.evolve(mc::Chain, μ) =  μ ⋅ mc.κ
DynamicIterators.evolve(mc::Chain) =  mc.μ

dyniterate(E::Chain, value) = dub(evolve(E, value))
dyniterate(E::Chain, ::Nothing) = dub(evolve(E))
Base.iterate(E::Chain) = dyniterate(E, nothing)
Base.iterate(E::Chain, value) = dyniterate(E, value)

function DynamicIterators.dyniterate(r::Chain, (x,rng)::Sample)
    μ = r.κ(x)
    y = rand(rng, μ)
    return (μ ↝ y), Sample(y, rng)
end

Base.IteratorSize(::Chain) = IsInfinite()
Base.IteratorSize(::Type{Chain}) = IsInfinite()

function Base.rand(rng::AbstractRNG, T::Type, chain::Chain)
    r = ResettableRNG(rng)
    return RealizedSamples(r, chain)
end

###############################################################################
# DynamicFor

# A `DynamicFor` is produced when `For` is called on a `DynamicIterator`.

@concrete terse struct DynamicFor{T,K,S} <: AbstractMeasure
    κ ::K
    iter :: S
end

Pretty.quoteof(r::DynamicFor) = :(DynamicFor($(Pretty.quoteof(r.κ)), $(Pretty.quoteof(r.iter))))

function DynamicFor(κ::K,iter::S) where {K,S}
    T = typeof(κ(first(iter)))
    DynamicFor{T,K,S}(κ,iter)
end

function Base.rand(rng::AbstractRNG, T::Type, df::DynamicFor)
    r = ResettableRNG(rng)
    return RealizedSamples(r, df)
end

function logdensity(df::DynamicFor, y)
    ℓ = 0.0
    for (xj, yj) in zip(df.iter, y)
        ℓ += logdensity(df.κ(xj), yj)
    end
    return ℓ
end

Base.eltype(::Type{D}) where {T,D<:DynamicFor{T}} = eltype(T)

Base.IteratorEltype(d::DynamicFor) = Base.HasEltype()

Base.IteratorSize(d::DynamicFor) = Base.IteratorSize(d.iter)

function Base.iterate(d::DynamicFor)
    (x,s) = iterate(d.iter)
    (d.κ(x), s)
end

function Base.iterate(d::DynamicFor, s)
    (x,s) = iterate(d.iter, s)
    (d.κ(x), s)
end

Base.length(d::DynamicFor) = length(d.iter)

For(f, r::Realized) = DynamicFor(f, copy(r))

function Base.rand(rng::AbstractRNG, dfor::DynamicFor)
    r = ResettableRNG(rng)
    return RealizedSamples(r, dfor)
end

function dyniterate(df::DynamicFor, st, args...)
    (val, state) = dyniterate(df.iter, st, args...)
    return (df.κ(val), state)
end

For(f, it::DynamicIterator) = DynamicFor(f, it)

For(f, it::DynamicFor) = DynamicFor(f, it)

function dyniterate(df::DynamicFor, state)
      ϕ = dyniterate(df.iter, state)
      ϕ === nothing && return nothing
      u, state = ϕ
      df.f(u), state
end

function Base.collect(r::Realized)
    next = iterate(r)
    isnothing(next) && return []
    (x,s) = next
    a = similar(r.iter, typeof(x))

    i = 1
    @inbounds a[i] = x
    while !isnothing(next)
        (x, s) = next
        @inbounds a[i] = x
        i += 1
        next = iterate(r, s)
    end
    return a
end

function testvalue(mc::Chain)
    μ = mc.μ
    κ = mc.κ
    rand(Chain(Dirac ∘ testvalue ∘ κ, (Dirac ∘ testvalue)(μ)))
end


text:
using ChebyshevQuantum
using Test

tol_val = 1e-10

@testset "test particle in box " begin

    H = ChebyshevQuantum.SolveEig.setup_ham(x->0.0, 40);
    #    vals, vects =  ChebyshevQuantum.SolveEig.solve(H)
    vals, vects = eigen( H[2:end-1,2:end-1]);

    @test isapprox(vals[1],   pi^2/8, atol=tol_val)
    @test isapprox(vals[2], 2^2*pi^2/8, atol=tol_val)
    @test isapprox(vals[3], 3^2*pi^2/8, atol=tol_val)
    @test isapprox(vals[4], 4^2*pi^2/8, atol=tol_val)
    @test isapprox(vals[5], 5^2*pi^2/8, atol=tol_val)

end

@testset "test particle in bigger box " begin

    H = ChebyshevQuantum.SolveEig.setup_ham(x->0.0, 40, a = 0, b = 10);
#    vals, vects =  ChebyshevQuantum.SolveEig.solve(H)
    vals, vects = eigen( H[2:end-1,2:end-1]);

    f = (2/10)^2

    @test isapprox(vals[1],   pi^2/8 * f, atol=tol_val)
    @test isapprox(vals[2], 2^2*pi^2/8 * f, atol=tol_val)
    @test isapprox(vals[3], 3^2*pi^2/8 * f, atol=tol_val)
    @test isapprox(vals[4], 4^2*pi^2/8 * f, atol=tol_val)
    @test isapprox(vals[5], 5^2*pi^2/8 * f, atol=tol_val)

end

tol_val = 1e-10

@testset "test particle in box N=100 " begin

    H = ChebyshevQuantum.SolveEig.setup_ham(x->0.0, 100);
#    vals, vects =  ChebyshevQuantum.SolveEig.solve(H)
    vals, vects = eigen( H[2:end-1,2:end-1]);

    @test isapprox(vals[1],   pi^2/8, atol=tol_val)
    @test isapprox(vals[2], 2^2*pi^2/8, atol=tol_val)
    @test isapprox(vals[3], 3^2*pi^2/8, atol=tol_val)
    @test isapprox(vals[4], 4^2*pi^2/8, atol=tol_val)
    @test isapprox(vals[5], 5^2*pi^2/8, atol=tol_val)

end

tol_val = 1e-9

@testset "test particle in finite box " begin

    #test edge detection

    function VV(x)
        if abs(x) >= 1.0
            return 10.0
        else
            return 0.0
        end
    end

#    H, S = ChebyshevQuantum.SolveEig.total_ham_piecewise(VV, 30, [-2.0, -1.0, 1.0, 2.0], d2=-0.5)
    #    H, S, lhs = ChebyshevQuantum.Op.setup(A2=-0.5, A0=VV,a=-2.0,b=2.0, N=30)

    vals, vects = ChebyshevQuantum.SolveEig.solve(V=VV, a=-2.0, b = 2.0)

    #    vals  =  eigvals(H, S)
#    vals = real.(vals)

    @test isapprox(vals[1],  0.819848145939759, atol=tol_val) #reference values from chebfun.m (see below)
    @test isapprox(vals[2],  3.222122379353262, atol=tol_val)

end

"""
matlab code using chebfun
tic
x = chebfun('x', [-2 2]);
V = 10.0*(abs(x)>=1);
L = chebop(-2,2);
L.op = @(x,u) -0.5*diff(u,2) + V*u;
L.bc = 0;
neigs = 12;
[EV,D] = eigs(L,neigs);
disp(diag(D))
toc
"""

nothing


text:
bitsof(::Type{T}) where {T<:Union{Signed, Unsigned}}  = sizeof(T) << 3
bitsof(x::T) where {T<:Union{Signed, Unsigned}} = bitsof(T)

for S in (:SafeSigned, :SafeUnsigned)
  @eval begin
    zero(::Type{T}) where T<:$S = safeint(zero(baseint(T)))
    one(::Type{T})  where T<:$S = safeint(one(baseint(T)))
    sizeof(::Type{T}) where T<:$S  = safeint(sizeof(baseint(T)))
    bitsof(::Type{T}) where T<:$S  = safeint(sizeof(baseint(T)) << 3)
    @inline zero(x::T) where T<:$S = zero(T)
    @inline one(x::T)  where T<:$S = one(T)
    @inline sizeof(x::T) where T<:$S = sizeof(T)
    @inline bitsof(x::T) where T<:$S = bitsof(T)
    ndigits0z(n::T) where T<:$S = ndigits0z(baseint(x)) # do not reconvert
    ndigits0z(n::T1, b::T2) where T1<:$S where T2<:SafeInteger = ndigits0z(baseint(x), baseint(b)) # do not reconvert
    ndigits0z(n::T1, b::T2) where T1<:$S where T2<:Integer = ndigits0z(baseint(x), b) # do not reconvert
    @inline count_zeros(x::T) where T<:$S = safeint(count_zeros(baseint(x)))
    @inline count_ones(x::T) where T<:$S = safeint(count_ones(baseint(x)))
    @inline leading_zeros(x::T) where T<:$S = safeint(leading_zeros(baseint(x)))
    @inline trailing_zeros(x::T) where T<:$S = safeint(trailing_zeros(baseint(x)))
    @inline leading_ones(x::T) where T<:$S = safeint(leading_ones(baseint(x)))
    @inline trailing_ones(x::T) where T<:$S = safeint(trailing_ones(baseint(x)))
  end
end

@inline signbit(x::T) where T<:SafeSigned = signbit(baseint(x))
@inline signbit(x::T) where T<:SafeUnsigned = false
@inline sign(x::T) where T<:SafeSigned = safeint(sign(baseint(x)))
@inline sign(x::T) where T<:SafeUnsigned = one(T)

function abs(x::T) where T<:SafeSigned
  x === typemin(x) && throw(OverflowError("cannot take `abs(typemin)`"))
  return safeint(abs(baseint(x)))
end
abs(x::T) where T<:SafeUnsigned = x

@inline function Base.:(~)(x::T) where T<:SafeSigned
    return safeint(~(baseint(x)))
end

@inline function Base.:(~)(x::T) where T<:SafeUnsigned
    return safeint(~(baseint(x)))
end

@inline function Base.:(-)(x::T) where T<:SafeSigned
  x === typemin(x) && throw(OverflowError("cannot negate typemin"))
  return safeint(-(baseint(x)))
end

@inline function Base.:(-)(x::T) where T<:SafeUnsigned
    throw(OverflowError("cannot negate unsigned type"))
end

@inline function copysign(x::T, y::T) where T<:SafeSigned
  return safeint(baseint( signbit(y) ? -abs(x) : abs(x) ))
end
@inline copysign(x::T, y::T) where T<:SafeUnsigned = x

@inline function flipsign(x::T, y::T) where T<:SafeSigned
  return safeint(baseint( signbit(y) ? -(x) : x ))
end
@inline flipsign(x::T, y::T) where T<:SafeUnsigned = x

typemin(::Type{T}) where {T<:SafeInteger} = safeint(typemin(baseint(T)))
typemax(::Type{T}) where {T<:SafeInteger} = safeint(typemax(baseint(T)))

widen(::Type{T}) where {T<:SafeInteger} = safeint(widen(baseint(T)))

# abs2

sqrtmax(::Type{T}) where {T<:Integer} = floor(T, sqrt(typemax(T)))

abs2max(::Type{T}) where {T<:Signed} = sqrtmax(T)

function abs2max(::Type{T}) where {T<:Unsigned}
    result = sqrtmax(T)
    if T(result*result) === zero(T)
       result -= one(T)
    end
    return result
end

for T in (:Int8, :Int16, :Int32, :Int64, :Int128,
          :UInt8, :UInt16, :UInt32, :UInt64, :UInt128)
  @eval maxabs2(::Type{$T}) = abs2max($T)
end

for T in (:SafeInt8, :SafeInt16, :SafeInt32, :SafeInt64, :SafeInt128,
          :SafeUInt8, :SafeUInt16, :SafeUInt32, :SafeUInt64, :SafeUInt128)
  @eval maxabs2(::Type{$T}) = abs2max(baseint($T))
end

function abs2(x::T) where T<:SafeSigned
  x > maxabs2(T) && throw(OverflowError("$x^2 exceeds $T"))
  return x * x
end

function abs2(x::T) where T<:SafeUnsigned
  x > maxabs2(T) && throw(OverflowError("$x^2 exceeds $T"))
  return x * x
end


text:
using Statistics
using ImageTransformations
using CoordinateTransformations
using Interpolations

###############################################################################
# Stacking routines

"""
    collapse(cube; method=median, fill=0, degree=Linear())
    collapse(cube, angles; method=:deweight, fill=0, degree=Linear())

Combine all the frames of a cube using `method`. If `angles` are provided, will use [`derotate`](@ref) before combining.

If `method` is `:deweight`, the method of _Bottom et al. 2017_ will be used in which the combined image will be the derotated weighted sum of the frames weighted by the temporal variance. Keyword arguments will be passed to [`derotate`](@ref).

### References
1. [Bottom et al. 2017 "Noise-weighted Angular Differential Imaging"](https://ui.adsabs.harvard.edu/abs/2017RNAAS...1...30B)

# Examples
```jldoctest
julia> X = ones(3, 3, 2);

julia> collapse(X)
3×3 Matrix{Float64}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

julia> collapse(X, [0, 90])
3×3 Matrix{Float64}:
 1.0  1.0  0.5
 1.0  1.0  1.0
 0.5  0.5  0.5

julia> collapse(X, [0, 90], fill=NaN)
3×3 Matrix{Float64}:
   1.0    1.0  NaN
   1.0    1.0    1.0
 NaN    NaN    NaN

```

# See Also
[`collapse!`](@ref)
"""
collapse(cube::AbstractArray{T,3}; method=median) where {T} = method(cube, dims = 3)[:, :, 1]

collapse(cube::AbstractArray{T,3}, angles::AbstractVector; method=:deweight, kwargs...) where T =
    method === :deweight ? _collapse_deweighted(cube, angles; kwargs...) :
               collapse(derotate(cube, angles; kwargs...); method = method)

"""
    collapse!(cube, angles; method=:deweight, fill=0)

An in-place version of the derotating [`collapse`](@ref). The only difference is in this version the cube will be derotated in-place.
"""
collapse!(cube::AbstractArray{T,3}, angles::AbstractVector; method=:deweight, kwargs...) where T =
    method === :deweight ? _collapse_deweighted!(cube, angles; kwargs...) :
               collapse(derotate!(cube, angles; kwargs...); method = method)

# deweight using Bottom et al. 2017
 _collapse_deweighted(cube::AbstractArray{T,3}, angles::AbstractVector; kwargs...) where T =
    _collapse_deweighted!(deepcopy(cube), angles; kwargs...)

function _collapse_deweighted!(cube::AbstractArray{T,3}, angles::AbstractVector; fill=zero(T), kwargs...) where T
    varframe = var(cube, dims=3)[:, :, 1]

    # have to check if no variance otherwise the returns will be NaN
    all(p -> p ≈ 0, varframe) && return mean(derotate!(cube, angles; fill=fill, kwargs...); dims=3)[:, :, 1]
    # create a cube from the variance of each pixel across time
    varcube = similar(cube)
    # derotate both cubes
    Threads.@threads for idx in axes(varcube, 3)
        cube[:, :, idx] .= derotate(view(cube, :, :, idx), angles[idx]; fill=fill, kwargs...)
        varcube[:, :, idx] .= derotate(varframe, angles[idx]; fill=fill, kwargs...)
    end

    # calculate weighted sum and replace NaNs with our fill value
    out = sum(cube ./ varcube, dims=3) ./ sum(inv.(varcube), dims=3)
    @. out[isnan(out)] = fill
    return out[:, :, 1]
end

"""
    flatten(cube)

Given a cube of size `(x, y, n)` returns a matrix with size `(x * y, n)` where each row is a flattened image from the cube.

# Examples
```jldoctest
julia> X = ones(2, 2, 3);

julia> flatten(X)
4×3 Matrix{Float64}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
```

# See Also
[`expand`](@ref)
"""
flatten(cube::AbstractArray{T,3}) where T = reshape(cube, size(cube, 1) * size(cube, 2), size(cube, 3))
flatten(mat::AbstractMatrix) = mat

"""
    expand(matrix)

Given a matrix of size `(z, n)`, returns a cube of size `(x, x, n)` where `x=√z`.

Will throw an error if `z` is not a perfect square.

# Examples
```jldoctest
julia> X = ones(4, 3)
4×3 Matrix{Float64}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

julia> expand(X)[:, :, 1]
2×2 Matrix{Float64}:
 1.0  1.0
 1.0  1.0
```

# See Also
[`flatten`](@ref)
"""
function expand(mat::AbstractMatrix)
    z, n = size(mat)
    x = sqrt(z)
    isinteger(x) || error("Array of size $((x, x, n)) is not compatible with input matrix of size $(size(mat)).")
    return reshape(mat, Int(x), Int(x), n)
end
expand(cube::AbstractArray{T,3}) where {T} = cube

###############################################################################
# Rotation routines

"""
    derotate!(cube, angles; fill=0, degree=Linear())

In-place version of [`derotate`](@ref) which modifies `cube`.
"""
function derotate!(cube::AbstractArray{T,3},
                   angles::AbstractVector;
                   fill=zero(T),
                   degree=Linear()) where T
    Threads.@threads for i in axes(cube, 3)
        frame = @view cube[:, :, i]
        frame .= imrotate(frame, -deg2rad(angles[i]), axes(frame); method=degree, fillvalue=fill)
    end
    return cube
end

"""
    derotate(frame, angle; fill=0, degree=Linear())

Rotates `frame` counter-clockwise by `angle`, given in degrees. This is merely a convenient wrapper around `ImageTransformations.imrotate`.
"""
function derotate(frame::AbstractMatrix{T},
                   angle;
                   fill=zero(T),
                   degree=Linear()) where T
    return imrotate(frame, -deg2rad(angle), axes(frame); method=degree, fillvalue=fill)
end

"""
    derotate(cube, angles; fill=0, degree=Linear())

Rotates an array using the given angles in degrees.

This will rotate frame `i` counter-clockwise. Any values outside the original axes will be replaced with `fill`. If the given angles are true parallactic angles, the resultant cube will have each frame aligned with top pointing North. `degree` is the corresponding `Interpolations.Degree` for the B-Spline used to subsample the pixel values.

# Examples
```jldoctest
julia> X = zeros(3, 3, 1); X[2, 1, 1] = 1;

julia> X[:, :, 1]
3×3 Matrix{Float64}:
 0.0  0.0  0.0
 1.0  0.0  0.0
 0.0  0.0  0.0

julia> derotate(X, [90])[:, :, 1]
3×3 Matrix{Float64}:
 0.0  0.0          0.0
 0.0  4.44089e-16  0.0
 0.0  0.0          0.0
```

# See Also
[`derotate!`](@ref)
"""
function derotate(cube::AbstractArray{T,3}, angles::AbstractVector; kwargs...) where T
    all(iszero, angles) && return cube
    return derotate!(copy(cube), angles; kwargs...)
end

#################################

"""
    shift_frame(frame, dx, dy; fill=0)
    shift_frame(frame, dpos; fill=0)

Shifts `frame` by `dx` and `dy` with bilinear interpolation. If necessary, empty indices will be filled with `fill`.

# Examples
```jldoctest
julia> shift_frame([0 0 0; 0 1 0; 0 0 0], 1, -1)
3×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 1.0  0.0  0.0

julia> shift_frame(ans, (-1, 1), fill=NaN)
3×3 Matrix{Float64}:
 NaN    0.0    0.0
 NaN    1.0    0.0
 NaN  NaN    NaN
```
"""
function shift_frame(frame::AbstractMatrix{T}, dx, dy; fill=zero(T)) where T
    tform = Translation(-dx, -dy)
    return warp(frame, tform, axes(frame); fillvalue=fill)
end
shift_frame(frame::AbstractMatrix{T}, dpos; fill=zero(T)) where T = shift_frame(frame, dpos...; fill=fill)

"""
    shift_frame(cube, dx, dy; fill=0)
    shift_frame(cube, dpos; fill=0)

Shift each frame of `cube` by `dx` and `dy`, which can be integers or vectors. The change in position can be given as a tuple, which can also be put into a vector to use across the cube. If a frame is shifted outside its axes, the empty indices will be filled with `fill`.

# See Also
[`shift_frame!`](@ref)
"""
shift_frame(cube::AbstractArray{T,3}, dx, dy; fill=zero(T)) where T = shift_frame!(copy(cube), dx, dy; fill=fill)
shift_frame(cube::AbstractArray{T,3}, dpos; fill=zero(T)) where T = shift_frame!(copy(cube), dpos; fill=fill)

"""
    shift_frame!(cube, dx, dy; fill=0)
    shift_frame!(cube, dpos; fill=0)

In-place version of [`shift_frame`](@ref) which modifies `cube`.
"""
function shift_frame!(cube::AbstractArray{T,3}, dx::Number, dy::Number; fill=zero(T)) where T
    @inbounds for idx in axes(cube, 3)
        frame = @view cube[:, :, idx]
        tform = Translation(-dx, -dy)
        frame .= warp(frame, tform, axes(frame); fillvalue=fill)
    end
    return cube
end

shift_frame!(cube::AbstractArray{T,3}, dpos::Tuple; fill=zero(T)) where T = shift_frame!(cube, dpos...; fill=fill)

function shift_frame!(cube::AbstractArray{T,3}, dx::AbstractVector, dy::AbstractVector; fill=zero(T)) where T
    @inbounds for idx in axes(cube, 3)
        frame = @view cube[:, :, idx]
        tform = Translation(-dx[idx], -dy[idx])
        frame .= warp(frame, tform, axes(frame); fillvalue=fill)
    end
    return cube
end

function shift_frame!(cube::AbstractArray{T,3}, dpos::AbstractVector{<:Tuple}; fill=zero(T)) where T
    @inbounds for idx in axes(cube, 3)
        frame = @view cube[:, :, idx]
        dx, dy = dpos[idx]
        tform = Translation(-dx, -dy)
        frame .= warp(frame, tform, axes(frame); fillvalue=fill)
    end
    return cube
end

#################################

"""
    crop(frame, size; center=center(frame), force=false)
    crop(cube, size; center=center(cube), force=false)

Crop a frame or cube to `size`. `size` can be a tuple or an integer, which will make a square crop. The indices will be relative to `center`. To avoid uneven (odd) cropping, we may change `size`. To disable this behavior, set `force` to `true`. To avoid allocations, consider [`cropview`](@ref).
"""
crop(input, size; kwargs...) = collect(cropview(input, size; kwargs...))

cropview(input, size::Integer; kwargs...) = cropview(input, (size, size); kwargs...)

"""
    cropview(cube::AbstractArray{T, 3}, size; center=center(cube), force=false)

Crop a frame to `size`, returning a view of the frame. `size` can be a tuple or an integer, which will make a square crop. The indices will be relative to `center`. To avoid uneven (odd) cropping, we may change `size`. To disable this behavior, set `force` to `true`.

# See Also
[`crop`](@ref)
"""
function cropview(cube::AbstractArray{T,3}, size::Tuple; center=center(cube)[1:2], force=false, verbose=true) where T
    frame_size = (Base.size(cube, 1), Base.size(cube, 2))
    out_size = force ? size : check_size(frame_size, size)
    out_size != size && verbose && @info "adjusted size to $out_size to avoid uneven (odd) cropping"
    wing = @. (out_size - 1) / 2
    _init = @. floor(Int, center - wing)
    _end = @. floor(Int, center + wing)
    return view(cube, _init[1]:_end[1], _init[2]:_end[2], :)
end

"""
    cropview(frame::AbstractMatrix, size; center=center(frame), force=false)

Crop a frame to `size`, returning a view of the frame. `size` can be a tuple or an integer, which will make a square crop. The indices will be relative to `center`. To avoid uneven (odd) cropping, we may change `size`. To disable this behavior, set `force` to `true`.

# See Also
[`crop`](@ref)
"""
function cropview(frame::AbstractMatrix, size::Tuple; center=center(frame), force=false, verbose=true)
    out_size = force ? size : check_size(Base.size(frame), size)
    out_size != size && verbose && @info "adjusted size to $out_size to avoid uneven (odd) cropping"
    wing = @. (out_size - 1) / 2
    _init = @. floor(Int, center - wing)
    _end = @. floor(Int, center + wing)
    return view(frame, _init[1]:_end[1], _init[2]:_end[2])
end

"""
    HCIToolbox.check_size(frame_size, crop_size)

Given two image shapes, will adjust the output size to make sure even amounts are clipped off each side.
"""
check_size(frame_size, crop_size) = map((f, c) -> iseven(f - c) ? c : c + 1, frame_size, crop_size)


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

"""
    AbstractSolver

A solver for geostatistical problems.
"""
abstract type AbstractSolver end

"""
    AbstractEstimationSolver

A solver for a geostatistical estimation problem.
"""
abstract type AbstractEstimationSolver <: AbstractSolver end

"""
    AbstractSimulationSolver

A solver for a geostatistical simulation problem.
"""
abstract type AbstractSimulationSolver <: AbstractSolver end

"""
    AbstractLearningSolver

A solver for a geostatistical learning problem.
"""
abstract type AbstractLearningSolver <: AbstractSolver end

"""
    solve(problem, solver)

Solve the `problem` with the `solver`.
"""
solve(::AbstractProblem, ::AbstractSolver) = @error "not implemented"

"""
    solve(problem, solver)

Solve the simulation `problem` with the simulation `solver`.

### Notes

Default implementation calls `solvesingle` in parallel.
"""
function solve(problem::SimulationProblem, solver::AbstractSimulationSolver)
  # sanity checks
  @assert variables(solver) ⊆ keys(variables(problem)) "invalid variables in solver"

  # optional preprocessing
  preproc = preprocess(problem, solver)

  # simulation loop
  results = []
  for covars in covariables(problem, solver)
    # simulate covariables
    reals = pmap(1:nreals(problem)) do _
      solvesingle(problem, covars, solver, preproc)
    end

    # rearrange realizations
    vnames = covars.names
    vtypes = [variables(problem)[var] for var in vnames]
    rdict  = Dict(vnames .=> [Vector{V}[] for V in vtypes])
    for real in reals
      for var in vnames
        push!(rdict[var], real[var])
      end
    end

    push!(results, rdict)
  end

  # merge results into a single dictionary
  pdomain = domain(problem)
  preals  = reduce(merge, results)

  SimulationSolution(pdomain, preals)
end

"""
    preprocess(problem, solver)

Preprocess the simulation `problem` once before generating each realization
with simulation `solver`.

### Notes

The output of the function is defined by the solver developer.
Default implementation returns nothing.
"""
preprocess(::SimulationProblem, ::AbstractSimulationSolver) = nothing

"""
    solvesingle(problem, covariables, solver, preproc)

Solve a single realization of `covariables` in the simulation `problem`
with the simulation `solver`, optionally using preprocessed data in `preproc`.

### Notes

By implementing this function instead of `solve`, the developer is
informing the framework that realizations generated with his/her
solver are indenpendent one from another. GeoStats.jl will trigger
the algorithm in parallel (if enough processes are available).
"""
solvesingle(::SimulationProblem, ::NamedTuple, ::AbstractSimulationSolver,
            ::Any) = @error "not implemented"

"""
    covariables(var, solver)

Return covariables associated with the variable `var`
in the `solver` and their respective parameters.
"""
covariables(::Symbol, ::AbstractSolver) = @error "not implemented"

"""
    covariables(problem, solver)

Return all covariables in the `solver` based on list of
variables in the `problem`.
"""
function covariables(problem::AbstractProblem, solver::AbstractSolver)
  vars = Set(keys(variables(problem)))

  result = []
  while !isempty(vars)
    # choose a variable from the problem
    var = pop!(vars)

    # find covariables of the variable
    covars = covariables(var, solver)

    # save covariables to result
    push!(result, covars)

    # update remaining variables
    for v in setdiff(covars.names, [var])
      pop!(vars, v)
    end
  end

  result
end

"""
    variables(solver)

Return flattened list of variable names in the `solver`.
"""
variables(::AbstractSolver) = @error "not implemented"

#------------------
# IMPLEMENTATIONS
#------------------
include("solvers/sequential_simulation.jl")
include("solvers/cookie_cutter_simulation.jl")
include("solvers/pointwise_learning.jl")


text:
"""
docstring_ext

Helper function to generate the docstring for a package.
"""
function docstring_ext(T; descr::String="")
    package_name = MLJBase.package_name(T)
    package_url  = MLJBase.package_url(T)
    model_name   = MLJBase.name(T)
    # the message to return
    message      = "$descr"
    message     *= "\n→ based on [$package_name]($package_url)."
    message     *= "\n→ do `@load $model_name pkg=\"$package_name\"` to use the model."
    message     *= "\n→ do `?$model_name` for documentation."
end

"""
metadata_pkg

Helper function to write the metadata for a package.
"""
function metadata_pkg(T; name::String="unknown", uuid::String="unknown", url::String="unknown",
                         julia::Union{Missing,Bool}=missing, license::String="unknown",
                         is_wrapper::Bool=false)
    ex = quote
        MLJBase.package_name(::Type{<:$T}) = $name
        MLJBase.package_uuid(::Type{<:$T}) = $uuid
        MLJBase.package_url(::Type{<:$T}) = $url
        MLJBase.is_pure_julia(::Type{<:$T}) = $julia
        MLJBase.package_license(::Type{<:$T}) = $license
        MLJBase.is_wrapper(::Type{<:$T}) = $is_wrapper
    end
    parentmodule(T).eval(ex)
end

"""
metadata_model

Helper function to write the metadata for a single model of a package (complements
[`metadata_ext`](@ref)).
"""
function metadata_model(T; input=Unknown, target=Unknown,
                         output=Unknown, weights::Bool=false,
                         descr::String="", path::String="")
    if isempty(path)
        path = "MLJModels.$(MLJBase.package_name(T))_.$(MLJBase.name(T))"
    end
    ex = quote
        MLJBase.input_scitype(::Type{<:$T}) = $input
        MLJBase.output_scitype(::Type{<:$T}) = $output
        MLJBase.target_scitype(::Type{<:$T}) = $target
        MLJBase.supports_weights(::Type{<:$T}) = $weights
        MLJBase.docstring(::Type{<:$T}) = MLJBase.docstring_ext($T; descr=$descr)
        MLJBase.load_path(::Type{<:$T}) = $path
    end
    parentmodule(T).eval(ex)
end


text:
ENV["QSG_RENDER_LOOP"] = "basic"
using CxxWrap # for safe_cfunction
using QML
using Observables
using GR

const qmlfile = joinpath(dirname(Base.source_path()), "qml_ex.qml")

x = randn(1000000)
y = randn(1000000)
nbins = Observable(30)

w, h = (600, 450)
zoom = Nothing

# Arguments here need to be the "reference types", hence the "Ref" suffix
function paint(p::CxxPtr{QPainter}, item::CxxPtr{JuliaPaintedItem})
  global w, h
  global zoom

  ENV["GKSwstype"] = 381
  ENV["GKSconid"] = split(repr(p.cpp_object), "@")[2]

  dev = device(p[])[]
  r = effectiveDevicePixelRatio(window(item[])[])
  w, h = width(dev) / r, height(dev) / r

  plt = gcf()
  plt[:size] = (w, h)

  if zoom === Nothing
    xmin, xmax, ymin, ymax = (-5, 5, -5, 5)
  elseif zoom != 0
    xmin, xmax, ymin, ymax = panzoom(0, 0, zoom)
  else
    xmin, xmax, ymin, ymax = inqwindow()
  end

  num_bins = Int64(round(nbins[]))
  hexbin(x, y, nbins=num_bins, xlim=(xmin, xmax), ylim=(ymin, ymax), title="nbins: $num_bins")

  return
end

function mousePosition(eventx, eventy, deltay)
  global zoom
  if deltay != 0
    zoom = deltay < 0 ? 1.02 : 1/1.02
  else
    zoom = 0
  end
  if w > h
    xn = eventx / w
    yn = (h - eventy) / w
  else
    xn = eventx / h
    yn = (h - eventy) / h
  end
  x, y = ndctowc(xn, yn)
  "($(round(x,digits=4)), $(round(y,digits=4)))"
end

load(qmlfile,
  paint_cfunction = @safe_cfunction(paint, Cvoid, (CxxPtr{QPainter}, CxxPtr{JuliaPaintedItem})),
  nbins = nbins
)

@qmlfunction mousePosition

exec()


text:
module DataScienceCourse

end


text:
# Code here taken from
# http://www.visiondummy.com/2014/04/draw-error-ellipse-representing-covariance-matrix/

##########################################################
function plot_ellipse(q::AbstractMvNormal, clr="b", label=""; alpha=1.0)
##########################################################
    plot_ellipse(mean(q), cov(q), clr, label; alpha=alpha)
end

##########################################################
function plot_ellipse(μ, Σ, clr="b", label=""; alpha=1.0)
##########################################################

    # calculate eigevanlues
    E = eigvals(Σ)
    @assert(E[1]<=E[2])

    # calculate eigenvectors, columns are the eigenvectors
    V = eigvecs(Σ)

    ϕ = atan(V[2,2], V[1,2]) # was atan2

    if(ϕ < 0.0)
        ϕ = ϕ + 2*pi;
    end

    chisquare_val = sqrt(2.41) # 70% confidence interval
    theta_grid    = LinRange(0.0, 2*pi, 500)
    a = chisquare_val*sqrt(maximum(E))
    b = chisquare_val*sqrt(minimum(E))

    # the ellipse in x and y coordinates
    ellipse_x_r  = a*cos.( theta_grid )
    ellipse_y_r  = b*sin.( theta_grid )

    # Define a rotation matrix
    R = [ cos.(ϕ) sin.(ϕ); -sin.(ϕ) cos.(ϕ) ];

    # let's rotate the ellipse to angle ϕ
    r_ellipse = [ellipse_x_r ellipse_y_r]  * R

    PyPlot.plot(μ[1], μ[2], "x", markersize=12, markeredgewidth=4, color=clr, alpha=alpha)
    PyPlot.plot(r_ellipse[:,1] .+ μ[1], r_ellipse[:,2] .+ μ[2], clr, linewidth=1, label=label, alpha=alpha)

    # plot ellipse axis (eigenvectors)
    a_pnts = [vec(LinRange(0.0, a, 500))  vec(LinRange(0.0, a, 500))*0.0    ]*R
    b_pnts = [vec(LinRange(0.0, b, 500))*0.0      vec(LinRange(0.0, b, 500))]*R
    PyPlot.plot(a_pnts[:,1] .+ μ[1], a_pnts[:,2] .+ μ[2], @sprintf("--%s",clr), linewidth=2, alpha=alpha)
    PyPlot.plot(b_pnts[:,1] .+ μ[1], b_pnts[:,2] .+ μ[2], @sprintf("--%s",clr), linewidth=2, alpha=alpha)

end


text:
module Types

export  AbstractProjectile, AP, AbstractMaterial, Impact, AbstractPenetrator, ShapeFactor, FragShapes, Fragment, Zone, Canon,Target1D,Target2D, dummyproj, Unkown, Properties, Projectile1D, Projectile2D, MaterialPenetrator, ST, AL, WF, DU, TA, TI, AT, MetallicParts, FH, MG, CI, CU, PB, TU, DO, PH, GE, UN, BN, LX,CP, SP, BG, PN, OK, NonMetallicParts, MaterialUnkown, MaterialHuman, Body
"""
AbstractPenetrator is the super type which contains the different types of penetrators
"""
abstract type AbstractPenetrator end
"""
AbstractTarget is the super type which contains all the different types of targets
"""
abstract type AbstractTarget end

"""
AbstractMaterial is the super type which contains all the different types of materials
"""
abstract type AbstractMaterial end

"""
AbstractProjectile is a subtype of AbstractPenetrator and contains all types of bullet
"""
abstract type AbstractProjectile <: AbstractPenetrator end
abstract type AbstractFragment <: AbstractPenetrator end

abstract type AbstractCritical <: AbstractTarget end

struct Properties
    bhn
    density
    bulk
    youngs
    hugoniot
    dynamic_yield
    static_ultimate
    plastic_strain
    plastic_stress
    static__yield
    static_shear
    dynamic_shear
    plastic_shear
    fracture
end

struct MetallicParts{x} <:AbstractMaterial
prop
end

struct NonMetallicParts{x} <:AbstractMaterial
prop
end

#struct MetallicPerson{x} <:AbstractMaterial
#nude::Bool
#end

struct MaterialPenetrator{x} <:AbstractMaterial
prop
end

struct MaterialUnkown{x} <:AbstractMaterial

end

struct MaterialHuman{x} <:AbstractMaterial

end

struct NonCriticalParts{x} <:AbstractTarget
    mat
end

mutable struct FirePower{x}  <: AbstractCritical
mat
hit
end
mutable struct Mobility{x} <: AbstractCritical
mat
hit
end
mutable struct Crew{x} <: AbstractCritical
#bodypart
mat
nude::Bool
hit
end
mutable struct Communication{x} <: AbstractCritical
mat
hit
end
mutable struct CatastrophicKill{x} <: AbstractCritical
mat
hit
end

struct ST end

struct WF end

struct DU end

struct AL end

struct TA end

struct TI end

struct AT end

struct FH end

struct MG end

struct DO end

struct PH end

struct CI end

struct CU end

struct PB end

struct TU end

struct UN end

struct BN end

struct LX end

struct CP end

struct SP end

struct BG end

struct GE end

struct PN end

struct OK end

mutable struct Parallelepiped <: AbstractFragment
    length
    width
    thickness
    end

mutable struct Cylinder <: AbstractFragment
    length
    diameter
    end

mutable struct Sphere <: AbstractFragment
    diameter
    end

mutable struct Round_nose <: AbstractFragment
    diameter
    length
    end

mutable struct Tapered <: AbstractFragment
    nose_diameter
    base_diameter
    nose_length
    base_length
    weight #this is used to compute the volume based on material data
end

mutable struct AP <: AbstractProjectile
    weight #this is used to compute the volume based on material data
    bhn
    nose_diameter
    base_diameter
    nose_length
    base_length
end

struct dummyproj <: AbstractPenetrator end

mutable struct Projectile1D{T<: Union{AbstractPenetrator,Nothing}, R<:Union{AbstractMaterial,Nothing}} <: AbstractPenetrator
    calibre::Float64
    type::T
    mat::R
    mass::Float64
    velocity_m::Float64 #magnitude
    e_v::Union{Array{Float64,1},Nothing} #velocity_vec #direction cosines of the unit vector directed along the initial shotline
    position::Float64
    spin::Union{Float64,Nothing}
    yaw::Union{Float64,Nothing}
    b_f::Union{Array{Float64,1},Nothing} #body_vec#b_1f, b_2f, b_3f
end
mutable struct Projectile2D{T<: Union{AbstractPenetrator,Nothing}, R<:Union{AbstractMaterial,Nothing}} <: AbstractPenetrator
    calibre::Float64
    type::T
    mat::R
    mass::Float64
    velocity::Array{Float64,1} #magnitude
    e_v::Union{Array{Float64,1},Nothing} #velocity_vec #direction cosines of the unit vector directed along the initial shotline
    position::Array{Float64,1}
    spin::Union{Float64,Nothing}
    yaw::Union{Float64,Nothing}
    b_f::Union{Array{Float64,1},Nothing} #body_vec#b_1f, b_2f, b_3f
    inertia::Union{Array{Float64,1},Nothing}
    tof::Union{Float64,Nothing}
end

mutable struct Target1D <:AbstractTarget
    size::Float64
    position::Float64
end

mutable struct Target2D <:AbstractTarget
    size::Float64
    position::Array{Float64,1}
    velocity::Array{Float64,1}
    #mass::Float64
    #density::Float64
end

mutable struct Canon
    tc::Float64
    lw::Float64
    X2w::Float64
end
#mutable struct component <: AbstractTarget
#    mat
#    critical::Bool
#end

#critical parts

 struct Ammo #<: AbstractFirePower
 end

 struct Tires #<: AbstractMobility
 end

 struct Diesel #<: AbstractCatastrophicKill
 end

 struct V8 #<: AbstractMobility
 end

 struct Radio #<: AbstractCommunication
 end

 #mutable struct target <: AbstractMaterial
#     type #::material
# end

struct Head #<: AbstractCrew
  # nude::Bool
end

struct Thorax #<: AbstractCrew
 #  nude::Bool
end

struct Abdomen #<: AbstractCrew
  # nude::Bool
end

struct Pelvis #<: AbstractCrew
    #nude::Bool
end

struct Arms #<: AbstractCrew
    #nude::Bool
end

struct Legs #<: AbstractCrew
    #nude::Bool
end

#non-critical parts
struct Hull end

struct Turret end
struct Unkown end
struct Body end

mutable struct Impact
    id::Int64
    x_position::Float64
    los::Float64
    obl::Float64
end

mutable struct Shotline
    id::Int64
    y_position::Float64
    z_position::Float64
    elements::Vector{Impact}
end

struct Zone
    angleBegin::Float64
    angleEnd::Float64
    μ_number::Float64
    σ_number::Float64
    μ_mass::Float64
    σ_mass::Float64
    μ_velocity::Float64
    σ_velocity::Float64
    density::Float64
end

mutable struct Fragment <:AbstractPenetrator
    rad::Float64
    ax::Float64
    mass::Float64
    density::Float64
    cd_sub::Float64
    Aₚ::Float64
    position::Array{Float64,1}
    velocity::Array{Float64,1}
    tof::Union{Float64,Nothing}

end

mutable struct FragShapes
    μ_cd::Float64
    σ_cd::Float64

end

mutable struct ShapeFactor
    μ::Float64
    σ::Float64
end

end