text:
mutable struct ParamType{Tdim, Tsol, Tres, Tmsh} <: AbstractParamType{Tdim}
  f::IOStream
  t::Float64  # current time value
  order::Int  # accuracy of elements (p=1,2,3...)

  q_vals::Array{Tsol, 1}  # resuable temporary storage for q variables at a node
  qg::Array{Tsol, 1}  # reusable temporary storage for boundary condition
  v_vals::Array{Tsol, 1}  # reusable storage for convert back to entropy vars.

  res_vals1::Array{Tres, 1}  # reusable residual type storage
  res_vals2::Array{Tres, 1}  # reusable residual type storage

  flux_vals1::Array{Tres, 1}  # reusable storage for flux values
  flux_vals2::Array{Tres, 1}  # reusable storage for flux values

  sat_vals::Array{Tres, 1}  # reusable storage for SAT term

  A0::Array{Tsol, 2}  # reusable storage for the A0 matrix
  A0inv::Array{Tsol, 2}  # reusable storage for inv(A0)
  A1::Array{Tsol, 2}  # reusable storage for a flux jacobian
  A2::Array{Tsol, 2}  # reusable storage for a flux jacobian

  A_mats::Array{Tsol, 3}  # reusable storage for flux jacobians

  Rmat1::Array{Tres, 2}  # reusable storage for a matrix of type Tres
  Rmat2::Array{Tres, 2}

  nrm::Array{Tmsh, 1}  # a normal vectora

  edgestab_gamma::Float64  # edge stabilization parameter
  # debugging options
  writeflux::Bool  # write Euler flux
  writeboundary::Bool  # write boundary data
  writeq::Bool # write solution variables
  use_edgestab::Bool  # use edge stabilization
  use_filter::Bool  # use filtering
  use_res_filter::Bool # use residual filtering

  filter_mat::Array{Float64, 2}  # matrix that performs filtering operation
  # includes transformations to/from modal representation

  use_dissipation::Bool  # use artificial dissipation
  dissipation_const::Float64  # constant used for dissipation filter matrix

  tau_type::Int  # type of tau to use for GLS stabilization

  krylov_itr::Int  # Krylov iteration number for iterative solve
  krylov_type::Int # 1 = explicit jacobian, 2 = jac-vec prod
  time::Timings

  function ParamType{Tdim, Tsol, Tres, Tmsh}(mesh, sbp, opts, order::Integer) where {Tdim, Tsol, Tres, Tmsh}
    # create values, apply defaults

    t = 0.0
    myrank = mesh.myrank
    f = open("log_$myrank.dat", "w")
    q_vals = Array{Tsol}(4)
    qg = Array{Tsol}(4)
    v_vals = Array{Tsol}(4)

    res_vals1 = Array{Tres}(4)
    res_vals2 = Array{Tres}(4)

    flux_vals1 = Array{Tres}(4)
    flux_vals2 = Array{Tres}(4)

    sat_vals = Array{Tres}(4)

    A0 = zeros(Tsol, 4, 4)
    A0inv = zeros(Tsol, 4, 4)
    A1 = zeros(Tsol, 4, 4)
    A2 = zeros(Tsol, 4, 4)
    A_mats = zeros(Tsol, 4, 4, Tdim)

    Rmat1 = zeros(Tres, 4, 4)
    Rmat2 = zeros(Tres, 4, 4)

    nrm = zeros(Tmsh, 2)

    edgestab_gamma = opts["edgestab_gamma"]

    # debugging options
    writeflux = opts[ "writeflux"]
    writeboundary = opts[ "writeboundary"]
    writeq = opts["writeq"]
    use_edgestab = opts["use_edgestab"]
    if use_edgestab println("edge stabilization enabled") end

    use_filter = opts["use_filter"]
    if use_filter println("solution variables filter enabled") end

    use_res_filter = opts["use_res_filter"]
    if use_res_filter println("residual filter enabled") end

    if use_filter || use_res_filter || opts["use_filter_prec"]
      filter_fname = opts["filter_name"]
      filter_mat = calcFilter(sbp, filter_fname, opts)
    else
      filter_mat = Array{Float64}(0,0)
    end

    use_dissipation = opts["use_dissipation"]
    if use_dissipation println("artificial dissipation enabled") end

    dissipation_const = opts["dissipation_const"]

    tau_type = opts["tau_type"]

    krylov_itr = 0
    krylov_type = 1 # 1 = explicit jacobian, 2 = jac-vec prod

    time = Timings()
    return new(f, t, order, q_vals, qg, v_vals, res_vals1, res_vals2, sat_vals, flux_vals1,
               flux_vals2, A0, A0inv, A1, A2, A_mats, Rmat1, Rmat2, nrm,
               edgestab_gamma, writeflux, writeboundary,
               writeq, use_edgestab, use_filter, use_res_filter, filter_mat,
               use_dissipation, dissipation_const, tau_type,
               krylov_itr, krylov_type,
               time)

  end   # end of ParamType function

end  # end type declaration

mutable struct EllipticData_{Tsol, Tres, Tdim, Tmsh} <: EllipticData{Tsol, Tres, Tdim}
  params::ParamType{Tdim, Tsol, Tres, Tmsh}
  comm::MPI.Comm
  commsize::Int
  myrank::Int
  #
  # Debug Use Only
  #
  xy_face::Array{Tsol, 4} # (x/y, L/R, numnodes, iface)
  #
  # End Debug
  #
  area_sum::Array{Tmsh, 1}
  w::Array{Tmsh, 2} # H/|J|, (numNodesPerElement, numEl)
  q::Array{Tsol, 3}
  q_irk::Array{Tsol, 4}
  q1::Array{Float64, 3}
  q2::Array{Float64, 3}
  q_face::Array{Tsol, 4} # (numDofs, left/right, numNodesPerFace, numFaces)
  # flux_parametric::Array{Tsol, 2}

  xflux_face::Array{Tres, 3}    # stores (u+ - u-)nx*, (numDofs, numNodes, numFaces)
  yflux_face::Array{Tres, 3}    # stores (u+ - u-)ny*
  flux_face ::Array{Tres, 3}    # stores σ⋅norm

  xflux_bndry::Array{Tres, 3}    # stores (u+ - u-)nx*, (numDofs, numNodes, numFaces)
  yflux_bndry::Array{Tres, 3}    # stores (u+ - u-)ny*
  flux_bndry ::Array{Tres, 3}    # stores σ⋅norm

  res::Array{Tres, 3}            # result of computation (numDofs, numNodesPerElem, numElems)
  res_irk::Array{Tres, 4}            # result of computation (numDofs, numNodesPerElem, numElems)
  res1::Array{Float64, 3}            # result of computation (numDofs, numNodesPerElem, numElems)
  res_vec::Array{Tres, 1}         # result of computation in vector form
  q_vec::Array{Tres,1}            # initial condition in vector form
  q_bndry::Array{Tsol, 3}         # store solution variables interpolated to

  shared_data::Array{SharedFaceData{Tsol}, 1}  # MPI send and receive buffers
  #q_face_send::Array{Array{Tsol, 3}, 1}    # send buffers for sending q values
  # to other processes
  #q_face_recv::Array{Array{Tsol, 3}, 1}    # recieve buffers for q values
  #flux_sharedface::Array{Array{Tres, 3}, 1}    # hold shared face flux
  # derivative
  src_func::SRCType    # functor for source term
  src::Array{Tsol, 3} # (numDofPerNode, numNodesPerElement, numEl)
  #
  # Declaration of variables exclusively used in elliptic problems
  #
  # flux_func::FluxType    # functor for the face flux
  diffusion_func::AbstractDiffn
  functional::AbstractFunctional
  calc_energy::AbstractFunctional
  energy::Array{Tsol, 1}
  # diffusion coefficients
  lambda::Array{Tsol, 5} # (d1, d2, ndof, numnodes, n_elemes)
  lambda_face::Array{Tsol, 6} # (d2, d1, ndof, L/R, numNodesPerFace, n_faces)
  lambda_bndry::Array{Tsol, 5} # (d1, d2, ndof, numNodesPerface, n_bndfaces)
  q_grad::Array{Tsol, 4}    # gradient of q, (dim, numDofPerNode, numNodesPerElem, numElems)
  #
  # TBD: Not sure if we store gradient of q at face cubature points
  #
  q_grad_face::Array{Tsol, 5} #(dim, numDofPerNode, left/right, numNodesPerFace, numFaces)
  q_grad_bndry::Array{Tsol, 4} #(dim, numDofPerNode, numNodesPerFace, numBndries)
  #
  # end declaration of variables exclusively used in elliptic problems
  #
  Minv::Array{Float64, 1}
  M::Array{Float64, 1}
  res_edge::Array{Tres, 4}       # edge based residual used for stabilization
  # (numdof, numNodesPerElement, numEl, numEdgesPerEl)
  flux_func::FluxType
  majorIterationCallback::Function
  # assembleSolution::Function
  disassembleSolution::Function
  multiplyA0inv::Function

  q_face_send::Array{Array{Tsol, 3}, 1}  # send buffers for sending q values to other processes
  q_face_recv::Array{Array{Tsol, 3}, 1}  # recieve buffers for q values

  nstages::UInt8
  istage::UInt8

  function EllipticData_{Tsol, Tres, Tdim, Tmsh}(mesh::AbstractMesh, sbp::AbstractOperator, opts) where {Tsol, Tres, Tdim, Tmsh}
    println("\nConstructing EllipticData object")
    println("  Tsol = ", Tsol)
    println("  Tres = ", Tres)
    println("  Tdim = ", Tdim)
    println("  Tmsh = ", Tmsh)

    eqn = new()

    eqn.comm = mesh.comm
    eqn.commsize = mesh.commsize
    eqn.myrank = mesh.myrank

    eqn.params = ParamType{Tdim, Tsol, Tres, Tmsh}(mesh, sbp, opts, mesh.order)

    numfacenodes = mesh.numNodesPerFace
    numfaces = mesh.numInterfaces
    numBndFaces = mesh.numBoundaryFaces
    numvars  = mesh.numDofPerNode
    #
    # volume variables
    #
    # TODO: switch the last two dimensions in all gradient variable
    # in order to keep consistency with lambda
    #
    eqn.w      = zeros(Tmsh, mesh.numNodesPerElement, mesh.numEl)
    eqn.q      = zeros(Tsol, numvars, sbp.numnodes, mesh.numEl)
    if haskey(opts, "TimeAdvance") && opts["TimeAdvance"] == "BDF2"
      eqn.q1      = zeros(Float64, numvars, sbp.numnodes, mesh.numEl)
      eqn.q2      = zeros(Float64, numvars, sbp.numnodes, mesh.numEl)
    end
    if haskey(opts, "TimeAdvance") && opts["TimeAdvance"] == "CN"
      eqn.q1      = zeros(Float64, numvars, sbp.numnodes, mesh.numEl)
    end
    eqn.res    = zeros(Tsol, numvars, sbp.numnodes, mesh.numEl)
    if haskey(opts, "TimeAdvance") && opts["TimeAdvance"] == "SDIRK4"
      # eqn.res_irk = zeros(Tsol, numvars, sbp.numnodes, mesh.numEl, 6)
      eqn.q_irk   = zeros(Tsol, numvars, sbp.numnodes, mesh.numEl, 6)
    end
    eqn.res1    = zeros(Float64, numvars, sbp.numnodes, mesh.numEl)
    eqn.src    = zeros(Tsol, numvars, sbp.numnodes, mesh.numEl)
    eqn.q_grad = zeros(Tsol, numvars, sbp.numnodes, mesh.numEl, Tdim)
    eqn.lambda = zeros(Tsol, Tdim, Tdim, numvars, sbp.numnodes, mesh.numEl)
    #
    # interface variables
    #
    eqn.q_face    = zeros(Tsol, numvars, 2, numfacenodes, numfaces)
    eqn.q_grad_face = zeros(Tsol, numvars, 2, numfacenodes, numfaces, Tdim)
    eqn.flux_face = zeros(Tsol, numvars, numfacenodes, numfaces)
    eqn.xflux_face = zeros(Tsol, numvars, numfacenodes, numfaces)
    eqn.yflux_face = zeros(Tsol, numvars, numfacenodes, numfaces)
    eqn.lambda_face = zeros(Tsol, Tdim, Tdim, numvars, 2, numfacenodes, numfaces)
    #
    # boundary variables
    #
    eqn.q_bndry    = zeros(Tsol, numvars, numfacenodes, numBndFaces)
    eqn.q_grad_bndry = zeros(Tsol, numvars, numfacenodes, numBndFaces, Tdim)
    eqn.flux_bndry = zeros(Tsol, numvars, numfacenodes, numBndFaces)
    eqn.xflux_bndry = zeros(Tsol, numvars, numfacenodes, numBndFaces)
    eqn.yflux_bndry = zeros(Tsol, numvars, numfacenodes, numBndFaces)
    eqn.lambda_bndry = zeros(Tsol, Tdim, Tdim, numvars, numfacenodes, numBndFaces)
    if mesh.isDG
      for i=1:mesh.npeers
        if opts["precompute_face_flux"]
          eqn.flux_sharedface[i] = zeros(Tres, mesh.numDofPerNode, numfacenodes,
                                         mesh.peer_face_counts[i])
        end
        # eqn.aux_vars_sharedface[i] = zeros(Tres, mesh.numDofPerNode,
                                        # numfacenodes, mesh.peer_face_counts[i])
      end
      eqn.shared_data = getSharedFaceData(Tsol, mesh, sbp, opts,
                                          opts["parallel_data"])
    else
      eqn.shared_data = Array{SharedFaceData}(0)
    end
#
    # reshape of q and res
    #
    eqn.q_vec = reshape(eqn.q, mesh.numDof)
    eqn.res_vec = reshape(eqn.res, mesh.numDof)

    # flux_func = FluxDict["SIPG"]
    eqn.majorIterationCallback = majorIterationCallback
    # eqn.assembleSolution = assembleSolution
    # eqn.disassembleSolution = disassembleSolution
    eqn.multiplyA0inv = multiplyA0inv

    eqn.Minv = calcMassMatrixInverse(mesh, sbp, eqn)
    eqn.M = calcMassMatrix(mesh, sbp, eqn)

    if opts["use_edge_res"]
      eqn.res_edge = zeros(Tres, mesh.numDofPerNode, sbp.numnodes, mesh.numEl,
                           mesh.numTypePerElement[2])
    else
      eqn.res_edge = zeros(Tres, 0, 0, 0, 0)
    end

    eqn.xy_face = zeros(Tsol, 2, 2, numfacenodes, numfaces)

    eqn.area_sum = zeros(Tmsh, mesh.numEl)
    calcWetArea(mesh, sbp, eqn)

    for el = 1 : mesh.numEl
      for n = 1 : mesh.numNodesPerElement
        eqn.w[n, el] = sbp.w[n]/mesh.jac[n, el]
      end
    end
    eqn.energy = zeros(Tsol, mesh.numDofPerNode)
    eqn.calc_energy = FunctionalDict["energy"]

    eqn.istage = 1
    eqn.nstages = 1
    if haskey(opts, "TimeAdvance") && opts["TimeAdvance"] == "SDIRK4"
      eqn.nstages = 5
    end

    # TODO: parallel related variables
    #
    eqn.q_face_send = Array{Array{Tsol, 3}}(1)
    eqn.q_face_recv = Array{Array{Tsol, 3}}(1)
    return eqn
  end # end function
end # end type

function calcWetArea(mesh::AbstractMesh{Tmsh},
                     sbp::AbstractOperator,
                     eqn::EllipticData{Tsol, Tres, Tdim}) where {Tmsh, Tsol, Tres, Tdim}
  nfaces = length(mesh.interfaces)
  nrm = zeros(Tmsh, mesh.numNodesPerFace, Tdim)
  area = zeros(Tmsh, mesh.numNodesPerFace)
  face_area = zero(Tmsh)
  sbpface = mesh.sbpface

  #
  # Compute the wet area of each element
  #
  for f = 1:nfaces
    face = mesh.interfaces[f]
    eL = face.elementL
    eR = face.elementR
    fL = face.faceL
    fR = face.faceR
    #
    # Compute the size of face
    face_area = 0.0
    for n=1:mesh.numNodesPerFace

      dxidx = sview(mesh.dxidx_face, :, :, n, f)
      # norm vector in reference element
      # nrm_xi = sview(sbp.facenormal, :, fL)
      # nrm_xi = sview(mesh.sbpface.normal, :, fL)
      # nrm[n,1] = dxidx[1, 1]*nrm_xi[1] + dxidx[2, 1]*nrm_xi[2]
      # nrm[n,2] = dxidx[1, 2]*nrm_xi[1] + dxidx[2, 2]*nrm_xi[2]
      # area[n] = sqrt(nrm_xy[n,1]*nrm_xy[n,1] + nrm_xy[n,2]*nrm_xy[n,2])

      nrm_xy = ro_sview(mesh.nrm_face, :, n, f)
      area[n] = norm(nrm_xy)
      face_area += sbpface.wface[n]*area[n]
    end
    eqn.area_sum[eL] += face_area
    eqn.area_sum[eR] += face_area
  end

  for bc = 1:mesh.numBC
    indx0 = mesh.bndry_offsets[bc]
    indx1 = mesh.bndry_offsets[bc+1] - 1

    for f = indx0:indx1
      face = mesh.bndryfaces[f].face
      elem = mesh.bndryfaces[f].element
      #
      # Compute the size of face
      face_area = 0.0
      for n=1:mesh.numNodesPerFace
        # dxidx = sview(mesh.dxidx_bndry, :, :, n, f)
        # # norm vector in reference element
        # nrm_xi = sview(mesh.sbpface.normal, :, face)
        # nrm[n,1] = dxidx[1, 1]*nrm_xi[1] + dxidx[2, 1]*nrm_xi[2]
        # nrm[n,2] = dxidx[1, 2]*nrm_xi[1] + dxidx[2, 2]*nrm_xi[2]
        # area[n] = sqrt(nrm[n,1]*nrm[n,1] + nrm[n,2]*nrm[n,2])
        nrm_xy = ro_sview(mesh.nrm_bndry, :, :, f)
        area[n] = norm(nrm_xy)
        face_area += sbpface.wface[n]*area[n]
      end
      eqn.area_sum[elem] += 2.0*face_area
    end
  end
  return nothing
end

function majorIterationCallback(itr::Integer,
                                mesh::AbstractMesh,
                                sbp::AbstractOperator,
                                eqn::AbstractEllipticData,
                                opts,
                                f::IO)
  if opts["write_vis"] && ((itr % opts["output_freq"]) == 0 || itr == 1)
    vals = real(eqn.q_vec)
    saveSolutionToMesh(mesh, vals)
    fname = string("solution_", itr)
    writeVisFiles(mesh, fname)
  end

  return nothing
end

function multiplyA0inv(mesh::AbstractMesh{Tmsh},
                       sbp::AbstractOperator,
                       eqn::EllipticData{Tsol, Tres, Tdim},
                       opts,
                       res::AbstractArray{Tsol, 3}) where {Tmsh, Tsol, Tdim, Tres}
  return nothing
end


text:
# These are various tasks to help with package development.

using Maker
using Glob

Maker.task("cleancov") do
    for fn in glob("*/*.cov")
        rm(fn)
    end
end

Maker.task("coverage", "cleancov") do
    cd("test") do
        run(`julia --code-coverage=user --inline=no runtests.jl`)
    end
    # coverage = Coverage.process_folder("src")
    # @show Coverage.get_summary(coverage)
end

Maker.task("tests") do
    Pkg.test("Maker")
end

Maker.task("localtests") do
    cd("$(Pkg.dir())/Maker/test") do
        include("runtests.jl")
    end
end

Maker.task("register") do
    Pkg.register("Maker")
end

Maker.task("tagminor") do
    Pkg.tag("Maker", :minor)
end

Maker.task("tagpatch") do
    Pkg.tag("Maker", :patch)
end

Maker.task("publish") do
    Pkg.publish()
end

Maker.task("pull") do
    run(`git pull origin $(Pkg.Git.branch())`)
end

Maker.task("push") do
    run(`git push origin $(Pkg.Git.branch())`)
end

make(ARGS)

ext: jl


text:
# Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)

function eponesweep!(epfields::EPFields{T}, epalg::EPAlg, epmat::EPMat, stat) where T
    @extract epfields : av va a d μ s siteflagave siteflagvar
    @extract epalg : alpha beta_vec minvar maxvar epsconv damp
    @extract epmat : αStS αStS_diag Σ lb ub αStb v

    # Parameters of the approximated join
    # We rebuild αStS and add D
    # αStS = αStS if diag(αStS) = αStS_diag
    # Σ^-1 = (αStS + D)
    Σ⁻¹ = αStS[diagind(αStS)] = αStS_diag + inv.(d)

    # Σ = inv(Σ^-1)
    stat[:elapsed_eponesweep_inv] = @elapsed inplaceinverse!(Σ, Σ⁻¹)

    # αStb = βF'y
    # v¯ = Σ(αStb + Da) (original ep)
    # mul!(v, Σ, (αStb + a./d))
    # v¯ = Σ(αStb + Da) + Σ * beta_vec (maxent-ep)
    v .= Σ * (αStb + a./d) + Σ * beta_vec

    errav = errva = errμ = errs = typemin(T)
    for i in eachindex(av)

        # Parameters of the normal part of the nth tilted
        newμ, news = newμs(Σ[i,i],a[i],d[i],v[i],lb[i],ub[i],minvar, maxvar)
        errμ = max(errμ, abs(μ[i]-newμ))
        errs = max(errs, abs(s[i]-news))
        μ[i] = newμ
        s[i] = news

        # Parameters of the whole tilted (truncated marginals)
        newave, newva = newav(s[i],μ[i],av[i],va[i],siteflagave[i],siteflagvar[i],lb[i],ub[i],minvar,maxvar)
        errav = max(errav,abs(av[i]-newave))
        errva = max(errva,abs(va[i]-newva))
        av[i] = newave
        va[i] = newva

        # Parameter of the prior
        newa,newb = matchmom(μ[i],s[i],av[i],va[i],minvar,maxvar)
        a[i] = damp * a[i] + (1.0-damp)*newa
        d[i] = damp * d[i] + (1.0-damp)*newb
    end

    return errav,errva,errμ,errs
end


text:
using JSON

# script to parse Brenda files -
function parse_brenda_file(path_to_brenda_file)

    # iniialize -
    brenda_records_array = Dict{String,Any}[]

    # process the file -
    open(path_to_brenda_file) do f

        # do stuff with the open file
        buffer = read(f, String)
        list_of_records = split(buffer,'!')

        # how many records do we have?
        number_of_records = length(list_of_records)
        for record_index = 1:number_of_records

            # initalize empty dictionary -
            record_dictionary = Dict{String,Any}()

            # convert the record to a string -
            record = string(list_of_records[record_index])

            # split into fields -
            field_array = split(record,'#')
            number_of_fields = length(field_array)
            for field_index = 1:number_of_fields

                # split into key*value pairs -
                key_value_array = split(field_array[field_index],'*')

                # split again if we have two elements -
                if (length(key_value_array) == 2)

                    # split into key,value pair
                    key = string(key_value_array[1])
                    value = string(key_value_array[2])

                    # cache the values -
                    record_dictionary[key] = value

                end
            end

            # add this record to the array -
            push!(brenda_records_array,record_dictionary)
        end
    end

    # return the array -
    return brenda_records_array
end

function main(path_to_model_files,path_to_json_file)

    # initialize -
    brenda_dictionary = Dict{String,Any}()
    brenda_records_array = Dict{String,Any}[]

    # brenda files end w/.dat
    for file in filter(x -> endswith(x, "dat"), readdir(path_to_model_files))

        # path to specific file -
        path_to_brenda_file = "$(path_to_model_files)/$(file)"

        # parse the brenda file,
        brenda_records_array = parse_brenda_file(path_to_brenda_file)

        # grab the first record, and get the ecnumber -
        ec_number = brenda_records_array[1]["ecNumber"]
        ec_number_key = "ec:$(ec_number)"

        # cache this records collection -
        brenda_dictionary[ec_number_key] = brenda_records_array
    end

    # write dictionary to JSON -
    file = open(path_to_json_file,"w")
    JSON.print(file,brenda_dictionary,6)
    close(file)
end

# build the dictionary -
main("../config/data/brenda_data","../config/data/ECN-EcoliCoreModel-Palsson-EcoSal-Plus-2013.json")


text:
# SINGLE NODE/CORE VERSION

import Pkg; Pkg.instantiate()

include("/users/yh31/scratch/projects/gigaword_64k/src/gigaword_64k.jl")

const path_nyt_dir = "/users/yh31/scratch/datasets/entity_linking/raw_data/gigaword/giga5/data/nyt_eng"

const path_intermed_data = "/users/yh31/scratch/datasets/entity_linking/raw_data/gigaword/giga_topk/intermediate_files"
const path_intermed_nyt = joinpath(path_intermed_data, "nyt")

process_part_of_tree(path_nyt_dir, path_intermed_nyt, 4)


text:
# This file is a part of SimilaritySearch.jl

"""
    append!(
        index::SearchGraph,
        db;
        neighborhood=Neighborhood(),
        parallel_block=1,
        parallel_minimum_first_block=parallel_block,
        callbacks=SearchGraphCallbacks(),
        pools=getpools(index)
    )

Appends all items in db to the index. It can be made in parallel or sequentially.

Arguments:

- `index`: the search graph index
- `neighborhood`: A [`Neighborhood`](@ref) object that specifies the kind of neighborhood that will be computed.
- `db`: the collection of objects to insert, an `AbstractDatabase` is the canonical input, but supports any iterable objects
- `parallel_block`: The number of elements that the multithreading algorithm process at once,
    it is important to be larger that the number of available threads but not so large since the quality of the search graph could degrade (a few times the number of threads is enough).
    If `parallel_block=1` the algorithm becomes sequential.
- `parallel_minimum_first_block`: The number of sequential appends before running parallel.
Note: Parallel doesn't trigger callbacks inside blocks.
- `callbacks`: A `SearchGraphCallbacks` object to be called after some insertions
    (specified by the `callbacks` object). These callbacks are used to maintain the algorithm
    in good shape after many insertions (adjust hyperparameters and the structure).
- `pools`: The set of caches used for searching.

Note 1: Callbacks are not executed inside parallel blocks
Note 2: Callbacks will be ignored if `callbacks=nothing`

"""
function Base.append!(
        index::SearchGraph,
        db;
        neighborhood=Neighborhood(),
        parallel_block=1,
        parallel_minimum_first_block=parallel_block,
        callbacks=SearchGraphCallbacks(),
        pools=getpools(index)
    )
    db = convert(AbstractDatabase, db)
    append!(index.db, db)

    parallel_block == 1 && return _sequential_append_loop!(index, neighborhood, callbacks, pools)

    n = length(index) + length(db)
    m = 0

    parallel_minimum_first_block = min(parallel_minimum_first_block, n)
    @inbounds while length(index) < parallel_minimum_first_block
        m += 1
        push_item!(index, db[m], neighborhood, false, callbacks, pools)
    end

    sp = length(index) + 1
    sp > n && return index

    resize!(index.links, n)
    _parallel_append_loop!(index, neighborhood, pools, sp, n, parallel_block, callbacks)
end

"""
    index!(index::SearchGraph; parallel_block=1, parallel_minimum_first_block=parallel_block, callbacks=SearchGraphCallbacks())

Indexes the already initialized database (e.g., given in the constructor method). It can be made in parallel or sequentially.
The arguments are the same than `append!` function but using the internal `index.db` as input.

"""
function index!(
        index::SearchGraph;
        neighborhood=Neighborhood(),
        parallel_block=1,
        parallel_minimum_first_block=parallel_block,
        callbacks=SearchGraphCallbacks(),
        pools=getpools(index)
    )
    @assert length(index) == 0 && length(index.db) > 0
    parallel_block == 1 && return _sequential_append_loop!(index, neighborhood, callbacks, pools)

    m = 0
    db = index.db
    n = length(db)

    parallel_minimum_first_block = min(parallel_minimum_first_block, n)
    @inbounds while length(index) < parallel_minimum_first_block
        m += 1
        push_item!(index, db[m], neighborhood, false, callbacks, pools)
    end

    sp = length(index) + 1
    sp > n && return index
    resize!(index.links, n)
    _parallel_append_loop!(index, neighborhood, pools, sp, n, parallel_block, callbacks)
end

function _sequential_append_loop!(index::SearchGraph, neighborhood::Neighborhood, callbacks, pools::SearchGraphPools)
    i = length(index)
    n = length(index.db)
    @inbounds while i < n
        i += 1
        push_item!(index, index.db[i], neighborhood, false, callbacks, pools)
    end

    index
end

function _connect_links(index, sp, ep)
    Threads.@threads for i in sp:ep
        @inbounds for id in index.links[i]
            lock(index.locks[id])
            try
                push!(index.links[id], i)
                # sat_should_push(index.links[id], index, index[i], i, -1.0) && push!(index.links[id], i)
            finally
                unlock(index.locks[id])
            end
        end
    end
end

function _parallel_append_loop!(index::SearchGraph, neighborhood::Neighborhood, pools::SearchGraphPools, sp, n, parallel_block, callbacks)
    while sp < n
        ep = min(n, sp + parallel_block)
        index.verbose && rand() < 0.01 && println(stderr, "appending chunk ", (sp=sp, ep=ep, n=n), " ", Dates.now())

        # searching neighbors
        # @show length(index.links), length(index.db), length(db), length(index.locks), length(index), sp, ep
        Threads.@threads for i in sp:ep
            @inbounds index.links[i] = find_neighborhood(index, index.db[i], neighborhood, pools)
        end

        # connecting neighbors
        _connect_links(index, sp, ep)

        # increasing locks => new items are enabled for searching (and reported by length so they can also be hints)
        resize!(index.locks, ep)
        for i in sp:ep
            @inbounds index.locks[i] = Threads.SpinLock()
        end

        # apply callbacks
        callbacks !== nothing && execute_callbacks(callbacks, index, sp, ep)
        sp = ep + 1
    end

    index
end

"""
    push!(
        index::SearchGraph,
        item;
        neighborhood=Neighborhood(),
        push_item=true,
        callbacks=SearchGraphCallbacks(),
        pools=getpools(index)
    )

Appends an object into the index. It accepts the same arguments that `push!` but assuming some default values.

Arguments:

- `index`: The search graph index where the insertion is going to happen
- `item`: The object to be inserted, it should be in the same space than other objects in the index and understood by the distance metric.
- `neighborhood`: A [`Neighborhood`](@ref) object that specifies the kind of neighborhood that will be computed.
- `push_item`: if `push_item=false` is an internal option, used by `append!` and `index!` (it avoids to insert `item` into the database since it is already inserted but not indexed)
- `callbacks`: The set of callbacks that are called whenever the index grows enough. Keeps hyperparameters and structure in shape.
- `pools`: The set of caches used for searching.

- Note: `callbacks=nothing` ignores the execution of any callback
"""
@inline function push!(
        index::SearchGraph,
        item;
        neighborhood=Neighborhood(),
        push_item=true,
        callbacks=SearchGraphCallbacks(),
        pools=getpools(index)
    )

    push_item!(index, item, neighborhood, push_item, callbacks, pools)
end

"""
    push_item!(
        index::SearchGraph,
        item,
        neighborhood,
        push_item,
        callbacks,
        pools
    )

Appends an object into the index

Arguments:

- `index`: The search graph index where the insertion is going to happen.
- `item`: The object to be inserted, it should be in the same space than other objects in the index and understood by the distance metric.
- `neighborhood`: A [`Neighborhood`](@ref) object that specifies the kind of neighborhood that will be computed.
- `push_item`: if `false` is an internal option, used by `append!` and `index!` (it avoids to insert `item` into the database since it is already inserted but not indexed).
- `callbacks`: The set of callbacks that are called whenever the index grows enough. Keeps hyperparameters and structure in shape.
- `pools`: The set of caches used for searching.

- Note: setting `callbacks` as `nothing` ignores the execution of any callback
"""
@inline function push_item!(
    index::SearchGraph,
    item,
    neighborhood::Neighborhood,
    push_item::Bool,
    callbacks::SearchGraphCallbacks,
    pools::SearchGraphPools
)
    neighbors = find_neighborhood(index, item, neighborhood, pools)
    push_neighborhood!(index, item, neighbors, callbacks; push_item)

    neighbors
end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/SimilaritySearch

max_forks_repo_name: sadit/SimilaritySearch


text:
using Test
using TensorKit
using TensorKitManifolds
using LinearAlgebra
using MERAKit
using Logging
using Random

function particle_number_operator(::Type{Vect[FibonacciAnyon]})
    V = Vect[FibonacciAnyon](FibonacciAnyon(:I) => 1, FibonacciAnyon(:τ) => 1)
    z = TensorMap(zeros, Float64, V ← V)
    z.data[FibonacciAnyon(:I)] .= 0.0
    z.data[FibonacciAnyon(:τ)] .= 1.0
    return z
end

function particle_number_operator(::Type{Vect[IsingAnyon]})
    V = Vect[IsingAnyon](
        IsingAnyon(:I) => 1,
        IsingAnyon(:σ) => 1,
        IsingAnyon(:ψ) => 1
    )
    z = TensorMap(zeros, Float64, V ← V)
    z.data[IsingAnyon(:I)] .= 0.0
    z.data[IsingAnyon(:σ)] .= 1.0
    z.data[IsingAnyon(:ψ)] .= 0.0
    return z
end

function particle_number_operator(::Type{Z2Space})
    V = Z2Space(0 => 1, 1 => 1)
    z = TensorMap(zeros, Float64, V ← V)
    z.data[ℤ₂(0)] .= 0.0
    z.data[ℤ₂(0)] .= 1.0
    return z
end

function particle_number_operator(::Type{ComplexSpace})
    V = ℂ^2
    z = TensorMap(zeros, Float64, V ← V)
    z.data[1,1] = 0.0
    z.data[2,2] = 1.0
    return z
end

function random_space(::Type{Vect[FibonacciAnyon]}, dlow=2, dhigh=6)
    dtotal = rand(dlow:dhigh)
    d0 = rand(1:dtotal-1)
    d1 = dtotal - d0
    V = Vect[FibonacciAnyon](FibonacciAnyon(:I) => d0, FibonacciAnyon(:τ) => d1)
    return V
end

function random_space(::Type{Vect[IsingAnyon]}, dlow=3, dhigh=6)
    @assert dlow > 2
    dtotal = rand(dlow:dhigh)
    d0 = rand(1:dtotal-2)
    d1 = rand(1:(dtotal-d0-1))
    d2 = dtotal - d0 - d1
    V = Vect[IsingAnyon](
        IsingAnyon(:I) => d0,
        IsingAnyon(:σ) => d1,
        IsingAnyon(:ψ) => d2
    )
    return V
end

function random_space(::Type{Z2Space}, dlow=2, dhigh=6)
    dtotal = rand(dlow:dhigh)
    d0 = rand(1:dtotal-1)
    d1 = dtotal - d0
    V = Z2Space(ℤ₂(0) => d0, ℤ₂(1) => d1)
    return V
end

function random_space(::Type{ComplexSpace}, dlow=1, dhigh=8)
    d = rand(dlow:dhigh)
    V = ℂ^d
    return V
end

"""
Generate vector spaces for a MERA of `M`, with `n` layers, and `S` vector space type.
"""
function random_layerspaces(::Type{S}, ::Type{M}, n, dlow=3, dhigh=6) where {S, M}
    width = causal_cone_width(M)
    V = random_space(S, dlow, dhigh)
    spaces = S[V]
    for i = 2:n
        V_prev = V
        V_prev_fusion = fuse(fill(V_prev, width)...)
        V = random_space(S, dlow, dhigh)
        # If a layer from V_prev to V would be of such a dimension that
        # it couldn't be isometric, try generating another V until
        # success.  Since V = V_prev is a valid choice, this will
        # eventually terminate.
        while infimum(V_prev_fusion, V) != V
            V = random_space(S, dlow, dhigh)
        end
        push!(spaces, V)
    end
    return tuple(spaces...)::NTuple{n,S}
end

"""
Generate vector spaces for the layer-internal indices of a MERA of `M`, that has
`extspaces` as its interlayer spaces.
"""
function random_internalspaces(extspaces, M)
    width = (M == TernaryMERA ? 3 :
    M in (BinaryMERA, ModifiedBinaryMERA) ? 2 : nothing)
    intspaces = ntuple(length(extspaces)) do i
        # The fuse just removes nested product state structure.
        Vext = fuse(extspaces[i])
        Vnext = i == length(extspaces) ? Vext : extspaces[i+1]
        sects = Dict()
        for s in sectors(Vext)
            d = dim(Vext, s)
            rand(Bool) && d > 1 && (d = d - 1)
            sects[s] = d
        end
        Vint = spacetype(Vext)(sects...)
        # Check that we didn't accidentally make the isometric part
        # non-isometric. If we did, just make Vint = Vext.
        Vint_fusion = fuse(fill(Vint, width)...)
        infimum(Vint_fusion, Vnext) != Vnext && (Vint = Vext)
        Vint
    end
    return intspaces
end

"""
Test ModifiedBinaryOp.
"""
function test_modifiedbinaryop(::Type{S}) where {S}
    V = random_space(S, 3)  # Use a minimum dimension of 3.
    # Unary operations
    a = 0.1
    unaryops = [
        copy,
        adjoint,
        one,
        (x) -> MERAKit.expand_support(x, 2),
        (x) -> a * x,
        (x) -> x * a,
        (x) -> x / a,
        (x) -> fill!(x, a),
        (x) -> rmul!(x, a),
        (x) -> lmul!(a, x),
    ]
    if isreal(sectortype(S))
        # These don't make sense complex sector types.
        push!(unaryops, imag)
        push!(unaryops, real)
    end
    for f in unaryops
        mid = TensorMap(randn, ComplexF64, V, V)
        gap = TensorMap(randn, ComplexF64, V, V)
        m = ModifiedBinaryOp(mid, gap)
        res_t = ModifiedBinaryOp(f(mid), f(gap))
        res_m = f(m)
        @test res_t == res_m
    end

    # Binary operations
    a = 0.1
    b = 2.7
    binaryops = (
        *,
        +,
        -,
        /,
        (x, y) -> axpby!(a, x, b, y),
        (x, y) -> axpy!(a, x, y),
        (x, y) -> mul!(x, a, y),
        (x, y) -> mul!(x, y, b),
    )
    for f in binaryops
        t1 = TensorMap(randn, ComplexF64, V, V)
        t2 = TensorMap(randn, ComplexF64, V, V)
        m1 = ModifiedBinaryOp(t1)
        m2 = ModifiedBinaryOp(t2)
        res_tt = ModifiedBinaryOp(f(t1, t2))
        res_mm = f(m1, m2)
        res_mt = f(m1, t2)
        res_tm = f(t1, m2)
        @test res_tt == res_mm
        @test res_tt == res_mt
        @test res_tt == res_tm

        mid1 = TensorMap(randn, ComplexF64, V, V)
        gap1 = TensorMap(randn, ComplexF64, V, V)
        mid2 = TensorMap(randn, ComplexF64, V, V)
        gap2 = TensorMap(randn, ComplexF64, V, V)
        m1 = ModifiedBinaryOp(mid1, gap1)
        m2 = ModifiedBinaryOp(mid2, gap2)
        res_t = ModifiedBinaryOp(f(mid1, mid2), f(gap1, gap2))
        res_m = f(m1, m2)
        @test res_t == res_m
    end

    # Ones that didn't fit the above patterns.
    mid1 = TensorMap(randn, ComplexF64, V, V)
    gap1 = TensorMap(randn, ComplexF64, V, V)
    mid2 = TensorMap(randn, ComplexF64, V, V)
    gap2 = TensorMap(randn, ComplexF64, V, V)
    m1 = ModifiedBinaryOp(mid1, gap1)
    m2 = ModifiedBinaryOp(mid2, gap2)
    res_t = (dot(mid1, mid2) + dot(gap1, gap2)) / 2
    res_m = dot(m1, m2)
    @test res_t == res_m
    res_t = (tr(mid1) + tr(gap1)) / 2
    res_m = tr(m1)
    @test res_t == res_m

    res_t = ModifiedBinaryOp(copy!(mid1, mid2), copy!(gap1, gap2))
    res_m = copy!(m1, m2)
    @test res_t == res_m
end

"""
Test type stability, and type stability only, of various methods.
"""
function test_type_stability(::Type{M}, ::Type{S}) where {M, S}
    layers = 4
    width = @inferred causal_cone_width(M)
    L = layertype(M)
    spaces = random_layerspaces(S, M, layers)
    intspaces = random_internalspaces(spaces, M)
    m = random_MERA(M, ComplexF64, spaces, intspaces; random_disentangler=true)

    V1, V2 = spaces[1:2]
    Vend = spaces[end]
    randomop1 = TensorMap(randn, ComplexF64, V1 ← V1)
    randomop1 = randomop1 + randomop1'
    randomop1 = convert(
        MERAKit.operatortype(m), MERAKit.expand_support(randomop1, causal_cone_width(M))
    )
    randomop2 = TensorMap(randn, ComplexF64, V1 ← V1)
    randomop2 = randomop2 + randomop2'
    randomop2 = convert(
        MERAKit.operatortype(m), MERAKit.expand_support(randomop2, causal_cone_width(M))
    )
    randomrho1 = TensorMap(randn, ComplexF64, Vend ← Vend)
    randomrho1 = randomrho1 + randomrho1'
    randomrho1 = convert(
        MERAKit.operatortype(m), MERAKit.expand_support(randomrho1, causal_cone_width(M))
    )
    randomrho2 = TensorMap(randn, ComplexF64, V2 ← V2)
    randomrho2 = randomrho2 + randomrho2'
    randomrho2 = convert(
        MERAKit.operatortype(m), MERAKit.expand_support(randomrho2, causal_cone_width(M))
    )

    nt = @inferred num_translayers(m)
    l1 = @inferred randomlayer(
        L, ComplexF64, V2, V1, intspaces[1]; random_disentangler=false
    )

    @inferred replace_layer(m, l1, 1)
    @inferred release_transitionlayer(m)
    @inferred projectisometric(m)
    @inferred outputspace(m, Inf)
    @inferred inputspace(m, Inf)
    @inferred MERAKit.environment(randomop1, l1, randomrho2)

    @inferred MERAKit.ascend(randomop1, m)
    @inferred ascended_operator(randomop1, m, 1)
    @inferred MERAKit.descend(randomrho1, m)
    @inferred MERAKit.thermal_densitymatrix(m, Inf)
    @inferred MERAKit.fixedpoint_densitymatrix(m)
    @inferred MERAKit.scale_invariant_operator_sum(randomop1, m, (;))
    @inferred densitymatrix(m, 1)

    # TODO Finish this part, and add @inferred checks for more functions, optimally all
    # functions that can reasonably be expected to be type stable.
    pars = (metric = :euclidean, precondition = false)
    pars = merge(MERAKit.default_pars, pars)
    mtan1 = gradient(randomop1, m, pars)
    mtan2 = gradient(randomop2, m, pars)
    ltan1 = getlayer(mtan1, 1)
    ltan2 = getlayer(mtan2, 1)
    @inferred (x -> tuple(x...))(ltan1)
    @inferred MERAKit.tensorwise_scale(ltan1, 0.1)
    @inferred MERAKit.tensorwise_sum(ltan1, ltan2)
    @inferred inner(l1, ltan1, ltan2)
    @inferred MERAKit.tensorwise_scale(mtan1, 0.1)
    @inferred MERAKit.tensorwise_sum(mtan1, mtan2)
    @inferred inner(m, mtan1, mtan2)
end

"""
For each layer, generate a random operator above and below it (not necessarily Hermitian),
and confirm that ascending the lower or descending the upper one both lead to the same
expectation value (trace of product).
"""
function test_ascend_and_descend(::Type{M}, ::Type{S}) where {M, S}
    layers = 4
    spaces = random_layerspaces(S, M, layers)
    intspaces = random_internalspaces(spaces, M)
    m = random_MERA(M, ComplexF64, spaces, intspaces; random_disentangler=true)
    width = causal_cone_width(M)

    for i in 1:layers
        Vin = inputspace(m, i)
        Vout = outputspace(m, i)
        upper_space = ⊗(fill(Vin, width)...)
        lower_space = ⊗(fill(Vout, width)...)
        randomop1 = TensorMap(randn, ComplexF64, upper_space ← upper_space)
        randomop2 = TensorMap(randn, ComplexF64, lower_space ← lower_space)
        down1 = MERAKit.descend(randomop1, m, i, i+1)
        up2 = MERAKit.ascend(randomop2, m, i+1, i)
        e1 = dot(down1, randomop2)
        e2 = dot(randomop1, up2)
        @test e1 ≈ e2
    end
end

"""
Test that the expectation value of the identity is 1.0, regardless of which layer we
evaluate it at.
"""
function test_expectation_of_identity(::Type{M}, ::Type{S}) where {M, S}
    layers = 4
    spaces = random_layerspaces(S, M, layers)
    intspaces = random_internalspaces(spaces, M)
    m = random_MERA(M, ComplexF64, spaces, intspaces; random_disentangler=true)
    V = outputspace(m, 1)
    eye = id(V)
    for i in 1:(layers+1)
        @test expect(eye, m, (;), 1, i) ≈ 1.0
    end
end

"""
Test that the expectation value of a random Hermitian operator does not depend on the layer
that we evaluate it at.
"""
function test_expectation_evalscale(::Type{M}, ::Type{S}) where {M, S}
    layers = 4
    spaces = random_layerspaces(S, M, layers)
    intspaces = random_internalspaces(spaces, M)
    m = random_MERA(M, ComplexF64, spaces, intspaces; random_disentangler=true)
    V = outputspace(m, 1)
    randomop = TensorMap(randn, ComplexF64, V ← V)
    randomop = (randomop + randomop')/2
    expectation = expect(randomop, m)
    for i in 1:(layers+1)
        expectation_i = expect(randomop, m, (;), 1, i)
        @test expectation_i ≈ expectation
    end
end

"""
Test that pseudoserializing and depseudoserializing back does not change the expectation
value of a random Hermitian operator.
"""
function test_pseudoserialization(::Type{M}, ::Type{S}) where {M, S}
    layers = 4
    spaces = random_layerspaces(S, M, layers)
    intspaces = random_internalspaces(spaces, M)
    m = random_MERA(M, ComplexF64, spaces, intspaces; random_disentangler=true)
    V = outputspace(m, 1)
    randomop = TensorMap(randn, ComplexF64, V ← V)
    randomop = (randomop + randomop')/2
    expectation = expect(randomop, m)
    ser = pseudoserialize(m)
    m_reco = depseudoserialize(ser...)
    new_expectation = expect(randomop, m_reco)
    @test new_expectation ≈ expectation
end

"""
Confirm that expanding bond dimensions does not change the expectation value of a random
Hermitian operator.
"""
function test_expand_bonddim(::Type{M}, ::Type{S}) where {M, S}
    num_layers = 4
    spaces = random_layerspaces(S, M, num_layers, 4)
    intspaces = random_internalspaces(spaces, M)
    m = random_MERA(M, ComplexF64, spaces, intspaces; random_disentangler=true)
    V = outputspace(m, 1)
    randomop = TensorMap(randn, ComplexF64, V ← V)
    randomop = (randomop + randomop')/2
    expectation = expect(randomop, m)
    sf = scalefactor(M)

    # Expand the interlayer spaces, check that expectation value is preserved.
    for i in 1:(num_layers-1)
        V = inputspace(m, i)
        expandable_sectors = sectors(V)
        Vint = fuse(⊗(fill(internalspace(m, i), sf)...))
        expandable_sectors = [s for s in expandable_sectors if dim(V, s) < dim(Vint, s)]
        if i == num_layers-1
            Vint = fuse(⊗(fill(internalspace(m, i+1), sf)...))
            expandable_sectors = [s for s in expandable_sectors if dim(V, s) < dim(Vint, s)]
        end
        newdims = Dict(s => dim(V, s) + 1 for s in expandable_sectors)
        m = expand_bonddim(m, i, newdims)
    end
    new_expectation = expect(randomop, m)
    @test new_expectation ≈ expectation

    # Expand the intralayer spaces, check that expectation value is preserved.
    for i in 1:(num_layers)
        V = internalspace(m, i)
        Vout = outputspace(m, i)
        newdims = Dict(
            s => dim(V, s) < dim(Vout, s) ? dim(V, s) + 1 : dim(V, s)
            for s in sectors(V)
        )
        m = expand_internal_bonddim(m, i, newdims)
    end
    new_expectation = expect(randomop, m)
    @test new_expectation ≈ expectation
end

"""
Confirm that releasing a does not change the expectation value of a random Hermitian
operator.
"""
function test_release_layer(::Type{M}, ::Type{S}) where {M, S}
    layers = 4
    spaces = random_layerspaces(S, M, layers)
    intspaces = random_internalspaces(spaces, M)
    m = random_MERA(M, ComplexF64, spaces, intspaces; random_disentangler=true)
    V = outputspace(m, 1)
    randomop = TensorMap(randn, ComplexF64, V ← V)
    randomop = (randomop + randomop')/2
    expectation = expect(randomop, m)
    m = release_transitionlayer(m)
    new_expectation = expect(randomop, m)
    @test new_expectation ≈ expectation
end

"""
Create a random MERA and operator, evaluate the expectation values, strip both of their
symmetry structure, and confirm that the expectation value hasn't changed.
"""
function test_remove_symmetry(::Type{M}, ::Type{S}) where {M, S}
    layers = 4
    spaces = random_layerspaces(S, M, layers)
    intspaces = random_internalspaces(spaces, M)
    m = random_MERA(M, ComplexF64, spaces, intspaces; random_disentangler=true)
    V = outputspace(m, 1)
    randomop = TensorMap(randn, ComplexF64, V ← V)
    randomop = (randomop + randomop')/2
    expectation = expect(randomop, m)
    m_nosym = remove_symmetry(m)
    randomop_nosym = remove_symmetry(randomop)
    expectation_nosym = expect(randomop_nosym, m_nosym)
    @test expectation ≈ expectation_nosym
end

"""
Create a random MERA and operator, evaluate the expectation values, strip both of their
symmetry structure, and confirm that the expectation value hasn't changed.
"""
function test_reset_storage(::Type{M}, ::Type{S}) where {M, S}
    layers = 4
    spaces = random_layerspaces(S, M, layers)
    intspaces = random_internalspaces(spaces, M)
    m = random_MERA(M, ComplexF64, spaces, intspaces; random_disentangler=true)
    V = outputspace(m, 1)
    randomop = TensorMap(randn, ComplexF64, V ← V)
    randomop = (randomop + randomop')/2
    expectation_orig = expect(randomop, m)

    m = replace_layer(m, getlayer(m, 2), 2)
    expectation_reset = expect(randomop, m)
    @test expectation_orig ≈ expectation_reset

    m = reset_storage(m)
    expectation_reset = expect(randomop, m)
    @test expectation_orig ≈ expectation_reset

    # TODO Consider writing tests that check that a) storage is correctly reset after
    # assigning new tensors, b) storage isn't needlessly reset when assigning tensors, i.e.
    # no unnecessary recomputation is done.
end

"""
Create a MERA that breaks the isometricity condition, by summing up two random MERAs.
Restore isometricity with projectisometric and projectisometric!, and confirm that the
expectation value of the identity is indeed 1 after this.
"""
function test_projectisometric(::Type{M}, ::Type{S}) where {M, S}
    layers = 4
    spaces = random_layerspaces(S, M, layers)
    intspaces = random_internalspaces(spaces, M)
    m1 = random_MERA(M, ComplexF64, spaces, intspaces; random_disentangler=true)
    m2 = random_MERA(M, ComplexF64, spaces, intspaces; random_disentangler=true)
    m = tensorwise_sum(m1, m2)
    V = outputspace(m, 1)
    eye = id(V)
    # Test that we've broken isometricity
    @test !(expect(eye, m) ≈ 1.0)
    # Test that the new MERA is isometric, but that the old one hasn't changed.
    m_new = projectisometric(m)
    @test expect(eye, m_new) ≈ 1.0
    @test !(expect(eye, m) ≈ 1.0)
    # Test that projectisometric! works too. Note that there's no guarantee that it modifies
    # the old m to be isometric. The ! simply gives it permission to mess with old objects.
    m = projectisometric!(m)
    @test expect(eye, m) ≈ 1.0
end

"""
Create a random MERA and compute its scaling dimensions. The numbers will be nonsense, but
at least we know the function runs.
"""
function test_scalingdimensions(::Type{M}, ::Type{S}) where {M, S}
    layers = 4
    spaces = random_layerspaces(S, M, layers)
    intspaces = random_internalspaces(spaces, M)
    m = random_MERA(M, ComplexF64, spaces, intspaces; random_disentangler=true)
    scalingdimensions(m)
end

"""
Test optimization on a Hamiltonian that is just the particle number operator We know what it
should converge to, and it should converge fast.
"""
function test_optimization(::Type{M}, ::Type{S}, method, precondition = false) where {M, S}
    layers = 3
    # eps is the threshold for how close we need to be to the actual ground state energy
    # to pass the test.
    eps = 1e-2
    dlow = 2
    dhigh = 3
    pars = (
        method = method,
        gradient_delta = 1e-6,
        maxiter = 1000,
        isometries_only_iters = 30,
        precondition = precondition,
        verbosity = 0,
        scaleinvariant_krylovoptions = (
            tol = 1e-8,
            verbosity = 0,
            maxiter = 20,
        ),
    )

    op = particle_number_operator(S)
    width = causal_cone_width(M)
    V = domain(op)
    ham = -1*⊗(fill(op, width)...)
    spaces = random_layerspaces(S, M, layers-1, dlow, dhigh)
    spaces = (V, spaces...)
    intspaces = random_internalspaces(spaces, M)
    m = random_MERA(M, ComplexF64, spaces, intspaces)
    m = minimize_expectation(m, ham, pars)
    expectation = expect(ham, m)
    @test abs(expectation + 1.0) < eps
    # TODO We used to check here that these entropies are small, but this doesn't hold for
    # anyons. Figure out what's a good check.
    entropies = densitymatrix_entropies(m)
end

"""
Test gradients and retraction.
"""
function test_gradient_and_retraction(::Type{M}, ::Type{S}, alg, metric) where {M, S}
    layers = 4
    spaces = random_layerspaces(S, M, layers)
    intspaces = random_internalspaces(spaces, M)
    morig = random_MERA(M, ComplexF64, spaces, intspaces; random_disentangler=true)

    width = causal_cone_width(M)
    V = outputspace(morig, 1)
    hamspace = ⊗(fill(V, width)...)
    ham = TensorMap(randn, ComplexF64, hamspace ← hamspace)
    ham = ham + ham'
    eye = id(V)

    pars = (metric = metric, precondition = false)
    pars = merge(MERAKit.default_pars, pars)

    fg(x) = (expect(ham, x), gradient(ham, x, pars))
    scale!(vec, beta) = tensorwise_scale(vec, beta)
    add!(vec1, vec2, beta) = tensorwise_sum(vec1, scale!(vec2, beta))

    # Retract along the trajectory generated by the gradient of ham, by alpha and by
    # alpha+delta.
    alpha = 2.7
    delta = 0.0001
    tanorig = fg(morig)[2]
    m1, tan1 = retract(morig, tanorig, alpha; alg=alg)
    m2, tan2 = retract(morig, tanorig, alpha+delta; alg=alg)

    # Check that the points we retracted to are normalized MERAs.
    @test expect(eye, m1) ≈ 1.0
    @test expect(eye, m2) ≈ 1.0

    # TODO If we implement getting a tangent vector that goes from m1 to m2, we should put
    # in a test that checks that that vector is roughly tan1*delta.

    # Get the energies and gradients at both points. Check that the energy difference
    # between them is the inner product of the gradient and the tangent.
    f1, g1 = fg(m1)
    f2, g2 = fg(m2)
    reco1 = (f2 - f1)/delta
    reco2 = (inner(m1, tan1, g1; metric=metric) + inner(m2, tan2, g2; metric=metric)) / 2.0
    @test isapprox(reco1, reco2; rtol=10*delta)
end

"""
Test vector transport.
"""
function test_transport(::Type{M}, ::Type{S}, alg, metric) where {M, S}
    layers = 4
    spaces = random_layerspaces(S, M, layers)
    intspaces = random_internalspaces(spaces, M)
    m = random_MERA(M, ComplexF64, spaces, intspaces; random_disentangler=true)

    width = causal_cone_width(M)
    V = outputspace(m, 1)
    hamspace = ⊗(fill(V, width)...)
    # Make three different random Hamiltonians.
    hams = [TensorMap(randn, ComplexF64, hamspace ← hamspace) for i in 1:3]
    hams = [ham + ham' for ham in hams]

    pars = (metric = metric, precondition = false)
    pars = merge(MERAKit.default_pars, pars)

    g1, g2, g3 = [gradient(ham, m, pars) for ham in hams]
    angle_pre = inner(m, g2, g3; metric=metric)
    # Transport g2 and g3 along the retraction by g1, by distance alpha.
    alpha = 2.7
    mt = retract(m, g1, alpha; alg=alg)[1]
    g2t = transport!(g2, m, g1, alpha, mt; alg=alg)
    g3t = transport!(g3, m, g1, alpha, mt; alg=alg)
    angle_post = inner(mt, g2t, g3t; metric=metric)

    # Check that the inner product has been preserved by the transport.
    @test angle_pre ≈ angle_post
end

function test_with_all_types(testfunc, meratypes, spacetypes, args...)
    for M in meratypes
        for S in spacetypes
            testfunc(M, S, args...)
        end
    end
end

Random.seed!(2)  # For reproducing the same tests again and again.
meratypes = (ModifiedBinaryMERA, BinaryMERA, TernaryMERA)
spacetypes = (ComplexSpace, Z2Space, Vect[FibonacciAnyon])

# Run the tests on different MERAs and vector spaces.
# Basics
@testset "ModifiedBinaryOp" begin
    for S in spacetypes
        test_modifiedbinaryop(S)
    end
end
@testset "Type stability" begin
    test_with_all_types(test_type_stability, meratypes, spacetypes)
end
@testset "Ascend and descend" begin
    test_with_all_types(test_ascend_and_descend, meratypes, spacetypes)
end
@testset "Expectation of identity" begin
    test_with_all_types(test_expectation_of_identity, meratypes, spacetypes)
end
@testset "Expectation at various evalscales" begin
    test_with_all_types(test_expectation_evalscale, meratypes, spacetypes)
end
@testset "Pseudoserialization" begin
    test_with_all_types(test_pseudoserialization, meratypes, spacetypes)
end
@testset "Expanding bond dimension" begin
    test_with_all_types(test_expand_bonddim, meratypes, spacetypes)
end
@testset "Releasing layers" begin
    test_with_all_types(test_release_layer, meratypes, spacetypes)
end
@testset "Removing symmetry" begin
    # This doesn't make sense for anyons.
    nonanonic_spacetypes = (
        ST for ST in spacetypes if BraidingStyle(sectortype(ST)) != Anyonic()
    )
    test_with_all_types(test_remove_symmetry, meratypes, nonanonic_spacetypes)
end
@testset "Reset storage" begin
    test_with_all_types(test_reset_storage, meratypes, spacetypes)
end
@testset "Projectisometric" begin
    test_with_all_types(test_projectisometric, meratypes, spacetypes)
end
@testset "Scaling dimensions" begin
    # TODO Make this work with anyons too.
    test_with_all_types(test_scalingdimensions, meratypes, (ComplexSpace, Z2Space))
end

# Manifold operations
@testset "Gradient and retraction, Cayley transform, canonical metric" begin
    test_with_all_types(
        test_gradient_and_retraction, meratypes, spacetypes, :cayley, :canonical
    )
end
@testset "Gradient and retraction, Cayley transform, Euclidean metric" begin
    test_with_all_types(
        test_gradient_and_retraction, meratypes, spacetypes, :cayley, :euclidean
    )
end
@testset "Gradient and retraction, exponential, canonical metric" begin
    test_with_all_types(
        test_gradient_and_retraction, meratypes, spacetypes, :exp, :canonical
    )
end
@testset "Gradient and retraction, exponential, Euclidean metric" begin
    test_with_all_types(
        test_gradient_and_retraction, meratypes, spacetypes, :exp, :euclidean
    )
end
@testset "Transport, cayley transform, canonical metric" begin
    test_with_all_types(test_transport, meratypes, spacetypes, :cayley, :canonical)
end
@testset "Transport, cayley transform, Euclidean metric" begin
    test_with_all_types(test_transport, meratypes, spacetypes, :cayley, :euclidean)
end
@testset "Transport, exponential, canonical metric" begin
    test_with_all_types(test_transport, meratypes, spacetypes, :exp, :canonical)
end
@testset "Transport, exponential, Euclidean metric" begin
    test_with_all_types(test_transport, meratypes, spacetypes, :exp, :euclidean)
end

# Optimization
@testset "Optimization E-V" begin
    test_with_all_types((M, S) -> test_optimization(M, S, :ev), meratypes, spacetypes)
end
@testset "Optimization LBFGS" begin
    test_with_all_types(
        (M, S) -> test_optimization(M, S, :lbfgs, false), meratypes, spacetypes
    )
end
@testset "Optimization LBFGS with preconditioning" begin
    test_with_all_types(
        (M, S) -> test_optimization(M, S, :lbfgs, true), meratypes, spacetypes
    )
end


text:
function simple_estimator(model::Potts, T::Real, Js::AbstractArray, _=nothing)
    nsites = numsites(model)
    nbonds = numbonds(model)
    invQ = 1.0/model.Q

    M = 0.0
    @inbounds for s in 1:nsites
        M += ifelse(model.spins[s]==1, 1.0-invQ, -invQ)
    end
    M /= nsites
    E = 0.0
    @inbounds for b in bonds(model)
        s1, s2 = source(b), target(b)
        E -= ifelse(model.spins[s1] == model.spins[s2], 1.0, 0.0) * Js[bondtype(b)]
    end
    E /= nsites

    res = Measurement()
    res["Magnetization"] = M
    res["|Magnetization|"] = abs(M)
    res["Magnetization^2"] = M^2
    res["Magnetization^4"] = M^4
    res["Energy"] = E
    res["Energy^2"] = E^2
    return res
end

@doc raw"""
    improved_estimator(model::Potts, T::Real, Js::AbstractArray, sw::SWInfo)

Returns the following observables as `Dict{String, Any}` using cluster information `sw`

# Observables
- `"Energy"`
    - energy density
- `"Energy^2"`
    - square of energy density
- `"Magnetization"`
    - magnetization density
- `"|Magnetization|"`
    - absolute value of magnetization density
- `"|Magnetization|^2"`
    - square of magnetization density
- `"|Magnetization|^4"`
    - quadruple of magnetization density
- `"Clustersize^2"`
    - ``\sum_c r_c^2 ``, where ``r_c`` is the size density of ``c``-th cluster
- `"Clustersize^4"`
    - ``\sum_c r_c^4 ``
- `"Clustersize^2 Clustersize^2"`
    - ``\sum_{c\ne c'} r_c^2 r_{c'}^2``
"""
function improved_estimator(model::Potts, T::Real, Js::AbstractArray, sw::SWInfo)
    nsites = numsites(model)
    Q = model.Q
    nc = numclusters(sw)
    invV = 1.0/nsites

    # magnetization
    I2 = (Q-1)/(Q*Q)
    I4 = (Q-1)*((Q-1)^3+1)/(Q^5)
    M = 0.0
    M2 = 0.0
    M4 = 0.0
    s2 = -1.0/Q
    s1 = 1.0+s2

    N4 = 0.0
    N2N2 = 0.0
    N2 = 0.0
    for (m,s) in zip(sw.clustersize, sw.clusterspin)
        M += (m*invV) * ifelse(s==1, s1, s2)
        m2 = (m*invV)^2
        N4 += m2*m2
        N2N2 += N2*m2
        N2 += m2
    end
    M4 = I4*N4 + 6*I2*I2*N2N2
    M2 = I2*N2

    # energy
    aJ = abs.(Js)
    mbeta = -1.0/T
    ns = sw.activated_bonds
    As = -aJ ./ expm1.(mbeta.*aJ)
    Ans = ns.*As
    E = 0.0
    E2 = 0.0
    for b in 1:numbondtypes(model)
        E2 += aJ[b]*Ans[b]
        E2 += Ans[b] * As[b]*(ns[b]-1)
        E2 += 2.0*Ans[b]*E
        E += Ans[b]
    end

    E *= -invV
    E2 *= invV*invV

    res = Measurement()
    res["Magnetization"] = M
    res["|Magnetization|"] = abs(M)
    res["Magnetization^2"] = M2
    res["Magnetization^4"] = M4
    res["Energy"] = E
    res["Energy^2"] = E2
    res["Clustersize^2"] = N2
    res["Clustersize^4"] = N4
    res["Clustersize^2 Clustersize^2"] = 2*N2N2

    return res
end

default_estimator(model::Potts, update) = ifelse(update==SW_update!, improved_estimator, simple_estimator)


text:
segy_file_1   = string(path, "/SampleFiles/SEGY/03.334.12.09.00.0362.1")
segy_file_2   = string(path, "/SampleFiles/Restricted/test_rev_1.segy")
segy_file_3   = string(path, "/SampleFiles/SEGY/04.322.18.06.46.92C3.1.coords1.segy")
segy_file_4   = string(path, "/SampleFiles/SEGY/04.322.18.06.46.92C3.1.coords2.segy")
segy_be_file  = string(path, "/SampleFiles/SEGY/02.104.00.00.08.7107.2")
segy_fpat     = string(path, "/SampleFiles/SEGY/02.104.00.00.08.7107.2")

printstyled("  SEG Y\n", color=:light_green)
printstyled("    PASSCAL/NMT SEG Y\n", color=:light_green)
SEG = read_data("passcal", segy_file_1, full=true)

printstyled("      header integrity\n", color=:light_green)

SEG = read_data("passcal", segy_file_1, full=true)
@test SEG.misc[1]["gain_const"] == 32
@test SEG.gain[1] ≈ SEG.misc[1]["scale_fac"]
@test isapprox(1.0/SEG.gain[1], 4.47021e-07/SEG.misc[1]["gain_const"], atol=eps(Float32))
@test SEG.fs[1] == 100.0 == 1.0e6 / SEG.misc[1]["delta"]
@test lastindex(SEG.x[1]) == 247698
@test SEG.misc[1]["trace_seq_line"] == 3
@test SEG.misc[1]["trace_seq_file"] == 3
@test SEG.misc[1]["event_no"] == 1
@test SEG.misc[1]["channel_no"] == 2
@test SEG.misc[1]["trace_id_code"] == 3
@test SEG.misc[1]["h_units_code"] == 2
@test SEG.misc[1]["nx"] == 32767
@test SEG.misc[1]["samp_rate"] == 10000
@test SEG.misc[1]["gain_type"] == 1
@test SEG.misc[1]["year"] == 2003
@test SEG.misc[1]["day"] ==  334
@test SEG.misc[1]["hour"] == 12
@test SEG.misc[1]["minute"] == 9
@test SEG.misc[1]["second"] == 0
@test SEG.misc[1]["ms"] == 5
@test SEG.misc[1]["time_code"] == 2
@test SEG.misc[1]["trigyear"] == 2003
@test SEG.misc[1]["trigday"] == 334
@test SEG.misc[1]["trighour"] == 12
@test SEG.misc[1]["trigminute"] == 9
@test SEG.misc[1]["trigsecond"] == 0
@test SEG.misc[1]["trigms"] == 5
@test SEG.misc[1]["data_form"] == 1
@test SEG.misc[1]["inst_no"] == 0x016a # 0362
@test strip(SEG.misc[1]["sensor_serial"]) == "UNKNOWN"
@test strip(SEG.misc[1]["station_name"]) == "362"

# Location
printstyled("      sensor position\n", color=:light_green)
h_sc = Float64(get(SEG.misc[1], "h_sc", 1.0))
h_sc = abs(h_sc)^(h_sc < 0.0 ? -1 : 1)
z_sc = Float64(get(SEG.misc[1], "z_sc", 1.0))
z_sc = abs(z_sc)^(z_sc < 0.0 ? -1 : 1)
x = get(SEG.misc[1], "rec_x", 0.0)
y = get(SEG.misc[1], "rec_y", 0.0)
z = get(SEG.misc[1], "rec_ele", 0.0)
# @test SEG.loc[1].lat == y*h_sc == 45.2896      # 45.2896 in wash.sta
# @test SEG.loc[1].lon == x*h_sc == -121.7915    # 121.791496 in wash.sta
# @test SEG.loc[1].el == z*z_sc == 1541.0       # 1541.0 in wash.sta

printstyled("      data integrity\n", color=:light_green)
@test Float64(SEG.misc[1]["max"]) == maximum(SEG.x[1]) == 396817
@test Float64(SEG.misc[1]["min"]) == minimum(SEG.x[1]) == -416512
# @test ≈(SEG.x[1][1:10], [47.0, 46.0, 45.0, 44.0, 51.0, 52.0, 57.0, 59.0, 40.0, 34.0])
@test length(SEG.x[1]) == SEG.misc[1]["num_samps"] == 247698

redirect_stdout(out) do
  segyhdr(segy_file_1, passcal=true)
end

printstyled("      big-endian support\n", color=:light_green)

SEG = read_data("passcal", segy_be_file, full=true, swap=true)
@test SEG.n == 1
@test SEG.id ==  ["...spn"] # lol, WHY BOB
@test isapprox(1.0/SEG.gain[1], 5.92346875e-8, atol=eps(Float32))
@test SEG.misc[1]["trigyear"] == SEG.misc[1]["year"] == 2002

printstyled("    wildcard support\n", color=:light_green)

SEG = read_data("passcal", segy_fpat, full=true, swap=true)
@test SEG.n == 1
@test Float64(SEG.misc[1]["max"]) == maximum(SEG.x[1]) == 49295.0
@test Float64(SEG.misc[1]["min"]) == minimum(SEG.x[1]) == -54454.0
@test ≈(SEG.x[1][1:5], [-615.0, -3994.0, -4647.0, -895.0, 190.0])
@test length(SEG.x[1]) == SEG.misc[1]["num_samps"] == 180027

if has_restricted
 printstyled("    SEG Y rev 1\n", color=:light_green)
  SEG = read_data("segy", segy_file_2)
  redirect_stdout(out) do
    segyhdr(segy_file_2)
  end
else
  printstyled("    Skipped SEG Y rev 1\n", color=:red)
end

printstyled("    Helper functions\n", color=:light_green)
printstyled("      auto_coords\n", color=:light_green)
sc = Int16[1,1]
lat0 = -39.4209
lon0 = -71.94061
ele = 2856.3
xy = [round(Int32, lon0*1.0e4), round(Int32, lat0*1.0e4)]
y,x = auto_coords(xy, sc)
@test latlon2xy(y,x) == xy

xy2 = latlon2xy(lat0, lon0)
auto_coords(xy2, sc)
@test isapprox([auto_coords(xy2, sc)...], [lat0, lon0], rtol=eps(Float32))

S1 = read_data(segy_file_3)
S2 = read_data(segy_file_4)
@test abs(S1.loc[1].lat-S2.loc[1].lat) < 1.0e-4
@test abs(S1.loc[1].lon-S2.loc[1].lon) < 1.0e-4
@test abs(S1.loc[1].el-S2.loc[1].el) < 1.0e-4
@test S1.x[1] == S2.x[1]

printstyled("      trid\n", color=:light_green)
@test trid(Int16(2), fs=S1.fs[1]) == "EHZ"
@test trid(Int16(2), fs=S1.fs[1], fc=1/30) == "HHZ"


text:
# not the most robust function, but works in practise
is_kwfunc(f) = (name = string(f); endswith(name, "##kw") || endswith(name, "##kw\""))
is_kwfunc(v::Variable) = is_kwfunc(v._op.val)


text:
module JungleHelperHawk

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/Hawk" Hawk(x::Integer, y::Integer, mainSpeed::Number=100.0, slowerSpeed::Number=80.0, sprite::String="")

const placements = Ahorn.PlacementDict(
    "Hawk (Jungle Helper)" => Ahorn.EntityPlacement(
        Hawk,
        "rectangle",
        Dict{String, Any}(
            "mainSpeed" => 100,
            "slowerSpeed" => 80
        )

    )
)

sprite = "JungleHelper/hawk/hold03"

function Ahorn.selection(entity::Hawk)
    x, y = Ahorn.position(entity)

    res = Ahorn.Rectangle[Ahorn.getSpriteRectangle(sprite, x, y)]
    return res
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Hawk)
    px, py = Ahorn.position(entity)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Hawk, room::Maple.Room)
    x, y = Ahorn.position(entity)
    Ahorn.drawSprite(ctx, sprite, x, y)
end

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
using Test
# using Revise
using PolynomialBasis

PB = PolynomialBasis

function allequal(v1, v2)
    return all(v1 .≈ v2)
end

function f(v)
    x = v[1]
    y = v[2]
    return x^3 + 3y^3 + 2x^2*y + 8x * y
end

function fx(v)
    x = v[1]
    y = v[2]
    return 3x^2 + 4x*y + 8y
end

function fxx(v)
    x = v[1]
    y = v[2]
    return 6*x+4*y
end

function fy(v)
    x = v[1]
    y = v[2]
    return 9y^2 + 2x^2 + 8x
end

function fyy(v)
    x = v[1]
    y = v[2]
    return 18y
end

function fxy(v)
    x = v[1]
    y = v[2]
    return 4x+8
end

function testhessian(v)
    hxx = fxx(v)
    hxy = fxy(v)
    hyy = fyy(v)
    return [hxx,hxy,hyy]
end

function testgradient(v)
    gx = fx(v)
    gy = fy(v)
    return [gx,gy]
end

hermitetp = PB.HermiteTensorProductBasis(2)
interp = PB.InterpolatingPolynomial(1,hermitetp)
points = [-1.0  -1.0   1.0   1.0
          -1.0   1.0  -1.0   1.0]
v = vec(mapslices(f,points,dims=1))
vx = vec(mapslices(fx,points,dims=1))
vy = vec(mapslices(fy,points,dims=1))
vxy = vec(mapslices(fxy,points,dims=1))

coeffs = zeros(16)
coeffs[[1,5,9,13]] .= v
coeffs[[2,6,10,14]] .= vy
coeffs[[3,7,11,15]] .= vx
coeffs[[4,8,12,16]] .= vxy

update!(interp,coeffs)

testp = range(-1.0,stop=1.0,length=4)'
p = PB.tensor_product_points(2,testp)

vals = mapslices(interp,p,dims=1)
testvals = mapslices(f,p,dims=1)
@test allequal(vals,testvals)

grads = mapslices(x->vec(gradient(interp,x)),p,dims=1)
testgrads = mapslices(testgradient,p,dims=1)
@test allequal(grads,testgrads)

hess = mapslices(x->vec(hessian(interp,x)),p,dims=1)
testhess = mapslices(testhessian,p,dims=1)
@test allequal(hess,testhess)


text:
struct Sphere <: Hitable
	center::Point3
	radius::Float64
	material::Material
end

export Sphere, Hit

struct Hit
	p::Point3
	normal::Vec3
	material::Material
	t::Float64
	front_face::Bool
	function Hit(s::Sphere, t, ray)
		p = at(ray, t)
		outward_normal = (p - s.center) / s.radius
		ff = dot(ray.direction, outward_normal)
		norm = ff < 0 ? outward_normal : -outward_normal
		new(p, norm, s.material, t, ff < 0)
	end
	Hit() = new()
end

function trace(sphere::Sphere, ray::Ray, t_min, t_max)
	oc = ray.origin - sphere.center
	a = magnitude(ray.direction)^2
	half_b = dot(oc, ray.direction)
	c = magnitude(oc)^2 - sphere.radius^2
	discriminant = half_b^2 - a*c
	if discriminant < 0
		return -1
	end

	sqrtd = sqrt(discriminant)
	root = (-half_b - sqrtd) / a
	if root < t_min || t_max < root
		root = (-half_b + sqrtd) / a
		if root < t_min || t_max < root
			return -1
		end
	end

	root
end


text:
using Unification
using Logging
using Test

include("test_bindings.jl")
include("test_unify.jl")

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


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "PulseAudio"
version = v"15.0.0"

short_version = "$(version.major).$(version.minor)"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://freedesktop.org/software/pulseaudio/releases/pulseaudio-$short_version.tar.gz", "a570b592351586541daf27b5e4b82555d6ac46bb6920eb847bcf5818e92f4c1e"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
apk update
apk add bash-completion doxygen gettext glib orc-compiler perl-xml-parser
# make sure meson can find everything
sed -i -e "s~c_args = .*~c_args = ['-I${includedir}', '-L${libdir}']~" ${MESON_TARGET_TOOLCHAIN}
# For some reason, librt fails to get linked correctly, so add a flag
sed -i -e "s~c_link_args = .*~c_link_args = ['-lrt']~" ${MESON_TARGET_TOOLCHAIN}
cd pulseaudio-*
# make rpath work with cross compilation
atomic_patch -p2 $WORKSPACE/srcdir/patches/rpath.patch
# disable fastmath
atomic_patch -p2 $WORKSPACE/srcdir/patches/fastmath.patch
# sys/capability.h doesn't seem to be workig on PowerPC
if [[ "${target}" == powerpc64le-* ]]; then
    atomic_patch -p2 $WORKSPACE/srcdir/patches/capabilities.patch
fi
mkdir build
cd build
# optional dependencies I can't build but might be useful
# Avahi Jack LIRC tdb WebRTC
meson ..  -Ddatabase="gdbm" --cross-file=${MESON_TARGET_TOOLCHAIN}
ninja
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(Sys.islinux, supported_platforms())

# The products that we will ensure are always built
products = [
    ExecutableProduct("pulseaudio", :pulseaudio),
    LibraryProduct("module-native-protocol-tcp", :module_native_protocol_tcp, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-combine", :module_combine, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-native-protocol-fd", :module_native_protocol_fd, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-oss", :module_oss, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-card-restore", :module_card_restore, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-cli-protocol-unix", :module_cli_protocol_unix, "lib/pulse-$short_version/modules"),
    LibraryProduct("libpulse-mainloop-glib", :libpulse_mainloop_glib),
    LibraryProduct("module-rescue-streams", :module_rescue_streams, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-filter-heuristics", :module_filter_heuristics, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-simple-protocol-tcp", :module_simple_protocol_tcp, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-raop-sink", :module_raop_sink, "lib/pulse-$short_version/modules"),
    LibraryProduct("libpulsecommon-$short_version", :libpulsecommon, "lib/pulseaudio"),
    ExecutableProduct("pasuspender", :pasuspender),
    LibraryProduct("libprotocol-cli", :libprotocol_cli, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-mmkbd-evdev", :module_mmkbd_evdev, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-position-event-sounds", :module_position_event_sounds, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-virtual-sink", :module_virtual_sink, "lib/pulse-$short_version/modules"),
    LibraryProduct("libpulsecore-$short_version", :libpulsecore, "lib/pulseaudio"),
    LibraryProduct("module-pipe-source", :module_pipe_source, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-null-sink", :module_null_sink, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-loopback", :module_loopback, "lib/pulse-$short_version/modules"),
    LibraryProduct("liboss-util", :liboss_util, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-bluez5-device", :module_bluez5_device, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-sine-source", :module_sine_source, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-always-sink", :module_always_sink, "lib/pulse-$short_version/modules"),
    LibraryProduct("libpulsedsp", :libpulsedsp, "lib/pulseaudio"),
    LibraryProduct("module-device-restore", :module_device_restore, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-equalizer-sink", :module_equalizer_sink, "lib/pulse-$short_version/modules"),
    LibraryProduct("libcli", :libcli, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-console-kit", :module_console_kit, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-role-cork", :module_role_cork, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-role-ducking", :module_role_ducking, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-tunnel-sink", :module_tunnel_sink, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-udev-detect", :module_udev_detect, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-detect", :module_detect, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-tunnel-source", :module_tunnel_source, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-stream-restore", :module_stream_restore, "lib/pulse-$short_version/modules"),
    LibraryProduct("libpulse-simple", :libpulse_simple),
    LibraryProduct("module-augment-properties", :module_augment_properties, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-rtp-send", :module_rtp_send, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-remap-source", :module_remap_source, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-sine", :module_sine, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-combine-sink", :module_combine_sink, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-intended-roles", :module_intended_roles, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-null-source", :module_null_source, "lib/pulse-$short_version/modules"),
    ExecutableProduct("pacmd", :pacmd),
    LibraryProduct("libbluez5-util", :libbluez5_util, "lib/pulse-$short_version/modules"),
    LibraryProduct("librtp", :librtp, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-bluetooth-discover", :module_bluetooth_discover, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-bluetooth-policy", :module_bluetooth_policy, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-echo-cancel", :module_echo_cancel, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-virtual-source", :module_virtual_source, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-switch-on-connect", :module_switch_on_connect, "lib/pulse-$short_version/modules"),
    LibraryProduct("libprotocol-simple", :libprotocol_simple, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-alsa-card", :module_alsa_card, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-tunnel-sink-new", :module_tunnel_sink_new, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-allow-passthrough", :module_allow_passthrough, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-http-protocol-unix", :module_http_protocol_unix, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-filter-apply", :module_filter_apply, "lib/pulse-$short_version/modules"),
    LibraryProduct("libalsa-util", :libalsa_util, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-cli-protocol-tcp", :module_cli_protocol_tcp, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-simple-protocol-unix", :module_simple_protocol_unix, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-pipe-sink", :module_pipe_sink, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-bluez5-discover", :module_bluez5_discover, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-native-protocol-unix", :module_native_protocol_unix, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-rtp-recv", :module_rtp_recv, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-alsa-sink", :module_alsa_sink, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-default-device-restore", :module_default_device_restore, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-virtual-surround-sink", :module_virtual_surround_sink, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-suspend-on-idle", :module_suspend_on_idle, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-rygel-media-server", :module_rygel_media_server, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-always-source", :module_always_source, "lib/pulse-$short_version/modules"),
    LibraryProduct("libpulse", :libpulse),
    LibraryProduct("libraop", :libraop, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-dbus-protocol", :module_dbus_protocol, "lib/pulse-$short_version/modules"),
    ExecutableProduct("pacat", :pacat),
    LibraryProduct("module-match", :module_match, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-alsa-source", :module_alsa_source, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-hal-detect", :module_hal_detect, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-tunnel-source-new", :module_tunnel_source_new, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-volume-restore", :module_volume_restore, "lib/pulse-$short_version/modules"),
    LibraryProduct("libprotocol-http", :libprotocol_http, "lib/pulse-$short_version/modules"),
    ExecutableProduct("pactl", :pactl),
    LibraryProduct("module-switch-on-port-available", :module_switch_on_port_available, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-ladspa-sink", :module_ladspa_sink, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-cli", :module_cli, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-device-manager", :module_device_manager, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-http-protocol-tcp", :module_http_protocol_tcp, "lib/pulse-$short_version/modules"),
    LibraryProduct("libprotocol-native", :libprotocol_native, "lib/pulse-$short_version/modules"),
    LibraryProduct("module-remap-sink", :module_remap_sink, "lib/pulse-$short_version/modules")
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="alsa_jll", uuid="45378030-f8ea-5b20-a7c7-1a9d95efb90e"))
    Dependency(PackageSpec(name="BlueZ_jll", uuid="471b5b61-da80-5748-8755-67d5084d21f2"))
    Dependency(PackageSpec(name="Check_jll", uuid="491db154-c145-5abe-9c32-446728d60cce"))
    Dependency(PackageSpec(name="Dbus_jll", uuid="ee1fde0b-3d02-5ea6-8484-8dfef6360eab"))
    Dependency(PackageSpec(name="eudev_jll", uuid="35ca27e7-8b34-5b7f-bca9-bdc33f59eb06"))
    Dependency(PackageSpec(name="FFTW_jll", uuid="f5851436-0d7a-5f13-b9de-f02708fd171a"))
    Dependency(PackageSpec(name="Gdbm_jll", uuid="54ca2031-c8dd-5cab-9ed4-295edde1660f"))
    Dependency(PackageSpec(name="Glib_jll", uuid="7746bdde-850d-59dc-9ae8-88ece973131d"); compat="2.68.1")
    Dependency(PackageSpec(name="GStreamer_jll", uuid="aaaaf01e-2457-52c6-9fe8-886f7267d736"))
    Dependency(PackageSpec(name="libsndfile_jll", uuid="5bf562c0-5a39-5b4f-b979-f64ac885830c"))
    Dependency(PackageSpec(name="libasyncns_jll", uuid="ed080073-db63-57db-a029-74e11ae80737"))
    Dependency(PackageSpec(name="libcap_jll", uuid="eef66a8b-8d7a-5724-a8d2-7c31ae1e29ed"))
    Dependency(PackageSpec(name="Libiconv_jll", uuid="94ce4f54-9a6c-5748-9c1c-f9c7231a4531"))
    Dependency(PackageSpec(name="Libtool_jll", uuid="a76c16ae-fb8f-5ff0-8826-da3b7a640f0b"))
    Dependency(PackageSpec(name="OpenSSL_jll", uuid="458c3c95-2e84-50aa-8efc-19380b2a3a95"))
    Dependency(PackageSpec(name="SBC_jll", uuid="da37f231-8920-5702-a09a-bdd970cb6ddc"))
    Dependency(PackageSpec(name="SoXResampler_jll", uuid="fbe68eb6-6641-54c6-99e3-f7c7c4d73a57"))
    Dependency(PackageSpec(name="SpeexDSP_jll", uuid="f2f9631b-9a4e-5b48-9975-88f638ec36a7"))

]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Wavelets
using Base.Test

# modified from Base.Test
function test_vecnorm_eq(va, vb, Eps, astr, bstr)
    if length(va) != length(vb)
        error("lengths of ", astr, " and ", bstr, " do not match: ",
              "\n  ", astr, " (length $(length(va))) = ", va,
              "\n  ", bstr, " (length $(length(vb))) = ", vb)
    end
    diff = vecnorm(va - vb)

    if !isnan(Eps) && !(diff <= Eps)
        sdiff = string("|", astr, " - ", bstr, "| <= ", Eps)
        error("assertion failed: ", sdiff,
              "\n  difference = ", diff, " > ", Eps)
    end
end
macro test_vecnorm_eq_eps(a, b, c)
    :(test_vecnorm_eq($(esc(a)), $(esc(b)), $(esc(c)), $(string(a)), $(string(b))))
end

include("util.jl")
include("transforms.jl")
include("threshold.jl")
include("plot.jl")

print("\ntesting: success\n")


text:
using FileIO
if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

@testset "FileIO" begin
	include("query.jl")
	include("loadsave.jl")
	include("error_handling.jl")
end

# make Travis fail when tests fail:


text:
"""
`ExponentiatedKernel([ρ=1])`
The exponentiated kernel is a Mercer kernel given by:
```
    κ(x,y) = exp(ρ²xᵀy)
```
"""
struct ExponentiatedKernel{Tr} <: Kernel{Tr}
    transform::Tr
end
@inline kappa(κ::ExponentiatedKernel, xᵀy::T) where {T<:Real} = exp(xᵀy)

metric(::ExponentiatedKernel) = DotProduct()


text:
function __init__()
    @require StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" include("init_StaticArrays.jl")
end


text:
function Pulsed(; kwargs...)
  cur_dict = Dict{Symbol, AbstractFloat}(
    :N_G => 0.85,
    :kappa => 1.8,
    :epsilon => 0.25,
    :beta_N => 0.026,
    :Q => 26,
    :f_D => 0.85,
    :tau_P => 1.5,
    :sigma_max => 650e6,
    :B_max => 23,
    :sigma_max_hat => 500e6,
    :B_max_hat => 23,
    :P_F => 500,
    :b => 1,
    :nu_n => 0.4,
    :nu_T => 1.1,
    :eta_T => 0.4,
    :eta_RF => 0.75
  )

  merge!(cur_dict, kwargs)

  Reactor(; cur_dict...)
end


text:
module Day16

using Advent: save_input
using Base: Fix1

function parse_file(filename)
    open(filename) do io
        ruledict = Dict{String, Vector{UnitRange}}()
        for line in eachline(io)
            line == "" && break
            rulename, rulevals = split(line, r": ")
            vals = split(rulevals, r" or ")
            ruledict[rulename] = [parse(Int, i):parse(Int, j) for (i, j) in split.(vals, "-")]
        end

        readline(io)
        your = parse.(Int, split(readline(io), ","))
        readline(io)
        readline(io)
        nearby = map(x -> parse.(Int, x), split.(eachline(io), ","))

        ruledict, your, nearby
    end
end

function solve_a(filename)
    ruledict, _, nearby = parse_file(filename)
    rules = Iterators.flatten(values(ruledict))
    sum(n for n in Iterators.flatten(nearby) if all(Fix1(!in, n), rules))
end

function recurse(dict, tocheck)
    isempty(tocheck) && return dict

    if any(isempty, values(dict)) || length(union(values(dict))) < length(dict)
        return empty(dict)
    end

    k = popfirst!(tocheck)
    for j in dict[k]
        save = deepcopy(dict)
        dict[k] = Set(j)
        for (k2, v2) in dict
            k == k2 || setdiff!(v2, dict[k])
        end
        res = recurse(dict, tocheck)
        isempty(res) || return res
        dict = save
    end
    return empty(dict)
end

function solve_b(filename)
    ruledict, your, _nearby = parse_file(filename)
    rules = Iterators.flatten(values(ruledict))
    nearby = filter(ticket -> all(n -> any(Fix1(in, n), rules), ticket), _nearby)

    possible = Dict(
        rulename => Set(
            Iterators.filter(1:length(your)) do j
                all(ticket -> any(Fix1(in, ticket[j]), rulevals), nearby)
            end
        )
        for (rulename, rulevals) in ruledict
    )

    tocheck = sort!(collect(keys(possible)); by=x -> length(possible[x]))
    actual = recurse(possible, tocheck)
    prod(your[only(v)] for (rule, v) in actual if startswith(rule, "departure"))
end

solve(filename = save_input(2020, 16)) = solve_a(filename), solve_b(filename)

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
# Tests related to validated integration

using TaylorModels
# using LinearAlgebra: norm
using Test

# const _num_tests = 1000
# const α_mid = TaylorModels.α_mid

setformat(:full)

# This is based on the examples from Buenger, Numer Algor (2018) 78:1001–1017,
# showing efectivity of shrink-wrapping
@testset "Testing `shrink_wrapping` 1" begin
    _order = 2
    set_variables("x", numvars=1, order=2*_order)
    x0 = IntervalBox(0..0, 1)
    dom = IntervalBox(-1..1, 1)
    for δ in 1/16:1/16:1
        rem = δ * Interval(-1, 1)
        p1 = TaylorModelN(TaylorN(1, order=_order), rem, x0, dom)
        p2 = deepcopy(p1)
        res = p1 * p2
        bound_prod = dom[1]^2 + δ*(2+δ)*dom[1]
        @test res(dom) ⊆ polynomial(res)(dom) + bound_prod
        @test remainder(res) ⊆ bound_prod
        q = [ res ]
        TaylorModels.shrink_wrapping!(q)
        @test q[1](dom) ⊆ res.pol(dom) + bound_prod
        @test remainder(q[1]) ⊆ bound_prod
        @test remainder(q[1]) == 0..0
    end
end

# This is based on the example of Sect 4 of M. Berz, K. Makino,
# in https://bt.pa.msu.edu/pub/papers/VIShrink06/VIShrink06.pdf
@testset "Testing `shrink_wrapping` 2" begin
    function two_state1!(v)
        x, y = v
        x_2 = x^2
        y_2 = y^2
        r_2 = x_2 + y_2
        s = sqrt(1 + r_2)
        v[1] = x * s
        v[2] = y * s
        return v
    end

    function two_state2!(v)
        x1, y1 = v
        x1_2 = x1^2
        y1_2 = y1^2
        r1_2 = x1_2 + y1_2
        r1 = sqrt( 1 + 4*r1_2 )
        s1 = 2 / (1+r1)
        ss = sqrt( s1 )
        v[1] = x1 * ss
        v[2] = y1 * ss
        return v
    end

    local δ = 0.05
    local ib = IntervalBox( -δ .. δ, 2)
    local mib = IntervalBox( 0..0, 2)

    # Diverges using naive Interval arithmetic methods
    ib0 = [ib...]
    for iter = 1:210 # iter is twice the number of iterates
        two_state1!(ib0)
        two_state2!(ib0)
    end
    @test minimum(diam.(ib0)) > 1e6

    # Taylor model controls grow and behaves as identity for after two iterates
    order = 10
    set_variables("x y", order=2*order)
    vm = [TaylorModelN(TaylorN(i, order=order), 0..0, mib, ib) for i = 1:2]
    vm0 = [deepcopy(vm)...]
    for iter = 1:1_000
        two_state1!(vm0)
        two_state2!(vm0)
    end
    @test maximum(remainder.(vm0)) < 1e-6 # 2.2e-7
    @test all(polynomial.(vm0) .== polynomial.(vm)) # only remainder grows!

    # Taylor model with shrink-wrapping after two iterates
    vm0 .= [deepcopy(vm)...]
    for iter = 1:1_000
        two_state1!(vm0)
        two_state2!(vm0)
        TaylorModels.shrink_wrapping!(vm0)
    end
    @test maximum(remainder.(vm0)) < 1e-8 #2.2e-13

    # Taylor model with shrink-wrapping after each iterate
    vm0 .= [deepcopy(vm)...]
    for iter = 1:1_000
        two_state1!(vm0)
        TaylorModels.shrink_wrapping!(vm0)
        two_state2!(vm0)
        TaylorModels.shrink_wrapping!(vm0)
    end
    @test maximum(remainder.(vm0)) < 2.2e-13
end


text:
function symlink_tree(dest::AbstractString, src::AbstractString; verbose::Bool=true)
    for (root, dirs, files) in walkdir(src)
        # Create all directories
        for d in dirs
            d_path = joinpath(root, d)
            dest_dir = joinpath(dest, relpath(root, src), d)
            if ispath(dest_dir) && !isdir(realpath(dest_dir))
                # We can't create a directory if the destination exists and
                # is not a directory or a symlink to a directory.
                error("Directory $(d) from artifact $(basename(src)) already exists on disk and is not a directory")
            end
            mkpath(dest_dir)
        end

        # Symlink all files
        for f in files
            src_file = joinpath(root, f)
            dest_file = joinpath(dest, relpath(root, src), f)
            if isfile(dest_file)
                # Ugh, destination file already exists.  If source and destination files
                # have the same size and SHA256 hash, just move on, otherwise issue a
                # warning.
                if filesize(src_file) == filesize(dest_file)
                    src_file_hash = open(io -> bytes2hex(sha256(io)), src_file, "r")
                    dest_file_hash = open(io -> bytes2hex(sha256(io)), dest_file, "r")
                    if src_file_hash == dest_file_hash
                        continue
                    end
                end

                # Find source artifact that this pre-existent destination file belongs to
                if verbose
                    dest_artifact_source = realpath(dest_file)
                    while occursin("artifacts", dest_artifact_source) && basename(dirname(dest_artifact_source)) != "artifacts"
                        dest_artifact_source = dirname(dest_artifact_source)
                    end
                    @warn("Symlink $(f) from artifact $(basename(src)) already exists in artifact $(basename(dest_artifact_source))")
                end
            else
                # If it's already a symlink, copy over the exact symlink target
                if islink(src_file)
                    symlink(readlink(src_file), dest_file)
                else
                    # Otherwise, point it at the proper location
                    symlink(relpath(src_file, dest_file), dest_file)
                end
            end
        end
    end
end

function unsymlink_tree(dest::AbstractString, src::AbstractString)
    for (root, dirs, files) in walkdir(src)
        # Unsymlink all symlinked directories, non-symlink directories will be culled in audit.
        for d in dirs
            dest_dir = joinpath(dest, relpath(root, src), d)
            if islink(dest_dir)
                rm(dest_dir)
            end
        end

        # Unsymlink all symlinked files
        for f in files
            dest_file = joinpath(dest, relpath(root, src), f)
            if islink(dest_file)
                rm(dest_file)
            end
        end
    end
end

function copy_tree(dest::AbstractString, src::AbstractString; verbose::Bool=true)
    for (root, dirs, files) in walkdir(src)
        # Create all directories
        for d in dirs
            d_path = joinpath(root, d)
            dest_dir = joinpath(dest, relpath(root, src), d)
            if ispath(dest_dir) && !isdir(realpath(dest_dir))
                # We can't create a directory if the destination exists and
                # is not a directory or a symlink to a directory.
                error("Directory $(d) from artifact $(basename(src)) already exists on disk and is not a directory")
            end
            mkpath(dest_dir)
        end

        # Copy all files
        for f in files
            src_file = joinpath(root, f)
            dest_file = joinpath(dest, relpath(root, src), f)
            if isfile(dest_file)
                # Ugh, destination file already exists.  If source and destination files
                # have the same size and SHA256 hash, just move on, otherwise issue a
                # warning.
                if filesize(src_file) == filesize(dest_file)
                    src_file_hash = open(io -> bytes2hex(sha256(io)), src_file, "r")
                    dest_file_hash = open(io -> bytes2hex(sha256(io)), dest_file, "r")
                    if src_file_hash == dest_file_hash
                        continue
                    end
                end

                # Find source artifact that this pre-existent destination file belongs to
                if verbose
                    @warn("File $(f) from $(dirname(src_file)) already exists in $(dest)")
                end
            else
                # If it's already a symlink, copy over the exact symlink target
                cp(src_file, dest_file)
            end
        end
    end
end


text:
using BinDeps
using CMakeWrapper

BinDeps.@setup

version = get(ENV, "LIBPLCTAG_VERSION",  "2.3.5")

plctag = library_dependency("libplctag", aliases = ["plctag"])
if Sys.iswindows()
	bindir = joinpath(BinDeps.depsdir(plctag), plctag.name)
	lib = joinpath(bindir, "Release", "plctag.dll")

	rm(bindir; force = true, recursive = true)

	provides(
		Binaries,
		URI("https://github.com/libplctag/libplctag/releases/download/v$(version)/libplctag_$(version)_windows_x86.zip"),
		plctag,
		installed_libpath = lib,
		unpacked_dir = ".",
	)
else
	srcdir = joinpath(BinDeps.srcdir(plctag), "libplctag-$(version)")
	blddir = joinpath(srcdir, "build")
	lib = joinpath(BinDeps.libdir(plctag), Sys.isapple() ? "libplctag.dylib" : "libplctag.so")

	rm(BinDeps.srcdir(plctag); force = true, recursive = true)
	rm(BinDeps.usrdir(plctag); force = true, recursive = true)

	provides(
		Sources,
		URI("https://github.com/libplctag/libplctag/archive/v$(version).tar.gz"),
		plctag,
		unpacked_dir = srcdir,
	)

	provides(
		BuildProcess,
		@build_steps(begin
			GetSources(plctag)
			CMakeBuild(
				srcdir = srcdir,
				builddir = blddir,
				prefix = BinDeps.usrdir(plctag),
				installed_libpath = [lib],
			)
		end),
		plctag,
	)

	open(joinpath(@__DIR__, "ab_server.jl"), "w+") do io
		println(io, "const AB_SERVER_BIN = $(repr(joinpath(blddir, "bin_dist", "ab_server")))")
	end
end

BinDeps.@install Dict(:plctag => :_plctag)


text:
using Test
using DFTK
import DFTK: Applyχ0Model

include("testcases.jl")

@testset "Compare different SCF algorithms (no spin, no temperature)" begin
    Ecut = 3
    n_bands = 6
    fft_size = [9, 9, 9]
    tol = 1e-7

    Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.psp))
    model = model_DFT(silicon.lattice, [Si => silicon.positions], [:lda_xc_teter93])
    basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)

    # Run nlsolve without guess
    ρ0 = RealFourierArray(basis)  # RFA of zeros
    ρ_nl = self_consistent_field(basis; ρ=ρ0, tol=tol).ρ.fourier

    # Run DM
    if mpi_nprocs() == 1  # Distributed implementation not yet available
        @testset "Direct minimization" begin
            ρ_dm = direct_minimization(basis; g_tol=tol).ρ.fourier
            @test maximum(abs.(ρ_dm - ρ_nl)) < sqrt(tol) / 10
        end
    end

    # Run other SCFs with SAD guess
    ρ0 = guess_density(basis, [Si => silicon.positions])
    for solver in (scf_nlsolve_solver(), scf_damping_solver(1.2), scf_anderson_solver(),
                   scf_CROP_solver())
        @testset "Testing $solver" begin
            ρ_alg = self_consistent_field(basis; ρ=ρ0, solver=solver, tol=tol).ρ.fourier
            @test maximum(abs.(ρ_alg - ρ_nl)) < sqrt(tol) / 10
        end
    end

    @testset "Enforce symmetry" begin
        ρ_alg = self_consistent_field(basis; ρ=ρ0, tol=tol, enforce_symmetry=true).ρ.fourier
        @test maximum(abs.(ρ_alg - ρ_nl)) < sqrt(tol) / 10
    end

    # Run other mixing with default solver (the others are too slow...)
    for mixing in (KerkerMixing(), SimpleMixing(), SimpleMixing(.5), DielectricMixing(εr=12),
                   KerkerDosMixing(), HybridMixing(), HybridMixing(εr=10, RPA=false))
        @testset "Testing $mixing" begin
            ρ_alg = self_consistent_field(basis; ρ=ρ0, mixing=mixing, tol=tol).ρ.fourier
            @test maximum(abs.(ρ_alg - ρ_nl)) < sqrt(tol) / 10
        end
    end
end

@testset "Compare different SCF algorithms (no spin, temperature)" begin
    Ecut = 3
    n_bands = 6
    fft_size = [9, 9, 9]
    tol = 1e-7

    Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.psp))
    model = model_DFT(silicon.lattice, [Si => silicon.positions], [:lda_xc_teter93],
                      temperature=0.01, smearing=Smearing.Gaussian())
    basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)

    # Reference: Default algorithm
    ρ0    = guess_density(basis)
    ρ_ref = self_consistent_field(basis, ρ=ρ0, tol=tol).ρ.fourier

    for mixing in (KerkerDosMixing(), HybridMixing(RPA=true), HybridMixing(εr=10, RPA=false), )
        @testset "Testing $mixing" begin
            ρ_mix = self_consistent_field(basis; ρ=ρ0, mixing=mixing, tol=tol).ρ.fourier
            @test maximum(abs.(ρ_mix - ρ_ref)) < sqrt(tol)
        end
    end
end

@testset "Compare different SCF algorithms (collinear spin, temperature)" begin
    Ecut = 11
    n_bands = 8
    fft_size = [13, 13, 13]
    tol = 1e-7

    Fe = ElementPsp(iron_bcc.atnum, psp=load_psp(iron_bcc.psp))
    magnetic_moments = [Fe => [4.0]]
    model = model_LDA(iron_bcc.lattice, [Fe => iron_bcc.positions],
                      temperature=0.01, magnetic_moments=magnetic_moments,
                      spin_polarization=:collinear)
    basis = PlaneWaveBasis(model, Ecut; fft_size=fft_size, kgrid=[3, 3, 3])
    ρspin0 = guess_spin_density(basis, magnetic_moments)

    # Reference: Default algorithm
    ρ0     = guess_density(basis)
    scfres = self_consistent_field(basis, ρ=ρ0, ρspin=ρspin0, tol=tol)
    ρspin_ref = scfres.ρspin.fourier
    ρ_ref     = scfres.ρ.fourier

    for mixing in (KerkerMixing(), KerkerDosMixing(α=1.0), DielectricMixing(εr=10),
                   HybridMixing(εr=10), χ0Mixing(χ0terms=[Applyχ0Model()], RPA=false),)
        @testset "Testing $mixing" begin
            scfres = self_consistent_field(basis; ρ=ρ0, ρspin=ρspin0, mixing=mixing, tol=tol)
            @test maximum(abs.(scfres.ρ.fourier     - ρ_ref    )) < sqrt(tol)
            @test maximum(abs.(scfres.ρspin.fourier - ρspin_ref)) < sqrt(tol)
        end
    end
end


text:
using FrankWolfe
using LinearAlgebra
using Test
using SparseArrays

@testset "Testing Blended Pairwise Conditional Gradients" begin
    f(x) = norm(x)^2
    function grad!(storage, x)
        @. storage = 2x
        return nothing
    end
    lmo_prob = FrankWolfe.ProbabilitySimplexOracle(4)
    x0 = FrankWolfe.compute_extreme_point(lmo_prob, zeros(10))
    res_bpcg = FrankWolfe.blended_pairwise_conditional_gradient(
        f,
        grad!,
        lmo_prob,
        x0,
        max_iteration=6000,
        line_search=FrankWolfe.Adaptive(),
        verbose=true,
        epsilon=3e-7,
    )
    res_afw = FrankWolfe.away_frank_wolfe(
        f,
        grad!,
        lmo_prob,
        x0,
        max_iteration=6000,
        line_search=FrankWolfe.Adaptive(),
        print_iter=100,
        verbose=true,
        epsilon=3e-7,
    )
    @test res_afw[3] ≈ res_bpcg[3]
    @test norm(res_afw[1] - res_bpcg[1]) ≈ 0 atol=1e-6
    res_bpcg2 = FrankWolfe.blended_pairwise_conditional_gradient(
        f,
        grad!,
        lmo_prob,
        x0,
        max_iteration=6000,
        line_search=FrankWolfe.Adaptive(),
        verbose=false,
        lazy=true
    )
    @test res_bpcg2[3] ≈ res_bpcg[3] atol=1e-5
    active_set_afw = res_afw[end]
    storage = copy(active_set_afw.x)
    grad!(storage, active_set_afw.x)
    epsilon = 1e-6
    @inferred FrankWolfe.afw_step(active_set_afw.x, storage, lmo_prob, active_set_afw, epsilon)
    @inferred FrankWolfe.lazy_afw_step(active_set_afw.x, storage, lmo_prob, active_set_afw, epsilon, epsilon)
end