text:
function calculate_costs_and_losses(input_data; equipment = "ohl")
    equipment_ = join([input_data["technology"],"_",equipment])
    if equipment_ == "ac_ohl"
        costs, equipment_details = ac_ohl_costs_and_losses(input_data)
    elseif equipment_ == "ac_ugc"
        costs, equipment_details = ac_ugc_costs_and_losses(input_data)
    elseif equipment_ == "dc_ohl"
        costs, equipment_details = dc_ohl_costs_and_losses(input_data)
    elseif equipment_ == "dc_ugc"
        costs, equipment_details = dc_ugc_costs_and_losses(input_data)
    end

    return costs, equipment_details
end

function ac_ohl_costs_and_losses(input_data)
    w = 2*pi*50 # frequency
    minimum_distance = sqrt((input_data["start_node"]["x"]-input_data["end_node"]["x"])^2 + (input_data["start_node"]["y"]-input_data["start_node"]["y"])^2)
    Avecorg=[16 25 50 70 95 120 150 170 185 210 230 240 265 300 340 380 435 450 490 550 560 680]
    A = Avecorg[end] # Start with highest cross-section
    C = -0.00003705 * input_data["voltages"]["ac_ohl"]^2 + 0.03257649 * input_data["voltages"]["ac_ohl"] + 5.965161265 # estimate capacitance in nF
    Imax = -0.00170712*A.^2 + 2.64505501*A + 97.70512024 # estimate maximum current
    Ic = w * (C*1e-9) * minimum_distance * (input_data["voltages"]["ac_ohl"] * 1e3) /sqrt(3) # Calculate capacitive current
    I = input_data["power_rating"] * 1e6 / (sqrt(3) * input_data["voltages"]["ac_ohl"] * 1e3) # Calculate transmission current
    n = ceil(I/(2*sqrt(Imax^2-(Ic/2)^2))) # Calculate minimum number of circuits
    nb = ceil(I/(n*sqrt(Imax^2-(Ic/2)^2))) # Calculate minumm number of bundles
    Imax1 = sqrt((I/(n*nb))^2+(Ic/2).^2) # Calculate actual current per conductor
    K1 = -0.00170712
    K2 = 2.64505501
    K3 = 97.70512024 - Imax1;
    A_new = min(679, real((1 / (2*K1) * (-K2 + sqrt(K2^2 - 4 * K1 * K3)))))
    Avec = Avecorg
    Aindex = A_new .< Avec
    if all(Aindex==0)
        A = Avecorg[end]
    else
        A = Avec[length(Avec) - sum(Aindex) + 1]
    end

    costs = Dict{String, Any}()
    costs["investment"] = n * ((60+0.4 * input_data["voltages"]["ac_ohl"] + 0.4 * A * nb^(1/4)) * 1e3)
    costs["installation"] = n * (500*1e3) * sqrt(input_data["voltages"]["ac_ohl"] / 400)
    equipment_details = Dict{String,Any}()
    equipment_details["cross_section"] = A
    equipment_details["circuits"] = n
    equipment_details["bundles"] = nb
    return costs, equipment_details
end

function ac_ugc_costs_and_losses(input_data)
    w = 2*pi*50 # frequency
    minimum_distance = 30;
    Avecorg_offshore = [70 95 120 150 185 240 300 400 500 630 800 1000 1200 1400]
    A_offshore = Avecorg_offshore[end]
    C_offshore = (-1.21*1e-7 * A_offshore^2 + 0.000286276 * A_offshore + 0.08212371) * 1.762251064 * input_data["voltages"]["ac_ugc"]^(-0.45531498) / 0.17
    Imax_offshore = -5.7698e-4 * A_offshore^2 + 1.1793 * A_offshore + 212.33
    Ic_max_offshore = sqrt(2) * Imax_offshore / 2
    maximum_distance = Ic_max_offshore / (w * (C_offshore * 1e-6) * (input_data["voltages"]["ac_ugc"] * 1e3 )/ sqrt(3))
    Ic_offshore = w * (C_offshore * 1e-6) * minimum_distance * (input_data["voltages"]["ac_ugc"] * 1e3) / sqrt(3)
    I_offshore = input_data["power_rating"] * 1e6 /(sqrt(3) * input_data["voltages"]["ac_ugc"] * 1e3)
    n_c_offshore = ceil(I_offshore / (sqrt(Imax_offshore^2-(Ic_offshore/2)^2)))
    Imax1_offshore = sqrt((I_offshore/n_c_offshore)^2+(Ic_offshore/2).^2)
    K1=5.7698e-4
    K2=1.1793;
    K3=212.33-Imax1_offshore;
    Anew_offshore = min(999,real((1/(2*K1).*(-K2+sqrt(K2.^2-4*K1.*K3)))))
    Avec_offshore = Avecorg_offshore
    Aindex = Anew_offshore .< Avec_offshore
    if all(Aindex==0)
        A_offshore = Avecorg_offshore[end]
    else
        A_offshore = Avec_offshore[length(Avec_offshore)-sum(Aindex)+1]
    end
    Q_offshore = sqrt(3) * Ic_offshore * input_data["voltages"]["ac_ugc"] / 1e3
    K1=5.8038;
    K2=0.044525;
    K3=0.0072;

    costs = Dict{String, Any}()

    costs["offshore"] = Dict{String, Any}()
    costs["offshore"]["investment"] = 1.3e6 * n_c_offshore * (K1+K2*exp(K3*(sqrt(3) * input_data["voltages"]["ac_ugc"] * Imax1_offshore)/1e6))/8.98+ n_c_offshore * Q_offshore *7e6 /(130*minimum_distance)
    costs["offshore"]["installation"] = 7e5 * (input_data["voltages"]["ac_ugc"]/400)^(1/2)*(n_c_offshore)^0.9

    equipment_details = Dict{String, Any}()
    equipment_details["offshore"] = Dict{String, Any}()
    equipment_details["offshore"]["cross_section"] = A_offshore
    equipment_details["offshore"]["circuits"] = n_c_offshore
    equipment_details["offshore"]["maximum_distance"] = maximum_distance
    equipment_details["offshore"]["circuits"] = n_c_offshore

    ## ONSHORE CABLES
    Avecorg = [300 400 500 630 800 1000 1200 1400 1600 2000 2500]
    A = Avecorg[end]
    C = 0.2
    Imax = 26.938 * A^(0.521);
    Ic= w * (C * 1e-6)*minimum_distance * (input_data["voltages"]["ac_ugc"] * 1e3) / sqrt(3)
    I = input_data["power_rating"] * 1e6 / (sqrt(3) * input_data["voltages"]["ac_ugc"] * 1e3)
    n_c = ceil(I / (sqrt(Imax^2 - (Ic/2)^2)))
    Imax1 = sqrt((I / n_c)^2 + (Ic/2).^2)
    Anew = (Imax1 / 26.938)^(1/0.521)
    Avec = Avecorg
    Aindex = Anew .< Avec
    if all(Aindex == 0)
        A = Avecorg[end]
    else
        A = Avec[length(Avec)-sum(Aindex)+1]
    end
    Q = sqrt(3) * Ic * input_data["voltages"]["ac_ugc"] / 1e3

    costs["onshore"] = Dict{String, Any}()
    costs["onshore"]["investment"] = 1e6 * ((input_data["voltages"]["ac_ugc"] / 400)^2+ 0.25 * A/2000)*n_c + n_c * Q * 3e6 /(130*minimum_distance) + 0.25e6*(input_data["voltages"]["ac_ugc"]/400)
    costs["onshore"]["installation"] = 1e6 * (input_data["voltages"]["ac_ugc"] / 400)^(1/2)*(n_c)^0.9
    costs["onshore"]["transition"] = 3 * n_c * 1e6 * (input_data["voltages"]["ac_ugc"]/400)^(1/2)

    equipment_details["onshore"] = Dict{String, Any}()
    equipment_details["onshore"]["cross_section"] = A
    equipment_details["onshore"]["circuits"] = n_c

    return costs, equipment_details
end

function dc_ohl_costs_and_losses(input_data)
    Avecorg = [16 25 50 70 95 120 150 170 185 210 230 240 265 300 340 380 435 450 490 550 560 680]
    A = maximum(Avecorg)
    Imax=(-0.00170712*A^2 + 2.64505501 * A +97.70512024) * 1.2
    I =  input_data["power_rating"] * 1e6 / (2 * input_data["voltages"]["dc_ohl"] * 1e3)
    n = ceil(I/(2*Imax))
    nb = ceil(I/(n*Imax))
    Imax1 = I / (n * nb)
    K1 = -0.00170712 * 1.2
    K2 = 2.64505501 * 1.2
    K3 = 97.70512024 * 1.2 - Imax1
    Anew = min(679, real((1/ (2 * K1) * (-K2 + sqrt(K2^2-4 * K1 * K3)))))
    Avec = Avecorg;
    Aindex = Anew .< Avec
    if all(Aindex == 0)
        A = Avecorg[end]
    else
        A = Avec[length(Avec)-sum(Aindex) + 1]
    end

    costs = Dict{String, Any}()
    costs["investment"] = n * (60+0.4 * input_data["voltages"]["dc_ohl"] + 0.4 * A * nb^(1/4)) * 1e3 * 2/3
    costs["installation"] = n * 2/3 * (500*1e3) * sqrt(input_data["voltages"]["dc_ohl"] / 320)
    costs["converter"] = (input_data["power_rating"] /1000) * (0.8 + 0.2 * sqrt(input_data["voltages"]["dc_ohl"]/320)) * 100e6
    equipment_details = Dict{String,Any}()
    equipment_details["cross_section"] = A
    equipment_details["circuits"] = n
    equipment_details["bundles"] = nb
    return costs, equipment_details
end

function dc_ugc_costs_and_losses(input_data)
    Avecorg = [95 120 150 185 240 300 400 500 630 800 1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 3000]
    A = maximum(Avecorg)
    Imax = 25.22218157 * A^(0.57255569)
    I = input_data["power_rating"] * 1e6/(2*input_data["voltages"]["dc_ugc"] *1e3)
    n = ceil(I/Imax)
    Imax1 = I /n
    Anew = (Imax1/25.22218157)^(1/0.57255569)
    Avec = Avecorg
    Aindex = Anew .< Avec
    if all(Aindex == 0)
        A = 3000
    else
        A = Avec[length(Avec)-sum(Aindex)+1]
    end

    costs = Dict{String, Any}()
    costs["offshore"] = Dict{String, Any}()
    costs["onshore"] = Dict{String, Any}()

    costs["onshore"]["converter"] = (input_data["power_rating"] /1000) * (0.8 + 0.2 * sqrt(input_data["voltages"]["dc_ugc"] / 320)) * 100e6
    costs["onshore"]["investment"] = n * ((0.0022 * input_data["voltages"]["dc_ugc"] -0.4338) * 1e6 + (0.4906 * input_data["voltages"]["dc_ugc"]^(-0.652)) * input_data["voltages"]["dc_ugc"] * Imax1 * 1e3)/8.94448
    costs["onshore"]["installation"] = n^0.9 * 700000 * (input_data["voltages"]["dc_ugc"]/320)^(1/2)
    costs["onshore"]["transition"] = 2 * n * 1e6 * (input_data["voltages"]["dc_ugc"]/320)^(1/2)

    costs["offshore"]["converter"] = (input_data["power_rating"] /1000) * (0.8 + 0.2 * sqrt(input_data["voltages"]["dc_ugc"]/320)) * 100e6
    costs["offshore"]["investment"] = n * ((0.0022 * input_data["voltages"]["dc_ugc"] - 0.4338) * 1e6 + (0.4906 * input_data["voltages"]["dc_ugc"]^(-0.652)) * input_data["voltages"]["dc_ugc"] * Imax1*1e3)/8.94448
    costs["offshore"]["installation"] = n^0.9 * 1e6 * (input_data["voltages"]["dc_ugc"]/320)^(1/2)

    equipment_details = Dict{String, Any}()
    equipment_details["onshore"] = Dict{String, Any}()
    equipment_details["onshore"]["cross_section"] = A
    equipment_details["onshore"]["circuits"] = n

    equipment_details["offshore"] = Dict{String, Any}()
    equipment_details["offshore"]["cross_section"] = A
    equipment_details["offshore"]["circuits"] = n
    equipment_details["offshore"]["maximum_distance"] = 5000 # dummy number
    return costs, equipment_details
end


text:
using JLD2
using FileIO
using ProgressMeter
using StatsBase

include("Simulation/enums.jl")

trajectory(times::AbstractVector{T} where T) = filter(isfinite, times) |> sort!
daily(times::AbstractArray{T} where T, max_days::Integer) = fit(Histogram, times, 0:(max_days) ).weights
daily(times::AbstractVector{T} where T) = daily(times, maximum(times))

const daily_trajectory = daily ∘ trajectory

if length(ARGS) < 1
  error("give path to the grid tiles and number of trajectories(optional)")
end

output_path = ARGS[1]
num_paths = length(ARGS) == 2 ? parse(Int, ARGS[2]) : 100

output_file = jldopen(output_path*"/daily_trajectories.jld2", "w", compress=true)
try
  hospitalization_progressions, death_progressions = load(output_path*"/run_params.jld2",
    "progressions/severe_symptoms_times",
    "progressions/death_times")

  @showprogress for i in 1:num_paths
    run_path = output_path*"/run_$i.jld2"

    if !isfile(run_path)
      @warn "could not load run file $run_path, stopping"
      break
    end

    try
      infection_times, detection_times, contact_kinds = load(run_path, "infection_times", "detection_times", "contact_kinds")

      max_infection_time = filter(!isnan, infection_times) |> maximum
      max_detection_time = filter(!isnan, detection_times) |> maximum
      max_days = min(max_infection_time, max_detection_time) |> floor |> Int

      trajectory_group = JLD2.Group(output_file, string(i))

      trajectory_group["daily_infections"] = daily( infection_times |> trajectory, max_days)
      trajectory_group["daily_detections"] = daily( detection_times |> trajectory, max_days)
      trajectory_group["daily_deaths"] = daily(infection_times.+death_progressions |> trajectory, max_days)
      trajectory_group["daily_hospitalizations"] = daily(infection_times.+hospitalization_progressions |> trajectory, max_days)
      for kind in values(instances(ContactKind))
        if kind == NoContact
          continue
        end
        trajectory_group["daily_" * lowercase(string(kind))] = daily(infection_times[contact_kinds.==Int(kind)], max_days)
      end
    catch ex
      @error "problem processing file " run_path ex
    end

  end
catch ex
  @error ex
finally
  close(output_file)
end

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ ddad81d2-184a-11eb-0a24-93759b4a0d17
using Pkg, DrWatson

# ╔═╡ ddd64660-184a-11eb-341e-a1896dcf506b
begin
	@quickactivate "StatisticalRethinkingTuring"
	using Turing
	using StatisticalRethinking
	Turing.setprogress!(false);
end

# ╔═╡ f771d004-1849-11eb-308a-fdd834f7e6b1
md"## Predict.jl"

# ╔═╡ ddd6c838-184a-11eb-3285-b704caa74767
begin
	N = 200
	df = DataFrame()
	df.weight = rand(Uniform(20, 90), N)
	f(x) = mean(df.weight) + 1.6x + rand(Normal(0, 10))
	df.height = f.(df.weight)
	Text(precis(df; io=String))
end

# ╔═╡ c4135506-2c09-11eb-2cba-a93a70db7651
md"##### Chains representation."

# ╔═╡ dddd83d0-184a-11eb-3b52-13e587ea984f
@model function ppl4_3(weights, heights)
    a ~ Normal(178, 20)
    b ~ LogNormal(0, 1)
    σ ~ Uniform(0, 50)
    μ = a .+ b .* weights
    for i in eachindex(heights)
        heights[i] ~ Normal(μ[i], σ)
    end
end

# ╔═╡ ddeb0d66-184a-11eb-0107-e5532d248439
begin
	m4_3t = ppl4_3(df.weight, df.height)
	chns4_3t = sample(m4_3t, NUTS(100, 0.65), N)
	Text(sprint(show, "text/plain", chns4_3t))
end

# ╔═╡ ddecc1a6-184a-11eb-2f79-bfe4c56ee921
md"##### DataFrame representation."

# ╔═╡ b6ba87d8-184c-11eb-2a47-5f926e636eef
begin
	chns4_3t_df = DataFrame(Array(chns4_3t, :parameters), names(chns4_3t, :parameters))
	Text(precis(chns4_3t_df; io=String))
end

# ╔═╡ 106a61ac-2c09-11eb-363c-9302f28e1808
md"##### Particles representation."

# ╔═╡ ddf925a4-184a-11eb-1556-f3b5938ba663
Particles(chns4_3t[[:a, :b, :σ]])

# ╔═╡ f8b39dc2-1875-11eb-1280-bf16e96a59c5
Particles(chns4_3t_df)

# ╔═╡ de03d602-184a-11eb-093e-a1747b291cc4
md"##### Prediction as a Chains object."

# ╔═╡ de131f84-184a-11eb-05ec-152d51bde0e2
begin
	x_test = [20, 40, 50, 60, 80]
	m_test = m4_3(x_test, fill(missing, length(x_test)))
	pred4_3t = Turing.Inference.predict(m_test, chns4_3t)
	Text(sprint(show, "text/plain", pred4_3t))
end

# ╔═╡ de14ddf8-184a-11eb-0c1d-b38c9a3d0f49
md"##### DataFrame representation of predictions."

# ╔═╡ de220e38-184a-11eb-2904-2fd0143ae1fe
begin
	pred4_3t_df = DataFrame(pred4_3t)
	Text(precis(pred4_3t_df; io=String))
end

# ╔═╡ de2b0542-184a-11eb-3bdf-59c61d532a2e
md"#### Get the mean predicted values."

# ╔═╡ de32cf16-184a-11eb-038a-1d7b097ff64f
begin
	pred_summary_df = DataFrame()
	pred_summary_df.parameters = names(pred4_3t)
	pred_summary_df.mean = mean(pred4_3t)[:, 2]
	pred_summary_df.std = vec(std(Array(group(pred4_3t, :heights)); dims = 1))
	pred_summary_df.lower = hpd(pred4_3t; alpha=0.11)[:, 2]
	pred_summary_df.upper = hpd(pred4_3t; alpha=0.11)[:, 3]
	pred_summary_df
end

# ╔═╡ de3a80c4-184a-11eb-221c-5b6199a3f0e1
begin
	scatter(df.weight, df.height, lab="Observations",leg=:bottomright)
	scatter!(x_test, pred_summary_df.mean,
    markersize=5, color=:red, lab="Predictions at x_test")
	x = range(minimum(df.weight), stop=maximum(df.weight), length=N)
	g(x) = mean(chns4_3t[:a]) + mean(chns4_3t[:b]) * x
	plot!(x, g.(x), lab="Regression line", leg=:bottomright)
end

# ╔═╡ de42317c-184a-11eb-0f91-87634b6f5007
md"## End of predict.jl"

# ╔═╡ Cell order:
# ╟─f771d004-1849-11eb-308a-fdd834f7e6b1
# ╠═ddad81d2-184a-11eb-0a24-93759b4a0d17
# ╠═ddd64660-184a-11eb-341e-a1896dcf506b
# ╠═ddd6c838-184a-11eb-3285-b704caa74767
# ╟─c4135506-2c09-11eb-2cba-a93a70db7651
# ╠═dddd83d0-184a-11eb-3b52-13e587ea984f
# ╠═ddeb0d66-184a-11eb-0107-e5532d248439
# ╟─ddecc1a6-184a-11eb-2f79-bfe4c56ee921
# ╠═b6ba87d8-184c-11eb-2a47-5f926e636eef
# ╟─106a61ac-2c09-11eb-363c-9302f28e1808
# ╠═ddf925a4-184a-11eb-1556-f3b5938ba663
# ╠═f8b39dc2-1875-11eb-1280-bf16e96a59c5
# ╟─de03d602-184a-11eb-093e-a1747b291cc4
# ╠═de131f84-184a-11eb-05ec-152d51bde0e2
# ╟─de14ddf8-184a-11eb-0c1d-b38c9a3d0f49
# ╠═de220e38-184a-11eb-2904-2fd0143ae1fe
# ╟─de2b0542-184a-11eb-3bdf-59c61d532a2e
# ╠═de32cf16-184a-11eb-038a-1d7b097ff64f
# ╠═de3a80c4-184a-11eb-221c-5b6199a3f0e1
# ╟─de42317c-184a-11eb-0f91-87634b6f5007


text:
## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 1 on the tetrahedron *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TETAHEDRON_G1N1 = ((SVector(0.25, 0.25, 0.25), 0.166666666666667),)

## 0 negative weights, 0 points outside of the tetrahedron,  total sum of the
## weights is 1/6

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 2 on the tetrahedron *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TETAHEDRON_G2N4 = (
    (SVector(0.138196601125, 0.138196601125, 0.138196601125), 0.0416666666666667),
    (SVector(0.585410196625, 0.138196601125, 0.138196601125), 0.0416666666666667),
    (SVector(0.138196601125, 0.585410196625, 0.138196601125), 0.0416666666666667),
    (SVector(0.138196601125, 0.138196601125, 0.585410196625), 0.0416666666666667))

## 0 negative weights, 0 points outside of the tetrahedron,  total sum of the
## weights is 1/6

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 3 on the tetrahedron *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TETAHEDRON_G3N5 = (
    (SVector(0.250000000000, 0.250000000000, 0.250000000000), -0.133333333333333),
    (SVector(0.166666666667, 0.166666666667, 0.166666666667), +0.075000000000000),
    (SVector(0.166666666667, 0.166666666667, 0.500000000000), +0.075000000000000),
    (SVector(0.166666666667, 0.500000000000, 0.166666666667), +0.075000000000000),
    (SVector(0.500000000000, 0.166666666667, 0.166666666667), +0.075000000000000))

## 1 negative weights, 0 points outside of the tetrahedron,  total sum of the
## weights is 1/6

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 4 on the tetrahedron *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TETAHEDRON_G4N11 = (
    (SVector(0.2500000000000, 0.2500000000000, 0.2500000000000), -0.0131555555555),
    (SVector(0.0714285714286, 0.0714285714286, 0.0714285714286), +0.0076222222222),
    (SVector(0.0714285714286, 0.0714285714286, 0.7857142857140), +0.0076222222222),
    (SVector(0.0714285714286, 0.7857142857140, 0.0714285714286), +0.0076222222222),
    (SVector(0.7857142857140, 0.0714285714286, 0.0714285714286), +0.0076222222222),
    (SVector(0.3994035761670, 0.3994035761670, 0.1005964238330), +0.0248888888888),
    (SVector(0.3994035761670, 0.1005964238330, 0.3994035761670), +0.0248888888888),
    (SVector(0.1005964238330, 0.3994035761670, 0.3994035761670), +0.0248888888888),
    (SVector(0.3994035761670, 0.1005964238330, 0.1005964238330), +0.0248888888888),
    (SVector(0.1005964238330, 0.3994035761670, 0.1005964238330), +0.0248888888888),
    (SVector(0.1005964238330, 0.1005964238330, 0.3994035761670), +0.0248888888888))

## 1 negative weights, 0 points outside of the tetrahedron,  total sum of the
## weights is 1/6

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 5 on the tetrahedron *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TETAHEDRON_G5N14 = (
    (SVector(0.0927352503109, 0.0927352503109, 0.0927352503109), 0.01224884051940),
    (SVector(0.7217942490670, 0.0927352503109, 0.0927352503109), 0.01224884051940),
    (SVector(0.0927352503109, 0.7217942490670, 0.0927352503109), 0.01224884051940),
    (SVector(0.0927352503109, 0.0927352503109, 0.7217942490670), 0.01224884051940),
    (SVector(0.3108859192630, 0.3108859192630, 0.3108859192630), 0.01878132095300),
    (SVector(0.0673422422101, 0.3108859192630, 0.3108859192630), 0.01878132095300),
    (SVector(0.3108859192630, 0.0673422422101, 0.3108859192630), 0.01878132095300),
    (SVector(0.3108859192630, 0.3108859192630, 0.0673422422101), 0.01878132095300),
    (SVector(0.4544962958740, 0.4544962958740, 0.0455037041256), 0.00709100346285),
    (SVector(0.4544962958740, 0.0455037041256, 0.4544962958740), 0.00709100346285),
    (SVector(0.0455037041256, 0.4544962958740, 0.4544962958740), 0.00709100346285),
    (SVector(0.4544962958740, 0.0455037041256, 0.0455037041256), 0.00709100346285),
    (SVector(0.0455037041256, 0.4544962958740, 0.0455037041256), 0.00709100346285),
    (SVector(0.0455037041256, 0.0455037041256, 0.4544962958740), 0.00709100346285))

## 0 negative weights, 0 points outside of the tetrahedron,  total sum of the
## weights is 1

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 6 on the tetrahedron *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TETAHEDRON_G6N24 = (
    (SVector(0.2146028712590, 0.2146028712590, 0.2146028712590), 0.006653791709700),
    (SVector(0.3561913862230, 0.2146028712590, 0.2146028712590), 0.006653791709700),
    (SVector(0.2146028712590, 0.3561913862230, 0.2146028712590), 0.006653791709700),
    (SVector(0.2146028712590, 0.2146028712590, 0.3561913862230), 0.006653791709700),
    (SVector(0.0406739585346, 0.0406739585346, 0.0406739585346), 0.001679535175883),
    (SVector(0.8779781243960, 0.0406739585346, 0.0406739585346), 0.001679535175883),
    (SVector(0.0406739585346, 0.8779781243960, 0.0406739585346), 0.001679535175883),
    (SVector(0.0406739585346, 0.0406739585346, 0.8779781243960), 0.001679535175883),
    (SVector(0.3223378901420, 0.3223378901420, 0.3223378901420), 0.009226196923950),
    (SVector(0.0329863295732, 0.3223378901420, 0.3223378901420), 0.009226196923950),
    (SVector(0.3223378901420, 0.0329863295732, 0.3223378901420), 0.009226196923950),
    (SVector(0.3223378901420, 0.3223378901420, 0.0329863295732), 0.009226196923950),
    (SVector(0.0636610018750, 0.0636610018750, 0.2696723314580), 0.008035714285717),
    (SVector(0.0636610018750, 0.2696723314580, 0.0636610018750), 0.008035714285717),
    (SVector(0.0636610018750, 0.0636610018750, 0.6030056647920), 0.008035714285717),
    (SVector(0.0636610018750, 0.6030056647920, 0.0636610018750), 0.008035714285717),
    (SVector(0.0636610018750, 0.2696723314580, 0.6030056647920), 0.008035714285717),
    (SVector(0.0636610018750, 0.6030056647920, 0.2696723314580), 0.008035714285717),
    (SVector(0.2696723314580, 0.0636610018750, 0.0636610018750), 0.008035714285717),
    (SVector(0.2696723314580, 0.0636610018750, 0.6030056647920), 0.008035714285717),
    (SVector(0.2696723314580, 0.6030056647920, 0.0636610018750), 0.008035714285717),
    (SVector(0.6030056647920, 0.0636610018750, 0.2696723314580), 0.008035714285717),
    (SVector(0.6030056647920, 0.0636610018750, 0.0636610018750), 0.008035714285717),
    (SVector(0.6030056647920, 0.2696723314580, 0.0636610018750), 0.008035714285717))

## 0 negative weights, 0 points outside of the tetrahedron,  total sum of the
## weights is 1

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 7 on the tetrahedron *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TETAHEDRON_G7N31 = (
    (SVector(0.50000000000000, 0.50000000000000, 0.00000000000000), +0.000970017636685),
    (SVector(0.50000000000000, 0.00000000000000, 0.50000000000000), +0.000970017636685),
    (SVector(0.00000000000000, 0.50000000000000, 0.50000000000000), +0.000970017636685),
    (SVector(0.00000000000000, 0.00000000000000, 0.50000000000000), +0.000970017636685),
    (SVector(0.00000000000000, 0.50000000000000, 0.00000000000000), +0.000970017636685),
    (SVector(0.50000000000000, 0.00000000000000, 0.00000000000000), +0.000970017636685),
    (SVector(0.25000000000000, 0.25000000000000, 0.25000000000000), +0.018264223466167),
    (SVector(0.07821319233030, 0.07821319233030, 0.07821319233030), +0.010599941524417),
    (SVector(0.07821319233030, 0.07821319233030, 0.76536042300900), +0.010599941524417),
    (SVector(0.07821319233030, 0.76536042300900, 0.07821319233030), +0.010599941524417),
    (SVector(0.76536042300900, 0.07821319233030, 0.07821319233030), +0.010599941524417),
    (SVector(0.12184321666400, 0.12184321666400, 0.12184321666400), -0.062517740114333),
    (SVector(0.12184321666400, 0.12184321666400, 0.63447035000800), -0.062517740114333),
    (SVector(0.12184321666400, 0.63447035000800, 0.12184321666400), -0.062517740114333),
    (SVector(0.63447035000800, 0.12184321666400, 0.12184321666400), -0.062517740114333),
    (SVector(0.33253916444600, 0.33253916444600, 0.33253916444600), +0.004891425263067),
    (SVector(0.33253916444600, 0.33253916444600, 0.00238250666074), +0.004891425263067),
    (SVector(0.33253916444600, 0.00238250666074, 0.33253916444600), +0.004891425263067),
    (SVector(0.00238250666074, 0.33253916444600, 0.33253916444600), +0.004891425263067),
    (SVector(0.10000000000000, 0.10000000000000, 0.20000000000000), +0.027557319224000),
    (SVector(0.10000000000000, 0.20000000000000, 0.10000000000000), +0.027557319224000),
    (SVector(0.10000000000000, 0.10000000000000, 0.60000000000000), +0.027557319224000),
    (SVector(0.10000000000000, 0.60000000000000, 0.10000000000000), +0.027557319224000),
    (SVector(0.10000000000000, 0.20000000000000, 0.60000000000000), +0.027557319224000),
    (SVector(0.10000000000000, 0.60000000000000, 0.20000000000000), +0.027557319224000),
    (SVector(0.20000000000000, 0.10000000000000, 0.10000000000000), +0.027557319224000),
    (SVector(0.20000000000000, 0.10000000000000, 0.60000000000000), +0.027557319224000),
    (SVector(0.20000000000000, 0.60000000000000, 0.10000000000000), +0.027557319224000),
    (SVector(0.60000000000000, 0.10000000000000, 0.20000000000000), +0.027557319224000),
    (SVector(0.60000000000000, 0.10000000000000, 0.10000000000000), +0.027557319224000),
    (SVector(0.60000000000000, 0.20000000000000, 0.10000000000000), +0.027557319224000))

## 4 negative weights, 0 points outside of the tetrahedron

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 8 on the tetrahedron *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TETAHEDRON_G8N43 = (
    (SVector(0.2500000000000, 0.2500000000000, 0.2500000000000), -0.020500188658667),
    (SVector(0.2068299316110, 0.2068299316110, 0.2068299316110), +0.014250305822867),
    (SVector(0.2068299316110, 0.2068299316110, 0.3795102051680), +0.014250305822867),
    (SVector(0.2068299316110, 0.3795102051680, 0.2068299316110), +0.014250305822867),
    (SVector(0.3795102051680, 0.2068299316110, 0.2068299316110), +0.014250305822867),
    (SVector(0.0821035883105, 0.0821035883105, 0.0821035883105), +0.001967033313133),
    (SVector(0.0821035883105, 0.0821035883105, 0.7536892350680), +0.001967033313133),
    (SVector(0.0821035883105, 0.7536892350680, 0.0821035883105), +0.001967033313133),
    (SVector(0.7536892350680, 0.0821035883105, 0.0821035883105), +0.001967033313133),
    (SVector(0.0057819505052, 0.0057819505052, 0.0057819505052), +0.000169834109093),
    (SVector(0.0057819505052, 0.0057819505052, 0.9826541484840), +0.000169834109093),
    (SVector(0.0057819505052, 0.9826541484840, 0.0057819505052), +0.000169834109093),
    (SVector(0.9826541484840, 0.0057819505052, 0.0057819505052), +0.000169834109093),
    (SVector(0.0505327400189, 0.0505327400189, 0.4494672599810), +0.004579683824467),
    (SVector(0.0505327400189, 0.4494672599810, 0.0505327400189), +0.004579683824467),
    (SVector(0.4494672599810, 0.0505327400189, 0.0505327400189), +0.004579683824467),
    (SVector(0.0505327400189, 0.4494672599810, 0.4494672599810), +0.004579683824467),
    (SVector(0.4494672599810, 0.0505327400189, 0.4494672599810), +0.004579683824467),
    (SVector(0.4494672599810, 0.4494672599810, 0.0505327400189), +0.004579683824467),
    (SVector(0.2290665361170, 0.2290665361170, 0.0356395827885), +0.005704485808683),
    (SVector(0.2290665361170, 0.0356395827885, 0.2290665361170), +0.005704485808683),
    (SVector(0.2290665361170, 0.2290665361170, 0.5062273449780), +0.005704485808683),
    (SVector(0.2290665361170, 0.5062273449780, 0.2290665361170), +0.005704485808683),
    (SVector(0.2290665361170, 0.0356395827885, 0.5062273449780), +0.005704485808683),
    (SVector(0.2290665361170, 0.5062273449780, 0.0356395827885), +0.005704485808683),
    (SVector(0.0356395827885, 0.2290665361170, 0.2290665361170), +0.005704485808683),
    (SVector(0.0356395827885, 0.2290665361170, 0.5062273449780), +0.005704485808683),
    (SVector(0.0356395827885, 0.5062273449780, 0.2290665361170), +0.005704485808683),
    (SVector(0.5062273449780, 0.2290665361170, 0.0356395827885), +0.005704485808683),
    (SVector(0.5062273449780, 0.2290665361170, 0.2290665361170), +0.005704485808683),
    (SVector(0.5062273449780, 0.0356395827885, 0.2290665361170), +0.005704485808683),
    (SVector(0.0366077495532, 0.0366077495532, 0.1904860419350), +0.002140519141167),
    (SVector(0.0366077495532, 0.1904860419350, 0.0366077495532), +0.002140519141167),
    (SVector(0.0366077495532, 0.0366077495532, 0.7362984589590), +0.002140519141167),
    (SVector(0.0366077495532, 0.7362984589590, 0.0366077495532), +0.002140519141167),
    (SVector(0.0366077495532, 0.1904860419350, 0.7362984589590), +0.002140519141167),
    (SVector(0.0366077495532, 0.7362984589590, 0.1904860419350), +0.002140519141167),
    (SVector(0.1904860419350, 0.0366077495532, 0.0366077495532), +0.002140519141167),
    (SVector(0.1904860419350, 0.0366077495532, 0.7362984589590), +0.002140519141167),
    (SVector(0.1904860419350, 0.7362984589590, 0.0366077495532), +0.002140519141167),
    (SVector(0.7362984589590, 0.0366077495532, 0.1904860419350), +0.002140519141167),
    (SVector(0.7362984589590, 0.0366077495532, 0.0366077495532), +0.002140519141167),
    (SVector(0.7362984589590, 0.1904860419350, 0.0366077495532), +0.002140519141167))

## 1 negative weights, 0 points outside of the tetrahedron

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
using BlockAverage
using Test

@testset "BlockAverage.jl" begin
    # Write your tests here.
end


text:
module TestBangBang

using Test
using BangBang
using InitialValues: InitialValue, hasinitialvalue

@testset begin
    @test hasinitialvalue(push!!)
    @test hasinitialvalue(append!!)
    @test push!!(InitialValue(push!!), 1) == [1]
    @test append!!(InitialValue(append!!), [1]) == [1]
    @test append!!([1], InitialValue(append!!)) == [1]
end

@testset "promote" begin
    init = InitialValue(+)
    @test eltype(push!!([init], 0)) == Union{typeof(init),Int}
end

end  # module


text:
using BinDeps

const version = "2.0.2"
const unpacked_dir = "optimpack-$version"

@BinDeps.setup

optimpack = library_dependency("libopk")

provides(Sources,
         URI("https://github.com/emmt/OptimPack/releases/download/v$version/optimpack-$version.tar.gz"),
         optimpack,
         unpacked_dir=unpacked_dir)

prefix = joinpath(BinDeps.depsdir(optimpack), "usr")
srcdir = joinpath(BinDeps.depsdir(optimpack), "src", unpacked_dir)
@unix_only libfilename = "libopk.so"
@osx_only libfilename = "libopk.dylib"
destlib = joinpath(prefix, "lib", libfilename)

provides(SimpleBuild,
         (@build_steps begin
             GetSources(optimpack)
             CreateDirectory(prefix)
             CreateDirectory(joinpath(prefix,"lib"))
             @build_steps begin
                 ChangeDirectory(srcdir)
                 FileRule(destlib,
                          @build_steps begin
                              `./configure --enable-shared --disable-static --prefix="$prefix"`
                              `make`
                              `make install`
                          end)
             end
         end),
         optimpack)

@BinDeps.install Dict(:libopk => :opklib)


text:
using BinaryBuilder

# Collection of sources required to build LibHealpixWrapper
sources = [
    "https://github.com/mweastwood/LibHealpixWrapper.git" =>
    "889645f510d7283f6ce3393eb6e2576c9d6f6be4",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd LibHealpixWrapper
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain
make
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    BinaryProvider.Linux(:i686, :glibc, :blank_abi),
    BinaryProvider.Linux(:aarch64, :glibc, :blank_abi),
    BinaryProvider.Linux(:powerpc64le, :glibc, :blank_abi),
    BinaryProvider.Linux(:armv7l, :glibc, :eabihf),
    BinaryProvider.MacOS(:x86_64, :blank_libc, :blank_abi),
    BinaryProvider.Linux(:x86_64, :glibc, :blank_abi)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libhealpixwrapper", :libhealpixwrapper)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaAstro/CFITSIOBuilder/releases/download/v3.440-1/build.jl",
    "https://github.com/mweastwood/LibHealpixBuilder/releases/download/v3.31-1/build.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "LibHealpixWrapper", sources, script, platforms, products, dependencies)

max_stars_repo_name: mweastwood/LibHealpixWrapperBuilder

max_issues_repo_name: mweastwood/LibHealpixWrapperBuilder

max_forks_repo_name: mweastwood/LibHealpixWrapperBuilder


text:
module MosaicPlutoPlugin
using HypertextLiteral

function MosaicWidget(schem; db="schematics")
    @htl("""

    """)

end

export MosaicWidget
end # module

max_stars_repo_name: JuliaComputing/MosaicPlutoPlugin

max_issues_repo_name: JuliaComputing/MosaicPlutoPlugin

max_forks_repo_name: JuliaComputing/MosaicPlutoPlugin


text:
# https://adventofcode.com/2021/day/1
using AdventOfCode
input = parse.(Int,readlines("2021/data/day_1.txt"))

function part_1(input)
    sum(>(0),diff(input))
end

function part_2(input)
    sum(>(0),input[begin+3:end].-input[begin:end-3])
end


text:
# tuples are ordered and immutable
t1 = (1, 2, 3, 4, "Richard")
println(typeof(t1))
println("the 2nd element of t1 is : ", t1[2])
println("a slice of 2nd to 4th elements of t1 is : ", t1[2:4])

# iterate over a tuple
println("*** iterate over a tuple ***")
for e in t1
    println(e)
end
println("******")

# you can check that a 'value' is 'in' a tuple
println("*** check a tuple for a value ***")
println(2 in t1)
println(9 in t1)
println("******")

# multi-dimensional tuples
println("*** multi-dimensional tuples ***")
t2 = ((1, 2), (3, 4))
println(typeof(t2))
println(t2[1][2])
println("******")

# named tuples
println("*** named tuples ***")
t3 = (a = 1, b = 2, name = "Richard")
println(typeof(t3))
println(t3)
println(t3.name)
println("******")

max_stars_repo_name: rpillar/JuliaExamples

max_issues_repo_name: rpillar/JuliaExamples

max_forks_repo_name: rpillar/JuliaExamples


text:
using Documenter
using FermionicOperatorStringAnalyzer

makedocs(sitename="FOSA Documentation")


text:
#=
    The Objective Sense
=#

MOI.canget(m::LinQuadOptimizer, ::MOI.ObjectiveSense) = true
MOI.get(m::LinQuadOptimizer,::MOI.ObjectiveSense) = m.obj_sense

MOI.canset(::LinQuadOptimizer, ::MOI.ObjectiveSense) = true
function MOI.set!(m::LinQuadOptimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    if sense == MOI.MinSense
        change_objective_sense!(m, :min)
        m.obj_sense = MOI.MinSense
    elseif sense == MOI.MaxSense
        change_objective_sense!(m, :max)
        m.obj_sense = MOI.MaxSense
    elseif sense == MOI.FeasibilitySense
        # we set the objective sense to :min, and the objective to 0.0
        change_objective_sense!(m, :min)
        unsafe_set!(m, MOI.ObjectiveFunction{Linear}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[],0.0))
        m.obj_type = AffineObjective
        m.obj_sense = MOI.FeasibilitySense
    else
        error("Sense $(sense) unknown.")
    end
end

#=
    The Objective Function
=#

function MOI.canset(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{F}) where F<:MOI.AbstractFunction
    if MOI.get(m, MOI.ObjectiveSense()) == MOI.FeasibilitySense
        # it doesn't make sense to set an objective for a feasibility problem
        return false
    end
    return F in supported_objectives(m)
end

function MOI.set!(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{MOI.SingleVariable}, var::MOI.SingleVariable)
    if m.obj_type == QuadraticObjective
        set_quadratic_objective!(m, Int[], Int[], Float64[])
    end
    m.obj_type = SingleVariableObjective
    m.single_obj_var = var.variable
    set_linear_objective!(m, [getcol(m, var.variable)], [1.0])
    m.objective_constant = 0
end

function MOI.set!(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{F}, objf::Linear) where F
    cobjf = MOIU.canonical(objf)
    unsafe_set!(m, MOI.ObjectiveFunction{Linear}(), cobjf)
end

"""
    unsafe_set!(m, ::MOI.ObjectiveFunction{F}, objective::Linear) where F

Sets a linear objective function without cannonicalizing `objective`.
"""
function unsafe_set!(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{F}, objf::Linear) where F
    if m.obj_type == QuadraticObjective
        # previous objective was quadratic...
        # zero quadratic part
        set_quadratic_objective!(m, Int[], Int[], Float64[])
    end
    m.obj_type = AffineObjective
    m.single_obj_var = nothing
    columns = [getcol(m, term.variable_index) for term in objf.terms]
    coefficients = [term.coefficient for term in objf.terms]
    set_linear_objective!(m, columns, coefficients)
    m.objective_constant = objf.constant
    nothing
end

function MOI.set!(m::LinQuadOptimizer, ::MOI.ObjectiveFunction, objf::Quad)
    m.obj_type = QuadraticObjective
    m.single_obj_var = nothing
    columns = [getcol(m, term.variable_index) for term in objf.affine_terms]
    coefficients = [term.coefficient for term in objf.affine_terms]
    set_linear_objective!(m, columns, coefficients)

    quadratic_columns_1 = [getcol(m, term.variable_index_1) for term in objf.quadratic_terms]
    quadratic_columns_2 = [getcol(m, term.variable_index_2) for term in objf.quadratic_terms]
    quadratic_coefficients = [term.coefficient for term in objf.quadratic_terms]
    ri, ci, vi = reduce_duplicates!(
        quadratic_columns_1,
        quadratic_columns_2,
        quadratic_coefficients
    )
    set_quadratic_objective!(m, ri, ci, vi)
    m.objective_constant = objf.constant
    nothing
end

#=
    Get the objective function
=#

MOI.canget(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{MOI.SingleVariable}) = m.obj_type == SingleVariableObjective
function MOI.get(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{MOI.SingleVariable})
    SingleVariable(get(m.single_obj_var))
end

MOI.canget(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{Linear}) = m.obj_type != QuadraticObjective
function MOI.get(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{Linear})
    variable_coefficients = zeros(length(m.variable_references))
    get_linear_objective!(m, variable_coefficients)
    terms = map(
        (v,c)->MOI.ScalarAffineTerm{Float64}(c,v),
        m.variable_references,
        variable_coefficients
    )
    Linear(terms, m.objective_constant)
end

MOI.canget(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{Quad}) = m.obj_type == QuadraticObjective
function MOI.get(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{Quad})
    variable_coefficients = zeros(length(m.variable_references))
    get_linear_objective!(m, variable_coefficients)
    affine_terms = map(
        (v,c)->MOI.ScalarAffineTerm{Float64}(c,v),
        m.variable_references,
        variable_coefficients
    )
    Q = get_quadratic_terms_objective(m)
    rows = rowvals(Q)
    vals = nonzeros(Q)
    nrows, ncols = size(Q)
    quadratic_terms = MOI.ScalarQuadraticTerm{Float64}[]
    sizehint!(quadratic_terms, length(vals))
    for i = 1:ncols
        for j in nzrange(Q, i)
            row = rows[j]
            val = vals[j]
            push!(quadratic_terms, MOI.ScalarQuadraticTerm{Float64}(val, m.variable_references[row], m.variable_references[i]))
        end
    end
    Quad(affine_terms, quadratic_terms, m.objective_constant)
end

#=
    Modify objective function
=#

function MOI.canmodify(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{F}, ::Type{MOI.ScalarCoefficientChange{Float64}}) where F
    if F <: MOI.ScalarQuadraticFunction
        return m.obj_type == QuadraticObjective
    elseif F <: MOI.ScalarAffineFunction
        return m.obj_type == AffineObjective
    end
    return false
end
function MOI.modify!(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{F}, chg::MOI.ScalarCoefficientChange{Float64}) where F
    if m.obj_type == SingleVariableObjective
        m.obj_type = AffineObjective
        m.single_obj_var = nothing
    end
    col = m.variable_mapping[chg.variable]
    change_objective_coefficient!(m, col, chg.new_coefficient)
end

function MOI.canmodify(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{F}, ::Type{MOI.ScalarConstantChange{Float64}}) where F
    return !(F == MOI.SingleVariable)
end
function MOI.modify!(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{F}, chg::MOI.ScalarConstantChange{Float64}) where F
    m.objective_constant = chg.new_constant
end


text:
"""
Load and plot simulations results on confidence bands' coverages
"""

#region import and models
using Pkg
Pkg.activate(".")
Pkg.instantiate()
using DrWatson
using DataFrames
using PyPlot
pygui(true)
using JLD2
using ScoreDrivenERGM
using ScoreDrivenERGM.DynNets
using Statistics
using Logging
using SharedArrays
using LinearAlgebra
using ProjUtilities

#endregion

@elapsed dfEst = collect_results( datadir("sims", "dgp&FIl_est"))
dfEst["modelTag"] = string.(dfEst["model"])

@elapsed df = dfEst

# @elapsed dfConf = collect_results( datadir("sims", "dgp&FIl_conf"))
# dfConf["modelTag"] = string.(dfConf["model"])

# @elapsed df = antijoin(dfEst, dfConf, on = [:allParDgpT, :modelTag])

df.dgpSettings

begin

# limitSample = 10
limitSample = nothing
# listParUncMethod = ["PB-MVN"] #"WHITE-MLE"
winsorProp=0
listParUncMethod = ["WHITE-MLE"]
#listParUncMethod = ["PB-MVN"] #"WHITE-MLE"
# m = listParUncMethod[1]
# estResRow = collect(eachrow(df))[1]
for estResRow in eachrow(df)

    processFlag = false
    if (estResRow.N in [100]) & (estResRow.T in [300, 3000])  &  (estResRow.model.scoreScalingType == "FISH_D") & (estResRow.dgpSettings.type == "AR")

        if  estResRow.dgpSettings.opt.B[1] == 1
            processFlag = true
            if haskey(estResRow.model.options, "integrated")
                if estResRow.model.options["integrated"]
                end
            end
        end
    end

    if processFlag
        for m in listParUncMethod

            if winsorProp!=0
                m = "$m$winsorProp"
            end

            timeConf = @elapsed avgCover, constInds, errInds, allConfBandsFiltPar, allConfBandsPar, allmvSDUnParEstCovWhite = average_coverages(estResRow, m; limitSample = limitSample, winsorProp=winsorProp)

            coverDict =  (;m, avgCover, constInds, errInds, allConfBandsFiltPar,  allConfBandsPar, allmvSDUnParEstCovWhite) |>DrWatson.ntuple2dict |> DrWatson.tostringdict

            coverDict = merge(DrWatson.tostringdict(estResRow), coverDict)
            delete!(coverDict, "path")
            delete!(coverDict, "modelTag")

            coverDict["dgp"] = coverDict["dgpSettings"]
            delete!(coverDict, "dgpSettings")
            coverDict["S"] = coverDict["nSample"]
            delete!(coverDict, "nSample")

            # loadPath = estResRow.path

            saveName = replace.( savename(coverDict, "jld2"; allowedtypes = (Real, String, Symbol, NamedTuple, Tuple, ScoreDrivenERGM.DynNets.SdErgm) ), r"[\"]" => "") #loadPath[findlast("\\", loadPath)[end]+1:end]

            timeSave = save( datadir("sims", "dgp&FIl_conf", saveName), coverDict)

            Logging.@info(" errRatio = $(mean(errInds)) , time = $timeConf")

        end
    end
end

end

# T = 100
# N = 100
# nSample = 100

# modelTag = string(DynNets.SdErgmDirBin0Rec0_pmle("HESS_D"))

# dgpSetting = DynNets.list_example_dgp_settings(DynNets.SdErgmDirBin0Rec0_mle()).dgpSetSDhigh

# res = filter([:modelTag, :T, :N, :dgpSettings, :nSample] => (m,t,n, d, s) -> all((m==modelTag, t==T, n==N, d == dgpSetting, s == nSample)), df)[1,:]

# res.allvEstSdResPar
# vEstSdUnParBootDist = mapslices(x -> DynNets.unrestrict_all_par(res.model, x), res.allvEstSdResPar, dims=1)

# # @show res.allvEstSdResPar
# covParBoot = cov(Utilities.drop_bad_un_estimates(vEstSdUnParBootDist)')

# covParBoot, minEigenVal = Utilities.make_pos_def(covParBoot)

# mvSDUnParEstCov = Symmetric(covParBoot)

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
# """
# 	print_results(ata_model; group_by_fs = false, results_folder = "results")

# # Description

# Print the features of the assembled tests.

# # Arguments

# - **`ata_model::AbstractModel`** : Required. The model built with `ATA` fuctions, `ata_model.design` matrix must be `IxT` or `nfsxT` if the items are grouped by friend sets.
# - **`group_by_fs`** : Optional. Default: `false`. Set to `true` if items have been grouped by friend sets by [`group_by_friends!`](#ATA.group_by_friends!-Tuple{ATA.Model}).
# - **`results_folder`** : Optional. Default: "results". The folder in which the output is stored.
# """
# function print_results(
#     ata_model::AbstractModel;
#     group_by_fs = false,
#     results_folder = "results",
# )
#     if !isdir(results_folder)
#         mkdir(results_folder)
#     else
#         println(
#             string(
#                 "There is already a folder with this name, files in ",
#                 results_folder,
#                 " will be overwritten."),
#             ),
#         )
#     end
#     if size(ata_model.output.design, 1) > 0
#         T = ata_model.settings.T
#         n_items = ata_model.settings.n_items
#         length_max = [ata_model.constraints[t].length_max for t = 1:T]
#         categories = ata_model.output.categories
#         n_fs = size(ata_model.settings.fs.items, 1)
#         if n_fs < ata_model.settings.n_items
#             n = n_fs * T
#         else
#             n = n_items * T
#         end

#         if ata_model.obj.name == "cc_maximin"
#             JLD2.@load "opt/IIF_CC.jld2" IIF
#             JLD2.@load "opt/ICF_CC.jld2" ICF
#             IIF_CC = copy(IIF)
#             ICF_CC = copy(ICF)
#         end
#         if ata_model.obj.name == "maximin" ||
#            ata_model.obj.name == "cc_maximin" ||
#            ata_model.obj.name == "minimax"
#             JLD2.@load "opt/IIF.jld2" IIF
#         end
#         if isfile("opt/ICF.jld2")
#             JLD2.@load "opt/ICF.jld2" ICF
#         end
#         if group_by_fs == true
#             design = reshape(ata_model.output.design, n_fs, T)
#             new_design = zeros(Float64, n_items, T)
#             for t = 1:T
#                 new_design[
#                     vcat(
#                         ata_model.settings.fs.items[findall(
#                             design[:, t] .== one(Float64),
#                         )]...,
#                     ),
#                     t,
#                 ] .= one(Float64)
#             end
#             design = new_design
#             DelimitedFiles.writedlm(string(results_folder, "/designItemLevel.csv"), design)
#         else
#             design = reshape(ata_model.output.design, n_items, T)
#         end
#         #overlap
#         olMatrixOUT = zeros(Int64, T, T)
#         for t = 1:T
#             for v = 1:T
#                 olMatrixOUT[t, v] = design[:, t]' * design[:, v]
#             end
#         end
#         #DelimitedFiles.writedlm(string(results_folder,"/olMatrixOUT.csv"),olMatrixOUT)
#         #TIF e ICF
#         if ata_model.obj.name == "maximin" ||
#            ata_model.obj.name == "cc_maximin" ||
#            ata_model.obj.name == "minimax"
#             if isfile("sim_pool.csv")
#                 sim_pool = CSV.read("sim_pool.csv", DataFrames.DataFrame)
#             else
#                 sim_pool = Float64[]
#             end
#         end
#         #save values
#         if ata_model.obj.name == "cc_maximin"
#             min = Float64[]
#             First = Float64[]
#             Median = Float64[]
#             Third = Float64[]
#             max = Float64[]
#             Alpha = Float64[]
#             estimated = Float64[]
#             if isfile("sim_pool.csv")
#                 true = Float64[]
#             end
#             Distributed.@sync for t = 1:T
#                 K = size(ata_model.obj.cores[t].points, 1)
#                 for k = 1:K
#                     IIF_t = (IIF_CC[t][k, :, :]' * design[:, t])
#                     min = vcat(min, minimum(IIF_t))
#                     First = vcat(First, StatsBase.quantile(IIF_t, [0.25]))
#                     Median = vcat(Median, StatsBase.quantile(IIF_t, [0.5]))
#                     Third = vcat(Third, StatsBase.quantile(IIF_t, [0.75]))
#                     max = vcat(max, maximum(IIF_t))
#                     Alpha =
#                         vcat(Alpha, StatsBase.quantile(IIF_t, ata_model.obj.cores[t].alpha))
#                     estimated = vcat(estimated, IIF[t][k, :]' * design[:, t])
#                     if isfile("sim_pool.csv")
#                         true = vcat(
#                             true,
#                             item_info(
#                                 sim_pool,
#                                 ata_model.obj.cores[t].points[k],
#                                 model = ata_model.settings.irt.model,
#                                 parametrization = ata_model.settings.irt.parametrization,
#                                 D = ata_model.settings.irt.D,
#                             )' * design[:, t],
#                         )
#                     end
#                 end
#             end
#             tif_at_theta_pts =
#                 DataFrames.DataFrame(hcat(min, First, Median, Third, max, Alpha, estimated))

#             if isfile("sim_pool.csv")
#                 tif_at_theta_pts[:true] = true
#                 DataFrames.rename(
#                     tif_at_theta_pts,
#                     Symbol.([
#                         "min",
#                         "0.25-Qle",
#                         "Median",
#                         "0.75-Qle",
#                         "max",
#                         string(ata_model.obj.cores[1].alpha, "-Qle"),
#                         "estimated",
#                         "true",
#                     ]),
#                 )
#             else
#                 DataFrames.rename(
#                     tif_at_theta_pts,
#                     Symbol.([
#                         "min",
#                         "0.25-Qle",
#                         "Median",
#                         "0.75-Qle",
#                         "max",
#                         string(ata_model.obj.cores[1].alpha, "-Qle"),
#                         "estimated",
#                     ]),
#                 )
#             end
#             CSV.write(string(results_folder, "/tif_at_theta_pts.csv"), tif_at_theta_pts)

#         elseif ata_model.obj.name == "maximin" || ata_model.obj.name == "minimax"
#             estimated = Float64[]
#             if isfile("sim_pool.csv")
#                 true = Float64[]
#             end
#             for t = 1:T
#                 K = size(ata_model.obj.cores[t].points, 1)
#                 for k = 1:K
#                     estimated = vcat(estimated, IIF[t][k, :]' * design[:, t])
#                     if isfile("sim_pool.csv")
#                         true = vcat(
#                             true,
#                             item_info(
#                                 sim_pool,
#                                 ata_model.obj.cores[t].points[k],
#                                 model = ata_model.settings.irt.model,
#                                 parametrization = ata_model.settings.irt.parametrization,
#                                 D = ata_model.settings.irt.D,
#                             )' * design[:, t],
#                         )
#                     end
#                 end
#             end
#             tif_at_theta_pts = DataFrames.DataFrame(estimated = estimated)
#             if isfile("sim_pool.csv")
#                 tif_at_theta_pts[:true] = true
#                 DataFrames.rename(tif_at_theta_pts, Symbol.(["estimated", "true"]))
#             end
#             CSV.write(string(results_folder, "/tif_at_theta_pts.csv"), tif_at_theta_pts)
#         end

#         #expected score
#         if isfile("opt/ICF.jld2")
#             es_print_irt = Vector{Vector{Float64}}(undef, T)
#             for t = 1:T
#                 es_print_irt[t] = zeros(size(ICF[t], 1))
#                 for k = 1:size(ICF[t], 1)
#                     es_print_irt[t][k] = ((ICF[t][k, :]'*design[:, t]))[1] / sum(design[:, t])
#                 end
#             end
#         end

#         #number of items
#         n = sum(design, dims = 1)
#         design_items = zeros(Int64, Int(maximum(n)), T)
#         #design_items=Array{String,2}(nothing,n_max,T)
#         for t = 1:T
#             design_items[collect(1:Int(n[t])), t] .= sort!(findall(design[:, t] .== 1))#sort!(CODEVar[findall(design[:,t].==1)])
#         end
#         open(string(results_folder, "/results.txt"), "w") do io
#             write(io, "Length")
#             write(io, "\r\n")
#             DelimitedFiles.writedlm(io, [1:T, Int.(n)'], "\t")
#             write(io, "\r\n")
#             write(io, "design")
#             write(io, "\r\n")
#             DelimitedFiles.writedlm(io, design_items, "\t")
#             write(io, "\r\n")
#             write(io, "Expected score")
#             write(io, "\r\n")
#             if isfile("opt/ICF.jld2")
#                 for t = 1:T
#                     DelimitedFiles.writedlm(
#                         io,
#                         vcat(t, round.(es_print_irt[t], digits = 3))',
#                         "\t",
#                     )
#                     write(io, "\r\n")
#                 end
#             end
#             if size(categories, 1) > 0
#                 write(io, "Categorical variables")
#                 write(io, "\r\n")
#                 for cats in categories
#                     write(io, cats)
#                     write(io, "\r\n")
#                     vals1 =
#                         setdiff(unique(ata_model.settings.bank[!, Symbol(cats)]), [missing])
#                     vals = copy(vals1)
#                     for t = 1:T
#                         valscat = Int64[]
#                         for val in vals1
#                             valscat_t =
#                                 ata_model.settings.bank[!, Symbol(cats)][design[
#                                     :,
#                                     t,
#                                 ].==1] .== val
#                             valscat_t = valscat_t[.!ismissing.(valscat_t)]
#                             if size(valscat_t, 1) > 0
#                                 valscat = vcat(valscat, sum(valscat_t))
#                             else
#                                 valscat = vcat(valscat, 0)
#                             end
#                         end
#                         vals = hcat(vals, valscat)
#                     end
#                     DelimitedFiles.writedlm(io, vals, "\t")
#                     #DelimitedFiles.writedlm(io,hcat(vals,[sum(skipmissing(ata_model.settings.bank[Symbol(cats)][design[:,t].==1].== i))   for i in vals,t=1:T]),", ")
#                     write(io, "\r\n")
#                 end
#                 write(io, "Sum variables")
#                 write(io, "\r\n")
#                 for var in unique(vcat([ata_model.constraints[t].sum_vars for t = 1:T]...))
#                     write(io, var)
#                     write(io, "\r\n")
#                     DelimitedFiles.writedlm(
#                         io,
#                         round.(
#                             [
#                                 sum(
#                                     skipmissing(
#                                         ata_model.settings.bank[!, var] .* design[:, t],
#                                     ),
#                                 ) for t = 1:T
#                             ],
#                             digits = 3,
#                         )',
#                         "\t",
#                     )
#                     write(io, "\r\n")
#                 end
#                 write(io, "\r\n")
#                 write(io, "Overlap matrix")
#                 write(io, "\r\n")
#                 olMatrixOUT =
#                     hcat(vcat("", string.(collect(1:T))), vcat(collect(1:T)', olMatrixOUT))
#                 DelimitedFiles.writedlm(io, olMatrixOUT, "\t")
#                 write(io, "\r\n")
#             end
#             write(io, "Item use")
#             write(io, "\r\n")
#             iu = [sum(design[i, :]) for i = 1:ata_model.settings.n_items]
#             ids = sortperm(iu; rev = true)
#             sort!(iu; rev = true)
#             for i = 1:ata_model.settings.n_items
#                 write(io, string(ids[i]))
#                 write(io, "\t")
#                 write(io, string(iu[i]))
#                 write(io, "\r\n")
#             end
#         end
#     else
#         println("No solution found")
#     end
# end


text:
############################################################################################
"""
$(TYPEDEF)

SMC kernel containing particles, dynamics and buffers for sampling.

# Fields
$(TYPEDFIELDS)
"""
struct SMCParticles{B<:AbstractAlgorithm,A<:ModelWrapper,D<:SMCBuffer}
    "SMC particles used for sampling."
    model::Vector{A}
    "Individual kernels associated to an SMC particle."
    kernel::Vector{B}
    "Weights struct for all SMC particles."
    weights::BaytesCore.ParameterWeights
    "Buffer values for allocation free evaluations."
    buffer::D
    function SMCParticles(
        model::Vector{A},
        algorithm::Vector{B},
        weights::BaytesCore.ParameterWeights,
        buffer::D,
    ) where {B<:AbstractAlgorithm,A<:ModelWrapper,D<:SMCBuffer}
        return new{B,A,D}(model, algorithm, weights, buffer)
    end
end

############################################################################################
"""
$(SIGNATURES)
Initialize SMC particles.

# Examples
```julia
```

"""
function SMCParticles(
    _rng::Random.AbstractRNG,
    JitterKernel::AbstractConstructor,
    objective::Objective,
    info::BaytesCore.SampleDefault,
    tune::SMCTune
)
    ## Initialize individual MCMC algorithm
    @unpack data, temperature = objective
    @unpack tagged, Ntuning = tune
    Nchains = tune.chains.Nchains
    ## Adjust observed data to current knowledge
    modelᵗᵉᵐᵖ = deepcopy(objective.model)
    ## Initialize individual Models ~ Assign copies for  Model and Algorithms ~ necessary in initiation so no pointer error
    modelᵛ = [
        ModelWrapper(deepcopy(modelᵗᵉᵐᵖ.val), modelᵗᵉᵐᵖ.info, modelᵗᵉᵐᵖ.id) for
        _ in Base.OneTo(Nchains)
    ]
    algorithmᵛ = [JitterKernel(_rng, modelᵛ[iter], data, temperature, info) for iter in Base.OneTo(Nchains)]
    ## Assign weights
    weights = BaytesCore.ParameterWeights(Nchains)
    ## Assign buffer for inplace resampling
    buffer = SMCBuffer(
        _rng, algorithmᵛ[1], modelᵗᵉᵐᵖ, data, Nchains
    )
    ## Loop through all models
    #Polyester.@batch per=thread minbatch=tune.batchsize
    Base.Threads.@threads for iter in eachindex(algorithmᵛ)
        ## Tune kernel
        propose!(_rng, algorithmᵛ[iter], modelᵛ[iter], data, temperature, BaytesCore.UpdateTrue())
        for _ in 2:Ntuning
            propose!(_rng, algorithmᵛ[iter], modelᵛ[iter], data, temperature, tune.capture)
        end
        ## Assign weight
        #!NOTE: Objective itself only stores logposterior function, no problem to call it as long as θ comes from modelᵛ.
        _objective = Objective(modelᵛ[iter], objective.data, tune.tagged, objective.temperature)
        buffer.cumweights[iter], weights.buffer[iter] = SMCweight(
            _rng,
            _objective,
            algorithmᵛ[iter],
            _objective.temperature * _objective(objective.model.val)
        )
    end
    ## Normalize weights
    weights(weights.buffer)
    ## Return container
    return SMCParticles(modelᵛ, algorithmᵛ, weights, buffer)
end

############################################################################################
function SMCDiagnostics(
    smc::SMCParticles, temperature::T, ESS::Float64, accepted::Bool, jittersteps::Int64, iter::Int64, generated::G
) where {T<:AbstractFloat, G}
    return SMCDiagnostics(
        BaytesCore.BaseDiagnostics(
            Statistics.mean(smc.buffer.cumweights),
            temperature,
            copy(smc.buffer.predictions),
            iter-1
        ),
        BaytesCore.weightedincrement(smc.weights),
        #!NOTE: There is not really any way around making a copy from buffer if no pointer issue for diagnostics
        copy(smc.buffer.cumweights),
        copy(smc.weights.ℓweightsₙ),
        deepcopy(smc.buffer.jitterdiagnostics),
        jittersteps,
        copy(smc.buffer.correlation.ρ),
        ESS,
        accepted,
        generated
    )
end

############################################################################################
#export
export SMCParticles, init, resample!, propagate!, weight!


text:
using Documenter, FeatureSelectors

makedocs(;
    modules = [FeatureSelectors],
    format = Documenter.HTML(prettyurls = false),
    pages = ["Home" => "index.md", "Module" => "module.md"],
    repo = "https://github.com/darrencl/FeatureSelectors.jl/blob/{commit}{path}#L{line}",
    sitename = "FeatureSelectors.jl",
    authors = "Darren Lukas <darrenc2995@gmail.com>",
    clean = true,
    doctest = true,
)

deploydocs(
    repo = "github.com/darrencl/FeatureSelectors.jl.git",
    target = "build",
    push_preview = true,
    devbranch = "master",
)


text:
module SimulacoesTemplate

  include("Data.jl")
  export Data

  include("./initial-point.jl")
  include("./initial-velocities.jl")
  export initial_point, initial_velocities

  include("./rpbc.jl")
  include("./upair.jl")
  include("./utotal.jl")
  export utotal

  include("./linear_index.jl")
  export linear_index

  include("./utotal_parallel.jl")
  export utotal
  export utotal_parallel1
  export utotal_parallel2

  include("./wrap.jl")
  include("./forcepair.jl")
  include("./energyandforce.jl")
  export energyandforce!

  include("./simulate.jl")
  export simulate

end