text:
using Revise
using PortfolioBedtest
import PortfolioBedtest: order, execute
using Dates
using BenchmarkTools
using Setfield

# Offensive assets
spyts = yahoo("SPY") |> x -> map(x -> TimedEvent(x.Date, (; spy = getproperty(x, Symbol("Adj Close")))), x);
efats = yahoo("EFA") |> x -> map(x -> TimedEvent(x.Date, (; efa = getproperty(x, Symbol("Adj Close")))), x);
eemts = yahoo("EEM") |> x -> map(x -> TimedEvent(x.Date, (; eem = getproperty(x, Symbol("Adj Close")))), x);
aggts = yahoo("AGG") |> x -> map(x -> TimedEvent(x.Date, (; agg = getproperty(x, Symbol("Adj Close")))), x);

# Defensive assets
lqdts = yahoo("LQD") |> x -> map(x -> TimedEvent(x.Date, (; lqd = getproperty(x, Symbol("Adj Close")))), x);
iefts = yahoo("IEF") |> x -> map(x -> TimedEvent(x.Date, (; ief = getproperty(x, Symbol("Adj Close")))), x);
shyts = yahoo("SHY") |> x -> map(x -> TimedEvent(x.Date, (; shy = getproperty(x, Symbol("Adj Close")))), x);

asts = glue(spyts, efats, eemts, aggts, lqdts, iefts, shyts)
asts = filter(x -> !any(isnan, x.event), asts)
masts = momentum_score(asts)
masts = filter(x -> !any(isnan, x.event), masts)

########################################
# Simulation part
########################################

# TODO: We can make a function, which accepts keys and translate them to indices
struct VAAStrategy{T1, T2}
    off::T1
    def::T2
end
function maxkv(v1, v2)
    i0 = v2[1]
    val0 = v1[i0]
    for i in 2:length(v2)
        i1 = v2[i]
        val1 = v1[i1]
        i0 = val1 > val0 ? i1 : i0
        val0 = val1 > val0 ? val1 : val0
    end
    return i0, val0
end
function minkv(v1, v2)
    i0 = v2[1]
    val0 = v1[i0]
    for i in 2:length(v2)
        i1 = v2[i]
        val1 = v1[i1]
        i0 = val1 < val0 ? i1 : i0
        val0 = val1 < val0 ? val1 : val0
    end
    return i0, val0
end

struct Order <: AbstractOrder
    id::Int
end

function order(strat::VAAStrategy, broker, signal, prices)
    _, v0 = minkv(signal, strat.off)
    i, v = v0 < 0 ? maxkv(signal, strat.def) : maxkv(signal, strat.off)
    return Order(i)
end

function execute(broker, order, prices)
    broker.shares[order.id] != 0 && return broker # no need to resell same share
    a = assets(broker, prices)
    nshares = floor(Int, a/prices[order.id])
    shares = NamedTuple{keys(prices)}(ntuple(i -> i == order.id ? nshares : 0, length(prices)))
    @set! broker.cash = a - nshares * prices[order.id]
    @set! broker.shares = shares
    return broker
end

broker = Broker(1_000_000.0, map(_ -> 0, asts[1].event))
strat = VAAStrategy((1, 2, 3, 4), (5, 6, 7))
signal = Signal(masts)
recorder = Recorder{Date, typeof((; shares = broker.shares, total = 10.0))}([])

# Here goes nothing...
sim(broker, strat, asts, signal, recorder)

########################################
# Plotting
########################################
using Plots

equity_curve = mutate.(x -> x.total, recorder.data)
plot(getfield.(equity_curve, :ts), getfield.(equity_curve, :event), legend = :topleft)

# Assets buy distribution

equities = map(recorder.data) do r
    mutate(r) do row
        for (k, v) in pairs(row.shares)
            v != 0 && return k
        end
        :none
    end
end

ks = (:none, keys(broker.shares)...)
findfirst(==(:shy), ks)
getfield.(equities, :event)

dts = getfield.(equities, :ts)
vals = getfield.(equities, :event)

g = scatter(dts[vals .== ks[2]], fill(1, count(==(ks[2]), vals)), legend = :outerright, markersize = 2, label = string(ks[2]), markerstrokewidth = 0.01)
for i in 3:length(ks)
    g = scatter!(dts[vals .== ks[i]], fill(i - 1, count(==(ks[i]), vals)), markersize = 2, label = string(ks[i]), markerstrokewidth = 0.01)
end
g

########################################
# Deprecated & experiments
########################################

# Move strategy to signals

# offts = mutate.(x -> (; spy = x.spy, efa = x.efa, eem = x.eem, agg = x.agg), masts)
# defts = mutate.(x -> (; lqd = x.lqd, ief = x.ief, shy = x.shy), masts)

# struct VAASignal{D, T1, T2}
#     off::Signal{D, T1}
#     def::Signal{D, T2}
# end

# function update(sig::VAASignal, ts)
#     off2 = update(sig.off, ts)
#     def2 = update(sig.def, ts)
#     VAASignal(off2, def2)
# end

# struct VAAStrategy end
# signal = VAASignal(Signal(offts), Signal(defts))


text:
const op_precedence = Dict(:+ => 1, :- => 1, :* => 2, :/ => 2, :^ => 3)

function deparse(ex::Expr, outer_precedence=0)
    if ex.head != :call
        return "$ex"
    end
    op = ex.args[1]
    args = ex.args[2:end]
    precedence = get(op_precedence, op, 0)
    if precedence == 0
        arg_list = join([deparse(arg) for arg in args], ", ")
        return "$op($arg_list)"
    end
    if length(args) == 1
        arg = deparse(args[1])
        return "$op$arg"
    end
    result = join([deparse(arg, precedence) for arg in args], " $op ")
    if precedence <= outer_precedence
        return "($result)"
    end
    return result
end

deparse(other, outer_precedence=0) = string(other)

ext: jl


text:
using Plots
using DSP
using FFTW
include("RSProject.jl")

function write_binary_file(data,filename::String,dtype::DataType=Matrix{ComplexF32})
  f = open(filename,"w");
  data = convert(dtype,data)
  write(f,data)
  close(f)
end

"""
    exportAndDisplay(u,a,iter,K,m,name)
Exports and displays PCFM algorith results.
  ...
# Arguments
- `u::Vector`: Window Function PSD.
- `a::Integer`: Log base. If zero use non log version.
- `iter::Integer`: Number of iterations.
- `K::Integer`: Oversampling factor.
- `m::Integer`: Size.
- `name::String`: Name of window function.
...
"""
function exportAndDisplay(u,a,iter,K,m,name)
    u[findall(<(-50), 10*log10.(u))] .= 10^-5
    u = abs.(u).^2
    result = RSProject.funPcfm(u,a,iter,K)
    (B,Bb,x) = RSProject.funPcfmHelper(m,K)
    s = exp.(im.*B*result)
    # write_binary_file(s,"data/"*name*".bin")
    sb = vcat(s, zeros(m-1,1))
    sbf =  fftshift(fft(sb))
    sbf = sbf ./maximum(abs.(sbf))
    display(plot(10*log10.(abs.(sbf).^2),ylim=(-50, 0)))
    display(plot!(10*log10.(u),ylim=(-50, 0)))
end
m = 512
K=3
a = 10
iter = 1000
#Guassian Window
u = gaussian((2*m-1,1),0.1; padding = 0, zerophase = false)
exportAndDisplay(u,a,iter,K,m,"Gaussian")
u[findall(<(-50), 10*log10.(u))] .= 10^-5
u = abs.(u).^2
time = @time begin
  result = RSProject.profm(sqrt.(u),iter)
end
#time2 = @time begin
  #result2 = RSProject.funPcfm(u,a,iter,K)
#end
#(B,Bb,x) = RSProject.funPcfmHelper(m,K)
#s = exp.(im.*B*result2)
#sb = vcat(s, zeros(m-1,1))
#sbf =  fftshift(fft(sb))
#sbf = sbf ./maximum(abs.(sbf))
resPlot = abs.(fftshift(fft(result)))
resPlot = resPlot./maximum(resPlot)
display(plot(abs.(resPlot).^2))
#display(plot!(abs.(sbf).^2))
#exportAndDisplay(u,a,iter,K,m,"Guassian")
#Hanning Window
#u = hanning((2*m-1,1); padding = 0, zerophase = false)
#exportAndDisplay(u,a,iter,K,m,"Hanning")
#Hamming Window
#u = hamming((2*m-1,1); padding = 0, zerophase = false)
#exportAndDisplay(u,a,iter,K,m,"Hamming")
#Tukey Window
# u = tukey((2*m-1,1), 0.5; padding = 0, zerophase = false)
# u[findall(<(-50), 10*log10.(u))] .= 10^-5
# u = abs.(u).^2
# time = @time begin
#   result = profm(sqrt.(u),iter)
# end
# time2 = @time begin
#   result2 = funPcfm(u,a,iter,K)
# end
# (B,Bb,x) = funPcfmHelper(m,K)
# s = exp.(im.*B*result2)
# sb = vcat(s, zeros(m-1,1))
# sbf =  fftshift(fft(sb))
# sbf = sbf ./maximum(abs.(sbf))
# resPlot = abs.(fftshift(fft(result)))
# resPlot = resPlot./maximum(resPlot)
# display(plot(abs.(resPlot).^2))

# display(plot!(abs.(sbf).^2))
# #exportAndDisplay(u,a,iter,K,m,"Tukey")
# #Rectangular Window

# #u = rect((2*m-1,1); padding = 0, zerophase = false)

# #exportAndDisplay(u,a,iter,K,m,"Rectangular")
# #Triangular Window
# #u = triang((2*m-1,1); padding = 0, zerophase = false)
#u = tukey((2*m-1,1), 0.5; padding = 0, zerophase = false)
#exportAndDisplay(u,a,iter,K,m,"Tukey")
#Rectangular Window
#u = rect((2*m-1,1); padding = 0, zerophase = false)
#exportAndDisplay(u,a,iter,K,m,"Rectangular")
#Triangular Window
#u = triang((2*m-1,1); padding = 0, zerophase = false)
# #exportAndDisplay(u,a,iter,K,m,"Triangular")


text:
using ITensors
using ITensors.ITensorNetworkMaps
using KrylovKit
using LinearAlgebra

include("utils.jl")

N = 4
s = siteinds("S=1/2", N; conserve_qns=true)
χ = 2
ψ = randomMPS(s, n -> isodd(n) ? "↑" : "↓"; linkdims=χ)
ℋ = OpSum()
for n in 1:(N - 1)
  ℋ .+= 0.5, "S+", n, "S-", n + 1
  ℋ .+= 0.5, "S-", n, "S+", n + 1
  ℋ .+= "Sz", n, "Sz", n + 1
end
H = MPO(ℋ, s)

T = [transfer_matrix(ψ, n) for n in 1:N]
Tᴴ = [transfer_matrix((H, ψ), n) for n in 1:N]
@show (T[1] * T[2] * T[3] * T[4] * ITensor(1))[] ≈ inner(ψ, ψ)
@show (ITensor(1) * T[1] * T[2] * T[3] * T[4])[] ≈ inner(ψ, ψ)
@show (ITensor(1) * (T[1] * T[2] * T[3] * T[4]))[] ≈ inner(ψ, ψ)
@show contract(T[1] * T[2] * T[3] * T[4])[] ≈ inner(ψ, ψ)

@show contract(prod(T))[] ≈ inner(ψ, ψ)
@show contract(prod(Tᴴ))[] ≈ inner(ψ, H, ψ)

Ttot = transfer_matrix(ψ, 1:N)
Tᴴtot = transfer_matrix((H, ψ), 1:N)
@show contract(Ttot)[] ≈ inner(ψ, ψ)

println("Time to contract transfer matrix")
@time contract(Tᴴtot)[]
println("Time to contract using `inner(ψ, H, ψ)`")
@time inner(ψ, H, ψ)
@show contract(Tᴴtot)[] ≈ inner(ψ, H, ψ)

v = cache(T)
vᴴ = cache(Tᴴ)

@show (v[2 => 3] * v[3 => 2])[] ≈ inner(ψ, ψ)
@show (v[1 => 2] * v[2 => 1])[] ≈ inner(ψ, ψ)

@show (vᴴ[2 => 3] * vᴴ[3 => 2])[] ≈ inner(ψ, H, ψ)
@show (vᴴ[1 => 2] * vᴴ[2 => 1])[] ≈ inner(ψ, H, ψ)

# One-site effective Hamiltonian
H² = hessian(H, vᴴ, 2)
ψ² = ψ[2]
@show (dag(ψ²) * H²(ψ²))[] ≈ inner(ψ, H, ψ)

# Two-site effective Hamiltonian
H²⁻³ = hessian(H, vᴴ, 2:3)
ψ²⁻³ = prod(ψ[2:3])
@show (dag(ψ²⁻³) * H²⁻³(ψ²⁻³))[] ≈ inner(ψ, H, ψ)

# Three-site effective Hamiltonian
H²⁻⁴ = hessian(H, vᴴ, 2:4)
ψ²⁻⁴ = prod(ψ[2:4])
@show (dag(ψ²⁻⁴) * H²⁻⁴(ψ²⁻⁴))[] ≈ inner(ψ, H, ψ)


text:
using EarthSphere

mylatdegree = Degree(45,0,0)
mylongdegree = Degree(25,0,0)
println(mylatdegree)
mylatlong = Earth_LatLong_Coord(mylatdegree,mylongdegree,"SW")
println(EarthSphere.sign_lat_long(mylatlong))
println(mylatlong)
my_sphere = convert(Spherical_Coord,mylatlong)
println(my_sphere)
myotherlat = Degree(45,0,0)
myotherlong = Degree(-25,0,0)
myotherlatlong = Earth_LatLong_Coord(myotherlat,myotherlong,"NW")
my_other_sphere = convert(Spherical_Coord,myotherlatlong)
mydistance = get_distance(my_sphere,my_other_sphere)
println(mydistance)


text:
function create_new_temporary_environment(cfg::Config,
                                          proj_or_mnfst::ProjectFileOrManifestFile)
    project_directory = mktempdir(; cleanup = true)
    _write_project_or_manifest(proj_or_mnfst; project_directory)

    env_without_project = copy(ENV)

    _delete_julia_env_vars!(env_without_project)

    env_without_project["JULIA_DEPOT_PATH"] = mktempdir(; cleanup = true)
    env_without_project["JULIA_PKG_SERVER"] = strip(cfg.pkg_server)

    env_with_project = copy(env_without_project)
    env_with_project["JULIA_PROJECT"] = project_directory

    temporary_environment = TemporaryEnvironment(;
        env_with_project,
        env_without_project,
    )

    return temporary_environment
end

function _write_project_or_manifest(proj_or_mnfst::ProjectFileOrManifestFile;
                                    project_directory::String)
    file_contents = proj_or_mnfst.contents
    full_file_path = joinpath(project_directory, proj_or_mnfst.name)
    open(full_file_path, "w") do io
        println(io, file_contents)
    end
    return nothing
end

function _delete_julia_env_vars!(env::AbstractDict)
    for k in keys(env)
        if _is_julia_env_var(k)
            delete!(env, k)
        end
    end
    return nothing
end

@inline function _is_julia_env_var(name::AbstractString)
    startswith(lowercase(strip(name)), "julia_") && return true
    startswith(lowercase(strip(name)), "julia-") && return true
    return false
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 51259b56-7567-11eb-16e3-63b248b060cd
begin
	using PlutoUI
	LocalResource("./assets/ice_cream_vs_whiskey.png")
end

# ╔═╡ abde3380-750f-11eb-119b-5dc493c19b15
begin
	using ML101
	using CSV
	using DataFrames
	using Dates
	using Plots
	using StatsPlots
	plotly()
end;

# ╔═╡ 62b389c8-755a-11eb-3c5e-0f650123135c
md"
# Linear Regression Analysis

#### Ice cream vs Whiskey

JingYu Ning
"

# ╔═╡ 9d2bad86-755a-11eb-1ac8-297ba703a65d
md"
## Data
"

# ╔═╡ 64ef9478-7679-11eb-0a56-7d92daf9b205
md"
**Load data**

This data is generated by `Google Trends`,
that shows the trends of ice cream and whiskey during `(2017, 2020]`.
"

# ╔═╡ da1eaa04-750f-11eb-185c-c11c2c913513
begin
	raw = CSV.read("../data/ice_cream_vs_whiskey_trends.csv", DataFrame)
	raw
end

# ╔═╡ c67bb8fc-755a-11eb-3338-a9cea1251f3c
md"
**Raw data visualization**

The plot bellow shows the trends of ice cream and whiskey.
"

# ╔═╡ 31120810-7510-11eb-36ec-db5fc38cbf89
begin
	@df raw plot(:Date, :IceCreamT, label="Ice Cream", lw=2)
	@df raw plot!(:Date, :WhiskeyT, label="Whiskey", lw=2)
	plot!(
		title="Trends",
		xlabel="Date",
		ylabel="Percentage",
		xlims=(raw.Date[1]-Month(2), raw.Date[end]+Month(2)),
		xticks=collect(raw.Date[1]:Month(6):raw.Date[end]+Week(1)),
		size=(650, 400)
	)
end

# ╔═╡ 10bd9ec6-755b-11eb-19f8-b71750df45b2
md"
**Preprocess**

Some important features is generated.

1. The average of temperature over the world is highly dependent to the date. Therefore, generate the `TempIndex` base on the proposition: $7/15\ is\ the\ hottest\ day,\ and\ 1/15\ is\ the\ coldest\ day$.

2. The trends may be dependent to the ways that temperature changes. so I labeled data by `IsTempInc` to determing the two pathes in the season cycle.

Finally, I normalize the trends data and `TempIndex` to 0-1.
"

# ╔═╡ e932d9fa-7513-11eb-1a99-5963f9568df3
begin
	df = copy(raw)
	function gen_temp_index(date::Date)
		if Month(date) == Month(1)
			temp_index = abs(date - Date(Year(date).value, 1, 15))
		else
			if dayofyear(date) <= dayofyear(Date(Year(date).value, 7, 15))
				temp_index = date - Date(Year(date).value, 1, 15)
			else
				temp_index = Date(Year(date).value+1, 1, 15) - date
			end
		end

		return temp_index.value
	end

	function gen_is_temp_inc(date::Date)
		Date(Year(date).value,1,15)<date<=Date(Year(date).value,7,15) ? true : false
	end

	df.TempIndex = map(gen_temp_index, df.Date) ./ 182.5
	df.IsTempInc = map(gen_is_temp_inc, df.Date)
	df[!, :IceCreamT] = df.IceCreamT ./ 100
	df[!, :WhiskeyT] = df.WhiskeyT ./ 100
	df
end

# ╔═╡ 3a04c9c0-767e-11eb-1760-c7beadaff557
md"
## Anaysis
"

# ╔═╡ 8f52c950-755b-11eb-0256-9767e13d82f4
md"
**Is trends Temperature dependent?**

According to the instinct, The trends of ice cream and whiskey may be hightly dependent to temperature.
"

# ╔═╡ a6adc136-755b-11eb-1a70-bd9f483b1bcf
begin
	@df df scatter(:TempIndex, :IceCreamT, label="Ice Cream")
	@df df scatter!(:TempIndex, :WhiskeyT, label="Whiskey")
	plot!(
		title="Trends",
		xlabel="Temperature Index (arb. unit)",
		ylabel="Percentage",
		xlims=(-5e-2, 1+5e-2),
		size=(650, 400)
	)
end

# ╔═╡ bb762efe-767e-11eb-20a8-d97c4e0c7e60
md"
**Linear Regression Model**

Construct two models `ice_cream_temp` and `whiskey_temp` to indicate the relation between temperature and trends of ice cream and the relation between temperature and trends of whiskey separately.
"

# ╔═╡ 7a4cc232-7564-11eb-3cfb-b727d76a46f5
begin
	ice_cream_temp = LinearRegressionModel(df, :IceCreamT, :TempIndex)
	ice_cream_temp.argv = [0.4, 0.4]
	fit!(ice_cream_temp, lr=1e-7, atol=8e-3)

	whiskey_temp = LinearRegressionModel(df, :WhiskeyT, :TempIndex)
	whiskey_temp.argv = [0.45, -0.1]
	fit!(whiskey_temp, lr=1e-7, atol=6e-3)
end

# ╔═╡ a158c064-7560-11eb-3a9c-c77cab156802
begin
	@df df scatter(:TempIndex, :IceCreamT, label="Ice Cream")
	@df df scatter!(:TempIndex, :WhiskeyT, label="Whiskey")
	plot!(
		df.TempIndex,
		x->predict(ice_cream_temp, [x]),
		label="Ice Cream Model",
		lw=5,
		color=ARGB(0, 0.5, 1)
	)
	plot!(
		df.TempIndex,
		x->predict(whiskey_temp, [x]),
		label="Whiskey Model",
		lw=5,
		color=ARGB(1, 0.5, 0)
	)
	plot!(
		title="Trends",
		xlabel="Temperature Index (arb. unit)",
		ylabel="Percentage",
		xlims=(-5e-2, 1+5e-2),
		size=(650, 400),
		legend=:topleft
	)
end

# ╔═╡ Cell order:
# ╟─62b389c8-755a-11eb-3c5e-0f650123135c
# ╟─51259b56-7567-11eb-16e3-63b248b060cd
# ╟─abde3380-750f-11eb-119b-5dc493c19b15
# ╟─9d2bad86-755a-11eb-1ac8-297ba703a65d
# ╟─64ef9478-7679-11eb-0a56-7d92daf9b205
# ╟─da1eaa04-750f-11eb-185c-c11c2c913513
# ╟─c67bb8fc-755a-11eb-3338-a9cea1251f3c
# ╟─31120810-7510-11eb-36ec-db5fc38cbf89
# ╟─10bd9ec6-755b-11eb-19f8-b71750df45b2
# ╟─e932d9fa-7513-11eb-1a99-5963f9568df3
# ╟─3a04c9c0-767e-11eb-1760-c7beadaff557
# ╟─8f52c950-755b-11eb-0256-9767e13d82f4
# ╟─a6adc136-755b-11eb-1a70-bd9f483b1bcf
# ╟─bb762efe-767e-11eb-20a8-d97c4e0c7e60
# ╠═7a4cc232-7564-11eb-3cfb-b727d76a46f5
# ╟─a158c064-7560-11eb-3a9c-c77cab156802


text:
Base.exponent(G::GroupsCore.Group) = exponent(conjugacy_classes(G))
Base.exponent(cclasses::AbstractArray) = lcm(order.(first.(cclasses)))
dixon_prime(G::GroupsCore.Group) = dixon_prime(order(G), exponent(G))

function dixon_prime(cclasses::AbstractArray)
    ordG = sum(length, cclasses)
    m = exponent(cclasses)
    return dixon_prime(ordG, m)
end

function dixon_prime(ordG::Integer, exponent::Integer)
    # `FiniteFields.GF{p}` needs `p` to be `Int` so no need to do the
    # computations of this function over `BigInt` if `ordG` is so we convert
    # to `Int`.
    p = 2 * convert(Int, isqrt(ordG))
    while true
        p = nextprime(p + 1)
        isone(p % exponent) && break # we need -1 to be in the field
    end
    return p
end

function common_esd(Ns, F::Type{FiniteFields.GF{q}}) where q
    @assert isprime(q)
    itr = iterate(Ns)
    @assert itr !== nothing
    esd = EigenSpaceDecomposition(F.(first(itr)))
    for N in Iterators.rest(Ns, last(itr))
        esd = refine(esd, F.(N, false))
        @debug N esd.eigspace_ptrs
        isdiag(esd) && return esd
    end
    return esd
end

function _multiplicities(
    chars::AbstractVector{<:Character{F}},
    cclasses = conjugacy_classes(first(chars)),
) where {F<:FiniteFields.GF}

    e = Int(exponent(cclasses))
    ie = inv(F(e))

    ω⁻¹ = inv(FiniteFields.rootofunity(F, e))
    ωs = Matrix{typeof(ω⁻¹)}(undef, e, e)
    Threads.@threads for k in 0:e-1
        ω⁻ᵏ = ω⁻¹^k
        for l in 0:e-1
            ωs[l+1, k+1] = ω⁻ᵏ^l
        end
    end

    multiplicities = zeros(Int, length(chars), length(cclasses), e)
    pmap = powermap(table(first(chars)))
    # pmap is a lazy object, computed on demand; to ensure its thread-safety
    # we precompute all its values here so in the loops below it is accessed read-only
    collect(pmap)

    for (i, χ) in enumerate(chars)
        Threads.@threads for j = 1:length(cclasses)
            for k = 0:e-1
                val = Int(ie * sum(χ[pmap[j, l]] * ωs[l+1, k+1] for l = 0:e-1))
                multiplicities[i, j, k+1] = val
            end
        end
    end

    return multiplicities
end


text:
"""
    load_from(::Type{T}, df::DataFrame) where T <: Any
Load a DataFrame `df` into a structure of type T.

!!! note

    This requires that all structure fieldnames are also DataFrame column names.
    Extra dataframe columns are acceptable, although that information will not be used.

# Arguments
- `::Type{T} where T <: Any`: Any DataType.
- `df::DataFrame`: The DataFrame storing the information to store as a DataType.

# Returns
- `x<:Any`: The DataType specified as an argument.
- `lst::Array{T} where T<:Any`: A list of elements of the DataType specified as an argument
    given a multi-row DataFrame.

# Example

```julia
df = DataFrame(from = ["State"], to = ["region"])
load_from(Rename, df)
```
"""
function load_from(::Type{T}, d::Array{Dict{Any,Any},1}) where T <: Any
    lst = if all(isarray.(T.types))
        ensurearray(load_from(T, convert_type(DataFrame, d)))
    else
        vcat(ensurearray(load_from.(T, d))...)
    end
    return size(lst)[1] == 1 ? lst[1] : lst
end

function load_from(::Type{T}, d::Dict{Any,Any}) where T <: Any
    # Fill the datatype with the values in the dictionary keys, ensuring correct t.
    (fields, types) = (string.(fieldnames(T)), T.types)
    d = _load_path(d)

    # Fill the datatype with the input.
    # if any(isarray.(types)) && !all(isarray.(types))
    if any(isarray.(types))
        # Restructure data into a list of inputs in the order and type required when
        # creating the datatype. Ensure that all array entries should, in fact, be arrays.
        inps = [_load_as_type(T, d[f], t) for (f, t) in zip(fields, types)]
        inpscorrect = isarray.(inps) .== isarray.(types)

        # If all inputs are of the correct structure, fill the data type.
        if all(inpscorrect)
            lst = [T(inps...)]
        # If some inputs are arrays when they shouldn't be, expand these into a new list of
        # dictionaries to create a list of datatypes, including all array values.
        # First, create a dictionary determining whether the entry needs to be split.
        # Then, split the dictionary into a list of arrays where necessary.
        else
            LEN = length(inps[findmax(.!inpscorrect)[2]])
            splitarray = Dict(fields[ii] => !inpscorrect[ii] for ii in 1:length(inps))
            lst = [Dict{Any,Any}(k => splitarray[k] ? d[k][ii] : d[k] for k in keys(d))
                for ii in 1:LEN]
            lst = ensurearray(load_from(T, lst))
        end
    else
        lst = ensurearray(load_from(T, convert_type(DataFrame, d)))
    end
    return size(lst)[1] == 1 ? lst[1] : lst
end

function load_from(::Type{T}, df::DataFrame) where T <: Any
    isempty(df) && @error("Cannot load datatype $T from an empty DataFrame.")

    (fields, types) = (fieldnames(T), T.types)

    # Print warning if DataFrame is missing required columns.
    missing_fields = setdiff(fields, propertynames(df))
    # if length(missing_fields) > 0
    #     @warn("DataFrame columns missing fields required to fill DataType $T" * missing_fields)
    # end

    df = df[:, ensurearray(fields)]
    # If all of the struct fields are arrays, we assume all DataFrame rows should be saved.
    if all(isarray.(T.types))
        inps = [_load_as_type(T, df[:, f], t) for (f, t) in zip(fields, types)]
        lst = [T(inps...)]
    else
        lst = [T((_load_as_type(T, row[f], t) for (f, t) in zip(fields, types))...)
            for row in eachrow(df)]
    end
    return size(lst)[1] == 1 ? lst[1] : lst
end

function load_from(::Type{Dict{T}}, df::DataFrame) where T <: Any
    return Dict(_inp_key(x) => x for x in load_from(T, df))
end

"""
    _load_path(d::Dict)
Edits directories containing a list of directories ending in a file name as one path.
"""
function _load_path(d::Dict)
    FILES = [".csv", ".xlsx", ".txt", ".map", ".set"]
    for (k, lst) in d
        if typeof(lst) .== Array{String,1}
            ii_file = [any(occursin.(FILES, x)) for x in lst]
            (ii_file[end] && .!any(ii_file[1:end - 1])) && (d[k] = joinpath(lst...))
        end
    end
    return d
end

"""
    _load_as_type(::Type{Any}, entry, type::DataType)
Converts an entry to the required DataType
"""
function _load_as_type(entry, type::DataType)
    entry = ensurearray(convert_type.(type, entry))
    (!isarray(type) && length(entry) == 1) && (entry = entry[1])
    return entry
end

_load_as_type(::Type{T}, entry, type::DataType) where T <: Any = _load_as_type(entry, type)
_load_as_type(::Type{Drop},    entry, type::Type{Any})    = _load_as_type(_load_case(entry), type)
_load_as_type(::Type{Rename},  entry, type::Type{Symbol}) = _load_as_type(_load_case(entry), type)
_load_as_type(::Type{Replace}, entry, type::Type{Any})    = _load_as_type(_load_case(entry), type)
_load_as_type(::Type{Operate}, entry, type::Type{Symbol}) = _load_as_type(_load_axis(entry), type)
_load_as_type(::Type{Parameter}, entry, type::Type{Array{Symbol,1}}) = _load_as_type(_load_index(entry), type)
_load_as_type(::Type{T}, entry::Missing, type::Type{String}) where T <: Any = _load_as_type(T, "", type)        # if we're reading in dataframe with missing values

"""
    _load_case(entry::AbstractString)
Standardizes string identifiers that indicate a case change (upper-to-lower or vice-versa)
for easier editing.
"""
function _load_case(entry::AbstractString)
    test = lowercase(entry)

    occursin("lower", test) && (entry = "lower")
    occursin("upper", test) && (entry = "uppercasefirst" == test ? "uppercasefirst" : "upper")
    occursin("titlecase", test) && (entry = "titlecase")

    "all" == test    && (entry = "all")
    "unique" == test && (entry = "unique")
    return entry
end

_load_case(entry::Any) = entry

"""
    _load_axis(entry::Any)
"""
function _load_axis(entry::AbstractString)
    entry = convert_type(String, entry)
    ("1" == entry || occursin("row", lowercase(entry))) && (entry = "row")
    ("2" == entry || occursin("col", lowercase(entry))) && (entry = "col")
    return entry
end

_load_axis(entry::Any) = _load_axis(convert_type.(String, entry))

"""
"""
function _load_index(entry::String)
    m = match(r"^[\[|\(](?<idx>.*)[\]|\)]$", entry)
    m !== nothing && (entry = m[:idx])
    return string.(split(entry, ","))
end

_load_index(entry::Any) = entry

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
__precompile__()

module Policy

# POLICY module; provides generic functions like gradloglikelihood, etc
# with Gaussian Linear Policy
# and Gaussian Radial Basis Function Policy
# types

using JLD
using ReverseDiff
using Distributions

abstract type AbstractPolicy end

# TODO may be able to do a getindex trick to address policy layers
struct GLP{T<:AbstractFloat} <: AbstractPolicy
   n::Integer
   m::Integer
   nparam::Integer

   theta::Vector{T}

   ranges::Dict{Symbol,Range}

   function GLP{T}(obsspace::Integer, actspace::Integer) where T<:AbstractFloat
      n = obsspace
      m = actspace
      d = m*n + m + m

      r=Dict{Symbol,Range}()
      r[:W] = (1:m*n)
      r[:b] = (1:m) + r[:W].stop
      r[:ls]= (1:m) + r[:b].stop
      @assert r[:ls].stop == d

      return new(n,m,d,zeros(T, d),r)
   end
end

mutable struct NN{T<:AbstractFloat} <: AbstractPolicy
   n::Integer
   m::Integer
   nparam::Integer
   nhidden::Integer
   nlayers::Integer

   theta::Vector{T}

   ranges::Dict{Symbol,Range}

   scratch1::Matrix{T}
   scratch2::Matrix{T}
   scratch3::Matrix{T}

   function NN{T}(obsspace::Integer, actspace::Integer,
                  nhidden::Integer, nlayers::Integer) where T<:AbstractFloat
      @assert nlayers == 2
      n = obsspace
      m = actspace
      nh = nhidden
      nlay = nlayers - 1 # skim off top layer
      d = (nh*n + nh) +
      (nlay*nh*nh) + nlay*nh +
      m*nh + m +
      m

      r = Dict{Symbol,Range}()
      r[:Win] = (1:n*nh)
      r[:bin] = (1:nh) + r[:Win].stop

      r[:Wh]  = (1:nh*nh*nlay) + r[:bin].stop
      r[:Wh1] = (1:nh*nh)      + r[:bin].stop
      #r[:Wh2] = (1:nh*nh)         + r[:Wh1].stop

      r[:bh]  = (1:nh*nlay) + r[:Wh].stop
      r[:bh1] = (1:nh)      + r[:Wh1].stop
      #r[:bh2] = (1:nh)         + r[:bh1].stop

      #r[:Wout]= (1:nh*m) + r[:bh2].stop
      r[:Wout]= (1:nh*m) + r[:bh1].stop
      r[:bout]= (1:m)    + r[:Wout].stop

      r[:ls]  = (1:m)    + r[:bout].stop
      @assert r[:ls].stop == d

      nn = new(n, m, d, nhidden, nlay, zeros(T, d), r,
               zeros(T, 1,1), zeros(T, 1,1), zeros(T, 1,1))

      Win, _, Wh, _, Wout, _, _ = split_theta!(nn)
      Win[:]  = sqrt(T(1.0)/(n)) * randn(T,nh*n)
      Wh[:]   = sqrt(T(1.0)/nh) * randn(T,nlay*nh*nh)
      Wout[:] = sqrt(T(1e-3)/(nh)) * randn(T,m*nh)

      return nn
   end
end

################################################################# Create Policy

function save_pol(p::AbstractPolicy, skip::Int, modelname::String, filename::String)
   jldopen(filename, "w") do file
      write(file, "policy", p)
      write(file, "skip", skip)
      write(file, "model", modelname)
   end
end

function loadpolicy(filename::String)
   if isdir(filename)
      filename *= "/policy.jld"
   end
   d = load(filename)
   return d["policy"], d["skip"], d["model"]
end

############################################################## Policy Functions

PolType = Union{Array{Float32}, Array{Float64}, ReverseDiff.TrackedArray}

function getls(p::AbstractPolicy)
   ls = reshape( view(p.theta, p.ranges[:ls]), p.m)
   return ls
end

#function loglikelihood{T<:PolType}(theta::T, p::AbstractPolicy,
#                                   features::T, actions::T)
#   N = size(features, 2) # N = T*numT
#   ll = Array{eltype(theta)}(1,N)
#   loglikelihood!(ll, theta, p, features, actions)
#   return ll
#end

function updatetheta{T<:AbstractFloat}(p::AbstractPolicy, newtheta::Vector{T})
   p.theta[:] = newtheta
end

# API:
# each policy needs the following functions:
# getmean!     : mean policy output c = Policy(input)
# getaction!   : policy output with noise stored in c
#              : c = Policy(input) + c * exp(logstd)
# gradLL!      : gradient of loglikelihood function wrt to policy params
# split_theta! : policy is stored as a vector; splits vector into usable chunks
# updatetheta  : updates vector of parameters (if different than copy)

include("./GLPolicy.jl")
include("./NNPolicy.jl")

end # module Policy


text:
using IntervalArithmetic
using Test

setformat(:standard, decorations=true, sigfigs=6)

setprecision(128)

@testset "Parse string to Interval" begin
    @test parse(Interval{Float64}, "1") == Interval(1, 1)
    @test parse(Interval{Float64}, "[1, 2]") == Interval(1, 2)
    @test parse(Interval{Float64}, "[-0x1.3p-1, 2/3]") == @interval(-0x1.3p-1, 2/3)
    @test parse(Interval{Float64}, "[1,+infinity]") == Interval(1.0, Inf)
    @test parse(Interval{Float64}, "[1.234e5,Inf]") == Interval(123400.0, Inf)
    @test parse(Interval{Float64}, "[]") == ∅
    @test parse(Interval{Float64}, "[,]") == entireinterval(Float64)
    @test parse(Interval{Float64}, "[ entire ]") == entireinterval(Float64)
    @test parse(Interval{Float64}, "3.56?1") == Interval(0x3.8ccccccccccccp+0, 0x3.91eb851eb8520p+0)
    @test parse(Interval{Float64}, "3.56?1e2") == Interval(355.0, 357.0)
    @test parse(Interval{Float64}, "3.560?2") == Interval(0x3.8ed916872b020p+0, 0x3.8fdf3b645a1ccp+0)
    @test parse(Interval{Float64}, "3.56?") == Interval(0x3.8e147ae147ae0p+0, 0x3.90a3d70a3d70cp+0)
    @test parse(Interval{Float64}, "3.560?2u") == Interval(0x3.8f5c28f5c28f4p+0, 0x3.8fdf3b645a1ccp+0)
    @test parse(Interval{Float64}, "-10?") == Interval(-10.5, -9.5)
    @test parse(Interval{Float64}, "-10?u") == Interval(-10.0, -9.5)
    @test parse(Interval{Float64}, "-10?12") == Interval(-22.0, 2.0)
    @test parse(Interval{Float64}, "0.0?d") == Interval(-0.05, 0.0)
    @test parse(Interval{Float64}, "2.5?d") == Interval(0x1.3999999999999p+1, 2.5)
    @test parse(Interval{Float64}, "0.000?5d") == Interval(-0.005, 0.0)
    @test parse(Interval{Float64}, "2.500?5d") == Interval(0x1.3f5c28f5c28f5p+1, 2.5)
    @test parse(Interval{Float64}, "2.5??d") == Interval(-Inf, 2.5)
    @test parse(Interval{Float64}, "2.500?5ue4") == Interval(0x1.86ap+14, 0x1.8768p+14)
    @test parse(Interval{Float64}, "2.500?5de-5") == Interval(0x1.a2976f1cee4d5p-16, 0x1.a36e2eb1c432dp-16)
    @test parse(Interval{Float64}, "3.1416?1") == Interval(0x3.24395810624dcp+0, 0x3.24467381d7dc0p+0)

    @test parse(Interval{BigFloat}, "1") == @biginterval(1)
    @test parse(Interval{BigFloat}, "[-0x1.3p-1, 2/3]") == @biginterval(-0x1.3p-1, 2/3)

    @test parse(Interval{Float64}, "[Empty]") == emptyinterval(Float64)
    @test parse(Interval{BigFloat}, "[Empty]") == emptyinterval(BigFloat)

    @test parse(Interval{Float64}, "3 ±  4") == (-1)..7
    @test parse(Interval{Float64}, "0.2 ± 0.1") == 0.1..0.3
    @test parse(Interval{BigFloat}, "0.2 ± 0.1") == big"0.2" ± big"0.1"
end

@testset "Parse string to DecoratedInterval" begin
    @test parse(DecoratedInterval{Float64}, "[3]") == DecoratedInterval(3)

    @test parse(DecoratedInterval{Float64}, "[3, 4]") == DecoratedInterval(3, 4)

    @test parse(DecoratedInterval{Float64}, "[3, 4]_dac") == DecoratedInterval(3, 4, dac)

    @test parse(DecoratedInterval{Float64}, "[ ]") == DecoratedInterval(∅, trv)

    @test parse(DecoratedInterval{Float64}, "[entire]") == DecoratedInterval(Interval(-Inf, Inf), dac)

    @test parse(DecoratedInterval{Float64}, "[,]") == DecoratedInterval(entireinterval(Float64), dac)
    @test parse(DecoratedInterval{Float64}, "3.56?1") == DecoratedInterval(Interval(0x3.8ccccccccccccp+0, 0x3.91eb851eb8520p+0), com)
    @test parse(DecoratedInterval{Float64}, "3.56?1e2") == DecoratedInterval(Interval(355.0, 357.0), com)
    @test parse(DecoratedInterval{Float64}, "3.560?2") == DecoratedInterval(Interval(0x3.8ed916872b020p+0, 0x3.8fdf3b645a1ccp+0), com)
    @test parse(DecoratedInterval{Float64}, "3.56?") == DecoratedInterval(Interval(0x3.8e147ae147ae0p+0, 0x3.90a3d70a3d70cp+0), com)
    @test parse(DecoratedInterval{Float64}, "3.56?e2") == DecoratedInterval(Interval(3.555e2, 3.565e2), com)
    @test parse(DecoratedInterval{Float64}, "3.560?2u") == DecoratedInterval(Interval(0x3.8f5c28f5c28f4p+0, 0x3.8fdf3b645a1ccp+0), com)
    @test parse(DecoratedInterval{Float64}, "-10?") == DecoratedInterval(Interval(-10.5, -9.5), com)
    @test parse(DecoratedInterval{Float64}, "-10?e2") == DecoratedInterval(Interval(-10.5e2, -9.5e2), com)
    @test parse(DecoratedInterval{Float64}, "-10?u") == DecoratedInterval(Interval(-10.0, -9.5), com)
    @test parse(DecoratedInterval{Float64}, "-10?12") == DecoratedInterval(Interval(-22.0, 2.0), com)
    @test parse(DecoratedInterval{Float64}, "-10??u") == DecoratedInterval(Interval(-10.0, Inf), dac)
    @test parse(DecoratedInterval{Float64}, "-10??") == DecoratedInterval(Interval(-Inf, Inf), dac)
    @test isnai(parse(DecoratedInterval{Float64}, "[nai]"))
    @test parse(DecoratedInterval{Float64}, "3.56?1_def") == DecoratedInterval(Interval(0x3.8ccccccccccccp+0, 0x3.91eb851eb8520p+0), def)
end


text:
#!/usr/bin/env julia6
using MyStats, ReadGlobal

function probdf(file::AbstractString,ntimess::AbstractString="3")
    ntimes = parse(Float64,ntimess)
    field = readfield(file)
    f = string(file,".info")
    if isfile(f)
        m,stand,_,_ = read_info(f)
    else
        m = tmean(field)
        stand = tstd(field,m)
    end
    writedlm(file*".pdf.txt",histstd(field,m,stand,ntimes,100,true))
    return 0
end

probdf(ARGS...)


text:
using Octavian
using Documenter

makedocs(;
    modules=[Octavian],
    authors="Mason Protter, Chris Elrod, Dilum Aluthge, and contributors",
    repo="https://github.com/JuliaLinearAlgebra/Octavian.jl/blob/{commit}{path}#L{line}",
    sitename="Octavian.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaLinearAlgebra.github.io/Octavian.jl",
        assets=String[],
    ),
    pages=[
        "Home"                => "index.md",
        "Getting Started"     => "getting-started.md",
        "Public API"          => "public-api.md",
        "Internals (Private)" => "internals.md",
    ],
    strict=true,
)

deploydocs(;
    repo="github.com/JuliaLinearAlgebra/Octavian.jl",
)


text:
@testset "chainrules" begin
    x = randn()
    z = x + exp(randn())
    y = z + exp(randn())
    test_frule(DistributionsAD.uniformlogpdf, x, y, z)
    test_rrule(DistributionsAD.uniformlogpdf, x, y, z)
end


text:
# TODO deprecate testshuffle
_checkErrorCCWNumerics(ccwl::Union{CommonConvWrapper{F},CommonConvWrapper{Mixture{N_,F,S,T}}}, testshuffle::Bool=false)  where {N_,F<:AbstractRelativeMinimize,S,T} = nothing

function _checkErrorCCWNumerics(ccwl::Union{CommonConvWrapper{F},CommonConvWrapper{Mixture{N_,F,S,T}}},
                                testshuffle::Bool=false)  where {N_,F<:AbstractRelativeRoots,S,T}
  #
  if ccwl.zDim < ccwl.xDim && !ccwl.partial || testshuffle || ccwl.partial
    error("<:AbstractRelativeRoots factors with less measurement dimensions than variable dimensions have been discontinued, easy conversion to <:AbstractRelativeMinimize is the better option.")
  elseif !( ccwl.zDim >= ccwl.xDim && !ccwl.partial )
    error("Unresolved numeric <:AbstractRelativeRoots solve case")
  end
  nothing
end

_perturbIfNecessary(fcttype::Union{F,<:Mixture{N_,F,S,T}},
                    len::Int=1,
                    perturbation::Real=1e-10 ) where {N_,F<:AbstractRelativeMinimize,S,T} = 0
#

_perturbIfNecessary(fcttype::Union{F,<:Mixture{N_,F,S,T}},
                    len::Int=1,
                    perturbation::Real=1e-10 ) where {N_,F<:AbstractRelativeRoots,S,T} = perturbation*randn(len)
#

# internal use only, and selected out from approxDeconv functions
_solveLambdaNumeric(fcttype::AbstractPrior,
                    objResX::Function,
                    residual::AbstractVector{<:Real},
                    u0::AbstractVector{<:Real},
                    islen1::Bool=false;
                    perturb::Real=1e-10 ) = u0
#

function _solveLambdaNumeric( fcttype::Union{F,<:Mixture{N_,F,S,T}},
                              objResX::Function,
                              residual::AbstractVector{<:Real},
                              u0::AbstractVector{<:Real},
                              islen1::Bool=false )  where {N_,F<:AbstractRelativeRoots,S,T}
  #

  #
  r = NLsolve.nlsolve( (res, x) -> res .= objResX(x), u0, inplace=true)

  #
  return r.zero
end

function _solveLambdaNumeric( fcttype::Union{F,<:Mixture{N_,F,S,T}},
                              objResX::Function,
                              residual::AbstractVector{<:Real},
                              u0::AbstractVector{<:Real},
                              islen1::Bool=false  )  where {N_,F<:AbstractRelativeMinimize,S,T}
                              # retries::Int=3 )
  #
  # wrt #467 allow residual to be standardize for Roots and Minimize and Parametric cases.
  r = if islen1
    Optim.optimize((x) -> (residual .= objResX(x); sum(residual.^2)), u0, Optim.BFGS() )
  else
    Optim.optimize((x) -> (residual .= objResX(x); sum(residual.^2)), u0)
  end

  #
  return r.minimizer
end

## ================================================================================================
## Heavy dispatch for all AbstractFactor / Mixture cases below
## ================================================================================================

# internal function to dispatch view on either vector or matrix, rows are dims and samples are columns
_viewdim1or2(other, ind1, ind2) = other
_viewdim1or2(arr::AbstractVector, ind1, ind2) = view(arr, ind2)
_viewdim1or2(arr::AbstractMatrix, ind1, ind2) = view(arr, ind1, ind2)

function _buildCalcFactorMixture( ccwl::CommonConvWrapper,
                                  _fmd_,
                                  smpid,
                                  measurement_,
                                  varParams )
  #
  CalcFactor( ccwl.usrfnc!, _fmd_, smpid,
              length(measurement_), measurement_, varParams)
end

function _buildCalcFactorMixture( ccwl::CommonConvWrapper{Mixture{N_,F,S,T}},
                                  _fmd_,
                                  smpid,
                                  measurement_,
                                  varParams ) where {N_,F <: FunctorInferenceType,S,T}
  #
  # just a passthrough similar to pre-v0.20
  CalcFactor( ccwl.usrfnc!.mechanics, _fmd_, smpid,
              length(measurement_), measurement_, varParams)
end

"""
    $SIGNATURES
Internal function to build lambda pre-objective function for finding factor residuals.

Notes
- Unless passed in as separate arguments, this assumes already valid in `cpt_`:
  - `cpt_.p`
  - `cpt_.activehypo`
  - `cpt_.factormetadata`
  - `ccwl.params`
  - `ccwl.measurement`

DevNotes
- TODO refactor relationship and common fields between (CCW, FMd, CPT, CalcFactor)
"""
function _buildCalcFactorLambdaSample(ccwl::CommonConvWrapper,
                                      smpid::Int,
                                      cpt_::ConvPerThread = ccwl.cpt[Threads.threadid()],
                                      target::AbstractVector = view(ccwl.params[ccwl.varidx], cpt_.p, smpid),
                                      measurement_ = ccwl.measurement,
                                      fmd_::FactorMetadata = cpt_.factormetadata  )
  #

  # build a view to the decision variable memory
  varParams = view(ccwl.params, cpt_.activehypo)

  # prepare fmd according to hypo selection
  # FIXME must refactor (memory waste)
  _fmd_ = FactorMetadata( view(fmd_.fullvariables, cpt_.activehypo),
                          view(fmd_.variablelist, cpt_.activehypo),
                          varParams, # view(fmd_.arrRef, cpt_.activehypo),
                          fmd_.solvefor,
                          fmd_.cachedata  )
  #
  # get the operational CalcFactor object
  cf = _buildCalcFactorMixture(ccwl, _fmd_, smpid, measurement_, varParams)
  # new dev work on CalcFactor
  # cf = CalcFactor(ccwl.usrfnc!, _fmd_, smpid,
  #                 length(measurement_), measurement_, varParams)
  #

  # reset the residual vector
  fill!(cpt_.res, 0.0) # Roots->xDim | Minimize->zDim

  # build static lambda
  unrollHypo! =  ()->cf( (_viewdim1or2.(measurement_, :, smpid))..., (view.(varParams, :, smpid))... )

  return unrollHypo!, target
end

"""
    $(SIGNATURES)

Solve free variable x by root finding residual function `fgr.usrfnc(res, x)`.  This is the
penultimate step before calling numerical operations to move actual estimates, which is
done by an internally created lambda function.

ccw.X must be set to memory ref the param[varidx] being solved, at creation of ccw

Notes
- Assumes `cpt_.p` is already set to desired X decision variable dimensions and size.
- Assumes only `ccw.particleidx` will be solved for
- small random (off-manifold) perturbation used to prevent trivial solver cases, div by 0 etc.
  - perturb is necessary for NLsolve cases, and smaller than 1e-10 will result in test failure
- Also incorporates the active hypo lookup

DevNotes
- TODO testshuffle is now obsolete, should be removed
- TODO perhaps consolidate perturbation with inflation or nullhypo
"""
function _solveCCWNumeric!( ccwl::Union{CommonConvWrapper{F},
                                        CommonConvWrapper{Mixture{N_,F,S,T}}};
                            perturb::Real=1e-10,
                            testshuffle::Bool=false  ) where {N_,F<:AbstractRelative,S,T}
  #
    # FIXME, move this check higher and out of smpid loop
  _checkErrorCCWNumerics(ccwl, testshuffle)

  #
  thrid = Threads.threadid()
  cpt_ = ccwl.cpt[thrid]
  smpid = cpt_.particleidx
  # cannot Nelder-Mead on 1dim, partial can be 1dim or more but being conservative.
  islen1 = length(cpt_.p) == 1 || ccwl.partial
  # islen1 = length(cpt_.X[:, smpid]) == 1 || ccwl.partial

  # build the pre-objective function for this sample's hypothesis selection
  unrollHypo!, target = _buildCalcFactorLambdaSample(ccwl, smpid, cpt_)

  # broadcast updates original view memory location
  ## using CalcFactor legacy path inside (::CalcFactor)
  _hypoObj = (x) -> (target.=x; unrollHypo!())

  # TODO small off-manifold perturbation is a numerical workaround only, make on-manifold requires RoME.jl #244
  # use all element dimensions : ==> 1:ccwl.xDim
  target .+= _perturbIfNecessary(getFactorType(ccwl), length(target), perturb)

  # do the parameter search over defined decision variables using Minimization
  retval = _solveLambdaNumeric(getFactorType(ccwl), _hypoObj, cpt_.res, cpt_.X[cpt_.p, smpid], islen1 )

  # Check for NaNs
  if sum(isnan.(retval)) != 0
    @error "$(ccwl.usrfnc!), ccw.thrid_=$(thrid), got NaN, smpid = $(smpid), r=$(retval)\n"
    return nothing
  end

  # insert result back at the correct variable element location
  cpt_.X[cpt_.p,smpid] .= retval

  nothing
end
# brainstorming
# should only be calling a new arg list according to activehypo at start of particle
# Try calling an existing lambda
# sensitive to which hypo of course , see #1024
# need to shuffle content inside .cpt.fmd as well as .params accordingly
#

#


text:
module potentialregistry

using PotentialCalculation
using TOML
using SHA

export CASnumber,
       PotentialRegistry,
       addpotential!,
       defaultregistry,
       listpotentials,
       loadpotential,
       saveregistry

struct CASnumber
    num::Tuple{Int,Int,Int}
    CASnumber(i,j,k) = new( (i,j,k) )
    function CASnumber(s::AbstractString)
        if length(s) > 4 && s[1:4] == "CAS:"
            parts = parse.(Int,split(s[5:end],"-"))
        else
            parts = parse.(Int,split(s,"-"))
        end
        @assert length(parts) == 3
        CASnumber(parts...)
    end
end

Base.show(io::IO, c::CASnumber) = print(io, "CAS: ", c.num[1], "-", c.num[2], "-", c.num[3])
Base.isequal(c1::CASnumber, c2::CASnumber) = c1.num == c2.num

"""
requiredkeys()

Returns a list of keys required for PotentialRegistry entry.
"""
function requiredkeys()
    return ["file", "keywords", "method", "basis", "CAS", "authors", "description"]
end

"""
checkformat(d) -> Bool

Checks if dictionary has information needed to make PotentialRegistry entry.
"""
function checkformat(d; checkfile=true)
    for k in requiredkeys()
        ! checkfile && continue
        if ! haskey(d, k)
            return false
        end
    end
    try
        for x in d["CAS"]
            isa(x, CASnumber) && continue
            CASnumber(x)
        end
    catch e
        @error "CAS numbers cannot be parsed"
        return false
    end
    return true
end

mutable struct PotentialRegistry
    "Registry TOML file"
    fname::String
    data::Dict
    "Keywords descriping potential and keys for `data`"
    keywords::Dict
    "CAS numbers of molecules that are in potentials and their key for `data`"
    cas::Dict

    function PotentialRegistry(fname::String)
        if ! isfile(fname)
            touch(fname)
        end
        d = TOML.parsefile(fname)
        for (k,v) in d
            if ! checkformat(v)
                error("Error in PotentialRegistry entry $k")
            end
        end
        tmpkeys = Dict()
        for (k,v) in d
            for x in v["keywords"]
                if haskey(tmpkeys, x)
                    push!(tmpkeys[x], k)
                else
                    push!(tmpkeys, x=>[k])
                end
            end
        end
        cas = Dict()
        for (k,v) in d
            for x in v["CAS"]
                t = CASnumber(x)
                if haskey(cas, t)
                    push!(cas[t], k)
                else
                    push!(cas, t=>[k])
                end
            end
        end
        new(fname, d, tmpkeys, cas)
    end
end

Base.length(p::PotentialRegistry) = length(p.data)
Base.keys(p::PotentialRegistry) = keys(p.data)
Base.haskey(p::PotentialRegistry, key) = haskey(p.data, key)
Base.values(p::PotentialRegistry) = values(p.data)
Base.show(io::IO, p::PotentialRegistry) = print(io, "Potential registry $(length(p)) potentials")

Base.getindex(p::PotentialRegistry,i) = p.data[i]
Base.getindex(p::PotentialRegistry,i::Integer) = p["$i"]

"""
saveregistry(p::PotentialRegistry)

Saves potentila registry to file given by `p.fname`
"""
function saveregistry(p::PotentialRegistry)
    open(p.fname,"w") do f
        TOML.print(f, p.data)
    end
    @info "Potential registry saved to $(p.fname)"
    return p
end

"""
loadpotential(p::PotentialRegistry,key)

Load potential
"""
function loadpotential(p::PotentialRegistry,key)
    fname = joinpath(splitdir(p.fname)[1], p[key]["file"])
    if haskey(p[key], "hash")
        s = open(fname,"r") do f
           bytes2hex(sha2_256(f))
        end
        if s != p[key]["hash"]
            error("Potential file ", fname ," has hash: ", s, " while ", p[key]["hash"],
                  " was expected")
        end
    end
    load_jld_data(fname)
end

"""
    defaultregistry() -> `PotentialRegistry`

Loads default registry
"""
function defaultregistry()
    PotentialRegistry( normpath(joinpath( @__DIR__, "../data", "Potentials.toml")))
end

"""
addpotential!(registry::PotentialRegistry, potential::Dict, registryentry::Dict)

Add potential to registry. Potential is saved to file in same directory that
registry file is and given name after default naming scheme.
"""
function addpotential!(registry::PotentialRegistry, potential::Dict, registryentry::Dict)
    entry = deepcopy(registryentry)
    ! checkformat(entry, checkfile=false) && error("Registry entry does not have the needed information")
    tmpkeys = deepcopy(registry.keywords)
    tmpcas = deepcopy(registry.cas)
    l = length(registry) + 1
    #Check that l is not reserved already
    while true
        if haskey(registry, "$l")
            l += 1
        else
            break
        end
    end
    if haskey(entry,"file")
        pfile = joinpath(splitdir(registry.fname)[1], entry["file"])
    else
        pfile = joinpath(splitdir(registry.fname)[1], "potential-$l.jld")
        entry["file"] = "potential-$l.jld"
    end
    isfile(pfile) && error("File $(pfile) already exists")
    for x in entry["keywords"]
        if haskey(tmpkeys, x)
            push!(tmpkeys[x], "$l")
        else
            push!(tmpkeys, x=>["$l"])
        end
    end
    for x in entry["CAS"]
        t = CASnumber(x)
        if haskey(tmpcas, t)
            push!(tmpcas[t], "$l")
        else
            push!(tmpcas, t=>["$l"])
        end
    end
    registry.cas = tmpcas
    registry.keywords = tmpkeys
    save_jld_data(pfile, potential)
    s = open(pfile, "r") do f
       bytes2hex(sha2_256(f))
    end
    entry["number of points"] = length(potential["Energy"])
    push!(entry, "hash"=>s)
    push!(registry.data, "$l"=>entry)

    #TODO Check potential for content

    saveregistry(registry)
    @info "All done!"
    return registry
end

function listpotentials(r::PotentialRegistry)
    println("\n[key]  =>  Description\n")
    for k in keys(r)
        println("[$k]  =>  ", r[k]["description"])
    end
end

end #module


text:
using NumericalMethodsforEngineers, IterativeSolvers

a = [16. 4. 8.; 4. 5. -4.; 8. -4. 22.]
b = [4., 2., 5.]

upper = cholesky(a).U
lower = upper'
lower |> display
y = lower \ b
c = upper \ y

tol = 1.0e-10
maxiters = 5
x = [(i=i, cg=cg(a, b; maxiter=i)) for i in 1:maxiters]
x|> display

println("\nSolution Vector:")
x[3].cg |> display
println()

@test round.(x[3].cg; digits=6) == round.(c; digits=6)


text:
export NavierStokes

"""
    struct NavierStokes{N} <: Application{N}
Defines the compressible Navier-Stokes equations:
∂ρ/∂t + ∇⋅(ρu) = 0
∂(ρu)/∂t + ∇⋅(ρu ⊗ u' + pI - τ) = 0
∂(ρE)/∂t + ∇⋅((ρE+p)u + q - τu) = 0
where ρ is the fluid density, u is the fluid velocity, E is the total energy,
p is the pressure, I is the N × N identity tensor, τ is the viscous stress
tensor
τ := μ (∇u + ∇u' - 2/3(∇⋅u)I),
with μ=1/Re the viscosity coefficient, and q is the heat flux
q := -μ/Pr ∇(E + p/ρ - 1/2||u||_2^2) = -γ/(Pr*Re) ∇(E+1/2||u||_2^2),
with Pr the Prandtl number (typically 0.72).
"""
struct NavierStokes{N} <: Application{N}
    # Constants TODO
    " Reynolds number "
    Re
    " Prandtl number "
    Pr
    NavierStokes{N}(Re,Pr) where {N} = new{N}(Re,Pr)
    function NavierStokes{N}(param_dict::Dict) where {N}
        # TODO - get all constants from dictionary
        Re = param_dict["Re"]
        Pr = param_dict["Pr"]
        new{N}(Re, Pr)
    end
end
is_second_order(app::NavierStokes) = true
nstates(app::NavierStokes{N}) where {N} = N+2
"""
    flux_c!(flux, app::NavierStokes{N}, u)
Computes the convective (inviscid) flux for the Navier-Stokes equations:
f_c([ρ,ρu,ρE]) = [ρu, ρ u⊗u'+pI, (ρE+p)u] ∈ ℜ^(nstates × N)
"""
function flux_c!(flux, app::NavierStokes{N}, u) where {N}
    ρ = u[1]
    ρu = @view u[2:1+N]
    ρE = u[2+N]
    γ = 7/5
    p = (γ-1)*(ρE-1/2*(ρu'*ρu)/ρ)
    flux[1,:] .= ρu
    flux[2:1+N,:] .= (ρu*ρu')./ρ
    for i = 1:N
        flux[1+i,i] += p
    end
    flux[2+N,:] .= ρu*(ρE+p)./ρ
    flux
end
"""
    flux_d!(flux, app::NavierStokes{N}, u, Du)
Computes the diffusive (viscous) flux for Navier-Stokes equations:
f_d([ρ,ρu,ρE]) = [0, τ, τu-q] ∈ ℜ^(nstates × N)
"""
function flux_d!(flux, app::NavierStokes{N}, state, Dstate) where {N}
    # TODO: Time with @view in different spots
    ρ = state[1]
    ρu = @view state[2:1+N]
    ρE = state[2+N]
    ∇ρ = @view Dstate[1,:]
    ∇ρu = @view Dstate[2:1+N,:]
    ∇ρE = @view Dstate[2+N,:]
    γ = 7/5
    u = ρu/ρ
    ∇u = (∇ρu-u*∇ρ')/ρ
    div_u = sum(∇u[LinearAlgebra.diagind(∇u)])
    ∇E = (∇ρE-(ρE/ρ)*∇ρ)/ρ

    flux[1,:] .= 0.0
    # Compute viscous stress tensor τ in-place
    τ = @view flux[2:1+N,:]
    τ .= ∇u + ∇u'
    for i = 1:N
        τ[i,i] -= 2/3*div_u
    end
    τ ./= app.Re
    # Compute heat flux vector in-place
    flux[2+N,:] .= γ/app.Pr/app.Re*(∇E-∇u'*u)
    flux[2+N,:] .+= τ*u
    flux
end
"""
    numerical_flux_c!(flux, app::NavierStokes{N}, uK, uN, normal_k) where {N}
TODO: Roe numerical flux
"""
function numerical_flux_c!(flux, app::NavierStokes{N}, uK, uN, normal_k) where {N}
    # TODO
    error("Not yet programmed")
end

max_stars_repo_name: NoseKnowsAll/DGToolkit

max_issues_repo_name: NoseKnowsAll/DGToolkit

max_forks_repo_name: NoseKnowsAll/DGToolkit


text:
#### = = = Benchmark Standards = = = = ####
# X = [range(-100., 100., length=30001);] #
# Y = similar(X)                          #
#### - - - - - - - - - - - - - - - - - ####

#### S-shaped functions ####
"""
    Generic logistic function
"""
logistic(x::Real, x0::Real, a::Real, k::Real, c::Real) = a / (1.0 + exp(k * (x0 - x))) + c
"""
Sigmoid S-shaped function.

```julia
>>> @btime sigmoid(\$Y, \$X)
158.625 μs (0 allocations: 0 bytes)
```
"""
sigmoid(x::Real) = inv(1.0 + exp(-x))
function sigmoid!(y::VecO{S}, x::VecI{T}) where {S<:Real,T<:Real}
    @simd for i in eachindex(y)
        @inbounds y[i] = sigmoid(x[i])
    end
    return nothing
end
"""
Algebraic S-shaped function.

```julia
>>> @btime algebraic(\$Y, \$X, 3.0)
30.127 μs (0 allocations: 0 bytes)
```
"""
algebraic(x::Real, a::Real) = x / (1.0 + abs(x / a))
function algebraic!(y::VecO{S}, x::VecI{T}, a::Real) where {S<:Real,T<:Real}
    @inbounds @simd for i in eachindex(y)
        y[i] = algebraic(x[i], a)
    end
    return nothing
end

#### Sampling functions ####
"""
    Duplicate a series of random numbers.

    Example #1:
    ```julia
    >>> @macroexpand @nrand(3)
    quote
        r1 = rand()
        r2 = rand()
        r3 = rand()
    end
    ```

    Example #2:
    ```julia
    >>> @macroexpand @nrand(2, 2)
    quote
        r11 = rand()
        r21 = rand()
        r12 = rand()
        r22 = rand()
    end
    ```
"""
macro nrand(N::Int, scale::Real=1.0)
    c = isone(scale) ? :(rand()) : :($scale * rand())
    a = Vector{Expr}(undef, N)
    @inbounds for i in 1:N
        a[i] = Expr(:(=), Symbol("r", i), c)
    end
    return Expr(:escape, Expr(:block, a...))
end
macro nrand(M::Int, N::Int, scale::Real=1.0)
    c = isone(scale) ? :(rand()) : :($scale * rand())
    a = Vector{Expr}(undef, M * N)
    @inbounds for j in 1:N, i in 1:M
        a[M*(j-1) + i] = Expr(:(=), Symbol("r", i, j), c)
    end
    return Expr(:escape, Expr(:block, a...))
end
"""
    Exclusive sampling from a given collection.

    params:
    -------
    * collection := the collection to sample from
    * except     := the one to be excluded
"""
function sample(collection, except)
    ret = rand(collection)
    while ret === except
        ret = rand(collection)
    end
    return ret
end
"""
    Perform a single step of Wolford algorithm
"""
function welford_step(μ::Real, s::Real, v::Real, c::Real)
    isone(c) && return v, zero(v)
    s = s * (c - 1)
    m = μ + (v - μ) / c
    s = s + (v - μ) * (v - m)
    μ = m
    return μ, s / (c - 1)
end


text:
using StochasticDiffEq#, Plots
Random.seed!(100)
prob = prob_sde_2Dlinear

## Solve and plot
println("Solve and Plot")
#Let the solver determine the initial stepsize for you!
sol =solve(prob,SRI())

TEST_PLOT && plot(sol,plot_analytic=true)
#gui()

#Make sure it does a good job
sol.t[2] > 1e-7


text:
################################################################
# collect the data for a team
################################################################
immutable Team
  name::ASCIIString
  team_id::Int
  bots::Vector{AbstractBot}
end

function Team(d::Dict)
  return Team(d["name"], Int(d["teamId"]), map(ClientAI.make_bot, d["bots"]))
end

###############################################################
#   data types that encapsulate network messages
###############################################################
abstract AbstractMessage

immutable ConnectedMsg <: AbstractMessage
  team_id::Int
  config::Config
end

function ConnectedMsg(d::Dict)
  return ConnectedMsg(Int(d["teamId"]),  Config(d["config"]))
end

immutable StartMsg <: AbstractMessage
  you::Team
  other_teams::Vector{Team}
end

function StartMsg(d::Dict)
  return StartMsg(Team(d["you"]),  map(Team, d["otherTeams"]))
end

immutable EndMsg <: AbstractMessage
  winner_team_id::Int
end

function EndMsg(d::Dict)
  winner = d["winnerTeamId"]
  if winner == nothing
    winner = -1
  else
    winner = Int(winner)
  end
  return EndMsg(winner)
end

type EventsMsg <: AbstractMessage
  round_id::Int
  you::Team
  other_teams::Vector{Team}
  events::Vector{AbstractEvent}
  config::Config
end

function EventsMsg(d::Dict)
  rid = Int(d["roundId"])
  you = Team(d["you"])
  config = Config(d["config"])
  others = map(Team, d["otherTeams"])
  events = convert(Vector{AbstractEvent}, map(ClientAI.make_event, d["events"]))
  return EventsMsg(rid, you, others, events, config)
end

max_stars_repo_name: Condimentum/spacetyckiting

max_issues_repo_name: Condimentum/spacetyckiting

max_forks_repo_name: Condimentum/spacetyckiting


text:
# function maprepend(a,b,i,j,k,l::Int64=1)
#   # 多维矢量数组的新插值方法
#   # a replace of prepend! for multi-dim vectors
#   # i,j,k is the dimension of Array a
#   # l,j,k is the dimension of Array b
#
#   if typeof(a)==Array{Array{T,1} where T,3}
#     c = Array{Vector}(i+l,j,k)
#   elseif typeof(a)==Array{Float64,3}
#     c = Array{Float64}(i+l,j,k)
#   else
#     print("Maprepend return an error message!")
#     print("The value input could not be handled!")
#     return false
#   end
#
#   for m in 1:l
#     for n in 1:j
#       for r in 1:k
#         c[m,n,r] = b[m,n,r]
#       end
#     end
#   end
#
#   for m in (l+1):(i+l)
#     for n in 1:j
#       for r in 1:k
#         c[m,n,r] = a[m-l,n,r]
#       end
#     end
#   end
#
#   return c
# end

# 坐标系转换矩阵
function systoro(m1,ψ) #固定坐标系转为旋转坐标系
  m2 = zeros(Float64,3)
  stor = [cos(ψ) sin(ψ) 0;
          -sin(ψ) cos(ψ) 0;
          0 0 1]
  # m2[1] = cos(ψ)*m1[1]+sin(ψ)*m1[2]
  # m2[2] = -sin(ψ)*m1[1]+cos(ψ)*m1[2]
  # m2[3] = m1[3]
  m2 = stor*m1
  return m2
end

function rotosys(m2,ψ) #旋转坐标系转为固定坐标系
  m1 = zeros(Float64,3)
  rtos = inv([cos(ψ) sin(ψ) 0;
              -sin(ψ) cos(ψ) 0;
              0 0 1])
  m1 = rtos*m2
  return m1
end

function rotobeta(mro,beta) # 旋转坐标系到挥舞坐标系
  mbeta = zeros(Float64,3)
  rtb = [cos(beta)  sin(beta) 0;
          0 1 1;
          -sin(beta)  0 cos(beta)]
  mbeta = rtb*mro
  return mbeta
end

function betatoro(mbeta,beta) # 挥舞坐标系到旋转坐标系
  mro = zeros(Float64,3)
  btr = inv([cos(beta)  sin(beta) 0;
              0 1 1;
              -sin(beta)  0 cos(beta)])
  mro = btr*mbeta
  return mro
end

function betatobe(mbeta,θ) # 从挥舞坐标系转换到叶素坐标系
  mbe = zeros(Float64,3)
  betobe = [1 0 0;
            0 -cos(θ) -sin(θ);
            0 sin(θ)  -cos(θ)]
  mbe = betobe*mbeta
  return mbe
end

function rotobe(mr,beta,θ) # 从旋转坐标系转化到叶素坐标系
  mbe = zeros(Float64,3)
  rtob = [cos(beta)   -cos(θ)*sin(beta) -sin(beta)sin(θ);
          0   -cos(θ) -sin(θ);
          -sin(beta)   cos(beta)*sin(θ)  -cos(beta)*cos(θ)]
  mbe = rtob*mr
  return mbe
end

function betoro(mbe,beta,θ) # 从叶素坐标系转换到旋转坐标系
  mr = zeros(Float64,3)
  btor = inv([cos(beta)   -cos(θ)*sin(beta) -sin(beta)sin(θ);
              0   -cos(θ) -sin(θ);
              -sin(beta)   cos(beta)*sin(θ) -cos(beta)*cos(θ)])
  mr = btor*mbe
  return mr
end

# 矢量夹角函数
function aoaang(vec_,x_=[1,0]) # 二维坐标系下，任意矢量与x轴夹角(-180deg~180deg)
#  x_ = [1,0]
  if vec_[2]>=x_[2] # 判断矢量与x轴夹角方向
    return -acos(dot(vec_,x_)/(norm(vec_)*norm(x_)))
  else
    return acos(dot(vec_,x_)/(norm(vec_)*norm(x_)))
  end
end

function cosvecang(a,b) #计算夹角余弦值
  if norm(a)<=1e-3||norm(b)<=1e-3
    return 0.0
  elseif norm(a+b)<1e-3||norm(a-b)<1e-3
    return cos(π)
  elseif norm(a-b)<1e-3
    return cos(0.0)
  else
    return dot(a,b)/(norm(a)*norm(b))
  end
end

max_stars_repo_name: dingzhw/RotorTrim

max_issues_repo_name: helicopterliu/RotorTrim

max_forks_repo_name: helicopterliu/RotorTrim


text:
# Types

"""
Type representing a multi-dimensional value
"""
struct MultiValue{S,T,N,L} <: Number
  array::SArray{S,T,N,L}
end

"""
Type representing a second-order tensor
"""
const TensorValue{D,T,L} = MultiValue{Tuple{D,D},T,2,L}

"""
Type representing a first-order tensor
"""
const VectorValue{D,T} = MultiValue{Tuple{D},T,1,D}

# Constructors (MultiValue)

function (::Type{MultiValue{S}})(x::Tuple) where S<:Tuple
  array = SArray{S}(x)
  MultiValue(array)
end

function (::Type{MultiValue{S}})(x::Tuple{}) where S<:Tuple
  s = """
  Unknown element type.

  Provide element type in the corresponding type parameter.
  Examples:
  MultiValue{Tuple{0,0},Int}()
  TensorValue{0,Int}()
  VectorValue{0,Int}()
  """
  error(s)
end

function (::Type{MultiValue{S,T}})(x::Tuple) where {S<:Tuple,T}
  array = SArray{S,T}(x)
  MultiValue(array)
end

function (::Type{MultiValue{S}})(x::Vararg) where S<:Tuple
  MultiValue{S}(x)
end

function (::Type{MultiValue{S,T}})(x::Vararg) where {S<:Tuple,T}
  MultiValue{S,T}(x)
end

function MultiValue(a::StaticArray{S,T}) where {S,T}
  MultiValue{S,T}(a.data)
end

# Constructors (TensorValue)

function (::Type{TensorValue{D}})(x::Tuple) where D
  S = Tuple{D,D}
  MultiValue{S}(x)
end

function (::Type{TensorValue{0}})()
  S = Tuple{0,0}
  MultiValue{S}()
end

function (::Type{TensorValue{D}})(x::Vararg) where D
  TensorValue{D}(x)
end

function (::Type{TensorValue{D,T}})(x::Tuple) where {D,T}
  S = Tuple{D,D}
  MultiValue{S,T}(x)
end

function (::Type{TensorValue{D,T}})(x::Vararg) where {D,T}
  TensorValue{D,T}(x)
end

@generated function TensorValue(arg::NTuple{DD,T}) where {T,DD}
  SQ = sqrt(DD)
  D = ceil(Int,SQ)
  @assert D == SQ
  :( TensorValue{$D,T}(arg)  )
end

function TensorValue(args::Vararg)
  TensorValue(args)
end

function TensorValue()
  S = Tuple{0,0}
  MultiValue{S}()
end

function TensorValue(a::StaticArray)
  TensorValue(a.data)
end

# Constructors (VectorValue)

function (::Type{VectorValue{D}})(x::Tuple) where D
  S = Tuple{D}
  MultiValue{S}(x)
end

function (::Type{VectorValue{D}})(x::Vararg) where D
  VectorValue{D}(x)
end

function (::Type{VectorValue{D,T}})() where {D,T}
  S = Tuple{D}
  MultiValue{S,T}()
end

function VectorValue(arg::NTuple{D,T}) where {D,T}
  VectorValue{D,T}(arg)
end

function (::Type{VectorValue{D,T}})(x::Vararg{Number,D}) where {T,D}
  VectorValue{D,T}(x)
end

function VectorValue(args::Vararg)
  VectorValue(args)
end

function VectorValue()
  S = Tuple{0}
  MultiValue{S}()
end

function VectorValue(a::StaticArray)
  VectorValue(a.data)
end

function VectorValue(a::SVector)
  MultiValue(a)
end

function VectorValue(a::MVector)
  MultiValue(a)
end

# Initializers

function zero(::Type{<:MultiValue{S,T,N,L}}) where {S,T,N,L}
  z = zero(SArray{S,T,N,L})
  MultiValue{S,T,N,L}(z)
end

function zero(::MultiValue{S,T,N,L}) where {S,T,N,L}
  zero(MultiValue{S,T,N,L})
end

function one(::Type{<:MultiValue{S,T,N,L}}) where {S,T,N,L}
  z = one(SArray{S,T,N,L})
  MultiValue{S,T,N,L}(z)
end

function one(::MultiValue{S,T,N,L}) where {S,T,N,L}
  one(MultiValue{S,T,N,L})
end

# Conversions

function convert(::Type{<:MultiValue{S,T,N,L}},a::StaticArray{S,T,N}) where {S,T,N,L}
  MultiValue(a)
end

function convert(
  ::Type{<:MultiValue{S,T,N,L}},a::AbstractArray{T,N}) where {S,T,N,L}
  b = convert(SArray{S,T,N,L},a)
  MultiValue(b)
end

function convert(::Type{<:MultiValue{S,T,N,L}},a::NTuple{L,T}) where {S,T,N,L}
  MultiValue(SArray{S,T}(a))
end

# Misc operations on the type itself

length(::Type{<: MultiValue{S,T,N,L} where {S,T,N}}  ) where L = L

function size(::Type{MultiValue{S,T,N,L}}) where {S,T,N,L}
  A = SArray{S,T,N,L}
  size(A)
end

function size(::Type{<:MultiValue{S}}) where S
  _s(Size(S))
end

@pure _s(s::Size{T}) where T = T

# Custom type printing

function show(io::IO,v::MultiValue)
  print(io,v.array.data)
end

function show(io::IO,::MIME"text/plain",v::MultiValue)
  print(io,typeof(v))
  print(io,v.array.data)
end

# Misc

mutable(::Type{MultiValue{S,T,N,L}}) where {S,T,N,L} = MArray{S,T,N,L}


text:
assignObj(scene, superObjType, obj, actPos) = nothing

function assignObj(scene::Scene, superObjType::SuperObjCollision, obj::Object3D, actPos::Int64)
    if canCollide(obj)
        push!(superObjType.superObj, obj)
    end
end

function assignObj(scene::Scene, superObjType::SuperObjMass, obj::Object3D, actPos::Int64)
    if featureHasMass(obj)
        push!(superObjType.superObj, obj)
    end
end

function assignObj(scene::Scene, superObjType::SuperObjMovable, obj::Object3D, actPos::Int64)
    if isMovable(obj)
        push!(superObjType.superObj, obj)
        append!(superObjType.superObj, reverse(obj.children) )
        obj.interactionManner.movablePos = actPos
        obj.interactionManner.actualPos  = actPos
        for child in obj.children
            child.interactionManner.movablePos = actPos
            child.interactionManner.actualPos  = actPos
        end
    end
end

function assignObj(scene::Scene, superObjType::SuperObjForce, obj::Object3D, actPos::Int64)
    if canCollide(obj) || hasJoint(obj) || hasForceElement(obj)
        push!(superObjType.superObj, obj)
    end
end

function assignObj(scene::Scene, superObjType::SuperObjVisu, obj::Object3D, actPos::Int64)
    renderer = Modia3D.renderer[1]
    if ( isVisible(obj, renderer) || !isnothing(obj.visualizationFrame)  ) && !hasJoint(obj) && !canCollide(obj) && !hasForceElement(obj) && !hasChildJoint(obj) #&& !hasCutJoint(obj) && !featureHasMass(obj)
        # if an Object3D is for visualization only it is stored in updateVisuElements
        if !(obj in scene.updateVisuElements)
            push!(scene.updateVisuElements, obj)
        end
    end
end

function assignAll(scene::Scene, superObj::SuperObjsRow, obj::Object3D, world::Object3D, actPos::Int64)
    names = fieldnames(typeof(superObj))
    for val in names
        tmp = getfield(superObj,val)
        obj.interactionManner.originPos = actPos
        obj.interactionManner.actualPos = actPos
        assignObj(scene,tmp,obj,actPos)
    end
    return nothing
end

function fillVisuElements!(scene::Scene, obj::Object3D, world::Object3D)
    if scene.options.enableVisualization || scene.exportAnimation
        if isNotCoordinateSystem(obj) && obj.visualizeFrame != Modia3D.False &&
        ((scene.options.visualizeFrames && obj.visualizeFrame == Modia3D.Inherited) || obj.visualizeFrame == Modia3D.True)
            name = Symbol(obj.path, ".", "visualizationFrame")
            if obj != world
                addObject3DVisualizationFrame!(obj, scene.autoCoordsys, name)
            else
                addObject3DVisualizationFrame!(obj, Shapes.CoordinateSystem(length=2*scene.options.defaultFrameLength), name)
            end
            append!(scene.allVisuElements, obj.visualizationFrame)
            # if an Object3D is for visualization only it is stored in updateVisuElements
            if !(obj in scene.updateVisuElements)
                push!(scene.updateVisuElements, obj)
            end
        end
        if isVisible(obj, Modia3D.renderer[1])  # todo: should not affect animation export; move to initializeVisualization?
            push!(scene.allVisuElements, obj)
        end
    end
    return nothing
end

function assignAccVelo(tree::Vector{Object3D}, obj::Object3D)
    # compute velocity of this object
    if canCollide(obj) && !hasChildJoint(obj)
        push!(tree, obj)
    end
    # compute acceleration of this object
    if hasChildJoint(obj) #|| hasCutJoint(obj)
        obj.computeAcceleration = true
        push!(tree, obj)
    end
end

assignAccVelo(tree, obj) = nothing


text:
##
using Test, StatsModels
using DataFrames
import Logging
using Unfold
include("test_utilities.jl")

Logging.global_logger(Logging.SimpleLogger(stdout, Logging.Debug))

data, evts = loadtestdata("testCase3") #
append!(data, zeros(1000))
data = reshape(data, 1, :)
data = vcat(data, data)
data = data .+ 1 * randn(size(data)) # we have to add minimal noise, else mixed models crashes.
data_missing = Array{Union{Missing,Number}}(undef, size(data))
data_missing .= data
data_missing[4500:4600] .= missing

categorical!(evts, :subject)
f = @formula 0 ~ 1 + condA + condB + (1 + condA + condB | subject)
#f  = @formula 0~1 + (1|subject)

# cut the data into epochs
# TODO This ignores subject bounds
data_e, times = Unfold.epoch(data = data, tbl = evts, τ = (-1.0, 1.9), sfreq = 10)
evts_e, data_e = Unfold.dropMissingEpochs(evts, data_e)

##

basisfunction = firbasis(τ = (-0.1, 0.3), sfreq = 10)

f = @formula 0 ~ 1 + condA + condB # 1
Xs = designmatrix(UnfoldLinearModel, f, evts_e)
ufModel_A = unfoldfit(UnfoldLinearModel, Xs, data_e)

basisfunction = firbasis(τ = (-0.1, 0.3), sfreq = 10)
Xs = designmatrix(UnfoldLinearModel, f, evts, basisfunction)
basisfunction = firbasis(τ = (-0.1, 0.5), sfreq = 10)
Xs2 = designmatrix(UnfoldLinearModel, f, evts, basisfunction)
Xs = Xs + Xs2

ufModel_B = unfoldfit(UnfoldLinearModel, Xs, data)

f = @formula 0 ~ 1 + condA + condB + (1 + condA | subject)
Xs = designmatrix(UnfoldLinearMixedModel, f, evts_e)

ufModel_C = unfoldfit(UnfoldLinearMixedModel, Xs, data_e)

basisfunction = firbasis(τ = (-0.1, 0.3), sfreq = 10)
Xs = designmatrix(UnfoldLinearMixedModel, f, evts, basisfunction)

ufModel_D = unfoldfit(UnfoldLinearMixedModel, Xs, data)

ufA = condense_long(ufModel_A, times)
ufB = condense_long(ufModel_B)
ufC = condense_long(ufModel_C, times)
ufD = condense_long(ufModel_D)

plot(ufA.colname_basis, ufA.estimate)
plot(ufC.colname_basis, ufC.estimate)
plot(ufB.colname_basis, ufB.estimate)
plot(ufD.colname_basis, ufD.estimate)

max_forks_repo_name: unfoldtoolbox/jlunfold


text:
#===============================================================================

    dispatch.jl

    dispatches API GET requests to the appropriate act() method

===============================================================================#

import JSON
import DataStructures.OrderedDict
import Memento: debug, warn, error

"""
Parse GET request, dispatch to the appropriate `act` method, and return a 2-Tuple of a
Boolean and a JSON string. Called by the Julia webserver in response to GET requests.
NB The API is documented in docs/juliaapi_new.txt.
"""
function dispatch(
    action_key      ::AbstractString,
    request_body    ::AbstractString;
    verify          ::Bool = false
)
    debug(logger, "at dispatch() with action \"$action_key\"")
    debug(logger, "request body: " * request_body)
	success = false
    result = try
        if !(action_key in keys(ACT))
            errmsg = "dispatch: $action_key action is not recognized"
            warn(logger, errmsg)
            throw(ErrorException(errmsg))
        end

        ## NB. DefaultDict and DefaultOrderedDict constructors sometimes don't work
        ## on OrderedDict (https://github.com/JuliaLang/DataStructures.jl/issues/205)
        #=
            Issue:
            By default the JSON parser converts Arrays and Dicts to type ::Any.
            It is important to annotate the type of data structures wherever
            they are known to avoid wasting time and mem`ory downstream.
            Since the data structures are specific to each action,
            this should generally be done in the generic act() methods.
        =#
		req_parsed = JSON.parse(request_body; dicttype = OrderedDict)

        ## else
        action = Val{ACT[action_key]}
        #
        production_env = (get(ENV, "JULIA_ENV", nothing) == PRODUCTION_MODE)
        @static if !production_env
            ## this code is hidden from the parser on the BeagleBone
            if verify
                verify_input = try
                    verify_request(action, req_parsed)
                catch
                    warn(logger, "data supplied with \"$action_key\"" *
                        "request is in the wrong format")
                end ## try
            end ## if verify
        end ## if !production_env
        #
        debug(logger, "dispatching to act() from dispatch()")
        response = act(action, req_parsed; out_format = pre_json_output)
        debug(logger, "response received from act() by dispatch()")
        debug(logger, repr(response))

        if isa(response, ArgumentError)
           throw(response)
        end ## if isa

        json_response = JSON.json(response)
        #
        @static if !production_env
            ## this code is hidden from the parser on the BeagleBone
            if verify
                verify_output = try
                    verify_response(action,
                        JSON.parse(json_response, dicttype = OrderedDict))
                catch
                    warn(logger, "data returned from \"$action_key\"" *
                        " request is in the wrong format")
                end ## try
            end ## if verify
        end ## if !production_env
        #
        ## return value
		success = true
        json_response
    catch err
		JSON.json(fail(logger, err; bt = true) |> out(pre_json_output))
    end ## try

    response_body = string(result)

    debug(logger, "returning from dispatch(): success=$success")

    return (success, response_body)
end ## dispatch()

#==============================================================================#

## unused function
## get keyword arguments from request
# function get_kw_from_req(key_vec ::AbstractVector, req ::Associative)
#     pair_vec = Vector{Pair}()
#     for key in key_vec
#         if haskey(req, key)
#             push!(pair_vec, Symbol(key) => req[key])
#         end ## if
#     end ## for
#     return OrderedDict(pair_vec)
# end

## obsolete function
## testing function: construct `request_body` from separate arguments
# function args2reqb(
#     action ::String,
#     exp_id ::Integer,
#     calib_info ::Union{Integer,OrderedDict};
#     stage_id ::Integer =0,
#     step_id ::Integer =0,
#     ramp_id ::Integer =0,
#     min_reliable_cyc ::Real =5,
#     baseline_method ::String ="sigmoid",
#     baseline_cyc_bounds ::AbstractVector =[],
#     guid ::String ="",
#     extra_args ::OrderedDict =OrderedDict(),
#     wdb ::String ="dflt", ## "handle", "dflt", "connect"
#     db_key ::String ="default", ## "default", "t1", "t2"
#     db_host ::String ="localhost",
#     db_usr ::String ="root",
#     db_pswd ::String ="",
#     db_name ::String ="chaipcr",
# )
#     reqb = OrderedDict{typeof(""), Any}("calibration_info" => calib_info)
#
#     if action == "amplification"
#         reqb["experiment_id"] = exp_id
#         reqb["min_reliable_cyc"] = min_reliable_cyc
#         reqb["baseline_cyc_bounds"] = baseline_cyc_bounds
#         if step_id != 0
#             reqb["step_id"] = step_id
#         elseif ramp_id != 0
#             reqb["ramp_id"] = ramp_id
#         # else
#         #     println("No step_id or ramp_id will be specified.")
#         end
#     elseif action == "meltcurve"
#         reqb["experiment_id"] = exp_id
#         reqb["stage_id"] = stage_id
#     elseif action == "analyze"
#         reqb["experiment_info"] = OrderedDict(
#             "id"   => exp_id,
#             "guid" => guid
#         )
#     else
#         error(logger, "unrecognized action")
#     end
#
#     for key in keys(extra_args)
#         reqb[key] = extra_args[key]
#     end
#
#     if wdb == "handle"
#         reqb["db_key"] = db_key
#     elseif wdb == "dflt"
#         nothing
#     elseif wdb == "connect"
#         reqb["db_host"] = db_host
#         reqb["db_usr"]  = db_usr
#         reqb["db_pswd"] = db_pswd
#         reqb["db_name"] = db_name
#     else
#         error(logger, "`wdb` must be one of the following: \"handle\", \"dflt\", \"connect\"")
#     end
#
#     return json(reqb)
# end ## args2reqb
#
#
## test: it works
# function test0()
#     println(guids)
# end

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
mutable struct PerfectBipartiteMatchingMunkresSolver <: PerfectBipartiteMatchingSolver
  weights_matrix::Matrix{Float64}

  function PerfectBipartiteMatchingMunkresSolver()
    return new(zeros(0, 0))
  end
end

supports_solve_budgeted_linear(::PerfectBipartiteMatchingMunkresSolver) = false
supports_solve_all_budgeted_linear(::PerfectBipartiteMatchingMunkresSolver) = false
has_lp_formulation(::PerfectBipartiteMatchingMunkresSolver) = false
approximation_ratio(::PerfectBipartiteMatchingMunkresSolver) = 1.0
approximation_term(::PerfectBipartiteMatchingMunkresSolver) = 0.0

function build!(solver::PerfectBipartiteMatchingMunkresSolver, n_arms::Int)
  solver.weights_matrix = zeros(n_arms, n_arms)
  nothing
end

function solve_linear(solver::PerfectBipartiteMatchingMunkresSolver, rewards::Dict{Tuple{Int, Int}, Float64})
  fill!(solver.weights_matrix, 0.0)
  for (k, v) in rewards
    solver.weights_matrix[k[1], k[2]] = v
  end
  solver.weights_matrix = maximum(solver.weights_matrix) .- solver.weights_matrix

  assignment = Munkres.munkres(solver.weights_matrix)
  return collect(enumerate(assignment))
end


text:
using BinDeps

@BinDeps.setup

rsvg = library_dependency("rsvg", aliases = ["librsvg", "librsvg-2.2", "librsvg-2-2", "librsvg-2", "librsvg-2.so.2"])
gio = library_dependency("gio", aliases = ["libgio-2.0", "libgio-2.0-0"])

@static if Sys.islinux() begin
        provides(AptGet, "librsvg2-2", rsvg)
        provides(AptGet, "libgio", gio)
    end
end

@static if Sys.isapple() begin
        using Homebrew
        provides(Homebrew.HB, "librsvg", [rsvg], os=:Darwin)
        provides(Homebrew.HB, "libgio", [gio], os=:Darwin)
    end
end

@static if Sys.iswindows() begin
        using WinRPM
        provides(WinRPM.RPM,"librsvg-2-2",rsvg,os = :Windows)
        provides(WinRPM.RPM,"glib2",gio,os = :Windows)
    end
end

@BinDeps.install Dict([
	(:rsvg => :librsvg),
    (:gio, :libgio),
])


text:
# Benchmark suite for Molly
# Run with something like:
#   using Molly, PkgBenchmark
#   results = benchmarkpkg(Molly, BenchmarkConfig(env=Dict("JULIA_NUM_THREADS" => 16)))
#   export_markdown(out_file, results)

using Molly
using BenchmarkTools
using CUDA

using Base.Threads
using DelimitedFiles

if nthreads() > 1
    @info "The parallel benchmarks will be run as Julia is running on $(nthreads()) threads"
else
    @warn "The parallel benchmarks will not be run as Julia is running on 1 thread"
end

if CUDA.functional()
    @info "The GPU benchmarks will be run as a CUDA-enabled device is available"
else
    @warn "The GPU benchmarks will not be run as a CUDA-enabled device is not available"
end

CUDA.allowscalar(false) # Check that we never do scalar indexing on the GPU

const SUITE = BenchmarkGroup(
    [],
    "interactions" => BenchmarkGroup(),
    "spatial"      => BenchmarkGroup(),
    "simulation"   => BenchmarkGroup(),
    "protein"      => BenchmarkGroup(),
)

c1 = SVector(1.0, 1.0, 1.0)u"nm"
c2 = SVector(1.4, 1.0, 1.0)u"nm"
a1 = Atom(charge=1.0u"q", σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1")
box_size = SVector(2.0, 2.0, 2.0)u"nm"
coords = [c1, c2]
s = Simulation(atoms=[a1, a1], coords=coords, box_size=box_size)
b1 = HarmonicBond(i=1, j=2, b0=0.6u"nm", kb=100_000.0u"kJ * mol^-1 * nm^-2")

SUITE["interactions"]["LennardJones force" ] = @benchmarkable force($(LennardJones()), $(c1), $(c2), $(a1), $(a1), $(box_size))
SUITE["interactions"]["LennardJones energy"] = @benchmarkable Molly.potential_energy($(LennardJones()), $(s), 1, 2)
SUITE["interactions"]["Coulomb force"      ] = @benchmarkable force($(Coulomb()), $(c1), $(c2), $(a1), $(a1), $(box_size))
SUITE["interactions"]["Coulomb energy"     ] = @benchmarkable Molly.potential_energy($(Coulomb()), $(s), 1, 2)
SUITE["interactions"]["HarmonicBond force" ] = @benchmarkable force($(b1), $(coords), $(s))
SUITE["interactions"]["HarmonicBond energy"] = @benchmarkable Molly.potential_energy($(b1), $(s))

SUITE["spatial"]["vector1D"] = @benchmarkable vector1D($(4.0u"nm"), $(6.0u"nm"), $(10.0u"nm"))
SUITE["spatial"]["vector"  ] = @benchmarkable vector($(SVector(4.0, 1.0, 1.0)u"nm"), $(SVector(6.0, 4.0, 3.0)u"nm"), $(SVector(10.0, 5.0, 3.5)u"nm"))

n_atoms = 400
mass = 10.0u"u"
box_size = SVector(6.0, 6.0, 6.0)u"nm"
temp = 1.0u"K"
starting_coords = placediatomics(n_atoms ÷ 2, box_size, 0.2u"nm", 0.2u"nm")
starting_velocities = [velocity(mass, temp) for i in 1:n_atoms]
starting_coords_f32 = [Float32.(c) for c in starting_coords]
starting_velocities_f32 = [Float32.(c) for c in starting_velocities]

function runsim(nl::Bool, parallel::Bool, gpu_diff_safe::Bool, f32::Bool, gpu::Bool)
    n_atoms = 400
    n_steps = 200
    mass = f32 ? 10.0f0u"u" : 10.0u"u"
    box_size = f32 ? SVector(6.0f0, 6.0f0, 6.0f0)u"nm" : SVector(6.0, 6.0, 6.0)u"nm"
    timestep = f32 ? 0.02f0u"ps" : 0.02u"ps"
    temp = f32 ? 1.0f0u"K" : 1.0u"K"
    simulator = VelocityVerlet()
    thermostat = NoThermostat()
    b0 = f32 ? 0.2f0u"nm" : 0.2u"nm"
    kb = f32 ? 10_000.0f0u"kJ * mol^-1 * nm^-2" : 10_000.0u"kJ * mol^-1 * nm^-2"
    bonds = [HarmonicBond(i=((i * 2) - 1), j=(i * 2), b0=b0, kb=kb) for i in 1:(n_atoms ÷ 2)]
    specific_inter_lists = (bonds,)

    neighbor_finder = NoNeighborFinder()
    cutoff = DistanceCutoff(1.0u"nm")
    general_inters = (LennardJones(nl_only=false, cutoff=cutoff),)
    if nl
        neighbor_finder = DistanceNeighborFinder(nb_matrix=trues(n_atoms, n_atoms), n_steps=10,
                                                    dist_cutoff=f32 ? 1.5f0u"nm" : 1.5u"nm")
        general_inters = (LennardJones(nl_only=true, cutoff=cutoff),)
    end

    if gpu
        coords = cu(deepcopy(f32 ? starting_coords_f32 : starting_coords))
        velocities = cu(deepcopy(f32 ? starting_velocities_f32 : starting_velocities))
        atoms = cu([Atom(charge=f32 ? 0.0f0u"q" : 0.0u"q", mass=mass, σ=f32 ? 0.2f0u"nm" : 0.2u"nm",
                            ϵ=f32 ? 0.2f0u"kJ * mol^-1" : 0.2u"kJ * mol^-1") for i in 1:n_atoms])
    else
        coords = deepcopy(f32 ? starting_coords_f32 : starting_coords)
        velocities = deepcopy(f32 ? starting_velocities_f32 : starting_velocities)
        atoms = [Atom(charge=f32 ? 0.0f0u"q" : 0.0u"q", mass=mass, σ=f32 ? 0.2f0u"nm" : 0.2u"nm",
                        ϵ=f32 ? 0.2f0u"kJ * mol^-1" : 0.2u"kJ * mol^-1") for i in 1:n_atoms]
    end

    s = Simulation(
        simulator=simulator,
        atoms=atoms,
        specific_inter_lists=specific_inter_lists,
        general_inters=general_inters,
        coords=coords,
        velocities=velocities,
        temperature=temp,
        box_size=box_size,
        neighbor_finder=neighbor_finder,
        thermostat=thermostat,
        timestep=timestep,
        n_steps=n_steps,
        gpu_diff_safe=gpu_diff_safe,
    )

    c = simulate!(s; parallel=parallel)
    return c
end

runs = [
    ("in-place"        , [false, false, false, false, false]),
    ("in-place NL"     , [true , false, false, false, false]),
    ("in-place f32"    , [false, false, false, true , false]),
    ("out-of-place"    , [false, false, true , false, false]),
    ("out-of-place f32", [false, false, true , true , false]),
]
if nthreads() > 1
    push!(runs, ("in-place parallel"   , [false, true , false, false, false]))
    push!(runs, ("in-place NL parallel", [true , true , false, false, false]))
end
if CUDA.functional()
    push!(runs, ("out-of-place gpu"    , [false, false, true , false, true ]))
    push!(runs, ("out-of-place gpu f32", [false, false, true , true , true ]))
end

for (name, args) in runs
    runsim(args...) # Run once for setup
    SUITE["simulation"][name] = @benchmarkable runsim($(args[1]), $(args[2]), $(args[3]), $(args[4]), $(args[5]))
end

data_dir = normpath(@__DIR__, "..", "data")
ff_dir = joinpath(data_dir, "force_fields")
openmm_dir = joinpath(data_dir, "openmm_6mrr")

ff = OpenMMForceField(joinpath.(ff_dir, ["ff99SBildn.xml", "tip3p_standard.xml", "his.xml"])...)

atoms, atoms_data, specific_inter_lists, general_inters, neighbor_finder, coords, box_size = setupsystem(
    joinpath(data_dir, "6mrr_equil.pdb"), ff)

n_steps = 25
timestep = 0.0005u"ps"
velocities = SVector{3}.(eachrow(readdlm(joinpath(openmm_dir, "velocities_300K.txt"))))u"nm * ps^-1"

s = Simulation(
    simulator=VelocityVerlet(),
    atoms=atoms,
    specific_inter_lists=specific_inter_lists,
    general_inters=general_inters,
    coords=coords,
    velocities=velocities,
    box_size=box_size,
    neighbor_finder=neighbor_finder,
    timestep=timestep,
)

simulate!(s, n_steps; parallel=true)
SUITE["protein"]["in-place NL parallel"] = @benchmarkable simulate!($(s), $(n_steps); parallel=true)


text:
#===
Below are functions copied from `utils.jl` in GraphRecipes
See: https://github.com/JuliaPlots/GraphRecipes.jl/blob/master/src/utils.jl
===#

# gets the alias replacement name to replace the attribute sym in plot_attributes
function _replacement_kwarg(sym, name, plot_attributes, attribute_aliases)
  replacement = name
  for alias in attribute_aliases[sym]
      if haskey(plot_attributes, alias)
          replacement = plot_attributes[alias]
      end
  end
  replacement
end

# adds the original keyword of all aliases into plot_attributes
# and sets it to the alias value
macro process_aliases(plot_attributes, attribute_aliases)
  ex = Expr(:block)
  attributes = getfield(__module__, attribute_aliases) |> keys
  ex.args = [Expr(:(=), :($(esc(plot_attributes))[$(Meta.quot(sym))]),
                  :($(esc(_replacement_kwarg))($(QuoteNode(sym)), $(esc(sym)),
                    $(esc(plot_attributes)), $(esc(attribute_aliases))))) for sym in attributes]
  ex
end
push!(_EXCLUDE_SYMBOLS, Symbol("@process_aliases"))

# called after aliases are processed via @process_aliases
# removes any alias keys from plot_attributes
function _remove_aliases!(sym, plot_attributes, attribute_aliases)
  for alias in attribute_aliases[sym]
      if haskey(plot_attributes, alias)
          delete!(plot_attributes, alias)
      end
  end
end

#===
New PowerPlots code below
===#

# Adds the given keyword args as entries into plot_attributes
function _convert_to_attributes!(plot_attributes, kwargs)
  for (var, val) in kwargs
    plot_attributes[var] = val
  end
end

function _attributes_as_variables(plot_attributes, kwargs)
  ex = Expr(:block)
  ex.args = [(haskey(kwargs, var) ? Expr(:(=), var, Meta.quot(kwargs[var]))
     : Expr(:(=), var, Meta.quot(val))) for (var, val) in plot_attributes]
  return ex
end

macro prepare_plot_attributes(kwargs)
  ex = quote
    _kwargs = $(esc(kwargs));
    plot_attributes = copy_default_attributes(); # get a copy of the default attributes
    _attributes_as_variables(plot_attributes, _kwargs) |> eval; # create variables from the default attributes
    _convert_to_attributes!(plot_attributes, _kwargs); # add kwargs into plot_attributes
    @process_aliases plot_attributes attribute_aliases; # process the aliases and add original attributes
    # remove alias names from plot_attributes
    for arg in keys(attribute_aliases)
      _remove_aliases!(arg, plot_attributes, attribute_aliases)
    end;
  end;
  push!(ex.args, Expr(:(=), esc(:plot_attributes), :plot_attributes)); # create external plot_attributes dict
  ex
end
push!(_EXCLUDE_SYMBOLS, Symbol("@prepare_plot_attributes"))

# #=== Example Usage ===#
# function powerplot(case; spring_constant=1e-3, kwargs...)
#   # Copy the line below at the start of the plot method
#   @prepare_plot_attributes(kwargs); # creates the plot_attributes dictionary

#   println(plot_attributes) # TODO remove this in real code (it's just for debug)

#   # Rest of plotting code...
# end

"Convert a color scheme `cs` into an array of string hex colors, usable by VegaLite"
function colorscheme2array(cs::ColorSchemes.ColorScheme)
    return a = ["#$(Colors.hex(c))" for c in cs]
end


text:
struct EchoReg{B} <: AbstractRegister{B} end
Base.show(io::IO, x::EchoReg) = print(io, "echo register")
EchoReg() = EchoReg{1}()

_snameof(x::IntrinsicRoutine) = string(routine_name(x))
_snameof(x::AbstractLocations) = string(x)

function execute(::typeof(Semantic.main), ::EchoReg, op::IntrinsicRoutine)
    @info "executing $op"
    return
end

function execute(::typeof(Semantic.gate), ::EchoReg, op::IntrinsicRoutine, loc::Locations)
    loc = sprint(print_locations, loc; context = :color => true)
    @info "executing $loc => $op"
    return
end

function execute(
    ::typeof(Semantic.ctrl),
    ::EchoReg,
    op::IntrinsicRoutine,
    loc::Locations,
    ctrl::CtrlLocations,
)
    loc = sprint(print_locations, loc; context = :color => true)
    ctrl = sprint(print_locations, ctrl; context = :color => true)
    @info "executing @ctrl $(ctrl) $loc => $op"
    return
end

function execute(::typeof(Semantic.measure), ::EchoReg, loc::Locations)
    loc = sprint(print_locations, loc; context = :color => true)
    @info "executing @measure $loc"
    return 0
end

function execute(::typeof(Semantic.barrier), ::EchoReg, loc::Locations)
    loc = sprint(print_locations, loc; context = :color => true)
    @info "executing @barrier $loc"
    return
end

"""
    TraceTape <: AbstractRegister{1}

Tracing all the intrinsic quantum operations.
"""
struct TraceTape <: AbstractRegister{1}
    inst::Vector{Any}
end

TraceTape() = TraceTape([])

function Base.:(==)(a::TraceTape, b::TraceTape)
    return a.inst == b.inst
end

Base.length(x::TraceTape) = length(x.inst)
Base.getindex(x::TraceTape, idx::Int) = x.inst[idx]

function Base.show(io::IO, tape::TraceTape)
    println(io, "YaoCompiler.trace:")
    nstmt = length(tape.inst)
    for i in 1:nstmt
        stmt = tape.inst[i]
        printstyled(io, nameof(stmt.args[1]); color = :light_blue)
        print(io, "\t"^2)

        for (i, each) in enumerate(stmt.args[2:end])
            print(io, each)

            if i != length(stmt.args[2:end])
                print(io, "\t")
            end
        end

        if i != nstmt
            println(io)
        end
    end
    return
end

function execute(stub::typeof(Semantic.main), r::TraceTape, op::IntrinsicRoutine)
    push!(r.inst, Expr(:call, stub, op))
    return
end

function execute(stub::typeof(Semantic.gate), r::TraceTape, op::IntrinsicRoutine, loc::Locations)
    push!(r.inst, Expr(:call, stub, op, loc))
    return
end

function execute(
    stub::typeof(Semantic.ctrl),
    r::TraceTape,
    op::IntrinsicRoutine,
    loc::Locations,
    ctrl::CtrlLocations,
)
    push!(r.inst, Expr(:call, stub, op, loc, ctrl))
    return
end

function execute(stub::typeof(Semantic.barrier), r::TraceTape, loc::Locations)
    push!(r.inst, Expr(:call, stub, loc))
    return
end

function execute(stub::typeof(Semantic.measure), r::TraceTape, locs::Locations)
    error("one should not trace programs contain @measure, since we cannot purify this hybrid program")
end

function trace_m(ex)
    ex isa Expr && ex.head === :call || error("expect a function call")
    tape = gensym(:tape)
    if ex.args[1] in (:gate, :ctrl)
        # execute(gate/ctrl, tape, spec, loc[, ctrl])
        return quote
            $tape = $TraceTape()
            $(Expr(:call, execute, GlobalRef(Semantic, ex.args[1]), tape, ex.args[2:end]...))
            $tape
        end
    else
        # execute(main, tape, spec)
        return quote
            $tape = $TraceTape()
            $(Expr(:call, execute, GlobalRef(Semantic, :main), tape, ex))
            $tape
        end
    end
end

function echo_m(ex)
    ex isa Expr && ex.head === :call || error("expect a function call")
    tape = gensym(:tape)
    if ex.args[1] in (:gate, :ctrl)
        return Expr(:call, execute, GlobalRef(Semantic, ex.args[1]), EchoReg(), ex.args[2:end]...)
    else
        return Expr(:call, execute, GlobalRef(Semantic, :main), EchoReg(), ex)
    end
end

macro trace(ex)
    esc(trace_m(ex))
end

macro echo(ex)
    esc(echo_m(ex))
end

# NOTE: we glue ArrayReg in a seperate package


text:
function get_files(dir::String)
    return collect(Base.Iterators.flatten([files for (_, _, files) in walkdir(dir)]))
end

function get_root_and_files(dir::AbstractString)
    allfiles = Matrix{Int}(undef, 0, 2)

    allfilesraw = [(root, files) for (root, _, files) in walkdir(dir)]

    for (objectdir, objects) in allfilesraw
        for object in objects
            allfiles = cat(allfiles, [objectdir object], dims = 1)
        end
    end

    return allfiles
end

function check_uniqueness(list::AbstractVector)
    processed = []
    for i in list
        i_parsed = split(i, '_')[1]
        if i_parsed ∈ processed
            println("We have a $(i_parsed) duplicate")
            continue
        end
        push!(processed, i_parsed)
    end
end

# check_uniqueness(get_files("/Users/jakeireland/Desktop/New original/"))
check_uniqueness(get_files("/Users/jakeireland/projects/GorillaSCExperiment/data/Original"))

function main()
    # selected_images = []
    # selected_images = get_files("/Users/jakeireland/projects/GorillaSCExperiment/data/chosen_object_images_1/")
    selected_images = vcat(get_files("/Users/jakeireland/projects/GorillaSCExperiment/data/chosen_object_images_1/"), get_files("/Users/jakeireland/projects/GorillaSCExperiment/data/chosen_object_images_2/"), get_files("/Users/jakeireland/projects/GorillaSCExperiment/data/chosen_object_images_3/"), get_files("/Users/jakeireland/projects/GorillaSCExperiment/data/chosen_object_images_4/"), get_files("/Users/jakeireland/projects/GorillaSCExperiment/data/chosen_object_images_5/"), get_files("/Users/jakeireland/projects/GorillaSCExperiment/data/chosen_object_images_6/"))
    allfiles = get_root_and_files("/Users/jakeireland/projects/GorillaSCExperiment/data/THINGS/")

    i = 0
    while i < 600
        randimageindex = rand(1:size(allfiles, 1))
        randimageinfo = allfiles[randimageindex, :]
        randimage = randimageinfo[2];
        # println(randimageinfo)
        if match(r"\.jpg", randimage).match != ""
            if randimage ∉ selected_images
                # println("passed")
                i += 1
                push!(selected_images, randimage)
                cp(joinpath(randimageinfo[1], randimage), joinpath("/Users/jakeireland/projects/GorillaSCExperiment/data/chosen_object_images_7", randimage))
            end
        end
    end

    return nothing
end

# main()

max_stars_repo_name: jakewilliami/GorillaSCExperiment

max_issues_repo_name: jakewilliami/GorillaSCExperiment

max_forks_repo_name: jakewilliami/GorillaSCExperiment