text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#
"""
$(TYPEDEF)

Linear matrix inequality cone parametrized by list of real symmetric or complex
Hermitian matrices `mats` of equal dimension.

    $(FUNCTIONNAME){T}(mats::Vector, use_dual::Bool = false)
"""
mutable struct LinMatrixIneq{T <: Real} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int
    side::Int
    As::Vector

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    hess_fact_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}
    use_hess_prod_slow::Bool
    use_hess_prod_slow_updated::Bool

    sumA::AbstractMatrix
    fact
    sumAinvAs::Vector

    function LinMatrixIneq{T}(
        As::Vector;
        use_dual::Bool = false,
        ) where {T <: Real}
        dim = length(As)
        @assert dim > 1
        side = 0
        for A_i in As
            if A_i isa AbstractMatrix
                if iszero(side)
                    side = size(A_i, 1)
                else
                    @assert size(A_i, 1) == side
                end
            end
            @assert ishermitian(A_i)
        end
        @assert side > 0
        # necessary to ensure linear independence of As (but not sufficient)
        @assert svec_length(side) >= dim
        @assert isposdef(first(As))
        cone = new{T}()
        cone.use_dual_barrier = use_dual
        cone.dim = dim
        cone.side = side
        cone.As = As
        return cone
    end
end

reset_data(cone::LinMatrixIneq) = (cone.feas_updated = cone.grad_updated =
    cone.hess_updated = cone.inv_hess_updated = cone.hess_fact_updated =
    cone.use_hess_prod_slow = cone.use_hess_prod_slow_updated = false)

get_nu(cone::LinMatrixIneq) = cone.side

function set_initial_point!(
    arr::AbstractVector,
    cone::LinMatrixIneq{T},
    ) where {T <: Real}
    arr .= 0
    arr[1] = 1
    return arr
end

lmi_fact(arr::AbstractSparseMatrix) = cholesky(Hermitian(arr),
    shift=false, check=false)
lmi_fact(arr::AbstractMatrix) = cholesky!(Hermitian(arr), check=false)

function update_feas(cone::LinMatrixIneq{T}) where {T <: Real}
    @assert !cone.feas_updated

    cone.sumA = sum(w_i * A_i for (w_i, A_i) in zip(cone.point, cone.As))
    cone.fact = lmi_fact(cone.sumA)
    cone.is_feas = isposdef(cone.fact)

    cone.feas_updated = true
    return cone.is_feas
end

function update_grad(cone::LinMatrixIneq{T}) where {T <: Real}
    @assert cone.is_feas

    L = cone.fact.L
    cone.sumAinvAs = [Hermitian(L \ (L \ A_i)', :U) for A_i in cone.As]
    @inbounds for (i, mat_i) in enumerate(cone.sumAinvAs)
        cone.grad[i] = -tr(mat_i)
    end

    cone.grad_updated = true
    return cone.grad
end

function update_hess(cone::LinMatrixIneq)
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)
    H = cone.hess.data
    sumAinvAs = cone.sumAinvAs

    @inbounds for i in 1:cone.dim, j in i:cone.dim
        H[i, j] = real(dot(sumAinvAs[i], sumAinvAs[j]'))
    end

    cone.hess_updated = true
    return cone.hess
end

function hess_prod_slow!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::LinMatrixIneq,
    )
    cone.use_hess_prod_slow_updated || update_use_hess_prod_slow(cone)
    @assert cone.hess_updated
    cone.use_hess_prod_slow || return hess_prod!(prod, arr, cone)

    @assert cone.grad_updated
    sumAinvAs = cone.sumAinvAs

    @inbounds for j in 1:size(arr, 2)
        j_mat = Hermitian(sum(arr[i, j] * sumAinvAs[i] for i in 1:cone.dim))
        for i in 1:cone.dim
            prod[i, j] = real(dot(j_mat, sumAinvAs[i]))
        end
    end

    return prod
end

function dder3(cone::LinMatrixIneq, dir::AbstractVector)
    @assert cone.grad_updated
    dder3 = cone.dder3
    sumAinvAs = cone.sumAinvAs

    dir_mat = sum(d_i * mat_i for (d_i, mat_i) in zip(dir, sumAinvAs))
    Z = Hermitian(dir_mat * dir_mat')
    @inbounds for i in 1:cone.dim
        dder3[i] = real(dot(Z, sumAinvAs[i]))
    end

    return dder3
end


text:
sig_test(line, char) = LanguageServer.textDocument_signatureHelp_request(LanguageServer.TextDocumentPositionParams(LanguageServer.TextDocumentIdentifier("testdoc"), LanguageServer.Position(line, char)), server, server.jr_endpoint)

def_test(line, char) = LanguageServer.textDocument_definition_request(LanguageServer.TextDocumentPositionParams(LanguageServer.TextDocumentIdentifier("testdoc"), LanguageServer.Position(line, char)), server, server.jr_endpoint)

ref_test(line, char) = LanguageServer.textDocument_references_request(LanguageServer.ReferenceParams(LanguageServer.TextDocumentIdentifier("testdoc"), LanguageServer.Position(line, char), missing, missing, LanguageServer.ReferenceContext(true)), server, server.jr_endpoint)

rename_test(line, char) = LanguageServer.textDocument_rename_request(LanguageServer.RenameParams(LanguageServer.TextDocumentIdentifier("testdoc"), LanguageServer.Position(line, char), missing, "newname"), server, server.jr_endpoint)

@testset "signatures" begin
    doc = settestdoc("""
    rand()
    Base.rand()
    func(arg) = 1
    func()
    struct T
        a
        b
    end
    T()
    struct S{R}
        a
        S() = new(1)
    end
    using Base:argtail
    argtail()
    S{R}()
    """)
    @test !isempty(sig_test(0, 5).signatures)
    @test !isempty(sig_test(1, 10).signatures)
    @test !isempty(sig_test(3, 5).signatures)
    @test !isempty(sig_test(8, 2).signatures)
    @test_broken !isempty(sig_test(15, 5).signatures)

    let sigs = LanguageServer.SignatureInformation[]
        LanguageServer.get_signatures(doc.cst[3].meta.binding, doc.cst.meta.scope, sigs, server)
        @test length(sigs) == 1
    end
    let sigs = LanguageServer.SignatureInformation[]
        LanguageServer.get_signatures(doc.cst[5].meta.binding, doc.cst.meta.scope, sigs, server)
        @test length(sigs) == 1
    end
    let sigs = LanguageServer.SignatureInformation[]
        LanguageServer.get_signatures(doc.cst[1][1].meta.ref, doc.cst.meta.scope, sigs, server)
        @test length(sigs) > 0
    end
    let sigs = LanguageServer.SignatureInformation[]
        LanguageServer.get_signatures(doc.cst[7].meta.binding, doc.cst.meta.scope, sigs, server)
        @test length(sigs) == 1
    end
    let sigs = LanguageServer.SignatureInformation[]
        LanguageServer.get_signatures(doc.cst[9][1].meta.ref, doc.cst.meta.scope, sigs, server)
        @test length(sigs) == 1
    end
end

@testset "definitions" begin
    settestdoc("""
    rand()
    func(arg) = 1
    func()
    Float64
    """)
    # @test !isempty(def_test(0, 3))
    @test !isempty(def_test(2, 3))
    @test !isempty(def_test(3, 3))
end

@testset "references" begin
    settestdoc("""
    func(arg) = 1
    func()
    """)
    @test length(ref_test(1, 2)) == 2
end

@testset "rename" begin
    settestdoc("""
    func(arg) = 1
    func()
    """)
    @test length(rename_test(0, 2).documentChanges[1].edits) == 2
end

@testset "get_file_loc" begin
    doc = settestdoc("""
    func(arg) = 1
    func()
    """)
    @test LanguageServer.get_file_loc(doc.cst.args[2].args[1]) == (doc, 14)
end

@testset "doc symbols" begin
    doc = settestdoc("""
    a = 1
    b = 2
    function func() end
    """)
    @test all(item.name in ("a", "b", "func") for item in LanguageServer.textDocument_documentSymbol_request(LanguageServer.DocumentSymbolParams(LanguageServer.TextDocumentIdentifier("testdoc"), missing, missing), server, server.jr_endpoint))
end

ext: jl


text:
export iauHfk5z
"""
Transform a Hipparcos star position into FK5 J2000.0, assuming
zero Hipparcos proper motion.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   rh            double    Hipparcos RA (radians)
   dh            double    Hipparcos Dec (radians)
   date1,date2   double    TDB date (Note 1)

Returned (all FK5, equinox J2000.0, date date1+date2):
   r5            double    RA (radians)
   d5            double    Dec (radians)
   dr5           double    FK5 RA proper motion (rad/year, Note 4)
   dd5           double    Dec proper motion (rad/year, Note 4)

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The proper motion in RA is dRA/dt rather than cos(Dec)*dRA/dt.

   3. The FK5 to Hipparcos transformation is modeled as a pure rotation
   and spin;  zonal errors in the FK5 catalogue are not taken into
   account.

   4. It was the intention that Hipparcos should be a close
   approximation to an inertial frame, so that distant objects have
   zero proper motion;  such objects have (in general) non-zero
   proper motion in FK5, and this function returns those fictitious
   proper motions.

   5. The position returned by this function is in the FK5 J2000.0
   reference system but at date date1+date2.

   6. See also iauFk52h, iauH2fk5, iauFk5zhz.

Called:
   iauS2c       spherical coordinates to unit vector
   iauFk5hip    FK5 to Hipparcos rotation and spin
   iauRxp       product of r-matrix and p-vector
   iauSxp       multiply p-vector by scalar
   iauRxr       product of two r-matrices
   iauTrxp      product of transpose of r-matrix and p-vector
   iauPxp       vector product of two p-vectors
   iauPv2s      pv-vector to spherical
   iauAnp       normalize angle into range 0 to 2pi

Reference:

   F.Mignard & M.Froeschle, 2000, Astron.Astrophys. 354, 732-739.

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauHfk5z(double rh, double dh, double date1, double date2,
#               double *r5, double *d5, double *dr5, double *dd5)

function iauHfk5z(rh::Real, dh::Real, date1::Real, date2::Real)
   ref_r5  = Ref{Float64}(0.0)
   ref_d5  = Ref{Float64}(0.0)
   ref_dr5 = Ref{Float64}(0.0)
   ref_dd5 = Ref{Float64}(0.0)

   ccall((:iauHfk5z, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Cdouble, Cdouble,
         Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
          convert(Float64, rh), convert(Float64, dh),
          convert(Float64, date1), convert(Float64, date2),
          ref_r5, ref_d5, ref_dr5, ref_dd5)

   return ref_r5[], ref_d5[], ref_dr5[], ref_dd5[]
end


text:
using Optim

using Distributions

# generate simple nonlinear dynamics
plant(x) = -x.^2

# Hyperparameters
K = 2  # number of Gaussians
d = 5  # dimension

batch_size = 100

# Generate some data
X = rand((batch_size, d))
V = plant(X)

function mixture(x, w, μ, Σ)
  normal_pdf(i) = pdf(MvNormal(μ[i, :], Σ[i, :]), x)
  return sum(weight*normal_pdf(i) for (i, weight) in enumerate(w))
end

function mse(G)
  w = G[:, 1]
  μ = G[:, 2:d+1]
  Σ = G[:, d+2:2*d+1]

  return sum(norm(V[i, :] - mixture(X[i, :], w, μ, Σ))^2 for i=1:batch_size)
end

# TODO non-unary callables in Optim?
# This may make some dumb assumptions
G = rand(K, 2*d+1)
optimize(mse, G)

# This doesn't work because nonlinear expressions don't allow vector data in JuMP :(
#=
using JuMP

model = Model()

@variable(model, W[1:K])  # Weights
@variable(model, μ[1:K, 1:d])  # Means
@variable(model, Σ[1:K, 1:d])  # Covariance diagonals

gaussian_pdf = @NLexpression(model, [i=1:batch_size, k=1:K],
                             exp(-(X[i, :] - μ[k, :])^2/(2*Σ[k, :]^2)/sqrt(2*π*Σ[k, :]^2)))

mixture = @NLexpression(model, [i=1:batch_size], sum([W[k]*(gaussian_pdf[i, k]) for k=1:K]))

@NLobjective(model, Min, sum([ norm(V[j, :] - mixture[j])^2 for j=1:batch_size]))
solve(model)
=#

max_stars_repo_name: nefthys/ModelLearning

max_issues_repo_name: nefthys/ModelLearning

max_forks_repo_name: nefthys/ModelLearning


text:
using CSyntax

@static if Sys.isapple()
    const VERSION_MAJOR = 4
    const VERSION_MINOR = 1
end

include(joinpath(@__DIR__, "glutils.jl"))

# init window
width, height = fb_width, fb_height = 640, 480
window = startgl(width, height)

glEnable(GL_DEPTH_TEST)
glDepthFunc(GL_LESS)

# vertex data
points = GLfloat[ 0.0,  0.5, 0.0,
                  0.5, -0.5, 0.0,
                 -0.5, -0.5, 0.0]

colors = GLfloat[ 1.0, 0.0, 0.0,
                  0.0, 1.0, 0.0,
                  0.0, 0.0, 1.0]

# create buffers located in the memory of graphic card
points_vbo = GLuint(0)
@c glGenBuffers(1, &points_vbo)
glBindBuffer(GL_ARRAY_BUFFER, points_vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(points), points, GL_STATIC_DRAW)

colors_vbo = GLuint(0)
@c glGenBuffers(1, &colors_vbo)
glBindBuffer(GL_ARRAY_BUFFER, colors_vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(colors), colors, GL_STATIC_DRAW)

# create VAO
vao = GLuint(0)
@c glGenVertexArrays(1, &vao)
glBindVertexArray(vao)
glBindBuffer(GL_ARRAY_BUFFER, points_vbo)
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0, C_NULL)
glBindBuffer(GL_ARRAY_BUFFER, colors_vbo)
glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, 0, C_NULL)
glEnableVertexAttribArray(0)
glEnableVertexAttribArray(1)

# load and compile shaders from file
vert_shader = createshader(joinpath(@__DIR__, "camera.vert"), GL_VERTEX_SHADER)
frag_shader = createshader(joinpath(@__DIR__, "camera.frag"), GL_FRAGMENT_SHADER)

# link program
shader_prog = createprogram(vert_shader, frag_shader)

# enable cull face
glEnable(GL_CULL_FACE)
glCullFace(GL_BACK)
glFrontFace(GL_CW)
# set background color to gray
glClearColor(0.2, 0.2, 0.2, 1.0)

# camera
near = 0.1            # clipping near plane
far = 100.0           # clipping far plane
fov = deg2rad(67)
aspect_ratio = width / height
# perspective matrix
range = tan(0.5*fov) * near
Sx = 2.0*near / (range * aspect_ratio + range * aspect_ratio)
Sy = near / range
Sz = -(far + near) / (far - near)
Pz = -(2.0*far*near) / (far - near)
proj_matrix = GLfloat[ Sx  0.0  0.0  0.0;
                      0.0   Sy  0.0  0.0;
                      0.0  0.0   Sz   Pz;
                      0.0  0.0 -1.0  0.0]
# view matrix
camera_speed = GLfloat(1.0)
camera_yaw_speed = GLfloat(10.0)
camera_pos = GLfloat[0.0, 0.0, 2.0]   # don't start at zero, or we will be too close
camera_yaw = GLfloat(0.0)             # y-rotation in degrees
trans_matrix = GLfloat[1.0 0.0 0.0 -camera_pos[1];
                       0.0 1.0 0.0 -camera_pos[2];
                       0.0 0.0 1.0 -camera_pos[3];
                       0.0 0.0 0.0            1.0]
rotation_y = GLfloat[cosd(-camera_yaw)  0.0  sind(-camera_yaw)  0.0;
                                   0.0  1.0                0.0  0.0;
                    -sind(-camera_yaw)  0.0  cosd(-camera_yaw)  0.0;
                                   0.0  0.0                0.0  1.0]
view_matrix = rotation_y * trans_matrix    # only rotate around the Y axis

view_matrix_loc = glGetUniformLocation(shader_prog, "view")
proj_matrix_loc = glGetUniformLocation(shader_prog, "proj")
glUseProgram(shader_prog)
glUniformMatrix4fv(view_matrix_loc, 1, GL_FALSE, view_matrix)
glUniformMatrix4fv(proj_matrix_loc, 1, GL_FALSE, proj_matrix)

let
updatefps = FPSCounter()
# render
previous = time()
while !GLFW.WindowShouldClose(window)
    updatefps(window)
    # clear drawing surface
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0, 0, GLFW.GetFramebufferSize(window)...)
    # drawing
    glUseProgram(shader_prog)
    glBindVertexArray(vao)
    glDrawArrays(GL_TRIANGLES, 0, 3)
    # check and call events
    GLFW.PollEvents()
    # move camera
    global camera_pos; global camera_speed; global camera_yaw
    current = time()
    elapsed = current - previous
    previous = current
    GLFW.GetKey(window, GLFW.KEY_A) && (camera_pos[1] -= camera_speed * elapsed)
    GLFW.GetKey(window, GLFW.KEY_D) && (camera_pos[1] += camera_speed * elapsed)
    GLFW.GetKey(window, GLFW.KEY_PAGE_UP) && (camera_pos[2] += camera_speed * elapsed)
    GLFW.GetKey(window, GLFW.KEY_PAGE_DOWN) && (camera_pos[2] -= camera_speed * elapsed)
    GLFW.GetKey(window, GLFW.KEY_W) && (camera_pos[3] -= camera_speed * elapsed)
    GLFW.GetKey(window, GLFW.KEY_S) && (camera_pos[3] += camera_speed * elapsed)
    GLFW.GetKey(window, GLFW.KEY_LEFT) && (camera_yaw += camera_yaw_speed * elapsed)
    GLFW.GetKey(window, GLFW.KEY_RIGHT) && (camera_yaw -= camera_yaw_speed * elapsed)
    # update view matrix
    trans_matrix .= GLfloat[1.0 0.0 0.0 -camera_pos[1];
                            0.0 1.0 0.0 -camera_pos[2];
                            0.0 0.0 1.0 -camera_pos[3];
                            0.0 0.0 0.0            1.0]
    rotation_y .= GLfloat[cosd(-camera_yaw)  0.0  sind(-camera_yaw) 0.0;
                                       0.0   1.0               0.0  0.0;
                         -sind(-camera_yaw)  0.0  cosd(-camera_yaw) 0.0;
                                       0.0   0.0               0.0  1.0]
    view_matrix .= rotation_y * trans_matrix
    glUniformMatrix4fv(view_matrix_loc, 1, GL_FALSE, view_matrix)
    # swap the buffers
    GLFW.SwapBuffers(window)
end
end # let

GLFW.DestroyWindow(window)

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
using Poisson2D
using JLD2
using Printf

# set problem parameters
a = 0;
b = 1;
c = 0;
d = 1;
n = 999;

kx = 5;
ky = 4;
f = (x,y)->((kx*π/(b-a))^2 + (ky * π/(d-c))^2) * sin(kx*π*(x-a)/(b-a)) * sin(ky*π*(y-c)/(d-c));

# set up and solve problem
problem = FDPoisson2DProblem(a,b,c,d,n,f);
assemble_system!(problem);
U = direct_solve_poisson2d(problem);

# name of file
fname = @sprintf("soln1_n%d.jld2", n);

# save to disk
jldsave(fname; problem, U);


text:
module JuliaZH

include("basedocs.jl")

end # module


text:
"""
    CubicRamanTensor() -> Tensor{4, 3}

    Returns a dimensionless tensor describing the Raman
    contribution for cubic crystals
"""
function CubicRamanTensor()
    δ = (i, j) -> i == j ? 1 : 0
    δ4 = (i, j, k, l) -> i == j && j == k && k == l ? 1 : 0
    f = (i, j, k, l) -> δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k) - 2*δ4(i,j,k,l)
    return Tensor{4, 3}(f)
end

"""
    CubicElectronicTensor(ρ) -> Tensor{4, 3}

    Returns a dimensionless tensor describing the electronic
    contribution for cubic crystals. ρ characterizes the nonlinear
    anisotropy ρ ≡ 3χ₁₁₂₂/χ₁₁₁₁
"""
function CubicElectronicTensor(ρ)
    δ = (i, j) -> i == j ? 1 : 0
    δ4 = (i, j, k, l) -> i == j && j == k && k == l ? 1 : 0
    f = (i, j, k, l) -> ρ/3*(δ(i,j)*δ(k,l) + δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) + (1-ρ)*δ4(i,j,k,l)
    return Tensor{4, 3}(f)
end

"""
    IsotropicTensor(λ, μ, ν) -> Tensor{4, 3}
    IsotropicTensor()        -> Tensor{4, 3}

    Returns a dimensionless tensor describing the
    an isotropic nonlinear contribution. λ, μ, and ν are
    parameters which control the coefficients, see
    https://farside.ph.utexas.edu/teaching/336L/Fluid/node252.html
"""
function IsotropicTensor(λ, μ, ν)
    δ = (i, j) -> i == j ? 1 : 0
    f = (i, j, k, l) -> (λ*δ(i,j)*δ(k,l) + μ*δ(i,k)*δ(j,l) + ν*δ(i,l)*δ(j,k))
    return Tensor{4, 3}(f)
end
IsotropicTensor() = IsotropicTensor(1,1,1)

const Diamond = Material(
                    ElectronicTensor(
                        CubicElectronicTensor(1.2),
                        1.3e-19m^2/W
                    ),
                    RamanTensor(
                        CubicRamanTensor(),
                        0.004ps,
                        5.7ps,
                        0.28
                    )
                )
const Silicon = Material(
                    ElectronicTensor(
                        CubicElectronicTensor(1.27),
                        3.8e-18m^2/W
                    ),
                    RamanTensor(
                        CubicRamanTensor(),
                        0.01ps,
                        3.0ps,
                        0.08
                    )
                )
const SiO2 = Material(
                    ElectronicTensor(
                        IsotropicTensor(),
                        2.7e-20m^2/W,
                    ),
                    RamanTensor(
                        IsotropicTensor(),
                        0.0122ps,
                        0.032ps,
                        0.18
                    )
                )


text:
include("Types.jl")
include("Interfaces.jl")
include("TypeParser.jl")
include("Base.jl")
include("DateTime.jl")
include("DateTime64.jl")
include("Date.jl")
include("Enum.jl")
include("FixedString.jl")
include("Tuple.jl")
include("UUID.jl")
include("Nullable.jl")
include("LowCardinality.jl")
include("Array.jl")
include("Nothing.jl")
include("SimpleAggregateFunction.jl")
include("IPv4.jl")
include("IPv6.jl")
include("Decimal.jl")


text:
module SpinQubits

    using LinearAlgebra, Statistics, Distributions, Test, Plots, TimerOutputs

    import Base: +, *, isequal, ==

    export calculateFidelities, saveFidelities, averager, plotter!, readmathematica, to

    const to = TimerOutput()

    include("spinors.jl")
    include("operators.jl")
    include("utils.jl")
    include("tensors.jl")
    include("IO.jl")

    @timeit to function calculateFidelities(L::Int64, β::Float64, disGam, sigmas, nReals::Int64, spacing::Float64; singlet=false, betas=[], verbose=true)

        # construct exponents; Scale στ to our units
        j0 = 1.0
        exponents = collect(range(0.0,3.0,step=spacing))
        δt = sigmas[3]*1e-9
        tShortestNichol = 10e-9
        δrel = δt/tShortestNichol
        tShortestMine = pi/*(4.0 * j0 * 10.0^exponents[end])
        sigmas[3] = δrel*tShortestMine

        nIterations = (maximum(sigmas) == 0) ? 1 : nReals

        # Initialize collections for building hamiltonian
        jtensor = getjtensor(L, β; betaArray=betas)
        numJs = Int(L*(L-1)/2) # n + (n-1) + (n-2) + ... = n(n+1)/2
        js = zeros(numJs)

        #sequence = singlet ? [collect(2:L-1); collect((L-1):-1:1); 1] : [collect(1:L-1);collect((L-1):-1:1)]
        sequence = [collect(1:L-1);collect((L-1):-1:1)] # Made sequence same for singlet and eigenstate

        # Find initial ket; pre allocate resulting kets after evolution
        initKet = getInitKet(L, singlet)
        finalKet::Vector{ComplexF64} = deepcopy(initKet)
        currentKet::Vector{ComplexF64} = deepcopy(initKet)

        # Pre allocation
        singleExpFidelities = zeros(nIterations)
        fidelities = zeros(length(exponents))
        D = zeros(2^L,2^L)
        R = zeros(2^L,2^L)
        ham::Matrix{ComplexF64} = zeros(2^L,2^L)
        diss::Matrix{ComplexF64} = Diagonal(zeros(2^L))
        UD::Matrix{ComplexF64} = zeros(2^L,2^L)
        U::Matrix{ComplexF64} = zeros(2^L,2^L)
        trueU::Matrix{ComplexF64} = zeros(2^L,2^L)
        j0s = fill(j0,(numJs, nIterations))
        γs = fill(disGam,nIterations)
        τs = zeros(nIterations)
        for expIndex in 1:length(exponents)

            exponent = exponents[expIndex]
            (verbose && expIndex % 25 == 0) ? println("Calculating ",expIndex,"th exponent out of ",length(exponents)) : nothing
            jSWAP = 10.0^exponent

            incorporateNoise!(j0s, γs, τs, sigmas, disGam, jSWAP, j0)

            for i in 1:nIterations
                currentKet .= initKet
                finalKet .= initKet

                for j in 1:(2^L)
                    diss[j,j] = (-im * (-im*γs[i]) * τs[i])
                end
                diagexp!(diss)

                for swapIndex in sequence
                    @timeit to "setting up mats" begin
                        zeros!(U)
                        zeros!(ham)
                        @timeit to "updating js and ham" begin
                            currentKet .= finalKet
                            js[:] .= @view j0s[:,i]
                            baseIndex = 0
                            for k in 1:L-1 # k is the interspin distance
                                if swapIndex <= L-k
                                    js[baseIndex + swapIndex] *= jSWAP/j0
                                end
                                if swapIndex > (k-1) && k > 1
                                    js[baseIndex + (swapIndex - (k-1))] *= jSWAP/j0
                                end
                                baseIndex += L-k
                            end
                            Ham!(ham,L,jtensor,js)
                        end
                        @timeit to "eigen" begin
                            eigenObject = eigen!(ham)
                            D .= Diagonal(eigenObject.values)

                            R .= eigenObject.vectors
                            UD .= (-im .* D .* τs[i])
                        end
                    end
                    @timeit to "matmuls" begin
                        diagexp!(UD)
                        mul!(ham,UD,transpose(R)) # ham here is just used as dummy memory space, not the hamiltonian
                        mul!(U,R,ham) # ham here is just used as dummy memory space, not the hamiltonian

                        mul!(trueU,U,diss) # ham here is just used as dummy memory space, not the hamiltonian
                        mul!(finalKet,trueU,currentKet)
                    end
                end # sequence
                singleExpFidelities[i] = abs2(initKet'*finalKet)
            end # average
            fidelities[expIndex] = mean(singleExpFidelities)
        end # exponents
        return 10 .^exponents,1 .-fidelities
    end
end


text:
#---------------------------------------------------------------------#
#This function computes the high-order grid & elements in 2D.
#Written by F.X. Giraldo on 5/2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#

using Plots
include("warp_grid.jl")
include("rotate_grid.jl")

function create_grid(Np,Npoin,Ne,Nboun,Nex,Ney,ξ,plot_grid,warp_grid,rotate_grid,DFloat)

    #Initialize Global Arrays
    coord=zeros(DFloat,2,Npoin)
    intma=zeros(Int64,Np,Np,Ne)
    bsido=zeros(Int64,4,Nboun)
    periodicity=zeros(Int64,Npoin,1)

    #Initialize Local Arrays
    node=zeros(Int64,Npoin,Npoin)

    #Set some constants
    xmin=-1
    xmax=+1
    ymin=-1
    ymax=+1
    dx=(xmax-xmin)/Nex
    dy=(ymax-ymin)/Ney
    N=Np-1
    nx=Nex*N + 1
    ny=Ney*N + 1

    #GENERATE COORD
    I=0
    jj=0
    for k=1:Ney
        y0=ymin + (k-1)*dy
        if (k == 1)
            l1=1
        else
            l1=2
        end
        for l=l1:Np
            y=( ξ[l]+1 )*dy/2 + y0
            jj=jj+1
            ii=0
            for i=1:Nex
                x0=xmin + (i-1)*dx
                if (i == 1)
                    j1=1
                else
                    j1=2
                end
                for j=j1:Np
                    ii=ii + 1
                    I=I + 1
                    x=( ξ[j]+1 )*dx/2 + x0
                    coord[1,I]=x
                    coord[2,I]=y
                    node[ii,jj]=I
                end #j

            end #i
        end #l
    end #k

    #GENERATE INTMA
    e=0
    for k=1:Ney
        for i=1:Nex
            e=e+1
            for l=1:Np
                jj=N*(k-1) + l
                for j=1:Np
                    ii=N*(i-1) + j
                    I=node[ii,jj]
                    intma[j,l,e]=I
                end #j
            end #l
        end #i
    end #k

    #Generate BSIDO
    b=0
    for i=1:Nex
        e=i
        b=b+1
        i1=(i-1)*(Np-1) + 1
        i2=(i-1)*(Np-1) + Np
        I1=node[i1,1]
        I2=node[i2,1]
        bsido[1,b]=I1
        bsido[2,b]=I2
        bsido[3,b]=e
        bsido[4,b]=6
    end

    #Right Boundary
    for i=1:Ney
        e=(Nex)*(i)
        b=b+1
        i1=(i-1)*(Np-1) + 1
        i2=(i-1)*(Np-1) + Np
        I1=node[nx,i1]
        I2=node[nx,i2]
        bsido[1,b]=I1
        bsido[2,b]=I2
        bsido[3,b]=e
        bsido[4,b]=6
    end

    #Top Boundary
    for i=Nex:-1:1
        e=Ne - (Nex - i)
        b=b+1
        i1=(i-1)*(Np-1) + Np
        i2=(i-1)*(Np-1) + 1
        I1=node[i1,ny]
        I2=node[i2,ny]
        bsido[1,b]=I1
        bsido[2,b]=I2
        bsido[3,b]=e
        bsido[4,b]=6
    end

    #Left Boundary
    for i=Ney:-1:1
        e=(Nex)*(i-1) + 1
        b=b+1
        i1=(i-1)*(Np-1) + Np
        i2=(i-1)*(Np-1) + 1
        I1=node[1,i1]
        I2=node[1,i2]
        bsido[1,b]=I1
        bsido[2,b]=I2
        bsido[3,b]=e
        bsido[4,b]=6
    end

    if (b != Nboun)
        println("Error in Create_Grid: b != Nboun =>  Nboun = ",Nboun," b = ",b)
        exit
    end

    #Periodicity
    for i=1:Npoin
        periodicity[i]=i;
    end

    #X-Periodicity
    for i=1:ny
        i1=node[1,i]
        i2=node[nx,i]
        periodicity[i2]=i1
    end

    #Y-Periodicity
    for i=1:nx
        i1=node[i,1]
        i2=node[i,ny]
        periodicity[i2]=periodicity[i1]
    end

    #Warp Grid
    if (warp_grid)
        warp_grid!(coord,Npoin,DFloat)
    end

    #Rotate Grid
    if (rotate_grid)
        rotate_grid!(coord,Npoin,DFloat)
    end

    #Plot Grid
    title_text=["Grid : Ne =  $Ne,  N = $N"]
    if (plot_grid)
        x=zeros(DFloat,5,1)
        y=zeros(DFloat,5,1)
        x[1]=-1; y[1]=-1
        x[2]=+1; y[2]=-1
        x[3]=+1; y[3]=+1
        x[4]=-1; y[4]=+1
        x[5]=-1; y[5]=-1
        p1=plot(x,y,linecolor="black",legend=false)

        for e=1:Ne
            for j=1:N
                for i=1:N
                    i1=intma[i,j,e]
                    i2=intma[i+1,j,e]
                    i3=intma[i+1,j+1,e]
                    i4=intma[i,j+1,e]
                    x[1]=coord[1,i1]; y[1]=coord[2,i1]
                    x[2]=coord[1,i2]; y[2]=coord[2,i2]
                    x[3]=coord[1,i3]; y[3]=coord[2,i3]
                    x[4]=coord[1,i4]; y[4]=coord[2,i4]
                    x[5]=coord[1,i1]; y[5]=coord[2,i1]
                    p1=plot!(p1,x,y,linecolor="red",legend=false)
                end
            end
            i1=intma[1,1,e]
            i2=intma[Np,1,e]
            i3=intma[Np,Np,e]
            i4=intma[1,Np,e]
            x[1]=coord[1,i1]; y[1]=coord[2,i1]
            x[2]=coord[1,i2]; y[2]=coord[2,i2]
            x[3]=coord[1,i3]; y[3]=coord[2,i3]
            x[4]=coord[1,i4]; y[4]=coord[2,i4]
            x[5]=coord[1,i1]; y[5]=coord[2,i1]
            p1=plot!(p1,x,y,linecolor="blue",legend=false)
        end
        display(plot(p1))
    end

    return(coord,intma,bsido,periodicity)
end


text:
# example to illustate how autoinit functions work, and used as a development script towards a standard unit test.

using Distributions
using IncrementalInference

import IncrementalInference: getSample

## MARKER START is a required block of code, but can be reviewed at the end of the tutorial,
# Run this but skip past these user defined functions for a quicker introduction======================
struct Prior{T} <: IncrementalInference.FunctorSingleton where T <: Distribution
  z::T
end
getSample(s::Prior, N::Int=1) = (rand(s.z,N), )
struct LinearOffset{T} <: IncrementalInference.FunctorPairwise where T <: Distribution
  z::T
end
getSample(s::LinearOffset, N::Int=1) = (rand(s.z,N), )
function (s::LinearOffset)(res::Array{Float64},
      idx::Int,
      meas::Tuple{Array{Float64, 1}},
      X1::Array{Float64,2},
      X2::Array{Float64,2}  )
  #
  res[1] = meas[1][idx] - (X2[1,idx] - X1[1,idx])
  nothing
end
# and a bi-modal conditional function
struct MultiModalOffset <: IncrementalInference.FunctorPairwise
  z::Vector{Distribution}
  c::Categorical
end
getSample(s::MultiModalOffset, N::Int=1) = (rand.(s.z, N)..., rand(s.c, N))
function (s::MultiModalOffset)(res::Array{Float64},
      idx::Int,
      meas::Tuple,
      X1::Array{Float64,2},
      X2::Array{Float64,2}
  )
  res[1] = meas[meas[end][idx]][idx] - (X2[1,idx] - X1[1,idx])
  nothing
end
## Define a model with these user defined functions ===============================================
## MARKER END

# Start with an empty graph
fg = emptyFactorGraph()

# add the first node
addNode!(fg, :x0, ContinuousScalar)

# this is unary (prior) factor and does not immediately trigger autoinit of :x0.
addFactor!(fg, [:x0], Prior(Normal(0,1)))

# writeGraphPdf(fg, file="/home/dehann/Downloads/fgx0.png")

# To visualize the factor graph structure, assuming GraphViz.jl is available.
# Graphs.plot(fg.g)
# Please find `writeGraphPdf` member definition at the end of this tutorial
# writeGraphPdf(fg, file="/home/dehann/Downloads/fgx01.png")

# Also consider automatic initialization of variables.
@show isInitialized(fg, :x0)
# Why is x0 not initialized?
# Since no other variable nodes have been 'connected to' (or depend) on :x0,
# and future intentions of the user are unknown, the initialization of :x0 is
# deferred until such criteria are met.

# Auto initialization of :x0 is triggered by connecting the next (probably uninitialized) variable node
addNode!(fg, :x1, ContinuousScalar)

# with a linear conditional belief to :x0
# P(Z | :x1 - :x0 ) where Z ~ Normal(10,1)
addFactor!(fg, [:x0, :x1], LinearOffset(Normal(10.0,1)))

# writeGraphPdf(fg, file="/home/dehann/Downloads/fgx01.png")

# x0 should not be initialized
@show isInitialized(fg, :x0)

# To plot the belief states of variables -- (assuming `RoMEPlotting` is available)
# Remember the first time executions are slow given required code compilation,
# and that future versions of these package will use more precompilation
# to reduce first execution running cost.
using RoMEPlotting

plotKDE(fg, :x0)

# Since no other variables 'are yet connected'/depend on :x1, it will not be initialized
@show isInitialized(fg, :x1)

# we can force all the variable nodes to initialize
ensureAllInitialized!(fg)

# now draw both :x0
plotKDE(fg, [:x0, :x1])

# add another node, but introduce more general beliefs
addNode!(fg, :x2, ContinuousScalar)

mmo = MultiModalOffset([Rayleigh(3); Uniform(30,55)], Categorical([0.4; 0.6]))
addFactor!(fg, [:x1, :x2], mmo)

# Graphs.plot(fg.g)
# writeGraphPdf(fg, file="/home/dehann/Downloads/fgx012.png")

# By again forcing the initialization of :x3 for illustration
ensureAllInitialized!(fg)

# the predicted marginal probability densities are
plotKDE(fg, [:x0, :x1, :x2])

# Now transmit this 'weird' multi-modal marginal belief through a another unimodal linear offset (conditional likelihood)
addNode!(fg, :x3, ContinuousScalar)

addFactor!(fg, [:x2, :x3], LinearOffset(Normal(-50, 1)))
# note, this addFactor step relies on :x2 being initialized and would have done so if we didn't call ensureAllInitialized! a few lines earlier.

# writeGraphPdf(fg, file="/home/dehann/Downloads/fgx0123.png")

ensureAllInitialized!(fg)
plotKDE(fg, [:x0, :x1, :x2, :x3])

lo3 = LinearOffset(Normal(40, 1))
addFactor!(fg, [:x3, :x0], lo3)

# writeGraphPdf(fg, file="/home/dehann/Downloads/fgx0123c.png")

plotKDE(fg, [:x0, :x1, :x2, :x3])

# Find global best likelihood solution (posterior belief)
# After defining the problem, we can find the 'minimum free energy' solution
tree = wipeBuildNewTree!(fg)
inferOverTree!(fg, tree)

# and look at the posterior belief, and notice which consensus modes stand out in the posterior
plotKDE(fg, [:x0, :x1, :x2, :x3])

# Helper function for graphing the factor graph structure (using GraphViz)

function writeGraphPdf(fgl::FactorGraph; file::String="fg.pdf")
  fgd = drawCopyFG(fgl)
  println("Writing factor graph file")
  fid = open("fg.dot","w+")
  write(fid,Graphs.to_dot(fgd.g))
  close(fid)
  fext = split(file, '.')[end]
  run(`dot fg.dot -T$fext -o $file`)
  Base.rm("fg.dot")
  nothing
end

# using Gadfly
# pl.guides[1] = Gadfly.Guide.xlabel("")
# push!(pl.guides, Gadfly.Guide.ylabel("density"))
#
# Gadfly.draw(PNG("/home/dehann/Downloads/plx012.png", 10cm, 7cm), pl)

## should complete and add to RoMEPlotting
# import KernelDensityEstimatePlotting: plot
# import Gadfly: plot
# import Graphs: plot
# import RoMEPlotting: plot
# function plot(fgl::FactorGraph, sym::Symbol; api::DataLayerAPI=IncrementalInference.dlapi)
#   PX = getKDE(getVert(fgl, sym, api=api))
#   plot(PX)
# end

#


text:
module BlindBake

using Distributed, Dates

include("def.jl")
include("imp.jl")

end # module


text:
contact_normal(obs::Obstacle) = obs.contact_face.outward_normal

RigidBodyDynamics.separation(obs::Obstacle, p::Point3D) = separation(obs.contact_face, p)

function contact_force(r::ContactResult)
    n = contact_normal(r.obs)
    D = contact_basis(r.obs)
    @framecheck(n.frame, D[1].frame)
    v = r.c_n .* n.v
    for i in eachindex(r.β)
        v += r.β[i] .* D[i].v
    end
    FreeVector3D(n.frame, r.scaling * v)
end

function append(v::AbstractVector{T}, x::T) where T
    result = similar(v, length(v) + 1)
    for i in 1:length(v)
        result[i] = v[i]
    end
    result[end] = x
    result
end

function add_contact_constraints(model::Model, μ, β, λ, c_n, D_transpose_times_v, separation_from_obstacle)

    @constraints model begin
        separation_from_obstacle >= 0 # (7)
        μ * c_n - sum(β) >= 0 # (9)
    end

    for d in D_transpose_times_v
        @constraint model λ + d >= 0 # (8)
    end

    @disjunction(model, separation_from_obstacle == 0, c_n == 0) # (10)
    for j in 1:length(D_transpose_times_v)
        d = D_transpose_times_v[j]
        λ_plus_d = AffExpr(append(d.vars, λ), append(d.coeffs, 1.0), d.constant)
        @disjunction(model, (λ_plus_d == 0), (β[j] == 0)) # (11)
    end
    @disjunction(model, (μ * c_n - sum(β) == 0), (λ == 0)) # (12)
end

function resolve_contact(xnext::LinearizedState, body::RigidBody, point::Point3D, obstacle::Obstacle, model::Model)
    mechanism = xnext.linearization_state.mechanism
    total_weight = mass(mechanism) * norm(mechanism.gravitational_acceleration)

    D = contact_basis(obstacle)
    β   = @variable(model, [1:length(D)], lowerbound=0, basename="β", upperbound=100)
    λ   = @variable(model, lowerbound=0, basename="λ", upperbound=100)
    c_n = @variable(model, lowerbound=0, basename="c_n", upperbound=100)
    world = default_frame(root_body(linearization_state(xnext).mechanism))

    separation_from_obstacle = linearized(xnext) do x
        separation(obstacle, transform(x, point, obstacle.contact_face.outward_normal.frame))
    end

    D_transpose_times_v = map(D) do d
        linearized(xnext) do x
            v = point_velocity(twist_wrt_world(x, body), transform_to_root(linearization_state(xnext), point.frame) * point)
            dot(transform_to_root(linearization_state(xnext), d.frame) * d,
                v)
        end
    end

    add_contact_constraints(model, obstacle.μ, β, λ, c_n, D_transpose_times_v, separation_from_obstacle)
    ContactResult(β, λ, c_n, point, obstacle, total_weight)
end

# function add_contact_constraints_nonsliding(model::Model, point, obstacle, β, λ, c_n, D, separation_from_obstacle, contact_velocity)
#     D_transpose_times_v = [dot(d, contact_velocity) for d in D]

#     @constraints model begin
#         separation_from_obstacle >= 0 # (7)
#         # λ .+ D_transpose_times_v .>= 0 # (8)
#         obstacle.μ * c_n .- sum(β) >= 0 # (9)
#     end

#     @disjunction(model,
#                  (&)(separation_from_obstacle == 0,
#                      [@?(contact_velocity.v[i] == 0) for i in 1:3]...),
#                  c_n == 0)
#     ContactResult(β, λ, c_n, point, obstacle)
# end

# function add_contact_constraints_sticking(model::Model, point, obstacle, β, λ, c_n, D, separation_from_obstacle, contact_velocity)

#     @constraints model begin
#         separation_from_obstacle <= 1e-3 # (7)
#         separation_from_obstacle >= -1e-3
#         contact_velocity.v .== 0
#         λ == 0 # (8)
#         obstacle.μ * c_n .- sum(β) >= 0 # (9)
#     end

#     ContactResult(β, λ, c_n, point, obstacle)
# end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

@testset "graph.plainhypergraph" begin
    g = plainhypergraph(6, [Set([1, 2, 3, 4]), Set([1, 4, 5]), Set([6])])
    @test issetequal(neighbors(g, 4)[2], [1, 4, 5])
    @test issetequal(getedge(g, 1), 1:4)
    @test hasedge(g, 2, 3)
    @test !hasedge(g, 2, 5)
    @test_throws ErrorException nodeattr(g, 1)
    @test_throws ErrorException edgeattr(g, 1)
    @test_throws ErrorException edgeattr(g, 3, 4)

    @test issetequal(nodeset(g), 1:6)
    @test issetequal(edgeset(g), 1:3)
    @test_throws ErrorException nodeattrs(g)
    @test_throws ErrorException edgeattrs(g)

    @test nodecount(g) == 6
    @test edgecount(g) == 3
    @test neighborcount(g, 2) == 1
    @test degree(g, 3) == 1

    addnode!(g)
    @test degree(g, 7) == 0
    addedge!(g, Set([4, 6]))
    @test degree(g, 4) == 3
    @test_throws ErrorException nodeattrtype(g)
    @test_throws ErrorException edgeattrtype(g)

    emp = plainhypergraph()
    @test nodecount(emp) == 0

end # graph.plaingraph


text:
@testset "Constrained" begin
    # EI Min
    f = EfficientGlobalOptimization.PressureVesselDesign()
    yₒ, xₒ = optimum(f)
    lb = lowerbounds(f)
    ub = upperbounds(f)
    cons = constrains(f)
    xs = sampling(lb, ub, 9, cons...)
    x0 = Vector{Float64}[]
    y0 = Vector{Float64}[]
    maxg(x) =  maximum(g->g(x) ,cons)
    for x in eachcol(xs)
        if maxg(x) > eps()
            continue
        end
        push!(x0, x)
        push!(y0, [f(x)])
    end
    @info "$(length(y0)) initial samples"
    krg = Kriging(hcat(x0...), hcat(y0...), MeanConst(5e3), Mat32Ard(zeros(4), 3.0))
    kbd = (fill(-10.0, 5), fill(10.0, 5))
    mbd = ([1e2], [1e4])
    tuner = KrigingTuner(kbd, mbd; global_every = 8, local_every = 2)
    ego = EGO(krg, tuner, lb, ub; constrains=cons)
    try
        for i = 1:64
            x, fx = acquire(ego)
            y = f(x)
            append!(ego, x, y)
        end
    catch e
        if isa(e, NoFeasibleInfill)
        else
            rethrow(e)
        end
    end
    yo, xo = optimum(ego)
    @test yo[1] ≈ f(xo)
    @test yo[1]/yₒ < 1.05

end


text:
module TestSubDataFrame
    using Test, DataFrames

    @testset "copy - SubDataFrame" begin
        df = DataFrame(x = 1:10, y = 1.0:10.0)
        sdf = view(df, 1:2, 1:1)
        @test sdf isa SubDataFrame
        @test copy(sdf) isa DataFrame
        @test sdf == copy(sdf)
    end

    @testset "view -- DataFrame" begin
        df = DataFrame(x = 1:10, y = 1.0:10.0)
        @test view(df, 1) == head(df, 1)
        @test view(df, UInt(1)) == head(df, 1)
        @test view(df, BigInt(1)) == head(df, 1)
        @test view(df, 1:2) == head(df, 2)
        @test view(df, vcat(trues(2), falses(8))) == head(df, 2)
        @test view(df, [1, 2]) == head(df, 2)
        @test view(df, 1, :x) == head(df[[:x]], 1)
        @test view(df, 1:2, :x) == head(df[[:x]], 2)
        @test view(df, vcat(trues(2), falses(8)), :x) == head(df[[:x]], 2)
        @test view(df, [1, 2], :x) == head(df[[:x]], 2)
        @test view(df, 1, 1) == head(df[[:x]], 1)
        @test view(df, 1:2, 1) == head(df[[:x]], 2)
        @test view(df, vcat(trues(2), falses(8)), 1) == head(df[[:x]], 2)
        @test view(df, [1, 2], 1) == head(df[[:x]], 2)
        @test view(df, 1, [:x, :y]) == head(df, 1)
        @test view(df, 1:2, [:x, :y]) == head(df, 2)
        @test view(df, vcat(trues(2), falses(8)), [:x, :y]) == head(df, 2)
        @test view(df, [1, 2], [:x, :y]) == head(df, 2)
        @test view(df, 1, [1, 2]) == head(df, 1)
        @test view(df, 1:2, [1, 2]) == head(df, 2)
        @test view(df, vcat(trues(2), falses(8)), [1, 2]) == head(df, 2)
        @test view(df, [1, 2], [1, 2]) == head(df, 2)
        @test view(df, 1, trues(2)) == head(df, 1)
        @test view(df, 1:2, trues(2)) == head(df, 2)
        @test view(df, vcat(trues(2), falses(8)), trues(2)) == head(df, 2)
        @test view(df, [1, 2], trues(2)) == head(df, 2)
        @test view(df, Integer[1, 2]) == head(df, 2)
        @test view(df, UInt[1, 2]) == head(df, 2)
        @test view(df, BigInt[1, 2]) == head(df, 2)
        @test view(df, Union{Int, Missing}[1, 2]) == head(df, 2)
        @test view(df, Union{Integer, Missing}[1, 2]) == head(df, 2)
        @test view(df, Union{UInt, Missing}[1, 2]) == head(df, 2)
        @test view(df, Union{BigInt, Missing}[1, 2]) == head(df, 2)
        @test view(df, :) == df
        @test view(df, :, :) == df
        @test view(df, 1, :) == head(df, 1)
        @test view(df, :, 1) == df[:, [1]]
        @test_throws MissingException view(df, [missing, 1])
    end

    @testset "view -- SubDataFrame" begin
        df = view(DataFrame(x = 1:10, y = 1.0:10.0), 1:10)
        @test view(df, 1) == head(df, 1)
        @test view(df, UInt(1)) == head(df, 1)
        @test view(df, BigInt(1)) == head(df, 1)
        @test view(df, 1:2) == head(df, 2)
        @test view(df, vcat(trues(2), falses(8))) == head(df, 2)
        @test view(df, [1, 2]) == head(df, 2)
        @test view(df, 1, :x) == head(df[[:x]], 1)
        @test view(df, 1:2, :x) == head(df[[:x]], 2)
        @test view(df, vcat(trues(2), falses(8)), :x) == head(df[[:x]], 2)
        @test view(df, [1, 2], :x) == head(df[[:x]], 2)
        @test view(df, 1, 1) == head(df[[:x]], 1)
        @test view(df, 1:2, 1) == head(df[[:x]], 2)
        @test view(df, vcat(trues(2), falses(8)), 1) == head(df[[:x]], 2)
        @test view(df, [1, 2], 1) == head(df[[:x]], 2)
        @test view(df, 1, [:x, :y]) == head(df, 1)
        @test view(df, 1:2, [:x, :y]) == head(df, 2)
        @test view(df, vcat(trues(2), falses(8)), [:x, :y]) == head(df, 2)
        @test view(df, [1, 2], [:x, :y]) == head(df, 2)
        @test view(df, 1, [1, 2]) == head(df, 1)
        @test view(df, 1:2, [1, 2]) == head(df, 2)
        @test view(df, vcat(trues(2), falses(8)), [1, 2]) == head(df, 2)
        @test view(df, [1, 2], [1, 2]) == head(df, 2)
        @test view(df, 1, trues(2)) == head(df, 1)
        @test view(df, 1:2, trues(2)) == head(df, 2)
        @test view(df, vcat(trues(2), falses(8)), trues(2)) == head(df, 2)
        @test view(df, [1, 2], trues(2)) == head(df, 2)
        @test view(df, Integer[1, 2]) == head(df, 2)
        @test view(df, UInt[1, 2]) == head(df, 2)
        @test view(df, BigInt[1, 2]) == head(df, 2)
        @test view(df, Union{Int, Missing}[1, 2]) == head(df, 2)
        @test view(df, Union{Integer, Missing}[1, 2]) == head(df, 2)
        @test view(df, Union{UInt, Missing}[1, 2]) == head(df, 2)
        @test view(df, Union{BigInt, Missing}[1, 2]) == head(df, 2)
        @test view(df, :) == df
        @test view(df, :, :) == df
        @test view(df, 1, :) == head(df, 1)
        @test view(df, :, 1) == df[:, [1]]
        @test_throws MissingException view(df, [missing, 1])
    end

    @testset "getproperty, setproperty! and propertynames" begin
        x = collect(1:10)
        y = collect(1.0:10.0)
        df = view(DataFrame(x = x, y = y), 2:6)

        @test Base.propertynames(df) == names(df)

        @test df.x == 2:6
        @test df.y == 2:6
        @test_throws KeyError df.z

        df.x = 1:5
        @test df.x == 1:5
        @test x == [1; 1:5; 7:10]
        df.y = 1
        @test df.y == [1, 1, 1, 1, 1]
        @test y == [1; 1; 1; 1; 1; 1; 7:10]
        @test_throws ErrorException df.z = 1:5
        @test_throws ErrorException df.z = 1
    end
end


text:
export CuStream

const CU_STREAM_DEFAULT = Cint(0)
const CU_STREAM_NON_BLOCKING = Cint(1)

mutable struct CuStream
    ptr::Ptr{Cvoid}

    function CuStream(flags::Cint=CU_STREAM_DEFAULT)
        ref = Ref{Ptr{Cvoid}}()
        @apicall :cuStreamCreate (Ptr{Ptr{Cvoid}},Cuint) ref flags
        s = new(ref[], getcontext())
        finalizer(s, destroy)
        s
    end
end

Base.unsafe_convert(::Type{Ptr{Cvoid}}, s::CuStream) = s.ptr

function destroy(s::CuStream)
    setcontext(s.ctx) do
        @apicall :cuStreamDestroy (Ptr{Cvoid},) s
    end
end

synchronize(s::CuStream) = @apicall :cuStreamSynchronize (Ptr{Cvoid},) s


text:
# The @benchmarkable macro combines a tuple of expressions into a function
# that implements the "benchmarkable" protocol. We assume that precompilation,
# when required, will be handled elsewhere; as such, we make no extra calls to
# the core expression that we want to benchmark.
#
# Arguments:
#
#     name: The name of the function that will be generated.
#
#     setup: An expression that will be executed before the core
#         expression starts executing.
#
#     core: The core expression that will be executed.
#
#     teardown: An expression that will be executed after the core
#         expression finishes executing.
#
# Generates:
#
#     A function expression that defines a new function of three-arguments,
#     which are:
#
#     (1): s::Samples: A Samples object in which benchmark data will be stored.
#
#     (2): n_samples::Integer: The number of samples that will be gathered.
#
#     (3): n_evals::Integer: The number of times the core expression will be
#         evaluated per sample.

macro benchmarkable(name, setup, core, teardown)
    # We only support function calls to be passed in `core`, but some syntaxes
    # have special parsing that will lower to a function call upon expansion
    # (e.g., A[i] or user macros). The tricky thing is that keyword functions
    # *also* have a special lowering step that we don't want (they are harder to
    # work with and would obscure some of the overhead of keyword arguments). So
    # we only expand the passed expression if it is not already a function call.
    expr = (core.head == :call) ? core : expand(core)
    expr.head == :call || throw(ArgumentError("expression to benchmark must be a function call"))
    f = expr.args[1]
    fargs = expr.args[2:end]
    nargs = length(expr.args)-1

    # Pull out the arguments -- both positional and keywords
    userargs = Any[]  # The actual expressions the user wrote
    args = Symbol[gensym("arg_$i") for i in 1:nargs] # The local argument names
    posargs = Symbol[] # The names that are positional arguments
    kws = Expr[]       # Names that are used in keyword arguments
    for i in 1:nargs
        if isa(fargs[i], Expr) && fargs[i].head == :kw
            push!(kws, Expr(:kw, fargs[i].args[1], args[i]))
            push!(userargs, fargs[i].args[2])
        else
            push!(posargs, args[i])
            push!(userargs, fargs[i])
        end
    end

    benchfn = gensym("bench")
    innerfn = gensym("inner")

    # Strategy: we create *three* functions:
    # * The outermost function is the entry point. It's simply a closure around
    #   the expressions the user passed in `core` as the arguments to the
    #   benchmarked function. This allows the arguments to be considered setup,
    #   which are evaluated in the correct scope. However, that means that
    #   within this outermost function, the arguments probably aren't
    #   concretely typed. This means that if we were to run the benchmarking
    #   function in this outermost function, we'd end up benchmarking dynamic
    #   dispatch most of the time.  So we introduce a function barrier here.
    # * The second level (`benchfn`) is the benchmarking loop.  Here is where
    #   the real work gets done.  However, if we were to call the benchmarked
    #   function directly here, it might get inlined.  And if it gets inlined,
    #   then LLVM can use optimizations that interact with the test loop itself.
    #   No longer are we simply testing the benchmarked function; we are testing
    #   the benchmark loops.  So in order to circumvent this, we introduce a
    #   third function that is explicitly marked `@noinline`
    # * It is within this third, `inner` function that we call the user's
    #   function that they want to benchmark. This means that all timings will
    #   include the overhead of at least one function call. But it also means
    #   that we can prevent LLVM from doing optimizations that are related to
    #   the benchmarking itself: it must always call the inner function in the
    #   benchmarking function (since at the mid-level it doesn't know what that
    #   function might do), and within the inner function it can only eliminate
    #   code that's unrelated to the return value (since it doesn't know what
    #   the caller might do).
    quote
        function $(esc(name))(
                s::Samples,
                n_samples::Integer,
                evaluations::Integer,
            )
            gc()
            $(benchfn)(s, n_samples, evaluations, $(map(esc, userargs)...))
        end
        function $(benchfn)(
                s::Samples,
                n_samples::Integer,
                evaluations::Integer,
                $(args...)
            )
            # Execute the setup expression exactly once
            $(esc(setup))

            # Generate n_samples by evaluating the core
            for _ in 1:n_samples
                # Store pre-evaluation state information
                stats = Base.gc_num()
                time_before = time_ns()

                # Evaluate the core expression n_evals times.
                for _ in 1:evaluations
                    out = $(esc(innerfn))($(args...))
                end

                # get time before comparing GC info
                elapsed_time = time_ns() - time_before

                # Compare post-evaluation state with pre-evaluation state.
                diff = Base.GC_Diff(Base.gc_num(), stats)
                bytes = diff.allocd
                allocs = diff.malloc + diff.realloc + diff.poolalloc + diff.bigalloc
                gc_time = diff.total_time

                # Append data for this sample to the Samples objects.
                push!(s, evaluations, elapsed_time, gc_time, bytes, allocs)
            end

            # Execute the teardown expression exactly once
            $(esc(teardown))

            # The caller receives all data via the mutated Results object.
            return
        end
        $(esc(:(@noinline function $innerfn($(args...))
            $f($(posargs...), $(kws...))
        end)))

        # "return" the outermost entry point as the final expression
        $(esc(name))
    end
end


text:
using GridWorlds
using Test
using Random
using ReinforcementLearningBase

ENVS = [EmptyGridWorld, FourRooms, GoToDoor, DoorKey, CollectGems, DynamicObstacles]
ENVS_RLBASE = [EmptyGridWorld, CollectGems]
ACTIONS = [TURN_LEFT, TURN_RIGHT, MOVE_FORWARD]

@testset "GridWorlds.jl" begin
    for env in ENVS
        @testset "$(env)" begin
            w = env()
            @test typeof(w.agent_pos) == CartesianIndex{2}
            @test typeof(w.agent.dir) <: LRUD
            @test size(w.world.world, 1) == length(w.world.objects)
            @test 1 ≤ w.agent_pos[1] ≤ size(w.world.world, 2)
            @test 1 ≤ w.agent_pos[2] ≤ size(w.world.world, 3)

            for _=1:1000
                w = w(rand(ACTIONS))
                @test 1 ≤ w.agent_pos[1] ≤ size(w.world.world, 2)
                @test 1 ≤ w.agent_pos[2] ≤ size(w.world.world, 3)
                @test w.world[WALL, w.agent_pos] == false
                view = get_agent_view(w)
                @test typeof(view) <: BitArray{3}
                @test size(view,1) == length(w.world.objects)
                @test size(view,2) == 7
                @test size(view,3) == 7
            end
        end
    end

    for env in ENVS_RLBASE
        @testset "$(env) RLBase API" begin
            w = env()

            @test get_reward(w) == 0.0
            @test get_terminal(w) == false
            @test get_state(w) == get_agent_view(w)

            π = RandomPolicy(w)
            total_reward = 0.0

            while !get_terminal(w)
                action = π(w)
                w(action)
                total_reward += get_reward(w)
            end

            if env == CollectGems
                @test total_reward == w.num_gem_init * w.gem_reward
            elseif env == EmptyGridWorld
                @test total_reward == w.goal_reward
            end
        end
    end
end


text:
# # Example: 1D Bose-Hubbard Model

# This is an example calculation finding the ground state of
# a 1D Bose-Hubbard chain with 6 particles in 6 lattice sites.
# The Julia run-able script is in [`scripts/BHM-example.jl`](../../scripts/BHM-example.jl).

# Firstly, we load all needed modules.
# `Rimu` for FCIQMC calculation;

using Rimu

# Now we define the physical problem:
# Setting the number of lattice sites `m = 6`;
# and the number of particles `n = 6`:
m = n = 6
# Generating a configuration that particles are evenly distributed:
aIni = near_uniform(BoseFS{n,m})
# where `BoseFS` is used to create a bosonic system.
# The Hamiltonian is defined based on the configuration `aIni`,
# with additional onsite interaction strength `u = 6.0`
# and the hopping strength `t = 1.0`:
Ĥ = HubbardReal1D(aIni; u = 6.0, t = 1.0)

# Now let's setup the Monte Carlo settings.
# The number of walkers to use in this Monte Carlo run:
targetwalkers = 1_000
# The number of time steps before doing statistics,
# i.e. letting the walkers to sample Hilbert and to equilibrate:
steps_equilibrate = 1_000
# And the number of time steps used for getting statistics,
# e.g. time-average of shift, projected energy, walker numbers, etc.:
steps_measure = 1_000

# Set the size of a time step
dτ = 0.001
# and we report QMC data every k-th step,
# setting `k = 1` means we record QMC data every step:
k = 1

# Now we prepare initial state and allocate memory.
# The initial address is defined above as `aIni = near_uniform(Ĥ)`.
# Define the initial number of walkers per rank:
nIni = 1
# Putting the `nIni` number of walkers into the initial address `aIni`
svec = DVec(aIni => nIni)
# Let's plant a seed for the random number generator to get consistent result:
Rimu.ConsistentRNG.seedCRNG!(17)

# Now let's setup all the FCIQMC strategies.

# Passing dτ and total number of time steps into params:
params = RunTillLastStep(dτ = dτ, laststep = steps_equilibrate + steps_measure)
# Strategy for updating the shift:
s_strat = DoubleLogUpdate(targetwalkers = targetwalkers, ζ = 0.08)
# Strategy for reporting info:
r_strat = ReportDFAndInfo(k = k, i = 100)
# Strategy for updating dτ:
t_strat = ConstantTimeStep()
# set up the calculation and reporting of the projected energy
# in this case we are projecting onto the starting vector,
# which contains a single configuration
post_step = ProjectedEnergy(Ĥ, copy(svec))

# Print out info about what we are doing:
println("Finding ground state for:")
println(Ĥ)
println("Strategies for run:")
println(params, s_strat)
println(t_strat)

# Finally, we can start the main FCIQMC loop:
df, state = lomc!(Ĥ,svec;
            params = params,
            laststep = steps_equilibrate + steps_measure,
            s_strat = s_strat,
            r_strat = r_strat,
            τ_strat = t_strat,
            post_step = post_step,
            threading = false, # only for reproducible runs
)
println("Writing data to disk...")

# Saving output data stored in `df` into a `.arrow` file which can be read in later:
save_df("fciqmcdata.arrow", df)

# Now let's look at the calculated energy from the shift:
# Loading the equilibrated data:
qmcdata = last(df,steps_measure)
using Rimu.StatsTools

# For the shift, it's easy to use `mean_and_se` from `Rimu.StatsTools`
(qmcShift,qmcShiftErr) = mean_and_se(qmcdata.shift)

# For the projected energy, it a bit more complicated as it's a ratio of two means:
r = ratio_of_means(qmcdata.hproj,qmcdata.vproj)
rwe = ratio_with_errs(r)

# Here we use the 95% CI for the lower and upper error bars:
(eProj,eProjErrLower,eProjErrUpper) = (rwe.ratio, rwe.err2_l, rwe.err2_u)

println("Energy from $steps_measure steps with $targetwalkers walkers:
         Shift: $qmcShift ± $qmcShiftErr;
         Projected Energy: $eProj ± ($eProjErrLower, $eProjErrUpper)")

# Finished !
println("Finished!")


text:
# Statistical utilities.

# Ambegaokar, V., & Troyer, M. (2010). Estimating errors reliably in Monte
# Carlo simulations of the Ehrenfest model. American Journal of Physics, 78(2),
# 150-157.

"""
    stderr_bin(dt::Float64, xs::AbstractVector)

Compute the standard error and autocorrelation time of `xs` by binning.

Stops when there are too few bins left. Returns the largest error attained and
the corresponding autocorrelation time in units of `dt`.
"""
function stderr_bin(dt::Float64, xs::AbstractVector)
    stderr_full = std(xs) / sqrt(length(xs))
    stderr_max = stderr_full

    while length(xs) >= 60
        xs = [0.5 * (xs[i] + xs[i+1]) for i in (1+length(xs)%2):2:length(xs)]
        stderr_max = max(stderr_max, std(xs) / sqrt(length(xs)))
    end

    autocorrelation_steps = 0.5 * ((stderr_max / stderr_full)^2 - 1)

    stderr_max, autocorrelation_steps * dt
end

"""
    aggregate(dt::Float64, xs::AbstractVector)

Perform a statistical analysis of `xs` for time step `dt`.
"""
function aggregate(dt::Float64, xs::AbstractVector)
    mean(xs), stderr_bin(dt, xs)...
end


text:
# # Car
#
# function car_dynamics!(ẋ::AbstractVector,x::AbstractVector,u::AbstractVector)
#     ẋ[1] = u[1]*cos(x[3])
#     ẋ[2] = u[1]*sin(x[3])
#     ẋ[3] = u[2]
#     return nothing
# end
# n,m = 3,2
#
# car = Model(car_dynamics!,n,m)

@with_kw struct DubinsCar <: AbstractModel
    radius::Float64 = 0.175  # [m] radius of a Roomba
end
Base.size(::DubinsCar) = 3,2

function dynamics(::DubinsCar,x,u)
    ẋ = @SVector [u[1]*cos(x[3]),
                  u[1]*sin(x[3]),
                  u[2]]
end

Base.position(::DubinsCar, x::SVector) = @SVector [x[1], x[2], 0.0]
orientation(::DubinsCar, x::SVector) = expm(x[3]*@SVector [0,0,1.])

# using Plots
# function Plots.plot(model::DubinsCar, Z::Traj; kwargs...)
#     x = [z.z[1] for z in Z]
#     y = [z.z[2] for z in Z]
#     Plots.plot(x, y, xlabel="x position", ylabel="y position", markershape=:circle;
#         kwargs...)
# end
# function Plots.plot!(model::DubinsCar, Z::Traj; kwargs...)
#     x = [z.z[1] for z in Z]
#     y = [z.z[2] for z in Z]
#     Plots.plot!(x, y, xlabel="x position", ylabel="y position", markershape=:circle;
#         kwargs...)
# end


text:
function JuAFEM.MixedDofHandler(grid::AdaptiveGrid{dim,T}) where {dim,T}
    JuAFEM.MixedDofHandler{dim,T,typeof(grid)}(FieldHandler[], JuAFEM.CellVector(Int[],Int[],Int[]), JuAFEM.CellVector(Int[],Int[],Int[]), JuAFEM.CellVector(Vec{dim,T}[],Int[],Int[]), JuAFEM.ScalarWrapper(false), grid, JuAFEM.ScalarWrapper(-1))
end

function JuAFEM.__close!(dh::MixedDofHandler{dim,T,AdaptiveGrid{dim,T}}) where {dim, T}

    @assert !JuAFEM.isclosed(dh)
    field_names = JuAFEM.getfieldnames(dh)  # all the fields in the problem
    numfields =  length(field_names)

    # Create dicts that stores created dofs
    # Each key should uniquely identify the given type
    vertexdicts = [Dict{Int, Array{Int}}() for _ in 1:numfields]
    edgedicts = [Dict{Tuple{Int,Int}, Array{Int}}() for _ in 1:numfields]
    facedicts = [Dict{Tuple{Int,Int}, Array{Int}}() for _ in 1:numfields]
    celldicts = [Dict{Int, Array{Int}}() for _ in 1:numfields]

    # Set initial values
    nextdof = 1  # next free dof to distribute

    append!(dh.cell_dofs.offset, zeros(getncells(dh.grid)))
    append!(dh.cell_dofs.length, zeros(getncells(dh.grid)))

    @debug "\n\nCreating dofs\n"
    for fh in dh.fieldhandlers
        # sort the cellset since we want to loop through the cells in a fixed order
        cellnumbers = sort(collect(fh.cellset))
        nextdof = JuAFEM._close!(
            dh,
            cellnumbers,
            JuAFEM.getfieldnames(fh),
            JuAFEM.getfielddims(fh),
            JuAFEM.getfieldinterpolations(fh),
            nextdof,
            vertexdicts,
            edgedicts,
            facedicts,
            celldicts)
    end
    dh.ndofs[] = maximum(dh.cell_dofs.values)
    dh.closed[] = true

    #Create cell_nodes and cell_coords (similar to cell_dofs)
    push!(dh.cell_nodes.offset, 1)
    push!(dh.cell_coords.offset, 1)
    for cell in dh.grid.cells
        for nodeid in JuAFEM.vertices(cell, dh.grid)
            push!(dh.cell_nodes.values, nodeid)
            push!(dh.cell_coords.values, dh.grid.nodes[nodeid].x)
        end
        push!(dh.cell_nodes.offset, length(dh.cell_nodes.values)+1)
        push!(dh.cell_coords.offset, length(dh.cell_coords.values)+1)
        push!(dh.cell_nodes.length, length(cell.nodes))
        push!(dh.cell_coords.length, length(cell.nodes))
    end

    return dh, vertexdicts, edgedicts, facedicts

end

function JuAFEM._close!(dh::JuAFEM.MixedDofHandler{dim, T, AdaptiveGrid{dim, T}}, cellnumbers, field_names, field_dims, field_interpolations, nextdof, vertexdicts, edgedicts, facedicts, celldicts) where {dim, T}

    ip_infos = JuAFEM.InterpolationInfo[]
    for interpolation in field_interpolations
        ip_info = JuAFEM.InterpolationInfo(interpolation)
        # these are not implemented yet (or have not been tested)
        @assert(ip_info.nvertexdofs <= 1)
        @assert(ip_info.nedgedofs <= 1)
        @assert(ip_info.nfacedofs <= 1)
        @assert(ip_info.ncelldofs <= 1)  # not tested but probably works
        push!(ip_infos, ip_info)
    end

    # loop over all the cells, and distribute dofs for all the fields
    for ci in cellnumbers
        dh.cell_dofs.offset[ci] = length(dh.cell_dofs.values)+1

        cell = dh.grid.cells[ci]
        cell_dofs = Int[]  # list of global dofs for each cell
        @debug "Creating dofs for cell #$ci"

        for fi in 1:length(field_names)
            @debug "\tfield: $(field_names[fi])"
            ip_info = ip_infos[fi]

            if ip_info.nvertexdofs > 0
                vertices = JuAFEM.vertices(cell, dh.grid)
                nextdof = JuAFEM.add_vertex_dofs(cell_dofs, cell, vertexdicts[fi], field_dims[fi], ip_info.nvertexdofs, nextdof, vertices)
            end

            if ip_info.nedgedofs > 0 && dim == 3 #Edges only in 3d
                edges = JuAFEM.edges(cell, dh.grid)
                nextdof = JuAFEM.add_edge_dofs(cell_dofs, cell, edgedicts[fi], field_dims[fi], ip_info.nedgedofs, nextdof, edges)
            end

            if ip_info.nfacedofs > 0 && (ip_info.dim == dim)
                faces = JuAFEM.faces(cell, dh.grid)
                nextdof = JuAFEM.add_face_dofs(cell_dofs, cell, facedicts[fi], field_dims[fi], ip_info.nfacedofs, nextdof, faces)
            end

            if ip_info.ncelldofs > 0
                nextdof = JuAFEM.add_cell_dofs(cell_dofs, ci, celldicts[fi], field_dims[fi], ip_info.ncelldofs, nextdof)
            end

        end
        # after done creating dofs for the cell, push them to the global list
        push!(dh.cell_dofs.values, cell_dofs...)
        dh.cell_dofs.length[ci] = length(cell_dofs)

        @debug "Dofs for cell #$ci:\n\t$cell_dofs"
    end # cell loop
    return nextdof
end

function JuAFEM.add_vertex_dofs(cell_dofs, cell, vertexdict, field_dim, nvertexdofs, nextdof, vertices)
    for vertex in vertices
        @debug "\tvertex #$vertex"
        nextdof, dofs = JuAFEM.get_or_create_dofs!(nextdof, field_dim, dict=vertexdict, key=vertex)
        push!(cell_dofs, dofs...)
    end
    return nextdof
end

function JuAFEM.add_face_dofs(cell_dofs, cell, facedict, field_dim, nfacedofs, nextdof, faces)
    @assert nfacedofs == 1 "Currently only supports interpolations with nfacedofs = 1"

    for face in faces
        sface = JuAFEM.sortface(face)
        @debug "\tface #$sface"
        nextdof, dofs = JuAFEM.get_or_create_dofs!(nextdof, field_dim, dict=facedict, key=sface)
        push!(cell_dofs, dofs...)
    end
    return nextdof
end

function JuAFEM.add_edge_dofs(cell_dofs, cell, edgedict, field_dim, nedgedofs, nextdof, edges)
    @assert nedgedofs == 1 "Currently only supports interpolations with nedgedofs = 1"
    for edge in edges
        sedge, dir = JuAFEM.sortedge(edge)
        @debug "\tedge #$sedge"
        nextdof, dofs = JuAFEM.get_or_create_dofs!(nextdof, field_dim, dict=edgedict, key=sedge)
        push!(cell_dofs, dofs...)
    end
    return nextdof
end

##not sure if needed
#function add_cell_dofs(cell_dofs, cell, celldict, field_dim, ncelldofs, nextdof)
#    for celldof in 1:ncelldofs
#        @debug "\tcell #$cell"
#        nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=celldict, key=cell)
#        push!(cell_dofs, dofs...)
#    end
#    return nextdof
#end

function JuAFEM.add!(ch::ConstraintHandler, dbc::Dirichlet, ncdofs::Vector)
    JuAFEM.dbc_check(ch, dbc)
    field_idx = JuAFEM.find_field(ch.dh, dbc.field_name)
    # Extract stuff for the field
    interpolation = JuAFEM.getfieldinterpolation(ch.dh, field_idx)#ch.dh.field_interpolations[field_idx]
    field_dim = JuAFEM.getfielddim(ch.dh, field_idx)#ch.dh.field_dims[field_idx] # TODO: I think we don't need to extract these here ...
    bcvalue = JuAFEM.getbcvalue(ch.dh, field_idx)
    JuAFEM._add!(ch, dbc, dbc.faces, interpolation, field_dim, JuAFEM.field_offset(ch.dh, dbc.field_name), bcvalue, ncdofs)
    return ch
end

function JuAFEM._add!(ch::ConstraintHandler, dbc::Dirichlet, bcfaces::Set{Tuple{Int,Int}}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::JuAFEM.BCValues, ncdofs::Vector)
    # calculate which local dof index live on each face
    # face `i` have dofs `local_face_dofs[local_face_dofs_offset[i]:local_face_dofs_offset[i+1]-1]
    local_face_dofs = Int[]
    local_face_dofs_offset = Int[1]
    for (i, face) in enumerate(JuAFEM.faces(interpolation))
        for fdof in face, d in 1:field_dim
            if d ∈ dbc.components # skip unless this component should be constrained
                push!(local_face_dofs, (fdof-1)*field_dim + d + offset)
            end
        end
        push!(local_face_dofs_offset, length(local_face_dofs) + 1)
    end
    JuAFEM.copy!!(dbc.local_face_dofs, local_face_dofs)
    JuAFEM.copy!!(dbc.local_face_dofs_offset, local_face_dofs_offset)

    # loop over all the faces in the set and add the global dofs to `constrained_dofs`
    constrained_dofs = Int[]
    #_celldofs = fill(0, ndofs_per_cell(ch.dh))
    for (cellidx, faceidx) in bcfaces
        _celldofs = fill(0, ndofs_per_cell(ch.dh, cellidx))
        celldofs!(_celldofs, ch.dh, cellidx) # extract the dofs for this cell
        r = local_face_dofs_offset[faceidx]:(local_face_dofs_offset[faceidx+1]-1)
        dofs = _celldofs[local_face_dofs[r]]
        for dof in dofs
            append!(constrained_dofs, dof - sum(ncdofs .<= dof )) # TODO: for-loop over r and simply push! to ch.prescribed_dofs
        end
        @debug println("adding dofs $(_celldofs[local_face_dofs[r]]) to dbc")
    end

    # save it to the ConstraintHandler
    push!(ch.dbcs, dbc)
    push!(ch.bcvalues, bcvalue)
    append!(ch.prescribed_dofs, constrained_dofs)
end

function JuAFEM.close!(ch::ConstraintHandler, ncdofs::Vector)
    fdofs = setdiff(1:ndofs(ch.dh), ch.prescribed_dofs)
    JuAFEM.copy!!(ch.free_dofs, fdofs)
    JuAFEM.copy!!(ch.prescribed_dofs, unique(ch.prescribed_dofs)) # for v0.7: unique!(ch.prescribed_dofs)
    sort!(ch.prescribed_dofs) # YOLO
    fill!(resize!(ch.values, length(ch.prescribed_dofs)), NaN)
    for i in 1:length(ch.prescribed_dofs)
        key = ch.prescribed_dofs[i]
        #ch.dofmapping[key-sum(ncdofs .<= key)] = i
        ch.dofmapping[key] = i
    end
    ch.closed[] = true
    return ch
end

# Updates the DBC's to the current time `time`
function JuAFEM.update!(ch::ConstraintHandler,ncdofs::Vector, time::Real=0.0)
    @assert ch.closed[]
    for (i,dbc) in enumerate(ch.dbcs)
        field_idx = JuAFEM.find_field(ch.dh, dbc.field_name)
        # Function barrier
        JuAFEM._update!(ch.values, dbc.f, dbc.faces, dbc.field_name, dbc.local_face_dofs, dbc.local_face_dofs_offset,
                 dbc.components, ch.dh, ch.bcvalues[i], ch.dofmapping, convert(Float64, time), ncdofs)
    end
end

# for faces
function JuAFEM._update!(values::Vector{Float64}, f::Function, faces::Set{Tuple{Int,Int}}, field::Symbol, local_face_dofs::Vector{Int}, local_face_dofs_offset::Vector{Int},
                  components::Vector{Int}, dh::JuAFEM.AbstractDofHandler, facevalues::JuAFEM.BCValues,
                  dofmapping::Dict{Int,Int}, time::T, ncdofs::Vector) where {T}

    dim = JuAFEM.getdim(dh.grid)
    _tmp_cellid = first(faces)[1]

    N = JuAFEM.nnodes_per_cell(dh.grid, _tmp_cellid)
    xh = zeros(Vec{dim, T}, N) # pre-allocate
    _celldofs = fill(0, ndofs_per_cell(dh, _tmp_cellid))

    for (cellidx, faceidx) in faces
        JuAFEM.cellcoords!(xh, dh, cellidx)
        JuAFEM.celldofs!(_celldofs, dh, cellidx) # update global dofs for this cell

        # no need to reinit!, enough to update current_face since we only need geometric shape functions M
        facevalues.current_face[] = faceidx

        # local dof-range for this face
        r = local_face_dofs_offset[faceidx]:(local_face_dofs_offset[faceidx+1]-1)
        counter = 1

        for location in 1:getnquadpoints(facevalues)
            x = spatial_coordinate(facevalues, location, xh)
            bc_value = f(x, time)
            @assert length(bc_value) == length(components)

            for i in 1:length(components)
                # find the global dof
                globaldof = _celldofs[local_face_dofs[r[counter]]]
                globaldof -= sum(ncdofs .<= globaldof)

                counter += 1

                dbc_index = dofmapping[globaldof]
                values[dbc_index] = bc_value[i]
                @debug println("prescribing value $(bc_value[i]) on global dof $(globaldof)")
            end
        end
    end
end

function facedofs(dh::JuAFEM.AbstractDofHandler, ele::Int, localface::Int)
    eledofs = celldofs(dh,ele)
    return eledofs[[JuAFEM.faces(Lagrange{2,RefCube,1}())[localface]...]], eledofs ## TODO decluster this
end

function find_local_coordinate(interpolation, cell_coordinates, global_coordinate)
    dim = length(global_coordinate)
    local_guess = zero(Vec{dim})
    n_basefuncs = getnbasefunctions(interpolation)
    max_iters = 10
    tol_norm = 1e-10
    for iter in 1:10
        if iter == max_iters
            error("did not find a local coordinate")
        end
        N = JuAFEM.value(interpolation, local_guess)

        global_guess = zero(Vec{dim})
        for j in 1:n_basefuncs
            global_guess += N[j] * cell_coordinates[j]
        end
        residual = global_guess - global_coordinate
        if norm(residual) <= tol_norm
            break
        end
        dNdξ = JuAFEM.derivative(interpolation, local_guess)
        J = zero(Tensor{2, 2})
        for j in 1:n_basefuncs
            J += cell_coordinates[j] ⊗ dNdξ[j]
        end
        local_guess -= inv(J) ⋅ residual
    end
    return local_guess
end


text:
using ReverseDiff: GradientTape, GradientConfig, gradient, gradient!, compile
using Turing: _hmc_step
using Turing
using HDF5, JLD
sv_data = load(TPATH*"/example-models/nips-2017/sv-data.jld.data")["data"]

@model sv_model(T, y) = begin
    ϕ ~ Uniform(-1, 1)
    σ ~ Truncated(Cauchy(0,5), 0, +Inf)
    μ ~ Cauchy(0, 10)
    # h = tzeros(Real, T)
    h = Vector{Real}(T)
    h[1] ~ Normal(μ, σ / sqrt(1 - ϕ^2))
    y[1] ~ Normal(0, exp.(h[1] / 2))
    for t = 2:T
      h[t] ~ Normal(μ + ϕ * (h[t-1] - μ) , σ)
      y[t] ~ Normal(0, exp.(h[t] / 2))
    end
  end

mf = sv_model(data=sv_data[1])
chain_nuts = sample(mf, HMC(0.05, 10), 2000)

println("mean of m: $(mean(chn[1000:end, :μ, :].value))")

# θ_dim = 1
# function lj_func(θ)
#   _lj = zero(Real)

#   s = 1

#   m = θ[1]
#   _lj += logpdf(Normal(0, sqrt(s)), m)

#   _lj += logpdf(Normal(m, sqrt(s)), 2.0)
#   _lj += logpdf(Normal(m, sqrt(s)), 2.5)

#   return _lj
# end

# neg_lj_func(θ) = -lj_func(θ)
# const f_tape = GradientTape(neg_lj_func, randn(θ_dim))
# const compiled_f_tape = compile(f_tape)

# function grad_func(θ)

#   inputs = θ
#   results = similar(θ)
#   all_results = DiffResults.GradientResult(results)

#   gradient!(all_results, compiled_f_tape, inputs)

#   neg_lj = all_results.value
#   grad, = all_results.derivs

#   return -neg_lj, grad

# end

# std = ones(θ_dim)
# θ = randn(θ_dim)
# lj = lj_func(θ)

# chn = []
# accept_num = 1

# totla_num = 5000
# for iter = 1:totla_num
#   push!(chn, θ)
#   θ, lj, is_accept, τ_valid, α = _hmc_step(θ, lj, lj_func, grad_func, 5, 0.05, std)
#   accept_num += is_accept
#   # if (iter % 50 == 0) println(θ) end
# end

# @show mean(chn), lj
# @show accept_num / totla_num


text:
module BitSAD

using LinearAlgebra
using DataStructures
using Random: MersenneTwister
using Setfield
using Ghost
using MacroTools: @capture
using LightGraphs, MetaGraphs
using Base: @kwdef
using Artifacts, LazyArtifacts

export AbstractBitstream
export SBit, SBitstream
# export DBit, DBitstream
export pos, neg
export decorrelate
export generate, generate!, estimate, estimate!, observe
# export SDM
export simulatable
export generatehw

include("tracing/trace.jl")

include("types/bitstream.jl")
include("types/sbitstream.jl")
# include("types/dbitstream.jl")
# include("modules/sdm.jl")

include("tracing/utilities.jl")
include("tracing/simulatable.jl")
include("tracing/netlist.jl")
include("tracing/circuit.jl")
include("tracing/transforms/insertrng.jl")
include("tracing/transforms/constantreplacement.jl")
include("tracing/transforms/constantreduction.jl")
include("tracing/hardware.jl")

include("hardware/utils.jl")
include("hardware/saddhandler.jl")
include("hardware/ssubhandler.jl")
include("hardware/smulthandler.jl")
include("hardware/sdivhandler.jl")
include("hardware/sfdivhandler.jl")
include("hardware/ssqrthandler.jl")
include("hardware/sl2normhandler.jl")
include("hardware/smaxhandler.jl")
include("hardware/transposehandler.jl")
include("hardware/reshapehandler.jl")
include("hardware/identityhandler.jl")
include("hardware/sbitstreamhandler.jl")
# include("hardware/daddhandler.jl")
# include("hardware/dsubhandler.jl")
# include("hardware/dmulthandler.jl")
# include("hardware/fxpaddhandler.jl")
# include("hardware/fxpsubhandler.jl")
# include("hardware/fxpmulthandler.jl")
# include("hardware/sdmhandler.jl")
# include("hardware/delaybufferhandler.jl")

function download_lib(dir = pwd())
    lib = artifact"verilog-stdlib"
    cp(lib, joinpath(dir, "verilog-stdlib"))
end

end # module


text:
"""
in the `positions` matrix, each column represents one particle and each row is one of
the coordinates (x, y, z, ...). If the input coordinates weren't in this format, it's
trivial to convert them to the desired matrix. A particle is counted if it's between `r`
and `r + dr` away from another particle.
"""
function radial_distribution(positions::Matrix{<:Real}, r::Real, dr::Real)
    inrange = 0
    for pos in eachcol(positions)
        distances = .√(sum(x -> x^2, positions .- pos, dims=1))
        inrange += count(x -> ((x <= r + dr) && (x >= r)), distances)
    end
    _, particle_count = size(positions)
    # average over all particles, then divide by the number of particles minus one because
    # the particle is never in range of "itself"
    return inrange / particle_count / (particle_count - 1)
end

max_stars_repo_name: slhshamloo/comp-phys

max_issues_repo_name: slhshamloo/comp-phys

max_forks_repo_name: slhshamloo/comp-phys


text:
##############################################################################
##
## Non Linear solver using Pseudo-Transient Continuation Method
##
##############################################################################

# Implicit time step
function implicit_timestep(G!, ypost, Δ; is_algebraic = fill(false, size(ypost)...), iterations = 100, verbose = true, method = :newton, autodiff = :forward, maxdist = sqrt(eps()), J0c = (nothing, nothing), y̲ = fill(-Inf, length(ypost)), ȳ = fill(Inf, length(ypost)), reformulation = :smooth)
    G_helper!(ydot, y) = (G!(ydot, y) ; ydot .-= .!is_algebraic .* (ypost .- y) ./ Δ)
    J0, colorvec = J0c
    if J0 == nothing
        result = nlsolve(G_helper!, ypost; iterations = iterations, show_trace = verbose, ftol = maxdist, method = method, autodiff = autodiff)
    else
        if autodiff == :forward
            jac_cache = ForwardColorJacCache(G_helper!, deepcopy(ypost); colorvec = colorvec, sparsity = J0)
            j_helper! = (J, y) -> forwarddiff_color_jacobian!(J, G_helper!, y, jac_cache)
        else
            j_helper! = (J, y) -> finite_difference_jacobian!(J, G_helper!, y; colorvec = colorvec)
        end
        if any(y̲ .!= -Inf) || any(ȳ .!= Inf)
            # using mcpsolve if lower/upper bounds are given
            result = mcpsolve(OnceDifferentiable(G_helper!, j_helper!, deepcopy(ypost), deepcopy(ypost), J0), y̲, ȳ, ypost; iterations = iterations, show_trace = verbose, ftol = maxdist, method = method, reformulation = reformulation)
        else
            result = nlsolve(OnceDifferentiable(G_helper!, j_helper!, deepcopy(ypost), deepcopy(ypost), J0), ypost; iterations = iterations, show_trace = verbose, ftol = maxdist, method = method)
        end
    end
    return result.zero, result.residual_norm
end

# Solve for steady state
function finiteschemesolve(G!, y0; Δ = 1.0, is_algebraic = fill(false, size(y0)...), iterations = 100, inner_iterations = 10, verbose = true, inner_verbose = false, method = :newton, autodiff = :forward, maxdist = sqrt(eps()), scale = 10.0, J0c = (nothing, nothing), minΔ = 1e-9, y̲ = fill(-Inf, length(y0)), ȳ = fill(Inf, length(y0)), reformulation = :smooth)
    ypost = y0
    ydot = zero(y0)
    G!(ydot, ypost)
    residual_norm = norm(ydot) / length(ydot)
    isnan(residual_norm) && throw("G! returns NaN with the initial value")
    if Δ == Inf
        ypost, residual_norm = implicit_timestep(G!, y0, Δ; is_algebraic = is_algebraic, verbose = verbose, iterations = iterations,  method = method, autodiff = autodiff, maxdist = maxdist, J0c = J0c, y̲ = y̲, ȳ = ȳ)
    else
        coef = 1.0
        oldresidual_norm = residual_norm
        iter = 0
        if verbose
            @printf "Iter   TimeStep   Residual\n"
            @printf "---- ---------- ----------\n"
        end
        while (iter < iterations) & (Δ >= minΔ) & (residual_norm > maxdist)
            iter += 1
            y, nlresidual_norm = implicit_timestep(G!, ypost, Δ; is_algebraic = is_algebraic, verbose = inner_verbose, iterations = inner_iterations, method = method, autodiff = autodiff, maxdist = maxdist, J0c = J0c, y̲ = y̲, ȳ = ȳ, reformulation = reformulation)
            G!(ydot, y)
            if any(y̲ .!= -Inf) || any(ȳ .!= Inf)
                mask = y̲ .+ eps() .<= y .<= ȳ .- eps() # only unconstrained ydot is relevant for residual_norm calculation
                residual_norm, oldresidual_norm = norm(ydot .* (mask))/sum(mask), residual_norm
            else
                residual_norm, oldresidual_norm = norm(ydot) / length(ydot), residual_norm
            end
            residual_norm = isnan(residual_norm) ? Inf : residual_norm
            if nlresidual_norm <= maxdist
                # if the implicit time step is correctly solved
                if verbose
                    @printf "%4d %8.4e %8.4e\n" iter Δ residual_norm
                end
                coef = (residual_norm <= oldresidual_norm) ? scale * coef : 1.0
                Δ = Δ * coef * oldresidual_norm / residual_norm
                ypost, y = y, ypost
            else
                # verbose && @show iter, Δ, NaN
                # if the implict time step is not solved
                # revert and diminish the time step
                coef = 1.0
                Δ = Δ / 10
                residual_norm = oldresidual_norm
            end
        end
    end
    verbose && ((residual_norm > maxdist) | (Δ < minΔ)) && @warn "Iteration did not converge"
    return ypost, residual_norm
end


text:
## Contains object for Delta Robot Params

"""
parameters for delta robot arm and delta robot models
"""
mutable struct DeltaRobotParams
    m1::Real # mass of first link (kg)
    m2::Real # mass of second link (kg)
    m3::Real # mass of end effector (kg)
    l1::Real # length of link 1 (m)
    l2::Real # length of link 2 (m)
    g::Real  # acceleration due to gravity (m/s^2)
    rb::Real # radius of base (m)
    rp::Real # radius of platform (m)
    d1::Real # damping between link 1 and base
    d2::Real # damping between link 1 and link 2 in phi2 direction
    d3::Real # damping between link 1 and link 2 in phi3 direction
    Kp::Real # proportional control constant
    Kd::Real # derivative control constant

    DeltaRobotParams(m1, m2, m3, l1, l2, g, rb=0.0, rp=0.0, d1=0.0, d2=0.0, d3=0.0, Kp=10.0, Kd=0.1) = new(m1, m2, m3, l1, l2, g, rb, rp, d1, d2, d3, Kp, Kd)
end

max_stars_repo_name: zborffs/Delta

max_issues_repo_name: zborffs/Delta

max_forks_repo_name: zborffs/Delta


text:
#=
regularized matrix regression problems

min 1/(2n) * ||Y - X * A||_fro^2 + lam_fro * ||A||_fro + lam_nuc * ||A||_nuc +
    lam_lass * ||A||_las + lam_glr * ||A||_glr + lamb_glc * ||A||_glc
- X is n x p
- Y is n x m
- A (variable) is p x m
- ||.||_fro is the Frobenius norm
- ||.||_nuc is the nuclear norm
- ||.||_las is the L1 norm
- ||.||_glr is the row group lasso penalty (L1 norm of row groups)
- ||.||_glc  is the column group lasso penalty (L1 norm of column groups)

TODO
- generalize for sparse Y,X but make sure qr factorization does not permute
=#

using SparseArrays

struct MatrixRegressionNative{T <: Real} <: ExampleInstanceNative{T}
    is_complex::Bool
    Y::Matrix
    X::Matrix
    lam_fro::Real # penalty on Frobenius norm
    lam_nuc::Real # penalty on nuclear norm
    lam_las::Real # penalty on l1 norm
    lam_glr::Real # penalty on penalty on row group l1 norm
    lam_glc::Real # penalty on penalty on column group l1 norm
end

function MatrixRegressionNative{T}(
    is_complex::Bool,
    n::Int,
    m::Int,
    p::Int,
    args...;
    A_max_rank::Int = div(m, 2) + 1,
    A_sparsity::Real = max(0.2, inv(sqrt(m * p))),
    Y_noise::Real = 0.01,
    ) where {T <: Real}
    @assert p >= m
    @assert 1 <= A_max_rank <= m
    @assert 0 < A_sparsity <= 1
    R = (is_complex ? Complex{T} : T)
    A_left = sprandn(R, p, A_max_rank, A_sparsity)
    A_right = sprandn(R, A_max_rank, m, A_sparsity)
    A = 10 * A_left * A_right
    X = randn(R, n, p)
    Y = Matrix{R}(X * A + Y_noise * randn(R, n, m))
    return MatrixRegressionNative{T}(is_complex, Y, X, args...)
end

function build(inst::MatrixRegressionNative{T}) where {T <: Real}
    (Y, X) = (inst.Y, inst.X)
    is_complex = inst.is_complex
    R = eltype(Y)
    @assert min(inst.lam_fro, inst.lam_nuc, inst.lam_las,
        inst.lam_glr, inst.lam_glc) >= 0
    (data_n, data_m) = size(Y)
    data_p = size(X, 2)
    @assert size(X, 1) == data_n
    @assert data_p >= data_m

    R_dim = (is_complex ? 2 : 1)
    data_pm = data_p * data_m
    data_nm = data_n * data_m
    # n is number of samples, m is number of responses, p is number of predictors
    # data_A (p * m) is matrix variable of coefficients, data_Y (n * m) is
    # response matrix, data_X (n * p) is design matrix
    model_n = 1 + R_dim * data_pm
    model_p = 0

    # loss function is 1/(2n) * ||Y - X * A||_fro^2
    # ||Y - X * A||^2 = tr((Y - X * A)' * (Y - X * A))
    # = tr((Y' - A' * X') * (Y - X * A))
    # = tr(A' * X' * X * A + Y' * Y - A' * X' * Y - Y' * X * A)
    # = tr(A' * X' * X * A) - 2 * tr(Y' * X * A) + tr(Y' * Y)
    # = sum(abs2, X * A) - 2 * dot(X' * Y, A) + sum(abs2, Y)
    # use EpiNormEucl for the sum(abs2, X * A) part
    # z >= ||Z||^2 is equivalent to (z, 1/2, vec(Z)) in EpiPerSquare is
    # equivalent to ((z + 1/2)/sqrt(2), (z - 1/2)/sqrt(2), vec(Z)) in EpiNormEucl
    if data_n > data_p
        # more samples than predictors: overdetermined
        # sum(abs2, X * A) = sum(abs2, cholesky(X' * X).U * A) =
        # sum(abs2, qr(X).R * A)
        # so can get a lower dimensional sum of squared terms
        Qhalf = qr(X).R
        # (Qhalf * A)_k,j = sum_k2 (Qhalf_k,k2 * A_k2,j)
        model_q = 2 + R_dim * data_pm
    else
        # fewer (or equal) samples than predictors
        # X * A is already low dimensional
        Qhalf = X
        # (X * A)_i,j = sum_k2 (X_i,k2 * A_k2,j)
        model_q = 2 + R_dim * data_nm
    end

    model_h = zeros(T, model_q)
    rtT2 = sqrt(T(2))
    model_h[1] = rtT2 / 4
    model_h[2] = -model_h[1]

    model_G = zeros(T, model_q, model_n)
    model_G[1:2, 1] .= -inv(rtT2)
    r_idx = 3
    for l in 1:min(data_p, data_n)
        c_idx = 2
        for j in 1:data_m
            if is_complex
                for k2 in 1:data_p
                    Qlk2 = Qhalf[l, k2]
                    model_G[r_idx, c_idx] = real(Qlk2)
                    model_G[r_idx, c_idx + 1] = -imag(Qlk2)
                    model_G[r_idx + 1, c_idx] = imag(Qlk2)
                    model_G[r_idx + 1, c_idx + 1] = real(Qlk2)
                    c_idx += 2
                end
                r_idx += 2
            else
                @. @views model_G[r_idx, c_idx:(c_idx + data_p - 1)] = Qhalf[l, :]
                r_idx += 1
                c_idx += data_p
            end
        end
    end

    cones = Cones.Cone{T}[Cones.EpiNormEucl{T}(model_q)]

    # put 1/(2n) * (epi(sum(abs2, X * A)) - 2 * dot(X' * Y, A)) in objective
    # ignore constant term 1/(2n) * sum(abs2, Y)
    model_c = zeros(T, model_n)
    mXpY = -X' * Y
    if is_complex
        @views Cones.vec_copyto!(model_c[2:end], vec(mXpY))
    else
        model_c[2:end] .= vec(mXpY)
    end
    model_c[1] = inv(T(2))
    model_c /= data_n

    # list of optional regularizers (group lasso handled separately below)
    reg_cone_dim = 1 + R_dim * data_pm
    lams = [
        (inst.lam_fro, Cones.EpiNormEucl{T}(reg_cone_dim)),
        (inst.lam_las, Cones.EpiNormInf{T, R}(reg_cone_dim, use_dual = true)),
        (inst.lam_nuc, Cones.EpiNormSpectral{T, R}(
            data_m, data_p, use_dual = true)),
        ]

    for (lam, cone) in lams
        if iszero(lam)
            continue # ignore regularizer because corresponding is zero
        end
        # obj term: lam * z
        # constraint: (z, data_A) in cone
        push!(model_c, lam)

        if cone isa Cones.EpiNormSpectral{T, R}
            # permute identity: need transpose of data_A since data_p >= data_m
            iden_mat = zeros(T, R_dim * data_pm, R_dim * data_pm)
            for j in 1:data_m, k in 1:data_p
                if is_complex
                    r_idx = (2 * j - 1) + 2 * (k - 1) * data_m
                    c_idx = (2 * k - 1) + 2 * (j - 1) * data_p
                    iden_mat[r_idx, c_idx] = 1
                    iden_mat[r_idx + 1, c_idx + 1] = 1
                else
                    r_idx = j + (k - 1) * data_m
                    c_idx = k + (j - 1) * data_p
                    iden_mat[r_idx, c_idx] = 1
                end
            end
        else
            iden_mat = Matrix(-one(T) * I, R_dim * data_pm, R_dim * data_pm)
        end
        model_G = T[
            model_G  zeros(T, model_q, 1);
            zeros(T, 1, model_n)  -one(T);
            zeros(T, R_dim * data_pm, 1)  iden_mat  zeros(T, R_dim * data_pm,
                model_n - R_dim * data_pm);
            ]
        append!(model_h, zeros(reg_cone_dim))

        model_n += 1
        model_q += reg_cone_dim

        push!(cones, cone)
    end

    # row group lasso regularizer (one lambda for all rows)
    if !iszero(inst.lam_glr)
        append!(model_c, fill(inst.lam_glr, data_p))

        q_glr = data_p * (1 + R_dim * data_m)
        G_glr = zeros(T, q_glr, model_n + data_p)
        r_idx = 1
        for k in 1:data_p
            G_glr[r_idx, model_n + k] = -1
            r_idx += 1
            for j in 1:data_m
                if is_complex
                    c_idx = 1 + (2 * k - 1) + 2 * (j - 1) * data_p
                    G_glr[r_idx, c_idx] = -1
                    r_idx += 1
                    G_glr[r_idx, c_idx + 1] = -1
                    r_idx += 1
                else
                    c_idx = 1 + k + (j - 1) * data_p
                    G_glr[r_idx, c_idx] = -1
                    r_idx += 1
                end
            end
        end

        append!(model_h, zeros(q_glr))
        model_G = T[
            model_G  zeros(model_q, data_p);
            G_glr;
            ]

        model_n += data_p
        model_q += q_glr

        append!(cones, Cones.Cone{T}[Cones.EpiNormEucl{T}(
            1 + R_dim * data_m) for k in 1:data_p])
    end

    # column group lasso regularizer (one lambda for all columns)
    if !iszero(inst.lam_glc)
        append!(model_c, fill(inst.lam_glc, data_m))

        q_glc = data_m * (1 + R_dim * data_p)
        G_glc = zeros(T, q_glc, model_n + data_m)
        r_idx = 1
        for j in 1:data_m
            G_glc[r_idx, model_n + j] = -1
            r_idx += 1
            for k in 1:data_p
                if is_complex
                    c_idx = 1 + (2 * k - 1) + 2 * (j - 1) * data_p
                    G_glc[r_idx, c_idx] = -1
                    r_idx += 1
                    G_glc[r_idx, c_idx + 1] = -1
                    r_idx += 1
                else
                    c_idx = 1 + k + (j - 1) * data_p
                    G_glc[r_idx, c_idx] = -1
                    r_idx += 1
                end
            end
        end

        append!(model_h, zeros(q_glc))
        model_G = T[
            model_G  zeros(model_q, data_m);
            G_glc;
            ]

        model_n += data_m
        model_q += q_glc

        append!(cones, Cones.Cone{T}[Cones.EpiNormEucl{T}(1 + R_dim * data_p)
            for k in 1:data_m])
    end

    model = Models.Model{T}(model_c, zeros(T, 0, model_n), zeros(T, 0),
        model_G, model_h, cones)
    return model
end

function test_extra(
    inst::MatrixRegressionNative{T},
    solve_stats::NamedTuple,
    solution::NamedTuple,
    ) where T
    @test solve_stats.status == Solvers.Optimal
    (solve_stats.status == Solvers.Optimal) || return

    # check objective value is correct
    (Y, X) = (inst.Y, inst.X)
    A_opt = zeros(eltype(Y), size(X, 2), size(Y, 2))
    A_len = length(A_opt) * (inst.is_complex ? 2 : 1)
    @views Cones.vec_copyto!(A_opt, solution.x[1 .+ (1:A_len)])
    loss = (sum(abs2, X * A_opt) / 2 - real(dot(X' * Y, A_opt))) / size(Y, 1)
    obj_result = loss +
        inst.lam_fro * norm(A_opt, 2) +
        inst.lam_nuc * sum(svd(A_opt).S) +
        inst.lam_las * norm(A_opt, 1) +
        inst.lam_glr * sum(norm, eachrow(A_opt)) +
        inst.lam_glc * sum(norm, eachcol(A_opt))
    tol = eps(T)^0.25
    @test solve_stats.primal_obj ≈ obj_result atol=tol rtol=tol
    return
end


text:
using HDF5
using Base.Test

@testset "external" begin

    # roughly following https://www.hdfgroup.org/ftp/HDF5/current/src/unpacked/examples/h5_extlink.c
    source_file = h5open(tempname(), "w")
    agroup = g_create(source_file, "agroup")
    target_file = h5open(tempname(), "w")
    target_group = g_create(target_file, "target_group")
    target_group["abc"]="abc"
    target_group["1"]=1
    target_group["1.1"]=1.1
    close(target_file)

    # create external link such that source_file["ext_link"] points to target_file["target_group"]
    # test both an HDF5File and an HDF5Group for first argument
    HDF5.create_external(source_file, "ext_link", target_file.filename, "target_group")
    HDF5.create_external(agroup, "ext_link", target_file.filename, "target_group")
    # write some things via the external link
    new_group = g_create(source_file["ext_link"], "new_group")
    new_group["abc"]="abc"
    new_group["1"]=1
    new_group["1.1"]=1.1

    # read things from target_group via exernal link created with HDF5File argument
    group = source_file["ext_link"]
    @test read(group["abc"])=="abc"
    @test read(group["1"])==1
    @test read(group["1.1"])==1.1
    # read things from target_group via the external link created with HDF5Group argument
    groupalt = source_file["agroup/ext_link"]
    @test read(groupalt["abc"])=="abc"
    @test read(groupalt["1"])==1
    @test read(groupalt["1.1"])==1.1
    close(source_file)

    ##### tests that should be included but don't work
    # when ggggggggg restarts julia and keeps track of target_file.filename,
    # these tests succeed
    # reopening the target_file crashes due to "file close degree doesn't match"
    # target_file = h5open(target_file.filename, "r")
    # group2 = target_file["target_group"]["new_group"]
    # @test read(group2["abc"])=="abc"
    # @test read(group2["1"])==1
    # @test read(group2["1.1"])==1.1

end


text:
## test expr creation
using ASTModels

@expr a b k
c = 2a
d = b*c

## test dictionary models

di = Dict("inputs"=>[a,b,k], "outputs"=>[c, k])
merge_models(di)

di1 = Dict("inputs"=>[a], "outputs"=>[c])
di2 = Dict("inputs"=>[c, b], "outputs"=>[d])
dm = merge_models(di1, di2)

# clone still does not make sense, as parameters are not cloned - We would have to generate new parameters....
# TODO only true inputs would have to be preserved ..
dc = clone(dm)
dc["inputs"] === dm["inputs"]
dc["outputs"] === dm["outputs"]
dc["outputs"] == dm["outputs"]

ti1 = TestTypeWithModel(di1)
ti2 = TestTypeWithModel(di2)
merge_models(ti1, ti2)

## test function compilation

# macroexpand(:(@ast_func h(a,b) = d))
@astfunc h(a,b) = d

anonym = astfunc([a,b], d)
anonym2 = @astfunc (a,b) -> d

d2 = clone([a,b], d)

## TODO test astfunc and clone !!!