text:
using Tmatrix
using Zygote

rx, rz = 1e-6, 1e-6
n_max = 1
k1_r = 1e7;
k1_i = 0.0;
k2_r = 1.5e7;
k2_i = 1e3;

k1 = Complex(k1_r, k1_i)
k2 = Complex(k2_r, k2_i)

n_θ_points = 10
θ_array = collect(LinRange(1e-6, pi, n_θ_points))
r_array = rx * ones(size(θ_array))
ϕ_array = zeros(size(θ_array))
rotationally_symmetric = true
symmetric_about_plane_perpendicular_z = false
BigFloat_precision = nothing

# testing out T-matrix calculation and automatic differentiation using <<T_matrix_SeparateRealImag_arbitrary_mesh>>
@time T_sph = calculate_Tmatrix_for_spheroid_SeparateRealImag(
    rx,
    rz,
    n_max,
    k1_r,
    k1_i,
    k2_r,
    k2_i;
    n_θ_points = n_θ_points,
    rotationally_symmetric = true,
    symmetric_about_plane_perpendicular_z = false,
);
println("T-matrix using spheroid <<calculate_Tmatrix_for_spheroid_SeparateRealImag>>")
display(T_sph)

@time T = Tmatrix.T_matrix_SeparateRealImag_arbitrary_mesh(
    n_max,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    rotationally_symmetric,
    symmetric_about_plane_perpendicular_z,
    BigFloat_precision,
);
println("T-matrix using an arbitrary mesh <<T_matrix_SeparateRealImag_arbitrary_mesh>>")
display(T)

@time ∂T = Zygote.jacobian(
    Tmatrix.T_matrix_SeparateRealImag_arbitrary_mesh,
    n_max,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    rotationally_symmetric,
    symmetric_about_plane_perpendicular_z,
    BigFloat_precision,
)
println("Jacobian for <<T_matrix_SeparateRealImag_arbitrary_mesh>> calcuated by Zygote:")
display(∂T)

# testing gradient descent optimization for <<sin>>
function target_to_be_minimized(t)
    return sin(t)
end

function ∂target_to_be_minimized(t)
    return Zygote.gradient(target_to_be_minimized, t)[1]
end

"""
learning_rate = 0.01
loss_here = 1e6
t_0 = 0.3
∂loss = 0
while loss_here > 0.01
    t_0 -= learning_rate * ∂loss
    loss_here = target_to_be_minimized(t_0)
    ∂loss = ∂target_to_be_minimized(t_0)
    println("loss_here = loss_here, ∂loss = ∂loss")
end
"""

# testing gradient descent optimization for <<T_matrix_SeparateRealImag_arbitrary_mesh>>
# TODO: use an optimization package instead, so that it has better stop criteria and tricks for faster convergence
rx, rz = 1e-6, 1e-6
n_max = 2
k1_r = 1e7;
k1_i = 0.0;
k2_r = 1.5e7;
k2_i = 1e3;

n_θ_points = 20
θ_array = collect(LinRange(1e-6, pi, n_θ_points))
r_array = rx * ones(size(θ_array))
ϕ_array = zeros(size(θ_array))
rotationally_symmetric = true
symmetric_about_plane_perpendicular_z = false
BigFloat_precision = nothing

function objective_function_normal_mesh_density(r_array)
    # println("I am <<objective_function>>, and the type of r_array is $(typeof(r_array))")
    # println("r_array = $r_array")
    ϕ_array = zeros(size(θ_array))
    T = Tmatrix.T_matrix_SeparateRealImag_arbitrary_mesh(
        n_max,
        k1_r,
        k1_i,
        k2_r,
        k2_i,
        (collect(r_array)),
        θ_array,
        ϕ_array,
        rotationally_symmetric,
        symmetric_about_plane_perpendicular_z,
        BigFloat_precision,
    )
    # return abs(T[1]) # if we need to minimize the first element of T-matrix

    # return sum(abs.(T)) # if we need to minimize the sum of all real and imag

    # return Tmatrix.get_OrentationAv_scattering_CrossSection_from_Tmatrix(T, Complex(k1_r, k1_i)) # if we need to minimize scattering cross section

    return 1 * Tmatrix.get_OrentationAv_emissivity_from_Tmatrix(
        T,
        Complex(k1_r, k1_i),
        Tmatrix.calculate_surface_area_of_axisymmetric_particle(r_array, θ_array),
    ) # if we need to minimize emissivity
end

function objective_function_more_mesh_density(r_array_input)
    r_θ_array = Tmatrix.double_mesh_density(r_array_input, θ_array)
    r_array_double_density = r_θ_array[:, 1]
    θ_array_double_density = r_θ_array[:, 2]
    ϕ_array_double_density = zeros(size(θ_array_double_density))

    # println("I am <<objective_function>>, and the type of r_array_double_density is $(typeof(r_array_double_density))")
    println("r_array_input = $r_array_input")
    println("r_array_double_density = $r_array_double_density")

    T = Tmatrix.T_matrix_SeparateRealImag_arbitrary_mesh(
        n_max,
        k1_r,
        k1_i,
        k2_r,
        k2_i,
        (collect(r_array_double_density)),
        θ_array_double_density,
        ϕ_array_double_density,
        rotationally_symmetric,
        symmetric_about_plane_perpendicular_z,
        BigFloat_precision,
    )
    # return abs(T[1]) # if we need to minimize the first element of T-matrix

    # return sum(abs.(T)) # if we need to minimize the sum of all real and imag

    # return Tmatrix.get_OrentationAv_scattering_CrossSection_from_Tmatrix(T, Complex(k1_r, k1_i)) # if we need to minimize scattering cross section

    return 1 * Tmatrix.get_OrentationAv_emissivity_from_Tmatrix(
        T,
        Complex(k1_r, k1_i),
        Tmatrix.calculate_surface_area_of_axisymmetric_particle(
            r_array_double_density,
            θ_array_double_density,
        ),
    )
end

function objective_function_more_mesh_density_using_Eps(r_array_input)
    if mesh_density_boost == "double"
        r_θ_array = Tmatrix.double_mesh_density(r_array_input, θ_array)
    elseif mesh_density_boost == "quadruple"
        r_θ_array = Tmatrix.quadruple_mesh_density(r_array_input, θ_array)
    elseif mesh_density_boost == "octuple"
        r_θ_array = Tmatrix.octuple_mesh_density(r_array_input, θ_array)
    elseif mesh_density_boost == "sexdecuple"
        r_θ_array = Tmatrix.sexdecuple_mesh_density(r_array_input, θ_array)
    else
        r_θ_array = r_array_input
    end

    r_array_double_density = r_θ_array[:, 1]
    θ_array_double_density = r_θ_array[:, 2]
    ϕ_array_double_density = zeros(size(θ_array_double_density))

    # println("I am <<objective_function>>, and the type of r_array_double_density is $(typeof(r_array_double_density))")
    println("r_array_input = $r_array_input")
    println("r_array_double_density = $r_array_double_density")

    T = Tmatrix.T_matrix_SeparateRealImag_arbitrary_mesh(
        n_max,
        wl_or_freq_input,
        input_unit,
        Eps_r_r_1,
        Eps_r_i_1,
        Mu_r_r_1,
        Mu_r_i_1,
        Eps_r_r_2,
        Eps_r_i_2,
        Mu_r_r_2,
        Mu_r_i_2,
        (collect(r_array_double_density)),
        θ_array_double_density,
        ϕ_array_double_density,
        rotationally_symmetric,
        symmetric_about_plane_perpendicular_z,
        BigFloat_precision,
    )
    # return abs(T[1]) # if we need to minimize the first element of T-matrix

    # return sum(abs.(T)) # if we need to minimize the sum of all real and imag

    # return Tmatrix.get_OrentationAv_scattering_CrossSection_from_Tmatrix(T, Complex(k1_r, k1_i)) # if we need to minimize scattering cross section
    if optical_property_tobe_optimized == "emissivity"
        k1_complex = Tmatrix.get_WaveVector(
            wl_or_freq_input;
            input_unit = input_unit,
            Eps_r = Complex(Eps_r_r_1, Eps_r_i_1),
            Mu_r = Complex(Mu_r_r_1, Mu_r_i_1),
        )

        return 1 * Tmatrix.get_OrentationAv_emissivity_from_Tmatrix(
            T,
            k1_complex,
            Tmatrix.calculate_surface_area_of_axisymmetric_particle(
                r_array_double_density,
                θ_array_double_density,
            ),
        ) # if we need to minimize emissivity
    elseif optical_property_tobe_optimized == "scattering_CrossSection"
        return Tmatrix.get_OrentationAv_scattering_CrossSection_from_Tmatrix(
            T,
            Complex(k1_r, k1_i),
        ) # if we need to minimize scattering cross section
    end
end

# objective_function = objective_function_normal_mesh_density
# objective_function = objective_function_more_mesh_density
objective_function = objective_function_more_mesh_density_using_Eps

function ∂objective_function(r_array)
    # println("I am <<∂objective_function>>, and the type of r_array is $(typeof(r_array))")
    # println("r_array = $r_array")
    return Zygote.gradient(objective_function, (collect(r_array)))[1]
end

using Plots

wl_or_freq_input = 1e-6;
input_unit = "m";
n_max = 1
Eps_r_r_1 = 1.0;
Eps_r_i_1 = 0.0;
Eps_r_r_2 = 1.5;
Eps_r_i_2 = 0.01;
Mu_r_r_1 = 1.0;
Mu_r_i_1 = 0.0;
Mu_r_r_2 = 1.0;
Mu_r_i_2 = 0.0;

optical_property_tobe_optimized = "scattering_CrossSection" # can be any member of this list: ["emissivity", "scattering_CrossSection"]
maximize_or_minimize = "maximize"
mesh_density_boost = "quadruple" # by how many we need to add more points between each two consecutive points defining the particle surface. Can be any of ["double", "quadruple", "octuple", "sexdecuple"]
n_max = 1
n_θ_points = 5
θ_array = collect(LinRange(1e-6, pi, n_θ_points))
r_array_initial, _ = Tmatrix.ellipsoid(0.1e-7, 0.11e-7, θ_array)
# r_array_initial = 0.1e-7 * ones(size(θ_array))
r_array = r_array_initial
ϕ_array = zeros(size(θ_array))
xyz_initial = vcat(
    Tmatrix.convert_coordinates_Sph2Cart.(
        r_array_initial,
        θ_array,
        zeros(size(r_array_initial)),
    )...,
)
rotationally_symmetric = true
symmetric_about_plane_perpendicular_z = false
BigFloat_precision = nothing

learning_rate = 0.5e-1
# learning_rate = 1e-2
global loss_here = -1e6
global ∂loss = 0
println("Starting optimization for T-matrix")
global n_iteration = 0
loss_array = []
loss_old = 9999999999999

directory_name_to_save_plots = "cache/iteration_particle_plots/$(maximize_or_minimize)_$(optical_property_tobe_optimized)_07_27_2021_$(mesh_density_boost)_mesh_density_learning_rate_$(learning_rate)__D_0.2e-7_nmax_$(n_max)"
mkpath(directory_name_to_save_plots)

while (n_iteration < 500) # (loss_here < 0)
    n_iteration += 1

    if maximize_or_minimize == "maximize"
        r_array = r_array .+ learning_rate .* ∂loss
    else
        r_array = r_array .- learning_rate .* ∂loss
    end

    loss_here = objective_function(r_array)
    ∂loss = ∂objective_function(r_array)
    append!(loss_array, loss_here)

    println()
    println("iteration #$n_iteration: loss_here = $loss_here, ∂loss = $∂loss")
    println("r_array = $r_array")

    println()
    println("loss_array=$loss_array")
    xyz = vcat(
        Tmatrix.convert_coordinates_Sph2Cart.(r_array, θ_array, zeros(size(r_array)))...,
    )
    p1 = plot(
        xyz_initial[:, 1],
        xyz_initial[:, 3],
        aspect_ratio = :equal,
        label = "initial particle",
    )
    p1 = plot!(
        xyz[:, 1],
        xyz[:, 3],
        aspect_ratio = :equal,
        label = "particle after $n_iteration iterations",
        title = "scattering cross section = $loss_here m",
    )
    p2 = plot(
        1:length(loss_array),
        loss_array,
        xlabel = "iteration #",
        ylabel = "scattering cross section (m)",
    )
    fig = plot(p1, p2, layout = (1, 2), size = (1200, 800))

    figure_name =
        joinpath(directory_name_to_save_plots, "particle_geom_iteration_$n_iteration.png")

    savefig(fig, figure_name)
    # if (loss_here - loss_old) < 1e-6; break; end
    loss_old = loss_here
end

# # trying JuMP with cos
# using JuMP
# using Ipopt
# model = Model(Ipopt.Optimizer)
# @variable(model, x >= 0)
# # @NLconstraint(model, sin(x)>=0)
# g(x::R) where {R <: Real} = cos(x)
# register(model, :g, 1, g; autodiff = true) # try registerNLFunction
# @NLobjective(model, Min, g(x))
# optimize!(model)
#
# ####################################################################
# # trying JuMP with objective_function
# function objective_function__(r_array...)
#     return objective_function(collect((r_array)))
# end
# function ∂objective_function__(r_array...)
#     return ∂objective_function(collect((r_array)))
# end
#
# using JuMP
# using Ipopt
#
# n_θ_points = 10
# θ_array = collect(LinRange(1e-6, pi, n_θ_points))
# r_array = rx * ones(size(θ_array))
# ϕ_array = zeros(size(θ_array))
# rotationally_symmetric = true
# symmetric_about_plane_perpendicular_z = false
# BigFloat_precision = nothing
# function zz()
#     model = Model(Ipopt.Optimizer)
#     # @variable(model, r[1:n_θ_points], lowerbound=zeros(n_θ_points), upperbound=1e-6*ones(n_θ_points))
#     @variable(model, r[1:n_θ_points])
#     for r_here in r
#         JuMP.set_lower_bound(r_here, 1e-7)
#         JuMP.set_upper_bound(r_here, 1e-6)
#     end
#     set_start_value.(r, 0.5e-7 * ones(n_θ_points))
#     # @NLconstraint(model, sin(x)>=0)
#     # g(x::R) where R <: Real = cos(x)
#     register(
#         model,
#         :objective_function__,
#         n_θ_points,
#         objective_function__,
#         ∂objective_function__;
#         autodiff = true,
#     )
#     @NLobjective(model, Min, objective_function__(r...))
#     optimize!(model)
#
#     return
# end
# zz()
#
# using Interpolations
# θ_array = collect((0:0.25:1) * pi)
# r_array_initial = 1e-8 .* ones(size(θ_array))
# θ_array_interp = collect(LinRange(0, pi, 10))
# function try_(r_array)
#     """
#     sitp = scale(
#         interpolate(r_array, BSpline(Quadratic(Reflect(OnCell())))),
#         θ_array
#     );
#     r_array_interp = (sitp(θ_array_interp))
#     """
#
#     """
#     III = BSplineInterpolation(r_array,θ_array,3,:ArcLen,:Average)
#     r_array_interp = III.(θ_array_interp)
#     """
#     return sum(r_array_interp)
# end
#
# function interpolate_curve(r_array)
#     sitp = scale(interpolate(r_array, BSpline(Quadratic(Reflect(OnCell())))), θ_array)
#     return (sitp(θ_array_interp))
# end


text:
using DitherPunk
using DitherPunk: gradient_image
using ReferenceTests

using ImageCore
using ImageCore: GenericGrayImage
using ImageInTerminal
using UnicodePlots

w = 200
h = 4 * 4 # multiple of 4 for unicode braille print
img, srgb = gradient_image(h, w)
println("Test image:")
imshow(srgb)
println()

## Run reference tests for deterministic algorithms
# using Dict for Julia 1.0 compatibility
algs_deterministic = Dict(
    # threshold methods
    "ConstantThreshold" => ConstantThreshold(),
    "ClosestColor" => ClosestColor(),
    # ordered dithering
    "Bayer" => Bayer(),
    "Bayer_l2" => Bayer(; level=2),
    "Bayer_l3" => Bayer(; level=3),
    "Bayer_l4" => Bayer(; level=4),
    "ClusteredDots" => ClusteredDots(),
    "CentralWhitePoint" => CentralWhitePoint(),
    "BalancedCenteredPoint" => BalancedCenteredPoint(),
    "Rhombus" => Rhombus(),
    "IM_checks" => IM_checks(),
    "IM_h4x4a" => IM_h4x4a(),
    "IM_h6x6a" => IM_h6x6a(),
    "IM_h8x8a" => IM_h8x8a(),
    "IM_h4x4o" => IM_h4x4o(),
    "IM_h6x6o" => IM_h6x6o(),
    "IM_h8x8o" => IM_h8x8o(),
    "IM_c5x5" => IM_c5x5(),
    "IM_c6x6" => IM_c6x6(),
    "IM_c7x7" => IM_c7x7(),
    # error error_diffusion
    "SimpleErrorDiffusion" => SimpleErrorDiffusion(),
    "FloydSteinberg" => FloydSteinberg(),
    "JarvisJudice" => JarvisJudice(),
    "Stucki" => Stucki(),
    "Burkes" => Burkes(),
    "Atkinson" => Atkinson(),
    "Sierra" => Sierra(),
    "TwoRowSierra" => TwoRowSierra(),
    "SierraLite" => SierraLite(),
    "Fan" => Fan93(),
    "ShiauFan" => ShiauFan(),
    "ShiauFan2" => ShiauFan2(),
    "FalseFloydSteinberg" => DitherPunk.FalseFloydSteinberg(),
    # Keyword arguments
    "FloydSteinberg_clamp_error" => FloydSteinberg(; clamp_error=false),
    "Bayer_invert_map" => Bayer(; invert_map=true),
)

for (name, alg) in algs_deterministic
    local img2 = copy(img)
    local d = dither(img2, alg)
    @test_reference "references/gradient/$(name).txt" Int.(d)
    @test eltype(d) == eltype(img)
    @test img2 == img # image not modified

    show(brailleprint(d; title=name)) # Visualize in terminal
    println()
end

## Algorithms with random output are currently only tested visually
algs_random = Dict(
    # threshold methods
    "white_noise_dithering" => WhiteNoiseThreshold(),
)

for (name, alg) in algs_random
    local img2 = copy(img)
    local d = dither(Gray{Bool}, img2, alg)
    @test eltype(d) == Gray{Bool}
    @test img2 == img # image not modified

    show(brailleprint(d; title=name)) # Visualize in terminal
    println()
end

## Test to_linear
img2 = copy(img)
alg = Bayer()
d = dither(img2, alg; to_linear=true)
@test_reference "references/gradient/Bayer_linear.txt" Int.(d)

## Test API
d = dither(img2, alg)

# Test setting output type
d2 = dither(Gray{Float16}, img2, alg)
@test eltype(d2) == Gray{Float16}
@test d2 == d
@test img2 == img # image not modified

# Inplace modify output image
out = zeros(Bool, size(img2)...)
d3 = dither!(out, img2, alg)
@test out == d # image updated in-place
@test d3 == d
@test eltype(out) == Bool
@test eltype(d3) == Bool
@test img2 == img # image not modified

# Inplace modify  image
d4 = dither!(img2, alg)
@test d4 == d
@test img2 == d # image updated in-place
@test eltype(d4) == eltype(img)
@test eltype(img2) == eltype(img)

ext: jl


text:
module JulieTest

include("JulieTest/context.jl")
include("JulieTest/is.jl")

export
  # reporter
  typeReport,

  #context
  describe, _describe, it, iit, _it,

  # @is
  @is,
  falsy, falsey, truthy, not,

  # comparision
  least, atleast, atLeast,
  most, ismost, isMost,
  above, isabove, isAbove,
  below, isbelow, isBelow,

  # array
  empty,

  # macros
  @p, @l, @d, @pp, @s

reporterLoaded = false

function runTests()
  global reporterLoaded
  reporterLoaded || reporter("spec")
  run()
end

function reporter(str::String)
  global reporterLoaded = true
  str == "" && return
  filepath = joinpath(Pkg.dir(), "JulieTest/src/reporter", str * ".jl")
  fliepath = isfile(filepath) ? filepath : abspath(str)
  JulieTest.eval(:(include($filepath)))
end

end # module test


text:
struct SpanningTreeInstance{T <: Real, O <: CombinatorialObjective} <:
       CombinatorialInstance
    graph::AbstractGraph{T}
    rewards::Dict{Edge{T}, Float64}
    objective::O

    function SpanningTreeInstance(
        graph::AbstractGraph{T},
        rewards::Dict{Edge{T}, Float64},
        objective::O=Maximise(),
    ) where {T <: Real, O <: CombinatorialObjective}
        return new{T, O}(graph, rewards, objective)
    end
end

function dimension(i::SpanningTreeInstance{T, O}) where {T, O}
    return ne(i.graph)
end

function copy(
    i::SpanningTreeInstance{T, O};
    graph::AbstractGraph{T}=i.graph,
    rewards::Dict{Edge{T}, Float64}=i.rewards,
    objective::CombinatorialObjective=i.objective,
) where {T, O <: CombinatorialObjective}
    return SpanningTreeInstance(graph, rewards, objective)
end

function reward(i::SpanningTreeInstance{T}, e::Edge{T}) where {T}
    if e in keys(i.rewards)
        return i.rewards[e]
    end
    return i.rewards[reverse(e)]
end

fastest_exact(::SpanningTreeInstance) = PrimAlgorithm()

# Solution.

struct SpanningTreeSolution{T, O} <: CombinatorialSolution
    instance::SpanningTreeInstance{T, O}
    variables::Vector{Edge{T}}
end

function make_solution(
    i::SpanningTreeInstance{T, O},
    tree::Dict{Edge{T}, Float64},
) where {T, O}
    tree_edges = Edge{T}[]
    for (k, v) in tree
        if v >= 0.5
            push!(tree_edges, k)
        end
    end

    return SpanningTreeSolution(i, tree_edges)
end

# Budgeted solution.

struct SimpleBudgetedSpanningTreeSolution{T, U} <: SingleMinBudgetedSolution
    instance::MinimumBudget{SpanningTreeInstance{T, Maximise}, U}
    variables::Vector{Edge{T}}

    function SimpleBudgetedSpanningTreeSolution(
        instance::MinimumBudget{SpanningTreeInstance{T, Maximise}, U},
        variables::Vector{Edge{T}},
    ) where {T, U}
        return new{T, U}(instance, variables)
    end

    function SimpleBudgetedSpanningTreeSolution(
        instance::MinimumBudget{SpanningTreeInstance{T, Maximise}, U},
    ) where {T, U}
        # No feasible solution.
        return new{T, U}(instance, edgetype(instance.instance.graph)[])
    end
end

struct BudgetedSpanningTreeLagrangianSolution{T, U} <: SingleMinBudgetedSolution
    # Used to store important temporary results from solving the Lagrangian dual.
    instance::MinimumBudget{SpanningTreeInstance{T, Maximise}, U}
    variables::Vector{Edge{T}}
    λ::Float64 # Optimum dual multiplier.
    value::Float64 # Optimum value of the dual problem (i.e. with penalised constraint).
    λmax::Float64 # No dual value higher than this is useful (i.e. they all yield the same solution).
end

struct BudgetedSpanningTreeDynamicProgrammingSolution{T, U} <:
       MultipleMinBudgetedSolution
    # Used to store important temporary results from dynamic programming.
    instance::MinimumBudget{SpanningTreeInstance{T, Maximise}, U}
    variables::Vector{Edge{T}}
    solutions::Dict{Int, Vector{Edge{T}}}
end

function make_solution(
    i::MinimumBudget{SpanningTreeInstance{T, O}, U},
    tree::Dict{Edge{T}, Float64},
) where {T, O, U}
    tree_edges = Edge{T}[]
    for (k, v) in tree
        if v >= 0.5
            push!(tree_edges, k)
        end
    end

    return SimpleBudgetedSpanningTreeSolution(i, tree_edges)
end


text:
##### Representation #####

@kwdef struct RepresentationSimpleHP
  width :: Int
  depth :: Int
  hiddenstate_shape :: Int
  use_batch_norm :: Bool = false
  batch_norm_momentum :: Float32 = 0.6f0
end

mutable struct RepresentationSimpleNetwork <: AbstractRepresentation
  gspec
  hyper
  common
end

function RepresentationNetwork(gspec::AbstractGameSpec, hyper::RepresentationSimpleHP)
  bnmom = hyper.batch_norm_momentum
  function make_dense(indim, outdim)
    if hyper.use_batch_norm
      Chain(
      Dense(indim, outdim, bias=false),
      BatchNorm(outdim, relu, momentum=bnmom))
    else
      Dense(indim, outdim, relu)
    end
  end
  indim = prod(GI.state_dim(gspec))
  outdim = hyper.hiddenstate_shape
  hsize = hyper.width
  hlayers(depth) = [make_dense(hsize, hsize) for i in 1:depth]
  if hyper.depth == -1 # somewhat unintuitive, jump from 1 to 3 layers #? depth-1
    common = Chain(
      flatten,
      # make_dense(indim, outdim)
      Dense(indim, outdim)
    )
  else
    common = Chain(
      flatten,
      make_dense(indim, hsize),
      hlayers(hyper.depth)...,
      make_dense(hsize, outdim)
    )
  end
  RepresentationSimpleNetwork(gspec, hyper, common)
end

### general SimpleNet used in Prediction and Dynamics
# TODO make constructor from vector of sizes
@kwdef struct SimpleNetHP_
  indim :: Int
  outdim :: Int
  width :: Int
  depth_common :: Int
  depth_vectorhead :: Int = 1
  depth_scalarhead :: Int = 1
  use_batch_norm :: Bool = false
  batch_norm_momentum :: Float32 = 0.6f0
end

mutable struct SimpleNet_
  gspec
  hyper
  common
  vectorhead
  scalarhead
end
# TODO remove gspec dependence
function SimpleNet_(gspec::AbstractGameSpec, hyper::SimpleNetHP_)
  bnmom = hyper.batch_norm_momentum
  function make_dense(indim, outdim)
    if hyper.use_batch_norm
      Chain(
      Dense(indim, outdim, bias=false),
      BatchNorm(outdim, elu, momentum=bnmom))
    else
      Dense(indim, outdim, elu)
    end
  end
  indim = hyper.indim
  outdim = hyper.outdim
  hsize = hyper.width
  hlayers(depth) = [make_dense(hsize, hsize) for i in 1:depth] #? 1:depth-1
  # common = depth_common == -1 ?
  #   flatten :
  #   Chain(
  #     flatten, #?
  #     make_dense(indim, hsize),
  #     hlayers(hyper.depth_common)...)
  # scalarhead = Chain(
  #   hlayers(hyper.depth_scalarhead)...,
  #   Dense(hsize, 1, tanh))
  # vectorhead = Chain(
  #   hlayers(hyper.depth_vectorhead)...,
  #   Dense(hsize, outdim),
  #   softmax)
  if hyper.depth_common == -1
    common = identity
    outcomm = indim
  else
    common = Chain(
      flatten, #? identity
      make_dense(indim, hsize),
      hlayers(hyper.depth_common)...)
    outcomm = hsize
  end
  if hyper.depth_scalarhead == -1
    scalarhead = Dense(outcomm, 1, tanh)
  else
    scalarhead = Chain(
      outcomm != hsize ? make_dense(outcomm, hsize) : identity,
      hlayers(hyper.depth_scalarhead)...,
      Dense(hsize, 1, tanh))
  end
  if hyper.depth_vectorhead == -1
    vectorhead = Chain(Dense(outcomm, outdim))
  else
    vectorhead = Chain(
      outcomm != hsize ? make_dense(outcomm, hsize) : identity,
      hlayers(hyper.depth_vectorhead)...,
      Dense(hsize, outdim)
      )
  end

  SimpleNet_(gspec, hyper, common, vectorhead, scalarhead)
end

### Dynamics ###

@kwdef struct DynamicsSimpleHP
  hiddenstate_shape :: Int
  width :: Int
  depth_common :: Int
  depth_vectorhead :: Int = 1 # depth state-head
  depth_scalarhead :: Int = 1 # depth reward-head
  use_batch_norm :: Bool = false
  batch_norm_momentum :: Float32 = 0.6f0
end

mutable struct DynamicsSimpleNetwork <: AbstractDynamics
  gspec
  hyper
  common
  statehead
  rewardhead
end

function DynamicsNetwork(gspec::AbstractGameSpec, hyper::DynamicsSimpleHP)
  indim = hyper.hiddenstate_shape + GI.num_actions(gspec)
  outdim = hyper.hiddenstate_shape
  simplenethyper = SimpleNetHP_(indim,
    outdim,
    hyper.width,
    hyper.depth_common,
    hyper.depth_vectorhead,
    hyper.depth_scalarhead,
    hyper.use_batch_norm,
    hyper.batch_norm_momentum)
  simplenet = SimpleNet_(gspec,simplenethyper)
  DynamicsSimpleNetwork(gspec, hyper, simplenet.common, simplenet.vectorhead, simplenet.scalarhead)
end

### Prediction ###

@kwdef struct PredictionSimpleHP
  hiddenstate_shape :: Int
  width :: Int
  depth_common :: Int
  depth_vectorhead :: Int = 1 # depth policy-head
  depth_scalarhead :: Int = 1 # depth value-head
  use_batch_norm :: Bool = false
  batch_norm_momentum :: Float32 = 0.6f0
end

mutable struct PredictionSimpleNetwork <: AbstractPrediction
  gspec
  hyper
  common
  policyhead
  valuehead
end

function PredictionNetwork(gspec::AbstractGameSpec, hyper::PredictionSimpleHP)
  indim = hyper.hiddenstate_shape
  outdim = GI.num_actions(gspec)
  simplenethyper = SimpleNetHP_(indim,
    outdim,
    hyper.width,
    hyper.depth_common,
    hyper.depth_vectorhead,
    hyper.depth_scalarhead,
    hyper.use_batch_norm,
    hyper.batch_norm_momentum)
  simplenet = SimpleNet_(gspec,simplenethyper)
  PredictionSimpleNetwork(gspec, hyper, simplenet.common, Chain(simplenet.vectorhead, softmax), simplenet.scalarhead)
end

# function evaluate_batch(nn::AbstractPrediction, batch) # is the same as resnet


text:
using CircGeometry, Test

@testset "write_circ test" begin
    vf = 0.4
    n_bodies = 400
    material = MaterialParameters(vf,n_bodies)

    radius = 1.5
    center = Point(-0.5,1.0)
    outline = OutlineCircle(radius,center)
    translate!(outline,1.0,-1.2)

    between_buffer = compute_between_buffer(outline,material)
    ps = generate_porous_structure(outline,material,between_buffer)

    translate!(ps,1.0,-1.2)
    write_circ("test-data/test_file.circ",ps)
    save_image("test-data/test_file.svg","test-data/test_file.circ",outline)
    save_image("test-data/test_file2.svg",ps)
end

@testset "csv_to_polygon test" begin
    plist = csv_to_polygon("test-data/car.txt")
    @test length(plist) == 19
    @test plist[5].x == 0.0639999744331845
end

@testset "read_in_circ test" begin
    radiusArr, xArr, yArr = CircGeometry.read_in_circ("test-data/circle.circ")
    @test length(radiusArr) == 400
    @test xArr[5] == -0.8808863280747441
end

# clean up for next test
rm("test-data/test_file.circ")
rm("test-data/test_file.svg")
rm("test-data/test_file2.svg")

max_stars_repo_name: pseastham/CircGeometry

max_issues_repo_name: pseastham/CircGeometry


text:
@info "IncrementalInference.jl is loading tools related to InteractiveUtils.jl."

# this requires InteractiveUtils

export getCurrentWorkspaceFactors, getCurrentWorkspaceVariables
export listTypeTree

"""
    $(SIGNATURES)

Return all factors currently registered in the workspace.
"""
function getCurrentWorkspaceFactors()
    return [
        InteractiveUtils.subtypes(AbstractPrior)...,
        InteractiveUtils.subtypes(AbstractRelativeRoots)...,
        InteractiveUtils.subtypes(AbstractRelativeMinimize)...];
end

"""
    $(SIGNATURES)

Return all variables currently registered in the workspace.
"""
function getCurrentWorkspaceVariables()
    return InteractiveUtils.subtypes(IIF.InferenceVariable);
end

function _listTypeTree(mytype, printlevel::Int)
  allsubtypes = InteractiveUtils.subtypes(mytype)
  for cursubtype in allsubtypes
    print("\t"^printlevel)
    println("|___",cursubtype)
    printlevel += 1
    _listTypeTree(cursubtype, printlevel)
    printlevel -= 1
  end
end

"""
    $SIGNATURES
List the types that inherit from `T`.

Notes
- from https://youtu.be/S5R8zXJOsUQ?t=1531
"""
function listTypeTree(T)
  println(T)
  _listTypeTree(T,0)
end


text:
import Base.@deprecate_binding

function floattype(::Type{T}) where {T <: Real}
    Base.depwarn("""
        In a future release, the fallback definition of `floattype` will throw a MethodError if it cannot return a type `<:AbstractFloat`.
        See the documentation on `floattype` for guidance on whether to define a custom `floattype(::Type{$T})` method.
        """, :floattype)
    return T
end

function typechar(::Type{X}) where {X}
    Base.depwarn("""
        `typechar` was deprecated since the prefix may not be a single character in the future.
        We recommend not using private functions, but if you need to, use `type_prefix` instead.
        """, :typechar)
    Char(string(type_prefix(X))[1])
end


text:
# Test transcription_model
@testset "transcription_model" begin
    # initialize model
    m = InfiniteModel()
    # test normal usage
    @test isa(transcription_model(m), Model)
    # test error
    set_optimizer_model(m, Model())
    @test_throws ErrorException transcription_model(m)
end

# Test build_optimizer_model!
@testset "build_optimizer_model!" begin
    # initialize model
    mockoptimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                             eval_objective_value=false)
    m = InfiniteModel(mockoptimizer)
    @infinite_parameter(m, 0 <= par <= 1, num_supports = 2)
    @infinite_parameter(m, 0 <= pars[1:2] <= 1, supports = [0, 1])
    @infinite_variable(m, 1 >= x(par) >= 0, Int)
    @infinite_variable(m, y(par, pars) == 2, Bin, start = 0)
    @point_variable(m, x(0), x0)
    @point_variable(m, y(0, [0, 0]), 0 <= y0 <= 1, Int)
    @hold_variable(m, 0 <= z <= 1, Bin)
    @hold_variable(m, w == 1, Int, start = 1)
    meas1 = support_sum(x - w, par)
    meas2 = integral(y, pars)
    @constraint(m, c1, x + par - z == 0)
    @constraint(m, c2, z + x0 >= -3)
    @constraint(m, c3, meas1 + z == 0)
    @BDconstraint(m, c4(par in [0.5, 1]), meas2 - 2y0 + x <= 1)
    @constraint(m, c5, meas2 == 0)
    @constraint(m, @deriv(x, par) == 0)
    @objective(m, Min, x0 + meas1)
    set_silent(m)
    set_time_limit_sec(m, 42)
    # test normal usage
    @test isa(build_optimizer_model!(m, Val(:TransData)), Nothing)
    @test optimizer_model_ready(m)
    @test num_variables(optimizer_model(m)) == 30
    @test time_limit_sec(optimizer_model(m)) == 42
end


text:
using OffsetArrays: OffsetArray
using Oceananigans.Utils
using Oceananigans.Architectures: device_event
using Oceananigans.Grids: architecture
using KernelAbstractions.Extras.LoopInfo: @unroll

#####
##### General halo filling functions
#####

fill_halo_regions!(::Nothing, args...) = nothing

"""
    fill_halo_regions!(fields::Union{Tuple, NamedTuple}, arch, args...)

Fill halo regions for each field in the tuple `fields` according to their boundary
conditions, possibly recursing into `fields` if it is a nested tuple-of-tuples.
"""
# Some fields have `nothing` boundary conditions, such as `FunctionField` and `ZeroField`.
fill_halo_regions!(c::OffsetArray, ::Nothing, args...; kwargs...) = nothing

for dir in (:west, :east, :south, :north, :bottom, :top)
    extract_bc = Symbol(:extract_, dir, :_bc)
    @eval begin
        @inline $extract_bc(bc) = bc.$dir
        @inline $extract_bc(bc::Tuple) = $extract_bc.(bc)
    end
end

# Finally, the true fill_halo!
const MaybeTupledData = Union{OffsetArray, NTuple{<:Any, OffsetArray}}

"Fill halo regions in ``x``, ``y``, and ``z`` for a given field's data."
function fill_halo_regions!(c::MaybeTupledData, boundary_conditions, loc, grid, args...; kwargs...)

    arch = architecture(grid)

    halo_tuple = permute_boundary_conditions(boundary_conditions)

    for task = 1:3
        barrier = device_event(arch)
        fill_halo_event!(task, halo_tuple, c, loc, arch, barrier, grid, args...; kwargs...)
    end

    return nothing
end

function fill_halo_event!(task, halo_tuple, c, loc, arch, barrier, grid, args...; kwargs...)
    fill_halo! = halo_tuple[1][task]
    bc_left    = halo_tuple[2][task]
    bc_right   = halo_tuple[3][task]
    event      = fill_halo!(c, bc_left, bc_right, loc, arch, barrier, grid, args...; kwargs...)
    wait(device(arch), event)
end

function permute_boundary_conditions(boundary_conditions)

    fill_halos! = [
        fill_west_and_east_halo!,
        fill_south_and_north_halo!,
        fill_bottom_and_top_halo!,
    ]

    boundary_conditions_array_left = [
        extract_west_bc(boundary_conditions),
        extract_south_bc(boundary_conditions),
        extract_bottom_bc(boundary_conditions)
    ]

    boundary_conditions_array_right = [
        extract_east_bc(boundary_conditions),
        extract_north_bc(boundary_conditions),
        extract_top_bc(boundary_conditions),
    ]

    perm = sortperm(boundary_conditions_array_left, lt=fill_first)
    fill_halos! = fill_halos![perm]
    boundary_conditions_array_left  = boundary_conditions_array_left[perm]
    boundary_conditions_array_right = boundary_conditions_array_right[perm]

    return (fill_halos!, boundary_conditions_array_left, boundary_conditions_array_right)
end

@inline validate_event(::Nothing) = NoneEvent()
@inline validate_event(event)     = event

#####
##### Halo filling order
#####

const PBCT = Union{PBC, NTuple{<:Any, <:PBC}}
const CBCT = Union{CBC, NTuple{<:Any, <:CBC}}

fill_first(bc1::PBCT, bc2)       = false
fill_first(bc1::CBCT, bc2)       = false
fill_first(bc1::PBCT, bc2::CBCT) = false
fill_first(bc1::CBCT, bc2::PBCT) = true
fill_first(bc1, bc2::PBCT)       = true
fill_first(bc1, bc2::CBCT)       = true
fill_first(bc1::PBCT, bc2::PBCT) = true
fill_first(bc1::CBCT, bc2::CBCT) = true
fill_first(bc1, bc2)             = true

#####
##### General fill_halo! kernels
#####

@kernel function _fill_west_and_east_halo!(c, west_bc, east_bc, loc, grid, args...)
    j, k = @index(Global, NTuple)
    _fill_west_halo!(j, k, grid, c, west_bc, loc, args...)
    _fill_east_halo!(j, k, grid, c, east_bc, loc, args...)
end

@kernel function _fill_south_and_north_halo!(c, south_bc, north_bc, loc, grid, args...)
    i, k = @index(Global, NTuple)
    _fill_south_halo!(i, k, grid, c, south_bc, loc, args...)
    _fill_north_halo!(i, k, grid, c, north_bc, loc, args...)
end

@kernel function _fill_bottom_and_top_halo!(c, bottom_bc, top_bc, loc, grid, args...)
    i, j = @index(Global, NTuple)
    _fill_bottom_halo!(i, j, grid, c, bottom_bc, loc, args...)
       _fill_top_halo!(i, j, grid, c, top_bc,    loc, args...)
end

#####
##### Tuple fill_halo! kernels
#####

@kernel function _fill_west_and_east_halo!(c::NTuple, west_bc, east_bc, loc, grid, args...)
    j, k = @index(Global, NTuple)
    ntuple(Val(length(west_bc))) do n
        Base.@_inline_meta
        @inbounds begin
            _fill_west_halo!(j, k, grid, c[n], west_bc[n], loc[n], args...)
            _fill_east_halo!(j, k, grid, c[n], east_bc[n], loc[n], args...)
        end
    end
end

@kernel function _fill_south_and_north_halo!(c::NTuple, south_bc, north_bc, loc, grid, args...)
    i, k = @index(Global, NTuple)
    ntuple(Val(length(south_bc))) do n
        Base.@_inline_meta
        @inbounds begin
            _fill_south_halo!(i, k, grid, c[n], south_bc[n], loc[n], args...)
            _fill_north_halo!(i, k, grid, c[n], north_bc[n], loc[n], args...)
        end
    end
end

@kernel function _fill_bottom_and_top_halo!(c::NTuple, bottom_bc, top_bc, loc, grid, args...)
    i, j = @index(Global, NTuple)
    ntuple(Val(length(bottom_bc))) do n
        Base.@_inline_meta
        @inbounds begin
            _fill_bottom_halo!(i, j, grid, c[n], bottom_bc[n], loc[n], args...)
               _fill_top_halo!(i, j, grid, c[n], top_bc[n],    loc[n], args...)
        end
    end
end

fill_west_and_east_halo!(c, west_bc, east_bc, loc, arch, dep, grid, args...; kwargs...) =
    launch!(arch, grid, :yz, _fill_west_and_east_halo!, c, west_bc, east_bc, loc, grid, args...; dependencies=dep, kwargs...)

fill_south_and_north_halo!(c, south_bc, north_bc, loc, arch, dep, grid, args...; kwargs...) =
    launch!(arch, grid, :xz, _fill_south_and_north_halo!, c, south_bc, north_bc, loc, grid, args...; dependencies=dep, kwargs...)

fill_bottom_and_top_halo!(c, bottom_bc, top_bc, loc, arch, dep, grid, args...; kwargs...) =
    launch!(arch, grid, :xy, _fill_bottom_and_top_halo!, c, bottom_bc, top_bc, loc, grid, args...; dependencies=dep, kwargs...)


text:
struct DynamicsData{T}
	m
	ip_dyn::InteriorPoint
	ip_jac::InteriorPoint
    z_subset_init::Vector{T}
	θ_params
	h::T
	diff_idx::Int
end

typeof(view(data.ip_dyn.z, collect(1:data.m.dim.q)))
typeof(view(data.ip_dyn.δz, collect(1:data.m.dim.q), collect(1:data.m.dim.q)))

function dynamics_data(m, h,
        r!, rz!, rθ!, rz, rθ;
        idx_ineq = collect(1:0),
        idx_soc = Vector{Int}[],
        z_subset_init = ones(size(rz)[1] - m.dim.q),
		θ_params = [],
		diff_idx = -1,
		dyn_opts =  InteriorPointOptions{Float64}(
						r_tol = 1.0e-8,
						κ_tol = 1.0e-4,
						κ_init = 0.1,
						diff_sol = false),
		jac_opts =  InteriorPointOptions{Float64}(
						r_tol = 1.0e-8,
						κ_tol = 1.0e-2,
						κ_init = 0.1,
						diff_sol = true))
    nz = size(rz)[1]
    nθ = size(rθ)[2]
	z = zeros(nz)
	θ = zeros(nθ)

	ip_dyn = interior_point(z, θ,
		idx_ineq = idx_ineq,
        idx_soc = idx_soc,
		r! = r!,
		rz! = rz!,
		rθ! = rθ!,
		rz = rz,
		rθ = rθ,
		opts = dyn_opts)

	ip_dyn.opts.diff_sol = false

	ip_jac = interior_point(z, θ,
		idx_ineq = idx_ineq,
        idx_soc = idx_soc,
		r! = r!,
		rz! = rz!,
		rθ! = rθ!,
		rz = rz,
		rθ = rθ,
		opts = jac_opts)

	ip_jac.opts.diff_sol = true

	DynamicsData(m, ip_dyn, ip_jac, z_subset_init, θ_params, h, diff_idx)
end

function f!(d::DynamicsData, q0, q1, u1, mode = :dynamics)
	ip = (mode == :dynamics ? d.ip_dyn : d.ip_jac)
	h = d.h

	ip.z .= copy([q1; d.z_subset_init])
	ip.θ .= copy([q0; q1; u1; h; d.θ_params])

	status = interior_point_solve!(ip)

	!status && (@warn "dynamics failure (res norm: $(norm(ip.r, Inf))) \n
		               z = $(ip.z), \n
					   θ = $(ip.θ)")
end

function f(d::DynamicsData, q0, q1, u1)
	f!(d, q0, q1, u1, :dynamics)
	return copy(d.ip_dyn.z[1:d.m.dim.q])
end

# function fq0(d::DynamicsData, q0, q1, u1)
# 	f!(d, q0, q1, u1, :jacobian)
# 	return copy(d.ip_jac.δz[1:d.m.dim.q, 1:d.m.dim.q])
# end
#
# function fq1(d::DynamicsData, q0, q1, u1)
# 	f!(d, q0, q1, u1, :jacobian)
# 	return copy(d.ip_jac.δz[1:d.m.dim.q, d.m.dim.q .+ (1:d.m.dim.q)])
# end

function fx1(d::DynamicsData, q0, q1, u1)
	f!(d, q0, q1, u1, :jacobian)
	return copy(d.ip_jac.δz[1:d.m.dim.q, 1:(2 * d.m.dim.q)])
end

function fu1(d::DynamicsData, q0, q1, u1)
	f!(d, q0, q1, u1, :jacobian)
	return copy(d.ip_jac.δz[1:d.m.dim.q, 2 * d.m.dim.q .+ (1:d.m.dim.u)])
end

struct ImplicitDynamics{I, T} <: Model{I, T}
    n::Int
    m::Int
    d::Int
	dynamics::DynamicsData
end

function fd(model::ImplicitDynamics{Midpoint, FixedTime}, x, u, w, h, t)
	nq = model.dynamics.m.dim.q
	q0 = x[1:nq]
	q1 = x[nq .+ (1:nq)]

	q2 = f(model.dynamics, q0, q1, u)

	return [q1; q2]
end

function fdx(model::ImplicitDynamics{Midpoint, FixedTime}, x, u, w, h, t)
	nq = model.dynamics.m.dim.q
	q0 = x[1:nq]
	q1 = x[nq .+ (1:nq)]
	dq2dx1 = fx1(model.dynamics, q0, q1, u)

	return [zeros(nq, nq) I; dq2dx1]
end

function fdu(model::ImplicitDynamics{Midpoint, FixedTime}, x, u, w, h, t)
	nq = model.dynamics.m.dim.q
	q0 = x[1:nq]
	q1 = x[nq .+ (1:nq)]
	dq2du1 = fu1(model.dynamics, q0, q1, u)
	return [zeros(nq, model.m); dq2du1]
end

# test
# data = dynamics_data(model, h,
# 	dyn_opts = InteriorPointOptions{Float64}(κ_tol = 1.0e-4, κ_init = 0.1),
# 	jac_opts = InteriorPointOptions{Float64}(κ_tol = 1.0e-2, κ_init = 0.1))

# f(data, q0, q1, zeros(m))
# fq0(data, q0, q1, zeros(m))
# fq1(data, q0, q1, zeros(m))
# fx1(data, q0, q1, zeros(m))
# fu1(data, q0, q1, zeros(m))

# model_implicit = ImplicitDynamics{Midpoint, FixedTime}(2 * model.dim.q, model.dim.u, 0, data)
# fd(model_implicit, [q0; q1], zeros(model_implicit.m), zeros(model_implicit.d), h, 1)
# fdx(model_implicit, [q0; q1], zeros(model_implicit.m), zeros(model_implicit.d), h, 1)
# fdu(model_implicit, [q0; q1], zeros(model_implicit.m), zeros(model_implicit.d), h, 1)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# The purpose of this script is to compare the filtering capability of the single snapshots
# estimates and the gas filter on a given

using StatsBase, JLD,  Hwloc,  DynNets, StaticNets, StatsFuns,Utilities

# STATIC ESTIMATES INCREASING N  ---------------------------------------------
#need to choose how to distribute the unconditional means. I want to stay
# in the case of chaterjee Diaconis Sly, hence have fixed maximum value for the
# parameters that does not increase with N
fractionParTV = 0.05;

Nrange = 50:50:100#100:150#800
Nsample = 5
degMin = 30
degMax = 30
function simulAndEstVarN(model::ErgmDirBin1,Nrange::StepRange{Int,Int},Nsample::Int,degMin::Int,degMax::Int)
    density=zeros(length(Nrange),3)
    dgpFits = Array{Array{Float64,1},2}(length(Nrange),3 )
    #generate multiple versions of the degree sequences for the dgp
     for n = 1:length(Nrange)
        println(n)
        Nn = Nrange[n]
        fits,degs =  dgpFitVarN(Nn,degMin,degMax;exponent = 1)
        density[n,1] = sum(degs[1:Nn])/(Nn^2-Nn)
        dgpFits[n,1] = fits
        fits,degs =  dgpFitVarN(Nn,degMin,Nn-degMin;exponent = 0.4 )
        density[n,2] = sum(degs[1:Nn])/(Nn^2-Nn)
        dgpFits[n,2] = fits
        fits,degs =  dgpFitVarN(Nn,degMin,Nn-degMin;exponent = 1 )
        density[n,3] = sum(degs[1:Nn])/(Nn^2-Nn)
        dgpFits[n,3] = fits
    end
    estFits = Array{Array{Float64,2},2}(length(Nrange),3 )
    @time for n = 1:length(Nrange)

        Nn = Nrange[n]
        for i=1:size(dgpFits)[2]
            # sample the dgpFits
            matSamp =  StaticNets.sampl(ErgmDirBin1(dgpFits[n,i]),Nsample;  parGroupsIO = dgpFits[n,i])
            degsSamp = [sumSq(matSamp,2);sumSq(matSamp,1)]
            #println(degsSamp   )
            tmpEst = zeros(2Nn,Nsample)
            for s=1:Nsample
                tmpEst[:,s],~ = StaticNets.estimate(fooErgmDirBin1,degIO = degsSamp[:,s])

            end
            estFits[n,i] = tmpEst
        end
            println(n)
    end
 return dgpFits,estFits
end

function simulAndEstVarN(model::ErgmDirBin1,dgpDegs::Array{<:Real,1},Nsample::Int)
    # if the degree sequence is given as input use that
    dgpFits,~,~ = meanFit ,~ = StaticNets.estimate(StaticNets.ErgmDirBin1(dgpDegs),degIO = dgpDegs  )
    N2 = length(dgpDegs)
    estFits = Array{Float64,2}(N2, Nsample )
    estFits = zeros(N2,Nsample)
    matSamp = StaticNets.sampl(ErgmDirBin1(dgpFits),Nsample;  parGroupsIO = dgpFits)
    degsSamp = [sumSq(matSamp,2);sumSq(matSamp,1)]
        for s=1:Nsample
            estFits[:,s],~ = StaticNets.estimate(fooErgmDirBin1,degIO = degsSamp[:,s])
            println(s)
        end
 return dgpFits,estFits
end

# sample and estimate networks of increasing size AND STORE DATA
#@time dgpFits,estFits = simulAndEstVarN(fooErgmDirBin1,Nrange,Nsample,degMin,degMax)

# Load data and plot results
save_fold = "./data/estimatesTest/asympTest/"

data = load(save_fold*"asympTestStaticLarge.jld")
estFits = data["estFits"]
dgpFits = data["dgpFits"]
N,Scal = size(dgpFits)
Sam = size(estFits[1,1])[2]
expDegsDgp = copy(dgpFits)
for n=1:N,scal=1:Scal expDegsDgp[n,scal] = expValStats(fooErgmDirBin1,dgpFits[n,scal]) end
expDegsEst = copy(estFits)
for n=1:N,scal=1:Scal,sam =1:Sam  expDegsEst[n,scal][:,sam] = expValStats(fooErgmDirBin1,estFits[n,scal][:,sam]) end
# plot the Degrees in the dgp
for n=1:N
    scal = 3
    plot(expDegsDgp[n,scal][1:midInd(expDegsDgp[n,scal])])
 end
 grid()
# plot densities
for scal=[1,3]
     dens = [sum(expDegsDgp[n,scal])/(length(expDegsDgp[n,scal])*(length(expDegsDgp[n,scal])-1)) for n=1:N]
     size = [length(expDegsDgp[n,scal]) for n=1:N]
     plot(size,dens ,".-")
  end
  grid()
  xlabel("number of nodes")
  ylabel("density")
#Plot mean relative bias over ensembles
for n=1:N
    scal = 3

 diffDegs = (expDegsDgp[n,scal] .- expDegsEst[n,scal])./expDegsDgp[n,scal]
  sampMenDiff = meanSq(diffDegs,2)
 plot(expDegsDgp[n,scal][1:midInd(expDegsDgp[n,scal])],n*0.2 + sampMenDiff[1:midInd(sampMenDiff)],".")
 end
 yticks(0.2:0.2:N*0.2)
 grid()
for n=1:N
     scal = 3
  diffDegs = (expDegsDgp[n,scal] .- expDegsEst[n,scal])./expDegsDgp[n,scal]
  sampStdDiff = squeeze(std(diffDegs,2),2)
  plot(expDegsDgp[n,scal][1:midInd(sampStdDiff)],n*0.2 + sampStdDiff[1:midInd(sampStdDiff)],".")
 end
  yticks(0.2:0.2:N*0.2)
  grid()
# load eMid data and run the estimates tests on realistically spaced fitnesses
LoadData =true
 halfPeriod=false
 fold_Path =  "/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/juliaFiles/"
 loadFilePartialName = "Weekly_eMid_Data_from_"
 halfPeriod? periodEndStr =  "2012_03_12_to_2015_02_27.jld": periodEndStr =  "2009_06_22_to_2015_02_27.jld"
 @load(fold_Path*loadFilePartialName*periodEndStr, AeMidWeekly_T,banksIDs,inactiveBanks, YeMidWeekly_T,weekInd,datesONeMid ,degsIO_T,strIO_T)
Sam = 500
degsIO_T = [sumSq(AeMidWeekly_T,2);sumSq(AeMidWeekly_T,1)]
#plot(sumSq(degsIO_T.>0,1)
meanDegsEmid = meanSq(degsIO_T,2)
dgpDegs =  round.(Int,meanDegsEmid)
dgpFits,estFits = simulAndEstVarN(fooErgmDirBin1,dgpDegs,Sam)
#plot empirical cumulative distribution
dgpDegDistr = ecdf(dgpDegs)
 plot(0:0.1:maximum(dgpDegs),dgpDegDistr(0:0.1:maximum(dgpDegs)),".-")
 grid()
 xticks(0:maximum(dgpDegs))
 xlabel("Degree")
 ylabel("emprical distribution")
expDegsDgp = expValStats(fooErgmDirBin1,dgpFits)
expDegsEst = copy(estFits)
for sam =1:Sam  expDegsEst[:,sam] = expValStats(fooErgmDirBin1,estFits[:,sam]) end
diff =  (expDegsDgp .- expDegsEst)
diffRel = diff./expDegsDgp
 diffRel[expDegsDgp.==0,:] = diff[expDegsDgp.==0,:]
# plot mean and std of relative error degrees estimates over the ensemble
fig, ax = subplots(2,1)
 ax[1,1][:plot](expDegsDgp,mean(diffRel,2),".")
 xticks(1:ceil(maximum(expDegsDgp)))
 ax[2,1][:plot](expDegsDgp,std(diffRel,2),".")
 ax[1,1][:grid]("on")
 ax[2,1][:grid]("on")
 ax[1,1][:set_xticks](0:ceil(maximum(expDegsDgp)))
 ax[1,1][:set_ylabel]("mean of relative error degs")
 ax[2,1][:set_xticks](0:ceil(maximum(expDegsDgp)))
 ax[2,1][:set_ylabel]("std of relative error degs")

indPlot = 7
plt[:hist](diffRel[indPlot,:],bins = 10,normed =true);
 grid()
 title("deg = $(round(Int,expDegsDgp[indPlot]))")
plt[:hist](diffRel[indmax(expDegsDgp),:],bins = 18,normed =true);
 grid()
 title("deg = $(round(Int,maximum(expDegsDgp)))")

 #

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
include("../examples/solenoid.jl")
using Zygote
using Test

s = [1.0,4.0,0.1]
x = rand(3)
g(x) = jacobian(x -> dstep(x,s), x)[1]
@show g(rand(3))

@time for i = 1:100
	x .= rand(3)
	du_ad = g(x)
end

@time for i = 1:100
	x .= rand(3)
	du_ana = dstep(x, s)
end
#@test du_ad ≈ du_ana atol=1.e-8

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
using BrownianDynamics
using OffsetArrays
using DifferentialEquations
using DiffEqDiffTools
using DSP
using FFTW
using Test
using SparseArrays
using LinearAlgebra

# write your own tests here
@testset "conv_mat" begin
    vv = [1.0;1.0+im;2.0;1.0*im;1.0;3.0+im]
    T = Complex{Float64}
    nu = length(vv)

    n = 2nu-1
    np2 = n > 1024 ? nextprod([2,3,5],n) : nextpow(2,n)
    upad = [vv;zeros(T,np2-nu)]

    p = plan_fft!(upad)
    ip = plan_ifft!(upad)

    @test conv_mat(vv) == conv_mat(p,ip,vv)
    @test isapprox(conv_mat(p,ip,vv)*vv,conv(vv,vv)[4:9])
    u = rand(6)
    v = rand(6)
    @test isapprox(conv_mat(u)*v,conv(u,v)[4:9])
end

# @testset "flux!" begin
#     # Test that the Jacobian is correct.
# end

## Initialize a spectral system.
nn = 64
nd2 = round(Int,nn/2)
p = OffsetArray{Complex{Float64}}(undef,-round(Int,nn/2):round(Int,nn/2-1))
p[:] = zeros(Complex{Float64},nn)
p[-2:2] = [0.0 0.0 1.0 0.0 0.0]
vv = OffsetArray{Complex{Float64}}(undef,-round(Int,nn/2):round(Int,nn/2-1))
vv[:] = zeros(Complex{Float64},nn)
# vv[-2:2] = [0.25 0.0 1.0 0.0 0.25]
vv[-2:2] = 2*[0.0 0.5im 0.0 -0.5im 0.0]
tt = OffsetArray{Complex{Float64}}(undef,-round(Int, nn/2):round(Int,nn/2-1))
tt[:] = zeros(Complex{Float64},nn)
tt[-1:1] = [0.0 1.0 0.0]

p = collect(p)
vv = collect(vv)
tt = collect(tt)
ϕ0 = [p;tt]

T = Complex{Float64}
nu = length(vv)
nv = length(vv)
n = nu+nv-1
np2 = n > 1024 ? nextprod([2,3,5],n) : nextpow(2,n)
upad = [vv;zeros(T,np2-nu)]
vpad = [vv;zeros(T,np2-nv)]

plan = plan_fft!(upad)
iplan = plan_ifft!(upad)
y = similar(upad)

A = 4.0
B = [1.0 100.0]
force = 4.0
T0 = 1.0

params = SpectralParameters{Float64}(plan,iplan,A,B |> vec,force,vv,T0)
tspan = (0.0,0.5)
prob = ODEProblem(spectral_rhs!,ϕ0,tspan,params)

@testset "spectral_rhs!" begin
    jac_diff = DiffEqDiffTools.finite_difference_jacobian((dϕ,ϕ) -> spectral_rhs!(dϕ,ϕ,params,0),ϕ0)
    jac = spzeros(Complex{Float64},2nn,2nn)
    jac = spectral_rhs!(Val{:jac},jac,ϕ0,params,0)
    @test norm(jac*ϕ0 - jac_diff*ϕ0)/nn < 1e-9

    sol = solve(prob,TRBDF2(autodiff=false),abstol=1e-10,reltol=1e-10)

    jac_diff = DiffEqDiffTools.finite_difference_jacobian((dϕ,ϕ) -> spectral_rhs!(dϕ,ϕ,params,0),sol[end])
    jac = spzeros(Complex{Float64},2nn,2nn)
    jac = spectral_rhs!(Val{:jac},jac,sol[end],params,0)
    @test norm(jac*sol[end] - jac_diff*sol[end])/nn < 1e-9
end

@testset "solve_steady_state" begin
    integrator = init(prob,TRBDF2(autodiff=false))
    u = solve_steady_state(integrator,params,steadytol=1e-10,maxiters=5)
    @test norm([spectral_rhs!(copy(ϕ0),u,params,0);u[nd2+1]-1])/nn < 1e-10
end


text:
function get_invest_values(u::String="*")
    invest_then = 0
    invest_now  = 0
    invest_paths = glob(joinpath(Pinguu.base_dir[2:end], "stats/users/" * u * "/invest.json"))
    for path in invest_paths
        tmp = JSON.parsefile(path)
        for (symbol, stat) in tmp
            invest_then += stat["value"]
            invest_now  += stat["value_now"]
        end
    end
    return invest_then, invest_now
end
return get_invest_values


text:
test_string = "ecl:gry pid:860033327 eyr:2020 hcl:#fffffd
byr:1937 iyr:2017 cid:147 hgt:183cm

iyr:2013 ecl:amb cid:350 eyr:2023 pid:028048884
hcl:#cfa07d byr:1929

hcl:#ae17e1 iyr:2013
eyr:2024
ecl:brn pid:760753108 byr:1931
hgt:179cm

hcl:#cfa07d eyr:2025 pid:166559648
iyr:2011 ecl:brn hgt:59in"

test_list = split(test_string, '\n')

list = readlines("day4_input.txt")

const PASSWDSET_NOCID = Set{String}(["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"])

#=
byr (Birth Year) - four digits; at least 1920 and at most 2002.
iyr (Issue Year) - four digits; at least 2010 and at most 2020.
eyr (Expiration Year) - four digits; at least 2020 and at most 2030.
hgt (Height) - a number followed by either cm or in:
If cm, the number must be at least 150 and at most 193.
If in, the number must be at least 59 and at most 76.
hcl (Hair Color) - a # followed by exactly six characters 0-9 or a-f.
ecl (Eye Color) - exactly one of: amb blu brn gry grn hzl oth.
pid (Passport ID) - a nine-digit number, including leading zeroes.
 =#

# answer 1: 196
# answer 2: ??
validator = Dict(
    "byr" => x -> length(x) == 4 && prod(['0' <= c <= '9' for c in x]) && 1920 ≤ parse(Int, x) ≤ 2002,
    "iyr" => x -> length(x) == 4 && prod(['0' <= c <= '9' for c in x]) && 2010 ≤ parse(Int, x) ≤ 2020,
    "eyr" => x -> length(x) == 4 && prod(['0' <= c <= '9' for c in x]) && 2020 ≤ parse(Int, x) ≤ 2030,
    "hgt" => x -> length(x) > 2 && x[end-1:end] in Set(["cm", "in"]) && ((x[end-1:end]== "cm" && prod(['0' <= c <= '9' for c in x[1:end-2]]) && 150 ≤ parse(Int, x[1:end-2]) ≤ 193 || (x[end-1:end] == "in" && prod(['0' <= c <= '9' for c in x[1:end-2]]) && 59 ≤ parse(Int, x[1:end-2]) ≤ 76))),
    "hcl" => x -> length(x) == 7 && x[1] == '#' && prod(['0' <= c <= 'f' for c in x[2:end]]),
    "ecl" => x -> x in Set(["amb", "blu", "brn", "gry", "grn", "hzl", "oth"]),
    "pid" => x -> length(x) == 9 && prod(['0' <= c <= '9' for c in x])
)

function validate_passinfo(list; part = 1)
    pass_keys = Set{String}(["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"])
    valid = 0
    nmax = length(list)
    info = ""
    for n in 1:nmax+1
        if n == nmax + 1 || list[n] == ""
            info_fields = split(info, [':', ' '])
            info_keys = Set(info_fields[1:2:end])
            if pass_keys ⊆ info_keys
                if part == 2
                    valid_keys = true
                    for j in 1:2:length(info_fields)
                        if info_fields[j] in pass_keys
                            result = validator[info_fields[j]](info_fields[j+1])
                            valid_keys *= validator[info_fields[j]](info_fields[j+1])
                        end
                    end
                    if valid_keys
                        valid += 1
                    end
                else
                    valid += 1
                end
            end
            info = ""
        else
            info *= list[n] * " "
        end
    end
    return valid
end

@show validate_passinfo(list, part=1)
@show validate_passinfo(list, part=2)

# A simpler validator after AaronRMatthis suggestions:

validator2 = Dict(
    "byr" => x -> length(x) == 4
        && all(isdigit, x) && "1920" ≤ x ≤ "2002",
    "iyr" => x -> length(x) == 4
        && all(isdigit, x) && "2010" ≤ x ≤ "2020",
    "eyr" => x -> length(x) == 4
        && all(isdigit, x) && "2020" ≤ x ≤ "2030",
    "hgt" => x -> length(x) > 2
        && ((x[end-1:end]== "cm"
            && all(isdigit, x[1:end-2]) && 150 ≤ parse(Int, x[1:end-2]) ≤ 193
            || (x[end-1:end] == "in" && all(isdigit, x[1:end-2]) && 59 ≤ parse(Int, x[1:end-2]) ≤ 76))),
    "hcl" => x -> length(x) == 7
        && x[1] == '#' && all(c->'0'<=c<='f',x[2:end]),
    "ecl" => x -> x in Set(["amb", "blu", "brn", "gry", "grn", "hzl", "oth"]),
    "pid" => x -> length(x) == 9 && all(isdigit, x)
)

nothing

# AaronRMatthis
#=
const dicts = map(split(read("day4_input.txt", String),"\n"=>" "), "  ") do single
    cols = split(single)
    d = Dict{String,String}()
    foreach(((key,val),)->(d[key] = val),split.(cols,':'))
    return d
end
 =#


text:
export Instrument

abstract type Instrument end

function Base.getproperty(instrument::Instrument, sym::Symbol)
    if sym === :tuning && !hasfield(typeof(instrument), sym)
        return EqualTemperedTuning(12)
    else # fallback to getfield
        return getfield(instrument, sym)
    end
end


text:
using SymEngine_jll: libsymengine

function __init__()
    __init_constants()
    __init_type_ids()
end

###########
## TYPES ##
###########

"""
    Expression <: Number

A symbolic expression.

```julia-repl
julia> expr = (Variable(:x) + 1)^2
(1 + x)^2

julia> Expression(2)
2

julia> Expression(Variable(:x))
x
```
"""
mutable struct Expression <: Number
    ptr::Ptr{Cvoid}

    function Expression()
        z = new(C_NULL)
        ccall((:basic_new_stack, libsymengine), Nothing, (Ref{Expression},), z)
        finalizer(free!, z)
        return z
    end

    function Expression(v::Ptr{Cvoid})
        z = new(v)
        return z
    end
end

Expression(ex::Expression) = ex
Expression(T) = convert(Expression, T)

free!(b::Expression) =
    ccall((:basic_free_stack, libsymengine), Nothing, (Ref{Expression},), b)

"""
    ExpressionRef

Basically an `Expression` with the difference that it internally only contains
a reference to an `Expression`. If the underlying reference is freed there is
no guarantee that `ExpressionRef` points to valid memory.
"""
struct ExpressionRef
    ptr::Ptr{Cvoid}

    ExpressionRef(ptr::Ptr{Cvoid}) = new(ptr)
end
ExpressionRef() = ExpressionRef(Ptr{Cvoid}())
ExpressionRef(ex::Expression) = ExpressionRef(ex.ptr)
ExpressionRef(x) = convert(ExpressionRef, x)

"""
    Variable(name::Union{String,Symbol}, indices...) <: Number

A data structure representing a variable.

```julia-repl
julia> Variable(:a)
a

julia> Variable(:x, 1)
x₁

julia> Variable(:x, 10, 5)
x₁₀₋₅
```

### Equality and ordering

Variables are identified by their `name` and `indices`.
That is, two variables are equal if and only if they have the same `name` and `indices`.

```julia-repl
julia> Variable(:a) == Variable(:a)
true

julia> Variable(:a, 1) == Variable(:a, 2)
false
```

Similarly, variables are first ordered lexicographically by their name and then by their indices.
```julia-repl
julia> Variable(:a, 1) < Variable(:a, 2)
true

julia> Variable(:a, 1) < Variable(:b, 1)
true

julia> a = [Variable(:a, i, j) for i in 1:2, j in 1:2]
2×2 Array{Variable,2}:
 a₁₋₁  a₁₋₂
 a₂₋₁  a₂₋₂

julia> sort(vec(a))
4-element Array{Variable,1}:
 a₁₋₁
 a₂₋₁
 a₁₋₂
 a₂₋₂
```
"""
struct Variable <: Number
    ex::Expression

    function Variable(ex::Expression)
        @assert class(ex) == :Symbol
        new(ex)
    end
    Variable(s::Union{String,Symbol}) = new(Expression(s))
end

name(v::Variable) = Symbol(to_string(v))

const Basic = Union{Expression,ExpressionRef,Variable}

Base.convert(::Type{Expression}, v::Variable) = v.ex
Base.convert(::Type{Expression}, e::ExpressionRef) = copy(e)
Base.convert(::Type{ExpressionRef}, e::Expression) = ExpressionRef(e)
Base.convert(::Type{ExpressionRef}, v::Variable) = ExpressionRef(v.ex)

Base.promote_rule(::Type{<:Basic}, ::Type{<:Number}) = Expression
Base.promote_rule(::Type{Expression}, ::Type{Variable}) = Expression
Base.promote_rule(::Type{Expression}, ::Type{ExpressionRef}) = Expression
Base.promote_rule(::Type{Variable}, ::Type{ExpressionRef}) = Expression

function to_string(x::Basic)
    b = ExpressionRef(x)
    if b.ptr == C_NULL
        return ""
    end
    a = ccall((:basic_str_julia, libsymengine), Cstring, (Ref{ExpressionRef},), b)
    string = unsafe_string(a)
    ccall((:basic_str_free, libsymengine), Nothing, (Cstring,), a)
    return string
end

Base.show(io::IO, b::Basic) = print(io, to_string(b))
function Base.show(io::IO, mime::MIME"text/latex", b::Basic)
    print(io, "\$\$ ")
    print(io, to_string(b))
    print(io, " \$\$")
end

function Base.hash(ex::Basic, h::UInt)
    h2 = ccall((:basic_hash, libsymengine), UInt, (Ref{ExpressionRef},), ex)
    Base.hash_uint(3h - h2)
end

function _copy(x::Basic)
    y = Expression()
    ccall(
        (:basic_assign, libsymengine),
        Nothing,
        (Ref{Expression}, Ref{ExpressionRef}),
        y,
        x,
    )
    y
end
Base.copy(x::Basic) = _copy(x)
Base.copy(x::Variable) = Variable(_copy(x))

function Base.:(==)(b1::Basic, b2::Basic)
    ccall(
        (:basic_eq, libsymengine),
        Int,
        (Ref{ExpressionRef}, Ref{ExpressionRef}),
        b1,
        b2,
    ) == 1
end

const EXPR_1 = Expression(C_NULL)
const EXPR_0 = Expression(C_NULL)
const EXPR_M1 = Expression(C_NULL)

Base.zero(::Basic) = zero(Expression)
Base.zero(::Type{<:Basic}) = Expression(0)
Base.iszero(x::Basic) = x == EXPR_0
Base.one(::Basic) = one(Expression)
Base.one(::Type{<:Basic}) = Expression(1)
Base.isone(x::Basic) = x == EXPR_1

is_one(x::Union{Number,Basic}) = isone(x)
is_one(x) = false
is_zero(x::Union{Number,Basic}) = iszero(x)
is_zero(x) = false
is_minus_one(x::Basic) = x == EXPR_M1
is_minus_one(x) = x == -1

for op in [:im, :π, :ℯ, :γ, :catalan]
    @eval begin
        const $(Symbol("__", op)) = Expression(C_NULL)
    end
end
const SYMENGINE_CONSTANTS = Dict{Expression,Irrational}()

macro init_constant(op, libnm)
    tup = (Base.Symbol("basic_const_$libnm"), libsymengine)
    alloc_tup = (:basic_new_stack, libsymengine)
    op_name = Symbol("__", op)
    c = quote
        ccall($alloc_tup, Nothing, (Ref{Expression},), $op_name)
        ccall($tup, Nothing, (Ref{Expression},), $op_name)
        finalizer(free!, $op_name)
        $(
            op != :im ? :(SYMENGINE_CONSTANTS[$op_name] = Base.MathConstants.$op) :
            :(nothing)
        )
    end
    c
end

function __init_constants()
    @init_constant im I
    @init_constant π pi
    @init_constant ℯ E
    @init_constant γ EulerGamma
    @init_constant catalan Catalan

    ccall((:basic_new_stack, libsymengine), Nothing, (Ref{Expression},), EXPR_1)
    ccall((:integer_set_si, libsymengine), Nothing, (Ref{Expression}, Int), EXPR_1, 1)
    finalizer(free!, EXPR_1)

    ccall((:basic_new_stack, libsymengine), Nothing, (Ref{Expression},), EXPR_0)
    ccall((:integer_set_si, libsymengine), Nothing, (Ref{Expression}, Int), EXPR_0, 0)
    finalizer(free!, EXPR_0)

    ccall((:basic_new_stack, libsymengine), Nothing, (Ref{Expression},), EXPR_M1)
    ccall((:integer_set_si, libsymengine), Nothing, (Ref{Expression}, Int), EXPR_M1, -1)
    finalizer(free!, EXPR_M1)
end

################
## arithmetic ##
################

## main ops
for (op, inplace, libnm) in [(:+, :add!, :add), (:-, :sub!, :sub), (:*, :mul!, :mul)]
    # $(Expr(:., :Base, QuoteNode(op)))
    @eval begin
        function $(inplace)(a::Expression, b1::Basic, b2::Basic)
            ccall(
                $((Symbol("basic_", libnm), libsymengine)),
                Nothing,
                (Ref{Expression}, Ref{ExpressionRef}, Ref{ExpressionRef}),
                a,
                b1,
                b2,
            )
            return a
        end
        function Base.$op(b1::Basic, b2::Basic)
            $inplace(Expression(), b1, b2)
        end
    end
end

function div!(a::Expression, b1::Basic, b2::Basic)
    # class(b2) ∈ NUMBER_TYPES || throw(ArgumentError("Cannot divide by " * string(b2)))
    ccall(
        (:basic_div, libsymengine),
        Nothing,
        (Ref{Expression}, Ref{ExpressionRef}, Ref{ExpressionRef}),
        a,
        b1,
        b2,
    )
    return a
end
Base.:(/)(b1::Basic, b2::Basic) = div!(Expression(), b1, b2)

Base.:(//)(x::Basic, y::Basic) = x / y
Base.:(//)(x::Basic, y::Integer) = x / Expression(y)

function Base.:(^)(x::Basic, k::Integer)
    a = Expression()
    ccall(
        (:basic_pow, libsymengine),
        Nothing,
        (Ref{Expression}, Ref{ExpressionRef}, Ref{Expression}),
        a,
        x,
        Expression(k),
    )
    return a
end
Base.:+(b::Basic) = b
Base.:-(b::Basic) = Expression(0) - b

# Functions
macro make_func(arg1, arg2)
    quote
        function $(esc(arg1))(b::Basic)
            a = Expression()
            ccall(
                ($(QuoteNode(arg2)), libsymengine),
                Base.Nothing,
                (Base.Ref{Expression}, Base.Ref{ExpressionRef}),
                a,
                b,
            )
            return a
        end
    end
end

@make_func symengine_expand basic_expand

function Base.sin(b::Basic)
    a = Expression()
    ccall((:basic_sin, libsymengine), Nothing, (Ref{Expression}, Ref{ExpressionRef}), a, b)
    return a
end
function Base.cos(b::Basic)
    a = Expression()
    ccall((:basic_cos, libsymengine), Nothing, (Ref{Expression}, Ref{ExpressionRef}), a, b)
    return a
end
Base.sincos(x::Basic) = (sin(x), cos(x))

function Base.log(b::ModelKit.Basic)
    a = Expression()
    ccall((:basic_log, libsymengine), Nothing, (Ref{Expression}, Ref{ExpressionRef}), a, b)
    return a
end
##############################
## conversion to Expression ##
##############################

function Base.convert(::Type{Expression}, s::String)
    a = Expression()
    ccall((:symbol_set, libsymengine), Nothing, (Ref{Expression}, Ptr{Int8}), a, s)
    return a
end
Base.convert(::Type{Expression}, s::Symbol) = convert(Expression, string(s))

Base.convert(::Type{Expression}, x::Irrational{:γ}) = __γ
Base.convert(::Type{Expression}, x::Irrational{:π}) = __π
Base.convert(::Type{Expression}, x::Irrational{:ℯ}) = __ℯ
Base.convert(::Type{Expression}, x::Irrational{:catalan}) = __catalan

@static if Int == Int64 && Sys.iswindows()
    # basic number types
    for (f, T) in
        [(:integer_set_si, Int32), (:integer_set_ui, UInt32), (:real_double_set_d, Float64)]
        @eval begin
            function Base.convert(::Type{Expression}, x::$T)
                a = Expression()
                ccall(($(QuoteNode(f)), libsymengine), Nothing, (Ref{Expression}, $T), a, x)
                return a
            end
        end
    end

    Base.convert(::Type{Expression}, x::Int64) = convert(Expression, BigInt(x))
    Base.convert(::Type{Expression}, x::UInt64) = convert(Expression, BigInt(x))
else
    # basic number types
    for (f, T) in
        [(:integer_set_si, Int), (:integer_set_ui, UInt), (:real_double_set_d, Float64)]
        @eval begin
            function Base.convert(::Type{Expression}, x::$T)
                a = Expression()
                ccall(($(QuoteNode(f)), libsymengine), Nothing, (Ref{Expression}, $T), a, x)
                return a
            end
        end
    end
    Base.convert(::Type{Expression}, x::Int32) = convert(Expression, convert(Int, x))
    Base.convert(::Type{Expression}, x::UInt32) = convert(Expression, convert(UInt, x))
end

function Base.convert(::Type{Expression}, x::BigInt)
    a = Expression()
    ccall((:integer_set_mpz, libsymengine), Nothing, (Ref{Expression}, Ref{BigInt}), a, x)
    return a
end

Base.convert(::Type{Expression}, x::Union{Float16,Float32}) =
    convert(Expression, convert(Float64, x))
Base.convert(::Type{Expression}, x::AbstractFloat) =
    convert(Expression, convert(BigFloat, x))
function Base.convert(::Type{Expression}, x::BigFloat)
    if (x.prec <= 53)
        return convert(Expression, Float64(x))
    else
        a = Expression()
        ccall(
            (:real_mpfr_set, libsymengine),
            Nothing,
            (Ref{Expression}, Ref{BigFloat}),
            a,
            x,
        )
        return a
    end
end

Base.convert(::Type{Expression}, x::Integer) = Expression(BigInt(x))
Base.convert(::Type{Expression}, x::Rational) =
    Expression(numerator(x)) / Expression(denominator(x))
Base.convert(::Type{Expression}, x::Complex) =
    Expression(real(x)) + Expression(imag(x)) * __im

Base.complex(x::Expression, y::Expression) = x + y * __im
Base.complex(x::Expression, y::Real) = x + y * __im
Base.complex(x::Real, y::Expression) = x + y * __im

################################
## Iterating over expressions ##
################################

mutable struct ExpressionSet
    ptr::Ptr{Cvoid}

    function ExpressionSet()
        z = new(ccall((:setbasic_new, libsymengine), Ptr{Cvoid}, ()))
        finalizer(free!, z)
        z
    end
end

function free!(x::ExpressionSet)
    if x.ptr != C_NULL
        ccall((:setbasic_free, libsymengine), Nothing, (Ptr{Cvoid},), x.ptr)
        x.ptr = C_NULL
    end
    nothing
end

Base.length(s::ExpressionSet) =
    ccall((:setbasic_size, libsymengine), Int, (Ptr{Cvoid},), s.ptr)

function Base.getindex(s::ExpressionSet, n::Int)
    result = Expression()
    ccall(
        (:setbasic_get, libsymengine),
        Nothing,
        (Ptr{Cvoid}, Int, Ref{Expression}),
        s.ptr,
        n - 1,
        result,
    )
    result
end

_variables(ex::Variable) = [ex]
function _variables(ex::Expression)
    syms = ExpressionSet()
    ccall(
        (:basic_free_symbols, libsymengine),
        Nothing,
        (Ref{ExpressionRef}, Ptr{Cvoid}),
        ex,
        syms.ptr,
    )
    sort!([Variable(syms[i]) for i = 1:length(syms)])
end

function differentiate(f::Basic, v::Variable)
    a = Expression()
    ret = ccall(
        (:basic_diff, libsymengine),
        Int,
        (Ref{Basic}, Ref{ExpressionRef}, Ref{ExpressionRef}),
        a,
        f,
        v,
    )
    return a
end
function differentiate(f::Basic, v::Variable, n::Int)
    n < 0 && throw(DomainError("n must be non-negative integer"))
    n == 0 && return f
    n == 1 && return differentiate(f, v)
    n > 1 && return differentiate(differentiate(f, v), v, n - 1)
end

# Get class of an Expression

const REAL_NUMBER_TYPES = [:Integer, :RealDouble, :Rational, :RealMPFR]
const COMPLEX_NUMBER_TYPES = [:Complex, :ComplexDouble, :ComplexMPC]
const NUMBER_TYPES = [
    REAL_NUMBER_TYPES
    COMPLEX_NUMBER_TYPES
]

function type_id(s::Basic)
    ccall((:basic_get_type, libsymengine), UInt, (Ref{ExpressionRef},), s)
end

function get_class_from_type_id(id::UInt)
    a = ccall((:basic_get_class_from_id, libsymengine), Ptr{UInt8}, (Int,), id)
    str = unsafe_string(a)
    ccall((:basic_str_free, libsymengine), Nothing, (Ptr{UInt8},), a)
    Symbol(str)
end

# prepopulate the dict
const TYPE_IDS = Dict{UInt,Symbol}()

function __init_type_ids()
    x = Expression("x")
    values = [
        x,
        Expression(1),
        Expression(0.5),
        Expression(2 // 3),
        Expression(big(0.5)),
        Expression(5 + 3im),
        Expression(0.5 + 0.2im),
        Expression(0.5 + big(0.2) * im),
        2x,
        x + 2,
        x^2,
    ]
    for v in values
        class(v)
    end
    nothing
end

function class(e::Basic)
    id = type_id(e)
    if haskey(TYPE_IDS, id)
        TYPE_IDS[id]
    else
        # add for futurue fast lookup
        cls = get_class_from_type_id(id)
        TYPE_IDS[id] = cls
        cls
    end
end

mutable struct ExprVec <: AbstractVector{ExpressionRef}
    ptr::Ptr{Cvoid}
    m::Union{Nothing,Ptr{ModelKit.ExpressionRef}}

    function ExprVec()
        ptr = ccall((:vecbasic_new, libsymengine), Ptr{Cvoid}, ())
        z = new(ptr, nothing)
        finalizer(free!, z)
        z
    end
end
function vec_set_ptr!(v::ExprVec)
    v.m = unsafe_load(Ptr{Ptr{ModelKit.ExpressionRef}}(v.ptr))
    v
end

function free!(x::ExprVec)
    if x.ptr != C_NULL
        ccall((:vecbasic_free, libsymengine), Nothing, (Ptr{Cvoid},), x.ptr)
        x.ptr = C_NULL
    end
end

Base.eltype(v::ExprVec) = ExpressionRef
Base.length(s::ExprVec) = ccall((:vecbasic_size, libsymengine), UInt, (Ptr{Cvoid},), s.ptr)
Base.size(s::ExprVec) = (length(s),)

function Base.getindex(v::ExprVec, n)
    @boundscheck checkbounds(v, n)
    if v.m === nothing
        vec_set_ptr!(v)
        unsafe_load(v.m, n)
    else
        unsafe_load(v.m, n)
    end
end

function Base.push!(v::ExprVec, x::Basic)
    ccall(
        (:vecbasic_push_back, libsymengine),
        Nothing,
        (Ptr{Cvoid}, Ref{ExpressionRef}),
        v.ptr,
        x,
    )
    v
end

args(ex::Basic) = args!(ExprVec(), ex)
function args!(vec::ExprVec, ex::Basic)
    ccall(
        (:basic_get_args, libsymengine),
        Nothing,
        (Ref{ExpressionRef}, Ptr{Cvoid}),
        ex,
        vec.ptr,
    )
    vec_set_ptr!(vec)
end

################################
## conversion from Expression ##
################################

@static if Int == Int64 && Sys.iswindows()
    function Base.convert(::Type{Int32}, n::Basic)
        @assert class(n) == :Integer
        ccall((:integer_get_si, libsymengine), Int32, (Ref{ExpressionRef},), n)
    end
    Base.convert(::Type{Int64}, n::Basic) = convert(Int64, convert(BigInt, n))
else
    function Base.convert(::Type{Int}, n::Basic)
        @assert class(n) == :Integer
        ccall((:integer_get_si, libsymengine), Int, (Ref{ExpressionRef},), n)
    end
    Base.convert(::Type{Int32}, n::Basic) = convert(Int32, convert(Int, n))
end

function Base.convert(::Type{BigInt}, c::Basic)
    a = BigInt()
    ccall(
        (:integer_get_mpz, libsymengine),
        Nothing,
        (Ref{BigInt}, Ref{ExpressionRef}),
        a,
        c,
    )
    return a
end

function Base.convert(::Type{BigFloat}, c::Basic)
    a = BigFloat()
    ccall(
        (:real_mpfr_get, libsymengine),
        Nothing,
        (Ref{BigFloat}, Ref{ExpressionRef}),
        a,
        c,
    )
    return a
end

function Base.convert(::Type{Float64}, c::Basic)
    return ccall((:real_double_get_d, libsymengine), Cdouble, (Ref{ExpressionRef},), c)
end

function Base.real(x::Basic)
    if class(x) ∈ COMPLEX_NUMBER_TYPES
        a = Expression()
        ccall(
            (:complex_base_real_part, libsymengine),
            Nothing,
            (Ref{Expression}, Ref{ExpressionRef}),
            a,
            x,
        )
        return a
    else
        return x
    end
end
function Base.imag(x::Basic)
    if class(x) ∈ COMPLEX_NUMBER_TYPES
        a = Expression()
        ccall(
            (:complex_base_imaginary_part, libsymengine),
            Nothing,
            (Ref{Expression}, Ref{ExpressionRef}),
            a,
            x,
        )
        return a
    else
        return Expression(0)
    end
end

function _numer_denom(x::Basic)
    a, b = Expression(), Expression()
    ccall(
        (:basic_as_numer_denom, libsymengine),
        Nothing,
        (Ref{Expression}, Ref{Expression}, Ref{ExpressionRef}),
        a,
        b,
        x,
    )
    return a, b
end

is_number(ex::Expression) = class(ex) in NUMBER_TYPES

"""
    to_number(x::Expression)

Tries to unpack the `Expression` `x` to a native number type.

```julia-repl
julia> x = to_number(Expression(2))
2

julia> typeof(x)
Int64
"""
function to_number(x::Basic)
    cls = class(x)

    if cls == :Integer
        n = convert(BigInt, x)
        if typemin(Int32) ≤ n ≤ typemax(Int32)
            return convert(Int32, n)
        elseif typemin(Int64) ≤ n ≤ typemax(Int64)
            return convert(Int64, n)
        elseif typemin(Int128) ≤ n ≤ typemax(Int128)
            return convert(Int128, n)
        else
            return n
        end
    elseif cls == :RealDouble
        return convert(Float64, x)
    elseif cls == :Rational
        a, b = _numer_denom(x)
        return to_number(a) // to_number(b)
    elseif cls == :RealMPFR
        return convert(BigFloat, x)
    elseif cls == :Constant
        return SYMENGINE_CONSTANTS[x]
    elseif cls in COMPLEX_NUMBER_TYPES
        a, b = reim(x)
        return complex(to_number(a), to_number(b))
    else
        return x
    end
end

function (::Type{T})(x::Basic) where {T<:Union{AbstractFloat,Rational,Complex,Integer}}
    y = to_number(x)
    if y isa Basic
        # Throw error to avoid endless loop
        # The T.name.name part is done the same way in Base
        throw(InexactError(T.name.name, T, x))
    else
        T(y)
    end
end

##########
## SUBS ##
##########

mutable struct ExpressionMap
    ptr::Ptr{Cvoid}

    function ExpressionMap()
        x = new(ccall((:mapbasicbasic_new, libsymengine), Ptr{Cvoid}, ()))
        finalizer(free!, x)
        x
    end
end

function free!(x::ExpressionMap)
    if x.ptr != C_NULL
        ccall((:mapbasicbasic_free, libsymengine), Nothing, (Ptr{Cvoid},), x.ptr)
        x.ptr = C_NULL
    end
end

ExpressionMap(args...) = ExpressionMap(ExpressionMap(), args...)
ExpressionMap(D::ExpressionMap) = D
function ExpressionMap(dict::Dict)
    c = ExpressionMap()
    for (key, value) in dict
        c[Expression(key)] = Expression(value)
    end
    return c
end
function ExpressionMap(
    D::ExpressionMap,
    (xs, ys)::Pair{<:AbstractArray{<:Basic},<:AbstractArray},
    args...,
)
    size(xs) == size(ys) ||
        throw(ArgumentError("Substitution arguments don't have the same size."))
    for (x, y) in zip(xs, ys)
        D[x] = y
    end
    ExpressionMap(D, args...)
end
function ExpressionMap(D::ExpressionMap, (x, y), args...)
    D[Expression(x)] = Expression(y)
    ExpressionMap(D, args...)
end

function Base.length(s::ExpressionMap)
    ccall((:mapbasicbasic_size, libsymengine), Int, (Ptr{Cvoid},), s.ptr)
end

Base.setindex!(s::ExpressionMap, v::Number, k::Basic) = setindex!(s, Expression(v), k)
function Base.setindex!(s::ExpressionMap, v::Basic, k::Basic)
    ccall(
        (:mapbasicbasic_insert, libsymengine),
        Nothing,
        (Ptr{Cvoid}, Ref{ExpressionRef}, Ref{ExpressionRef}),
        s.ptr,
        k,
        v,
    )
    v
end

function subs(ex::Basic, d::ExpressionMap)
    s = Expression()
    ccall(
        (:basic_subs, libsymengine),
        Nothing,
        (Ref{Expression}, Ref{ExpressionRef}, Ptr{Cvoid}),
        s,
        ex,
        d.ptr,
    )
    return s
end

function cse(exprs::Vector{Expression})
    vec = ExprVec()
    for ex in exprs
        push!(vec, ex)
    end

    replacement_syms = ExprVec()
    replacement_exprs = ExprVec()
    reduced_exprs = ExprVec()

    ccall(
        (:basic_cse, libsymengine),
        Nothing,
        (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        replacement_syms.ptr,
        replacement_exprs.ptr,
        reduced_exprs.ptr,
        vec.ptr,
    )

    # The symengine cse alg also considers in (x - a) * (y - a) the -a as a cse.
    # We don't want these since we can handle them with no overhead directly,
    # so we remove these again.
    bad_subs = ExpressionMap()
    substs = Dict{Expression,Expression}()
    for (sᵢ, exᵢ) in zip(replacement_syms, replacement_exprs)
        if class(exᵢ) == :Mul
            exs = args(exᵢ)
            if length(exs) == 2 && (is_minus_one(exs[1]) || is_one(exs[1]))
                bad_subs[sᵢ] = exᵢ
            else
                substs[sᵢ] = subs(exᵢ, bad_subs)
            end
        else
            substs[sᵢ] = subs(exᵢ, bad_subs)
        end
    end

    map(e -> subs(e, bad_subs), reduced_exprs), substs
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

mutable struct IterationChannel{D} <: AbstractChannel{D}
    data::Dict{Int,D}
    cond_take::Condition
    IterationChannel(data::Dict{Int,D}) where D = new{D}(data, Condition())
end

function put!(channel::IterationChannel, t, x)
    channel.data[t] = copy(x)
    notify(channel.cond_take)
    return channel
end

function take!(channel::IterationChannel, t)
    x = fetch(channel, t)
    delete!(channel.data, t)
    return x
end

isready(channel::IterationChannel) = length(channel.data) > 1
isready(channel::IterationChannel, t) = haskey(channel.data, t)

function fetch(channel::IterationChannel, t)
    wait(channel, t)
    return channel.data[t]
end

function wait(channel::IterationChannel, t)
    while !isready(channel, t)
        wait(channel.cond_take)
    end
end

mutable struct RunningAverage{D}
    average::D
    data::Vector{D}
    buffer::Dict{Int,D}

    RunningAverage(average::D, data::Vector{D}) where D = new{D}(average, data, Dict{Int,D}())
end

function average(running_average::RunningAverage)
    return running_average.average
end

function subtract!(running_average::RunningAverage, i::Integer)
    running_average.buffer[i] = copy(running_average.data[i])
end

function add!(running_average::RunningAverage, i::Integer, π::AbstractFloat)
    running_average.average -= π * running_average.buffer[i]
    running_average.average += π * running_average.data[i]
    return running_average
end

function add!(running_average::RunningAverage{D}, i::Integer, x::D, π::AbstractFloat) where D
    running_average.average -= π * running_average.buffer[i]
    running_average.average += π * x
    running_average.data[i] = copy(x)
    return running_average
end

IteratedValue{T <: AbstractFloat} = RemoteChannel{IterationChannel{T}}
RemoteRunningAverage{D} = RemoteChannel{Channel{RunningAverage{D}}}
RemoteIterates{A <: AbstractArray} = RemoteChannel{IterationChannel{A}}


text:
# fc = fcut($desc, $rθ)

using NBodyIPs, StaticArrays, BenchmarkTools, JuLIP, Test, Profile,
      LinearAlgebra
using NBodyIPs.Polys: NBPoly, StNBPoly
using NBodyIPs: evaluate

TRANSFORM = "r -> (2.9/r)^3"
rcut = 7.5
desc = BondLengthDesc(TRANSFORM, (:cos, (0.66*rcut), rcut))
r = 9.5 * (@SVector rand(6))
x = (2.9./r).^3
@btime NBodyIPs.fcut($desc, $r)
@btime NBodyIPs.invariants($desc, $x)

@which NBodyIPs.fcut(desc, r)

@btime NBodyIPs.fcut($(desc.cutoff), $x)
@which NBodyIPs.fcut(desc.cutoff, x)
@btime NBodyIPs.fcut($(desc.cutoff), 4.567)
@which NBodyIPs.fcut(desc.cutoff, 4.567)

fff = desc.cutoff.f
@btime $fff(4.567)

f1, df1, rc = let rc1=3.456, rc2=7.5
   r -> JuLIP.Potentials.coscut(r, rc1, rc2), r -> JuLIP.Potentials.coscut_d(r, rc1, rc2), rc2
end
@btime $f1(4.567)

f1w = JuLIP.Potentials.F64fun(f1)
@btime $f1w(4.567)
@btime $(fff.obj.x)(4.567)
fffw = JuLIP.Potentials.F64fun(fff.obj.x)
@btime $fffw(4.567)

cutoff1 = NBodyIPs.Cutoff(desc.cutoff.sym,
                 desc.cutoff.params,
                 JuLIP.Potentials.F64fun(f1),
                 JuLIP.Potentials.F64fun(df1),
                 cutoff(desc))
@btime $(cutoff1.f)(4+rand())
@btime NBodyIPs.fcut($cutoff1, 4+rand())
@btime NBodyIPs.fcut($(desc.cutoff), 4+rand())
@btime NBodyIPs.fcut($desc, 4+rand())
@btime $(desc.cutoff.f.obj.x)(4+rand())
@btime $(JuLIP.Potentials.F64fun(desc.cutoff.f.obj.x))(4+rand())

x = 4 .+ (@SVector rand(6)) .* 2
@btime NBodyIPs.fcut_d($cutoff1, $x)
@btime NBodyIPs.fcut($cutoff1, $x)
@btime NBodyIPs.fcut($desc, $x)
@btime NBodyIPs.fcut_d($desc, $x)

@btime NBodyIPs.transform($desc, 4+rand())

##
# # 3-body potential
# B3 = nbpolys(3, desc, 8)
# c = rand(length(B3))
# V3 = NBPoly(B3, c, desc)
# V3sp = StNBPoly(V3)
#
# (R, ii, J) = (SArray{Tuple{3},Float64,1,3}[[1.58988, 1.58644, 1.58548], [3.16204, -0.00261097, -0.00218217], [-4.73656, 1.58054, 1.57758], [4.74344, 1.58054, 1.57758], [-3.15081, 0.00344743, 0.000898426], [-1.5723, 1.58215, 1.57658], [0.00678987, 3.16114, 0.00410215], [1.58579, -4.74039, 1.58146], [1.58579, 4.73961, 1.58146], [3.1617, 3.15795, 0.00452045], [-3.15924, 3.1666, 0.00571087], [-1.57802, -4.74068, 1.5814], [-1.57802, 4.73932, 1.5814], [0.00726609, -3.16374, 0.00326648], [1.5791, -1.57564, 1.57527], [3.15877, -3.15342, 0.00183828], [-4.73569, -1.58259, 1.58014], [4.74431, -1.58259, 1.58014], [-3.1615, -3.16062, 0.00499946], [-1.57353, -1.57909, 1.57556], [0.00953652, 0.00570647, 3.1631], [1.58575, 1.58157, -4.73906], [1.58575, 1.58157, 4.74094], [3.16764, 0.00644966, 3.16233], [-3.1569, -0.00360474, 3.15531], [-1.57745, 1.58326, -4.7442], [-1.57745, 1.58326, 4.7358], [-0.00115697, 3.16279, 3.15497], [3.1669, 3.15745, 3.16212], [-3.15976, 3.16027, 3.15724], [0.000180625, -3.158, 3.15627], [1.5873, -1.57522, -4.73615], [1.5873, -1.57522, 4.74385], [3.16105, -3.15639, 3.15749], [-3.16145, -3.15836, 3.16352], [-1.58108, -1.57531, -4.74534], [-1.58108, -1.57531, 4.73466], [0.00234916, -0.00279929, -3.15585], [1.58104, 1.5818, -1.57899], [3.16084, 0.00199752, -3.15851], [-4.73986, 1.5849, -1.58027], [4.74014, 1.5849, -1.58027], [-3.15168, 0.00143851, -3.15614], [-1.57103, 1.58238, -1.58417], [0.00348743, 3.1621, -3.15689], [1.5891, -4.74034, -1.57682], [1.5891, 4.73966, -1.57682], [3.16149, 3.16205, -3.15535], [-3.15445, 3.16553, -3.16449], [-1.57679, -4.73654, -1.585], [-1.57679, 4.74346, -1.585], [0.00439256, -3.15483, -3.15568], [1.58515, -1.58148, -1.57721], [3.16848, -3.16355, -3.15797], [-4.73615, -1.57912, -1.58197], [4.74385, -1.57912, -1.58197], [-3.15364, -3.15896, -3.16494], [-1.5789, -1.58257, -1.58277]], 1, [57, 58])
# J = SVector(J...)
#
# @btime NBodyIPs.evaluate($V3sp, $desc, $R, $ii, $J)
#
# @code_warntype NBodyIPs.evaluate(V3sp, desc, R, ii, J)
#
# # function evaluate(V::NBodyFunction{N},
# #                   desc::NBSiteDescriptor,
# #                   Rs::AbstractVector{JVec{T}},
# #                   i::Int,
# #                   J::SVector{K, Int}) where {N, T, K}
# #    # get the physical descriptor: bond-lengths (+ bond-angles)
# #    rθ = ricoords(desc, Rs, J)
# #    # check whether to skip this N-body term?
# #    skip_simplex(desc, rθ) && return zero(T)
# #    # compute the cutoff (and skip this site if the cutoff is zero)
# #    fc = fcut(desc, rθ)
# #    # @btime fcut($desc, $rθ)
# #    # @show typeof(fc), fc
# #    fc == 0 && return zero(T)
# #    # compute the invariants (this also applies the transform)
# #    II = invariants(desc, rθ)
# #    # @btime invariants($desc, $rθ)
# #    # @show typeof(II), II
# #    # exit()
# #    # evaluate the inner potential function (e.g. polynomial)
# #    return evaluate_I(V, II) * fc
# # end
#
# # @btime NBodyIPs.ricoords($desc, $R, $J) # ok
# rθ = NBodyIPs.ricoords(desc, R, J)
# # @btime NBodyIPs.skip_simplex($desc, $rθ)  # ok
# fc = NBodyIPs.fcut(desc, rθ)
# @btime NBodyIPs.fcut($desc, $rθ)
# II = NBodyIPs.invariants(desc, rθ)
# @btime NBodyIPs.invariants($desc, $rθ)
# @btime NBodyIPs.evaluate_I($V3sp, $II)
#
# x = @SVector rand(3)
# @btime NBodyIPs.invariants($desc, $x)

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
function ensure_seeds_dimensions(model::Model{D, T}, initial_points::Vector{Vector{T}}) where {D, T}
    # Querying number of unknowns
    n_psi = dim_unknowns(model)

    n_initial_points = length(initial_points)

    for (i, initial_point) in enumerate(initial_points)
        if length(initial_point) != n_psi
            error("Seed $i has $(length(initial_point)) elements and the model has $n_psi unknowns.")
        end
    end
    return
end

function create_initial_points(model::Model{D, T}, n_initial_points::Int, LB::T, UB::T) where {D, T}
    # Querying number of unknowns
    n_psi = dim_unknowns(model)

    initial_points = Vector{Vector{T}}(undef, n_initial_points)

    # Generate initial values in [-1e3, 1e3]
    for i = 1:n_initial_points
        initial_points[i] = rand(Distributions.Uniform(LB, UB), n_psi)
    end
    return initial_points
end

function optimize(func::Optim.TwiceDifferentiable, opt_method::AbstractOptimizationMethod{T},
                  optimizer::Optim.AbstractOptimizer, verbose::Int, i::Int, time_limit_sec::Int) where T

    return Optim.optimize(func, opt_method.initial_points[i], optimizer,
                                                     Optim.Options(f_tol = opt_method.f_tol,
                                                                   g_tol = opt_method.g_tol,
                                                                   iterations = opt_method.iterations,
                                                                   show_trace = show_trace(verbose),
                                                                   time_limit = time_limit_sec))
end

function optimize(func::Optim.TwiceDifferentiable, opt_method::AbstractOptimizationMethod{T},
                  cons::Optim.TwiceDifferentiableConstraints,
                  optimizer::Optim.AbstractOptimizer,
                  verbose::Int, i::Int, time_limit_sec::Int) where T

    return Optim.optimize(func, cons, opt_method.initial_points[i], optimizer,
                                                     Optim.Options(f_tol = opt_method.f_tol,
                                                                   g_tol = opt_method.g_tol,
                                                                   iterations = opt_method.iterations,
                                                                   show_trace = show_trace(verbose),
                                                                   time_limit = time_limit_sec))
end

show_trace(verbose::Int) =  (verbose == 3 ? true : false)


text:
# # Example 1 : Regular Grid
#
# - two-dimensional fields on a regular `longitude,latitude` grid are retrieved from binary files
# - they are then written to a new netcdf file using `NCTiles.write()`
#
# This example illustrates the use of either `NCDatasets.jl` or `NetCDF.jl` as the backend.

using NCTiles, NCDatasets, NetCDF

# ### File Paths
inputs=NCTiles.NCTILES_TESTCASES
NCTiles.ensure_testcases_installed()

outputs = joinpath(tempdir(),"NCTILES_TESTCASES_OUTPUT/")
if ~ispath(outputs); mkpath(outputs); end

selectfields = ["ETAN"]
indir = joinpath(inputs,"diags_interp",selectfields[1])
fnames = joinpath.(Ref(indir),filter(x -> occursin(".data",x), readdir(indir)))

savedir = joinpath(outputs,"ex1")
if ~ispath(savedir); mkpath(savedir); end

# ### Grid specifications & Metadata

prec = Float32
README = readlines(joinpath(inputs,"README"))
lon=-179.75:0.5:179.75; lat=-89.75:0.5:89.75;
n1,n2 = (length(lon),length(lat))
tim_units = "days since 1992-01-01"
tim = vec([14.0 45.0 74.0])
longname = "Surface Height Anomaly"
units = "m"

# ### 1. Using `NCDatasets.jl` As The Backend

# Define Dimensions as `NCvar`s
dims = [NCvar("lon_c","degrees_east",size(lon),lon,Dict("long_name" => "longitude"),NCDatasets),
        NCvar("lat_c","degrees_north",size(lat),lat,Dict("long_name" => "latitude"),NCDatasets),
        NCvar("tim",tim_units,Inf,tim,Dict(("long_name" => "tim","standard_name" => "time")),NCDatasets)
        ]

# Define Variable Array as `BinData`
fielddata = BinData(fnames,prec,(n1,n2))
field = NCvar(selectfields[1],units,dims,fielddata,Dict("long_name" => longname),NCDatasets)

# Create the NetCDF file and populate with dimensions and variable
write(field,joinpath(savedir,"ex1_NCDatasets.nc"),README=README)

# ### 2. Using `NetCDF.jl` As The Backend

# Define Dimensions as `NCvar`s
dims = [NCvar("lon_c","degrees_east",size(lon),lon,Dict("long_name" => "longitude"),NetCDF),
        NCvar("lat_c","degrees_north",size(lat),lat,Dict("long_name" => "latitude"),NetCDF),
        NCvar("tim",tim_units,Inf,tim,Dict(("long_name" => "time","standard_name" => "time")),NetCDF)
        ]

# Define Variable Array as `BinData`
fielddata = BinData(fnames,prec,(n1,n2))
field = NCvar(selectfields[1],units,dims,fielddata,Dict("long_name" => longname),NetCDF)

write(field,joinpath(savedir,"ex1_NetCDF.nc"),README=README)

# _Note: the `write` function is a shorthand for_
#
# ```
# # Create the NetCDF file and populate with dimension and field info, as well as dimension data
# ncfile,fldvar,dimlist = createfile(joinpath(savedir,"ex1_NetCDF.nc"),field,README)
#
# # Add field data
# addData(fldvar,field)
#
# # Close the file-  only needed for NetCDF v < 0.10
# finalize(ncfile)
# ```


text:
using SIMD

@inline function vdot(xs, ys, ::Type{Vec{N,T}} = Vec{4,Float64}) where {N,T}
    @assert length(ys) == length(xs)
    @assert length(xs) % N == 0
    lane = VecRange{N}(0)
    vacc = zero(Vec{N, T})
    @inbounds for i in 1:N:length(xs)
        vacc = muladd(xs[lane + i], ys[lane + i], vacc)
    end
    return sum(vacc)
end


text:
# import Pkg
# Pkg.add("DataFrames")
# Pkg.add("CSV")
# Pkg.add("Distributions")
# Pkg.add("Random")
# Pkg.add("StatsBase")
# Pkg.add("SpecialFunctions")
# Pkg.add("Plots")

using DataFrames
using CSV
using Distributions
using Random
using StatsBase
using SpecialFunctions
using Plots

#testing protocol definitions
const PCR_mass_protocol, LFD_mass_protocol = 1,2
#LFD sens definitions
const SC_data_2021, SC_data_Oct, porton_down = 1,2,3
#Infectivity model definitions
const ke_inf_model_no, flat_inf_model_no, linear_inf_model_no = 1,2,3
#VL model definitions
const ke_model_no, kissler_model_no = 1,2
#options definitons for kissler model
const V0_model_opt, PVT_model_opt = 1,2
#options for infectivity model
const marks_peakinf_opt, unif_peakinf_opt = 1,2

#default options
VL_model = ke_model_no         #choice of viral load model (Ke et al is default)
LFD_model = SC_data_2021         #choice of LFD sensitivity model (social care binned 2021 data is default)
PCR_sens_max = 0.83           #max PCR sensitivity (Ferretti et al 2021)
Inf_model = ke_inf_model_no    #infectivity model
p_asymp = 0.5                  #asymptomatic fraction
onset_opt = V0_model_opt        #IF using kissler VL model, option for how to generate onset time
peak_inf_opt = marks_peakinf_opt   #IF using flat or linear inf model, what function is used for peak
PCR_TaT_scale = 1.0    #scale PCR delays by this factor

#===================== Symptoms onset timing =====================#

const symp_beta = 4.84 / 2.6^2
const symp_alpha = 4.84 * symp_beta

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

include("Kissler_VL_model.jl")
include("Ke_VL_inf_models.jl")
include("testing_models.jl")
include("Alt_inf_models.jl")

"""
    generate_asymptomatic()

Generate and return a boolean value determining whether
somebody will develop symptoms.

## Returns:
`Bool` = whether or not the person in asymptomatic (true =
          asymptomatic, false = symptoms)
"""
function generate_asymptomatic()
    #randomly generate if people are asymptomatic
    return (rand() < p_asymp)
end


text:
include("utils.jl")

@testset "Euclidean" begin
    M = Manifolds.Euclidean(3)
    types = [Vector{Float64},
             SizedVector{3, Float64},
             MVector{3, Float64},
             Vector{Float32},
             SizedVector{3, Float32},
             MVector{3, Float32},
             Vector{Double64},
             MVector{3, Double64},
             SizedVector{3, Double64}]
    for T in types
        @testset "Type $T" begin
            pts = [convert(T, [1.0, 0.0, 0.0]),
                   convert(T, [0.0, 1.0, 0.0]),
                   convert(T, [0.0, 0.0, 1.0])]
            test_manifold(M,
                          pts,
                          test_reverse_diff = isa(T, Vector),
                          test_project_tangent = true,
                          test_musical_isomorphisms = true,
                          test_vector_transport = true,
                          test_mutating_rand = isa(T, Vector),
                          point_distributions = [Manifolds.projected_distribution(M, Distributions.MvNormal(zero(pts[1]), 1.0))],
                          tvector_distributions = [Manifolds.normal_tvector_distribution(M, pts[1], 1.0)])
        end
    end
end