text:
### Edge and Gradient related tests ###

EPS = 1e-14

## Checkerboard array, used to test image gradients

let
    white{T}(::Type{T}) = one(T)
    black{T}(::Type{T}) = zero(T)
    white{T<:Unsigned}(::Type{T}) = typemax(T)
    black{T<:Unsigned}(::Type{T}) = typemin(T)

    global checkerboard
    function checkerboard{T}(::Type{T}, sq_width::Integer, count::Integer)
        wh = fill(white(T), (sq_width,sq_width))
        bk = fill(black(T), (sq_width,sq_width))
        bw = [wh bk; bk wh]
        vert = repmat(bw, (count>>1), 1)
        isodd(count) && (vert = vcat(vert, [wh bk]))
        cb = repmat(vert, 1, (count>>1))
        isodd(count) && (cb = hcat(cb, vert[:,1:sq_width]))
        cb
    end

    checkerboard(sq_width::Integer, count::Integer) = checkerboard(Uint8, sq_width, count)
end

SZ=5

cb_array    = checkerboard(SZ,3)
cb_image_xy = grayim(cb_array)
cb_image_yx = grayim(cb_array)
cb_image_yx["spatialorder"] = ["y","x"]

for method in ["sobel", "prewitt", "ando3", "ando4", "ando5", "ando4_sep", "ando5_sep"]
    ## Checkerboard array

    (agx, agy) = imgradients(cb_array, method)
    amag = magnitude(agx, agy)
    agphase = phase(agx, agy)
    @assert (amag, agphase) == magnitude_phase(agx, agy)

    @assert agx[1,SZ]   < 0.0   # white to black transition
    @assert agx[1,2*SZ] > 0.0   # black to white transition
    @assert agy[SZ,1]   < 0.0   # white to black transition
    @assert agy[2*SZ,1] > 0.0   # black to white transition

    # Test direction of increasing gradient
    @assert cos(agphase[1,SZ])   - (-1.0) < EPS   # increasing left  (=  pi   radians)
    @assert cos(agphase[1,2*SZ]) -   1.0  < EPS   # increasing right (=   0   radians)
    @assert sin(agphase[SZ,1])   -   1.0  < EPS   # increasing up    (=  pi/2 radians)
    @assert sin(agphase[2*SZ,1]) - (-1.0) < EPS   # increasing down  (= -pi/2 radians)

    # Test that orientation is perpendicular to gradient
    aorient = orientation(agx, agy)
    @assert all((cos(agphase).*cos(aorient) .+ sin(agphase).*sin(aorient) .< EPS) |
                ((agphase .== 0.0) & (aorient .== 0.0)))  # this part is where both are
                                                          # zero because there is no gradient

    ## Checkerboard Image with row major order

    (gx, gy) = imgradients(cb_image_xy, method)
    mag = magnitude(gx, gy)
    gphase = phase(gx, gy)
    @assert (mag, gphase) == magnitude_phase(gx, gy)

    @assert gx[SZ,1]   < 0.0   # white to black transition
    @assert gx[2*SZ,1] > 0.0   # black to white transition
    @assert gy[1,SZ]   < 0.0   # white to black transition
    @assert gy[1,2*SZ] > 0.0   # black to white transition

    @assert cos(gphase[SZ,1])   - (-1.0) < EPS   # increasing left  (=  pi   radians)
    @assert cos(gphase[2*SZ,1]) -   1.0  < EPS   # increasing right (=   0   radians)
    @assert sin(gphase[1,SZ])   -   1.0  < EPS   # increasing up    (=  pi/2 radians)
    @assert sin(gphase[1,2*SZ]) - (-1.0) < EPS   # increasing down  (= -pi/2 radians)

    # Test that orientation is perpendicular to gradient
    orient = orientation(gx, gy)
    @assert all((cos(gphase).*cos(orient) .+ sin(gphase).*sin(orient) .< EPS) |
                ((gphase .== 0.0) & (orient .== 0.0)))  # this part is where both are
                                                        # zero because there is no gradient

    ## Checkerboard Image with column-major order

    (gx, gy) = imgradients(cb_image_yx, method)
    mag = magnitude(gx, gy)
    gphase = phase(gx, gy)
    @assert (mag, gphase) == magnitude_phase(gx, gy)

    @assert gx[1,SZ]   < 0.0   # white to black transition
    @assert gx[1,2*SZ] > 0.0   # black to white transition
    @assert gy[SZ,1]   < 0.0   # white to black transition
    @assert gy[2*SZ,1] > 0.0   # black to white transition

    # Test direction of increasing gradient
    @assert cos(gphase[1,SZ])   - (-1.0) < EPS   # increasing left  (=  pi   radians)
    @assert cos(gphase[1,2*SZ]) -   1.0  < EPS   # increasing right (=   0   radians)
    @assert sin(gphase[SZ,1])   -   1.0  < EPS   # increasing up    (=  pi/2 radians)
    @assert sin(gphase[2*SZ,1]) - (-1.0) < EPS   # increasing down  (= -pi/2 radians)

    # Test that orientation is perpendicular to gradient
    orient = orientation(gx, gy)
    @assert all((cos(gphase).*cos(orient) .+ sin(gphase).*sin(orient) .< EPS) |
                ((gphase .== 0.0) & (orient .== 0.0)))  # this part is where both are
                                                        # zero because there is no gradient

end

# Create an image with white along diagonals -2:2 and black elsewhere
m = zeros(Uint8, 20,20)
for i = -2:2; m[diagind(m,i)] = 0xff; end

m_xy = grayim(m')
m_yx = grayim(m)
m_yx["spatialorder"] = ["y","x"]

for method in ["sobel", "prewitt", "ando3", "ando4", "ando5", "ando4_sep", "ando5_sep"]
    ## Diagonal array

    (agx, agy) = imgradients(m, method)
    amag = magnitude(agx, agy)
    agphase = phase(agx, agy)
    @assert (amag, agphase) == magnitude_phase(agx, agy)

    @assert agx[7,9]  < 0.0   # white to black transition
    @assert agx[10,8] > 0.0   # black to white transition
    @assert agy[10,8] < 0.0   # white to black transition
    @assert agy[7,9]  > 0.0   # black to white transition

    # Test direction of increasing gradient
    @assert abs(agphase[10,8] -    pi/4 ) < EPS   # lower edge (increasing up-right  =   pi/4 radians)
    @assert abs(agphase[7,9]  - (-3pi/4)) < EPS   # upper edge (increasing down-left = -3pi/4 radians)

    # Test that orientation is perpendicular to gradient
    aorient = orientation(agx, agy)
    @assert all((cos(agphase).*cos(aorient) .+ sin(agphase).*sin(aorient) .< EPS) |
                ((agphase .== 0.0) & (aorient .== 0.0)))  # this part is where both are
                                                          # zero because there is no gradient

    ## Diagonal Image, row-major order

    (gx, gy) = imgradients(m_xy, method)
    mag = magnitude(gx, gy)
    gphase = phase(gx, gy)
    @assert (mag, gphase) == magnitude_phase(gx, gy)

    @assert gx[9,7]  < 0.0   # white to black transition
    @assert gx[8,10] > 0.0   # black to white transition
    @assert gy[8,10] < 0.0   # white to black transition
    @assert gy[9,7]  > 0.0   # black to white transition

    # Test direction of increasing gradient
    @assert abs(gphase[8,10] -    pi/4 ) < EPS   # lower edge (increasing up-right  =   pi/4 radians)
    @assert abs(gphase[9,7]  - (-3pi/4)) < EPS   # upper edge (increasing down-left = -3pi/4 radians)

    # Test that orientation is perpendicular to gradient
    orient = orientation(gx, gy)
    @assert all((cos(gphase).*cos(orient) .+ sin(gphase).*sin(orient) .< EPS) |
                ((gphase .== 0.0) & (orient .== 0.0)))  # this part is where both are
                                                        # zero because there is no gradient

    ## Diagonal Image, column-major order

    (gx, gy) = imgradients(m_yx, method)
    mag = magnitude(gx, gy)
    gphase = phase(gx, gy)
    @assert (mag, gphase) == magnitude_phase(gx, gy)

    @assert gx[7,9]  < 0.0   # white to black transition
    @assert gx[10,8] > 0.0   # black to white transition
    @assert gy[10,8] < 0.0   # white to black transition
    @assert gy[7,9]  > 0.0   # black to white transition

    # Test direction of increasing gradient
    @assert abs(gphase[10,8] -    pi/4 ) < EPS   # lower edge (increasing up-right  =   pi/4 radians)
    @assert abs(gphase[7,9]  - (-3pi/4)) < EPS   # upper edge (increasing down-left = -3pi/4 radians)

    # Test that orientation is perpendicular to gradient
    orient = orientation(gx, gy)
    @assert all((cos(gphase).*cos(orient) .+ sin(gphase).*sin(orient) .< EPS) |
                ((gphase .== 0.0) & (orient .== 0.0)))  # this part is where both are
                                                        # zero because there is no gradient
end

# Nonmaximal suppression

function thin_edges(img)
    # Get orientation
    gx,gy = imgradients(img)
    orient = phase(gx,gy)

    # Do NMS thinning
    thin_edges_nonmaxsup_subpix(img, orient, radius=1.35)
end

function nms_test_horiz_vert(img, which)
    ## which = :horizontal or :vertical

    # Do NMS thinning
    t,s = thin_edges(img)

    # Calc peak location by hand

    # Interpolate values 1.35 pixels left and right
    # Orientation is zero radians -> to the right
    v1 = 6 - 0.35   # slope on right is 1
    v2 = 5 - 0.35*2 # slope on left is 2
    c = 7.0         # peak value

    # solve v = a*r^2 + b*r + c
    a = (v1 + v2)/2 - c
    b = a + c - v2
    r = -b/2a

    @assert abs(r - 1/6) < EPS

    # Location and value at peak
    peakloc = r*1.35 + 3
    peakval = a*r^2 + b*r + c

    transposed = spatialorder(img)[1] == "x"
    horizontal = which == :horizontal

    test_axis1 = transposed $ !horizontal

    @assert test_axis1 ? all(t[:,[1,2,4,5]] .== 0) : all(t[[1,2,4,5],:] .== 0)
    @assert test_axis1 ? all(t[:,3]   .== peakval) : all(t[3,:]   .== peakval)
    @assert test_axis1 ? all(s[:,[1,2,4,5]] .== zero(Base.Graphics.Point)) :
                         all(s[[1,2,4,5],:] .== zero(Base.Graphics.Point))

    if transposed
        if which == :horizontal
            @assert     [pt.x for pt in s[:,3]]  == [1:5]
            @assert all([pt.y for pt in s[:,3]] .== peakloc)
        else
            @assert all([pt.x for pt in s[3,:]] .== peakloc)
            @assert     [pt.y for pt in s[3,:]]  == [1:5]
        end
    else
        if which == :horizontal
            @assert     [pt.x for pt in s[3,:]]  == [1:5]
            @assert all([pt.y for pt in s[3,:]] .== peakloc)
        else
            @assert all([pt.x for pt in s[:,3]] .== peakloc)
            @assert     [pt.y for pt in s[:,3]]  == [1:5]
        end
    end
end

# Test image: vertical edge
m = [3.0  5.0  7.0  6.0  5.0
     3.0  5.0  7.0  6.0  5.0
     3.0  5.0  7.0  6.0  5.0
     3.0  5.0  7.0  6.0  5.0
     3.0  5.0  7.0  6.0  5.0]

m_xy = grayim(m')
m_yx = grayim(m)
m_yx["spatialorder"] = ["y","x"]

nms_test_horiz_vert(m, :vertical)
nms_test_horiz_vert(m_xy, :vertical)
nms_test_horiz_vert(m_yx, :vertical)

# Test image: horizontal edge
m = m'
m_xy = grayim(m')
m_yx = grayim(m)
m_yx["spatialorder"] = ["y","x"]

nms_test_horiz_vert(m, :horizontal)
nms_test_horiz_vert(m_xy, :horizontal)
nms_test_horiz_vert(m_yx, :horizontal)

function nms_test_diagonal(img)
    # Do NMS thinning
    t,s = thin_edges(img)

    # Calc peak location by hand

    # Interpolate values 1.35 pixels up and left, down and right
    # using bilinear interpolation
    # Orientation is π/4 radians -> 45 degrees up
    fr = 1.35*cos(π/4)
    lower = (7 + fr*(6-7))
    upper = (6 + fr*(5-6))
    v1 = lower + fr*(upper-lower)

    lower = (7 + fr*(5-7))
    upper = (5 + fr*(3-5))
    v2 = lower + fr*(upper-lower)

    c = 7.0         # peak value

    # solve v = a*r^2 + b*r + c
    a = (v1 + v2)/2 - c
    b = a + c - v2
    r = -b/2a

    @assert (r - 1/6) < EPS

    transposed = spatialorder(img)[1] == "x"

    # Location and value at peak

    x_peak_offset, y_peak_offset = r*fr, -r*fr
    peakval = a*r^2 + b*r + c

    @assert all(diag(data(t))[2:4] .== peakval)  # Edge pixels aren't interpolated here
    @assert all(t - diagm(diag(data(t))) .== 0)

    diag_s = copy(s, diagm(diag(data(s))))
    @assert s == diag_s

    @assert all([pt.x for pt in diag(data(s))[2:4]] - ([2:4] + x_peak_offset) .< EPS)
    @assert all([pt.y for pt in diag(data(s))[2:4]] - ([2:4] + y_peak_offset) .< EPS)

end

# Test image: diagonal edge
m = [7.0  6.0  5.0  0.0  0.0
     5.0  7.0  6.0  5.0  0.0
     3.0  5.0  7.0  6.0  5.0
     0.0  3.0  5.0  7.0  6.0
     0.0  0.0  3.0  5.0  7.0]

m_xy = grayim(m')
m_yx = grayim(m)
m_yx["spatialorder"] = ["y","x"]

nms_test_diagonal(m)
nms_test_diagonal(m_xy)
nms_test_diagonal(m_yx)


text:
module EvoTrees

export init_evotree, grow_evotree!, grow_tree, predict, fit_evotree,
    EvoTreeRegressor, EvoTreeCount, EvoTreeClassifier, EvoTreeGaussian,
    EvoTreeRModels, importance

using Statistics
using Base.Threads: @threads
using StatsBase: sample, quantile
import StatsBase: predict
using Random: seed!
using StaticArrays
using Distributions
using CategoricalArrays
import MLJBase
# import MLJ

include("models.jl")
include("structs.jl")
include("loss.jl")
include("eval.jl")
include("predict.jl")
include("find_split.jl")
include("fit.jl")
include("importance.jl")
include("MLJ.jl")

end # module

ext: jl


text:
module QPointCellFieldsTests

using Test
using FillArrays
using Gridap.Helpers
using Gridap.Fields
using Gridap.CellData
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.Integration

using Gridap.Fields: MockField, MockBasis, OtherMockBasis

ndofs = 4
c = fill(1.0,ndofs)
f = OtherMockBasis{2}(ndofs)

l = 10
cl = fill(c,l)
fl = Fill(f,l)
ϕl = lincomb(fl,cl)

degree = 3
quad = CellQuadrature(degree,[QUAD],Fill(1,l))
q = get_coordinates(quad)

s_q = [i*ones(size(qi)) for (i,qi) in enumerate(q)]
a = CellData.ArrayOfEvaluatedFields(s_q,q)
test_array_of_fields(a,q,s_q)

v = VectorValue(3.0,4.0)
cf = QPointCellField(v,ϕl,quad)
cf_q = evaluate(cf,q)
test_cell_field(cf,q,cf_q)

cf = CellField(v,ϕl,quad)
cf_q = evaluate(cf,q)
test_cell_field(cf,q,cf_q)

cf = CellField(1.0,ϕl,quad)
@test sum(integrate(cf,ϕl,quad)) ≈ 213.33333333333323

a = CellField(1.0,ϕl,quad)
b = CellField(1.0,ϕl,quad)
c = CellField(1.0,ϕl,quad)

function updater(a,b,c)
  b_new = 2*b
  c_new = 4*c
  true, b_new, c_new
end

@test evaluate(a,q)[1] == fill(1.0,4)
@test evaluate(b,q)[1] == fill(1.0,4)
@test evaluate(c,q)[1] == fill(1.0,4)

update_state_variables!(updater,quad,a,b,c)
@test evaluate(a,q)[1] == fill(1.0,4)
@test evaluate(b,q)[1] == fill(2.0,4)
@test evaluate(c,q)[1] == fill(4.0,4)

update_state_variables!(updater,quad,a,b,c)
@test evaluate(a,q)[1] == fill(1.0,4)
@test evaluate(b,q)[1] == fill(4.0,4)
@test evaluate(c,q)[1] == fill(16.0,4)

end # module


text:
module FunctionZeros
import SpecialFunctions
import Roots

export besselj_zero, bessely_zero

besselj_zero_asymptotic(nu, n) = bessel_zero_asymptotic(nu, n, 1)

"""
    bessel_zero_asymptotic(nu, n, kind=1)

Asymptotic formula for the `n`th zero of the the Bessel J (Y) function of order `nu`.
`kind == 1 (2)` for Bessel function of the first (second) kind, J (Y).
"""
function bessel_zero_asymptotic(nu_in::Real, n::Integer, kind=1)
    nu = abs(nu_in)
    if kind == 1
        beta = MathConstants.pi * (n + nu / 2 - 1//4)
    else # kind == 2
        beta = MathConstants.pi * (n + nu / 2 - 3//4)
    end
    delta = 8 * beta
    mu = 4 * nu^2
    mup2 = mu * mu
    mup3 = mup2 * mu
    deltap2 = delta * delta
    deltap3 = deltap2 * delta
    deltap4 = deltap2 * deltap2
    deltap6 = deltap3 * deltap3
    t1 = 1
    t2 = 4 * (7 * mu - 31) / (3 * deltap2)
    t3 = 32 * (83 * mup2 - 982 * mu + 3779) / (15 * deltap4)
    t4 = 64 * (6949 * mup3 - 153855 * mup2 + 1585743 * mu - 6277237) /
        (105 * deltap6)
    zero_asymp = beta - ((mu - 1) / delta) * (t1 + t2 + t3 + t4)
    return zero_asymp
end

# Use the asymptotic values as starting values.
# These find the correct zeros even for n = 1,2,...
# Order 0 is 6 times slower and 50-100 times less accurate
# than higher orders, with other parameters constant.
"""
    besselj_zero(nu, n; order=2)

`n`th zero of the Bessel J function of order `nu`,
for `n` = `1,2,...`.

`order` is passed to the function `Roots.fzero`.
"""
besselj_zero(nu, n; order=2) = Roots.fzero((x) -> SpecialFunctions.besselj(nu, x),
                                           bessel_zero_asymptotic(nu, n, 1); order=order)

"""
    bessely_zero(nu, n; order=2)

`n`th zero of the Bessel Y function of order `nu`,
for `n` = `1,2,...`.

`order` is passed to the function `Roots.fzero`.
"""
bessely_zero(nu, n; order=2) = Roots.fzero((x) -> SpecialFunctions.bessely(nu, x),
                                           bessel_zero_asymptotic(nu, n, 2); order=order)

end # module FunctionZeros


text:
#script for pre-processing data. Ends with giving Jacard similarity between
#ground-truth world states and the mode inferred world state, and 95% CI on that

using CSV
using DataFrames
using JSON
using Pipe: @pipe
using MetaGen
using Bootstrap

include("helper_function.jl")

#path = "../Data/20particles_threshold64_18945877/"
path = "../../scratch_work_07_16_21/09_18/shuffle_0/"

online_data = CSV.read(path * "online_ws.csv", DataFrame; delim = "&")
retro_data = CSV.read(path * "retro_ws.csv", DataFrame; delim = "&")
lesioned_data = CSV.read(path * "lesioned_ws.csv", DataFrame; delim = "&")

################################################################################
#could equally use input or output dictionary
dict = @pipe path * "output.json" |> open |> read |> String |> JSON.parse
#dict = @pipe "../../scratch_work_07_16_21/08_20/data_labelled.json" |> open |> read |> String |> JSON.parse

num_videos = 100
num_frames = 20

num_particles = 100

num_training_videos = 50

ground_truth_world_states = get_ground_truth(dict, num_videos)

################################################################################
online_world_states = new_parse_data(online_data, num_training_videos, num_particles)
retrospective_world_states = new_parse_data(retro_data, num_videos, num_particles)
lesioned_world_states = new_parse_data(lesioned_data, num_videos, num_particles)

################################################################################
#get just object categories
ground_truth_categories = extract_category(ground_truth_world_states)
online_categories = extract_category(online_world_states, num_training_videos, num_particles)
retrospective_categories = extract_category(retrospective_world_states, num_videos, num_particles)
lesioned_categories = extract_category(lesioned_world_states, num_videos, num_particles)

#jaccard similarity
sim_online = jaccard_similarity(ground_truth_categories, online_categories, num_particles)
sim_retrospective = jaccard_similarity(ground_truth_categories, retrospective_categories, num_particles)
sim_lesioned = jaccard_similarity(ground_truth_categories, lesioned_categories, num_particles)

#get CI
sim_online, sim_online_lower_ci, sim_online_upper_ci = add_confidence_interval(sim_online, 1000, 0.95)
sim_retrospective, sim_retrospective_lower_ci, sim_retrospective_upper_ci = add_confidence_interval(sim_retrospective, 1000, 0.95)
sim_lesioned, sim_lesioned_lower_ci, sim_lesioned_upper_ci = add_confidence_interval(sim_lesioned, 1000, 0.95)

################################################################################
#make new dataframe. just has stuff on the similarity.

new_df = DataFrame(video = 1:num_videos,
    sim_online = vcat(sim_online, fill(NaN, num_videos - num_training_videos)),
    sim_online_lower_ci = vcat(sim_online_lower_ci, fill(NaN, num_videos - num_training_videos)),
    sim_online_upper_ci = vcat(sim_online_upper_ci, fill(NaN, num_videos - num_training_videos)),
    sim_retrospective = sim_retrospective,
    sim_retrospective_lower_ci = sim_retrospective_lower_ci,
    sim_retrospective_upper_ci = sim_retrospective_upper_ci,
    sim_lesioned = sim_lesioned,
    sim_lesioned_lower_ci = sim_lesioned_lower_ci,
    sim_lesioned_upper_ci = sim_lesioned_upper_ci,
    )
#new_df = DataFrame(video = 1:num_videos, sim_online = sim_online)
CSV.write(path * "similarity3D.csv", new_df)

################################################################################
#repeat Jaccard sim in 2D

include("scripts/useful_functions.jl")

params = Video_Params(n_possible_objects = 5)

sim_online, sim_online_lower_ci, sim_online_upper_ci = jacccard_sim_2D(num_training_videos,
    num_frames, params, dict, ground_truth_world_states, online_world_states,
    1000, 0.95)

sim_retrospective, sim_retrospective_lower_ci, sim_retrospective_upper_ci = jacccard_sim_2D(num_videos,
    num_frames, params, dict, ground_truth_world_states, retrospective_world_states,
    1000, 0.95)

sim_lesioned, sim_lesioned_lower_ci, sim_lesioned_upper_ci = jacccard_sim_2D(num_videos,
    num_frames, params, dict, ground_truth_world_states, lesioned_world_states,
    1000, 0.95)

sim_NN_fitted = jacccard_sim_2D(num_videos, num_frames, params, dict, ground_truth_world_states, threshold = 0.09)

sim_NN_input = jacccard_sim_2D(num_videos, num_frames, params, dict, ground_truth_world_states, threshold = 0.0, top_n = 5)

new_df = DataFrame(video = 1:num_videos,
    sim_online = vcat(sim_online, fill(NaN, num_videos - num_training_videos)),
    sim_online_lower_ci = vcat(sim_online_lower_ci, fill(NaN, num_videos - num_training_videos)),
    sim_online_upper_ci = vcat(sim_online_upper_ci, fill(NaN, num_videos - num_training_videos)),
    sim_retrospective = sim_retrospective,
    sim_retrospective_lower_ci = sim_retrospective_lower_ci,
    sim_retrospective_upper_ci = sim_retrospective_upper_ci,
    sim_lesioned = sim_lesioned,
    sim_lesioned_lower_ci = sim_lesioned_lower_ci,
    sim_lesioned_upper_ci = sim_lesioned_upper_ci,
    sim_NN_fitted = sim_NN_fitted,
    sim_NN_input = sim_NN_input)

CSV.write(path * "similarity2D.csv", new_df)

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
mutable struct MetaDiGraph{T<:Integer,U<:Real} <: AbstractMetaGraph{T}
    graph::SimpleDiGraph{T}
    vprops::Dict{T,PropDict}
    eprops::Dict{SimpleEdge{T},PropDict}
    gprops::PropDict
    weightfield::Symbol
    defaultweight::U
    metaindex::MetaDict
    indices::Set{Symbol}
end

function MetaDiGraph(x, weightfield::Symbol, defaultweight::U) where U <: Real
    T = eltype(x)
    g = SimpleDiGraph(x)
    vprops = Dict{T,PropDict}()
    eprops = Dict{SimpleEdge{T},PropDict}()
    gprops = PropDict()
    metaindex = MetaDict()
    idxs = Set{Symbol}()

    MetaDiGraph(g, vprops, eprops, gprops,
        weightfield, defaultweight,
        metaindex, idxs)
end

MetaDiGraph() = MetaDiGraph(SimpleDiGraph())
MetaDiGraph{T,U}() where T <: Integer where U <: Real = MetaDiGraph(SimpleDiGraph{T}(), one(U))
MetaDiGraph{T,U}(x::Integer) where T <: Integer where U <: Real = MetaDiGraph(T(x), :weight, U(1.0))
MetaDiGraph(x) = MetaDiGraph(x, :weight, 1.0)
MetaDiGraph(x, weightfield::Symbol) = MetaDiGraph(x, weightfield, 1.0)
MetaDiGraph(x, defaultweight::Real) = MetaDiGraph(x, :weight, defaultweight)

# converts MetaDiGraph{Int,Float64} to MetaDiGraph{UInt8, Float32}
function MetaDiGraph{T,U}(g::MetaDiGraph) where T<:Integer where U<:Real
    newg = SimpleDiGraph{T}(g.graph)
    return MetaDiGraph(newg, U(g.defaultweight))
end

function MetaDiGraph{T,U}(g::SimpleDiGraph, weightfield::Symbol=:weight, defaultweight::Real=1.0) where T<:Integer where U<:Real
    newg = SimpleDiGraph{T}(g)
    return MetaDiGraph(newg, weightfield, U(defaultweight))
end

function MetaDiGraph{T,U}(g::SimpleDiGraph, defaultweight::Real) where T<:Integer where U<:Real
    newg = SimpleDiGraph{T}(g)
    return MetaDiGraph(newg, :weight, U(defaultweight))
end

SimpleDiGraph(g::MetaDiGraph) = g.graph

is_directed(::Type{MetaDiGraph}) = true
is_directed(::Type{MetaDiGraph{T,U}}) where T where U = true
is_directed(g::MetaDiGraph) = true

weighttype(g::MetaDiGraph{T,U}) where T where U = U

props(g::MetaDiGraph, e::SimpleEdge) = get(g.eprops, e, PropDict())

function set_props!(g::MetaDiGraph, e::SimpleEdge, d::Dict)
    if has_edge(g, e)
        if !_hasdict(g, e)
            g.eprops[e] = d
        else
            merge!(g.eprops[e], d)
        end
        return true
    end
    return false
end

zero(g::MetaDiGraph{T,U}) where T where U = MetaDiGraph{T,U}(SimpleDiGraph{T}())


text:
#=
Given an undirected graph represented as an adjacency matrix and an integer k, write a function to determine whether each vertex in the graph can be colored such that no two adjacent vertices share the same color using at most k colors.
=#

using Test
include("Solutions/problem56_color_graph_k_colors.jl")

graph = [0 1 1 0;1 0 0 0;1 0 0 1;0 0 1 0]

@test colorable(graph, 1) == false
@test colorable(graph, 2)
@test colorable(graph, 3)

graph = [0 1 1 0 1;1 0 0 0 1;1 0 0 1 0;0 0 1 0 1;1 1 0 1 0]

@test colorable(graph, 1) == false
@test colorable(graph, 2) == false
@test colorable(graph, 3)

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using Base
using Test
using C2Julia
using C2Julia:value
using C2Julia.stdlib
const float = C2Julia.float

@testset "CNumber type promotion" begin
    @test promote_type(int, char) == int
    @test promote_type(char, unsigned_char) == int
    @test promote_type(unsigned_int, long) == unsigned_int
    @test promote_type(long_long, float) == float
end

@testset "CNumber type conversion" begin
    @test convert(int, float(3)) === int(3)
    @test convert(int, float(3)) == 3
    @test convert(float, int(3)) === float(3)
    @test convert(float, int(3)) == 3
    @test convert(long, double(1e9)) |> value == 1e9
    @test convert(unsigned_char, int(1000)) |> value == 1000 % Cuchar
    @test convert(Cint, int(2000)) == 2000
    @test convert(Float32, int(2000)) == 2000
end

@testset "CFunction" begin
    @C function foo(a::int, b::float)::int
        return @bool(a) ? a + b : b
    end
    @test foo(2, 2) === int(4)
    @test foo(3, 2) === int(5)
    @test foo(0, 2) === int(2)
end

mutable struct S1
    a::int
    b::double
end

@testset "@pointer" begin
    s = S1(2, 3)
    @test @pointer(s.a)[] == 2
    a::Pointer{int} = C2Julia.malloc(sizeof(int) * 10)
    a[] = 1
    a[1] = 2
    @test @pointer(a[0])[] == 1
    @test @pointer(a[1])[] == 2
    sscanf("1234", "%d", @pointer(a[1]))
    @test a[1] == 1234
end

@testset "FieldRef" begin
    function add_one(x)
        @post_inc x[]
    end
    x = int(123)
    s = S1(x, 2.3)
    @test s.a == x
    add_one(@pointer s.a)
    @test s.a == x + 1
end

@testset "Postfix operators" begin
    orig = 234
    x::int = orig
    y::int = @post_dec x
    @test y == orig
    @test x == orig - 1
end

to_string(buf) = buf |> pointer |> unsafe_string

@testset "printf" begin
    buf = Vector{UInt8}(undef, 128)
    hello = "hello, world!"
    sprintf(buf, hello)
    @test to_string(buf) == hello
    sprintf(buf, "%s", hello)
    @test to_string(buf) == hello
    sprintf(buf, "%f", 123.)
    @test parse(Float32, to_string(buf)) == 123.
    sprintf(buf, "%d", 123)
    @test parse(Int, to_string(buf)) == 123
    sprintf(buf, "%f", double(123.))
    @test parse(Float32, to_string(buf)) == 123.
    sprintf(buf, "%d", int(123))
    @test parse(Int, to_string(buf)) == 123
end

@testset "scanf" begin
    x = Ref{Cint}()
    sscanf("1234", "%d", x)
    @test x[] == 1234
    x = Ref{int}()
    sscanf("1234", "%d", x)
    @test x[] == 1234
    x = Ref{Cfloat}()
    sscanf("1234", "%f", x)
    @test x[] == 1234
    x = Ref{Cdouble}()
    sscanf("1234", "%lf", x)
    @test x[] == 1234
    buf = Vector{UInt8}(undef, 128)
    sscanf("\n1234", "%s", buf)
    @test to_string(buf) == "1234"
end


text:
hex_string(' ')

escape("t est\n")

u=URI("http://google.com")

sprint(writemime, mime("text/html"), u)


text:
module epstate

export EPState, SingleSpinState, MultiSpinState, spin_configs

"""
Abstract base class for a configuration of an exclusion process.

All deriving state classes have to define the field `spins` that encodes the
spin configuration.
"""
abstract type EPState end

"""
    SingleSpinState(state::Vector{Bool})

State of a single process. The `spins` field gives spin configuration as Booleans.
"""
struct SingleSpinState <: EPState
    spins::BitVector

    SingleSpinState(init_spins::BitVector) = new(init_spins)

    SingleSpinState(n_sites::Int) = new(bitrand(n_sites))

    function SingleSpinState(n_sites::Int, n_particles::Int)
        if n_particles > n_sites
            throw(ArgumentError("n_partices must be <= n_sites"))
        end
        new(random_spins(n_sites, n_particles))
    end

end

"""
    MultiSpinState(state::Vector{UInt})

State of several processes (up to 64). The `spins` field gives spin
configurations, encoded in the bitstring forming `UInt`.
"""
struct MultiSpinState <: EPState
    spins::Vector{UInt}

    MultiSpinState(init_spins::Vector{UInt}) = new(init_spins)

    MultiSpinState(n_sites::Int) = new(rand(UInt, n_sites))

    function MultiSpinState(n_sites::Int, n_particles::Int)
        if n_particles > n_sites
            throw(ArgumentError("n_partices must be <= n_sites"))
        end
        new(random_multispins(n_sites, n_particles))
    end

end

"""
Returns vector of `UInt64`s corresponding to 64 random configurations with fixed particle number.
"""
function random_multispins(n_sites::Int, n_particles::Int)

    configs = hcat([random_spins(n_sites, n_particles) for _ in 1:64]...)

    integer_configs = []

    for j in 1:n_sites
        bit_string = join(Vector{Int}(configs[j,:]))
        append!(integer_configs, parse(UInt64, bit_string, 2))
    end

    return integer_configs
end

"""
Returns `BitArray` with random configuration of fixed particle number
"""
random_spins(n_sites::Int, n_particles::Int) = shuffle([trues(n_particles);falses(n_sites - n_particles)])

"""
Converts `UInt64` to 64 bit `BitArray`
"""
uint2spins(ms::UInt64) = BitArray(split(bin(ms, 64), "") .== "1")

"""
Returns 2D array giving the 64 configurations encoded in a `MultiSpinState`
"""
spin_configs(mss::MultiSpinState) = hcat([uint2spins(ms) for ms in mss.spins]...)

end

max_stars_repo_name: AustenLamacraft/ExclusionProcesses

max_issues_repo_name: AustenLamacraft/ExclusionProcesses

max_forks_repo_name: AustenLamacraft/ExclusionProcesses


text:
function debugprint(ssa, cinfo, binfo)
    display(cinfo)
    debugprintssa(ssa, cinfo)
    debugprintblockinfo(binfo)
end
printslots(ci) = [println("_$i:  $(ci.slotnames[i])  $(ci.slottypes[i])") for i = 1:length(ci.slotnames)];
printssa(ssa) = [println("%$i:  ", itemname(ssa[i])) for i = 1:length(ssa)];
printssatypes(ssa,ci) = [println("$(ci.ssavaluetypes[i])\t%$i:  ", itemname(ssa[i])) for i = 1:length(ci.ssavaluetypes)];
# printssatypes(ci) = [println("%$i:  $(ci.ssavaluetypes[i])") for i = 1:length(ci.ssavaluetypes)];
itemname(item) = isa(item, GlobalRef) ? item.name : item
itemname(item::Array) = itemname.(item)

function debugprintssa(ssa, ci)
    println("\nSlots:")
    printslots(ci)
    println("\nSSA:")
    printssatypes(ssa, ci)
    println()
end

function debugprintblockinfo(binfo)
    println("Inferred blocktree (parents of each ssa index):")
    for (i, v) in enumerate(binfo.parents)
        if length(v) > 0
            println(i, ": ", v)
        else
            println(i, ": ")
        end
    end
    println()
end


text:
# Only run coverage from linux nightly build on travis.
get(ENV, "TRAVIS_OS_NAME", "")       == "linux"   || exit()
get(ENV, "TRAVIS_JULIA_VERSION", "") == "nightly" || exit()

using Coverage

cd(joinpath(dirname(@__FILE__), "..")) do
    const fileCoverage = process_folder();
    Coveralls.submit(fileCoverage);
    Codecov.submit(fileCoverage);
end


text:
using Dorayaki
using Random
using Plots
using BenchmarkTools

Random.seed!(0)

const lr = 0.1
const iters = 100

x = rand(1, 100)
y = 5 .+ 2 .* x .+ rand(1, 100)

function main(x, y)
    W = Var(zeros(1, 1))
    b = Var(zeros(1))

    predict(x) = W * x + b

    for i in 1:iters
        y_pred = predict(x)
        loss = mean_squared_error(y, y_pred)

        cleargrad!(W)
        cleargrad!(b)

        gradient!(loss)
        W.data .-= lr .* W.grad.data
        b.data .-= lr .* b.grad.data
        # println(W.data, b.data, loss.data)
    end
    return W.data[1], b.data[1]
end

ret = @benchmark main(x, y)
#=
f(a) = ret[1] * a + ret[2]
fig = plot(x, y, st=:scatter)
plot!(fig, f) =#


text:
export Layout, LForwardReference, LFixedVector, LVarVector
export is_layout_fixed, is_layout_variable, simple_size, total_size, Cserialize

"""
Layout

All structs used to describe the memory layout (of a C-data structure) need to be
subtypes of this.
Some controlling objects used in such templates to describe vectors and pointers
have also this type.
A `Layout` structure and a memory pointer are needed to construct an `CAccessor` object.
"""
abstract type Layout end

abstract type LVector{T} <: Layout end

Base.eltype(::Type{<:LVector{T}}) where T = T

# Layout Elements
"""
    LFixedVector{T,N}

Denote a fixed size vector with element type `T` and size `N`.
"""
struct LFixedVector{T,N} <: LVector{T}
    p::NTuple{N,T}
end
get_length(::Type{LFixedVector{T,N}}, ::Any) where {T,N} = N
Base.eltype(::Type{LFixedVector{T,N}}) where {T,N} = T

"""
    LVarVector{T,F}

Denote a variable length vector with element type `T` in a template.
`F` is a function, which calculates the length of the vector, given the
accessor object containing the vector.

Example:
    struct A <: Layout
        len::Int
        vec::NVarVector{Float64, (x) -> x.len}
    end

"""
struct LVarVector{T,F}  <: LVector{T}
    p::NTuple{1,T}
end

function get_length(::Type{LVarVector{T,F}}, x) where {T,F}
    F isa Symbol ? getproperty(x, F) :
    F isa Integer ? getindex(x, F) :
    F isa Function ? F(x) :
    0
end
function set_length!(::Type{LVarVector{T,F}}, t::Tuple{U,Ptr}, v) where {T,F,U}
    S, p = t
    ix = 1
    G = Nothing
    if F isa Symbol
        # setproperty!(x, v, F)
        S <: Layout || return v
        i = findfirst(isequal(F), fieldnames(S))
        i === nothing && return v
        G = fieldtype(S, i)
        p += fieldoffset(S, i)
    elseif F isa Integer
        # setindex!(x, v, F)
        S <: LFixedVector || return v
        G = eltype(S)
        ix = Int(F)
        ix <= 0 && return v
    end
    G === Nothing && return v
    vv = G(v)
    unsafe_store!(Ptr{G}(p), vv, ix)
    vv
end

Base.eltype(::Type{LVarVector{T,F}}) where {T,F} = T

struct LForwardReference{M,L} <: Layout
    p::Ptr{Nothing}
end
Base.eltype(::Type{LForwardReference{M,L}}) where {M,L} = M.name.module.eval(L)

const TEMPLATE_FIXED = true
const TEMPLATE_VAR = false
"""
    is_template_variable(type)

Has the layout described by `type` a variable size
(for example variable sized vector in last field of a struct)?
"""
is_layout_variable(T::Type, deep::Bool=false) = !is_layout_fixed(T, deep)

"""
    is_template_fixed(type)

Has the layout described by `type` a fixed size.
"""
is_layout_fixed(T::Type, deep::Bool=false) = is_layout_fixed(T, deep, Dict())
function is_layout_fixed(::Type{T}, deep::Bool, dup) where T
    isprimitivetype(T) || throw(ArgumentError("$T is not a supported layout type"))
    TEMPLATE_FIXED
end
function is_layout_fixed(::Type{S}, deep::Bool, dup) where {T,S<:Ptr{T}}
    T <: Ptr && throw(ArgumentError("$S is not a supported layout type"))
    get!(dup, S) do
        dup[S] = TEMPLATE_FIXED
        d = is_layout_fixed(T, deep, dup)
        deep ? d : TEMPLATE_FIXED
    end
end
function is_layout_fixed(::Type{S}, deep::Bool, dup) where {S<:LForwardReference}
    is_layout_fixed(Ptr{eltype(S)}, deep, dup)
end
function is_layout_fixed(::Type{S}, deep::Bool, dup) where {T,N,S<:LFixedVector{T,N}}
    get!(dup, S) do
        dup[S] = TEMPLATE_FIXED
        k = is_layout_fixed(T, deep, dup)
        if N > 1 && k == TEMPLATE_VAR
            throw(ArgumentError("$S with variable length elements"))
        end
        N == 0 ? TEMPLATE_FIXED : k
    end
end
function is_layout_fixed(::Type{S}, deep::Bool, dup) where {T,S<:LVarVector{T}}
    get!(dup, S) do
        dup[S] = TEMPLATE_VAR
        is_layout_fixed(T, deep, dup)
        TEMPLATE_VAR
    end
end
function is_layout_fixed(::Type{T}, deep::Bool, dup) where {T<:Layout}
    get!(dup, T) do
        k = dup[T] = TEMPLATE_FIXED
        if !isbitstype(T)
            text = isconcretetype(T) ? "bits" : "concrete"
            throw(ArgumentError("$T is not a $text type struct"))
        end
        fields = fieldnames(T)
        n = length(fields)
        for i = 1:n
            f = fields[i]
            F = fieldtype(T, f)
            k = is_layout_fixed(F, deep, dup)
            if i < n && k == TEMPLATE_VAR
                throw(ArgumentError("$F has variable length in '$T.$f' - not last field"))
            end
        end
        k
    end
end

function align(p::Integer, s::Integer=sizeof(Ptr)) # s must be 2^n
    t = s - 1
    (p + t )  & ~t
end

function align(p::Integer, ::Type{T}) where T
    align(p, Base.aligned_sizeof(T))
end

simple_size(T::Type, veclens) = blength(T, veclens, Val(false))
total_size(T::Type, veclens) = blength(T, veclens, Val(true))

function blength(::Type{T}, veclens, v::Val{P}) where {P,F,N,T<:LFixedVector{F,N}}
    s = sizeof(T)
    j = 0
    for _ = 1:N
        j, s = blength_helper(F, veclens, j, s, v, T)
    end
    if j < length(veclens)
        throw(ArgumentError("too many variable length specifiers for $T only $j are needed"))
    end
    s
end

function blength(::Type{T}, veclens, v::Val{P}) where {P,S,T<:LVarVector{S}}
    isempty(veclens) && return 0
    n = first(veclens)
    n == 0 && return 0
    n < 0 && throw(ArgumentError("negative vector length '$n' not allowed"))
    vl(i) = i < length(veclens) ? veclens[i+1] : ()
    sum(blength(S, vl(i), v) for i = 1:n)
end

blength(::Type{Ptr{T}}, veclens, v::Val{P}) where {P,T} = P ? blength(T, veclens, v) : 0

function blength(::Type{T}, veclens, v::Val{P}) where {P,T<:Layout}
    s = sizeof(T)
    j = 0
    for i = 1:fieldcount(T)
        F = fieldtype(T, i)
        j, s = blength_helper(F, veclens, j, s, v, T)
    end
    if j < length(veclens)
        throw(ArgumentError("too many variable length specifiers for $T- only $j are needed"))
    end
    s
end

blength(::Type{T}, veclens, ::Val) where T = sizeof(T)

function blength_helper(::Type{F}, veclens, j, s, v::Val{P}, T) where {P,F}
    if is_layout_variable(F, true)
        j += 1
        if j > length(veclens)
            throw(ArgumentError("not enough variable length specifiers for $T"))
        end
        vl = veclens[j]
    else
        vl = ()
    end
    al = Base.datatype_alignment(F)
    s = align(s, al)
    s += F <: LVarVector ? blength(F, vl, v) :
        P && F <: Union{Ptr,LForwardReference} ? blength(eltype(F), vl, v) : 0

    j, s
end

"""
    Cserialize(::Type{T}, source::Any)

Convert the julia object `source` into a byte vector to be used in C.
The process is controlled by the layout type recursively.

The resulting vector contains only data described in `T`.
The field, vector element or bit data required by `T` are taken from `source`
if available in a corresponding part. Other data are filled with 0-bytes.

If `T` is a structure, corresponding fields in source are by the same name.
If `T` is a vector, corresponding elements in source are by the same index.
If `T` is a `Ptr{S}`, the space for a pointer is reserved and filled with
the offset integer (of same size), while the `S` object is appended at the end
of the serialization stream.

Finally all offset integers are replaced by actual pointers.
"""
function Cserialize(::Type{T}, src) where T
    buf = UInt8[]
    rea = Int[]
    off = Cserialize!(T, src, buf, 0, rea)
    resize!(buf, off)
    relocate!(buf, rea)
end

function Cserialize!(::Type{T}, src, buf, off, rea::Vector{Int}) where T
    ctx = Tuple{Integer,Type,Any}[]
    noff = _Cserialize!(T, src, buf, off, rea, ctx)
    for (poff, F, src) in ctx
        noff = align(noff)
        ensure!(buf, poff, sizeof(Ptr))
        p = pointer(buf) + poff
        q = noff
        Base.unsafe_store!(Ptr{Ptr{F}}(p), Ptr{F}(q))
        noff = Cserialize!(F, src, buf, noff, rea)
    end
    align(noff)
end

function _Cserialize!(::Type{T}, src::S, buf::Vector{UInt8}, off::Integer, rea, ctx) where {T<:Layout,S}
    ensure!(buf, off, sizeof(T))
    noff = off
    for i = 1:fieldcount(T)
        F = fieldtype(T, i)
        f = fieldname(T, i)
        x = fieldoffset(T, i)
        if hasproperty(src, f)
            noff = _Cserialize!(F, getproperty(src, f), buf, off + x, rea, ctx)
        elseif F <: LVarVector
            throw(ArgumentError("need src vector $f to determine length"))
        else
            noff = off + x + Base.aligned_sizeof(F)
        end
    end
    noff
end
function _Cserialize!(::Type{<:LFixedVector{T,N}}, src::AbstractVector, buf::Vector{UInt8}, off::Integer, rea, ctx) where {T,N}
    as = Base.aligned_sizeof(T)
    ensure!(buf, off, as * N)
    n = length(src)
    for i = 1:min(N, n)
        off = _Cserialize!(T, src[i], buf, off, rea, ctx)
        off = align(off, T)
    end
    if N > n
        off += as * (N - n)
    end
    off
end
function _Cserialize!(::Type{S}, src::AbstractVector, buf::Vector{UInt8}, off::Integer, rea, ctx) where {T,F,S<:LVarVector{T,F}}
    for i = 1:length(src)
        off = _Cserialize!(T, src[i], buf, off, rea, ctx)
        off = align(off, T)
    end
    off
end

function _Cserialize!(::Type{Cstring}, src::Nothing, buf::Vector{UInt8}, off::Integer, rea, ctx)
    alignptr(off)
end
function _Cserialize!(::Type{Cstring}, src, buf::Vector{UInt8}, off::Integer, rea, ctx)
    s = string(src)
    n = length(s) + 1
    v = codeunits(s)
    pushall!(rea, ctx, off, LFixedVector{UInt8,n}, v)
    alignptr(off)
end

function _Cserialize!(::Type{T}, src, buf::Vector{UInt8}, off::Integer, rea, ctx) where T
    if isbitstype(T)
        s = Base.aligned_sizeof(T)
        ensure!(buf, off, s)
        p = pointer(buf) + off
        Base.unsafe_store!(Ptr{T}(p), convert(T, src))
        off + s
    else
        throw(ArgumentError("cannot serialize type $T"))
    end
end

function _Cserialize!(::Type{Ptr{T}}, src::Nothing, buf::Vector{UInt8}, off::Integer, rea, ctx) where T
    alignptr(off)
end
function _Cserialize!(::Type{Ptr{T}}, src, buf::Vector{UInt8}, off::Integer, rea, ctx) where T
    pushall!(rea, ctx, off, T, src)
    alignptr(off)
end

alignptr(off) = align(off + sizeof(Ptr))

"""
    pushall!(relocs::Vector{Int}, ctx::Vector{Tuple}, offset, Type, value}

push! the `offset` to `relocs` and the tuple `(offset, Type, value)` in `ctx`.
The `relocs` are finally used to replace offset values by pointers.
The `ctx` is used push back processing for later serializing.
"""
function pushall!(rea::Vector{<:Integer}, ctx::Vector{<:Tuple}, off, T, v)
    push!(rea, off)
    push!(ctx, (off, T, v))
end

"""
    ensure!(buf::Vector, off, size)

Ensure that the size of `buf` is at least `off + size` by maybe resizing `buf`.
Added space is filled with zero bytes.
"""
function ensure!(buf::Vector{UInt8}, off::Integer, siz::Integer)
    n = sizeof(buf)
    m = off + siz
    if n < m
        resize!(buf, m)
        for i = n+1:m
            buf[i] = 0
        end
    end
    buf
end

"""
    relocate!(buffer::Vector, offsets)

In vector `buffer`, at the byte offsets stored in `offsets`, offset values (into buffer) are
stored as `Int` values. The are replaced by `Ptr` values into the data area of `buffer`.

It is essential, that the data area is not changed after this process, that means no
`resize!`, `push!`, etc. are allowed after this final fix of pointer values to be used in
C-calls.
"""
function relocate!(buf::AbstractVector, rea::AbstractVector{<:Integer})
    p0 = pointer(buf)
    for off in rea
        p = p0 + off
        q = unsafe_load(Ptr{UInt}(p)) + p0
        unsafe_store!(Ptr{Ptr{UInt8}}(p), q)
    end
    buf
end


text:
using Query
using DataFrames

df = DataFrame(name=["John", "Sally", "Kirk"], age=[23., 42., 59.], children=[3,5,2])

x = @from i in df begin
    @let count = length(i.name)
    @let kids_per_year = i.children / i.age
    @where count > 4
    @select {Name=i.name, Count=count, KidsPerYear=kids_per_year}
    @collect DataFrame
end

println(x)


text:
export residual, residual!, jac_residual, jac_structure_residual, jac_structure_residual!
export jac_coord_residual!, jac_coord_residual, jprod_residual, jprod_residual!
export jtprod_residual, jtprod_residual!, jac_op_residual, jac_op_residual!
export hess_residual, hess_structure_residual, hess_structure_residual!
export hess_coord_residual!, hess_coord_residual, jth_hess_residual
export hprod_residual, hprod_residual!, hess_op_residual, hess_op_residual!

mutable struct NLSProjector{T, S} <: AbstractNLSModel{T, S}
  model::AbstractNLSModel{T, S}
  meta::NLPModelMeta{T, S}
  nls_meta::NLSMeta{T, S}
  counters::NLSCounters
  d::S
  function NLSProjector(
    model::AbstractNLSModel{T, S},
    d::S;
    name = "Projection over $(model.meta.name)",
  ) where {T, S}
    nvar = length(d)
    x0 = d #fill!(S(undef, nvar), zero(T))
    meta = NLPModelMeta{T, S}(nvar, x0 = x0, name = name)
    nls_meta = NLSMeta{T, S}(nvar, nvar, nnzj = nvar, nnzh = 0)
    return new{T, S}(model, meta, nls_meta, NLSCounters(), d)
  end
end

function NLPModels.residual!(model::NLSProjector, x, Fx)
  Fx .= x - model.d
  return Fx
end

function NLPModels.jac_structure_residual!(
  model::NLSProjector,
  rows::AbstractVector{<:Integer},
  cols::AbstractVector{<:Integer},
)
  @lencheck model.nls_meta.nnzj rows cols
  for i = 1:(model.meta.nvar)
    rows[i] = i
    cols[i] = i
  end
  return rows, cols
end

function NLPModels.jac_coord_residual!(
  model::NLSProjector{T, S},
  x::AbstractVector,
  vals::AbstractVector,
) where {T, S}
  @lencheck model.meta.nvar x
  @lencheck model.nls_meta.nnzj vals
  increment!(model, :neval_jac_residual)
  vals .= one(T)
  return vals
end

function NLPModels.jprod_residual!(
  model::NLSProjector,
  x::AbstractVector,
  v::AbstractVector,
  Jv::AbstractVector,
)
  @lencheck model.meta.nvar x v
  @lencheck model.nls_meta.nequ Jv
  increment!(model, :neval_jprod_residual)
  Jv .= v
  return Jv
end

function NLPModels.jtprod_residual!(
  model::NLSProjector,
  x::AbstractVector,
  v::AbstractVector,
  Jtv::AbstractVector,
)
  @lencheck model.meta.nvar x Jtv
  @lencheck model.nls_meta.nequ v
  increment!(model, :neval_jtprod_residual)
  Jtv .= v
  return Jtv
end

function NLPModels.hess_structure_residual!(::NLSProjector, rows, cols)
  return rows, cols
end

function NLPModels.hess_coord_residual!(::NLSProjector, x, v, vals)
  return vals
end

function NLPModels.hprod_residual!(::NLSProjector, x, i, v, Hiv)
  return Hiv
end

for meth in [
  :cons!,
  :jth_congrad!,
  :jth_sparse_congrad,
  :jac_structure!,
  :jac_coord!,
  :jprod!,
  :jtprod!,
  :jth_hess_coord!,
  :jth_hprod!,
  :ghjvprod!,
  :hess_structure!,
  :hess_coord!,
  :hprod!,
]
  premeth = Symbol("NLPModels.", meth)
  @eval begin
    $premeth(model::NLSProjector, args...; kwargs...) = $meth(model.model, args...; kwargs...)
  end
end


text:
module Resonators
using AxisArrays
using InstrumentControl
using InstrumentControl.VNA
using InstrumentControl.E5071C

import InstrumentControl.measure

export search_sidecoupled

function measure(x::VNA.FSweep)

    old_timeout = x.ins[Timeout]
    old_avgtrig = x.ins[AveragingTrigger]

    x.ins[Timeout] = 10000
    x.ins[NumTraces] = 3
    x.ins[VNA.Graphs] = [1 2; 1 3]
    x.ins[VNA.Format] = :PolarComplex
    x.ins[VNA.Format, 1, 2] = :LogMagnitude
    x.ins[VNA.Format, 1, 3] = :Phase
    x.ins[VNA.Parameter, 1, 2] = :S21
    x.ins[VNA.Parameter, 1, 3] = :S21
    x.ins[VNA.Parameter] = :S21
    x.ins[TriggerSource] = :BusTrigger

    if x.reject > 0
        old_avg = x.ins[Averaging]
        x.ins[Averaging] = false
        for i in 1:x.reject
            trig1(x.ins)
            sleep(sweeptime(x.ins))
            opc(x.ins)
        end
        x.ins[Averaging] = old_avg
    end

    x.ins[Averaging] && (x.ins[AveragingTrigger] = true)

    trig1(x.ins)
    sleep(sweeptime(x.ins))
    opc(x.ins)

    npts = x.ins[NumPoints]
    sparams = [:S11, :S21]
    result = AxisArray(Array{Complex{Float64}}(npts, length(sparams)),
        Axis{:f}(stimdata(x.ins)), Axis{:sparam}(sparams))
    for s in sparams
        x.ins[VNA.Parameter] = s
        result[Axis{:sparam}(s)] = VNA.data(x.ins, :PolarComplex)
    end
    result = transpose(result)

    autoscale(x.ins,1,1)
    autoscale(x.ins,1,2)
    autoscale(x.ins,1,3)

    # respect previous settings
    x.ins[Timeout] = old_timeout
    x.ins[AveragingTrigger] = old_avgtrig

    result
end

function search_sidecoupled(ins::InstrumentVNA, startf::Real, stopf::Real, cutoff=3, span=5e5)
    startf > stopf && error("Start frequency > stop frequency.")

    # Calculate frequency steps
    stepf = stopf/span
    numpts = length(startf:stepf:stopf)

    ins[TriggerSource] = :BusTrigger
    ins[FrequencyStart] = startf
    ins[FrequencyStop] = stopf
    ins[NumPoints] = numpts
    ins[Averaging] = false

    ins[VNA.Graphs] = [1]
    ins[NumTraces] = 1
    ins[VNA.Format, 1, 1] = :LogMagnitude
    ins[VNA.Parameter] = :S21
    ti = ins[SweepTime]

    trig1(ins)
    sleep(ti)
    opc(ins)

    # shotgun(ins, markers)
    markers = 1:nummarkers(ins)
    for m in markers
        ins[VNA.Marker, m] = true
        ins[VNA.MarkerX, m] = startf
    end

    peaks = Float64[]
    for (i,mi) in enumerate(markers)
        println(i)
        mx = search(ins, MarkerSearch(:RightPeak, 1, 1, mi, cutoff, VNA.Negative()))
        println(mx)
        isnan(mx) && break
        push!(peaks, mx)
        for mj in (i+1):nummarkers(ins)
            ins[VNA.MarkerX, mj] = mx+span
        end
    end

    for m in (length(peaks)+1):nummarkers(ins)
        ins[VNA.Marker, m] = false
    end

    # isnan(mx[1]) && (ins[VNA.Marker, 1] = false)
    # for m in 2:length(markers)
    #     (mx[m] == mx[m-1] || isnan(mx[m])) && (ins[VNA.Marker, m] = false)
    # end

    # fs = Float64[]
    # for m in markers
    #     if ins[VNA.Marker, m] == true
    #         count += 1
    #         push!(fs, mx[m])
    #     end
    # end

    # if count == 0
    if length(peaks) == 0
        warn("No peaks were found.")
    end

    autoscale(ins, 1, 1)
    peaks

    # count == 2 && configure(ins, Graphs, [1 2])
    # 3 <= count <= 4 && configure(ins, Graphs, [1 2; 3 4])
    # 5 <= count <= 6 && configure(ins, Graphs, [1 2 3; 4 5 6])
    # 7 <= count <= 8 && configure(ins, Graphs, [1 2 3 4; 5 6 7 8])
    # count == 9 && configure(ins, Graphs, [1 2 3; 4 5 6; 7 8 9])
    # 10 <= count <= 12 && configure(ins, Graphs, [1 2 3 4; 5 6 7 8; 9 10 11 12])
    # 13 <= count && configure(ins, Graphs, [1 2 3 4; 5 6 7 8; 9 10 11 12])
end

end # module


text:
using Test

using Unitful, UnitfulAtomic

using EquationsOfState
using EquationsOfState.Collections
using EquationsOfState.NonlinearFitting

@testset "Test fitting energy with different element types" begin
    result = BirchMurnaghan3rd(
        0.0057009512119028044,
        103.58772269057364,
        -144.45152457521132,
        -40.31992619868024,
    ) |> Collections.fieldvalues
    @test isapprox(
        lsqfit(
            EnergyForm(),
            BirchMurnaghan3rd(1, 2, 3.0, 0),
            [1, 2, 3, 4, 5],
            [5, 6, 9, 8, 7],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-5,
    )
    @test isapprox(
        lsqfit(
            EnergyForm(),
            BirchMurnaghan3rd(1, 2, 3, 0),
            [1, 2, 3, 4, 5.0],
            [5, 6, 9, 8, 7],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-5,
    )
    @test isapprox(
        lsqfit(
            EnergyForm(),
            BirchMurnaghan3rd(1, 2, 3.0, 0),
            [1, 2, 3, 4, 5],
            [5, 6, 9, 8, 7.0],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-5,
    )
    @test isapprox(
        lsqfit(
            EnergyForm(),
            BirchMurnaghan3rd(1, 2, 3, 0),
            [1, 2, 3, 4, 5],
            [5, 6, 9, 8, 7],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-5,
    )
    @test isapprox(
        lsqfit(
            EnergyForm(),
            BirchMurnaghan3rd(big(1), 2, big(3), 0),
            [1, 2, 3, 4, 5],
            [5, 6, 9, 8, 7],
        ) |> Collections.fieldvalues,
        big.(result);
        atol = 1e-4,
    )
    @test isapprox(
        lsqfit(
            EnergyForm(),
            BirchMurnaghan3rd(big(1), 2, 3, 0),
            BigFloat[1, 2, 3, 4, 5],
            [5, 6, 9, 8, 7],
        ) |> Collections.fieldvalues,
        big.(result);
        atol = 1e-4,
    )
    @test isapprox(
        lsqfit(
            EnergyForm(),
            BirchMurnaghan3rd(big(1.0), 2, 3, 0),
            [1, 2, 3, 4, 5],
            BigInt[5, 6, 9, 8, 7],
        ) |> Collections.fieldvalues,
        big.(result);
        atol = 1e-4,
    )
end

@testset "Test fitting pressure with different element types" begin
    result = BirchMurnaghan3rd(
        1.1024687826597717,
        29.30861698140365,
        12.689089871112746,
        0.0,
    ) |> Collections.fieldvalues
    @test isapprox(
        lsqfit(
            PressureForm(),
            BirchMurnaghan3rd(1, 2, 3.0, 0),
            [1, 2, 3, 4, 5],
            [5, 6, 9, 8, 7],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-6,
    )
    @test isapprox(
        lsqfit(
            PressureForm(),
            BirchMurnaghan3rd(1, 2, 3, 0),
            [1, 2, 3, 4, 5.0],
            [5, 6, 9, 8, 7],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-6,
    )
    @test isapprox(
        lsqfit(
            PressureForm(),
            BirchMurnaghan3rd(1, 2, 3.0, 0),
            [1, 2, 3, 4, 5],
            [5, 6, 9, 8, 7.0],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-6,
    )
    @test isapprox(
        lsqfit(
            PressureForm(),
            BirchMurnaghan3rd(1, 2, 3, 0),
            [1, 2, 3, 4, 5],
            [5, 6, 9, 8, 7],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-6,
    )
    @test isapprox(
        lsqfit(
            PressureForm(),
            BirchMurnaghan3rd(1, 2, big(3), 0),
            [1, 2, 3, 4, 5],
            BigInt[5, 6, 9, 8, 7],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-6,
    )
    @test isapprox(
        lsqfit(
            PressureForm(),
            BirchMurnaghan3rd(1, 2, big(3.0), 0),
            [1, 2, 3, 4, 5],
            [5, 6, 9, 8, 7],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-6,
    )
    @test isapprox(
        lsqfit(
            PressureForm(),
            BirchMurnaghan3rd(1, 2, big(3), 0),
            [big(1), 2, 3, 4, 5],
            [5, 6, 9, 8, 7],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-6,
    )
end

@testset "Test fitting bulk modulus with different element types" begin
    result = BirchMurnaghan3rd(
        7.218928431312577,
        5.007900469653902,
        4.06037725509478,
        0.0,
    ) |> Collections.fieldvalues
    @test isapprox(
        lsqfit(
            BulkModulusForm(),
            BirchMurnaghan3rd(1, 2, 3.0, 0),
            [1, 2, 3, 4, 5],
            [5, 6, 9, 8, 7],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-5,
    )
    @test isapprox(
        lsqfit(
            BulkModulusForm(),
            BirchMurnaghan3rd(1, 2, 3, 0),
            [1, 2, 3, 4, 5.0],
            [5, 6, 9, 8, 7],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-5,
    )
    @test isapprox(
        lsqfit(
            BulkModulusForm(),
            BirchMurnaghan3rd(1, 2, 3.0, 0),
            [1, 2, 3, 4, 5],
            [5, 6, 9, 8, 7.0],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-5,
    )
    @test isapprox(
        lsqfit(
            BulkModulusForm(),
            BirchMurnaghan3rd(1, 2, 3, 0),
            [1, 2, 3, 4, 5],
            [5, 6, 9, 8, 7],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-5,
    )
    @test isapprox(
        lsqfit(
            BulkModulusForm(),
            BirchMurnaghan3rd(1, 2, 3, big(0)),
            [1, 2, 3, 4, 5],
            [5, 6, 9, 8, 7],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-5,
    )
    @test isapprox(
        lsqfit(
            BulkModulusForm(),
            BirchMurnaghan3rd(1, 2, 3, big(0)),
            [1, 2, 3, big(4.0), 5],
            [big(5), 6, 9, 8, 7.0],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-5,
    )
    @test isapprox(
        lsqfit(
            BulkModulusForm(),
            BirchMurnaghan3rd(1, 2, 3, big(0)),
            [1, 2, 3, 4, 5],
            [big(5), 6, 9, 8, 7.0],
        ) |> Collections.fieldvalues,
        result;
        atol = 1e-5,
    )
end

# Data in the following tests are from
# https://github.com/materialsproject/pymatgen/blob/1f0957b8525ddc7d12ea348a19caecebe6c7ff34/pymatgen/analysis/tests/test_eos.py
@testset "Test data from Pymatgen" begin
    volumes = [
        25.987454833,
        26.9045702104,
        27.8430241908,
        28.8029649591,
        29.7848370694,
        30.7887887064,
        31.814968055,
        32.8638196693,
        33.9353435494,
        35.0299842495,
        36.1477417695,
        37.2892088485,
        38.4543854865,
        39.6437162376,
        40.857201102,
        42.095136449,
        43.3579668329,
        44.6456922537,
        45.9587572656,
        47.2973100535,
        48.6614988019,
        50.0517680652,
        51.4682660281,
        52.9112890601,
        54.3808371612,
        55.8775030703,
        57.4014349722,
        58.9526328669,
    ]
    energies = [
        -7.63622156576,
        -8.16831294894,
        -8.63871612686,
        -9.05181213218,
        -9.41170988374,
        -9.72238224345,
        -9.98744832526,
        -10.210309552,
        -10.3943401353,
        -10.5427238068,
        -10.6584266073,
        -10.7442240979,
        -10.8027285713,
        -10.8363890521,
        -10.8474912964,
        -10.838157792,
        -10.8103477586,
        -10.7659387815,
        -10.7066179666,
        -10.6339907853,
        -10.5495538639,
        -10.4546677714,
        -10.3506386542,
        -10.2386366017,
        -10.1197772808,
        -9.99504030111,
        -9.86535084973,
        -9.73155247952,
    ]
    @test isapprox(
        lsqfit(EnergyForm(), BirchMurnaghan3rd(40, 0.5, 4, 0), volumes, energies) |>
        Collections.fieldvalues,
        BirchMurnaghan3rd(
            40.98926572528106,
            0.5369258245417454,
            4.178644235500821,
            -10.842803908240892,
        ) |> Collections.fieldvalues,
    )
    @test isapprox(
        lsqfit(EnergyForm(), Murnaghan(41, 0.5, 4, 0), volumes, energies) |>
        Collections.fieldvalues,
        Murnaghan(
            41.13757930387086,
            0.5144967693786603,
            3.9123862262572264,
            -10.836794514626673,
        ) |> Collections.fieldvalues,
    )
    @test isapprox(
        lsqfit(EnergyForm(), PoirierTarantola3rd(41, 0.5, 4, 0), volumes, energies) |>
        Collections.fieldvalues,
        PoirierTarantola3rd(
            40.86770643373908,
            0.5667729960804602,
            4.331688936974368,
            -10.851486685041658,
        ) |> Collections.fieldvalues,
    )
    @test isapprox(
        lsqfit(EnergyForm(), Vinet(41, 0.5, 4, 0), volumes, energies) |>
        Collections.fieldvalues,
        Vinet(
            40.916875663779784,
            0.5493839425156859,
            4.3051929654936885,
            -10.846160810560756,
        ) |> Collections.fieldvalues,
    )
    # 'deltafactor': {'b0': 0.5369258245611414,
#             'b1': 4.178644231924639,
#             'e0': -10.842803908299294,
#             'v0': 40.989265727927936},
# 'numerical_eos': {'b0': 0.5557257614101998,
#             'b1': 4.344039148405489,
#             'e0': -10.847490826530702,
#             'v0': 40.857200064982536},
# }
end

@testset "Test Mg dataset" begin
    mp153_volumes = [
        16.69182365,
        17.25441763,
        17.82951915,
        30.47573817,
        18.41725977,
        29.65211363,
        28.84346369,
        19.01777055,
        28.04965916,
        19.63120886,
        27.27053682,
        26.5059864,
        20.25769112,
        25.75586879,
        20.89736201,
        25.02003097,
        21.55035204,
        24.29834347,
        22.21681221,
        23.59066888,
        22.89687316,
    ]

    mp153_energies = [
        -1.269884575,
        -1.339411225,
        -1.39879471,
        -1.424480995,
        -1.44884184,
        -1.45297499,
        -1.4796246,
        -1.49033594,
        -1.504198485,
        -1.52397006,
        -1.5264432,
        -1.54609291,
        -1.550269435,
        -1.56284009,
        -1.569937375,
        -1.576420935,
        -1.583470925,
        -1.58647189,
        -1.591436505,
        -1.592563495,
        -1.594347355,
    ]

    mp153_known_energies_vinet = [
        -1.270038831,
        -1.339366487,
        -1.398683238,
        -1.424556061,
        -1.448746649,
        -1.453000456,
        -1.479614511,
        -1.490266797,
        -1.504163502,
        -1.523910268,
        -1.526395734,
        -1.546038792,
        -1.550298657,
        -1.562800797,
        -1.570015274,
        -1.576368392,
        -1.583605186,
        -1.586404575,
        -1.591578378,
        -1.592547954,
        -1.594410995,
    ]

    fitted_eos = lsqfit(EnergyForm(), Vinet(23, 0.5, 4, -2), mp153_volumes, mp153_energies)
    @test isapprox(
        fitted_eos |> Collections.fieldvalues,
        Vinet(
            22.95764559358769,
            0.2257091141420788,
            4.060543387224629,
            -1.5944292606251582,
        ) |> Collections.fieldvalues,
    )
    @test isapprox(
        map(apply(EnergyForm(), fitted_eos), mp153_volumes),
        mp153_known_energies_vinet;
        atol = 1e-5,
    )
end

@testset "Test Si dataset" begin
    mp149_volumes = [
        15.40611854,
        14.90378698,
        16.44439516,
        21.0636307,
        17.52829835,
        16.98058208,
        18.08767363,
        18.65882487,
        19.83693435,
        15.91961152,
        22.33987173,
        21.69548924,
        22.99688883,
        23.66666322,
        20.44414922,
        25.75374305,
        19.24187473,
        24.34931029,
        25.04496106,
        27.21116571,
        26.4757653,
    ]

    mp149_energies = [
        -4.866909695,
        -4.7120965,
        -5.10921253,
        -5.42036228,
        -5.27448405,
        -5.200810795,
        -5.331915665,
        -5.3744186,
        -5.420058145,
        -4.99862686,
        -5.3836163,
        -5.40610838,
        -5.353700425,
        -5.31714654,
        -5.425263555,
        -5.174988295,
        -5.403353105,
        -5.27481447,
        -5.227210275,
        -5.058992615,
        -5.118805775,
    ]

    mp149_known_energies_vinet = [
        -4.866834585,
        -4.711786499,
        -5.109642598,
        -5.420093739,
        -5.274605844,
        -5.201025714,
        -5.331899365,
        -5.374315789,
        -5.419671568,
        -4.998827503,
        -5.383703409,
        -5.406038887,
        -5.353926272,
        -5.317484252,
        -5.424963418,
        -5.175090887,
        -5.403166824,
        -5.275096644,
        -5.227427635,
        -5.058639193,
        -5.118654229,
    ]

    fitted_eos = lsqfit(EnergyForm(), Vinet(20, 0.5, 4, -5), mp149_volumes, mp149_energies)
    @test isapprox(
        fitted_eos |> Collections.fieldvalues,
        Vinet(
            20.446696754873944,
            0.5516638521306302,
            4.324373909783161,
            -5.424963389876503,
        ) |> Collections.fieldvalues,
    )
    @test isapprox(
        map(apply(EnergyForm(), fitted_eos), mp149_volumes),
        mp149_known_energies_vinet;
        atol = 1e-5,
    )
end

@testset "Test Ti dataset" begin
    mp72_volumes = [
        12.49233296,
        12.91339188,
        13.34380224,
        22.80836212,
        22.19195533,
        13.78367177,
        21.58675559,
        14.23310328,
        20.99266009,
        20.4095592,
        14.69220297,
        19.83736385,
        15.16106697,
        19.2759643,
        15.63980711,
        18.72525771,
        16.12851491,
        18.18514127,
        16.62729878,
        17.65550599,
        17.13626153,
    ]

    mp72_energies = [
        -7.189983803,
        -7.33985647,
        -7.468745423,
        -7.47892835,
        -7.54945107,
        -7.578012237,
        -7.61513166,
        -7.66891898,
        -7.67549721,
        -7.73000681,
        -7.74290386,
        -7.77803379,
        -7.801246383,
        -7.818964483,
        -7.84488189,
        -7.85211192,
        -7.87486651,
        -7.876767777,
        -7.892161533,
        -7.892199957,
        -7.897605303,
    ]

    mp72_known_energies_vinet = [
        -7.189911138,
        -7.339810181,
        -7.468716095,
        -7.478678021,
        -7.549402394,
        -7.578034391,
        -7.615240977,
        -7.669091347,
        -7.675683891,
        -7.730188653,
        -7.74314028,
        -7.778175824,
        -7.801363213,
        -7.819030923,
        -7.844878053,
        -7.852099741,
        -7.874737806,
        -7.876686864,
        -7.891937429,
        -7.892053535,
        -7.897414664,
    ]

    fitted_eos = lsqfit(EnergyForm(), Vinet(17, 0.5, 4, -7), mp72_volumes, mp72_energies)
    @test isapprox(
        fitted_eos |> Collections.fieldvalues,
        Vinet(
            17.13223026131245,
            0.7029766224730147,
            3.6388077563621812,
            -7.897414959124461,
        ) |> Collections.fieldvalues,
    )
    @test isapprox(
        map(apply(EnergyForm(), fitted_eos), mp72_volumes),
        mp72_known_energies_vinet;
        atol = 1e-5,
    )
end

@testset "`Test w2k-lda-na.dat` from `Gibbs2`" begin
    data = [
        159.9086 -323.4078898
        162.5738 -323.4089153
        165.2389 -323.4098546
        167.9041 -323.410722
        170.5692 -323.4115195
        173.2344 -323.4122481
        175.8995 -323.4129189
        178.5647 -323.413528
        181.2298 -323.4140871
        183.8949 -323.4145889
        186.5601 -323.4150471
        189.2252 -323.415459
        191.8904 -323.4158302
        194.5555 -323.4161579
        197.2207 -323.4164498
        199.8858 -323.4167071
        202.551 -323.4169305
        205.2161 -323.4171194
        207.8812 -323.4172809
        210.5464 -323.4174144
        213.2115 -323.4175216
        215.8767 -323.4176029
        218.5418 -323.417661
        221.207 -323.4176975
        223.8721 -323.41771
        226.5373 -323.4177051
        229.2024 -323.417682
        231.8675 -323.4176375
        234.5327 -323.417579
        237.1978 -323.4175048
        239.863 -323.4174142
        242.5281 -323.4173101
        245.1933 -323.4171922
        247.8584 -323.4170611
        250.5236 -323.4169184
        253.1887 -323.4167647
        255.8538 -323.4166002
        258.519 -323.4164244
        261.1841 -323.4162386
        263.8493 -323.4160446
        266.5144 -323.4158421
        269.1796 -323.4156312
        271.8447 -323.4154125
        274.5098 -323.4151861
        277.175 -323.4149528
        279.8401 -323.4147131
        282.5053 -323.414467
        285.1704 -323.414215
        287.8356 -323.4139583
        290.5007 -323.4136953
        293.1659 -323.4134285
        295.831 -323.4131559
        298.4961 -323.4128797
        301.1613 -323.4125984
        303.8264 -323.4123147
        306.4916 -323.4120269
        309.1567 -323.411736
        311.8219 -323.4114399
        314.487 -323.4111421
        317.1522 -323.4108418
        319.8173 -323.4105393
    ]

    @testset "without unit" begin
        volumes = data[:, 1]
        energies = data[:, 2]
        @test isapprox(
            lsqfit(EnergyForm(), Murnaghan(224, 0.006, 4, -323), volumes, energies) |>
            Collections.fieldvalues,
            Murnaghan(224.501825, 0.00060479524074699499, 3.723835, -323.417686) |>
            Collections.fieldvalues;
            atol = 1e-5,
        )
        @test isapprox(
            lsqfit(
                EnergyForm(),
                BirchMurnaghan2nd(224, 0.0006, -323),
                volumes,
                energies,
            ) |> Collections.fieldvalues,
            BirchMurnaghan2nd(
                223.7192539523166,
                0.0006268341030294977,
                -323.4177121144877,
            ) |> Collections.fieldvalues;
            atol = 1e-3,
        )
        @test lsqfit(
            EnergyForm(),
            BirchMurnaghan3rd(224, 0.0006, 4, -323),
            volumes,
            energies,
        ) |> Collections.fieldvalues ≈ BirchMurnaghan3rd(
            224.444565,
            0.00062506191050572675,
            3.740369,
            -323.417714,
        ) |> Collections.fieldvalues
        @test isapprox(
            lsqfit(
                EnergyForm(),
                BirchMurnaghan4th(224, 0.0006, 4, -5460, -323),
                volumes,
                energies,
            ) |> Collections.fieldvalues,
            BirchMurnaghan4th(
                224.45756238103118,
                0.0006229382380931005,
                3.730991532958105,
                -5322.696706065215,
                -323.4177113158582,
            ) |> Collections.fieldvalues;
            rtol = 1e-6,
        )
        @test lsqfit(EnergyForm(), Vinet(224, 0.0006, 4, -323), volumes, energies) |>
              Collections.fieldvalues ≈ Vinet(
            224.45278665796354,
            0.0006313500637481759,
            3.7312381477678853,
            -323.4177229576912,
        ) |> Collections.fieldvalues
        @test isapprox(
            lsqfit(
                EnergyForm(),
                PoirierTarantola3rd(100, 0.0006, 3.7, -323),
                volumes,
                energies,
            ) |> Collections.fieldvalues,
            PoirierTarantola3rd(224.509208, 0.000635892264159838, 3.690448, -323.41773) |>
            Collections.fieldvalues;
            atol = 1e-5,
        )
        # @test lsqfit(EnergyForm(), PoirierTarantola4th(220, 0.0006, 3.7, -5500, -323), volumes, energies; lower = Float64[220, 0, 3, -6000, -400], upper = Float64[300, 0.01, 5, -5000, -300]) ≈ PoirierTarantola4th(224.430182, 0.0006232241765069493, 3.758360, -5493.859729817176, -323.417712)
    end # testset

    @testset "with units" begin
        volumes = data[:, 1] .* u"bohr^3"
        energies = data[:, 2] .* u"Ry"
        @test isapprox(
            ustrip.(lsqfit(
                EnergyForm(),
                Murnaghan(
                    224.445371 * u"bohr^3",
                    9.164446 * u"GPa",
                    3.752432,
                    -161.708856 * u"hartree",
                ),
                volumes,
                energies,
            ) |> Collections.fieldvalues),
            ustrip.(BirchMurnaghan3rd(
                224.501825 * u"bohr^3",
                8.896845 * u"GPa",
                3.723835,
                -161.708843 * u"hartree",
            ) |> Collections.fieldvalues);
            atol = 1e-5,
        )
        @test ustrip.(lsqfit(
            EnergyForm(),
            BirchMurnaghan3rd(
                224.445371 * u"bohr^3",
                9.164446 * u"GPa",
                3.752432,
                -161.708856 * u"hartree",
            ),
            volumes,
            energies,
        ) |> Collections.fieldvalues) ≈ ustrip.(BirchMurnaghan3rd(
            224.444565 * u"bohr^3",
            9.194978 * u"GPa",
            3.740369,
            -161.708857 * u"hartree",
        ) |> Collections.fieldvalues)
        # @test ustrip.(
        #     lsqfit(EnergyForm(), BirchMurnaghan4th(224.445371u"bohr^3", 9.164446u"GPa", 3.752432, -0.371174u"1/GPa", -161.708856u"hartree"), volumes, energies) |> Collections.fieldvalues
        # ) ≈
        # ustrip.(
        #     BirchMurnaghan4th(224.457562u"bohr^3", 9.163736u"GPa", 3.730992, -0.361830u"1/GPa", -161.708856u"hartree") |> Collections.fieldvalues
        # )
        # Non-linear fitting: BM4: 4th order Birch-Murnaghan EOS
        # Parameters (5) start / converged
        # E0 (Hy)               -161.708856      -161.708856
        # V0 (bohr^3)            224.445371       224.457562
        # B0 (GPa)                 9.164446         9.163736
        # B1p                      3.752432         3.730992
        # B2p (1/GPa)             -0.371174        -0.361830
    end
end

@testset "`Test w2k-lda-k.dat` from `Gibbs2`" begin
    data = [
        243.6753 -1201.1776527
        246.941 -1201.1792783
        250.2358 -1201.180845
        253.5598 -1201.1823535
        256.9131 -1201.1838048
        260.2958 -1201.1857092
        263.7081 -1201.186542
        267.1501 -1201.1878269
        270.6219 -1201.1890616
        274.1236 -1201.1906847
        277.6555 -1201.1913795
        281.2175 -1201.1924633
        284.8099 -1201.1935021
        288.4327 -1201.194493
        292.0861 -1201.1954413
        295.7703 -1201.1963445
        299.4853 -1201.1975477
        303.2312 -1201.1980347
        307.0083 -1201.1988083
        310.8166 -1201.1995541
        314.6563 -1201.2002575
        318.5275 -1201.2009258
        322.4303 -1201.2015585
        326.3648 -1201.2021556
        330.3313 -1201.2028674
        334.3297 -1201.2032521
        338.3603 -1201.2037525
        342.4231 -1201.2042232
        346.5183 -1201.2046643
        350.6461 -1201.2050747
        354.8065 -1201.2054579
        358.9997 -1201.2058132
        363.2257 -1201.206143
        367.4849 -1201.2064473
        371.7771 -1201.2067259
        376.1027 -1201.20698
        380.4617 -1201.2072114
        384.8543 -1201.2074353
        389.2805 -1201.207606
        393.7406 -1201.2077761
        398.2345 -1201.2079159
        402.7626 -1201.2080394
        407.3248 -1201.2081438
        411.9214 -1201.2082292
        416.5524 -1201.2082962
        421.218 -1201.2083459
        425.9183 -1201.2083783
        430.6534 -1201.2083929
        435.4235 -1201.2084741
        440.2287 -1201.2083688
        445.0692 -1201.2083404
        449.9449 -1201.2082935
        454.8562 -1201.208232
        459.803 -1201.2081564
        464.7856 -1201.2081173
        469.8041 -1201.2079657
        474.8586 -1201.2078516
        479.9492 -1201.2077239
        485.076 -1201.2075854
        490.2392 -1201.2074633
        495.4389 -1201.207274
        500.6753 -1201.207102
        505.9484 -1201.2069354
        511.2585 -1201.2067294
        516.6055 -1201.2065235
        521.9897 -1201.2063119
        527.4112 -1201.2060909
        532.8701 -1201.2058612
        538.3665 -1201.2056224
        543.9006 -1201.2053821
        549.4725 -1201.2051214
        555.0823 -1201.2048592
        560.7302 -1201.2045886
        566.4163 -1201.2043109
        572.1406 -1201.2039935
        577.9034 -1201.2037374
        583.7048 -1201.2034399
        589.5448 -1201.2031357
        595.4237 -1201.2028264
        601.3416 -1201.2025196
        607.2985 -1201.2021901
        613.2946 -1201.2018636
        619.3301 -1201.201532
        625.405 -1201.2011948
        631.5196 -1201.200853
        637.6738 -1201.2005083
        643.868 -1201.2001563
        650.1021 -1201.199944
        656.3763 -1201.1994441
        662.6907 -1201.1990811
        669.0455 -1201.1987147
        675.4408 -1201.1983441
        681.8768 -1201.1978577
        688.3535 -1201.1975937
        694.871 -1201.1970935
        701.4296 -1201.1968292
        708.0294 -1201.1963133
        714.6703 -1201.1960557
        721.3527 -1201.1956647
        728.0767 -1201.1952707
    ]

    @testset "without unit" begin
        volumes = data[:, 1]  # unit: bohr^3
        energies = data[:, 2]  # unit: Rydberg
        @test isapprox(
            lsqfit(EnergyForm(), Murnaghan(224, 0.0006, 4, -323), volumes, energies) |>
            Collections.fieldvalues,
            Murnaghan(
                435.05782299050884,
                0.00028297159355249787,
                3.5705032675000785,
                -1201.2082739321822,
            ) |> Collections.fieldvalues;
            atol = 1e-3,
        )
        @test isapprox(
            lsqfit(
                EnergyForm(),
                BirchMurnaghan2nd(224, 0.0006, -323),
                volumes,
                energies,
            ) |> Collections.fieldvalues,
            BirchMurnaghan2nd(
                430.10027687726716,
                0.000302451215462375,
                -1201.2083221436026,
            ) |> Collections.fieldvalues;
            atol = 1e-3,
        )
        @test isapprox(
            lsqfit(
                EnergyForm(),
                BirchMurnaghan3rd(224, 0.0006, 4, -323),
                volumes,
                energies,
            ) |> Collections.fieldvalues,
            BirchMurnaghan3rd(
                432.67139080209046,
                0.00030508544859901674,
                3.7894868450211923,
                -1201.2083959943404,
            ) |> Collections.fieldvalues;
            atol = 1e-3,
        )
        @test isapprox(
            lsqfit(
                EnergyForm(),
                BirchMurnaghan4th(432, 0.0003, 3.8, -11773, -1201),
                volumes,
                energies,
            ) |> Collections.fieldvalues,
            BirchMurnaghan4th(
                432.8012195854224,
                0.0003041889904166284,
                3.774020919355492,
                -11773.192574765615,
                -1201.2083912308235,
            ) |> Collections.fieldvalues;
            rtol = 1e-5,
        )
        @test isapprox(
            lsqfit(EnergyForm(), Vinet(432, 0.0003, 3.8, -1201), volumes, energies) |>
            Collections.fieldvalues,
            Vinet(
                432.04609865398015,
                0.0003137631070690569,
                3.837666939407128,
                -1201.2084453225773,
            ) |> Collections.fieldvalues;
            atol = 1e-3,
        )
    end # testset

    @testset "with units" begin
        volumes = data[:, 1] .* u"bohr^3"
        energies = data[:, 2] .* u"Ry"
        fitted_eos = lsqfit(
            EnergyForm(),
            BirchMurnaghan3rd(224 * u"bohr^3", 0.0006 * u"Ry/bohr^3", 4, -323 * u"Ry"),
            volumes,
            energies,
        )
        @test ustrip.(fitted_eos |> Collections.fieldvalues) ≈ ustrip.(BirchMurnaghan3rd(
            432.6713907942206 * u"bohr^3",
            0.00030508544829126676 * u"Ry/bohr^3",
            3.789486849598267,
            -1201.208395994332 * u"Ry",
        ) |> Collections.fieldvalues)
        @test ustrip.(lsqfit(
            EnergyForm(),
            BirchMurnaghan3rd(224 * u"bohr^3", 10 * u"GPa", 3.75, -161 * u"hartree"),
            volumes,
            energies,
        ) |> Collections.fieldvalues) ≈ ustrip.(BirchMurnaghan3rd(
            432.671390388525 * u"bohr^3",
            4.487961877912739 * u"GPa",
            3.7894868798185493,
            -600.6041979971637 * u"hartree",
        ) |> Collections.fieldvalues)
    end
end

@testset "`Test mgo-sety1.dat` from `Gibbs2`" begin
    data = [
        71.51955637984702 -34.00775398
        72.0783029140629 -34.0190082
        72.63704944827879 -34.02992647
        73.19579598249467 -34.04052022
        73.75454251671056 -34.05079922
        74.31328905092643 -34.06077277
        74.87203558514231 -34.0704506
        75.4307821193582 -34.07984653
        75.98952865357408 -34.08896085
        76.54827518778997 -34.0978042
        77.10702172200585 -34.10638855
        77.66576825622174 -34.11471583
        78.22451479043761 -34.12279583
        78.7832613246535 -34.13063609
        79.34200785886938 -34.13824518
        79.90075439308526 -34.1456268
        80.45950092730115 -34.15279001
        81.01824746151703 -34.15973876
        81.57699399573292 -34.16647895
        82.1357405299488 -34.1730202
        82.69448706416469 -34.17936412
        83.25323359838056 -34.18551683
        83.81198013259645 -34.19148806
        84.37072666681233 -34.19727555
        84.92947320102822 -34.20288898
        85.4882197352441 -34.20833138
        86.04696626945999 -34.21360745
        86.60571280367587 -34.21872393
        87.16445933789174 -34.2236828
        87.72320587210763 -34.22849177
        88.28195240632351 -34.23315234
        88.8406989405394 -34.23766788
        89.39944547475528 -34.24204214
        89.95819200897117 -34.24627838
        90.51693854318705 -34.25038459
        91.07568507740294 -34.25435835
        91.63443161161882 -34.25820545
        92.19317814583471 -34.26193399
        92.75192468005058 -34.26553901
        93.31067121426646 -34.2690287
        93.86941774848235 -34.27240614
        94.42816428269823 -34.27566932
        94.98691081691412 -34.27882789
        95.54565735112999 -34.28188043
        96.10440388534587 -34.28483105
        96.66315041956176 -34.28768098
        97.22189695377764 -34.29043394
        97.78064348799353 -34.29309257
        98.33939002220941 -34.29565801
        98.8981365564253 -34.29813575
        99.45688309064118 -34.30052413
        100.01562962485707 -34.30282697
        100.57437615907295 -34.30504676
        101.13312269328884 -34.30718589
        101.69186922750471 -34.30924656
        102.2506157617206 -34.31123046
        102.80936229593648 -34.31313623
        103.36810883015237 -34.31496818
        103.92685536436825 -34.31673139
        104.48560189858412 -34.31842177
        105.0443484328 -34.32004585
        105.60309496701589 -34.32160156
        106.16184150123178 -34.32309155
        106.72058803544766 -34.32451965
        107.27933456966355 -34.32588527
        107.83808110387943 -34.32719027
        108.39682763809532 -34.32844001
        108.9555741723112 -34.32962806
        109.51432070652709 -34.3307603
        110.07306724074297 -34.33183713
        110.63181377495884 -34.33286456
        111.19056030917473 -34.33383534
        111.74930684339061 -34.33475519
        112.3080533776065 -34.33562709
        112.86679991182238 -34.33644925
        113.42554644603825 -34.33722484
        113.98429298025414 -34.3379529
        114.54303951447002 -34.33863669
        115.10178604868591 -34.33927612
        115.6605325829018 -34.33987151
        116.21927911711768 -34.34042433
        116.77802565133356 -34.3409375
        117.33677218554945 -34.34140944
        117.89551871976533 -34.3418419
        118.45426525398122 -34.34223621
        119.0130117881971 -34.34259266
        119.57175832241298 -34.34291588
        120.13050485662886 -34.34320277
        120.68925139084475 -34.34345231
        121.24799792506063 -34.34366646
        121.80674445927652 -34.34384815
        122.36549099349239 -34.34399664
        122.92423752770827 -34.3441149
        123.48298406192416 -34.34420115
        124.04173059614004 -34.34425594
        124.60047713035593 -34.34428195
        125.15922366457181 -34.34427761
        125.7179701987877 -34.34424807
        126.27671673300358 -34.34418868
        126.83546326721947 -34.34410162
        127.39420980143535 -34.34398884
        127.95295633565124 -34.34385192
        128.51170286986712 -34.34368896
        129.07044940408298 -34.3434993
        129.6291959382989 -34.34328585
        130.18794247251475 -34.34305126
        130.74668900673066 -34.34279123
        131.30543554094652 -34.34250948
        131.8641820751624 -34.34220541
        132.4229286093783 -34.34187944
        132.98167514359417 -34.34153303
        133.54042167781006 -34.34116878
        134.09916821202594 -34.34078072
        134.65791474624183 -34.34037385
        135.2166612804577 -34.33994687
        135.7754078146736 -34.33950111
        136.33415434888948 -34.33903726
        136.89290088310537 -34.33855599
        137.45164741732123 -34.33805699
        138.01039395153714 -34.33754037
        138.569140485753 -34.33700777
        139.1278870199689 -34.33645907
        139.68663355418477 -34.33589553
        140.24538008840068 -34.33531772
        140.80412662261654 -34.33472332
        141.36287315683242 -34.33411208
        141.9216196910483 -34.33348704
        142.4803662252642 -34.33284619
        143.03911275948008 -34.33219292
    ]

    @testset "without unit" begin
        volumes = data[:, 1]  # unit: bohr^3
        energies = data[:, 2]  # unit: Rydberg
        @test isapprox(
            lsqfit(EnergyForm(), Murnaghan(110, 0.01, 4, -34), volumes, energies) |>
            Collections.fieldvalues,
            Murnaghan(
                124.88539638285143,
                0.012047999390789954,
                3.505337379827799,
                -34.34458025509868,
            ) |> Collections.fieldvalues;
            atol = 1e-3,
        )
        @test isapprox(
            lsqfit(EnergyForm(), BirchMurnaghan2nd(124, 0.012, -34), volumes, energies) |>
            Collections.fieldvalues,
            BirchMurnaghan2nd(
                124.60346122403192,
                0.0119478059177848,
                -34.344520503316495,
            ) |> Collections.fieldvalues;
            atol = 1e-3,
        )
        @test isapprox(
            lsqfit(
                EnergyForm(),
                BirchMurnaghan3rd(110, 0.01, 4, -34),
                volumes,
                energies,
            ) |> Collections.fieldvalues,
            BirchMurnaghan3rd(
                124.82366127014902,
                0.011559181270548115,
                4.115715176896173,
                -34.344311191717864,
            ) |> Collections.fieldvalues;
            atol = 1e-3,
        )
        @test isapprox(
            lsqfit(
                EnergyForm(),
                BirchMurnaghan4th(124, 0.01, 4, -5300, -34),
                volumes,
                energies,
            ) |> Collections.fieldvalues,
            BirchMurnaghan4th(
                124.82134004287863,
                0.011506626229629386,
                4.171189921447817,
                -373.8291670908353,
                -34.34428288861632,
            ) |> Collections.fieldvalues;
            atol = 1e-3,
        )
        @test isapprox(
            lsqfit(EnergyForm(), Vinet(124, 0.01, 4, -34), volumes, energies) |>
            Collections.fieldvalues,
            Vinet(
                124.78343088049905,
                0.011244915521226076,
                4.5244363120324955,
                -34.344139097248,
            ) |> Collections.fieldvalues;
            atol = 1e-3,
        )
    end # testset

    @testset "with units" begin
        volumes = data[:, 1] .* u"angstrom^3"
        energies = data[:, 2] .* u"Ry"
        fitted_eos = lsqfit(
            EnergyForm(),
            BirchMurnaghan3rd(110 * u"angstrom^3", 0.01 * u"Ry/angstrom^3", 4, -34 * u"Ry"),
            volumes,
            energies,
        )
        @test ustrip.(fitted_eos |> Collections.fieldvalues) ≈ ustrip.(BirchMurnaghan3rd(
            124.82366127026123 * u"angstrom^3",
            0.011559181270413301 * u"Ry/angstrom^3",
            4.1157151769280995,
            -34.34431119171783 * u"Ry",
        ) |> Collections.fieldvalues)
        @test ustrip.(lsqfit(
            EnergyForm(),
            BirchMurnaghan3rd(124 * u"angstrom^3", 20 * u"GPa", 4, -17 * u"hartree"),
            volumes,
            energies,
        ) |> Collections.fieldvalues) ≈ ustrip.(BirchMurnaghan3rd(
            124.82366126764893 * u"angstrom^3",
            25.19753977636424 * u"GPa",
            4.115715175957134,
            -17.172155595859646 * u"hartree",
        ) |> Collections.fieldvalues)
    end
end

@testset "`Test test01a.dat` from `Gibbs2`" begin
    data = [
        31.5176 -14.5664202
        32.4533 -14.5804702
        33.4073 -14.5938196
        34.3799 -14.6065011
        35.3711 -14.618562
        36.3812 -14.6300054
        37.4103 -14.640879
        38.4587 -14.6512193
        39.5265 -14.6610287
        40.6139 -14.670344
        41.721 -14.6792001
        42.8481 -14.687601
        43.9952 -14.695594
        45.1627 -14.703158
        46.3507 -14.7103533
        47.5593 -14.717181
        48.7887 -14.723658
        50.0392 -14.7298098
        51.3108 -14.7356349
        52.6038 -14.7411619
        53.9183 -14.7463919
        55.2545 -14.7513556
        56.6127 -14.7560566
        57.9929 -14.7605144
        59.3954 -14.7647283
        60.8203 -14.7687083
        62.2678 -14.7724766
        63.7381 -14.7760428
        65.2313 -14.7793983
        66.7478 -14.7825746
        68.2875 -14.7855695
        69.8507 -14.7883851
        71.4376 -14.7910427
        73.0484 -14.7935382
        74.6832 -14.7958923
        76.3421 -14.7980946
        78.0255 -14.8001569
        79.7335 -14.8020929
        81.4662 -14.8039003
        83.2238 -14.8055857
        85.0065 -14.8071605
        86.8145 -14.8086231
        88.648 -14.8099792
        90.5071 -14.8112335
        92.392 -14.8123882
        94.3029 -14.8134524
        96.2399 -14.8144362
        98.2033 -14.8153196
        100.1932 -14.8161318
        102.2099 -14.8168626
        104.2534 -14.8175204
        106.3239 -14.8181053
        108.4217 -14.8186225
        110.547 -14.8190725
        112.6998 -14.8194614
        114.8803 -14.8197917
        117.0889 -14.8200626
        119.3255 -14.820278
        121.5905 -14.8204398
        123.8839 -14.8205529
        126.206 -14.8206173
        128.5569 -14.8206353
        130.9369 -14.8206092
        133.346 -14.8205405
        135.7845 -14.8204307
        138.2526 -14.820282
        140.7504 -14.8200967
        143.2781 -14.8198746
        145.8358 -14.8196193
        148.4239 -14.819331
        151.0423 -14.8190113
        153.6914 -14.818666
        156.3713 -14.8182878
        159.0822 -14.8178823
        161.8242 -14.8174505
        164.5975 -14.816994
        167.4024 -14.8165098
        170.2389 -14.8160064
        173.1073 -14.8154811
        176.0077 -14.8149353
        178.9404 -14.8143691
        181.9054 -14.8137844
        184.903 -14.8131807
        187.9334 -14.8125612
        190.9967 -14.8119241
        194.0931 -14.8112717
        197.2228 -14.8106045
        200.386 -14.8099231
        203.5828 -14.8092275
        206.8134 -14.8085189
        210.078 -14.8077991
        213.3768 -14.807068
        216.71 -14.8063242
        220.0777 -14.8055711
        223.4801 -14.8048086
        226.9174 -14.8040364
        230.3897 -14.8032545
        233.8973 -14.802466
        237.4403 -14.8016672
        241.0189 -14.8008634
        244.6333 -14.8000517
        248.2837 -14.7992338
        251.9701 -14.7984082
    ]

    @testset "without unit" begin
        volumes = data[:, 1]  # unit: bohr^3
        energies = data[:, 2]  # unit: Rydberg
        @test isapprox(
            lsqfit(EnergyForm(), Murnaghan(132, 0.01, 3.68, -14), volumes, energies) |>
            Collections.fieldvalues,
            Murnaghan(
                132.9174710492377,
                0.000997071972650323,
                2.6506133298092553,
                -14.821067528760938,
            ) |> Collections.fieldvalues;
            atol = 1e-3,
        )
        @test isapprox(
            lsqfit(EnergyForm(), BirchMurnaghan2nd(132, 0.01, -14), volumes, energies) |>
            Collections.fieldvalues,
            BirchMurnaghan2nd(
                130.7191505712864,
                0.000706623449967504,
                -14.817402887854884,
            ) |> Collections.fieldvalues;
            atol = 1e-3,
        )
        @test isapprox(
            lsqfit(
                EnergyForm(),
                BirchMurnaghan3rd(128, 0.03, 4, -14),
                volumes,
                energies,
            ) |> Collections.fieldvalues,
            BirchMurnaghan3rd(
                126.49515516259525,
                0.0010084167615290376,
                3.679199350825455,
                -14.820223865421232,
            ) |> Collections.fieldvalues;
            atol = 1e-3,
        )
        @test isapprox(
            lsqfit(
                EnergyForm(),
                BirchMurnaghan4th(128, 0.03, 4, -320, -14),
                volumes,
                energies,
            ) |> Collections.fieldvalues,
            BirchMurnaghan4th(
                127.67097934611786,
                0.001041910691949355,
                3.463365305040609,
                -3214.0364009441837,
                -14.820651678910448,
            ) |> Collections.fieldvalues;
            atol = 1e-3,
        )
        @test isapprox(
            lsqfit(EnergyForm(), Vinet(128, 0.03, 4, -14), volumes, energies) |>
            Collections.fieldvalues,
            Vinet(
                124.71725873851614,
                0.001064866793589798,
                3.905491499988182,
                -14.820340399744088,
            ) |> Collections.fieldvalues;
            atol = 1e-3,
        )
    end # testset

    @testset "with units" begin
        volumes = data[:, 1] .* u"bohr^3"
        energies = data[:, 2] .* u"Ry"
        fitted_eos = lsqfit(
            EnergyForm(),
            BirchMurnaghan3rd(128 * u"bohr^3", 0.03 * u"Ry/bohr^3", 4, -14 * u"Ry"),
            volumes,
            energies,
        )
        @test ustrip.(fitted_eos |> Collections.fieldvalues) ≈ ustrip.(BirchMurnaghan3rd(
            126.49515516270048 * u"bohr^3",
            0.001008416761528213 * u"Ry/bohr^3",
            3.6791993508231235,
            -14.820223865421246 * u"Ry",
        ) |> Collections.fieldvalues)
        @test ustrip.(lsqfit(
            EnergyForm(),
            BirchMurnaghan3rd(128 * u"bohr^3", 44 * u"GPa", 4, -14 * u"hartree"),
            volumes,
            energies,
        ) |> Collections.fieldvalues) ≈ ustrip.(BirchMurnaghan3rd(
            126.49515516262232 * u"bohr^3",
            14.834322684855008 * u"GPa",
            3.679199350823988,
            -7.410111932710624 * u"hartree",
        ) |> Collections.fieldvalues)
    end
end