text:
module LP

import ..FileFormats
import VecMathOptInterface
const MOI = VecMathOptInterface

# Julia 1.6 removes Grisu from Base. Previously, we went
#   print_shortest(io, x) = Base.Grisu.print_shortest(io, x)
# To avoid adding Grisu as a dependency, use the following printing heuristic.
# TODO(odow): consider printing 1.0 as 1.0 instead of 1, i.e., without the
# rounding branch.
function print_shortest(io::IO, x::Real)
    x_int = round(Int, x)
    if isapprox(x, x_int)
        print(io, x_int)
    else
        print(io, x)
    end
    return
end

MOI.Utilities.@model(
    Model,
    (MOI.ZeroOne, MOI.Integer),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
    (),
    (),
    (),
    (MOI.ScalarAffineFunction,),
    (),
    ()
)
function MOI.supports(
    ::Model{T},
    ::MOI.ObjectiveFunction{<:MOI.ScalarQuadraticFunction{T}},
) where {T}
    return false
end

struct Options
    maximum_length::Int
    warn::Bool
end

function get_options(m::Model)
    default_options = Options(255, false)
    return get(m.ext, :LP_OPTIONS, default_options)
end

"""
    Model(; kwargs...)

Create an empty instance of FileFormats.LP.Model.

Keyword arguments are:

 - `maximum_length::Int=255`: the maximum length for the name of a variable.
   lp_solve 5.0 allows only 16 characters, while CPLEX 12.5+ allow 255.

 - `warn::Bool=false`: print a warning when variables or constraints are renamed.
"""
function Model(; maximum_length::Int = 255, warn::Bool = false)
    model = Model{Float64}()
    options = Options(maximum_length, warn)
    model.ext[:LP_OPTIONS] = options
    return model
end

function Base.show(io::IO, ::Model)
    print(io, "A .LP-file model")
    return
end

# ==============================================================================
#
#   Base.write
#
# ==============================================================================

const START_REG = r"^([\.0-9eE])"
const NAME_REG = r"([^a-zA-Z0-9\!\"\#\$\%\&\(\)\/\,\.\;\?\@\_\`\'\{\}\|\~])"

function write_function(
    io::IO,
    model::Model,
    func::MOI.VariableIndex,
    variable_names::Dict{MOI.VariableIndex,String},
)
    print(io, variable_names[func])
    return
end

function write_function(
    io::IO,
    model::Model,
    func::MOI.ScalarAffineFunction{Float64},
    variable_names::Dict{MOI.VariableIndex,String},
)
    is_first_item = true
    if !(func.constant ≈ 0.0)
        print_shortest(io, func.constant)
        is_first_item = false
    end
    for term in func.terms
        if !(term.coefficient ≈ 0.0)
            if is_first_item
                print_shortest(io, term.coefficient)
                is_first_item = false
            else
                print(io, term.coefficient < 0 ? " - " : " + ")
                print_shortest(io, abs(term.coefficient))
            end

            print(io, " ", variable_names[term.variable])
        end
    end
    return
end

function write_constraint_suffix(io::IO, set::MOI.LessThan)
    print(io, " <= ")
    print_shortest(io, set.upper)
    println(io)
    return
end

function write_constraint_suffix(io::IO, set::MOI.GreaterThan)
    print(io, " >= ")
    print_shortest(io, set.lower)
    println(io)
    return
end

function write_constraint_suffix(io::IO, set::MOI.EqualTo)
    print(io, " = ")
    print_shortest(io, set.value)
    println(io)
    return
end

function write_constraint_suffix(io::IO, set::MOI.Interval)
    print(io, " <= ")
    print_shortest(io, set.upper)
    println(io)
    return
end

function write_constraint_prefix(io::IO, set::MOI.Interval)
    print_shortest(io, set.lower)
    print(io, " <= ")
    return
end

write_constraint_prefix(io::IO, set) = nothing

function write_constraint(
    io::IO,
    model::Model,
    index::MOI.ConstraintIndex,
    variable_names::Dict{MOI.VariableIndex,String};
    write_name::Bool = true,
)
    func = MOI.get(model, MOI.ConstraintFunction(), index)
    set = MOI.get(model, MOI.ConstraintSet(), index)
    if write_name
        print(io, MOI.get(model, MOI.ConstraintName(), index), ": ")
    end
    write_constraint_prefix(io, set)
    write_function(io, model, func, variable_names)
    return write_constraint_suffix(io, set)
end

const SCALAR_SETS = (
    MOI.LessThan{Float64},
    MOI.GreaterThan{Float64},
    MOI.EqualTo{Float64},
    MOI.Interval{Float64},
)

function write_sense(io::IO, model::Model)
    if MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        println(io, "maximize")
    else
        println(io, "minimize")
    end
    return
end

function write_objective(
    io::IO,
    model::Model,
    variable_names::Dict{MOI.VariableIndex,String},
)
    print(io, "obj: ")
    obj_func_type = MOI.get(model, MOI.ObjectiveFunctionType())
    obj_func = MOI.get(model, MOI.ObjectiveFunction{obj_func_type}())
    write_function(io, model, obj_func, variable_names)
    println(io)
    return
end

"""
    Base.write(io::IO, model::FileFormats.LP.Model)

Write `model` to `io` in the LP file format.
"""
function Base.write(io::IO, model::Model)
    options = get_options(model)
    FileFormats.create_unique_names(
        model,
        warn = options.warn,
        replacements = [
            s -> match(START_REG, s) !== nothing ? "_" * s : s,
            s -> replace(s, NAME_REG => "_"),
            s -> s[1:min(length(s), options.maximum_length)],
        ],
    )
    variable_names = Dict{MOI.VariableIndex,String}(
        index => MOI.get(model, MOI.VariableName(), index) for
        index in MOI.get(model, MOI.ListOfVariableIndices())
    )
    write_sense(io, model)
    write_objective(io, model, variable_names)
    println(io, "subject to")
    for S in SCALAR_SETS
        for index in MOI.get(
            model,
            MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64},S}(),
        )
            write_constraint(
                io,
                model,
                index,
                variable_names;
                write_name = true,
            )
        end
    end
    println(io, "Bounds")
    free_variables = Set(keys(variable_names))
    for S in SCALAR_SETS
        for index in
            MOI.get(model, MOI.ListOfConstraintIndices{MOI.VariableIndex,S}())
            delete!(free_variables, MOI.VariableIndex(index.value))
            write_constraint(
                io,
                model,
                index,
                variable_names;
                write_name = false,
            )
        end
    end
    for index in MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.ZeroOne}(),
    )
        delete!(free_variables, MOI.VariableIndex(index.value))
    end
    for variable in sort(collect(free_variables), by = x -> x.value)
        println(io, variable_names[variable], " free")
    end
    for (S, str_S) in [(MOI.Integer, "General"), (MOI.ZeroOne, "Binary")]
        indices =
            MOI.get(model, MOI.ListOfConstraintIndices{MOI.VariableIndex,S}())
        if length(indices) > 0
            println(io, str_S)
            for index in indices
                write_function(
                    io,
                    model,
                    MOI.get(model, MOI.ConstraintFunction(), index),
                    variable_names,
                )
                println(io)
            end
        end
    end
    println(io, "End")
    return
end

# ==============================================================================
#
#   Base.read!
#
# ==============================================================================

function Base.read!(io::IO, model::Model)
    return error("read! is not implemented for LP files.")
end

end


text:
using Test, ReferenceTests, BSON

# make sure to turn off GPU usage, at least for Github Actions

include("part1.jl")

ext: jl


text:
using IncrementalInference
using Test

# import IncrementalInference: buildCliquePotentials

global N=100
global fg = initfg()

# doors = [-100.0;0.0;100.0;300.0]'
# cov = [3.0]

addVariable!(fg,:x1, ContinuousScalar, N=N)
addFactor!(fg, [:x1;], Prior(Normal()))

addVariable!(fg,:x2, ContinuousScalar, N=N)
addFactor!(fg,[:x1, :x2], LinearConditional(Normal()))

addVariable!(fg, :x3, ContinuousScalar, N=N)
addFactor!(fg,[:x2,:x3],LinearConditional(Normal()))

addVariable!(fg, :l1, ContinuousScalar, N=N)
addFactor!(fg, [:x1,:l1], LinearConditional(Normal()) )
addFactor!(fg, [:x2,:l1], LinearConditional(Normal()) )

addVariable!(fg, :l2, ContinuousScalar, N=N)
addFactor!(fg, [:x3,:l2], LinearConditional(Normal()))

@testset "test building tree native" begin

global fg

resetFactorGraphNewTree!(fg)

# p = getEliminationOrder(fg, ordering=:qr)
p = [:l1;:l2;:x1;:x2;:x3]

println()
global fge = deepcopy(fg)
println("Building Bayes net...")
buildBayesNet!(fge, p)

global tree = emptyBayesTree()
buildTree!(tree, fge, p)

@test num_vertices(tree.bt) == 3

# Michael reference -- x2->x1, x2->x3, x2->x4, x2->l1, x4->x3, l1->x3, l1->x4
# if false
#   println("Bayes Tree")
#   # Graphs.plot(tree.bt)
#   fid = open("bt.dot","w+")
#   write(fid,Graphs.to_dot(tree.bt))
#   close(fid)
#   run(`dot bt.dot -Tpdf -o bt.pdf`)
# end

println("Find potential functions for each clique")
cliq = tree.cliques[1] # start at the root
buildCliquePotentials(fg, tree, cliq); # fg does not have the marginals as fge does

end

# TODO -- add testing to ensure this is the correct tree!

@warn "add test tree verification"

# run(`evince bt.pdf`)

@testset "build tree from ordering" begin

global fg

resetFactorGraphNewTree!(fg)
vo = getEliminationOrder(fg)
tree = buildTreeFromOrdering!(fg, vo)

end

#


text:
# if Sys.islinux()
#     @testset "mpi_SparseTensor" begin
#         A = sprand(10,10,0.3)
#         B = mpi_SparseTensor(A)
#         @test run(sess, B) ≈ A

#         B = mpi_SparseTensor(A+5I)
#         ADCME.options.mpi.solver = "GMRES"
#         g = rand(10)
#         @test run(sess, B\g) ≈ (A+5I)\g
#     end
# end


text:
"""
The pinched hemisphere benchmark for the configuration with a 18 deg hole at the
top

The reference below states:

The spherical shell shown in Fig. 9 is our proposed doubly-curved shell
problem. Note that the equator is a free edge so that the problem
represents a hemisphere with four point loads alternating in sign at 90 °
intervals on the equator. The hole at the top has been introduced to
avoid the use of triangles near the axis of revolution. Convergence can
be studied by varying  mesh size. Both membrane and bending strains
contribute significantly to the radial displacement at the load point. A
theoretical value of the displacement under load has been computed for a
slightly different configuration [7] in which the hole at the axis is
closed.

Macneal RH, Harder RL (1985) A proposed standard set of problems to test
finite element accuracy. Finite Elements in Analysis and Design 1: 3-20.

Regarding the convergence: refer to the article
Performance of the MITC3+ and MITC4+ shell elements in widely-used
benchmark problems
Yeongbin Ko a, Youngyu Lee b, Phill-Seung Lee a,⇑, Klaus-Jürgen Bathe

The drilling degrees of freedom stiffness has an effect on the convergence for
coarser meshes.

Abaqus KL eight-node curved element gives this sequence (element size changed
with a factor of two): 9.345, 9.353, 9.357. Extrapolated value 9.3597. (/100)
"""
module hemisphere_open_examples

using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FESetShellQ4Module: FESetShellQ4
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using FinEtools.MeshExportModule.VTKWrite: vtkwrite

using Infiltrator

function spherical!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    r = vec(XYZ);
    csmatout[:, 3] .= vec(r)/norm(vec(r))
    csmatout[:, 2] .= (0.0, 0.0, 1.0)
    cross3!(view(csmatout, :, 1), view(csmatout, :, 2), view(csmatout, :, 3))
    csmatout[:, 1] .= vec(view(csmatout, :, 1))/norm(vec(view(csmatout, :, 1)))
    cross3!(view(csmatout, :, 2), view(csmatout, :, 3), view(csmatout, :, 1))
    return csmatout
end

function _execute(n = 8, visualize = true, exact_normals = false, drilling_stiffness_multiplier = 1.0)
    E = 6.825e7;
    nu = 0.3;
    thickness  =  0.04;
    # analytical solution for the vertical deflection under the load
    analyt_sol = 0.093;
    R = 10.0;
    formul = FEMMShellT3FFModule

    tolerance = R/n/100
    fens, fes = Q4block(90.0, 70.0, n, n)
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        phi, psi = fens.xyz[i, 1:2]
        fens.xyz[i, :] .= (cos(psi/180*pi) .* cos(phi/180*pi)*R, cos(psi/180*pi) .* sin(phi/180*pi)*R, sin(psi/180*pi) .* R)
    end
    fens, fes = Q4toT3(fens, fes)

    vtkwrite("geom.vtu", fens, fes)

    ocsys = CSys(3, 3, spherical!)

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

    # Report
    @info "Mesh: $n elements per side"

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    if exact_normals
        femm = formul.make(IntegDomain(fes, TriRule(1), thickness), ocsys, mater)
    else
        femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    end
    femm.drilling_stiffness_scale = 0.1 * drilling_stiffness_multiplier
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # plane of symmetry perpendicular to X
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [1,5,6]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to Y
    l1 = selectnode(fens; box = Float64[-Inf Inf 0 0 -Inf Inf], inflate = tolerance)
    for i in [2,4,6]
        setebc!(dchi, l1, true, i)
    end
    # one point from the list
    for i in [3]
        setebc!(dchi, l1[1:1], true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Load
    nl = selectnode(fens; box = Float64[0 0 R R 0 0], inflate = tolerance)
    loadbdry = FESetP1(reshape(nl, 1, 1))
    lfemm = FEMMBase(IntegDomain(loadbdry, PointRule()))
    fi = ForceIntensity(FFlt[0, -1, 0, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 3);
    nl = selectnode(fens; box = Float64[R R 0 0 0 0], inflate = tolerance)
    loadbdry = FESetP1(reshape(nl, 1, 1))
    lfemm = FEMMBase(IntegDomain(loadbdry, PointRule()))
    fi = ForceIntensity(FFlt[1, 0, 0, 0, 0, 0]);
    F += distribloads(lfemm, geom0, dchi, fi, 3);

    # @infiltrate
    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    resultpercent =  dchi.values[nl, 1][1]*100
    @info "Solution: $(round(resultpercent/analyt_sol, digits = 4))%"

    # Generate a graphical display of resultants

    # Visualization
    if visualize

        scalars = []
        for nc in 1:3
            fld = fieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
            push!(scalars, ("m$nc", fld.values))
        end
        vtkwrite("hemisphere_open-$n-m.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
        scalars = []
        for nc in 1:3
            fld = fieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
            push!(scalars, ("n$nc", fld.values))
        end
        vtkwrite("hemisphere_open-$n-n.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
        scalars = []
        for nc in 1:2
            fld = fieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
            push!(scalars, ("q$nc", fld.values))
        end
        vtkwrite("hemisphere_open-$n-q.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])

        vtkwrite("hemisphere_open-$n-normals.vtu", fens, fes; vectors = [("normals", femm._normals[:, 1:3])])

        scattersysvec!(dchi, (R/4)/maximum(abs.(U)).*U)
        update_rotation_field!(Rfield0, dchi)
        plots = cat(plot_space_box([[0 0 -R]; [R R R]]),
        #plot_nodes(fens),
            plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
            dims = 1)
        pl = render(plots)
    end
    return resultpercent/analyt_sol
end

function test_convergence(ns = [4, 8, 16,] )
    @info "Hemisphere with opening"

    results = Float64[]
    for n in ns
        v = _execute(n, true)
        push!(results, v)
    end
    return ns, results
end

function test_convergence_normals(ns = [4, 8, 16,], exact_normals = false, drilling_stiffness_multiplier =  1.0)
    @info "Hemisphere with opening"

    results = Float64[]
    for n in ns
        v = _execute(n, false, exact_normals, drilling_stiffness_multiplier)
        push!(results, v)
    end
    return ns, results
end

end # module

using PGFPlotsX

let
    using .hemisphere_open_examples
    ns, results = hemisphere_open_examples.test_convergence()

    objects = []

    Allman = [
    20.33
    37.68
    87.11
    ]
    Cook_flat_stiffened = [
    66.22
    94.55
    98.22
    ]
    Providas_Kattis = [
    99.90
    99.06
    98.37
    ]

    all_results = [("Present", results, "*"), ("Allman", Allman, "x"), ("Providas, Kattis", Providas_Kattis, "triangle"), ("Cook, flat, stiffened", Cook_flat_stiffened, "square")]

    for r in  all_results
        @pgf p = PGFPlotsX.Plot(
        {
        color = "black",
        line_width  = 0.7,
        style = "solid",
        mark = "$(r[3])"
        },
        Coordinates([v for v in  zip(ns, r[2])])
        )
        push!(objects, p)
        push!(objects, LegendEntry("$(r[1])"))
    end

    @pgf ax = Axis(
    {
    xlabel = "Number of Elements / side [ND]",
    ylabel = "Normalized Displacement [ND]",
        # xmin = range[1],
        # xmax = range[2],
    xmode = "linear",
    ymode = "linear",
    yminorgrids = "true",
    grid = "both",
    legend_style = {
    at = Coordinate(0.5, 1.05),
    anchor = "south",
    legend_columns = -1
    },
    },
    objects...
    )

    display(ax)
    pgfsave("hemisphere_open-convergence.pdf", ax)
end

let
    test_convergence_normals = hemisphere_open_examples.test_convergence_normals
    ns = [4, 8, 16, 32, 64, 128]

    objects = []

    ns, results = test_convergence_normals(ns, false, 0.1)
    @pgf p = PGFPlotsX.Plot(
    {
    color = "black",
    line_width  = 0.7,
    style = "solid",
    mark = "x"
    },
    Coordinates([v for v in  zip(ns, results)])
    )
    push!(objects, p)
    push!(objects, LegendEntry("0.1"))

    ns, results = test_convergence_normals(ns, false, 1.0)
    @pgf p = PGFPlotsX.Plot(
    {
    color = "black",
    line_width  = 0.7,
    style = "solid",
    mark = "o"
    },
    Coordinates([v for v in  zip(ns, results)])
    )
    push!(objects, p)
    push!(objects, LegendEntry("1.0"))

    ns, results = test_convergence_normals(ns, false, 100.0)
    @pgf p = PGFPlotsX.Plot(
    {
    color = "black",
    line_width  = 0.7,
    style = "solid",
    mark = "diamond"
    },
    Coordinates([v for v in  zip(ns, results)])
    )
    push!(objects, p)
    push!(objects, LegendEntry("100.0"))

    @pgf ax = Axis(
    {
    xlabel = "Number of Elements / side [ND]",
    ylabel = "Normalized Displacement [ND]",
        # xmin = range[1],
        # xmax = range[2],
    xmode = "linear",
    ymode = "linear",
    yminorgrids = "true",
    grid = "both",
    legend_style = {
    at = Coordinate(0.5, 1.05),
    anchor = "south",
    legend_columns = -1
    },
    },
    objects...
    )

    display(ax)
    pgfsave("hemisphere_open-normals-approximate-convergence.pdf", ax)
end

let
    test_convergence_normals = hemisphere_open_examples.test_convergence_normals
    ns = [4, 8, 16, 32, 64, 128]

    objects = []

    ns, results = test_convergence_normals(ns, true, 0.1)
    @pgf p = PGFPlotsX.Plot(
    {
    color = "black",
    line_width  = 0.7,
    style = "solid",
    mark = "x"
    },
    Coordinates([v for v in  zip(ns, results)])
    )
    push!(objects, p)
    push!(objects, LegendEntry("0.1"))

    ns, results = test_convergence_normals(ns, true, 1.0)
    @pgf p = PGFPlotsX.Plot(
    {
    color = "black",
    line_width  = 0.7,
    style = "solid",
    mark = "o"
    },
    Coordinates([v for v in  zip(ns, results)])
    )
    push!(objects, p)
    push!(objects, LegendEntry("1.0"))

    ns, results = test_convergence_normals(ns, true, 100.0)
    @pgf p = PGFPlotsX.Plot(
    {
    color = "black",
    line_width  = 0.7,
    style = "solid",
    mark = "diamond"
    },
    Coordinates([v for v in  zip(ns, results)])
    )
    push!(objects, p)
    push!(objects, LegendEntry("100.0"))

    @pgf ax = Axis(
    {
    xlabel = "Number of Elements / side [ND]",
    ylabel = "Normalized Displacement [ND]",
        # xmin = range[1],
        # xmax = range[2],
    xmode = "linear",
    ymode = "linear",
    yminorgrids = "true",
    grid = "both",
    legend_style = {
    at = Coordinate(0.5, 1.05),
    anchor = "south",
    legend_columns = -1
    },
    },
    objects...
    )

    display(ax)
    pgfsave("hemisphere_open-normals-exact-convergence.pdf", ax)
end


text:
import CVRPSEP_jll
const LIB_CVRPSEP = CVRPSEP_jll.libcvrpsep

include("types.jl")
include("cnstrmgr.jl")
include("rounded_capacity.jl")
include("strengthened_comb.jl")
include("framed_capacity.jl")
include("multistar.jl")
include("hypotour.jl")

function input_conversion(demand::Vector{Int}, edge_tail::Vector{Int}, edge_head::Vector{Int}, edge_x::Vector{Float64})
        n_edges = length(edge_tail)

        # In CVRPSEP, customers are 1, 2, ..., n_customers.
        # Since C starts with index 0, we need to fill the first element.
        # The input vector `demand` should already contain
        _demand = demand
        n_customers = length(demand) - 1

        # make copies
        _edge_tail = copy(edge_tail)
        _edge_head = copy(edge_head)
        _edge_x = copy(edge_x)

        # In Julia,   customers are 2, 3, ..., n_customers+1
        # In CVRPSEP, customers are 1, 2, ..., n_customers.
        _edge_tail .-= 1
        _edge_head .-= 1

        # In Julia,   depot is 1
        # In CVRPSEP, depot is n_customers+1
        _edge_tail[_edge_tail .== 0] .= n_customers + 1
        _edge_head[_edge_head .== 0] .= n_customers + 1

        # In CVRPSEP, edges with positive x values are numbered 1, 2, ..., n_edges
        # a dummy zero is necessary to fill the first element
        _edge_tail = [-1; _edge_tail]
        _edge_head = [-1; _edge_head]
        _edge_x = [-1; _edge_x]

        @assert length(_edge_tail) == n_edges + 1

        return n_customers, _demand, n_edges, _edge_tail, _edge_head, _edge_x
end


text:
module TestFold

include("../preamble.jl")

@testset "$fold" for fold in [foldl, reduce, dreduce]
    @testset "IdentityTransducer" begin
        @test fold(+, IdentityTransducer(), 1:10) == sum(1:10)
    end
end

end  # module


text:
using Core.Compiler: LineInfoNode
import Base: push!, insert!, getindex, setindex!, iterate, length

# We have our own versions of these in order to
# (1) be more robust to Base IR changes, and
# (2) make sure that mistakes/bugs do not cause bad LLVM IR.

struct Undefined end
const undef = Undefined()

"""
    Variable(N)
    var(N)

Represents an SSA variable. Primarily used as an index into `IR` objects.
"""
struct Variable
  id::Int
end

var(id::Integer) = Variable(id)

struct Branch
  condition::Any
  block::Int
  args::Vector{Any}
end

Branch(br::Branch; condition = br.condition,
                   block = br.block, args = br.args) =
  Branch(condition, block, args)

isreturn(b::Branch) = b.block == 0 && length(b.args) == 1
isconditional(b::Branch) = b.condition != nothing

Base.:(==)(a::Branch, b::Branch) =
  (a.condition, a.block, a.args) == (b.condition, b.block, b.args)

arguments(b::Branch) = b.args

const unreachable = Branch(nothing, 0, [])

"""
    Statement(expr; type, line)

Represents a single statement in the IR. The `expr` is a non-nested Julia
expression (`Expr`). `type` represents the return type of the statement; in most
cases this can be ignored and defaults to `Any`. `line` represents the source
location of the statement; it is an integer index into the IR's line table.

As a convenience, if `expr` is already a statement, the new statement will
inherit its type and line number.
"""
struct Statement
  expr::Any
  type::Any
  line::Int
end

Statement(x; expr = x, type = Any, line = 0) =
  Statement(expr, type, line)

Statement(x::Statement; expr = x.expr, type = x.type, line = x.line) =
  Statement(expr, type, line)

const stmt = Statement

struct BasicBlock
  stmts::Vector{Statement}
  args::Vector{Any}
  argtypes::Vector{Any}
  branches::Vector{Branch}
end

BasicBlock(stmts = []) = BasicBlock(stmts, [], [], Branch[])

Base.copy(bb::BasicBlock) = BasicBlock(copy(bb.stmts), copy(bb.args), copy(bb.argtypes), copy(bb.branches))

branches(bb::BasicBlock) = bb.branches
arguments(bb::BasicBlock) = bb.args

"""
    IR()
    IR(metadata; slots = false)

Represents a fragment of SSA-form code.

IR can be constructed from scratch, but more usually an existing Julia method is
used as a starting point (see [`meta`](@ref) for how to get metadata for a
method). The `slots` argument determines whether the IR preserves mutable
variable slots; by default, these are converted to SSA-style variables.

As a shortcut, IR can be constructed directly from a type signature, e.g.

    julia> IR(typeof(gcd), Int, Int)
    1: (%1, %2, %3)
      %4 = %2 == 0
      br 4 unless %4
    2: ...
"""
struct IR
  defs::Vector{Tuple{Int,Int}}
  blocks::Vector{BasicBlock}
  lines::Vector{LineInfoNode}
  meta::Any
end

IR(; meta = nothing) = IR([],[BasicBlock()],[],meta)
IR(lines::Vector{LineInfoNode}; meta = nothing) = IR([],[BasicBlock()],lines,meta)

Base.copy(ir::IR) = IR(copy(ir.defs), copy.(ir.blocks), copy(ir.lines), ir.meta)

length(ir::IR) = sum(x -> x[2] > 0, ir.defs)

function block!(ir::IR, i = length(blocks(ir))+1)
  insert!(ir.blocks, i, BasicBlock())
  if i != length(blocks(ir))
    for b in blocks(ir), i = 1:length(branches(b))
      br = branches(b)[i]
      br.block >= i && (branches(b)[i] = Branch(br, block = br.block+1))
    end
  end
  for (ii, (b, j)) = enumerate(ir.defs)
    b >= i && (ir.defs[ii] = (b+1, j))
  end
  return block(ir, i)
end

struct Block
  ir::IR
  id::Int
end

basicblock(b::Block) = b.ir.blocks[b.id]
branches(b::Block) = branches(basicblock(b))
arguments(b::Block) = arguments(basicblock(b))
arguments(ir::IR) = arguments(block(ir, 1))

canbranch(bb::Block) = length(branches(bb)) == 0 || isconditional(branches(bb)[end])

isreturn(b::Block) = any(isreturn, branches(b))

function explicitbranch!(b::Block)
  b.id == 1 && return
  a = block(b.ir, b.id-1)
  if all(isconditional, branches(a))
    branch!(a, b.id)
  end
  return
end

explicitbranch!(ir::IR) = foreach(explicitbranch!, blocks(ir))

function branches(b::Block, c::Block)
  c.id == b.id+1 && explicitbranch!(c)
  filter(br -> br.block == c.id, branches(b))
end

branches(b::Block, c::Integer) = branches(b, block(b.ir, c))

"""
    returnvalue(block)

Retreive the return value of a block.

    julia> f(x) = 3x + 2;

    julia> IRTools.block(@code_ir(f(1)), 1)
    1: (%1, %2)
      %3 = 3 * %2
      %4 = %3 + 2
      return %4

    julia> IRTools.returnvalue(ans)
    %4
"""
function returnvalue(b::Block)
  isreturn(branches(b)[end]) || error("Block does not return")
  return branches(b)[end].args[1]
end

"""
    argument!(block, [value, type])

Create a new argument for the given block / IR fragment, and return the variable
representing the argument.

    julia> ir = IR();

    julia> argument!(ir)
    %1

    julia> ir
    1: (%1)

The `at` keyword argument can be used to specify where the new argument should
go; by default it is appended to the end of the argument list.

If there are branches to this block, they will be updated to pass `value`
(`nothing` by default) as an argument.
"""
function argument!(b::Block, value = nothing, type = Any; insert = true, at = length(arguments(b))+1)
  if at < length(arguments(b))
    for i = 1:length(b.ir.defs)
      (c, j) = b.ir.defs[i]
      c == b.id && -j >= at && (b.ir.defs[i] = (c, j-1))
    end
  end
  push!(b.ir.defs, (b.id, -at))
  arg = var(length(b.ir.defs))
  insert!(arguments(b), at, arg)
  insert!(basicblock(b).argtypes, at, type)
  if insert
    explicitbranch!(b)
    for c in blocks(b.ir), br in branches(c)
      br.block == b.id && insert!(arguments(br), at, value)
    end
  end
  return arg
end

argument!(ir::IR, a...; kw...) =
  argument!(block(ir, 1), nothing, a...; kw..., insert = false)

function emptyargs!(b::Block)
  empty!(arguments(b))
  for c in blocks(b.ir), br in branches(c)
    br.block == b.id && empty!(arguments(br))
  end
  return
end

function deletearg!(b::Block, i)
  arg = arguments(b)[i]
  deleteat!(arguments(b), i)
  for c in blocks(b.ir), br in branches(c)
    br.block == b.id && deleteat!(arguments(br), i)
  end
  b.ir.defs[arg.id] = (-1, -1)
  return
end

function deletearg!(b::Block, i::AbstractVector)
  for i in sort(i, lt = >)
    deletearg!(b, i)
  end
end

deletearg!(ir::IR, i) = deletearg!(block(ir, 1), i)

block(ir::IR, i) = Block(ir, i)
blocks(ir::IR) = [block(ir, i) for i = 1:length(ir.blocks)]

function blockidx(ir::IR, x::Variable)
  b, i = get(ir.defs, x.id, (-1, -1))
  i > 0 || error("No such variable $x")
  block(ir, b), i
end

"""
    haskey(ir, var)

Check whether the variable `var` was defined in `ir`.

    julia> f(x) = 3x + 2;

    julia> ir = @code_ir f(1)
    1: (%1, %2)
      %3 = 3 * %2
      %4 = %3 + 2
      return %4

    julia> haskey(ir, var(3))
    true

    julia> haskey(ir, var(7))
    false
"""
function Base.haskey(ir::IR, x::Variable)
  b, i = get(ir.defs, x.id, (-1, -1))
  return i > 0
end

getindex(b::Block, i::Integer) = basicblock(b).stmts[i]
getindex(b::Block, i::Variable) = b.ir[i]
setindex!(b::Block, x::Statement, i::Integer) = (basicblock(b).stmts[i] = x)
setindex!(b::Block, x, i::Integer) = (b[i] = Statement(b[i], expr = x))

branch(block::Integer, args...; unless = nothing) =
  Branch(unless, block, Any[args...])

branch(block::Block, args...; kw...) = branch(block.id, args...; kw...)

function branch!(b::Block, block, args...; unless = nothing)
  brs = branches(b)
  unless === nothing && deleteat!(brs, findall(br -> br.condition === nothing, brs))
  args = map(a -> push!(b, a), args)
  push!(brs, branch(block, args...; unless = unless))
  return b
end

function branch!(ir::IR, args...; kw...)
  branch!(blocks(ir)[end], args...; kw...)
  return ir
end

return!(ir, x) = branch!(ir, 0, x)

function getindex(ir::IR, i::Variable)
  b, i = blockidx(ir, i)
  return b[i]
end

Base.get(ir::IR, i::Variable, default) = haskey(ir, i) ? ir[i] : default

function setindex!(ir::IR, x, i::Variable)
  b, i = blockidx(ir, i)
  b[i] = x
end

function Base.delete!(ir::IR, i::Variable)
  ir[i] = nothing
  ir.defs[i.id] = (-1, -1)
  return ir
end

length(b::Block) = count(x -> x[1] == b.id, b.ir.defs)

function successors(b::Block)
  brs = basicblock(b).branches
  succs = Int[br.block for br in brs if br.block > 0]
  all(br -> br.condition != nothing, brs) && push!(succs, b.id+1)
  return [block(b.ir, succ) for succ in succs]
end

predecessors(b::Block) = [c for c in blocks(b.ir) if b in successors(c)]

"""
    keys(ir)

Return the variable keys for all statements defined in `ir`.

    julia> f(x) = 3x + 2;

    julia> ir = @code_ir f(1)
    1: (%1, %2)
      %3 = 3 * %2
      %4 = %3 + 2
      return %4

    julia> keys(ir)
    2-element Array{IRTools.Variable,1}:
     %3
     %4
"""
Base.keys(b::Block) = first.(sort([Variable(i) => v for (i, v) in enumerate(b.ir.defs) if v[1] == b.id && v[2] > 0], by = x->x[2]))

function iterate(b::Block, (ks, i) = (keys(b), 1))
  i > length(ks) && return
  return (ks[i]=>b.ir[ks[i]], (ks, i+1))
end

Base.keys(ir::IR) = first.(sort([Variable(i) => v for (i, v) in enumerate(ir.defs) if v[2] > 0], by = x->x[2]))

function iterate(ir::IR, (ks, i) = (keys(ir), 1))
  i > length(ks) && return
  return (ks[i]=>ir[ks[i]], (ks, i+1))
end

applyex(f, x) = x
applyex(f, x::Expr) =
  Expr(x.head, [x isa Expr ? f(x) : x for x in x.args]...)
applyex(f, x::Statement) = Statement(x, expr = applyex(f, x.expr))

"""
    push!(ir, x)

Append the statement or expression `x` to the IR or block `ir`, returning the
new variable. See also [`pushfirst!`](@ref), [`insert!`](@ref).

    julia> ir = IR();

    julia> x = argument!(ir)
    %1

    julia> push!(ir, xcall(:*, x, x))
    %2

    julia> ir
    1: (%1)
      %2 = %1 * %1
"""
function push!(b::Block, x::Statement)
  x = applyex(a -> push!(b, Statement(x, expr = a)), x)
  x = Statement(x)
  push!(basicblock(b).stmts, x)
  push!(b.ir.defs, (b.id, length(basicblock(b).stmts)))
  return Variable(length(b.ir.defs))
end

push!(b::Block, x) = push!(b, Statement(x))

push!(b::Block, x::Variable) = x

# TODO make this work on nested Exprs.
function insert!(b::Block, idx::Integer, x)
  insert!(basicblock(b).stmts, idx, Statement(x))
  for i = 1:length(b.ir.defs)
    c, j = b.ir.defs[i]
    if c == b.id && j >= idx
      b.ir.defs[i] = (c, j+1)
    end
  end
  push!(b.ir.defs, (b.id, idx))
  return Variable(length(b.ir.defs))
end

Base.pushfirst!(b::Block, x) = insert!(b, 1, x)

push!(ir::IR, x) = push!(block(ir, length(ir.blocks)), x)

"""
    pushfirst!(ir, x)

Insert the expression or statement `x` into the given IR or block at the
beginning, returning the new variable. See also [`push!`](@ref),
[`insert!`](@ref).

    julia> f(x) = 3x + 2
    f (generic function with 1 method)

    julia> ir = @code_ir f(1)
    1: (%1, %2)
      %3 = 3 * %2
      %4 = %3 + 2
      return %4

    julia> pushfirst!(ir, :(println("hello, world")))
    %5

    julia> ir
    1: (%1, %2)
      %5 = println("hello, world")
      %3 = 3 * %2
      %4 = %3 + 2
      return %4
"""
Base.pushfirst!(ir::IR, x) = pushfirst!(block(ir, 1), x)

"""
    insert!(ir, v, x)

Insert the expression or statement `x` into the given IR, just before the
variable `v` is defined, returning the new variable for `x`. See also
[`insertafter!`](@ref).

    julia> f(x) = 3x + 2
    f (generic function with 1 method)

    julia> ir = @code_ir f(1)
    1: (%1, %2)
      %3 = 3 * %2
      %4 = %3 + 2
      return %4

    julia> insert!(ir, IRTools.var(4), :(println("hello, world")))
    %5

    julia> ir
    1: (%1, %2)
      %3 = 3 * %2
      %5 = println("hello, world")
      %4 = %3 + 2
      return %4
"""
function insert!(ir::IR, i::Variable, x; after = false)
  if after && ir.defs[i.id][2] < 0
    pushfirst!(block(ir, ir.defs[i.id][1]), x)
  else
    b, i = blockidx(ir, i)
    insert!(b, i+after, x)
  end
end

"""
    insertafter!(ir, v, x)

Insert the expression or statement `x` into the given IR, just before the
variable `v` is defined, returning the new variable for `x`. See also
[`insert!`](@ref).

    julia> f(x) = 3x + 2
    f (generic function with 1 method)

    julia> ir = @code_ir f(1)
    1: (%1, %2)
      %3 = 3 * %2
      %4 = %3 + 2
      return %4

    julia> IRTools.insertafter!(ir, IRTools.var(4), :(println("hello, world")))
    %5

    julia> ir
    1: (%1, %2)
      %3 = 3 * %2
      %4 = %3 + 2
      %5 = println("hello, world")
      return %4
"""
insertafter!(ir, i, x) = insert!(ir, i, x, after=true)

"""
    empty(ir)

Create an empty IR fragment based on the given IR. The line number table and
any metadata are preserved from the original IR.

    julia> ir = empty(@code_ir gcd(10, 5))
    1:

    julia> ir.meta
    Metadata for gcd(a::T, b::T) where T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8} in Base at intfuncs.jl:31
"""
Base.empty(ir::IR) = IR(copy(ir.lines), meta = ir.meta)

function Base.permute!(ir::IR, perm::AbstractVector)
  explicitbranch!(ir)
  permute!(ir.blocks, perm)
  iperm = invperm(perm)
  for v = 1:length(ir.defs)
    b, i = ir.defs[v]
    b == -1 && continue
    ir.defs[v] = (iperm[b], i)
  end
  for b in blocks(ir), i = 1:length(branches(b))
    branches(b)[i].block > 0 || continue
    br = branches(b)[i]
    branches(b)[i] = Branch(br, block = iperm[br.block])
  end
  return ir
end

# Pipe

struct NewVariable
  id::Int
end

"""
    Pipe(ir)

In general, it is not efficient to insert statements into IR; only appending is
fast, for the same reason as with `Vector`s.

For this reason, the `Pipe` construct makes it convenient to incrementally build
an new IR fragment from an old one, making efficient modifications as you go.

The general pattern looks like:

    pr = IRTools.Pipe(ir)
    for (v, st) in pr
      # do stuff
    end
    ir = IRTools.finish(pr)

Iterating over `pr` is just like iterating over `ir`, except that within the
loop, inserting and deleting statements in `pr` around `v` is efficient. Later,
`finish(pr)` converts it back to a normal IR fragment (in this case just a plain
copy of the original).
"""
mutable struct Pipe
  from::IR
  to::IR
  map::Dict{Any,Any}
  var::Int
end

var!(p::Pipe) = NewVariable(p.var += 1)

substitute!(p::Pipe, x, y) = (p.map[x] = y; x)
substitute(p::Pipe, x::Union{Variable,NewVariable}) = p.map[x]
substitute(p::Pipe, x) = get(p.map, x, x)
substitute(p::Pipe) = x -> substitute(p, x)

function Pipe(ir)
  p = Pipe(ir, IR(copy(ir.lines), meta = ir.meta), Dict(), 0)
  for (x, T) in zip(p.from.blocks[1].args, p.from.blocks[1].argtypes)
    y = argument!(blocks(p.to)[end], nothing, T, insert = false)
    substitute!(p, x, y)
  end
  return p
end

function pipestate(ir::IR)
  ks = sort([Variable(i) => v for (i, v) in enumerate(ir.defs) if v[2] > 0], by = x->x[2])
  [first.(filter(x -> x[2][1] == b, ks)) for b = 1:length(ir.blocks)]
end

function iterate(p::Pipe, (ks, b, i) = (pipestate(p.from), 1, 1))
  if i == 1 && b != 1
    for (x, T) in zip(p.from.blocks[b].args, p.from.blocks[b].argtypes)
      y = argument!(blocks(p.to)[end], nothing, T, insert = false)
      substitute!(p, x, y)
    end
  end
  if i > length(ks[b])
    for br in branches(block(p.from, b))
      push!(p.to.blocks[end].branches, map(substitute(p), br))
    end
    b == length(ks) && return
    block!(p.to)
    return iterate(p, (ks, b+1, 1))
  end
  v = ks[b][i]
  st = p.from[v]
  substitute!(p, v, push!(p.to, prewalk(substitute(p), st)))
  ((v, st), (ks, b, i+1))
end

finish(p::Pipe) = p.to

islastdef(ir::IR, v::Variable) =
  v.id == length(ir.defs) &&
  ir.defs[v.id] == (length(ir.blocks), length(ir.blocks[end].stmts))

setindex!(p::Pipe, x, v) = p.to[substitute(p, v)] = prewalk(substitute(p), x)

function Base.push!(p::Pipe, x)
  tmp = var!(p)
  substitute!(p, tmp, push!(p.to, prewalk(substitute(p), x)))
  return tmp
end

function Base.pushfirst!(p::Pipe, x)
  tmp = var!(p)
  substitute!(p, tmp, pushfirst!(p.to, prewalk(substitute(p), x)))
  return tmp
end

function Base.delete!(p::Pipe, v)
  v′ = substitute(p, v)
  delete!(p.map, v)
  if islastdef(p.to, v′)
    pop!(p.to.defs)
    pop!(p.to.blocks[end].stmts)
  else
    delete!(p.to, v′)
  end
end

function insert!(p::Pipe, v, x; after = false)
  v′ = substitute(p, v)
  x = prewalk(substitute(p), x)
  tmp = var!(p)
  if islastdef(p.to, v′) # we can make this case efficient by renumbering
    if after
      substitute!(p, tmp, push!(p.to, x))
    else
      substitute!(p, v, push!(p.to, p.to[v′]))
      p.to[v′] = Statement(x)
      substitute!(p, tmp, v′)
    end
  else
    substitute!(p, tmp, insert!(p.to, v′, x, after = after))
  end
  return tmp
end

argument!(p::Pipe, a...; kw...) =
  substitute!(p, var!(p), argument!(p.to, a...; kw...))

function branch!(ir::Pipe, b, args...; kw...)
  args = map(a -> postwalk(substitute(ir), a), args)
  branch!(blocks(ir.to)[end], b, args...; kw...)
  return ir
end


text:
include(joinpath(@__DIR__, "pedagogy_spatial.jl"))
using Random

β₀ = 0.1
r = 0.05  # interaction range = 2*radius
b = β₀ / (π * (r/2)^2)

c, w = [0.0, 0.0], 1.0
win = PRS.SquareWindow(c, w)

hc = PRS.HardCorePointProcess(b, r, win)

rng = Random.MersenneTwister(123)
@time sample = PRS.generate_sample_prs(rng, hc)
# @time sample = PRS.generate_sample_grid_prs(rng, hc)

p = plot(hc, sample; title="")


text:
using SumOfSquares
using MultivariateMoments
using MultivariatePolynomials
using PolyJuMP

"""
Finds the initial order of the hierarchy; i.e what degrees are required
to compute: Eμ[p]; M[gμ] ∀g ∈ Sy; and positivstellensatz of Eμ[p] on Sx.
"""
function min_order(
        p::AbstractPolynomial,
        Sx::AbstractSemialgebraicSet,
        Sy::AbstractSemialgebraicSet,
)
        round_up_even(i::Integer) = (i + 1) & ~(1)
        deg(s::AbstractBasicSemialgebraicSet) = maxdegree.(inequalities(s))
        deg(s::AbstractAlgebraicSet) = maxdegree.(equalities(s))
        deg(p::AbstractPolynomial) = maxdegree(p)

        y = variables(Sy)
        x = variables(Sx)
        px = subs(p, y => ones(length(y)))
        py = subs(p, x => ones(length(x)))

        round_up_even(maximum([deg(px); deg(py); deg(Sy); deg(Sx)]))
end

"""
Tries to solve a zero-sum two-player polynomial game on semialgebraic sets,
using an iteration of the moment--sos hierarchy.

The inputs are: a payoff function p, semialgebraic strategy sets Sx and Sy, and
an "iteration", where "minimal order + 2*iteration" defines the actual order.

Returns: the payoff and the optimal strategies of both players.

Lasserre, Jean-Bernard. (2004). Global Optimization With Polynomials And The
Problem Of Moments. SIAM Journal on Optimization. 11. 10.1137/S1052623400366802.
P. A. Parrilo (2006). Polynomial games and sum of squares optimization In
Proceedings of the 45th IEEE Conference on Decision and Control (pp. 2855-2860).
"""
function solve_game(
        p::AbstractPolynomial,
        Sx::AbstractSemialgebraicSet,
        Sy::AbstractSemialgebraicSet,
        optimizer_factory;
        iteration::Unsigned = 0x0,
)
        order = min_order(p, Sx, Sy) + iteration * 2

        m = SOSModel(optimizer_factory)

        @variable m α
        @objective m Min α
        @variable m μ Moments(order, Sy)
        @constraint(m, q, linearize(μ, p) <= α, domain = Sx, maxdegree = order)

        set_silent(m)
        optimize!(m)

        if !available_result(m)
                return JuMP.termination_status(m), nothing, nothing
        end

        xmat = moment_matrix(q)
        ymat = moment_matrix(μ)
        value(α), extractatoms(xmat, 1e-4), extractatoms(ymat, 1e-4)
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# Defines the NonlinearExpression, BufferedNonlinearFunction used in
# constructing relaxations of nonlinear functions along with a number of
# helper functions including an Evaluator structure and: set_node_flag!
# set_node!, set_reference_point!, retrieve_node, prior_eval
# copy_subexpression_value!, eliminate_fixed_variables!
#############################################################################

const DEBUG_NL = false

include(joinpath(@__DIR__, "register_special.jl"))
include(joinpath(@__DIR__, "graph", "graph.jl"))
include(joinpath(@__DIR__, "interval", "interval.jl"))
include(joinpath(@__DIR__, "composite_relax", "composite_relax.jl"))
include(joinpath(@__DIR__, "apriori_relax", "apriori_relax.jl"))

@enum(RelaxType, STD_RELAX, MC_AFF_RELAX, MC_ENUM_RELAX)

_set_has_value!(d, v) = nothing

"""
$(TYPEDEF)

Stores a general quadratic function with a buffer.
"""
mutable struct NonlinearExpression{V,N,T<:RelaxTag} <: AbstractEAGOConstraint
    g::DirectedTree
    relax_cache::RelaxCache{V,N,T}
    has_value::Bool
    last_reverse::Bool
    lower_bound::Float64
    upper_bound::Float64
    grad_sparsity::Vector{Int}
end
function NonlinearExpression()
    g = DirectedTree()
    c = RelaxCache{MC{1,NS},1,NS}()
    return NonlinearExpression{MC{1,NS},1,NS}(g, c, false, false, -Inf, Inf, Int[])
end

relax_info(s::Relax, n::Int, t::T) where T = MC{n,T}
function NonlinearExpression!(aux_info, rtype::S, sub::Union{JuMP._SubexpressionStorage,JuMP._FunctionStorage},
                              b::MOI.NLPBoundsPair, sub_sparsity::Dict{Int,Vector{Int}},
                              subexpr_indx::Int,
                              subexpr_linearity::Vector{JuMP._Derivatives.Linearity},
                              op::OperatorRegistry, parameter_values,
                              tag::T, use_apriori_flag::Bool; is_sub::Bool = false) where {S,T}
    g = DirectedTree(aux_info, sub, op, sub_sparsity, subexpr_linearity, parameter_values, is_sub, subexpr_indx)
    grad_sparsity = sparsity(g, 1)
    n = length(grad_sparsity)
    V = relax_info(rtype, n, tag)
    c = RelaxCache{V,n,T}()
    c.use_apriori_mul = use_apriori_flag
    initialize!(c, g)
    return NonlinearExpression{V,n,T}(g, c, false, false, b.lower, b.upper, grad_sparsity)
end

@inline has_value(d::NonlinearExpression) = d.has_value
@inline dep_subexpr_count(d::NonlinearExpression) = dep_subexpr_count(d.g)
@inline set_has_value!(d::NonlinearExpression, v::Bool) = (d.has_value = v; return )
function _set_last_reverse!(d::NonlinearExpression{V, N, T}, v::Bool) where {V,N,T<:RelaxTag}
    d.last_reverse = v;
    return
end
function set_variable_storage!(d::NonlinearExpression, v::VariableValues{S}) where S<:Real
    d.relax_cache.ic.v = v
    return
end
@inbounds sparsity(d::NonlinearExpression) = sparsity(d.g, 1)
@inbounds set(d::NonlinearExpression{V,N,T}) where {V,N,T<:RelaxTag} = set(d.relax_cache, 1)
@inbounds info(d::NonlinearExpression{V,N,T}) where {V,N,T<:RelaxTag} = info(d.relax_cache, 1)
@inbounds num(d::NonlinearExpression{V,N,T}) where {V,N,T<:RelaxTag} = num(d.relax_cache, 1)
@inbounds is_num(d::NonlinearExpression) = is_num(d.relax_cache, 1)
var_num(d::NonlinearExpression{V,N,T}) where {V,N,T<:RelaxTag} = N

mc_type(rc::NonlinearExpression{V,N,T}) where {V,N,T<:RelaxTag} = MC{N,T}

"""
$(TYPEDEF)

Stores a general nonlinear function with a buffer represented by the sum of a tape
and a scalar affine function.
"""
mutable struct BufferedNonlinearFunction{V,N,T<:RelaxTag} <: AbstractEAGOConstraint
    ex::NonlinearExpression{V,N,T}
    saf::SAF
end
function BufferedNonlinearFunction()
    ex = NonlinearExpression()
    saf = SAF(SAT[], 0.0)
    return BufferedNonlinearFunction{MC{1,NS},1,NS}(ex, saf)
end

function BufferedNonlinearFunction(aux_info, rtype::RELAX_ATTRIBUTE, f::JuMP._FunctionStorage, b::MOI.NLPBoundsPair,
                                   sub_sparsity::Dict{Int,Vector{Int}},
                                   subexpr_lin::Vector{JuMP._Derivatives.Linearity},
                                   op::OperatorRegistry, parameter_values,
                                   tag::T, use_apriori_flag::Bool) where T <: RelaxTag

    ex = NonlinearExpression!(aux_info, rtype, f, b, sub_sparsity, -1, subexpr_lin, op, parameter_values, tag, use_apriori_flag)
    n = length(sparsity(ex.g, 1))
    saf = SAF(SAT[SAT(0.0, VI(i)) for i = 1:n], 0.0)
    V = relax_info(rtype, n, tag)
    return BufferedNonlinearFunction{V,n,T}(ex, saf)
end

function expand_sv!(out::Vector{Float64}, n::Int, m::Int, vs::Vector{Int}, gs::Vector{Int}, x::SVector{N,T}) where {N,T}
    k = 1
    for q = 1:m
        i = @inbounds vs[q]
        for j = k:n
            if i == @inbounds gs[j]
                @inbounds out[j] = x[q]
                k = j
                break
            end
        end
    end
    nothing
end
function _load_subexprs!(d::RelaxCache{V,N,T}, g, subexpressions, dep_subexprs) where {V,N,T<:RelaxTag}
    gs = sparsity(g, 1)
    for (i,ds) in enumerate(dep_subexprs)
        s = subexpressions[ds]
        if is_num(s)
            store_subexpression_num!(d, num(s), i)
        else
            vs = sparsity(s)
            v = set(s)
            m = var_num(s)
            expand_sv!(d._cv_grad_buffer, N, m, vs, gs, v.cv_grad)
            expand_sv!(d._cc_grad_buffer, N, m, vs, gs, v.cc_grad)
            cvg = SVector{N,Float64}(d._cv_grad_buffer)
            ccg = SVector{N,Float64}(d._cc_grad_buffer)
            store_subexpression_set!(d, MC{N,T}(v.cv, v.cc, v.Intv, cvg, ccg, false), i)
        end
    end
    return
end

@inline _set_last_reverse!(d::BufferedNonlinearFunction{V,N,T}, v::Bool) where {V,N,T<:RelaxTag} = _set_last_reverse!(d.ex, v)
function set_variable_storage!(d::BufferedNonlinearFunction{V,N,T}, v::VariableValues{Float64}) where {V,N,T<:RelaxTag}
    set_variable_storage!(d.ex, v)
end

has_value(d::BufferedNonlinearFunction) = has_value(d.ex)
dep_subexpr_count(d::BufferedNonlinearFunction) = dep_subexpr_count(d.ex)
set_has_value!(d::BufferedNonlinearFunction, v::Bool) = set_has_value!(d.ex, v)
sparsity(d::BufferedNonlinearFunction) = sparsity(d.ex)
set(d::BufferedNonlinearFunction{V,N,T}) where {V,N,T<:RelaxTag} = set(d.ex)
num(d::BufferedNonlinearFunction{V,N,T}) where {V,N,T<:RelaxTag} = num(d.ex)
lower_bound(d::BufferedNonlinearFunction{V,N,T}) where {V,N,T<:RelaxTag} = d.ex.lower_bound
upper_bound(d::BufferedNonlinearFunction{V,N,T}) where {V,N,T<:RelaxTag} = d.ex.upper_bound

# returns the interval bounds associated with the set
interval(d::BufferedNonlinearFunction{V,N,T}) where {V,N,T<:RelaxTag} = Interval{Float64}(set(d))
is_num(d::BufferedNonlinearFunction) = is_num(d.ex)

mc_type(rc::BufferedNonlinearFunction{V,N,T}) where {V,N,T<:RelaxTag} = MC{N,T}

"""
    Evaluator

MOI.AbstractNLPEvaluator for calculating relaxations of nonlinear terms.

Checks that the resulting value should be a number...

$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct Evaluator <: MOI.AbstractNLPEvaluator
    user_operators::OperatorRegistry = OperatorRegistry()
    has_user_mv_operator::Bool = false
    num_mv_buffer::Vector{Float64} = Float64[]
    parameter_values::Vector{Float64} = Float64[]
    node::NodeBB = NodeBB()
    variable_values::VariableValues{Float64} = VariableValues{Float64}()
    subgrad_tighten::Bool = false
    reverse_subgrad_tighten::Bool = false
    subexpressions::Vector{NonlinearExpression} = NonlinearExpression[]
    subexpressions_eval::Vector{Bool}           = Bool[]
    is_post::Bool = false
    is_intersect::Bool = false
    is_first_eval::Bool = false
    interval_intersect::Bool = false
    subgrad_tol::Float64 = 1E-10
    relax_type::RelaxType                       = STD_RELAX
    pass_number::Int       = 0
end
set_variable_values!(d::Evaluator, v) = d.variable_values = v

"""
$(FUNCTIONNAME)

Sets the current node in the Evaluator structure.
"""
function set_node!(d::Evaluator, n::NodeBB)
    d.node = NodeBB(n)
    for i = 1:length(n)
        vi = d.variable_values.node_to_variable_map[i]
        d.variable_values.lower_variable_bounds[vi] = n.lower_variable_bounds[i]
        d.variable_values.upper_variable_bounds[vi] = n.upper_variable_bounds[i]
    end
    fill!(d.subexpressions_eval, false)
    d.is_first_eval = true
    return nothing
end

function retrieve_node(d::Evaluator)
    n = d.node
    nv_map = d.variable_values.node_to_variable_map
    return NodeBB(copy(d.variable_values.lower_variable_bounds[nv_map]),
                  copy(d.variable_values.upper_variable_bounds[nv_map]),
                  copy(n.is_integer), n.continuous,
                  n.lower_bound, n.upper_bound, n.depth, n.cont_depth, n.id,
                  n.branch_direction, n.last_branch, n.branch_extent)
end
@inline function _get_x!(::Type{BranchVar}, out::Vector{Float64}, d::Evaluator)
    return _get_x!(BranchVar, out, d.variable_values)
end
prior_eval(d::Evaluator, i::Int) = d.subexpressions_eval[i]

#=
Assumes the sparsities are sorted...
=#
function copy_subexpression_value!(k::Int, op::Int, subexpression::NonlinearExpression{V,MC{N1,T}},
                                   numvalued::Vector{Bool}, numberstorage::Vector{S}, setstorage::Vector{MC{N2,T}},
                                   cv_buffer::Vector{S}, cc_buffer::Vector{S},
                                   func_sparsity::Vector{Int}) where {V, N1, N2, S, T <: RelaxTag}

    # fill cv_grad/cc_grad buffers
    sub_sparsity = subexpression.grad_sparsity
    sset = subexpression.setstorage[1]
    fill!(cv_buffer, zero(S))
    fill!(cc_buffer, zero(S))

    sub_sparsity_count = 1
    subs_index = @inbounds sub_sparsity[1]
    for i = 1:length(func_sparsity)
        func_index = @inbounds func_sparsity[i]
        if func_index === subs_index
            @inbounds cv_buffer[i] = sset.cv_grad[sub_sparsity_count]
            @inbounds cc_buffer[i] = sset.cc_grad[sub_sparsity_count]
            sub_sparsity_count += 1
            subs_index = @inbounds sub_sparsity[sub_sparsity_count]
        end
    end

    cv_grad = SVector(cv_buffer)
    cc_grad = SVector(cc_buffer)

    setstorage[k] = MC{N1,T}(sset.cv, sset.cc, sset.Intv, cv_grad, cc_grad, sset.cnst)

    return nothing
end

function eliminate_fixed_variables!(f::NonlinearExpression{V,N,T}, v::Vector{VariableInfo}) where {V,N,T<:RelaxTag}
    num_constants = length(f.const_values)
    indx_to_const_loc = Dict{Int,Int}()
    for i = 1:length(expr.nd)
        nd = @inbounds expr.nd[i]
        if nd.nodetype === JuMP._Derivatives.VARIABLE
            indx = nd.index
            if v[indx].is_fixed
                if haskey(indx_to_const_loc, indx)
                    const_loc = indx_to_const_loc[indx]
                    expr.nd[i] = NodeData(JuMP._Derivatives.VALUE, const_loc, nd.parent)
                else
                    push!(const_values, v[indx].lower_bound)
                    num_constants += 1
                    indx_to_const_loc[indx] = num_constants
                    expr.nd[i] = NodeData(nd.nodetype, num_constants, nd.parent)
                end
                f.isnumber[i] = true
            end
        end
    end
    return nothing
end

eliminate_fixed_variables!(f::BufferedNonlinearFunction{N,T}, v::Vector{VariableInfo}) where {N,T<:RelaxTag} = eliminate_fixed_variables!(f.ex, v)
f_init_prop!(t, g::DAT, c::RelaxCache, flag::Bool) = flag ? f_init!(t, g, c) : fprop!(t, g, c)
function forward_pass!(z::Evaluator, d::NonlinearExpression{V,N,T}) where {V,N,T<:RelaxTag}
    #println("ran forward pass")
    b = d.relax_cache
    update_box_and_pnt!(b.ic.v, z.variable_values, z.is_first_eval)
    if b.use_apriori_mul
        s = sparsity(d)
        v = b.ic.v
        x = v.x
        for j in s
            if isone(z.pass_number)
                v.x0[j] = x[j]
            end
        end
        x0 = v.x0
        isempty(b.dp) && (b.dp = zeros(length(x));)
        isempty(b.dP) && (b.dP = zeros(Interval{Float64}, length(x));)
        isempty(b.p_rel) && (b.p_rel = zeros(length(x));)
        isempty(b.p_diam) && (b.p_diam  = zeros(length(x));)
        for j in s
            l = lbd(b, j)
            u = ubd(b, j)
            b.dp[j] = x[j] - x0[j]
            b.dP[j] = Interval(l, u) - x0[j]
            b.p_rel[j] = (x[j] - 0.5*(u + l))/(0.5*(u - l))
            b.p_diam[j] = u - l
        end
    end
    for i = 1:dep_subexpr_count(d)
        j = d.g.dependent_subexpressions[i]
        forward_pass!(z, z.subexpressions[j])
    end
    _load_subexprs!(d.relax_cache, d.g, z.subexpressions, d.g.dependent_subexpressions)
    (z.relax_type == STD_RELAX) && (return f_init_prop!(Relax(), d.g, d.relax_cache, z.is_first_eval))
    (z.relax_type == MC_AFF_RELAX) && (return f_init_prop!(RelaxAA(d.grad_sparsity), d.g, d.relax_cache, z.is_first_eval))
    return f_init_prop!(RelaxMulEnum(d.grad_sparsity), d.g, d.relax_cache, z.is_first_eval)
end

function forward_pass!(x::Evaluator, d::BufferedNonlinearFunction{V,N,T}) where {V,N,T<:RelaxTag}
    forward_pass!(x, d.ex)
    _set_has_value!(d, true)
    _set_last_reverse!(d, false)
    return
end

rprop!(::Relax, x::Evaluator, d::NonlinearExpression{V,N,T}) where {V,N,T<:RelaxTag} = rprop!(Relax(), d.g, d.relax_cache)
function rprop!(::Relax, x::Evaluator, d::BufferedNonlinearFunction{V,N,T}) where {V,N,T<:RelaxTag}
    _set_last_reverse!(d, true)
    return rprop!(Relax(), x, d.ex)
end

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
"""
    env_operator(env, channels)

Saves/Loads/Updates the search environment `env`.
Communication with the search server i.e. getting the command
and its arguments and sending back a new environment is done
via `channels`.
"""
function env_operator(env, channels)
    in_channel, out_channel = channels
    _env = env
    while true
        sleep(DEFAULT_ENVOP_SLEEP_INTERVAL)
        opdict = try
            JSON.parse(take!(in_channel))
        catch
            Dict{String, String}()
        end
        cmd = Symbol(get(opdict, "cmd", ""))
        cmd_argument = get(opdict, "cmd_argument", "")
        if cmd === :save
            try
                serialize(cmd_argument, _env)
                @info "• Environment successfully saved (serialized) in $cmd_argument."
            catch
                @warn "Could not save (serialize) environment in $cmd_argument."
            end
        elseif cmd === :load
            try
                _env = deserialize(cmd_argument)
                @info "• Environment successfully loaded (deserialized) from $cmd_argument."
            catch
                @warn "Could not load (deserialize) environment from $cmd_argument."
            end
        elseif cmd === :reindex
            try
                env_config = parse_configuration(env.config_path)
                _dbdata = env_config.data_loader()
                db_check_id_key(_dbdata, env_config.id_key)

                # Selectively reload
                new_searchers = similar(env.searchers)  # initialize
                cnt = 0
                for i in eachindex(env.searchers)
                    if cmd_argument == "*" || isequal(id(env.searchers[i]), cmd_argument)
                        new_searchers[i] = build_searcher(_dbdata,
                                                          env.embedders,
                                                          env_config.searcher_configs[i];
                                                          id_key=env_config.id_key)
                        cnt+= 1
                    else
                        new_searchers[i] = env.searchers[i]
                    end
                end
                _env = SearchEnv(_dbdata,
                                 env_config.id_key,
                                 env.sampler,
                                 env.embedders,
                                 new_searchers,
                                 env_config.config_path)
                @info "• Updated $cnt searcher(s) in the environment."
            catch
                @warn "Could not update environment (searchers=$cmd_argument)."
            end
        else
            @warn "Unknown environment operation command $cmd"
        end
        put!(out_channel, _env)
    end
    return nothing
end


text:
@testset "resample.jl" begin
    println("             > resample.jl")
    @testset "benchmark" begin
        m_out = 1024
        m_in  = 1024
        state = BinomialState(128, m_out, m_in)
        model = BinomialGridModel(
            m_out,
            1:5,
            LinRange(0.05,0.95,5),
            LinRange(0.1,2,5),
            LinRange(0.05,2,5),
            LinRange(0.05,2,5)
        )

        u = CUDA.rand(m_out)

        if RUN_BENCHMARKS
            println("")
            println("Benchmarking function outer_resample!: should take about 300μs")
            display(@benchmark outer_resample!($state, $model, $u))
            println("")
        end
    end

    @testset "correctness of values" begin
        using BinomialSynapses: outer_indices!

        for M_out in [4, 8, 16, 32, 64, 128, 256, 512, 1024]
            u = rand(M_out)
            v = cumsum(u)
            outer_indices!(u)
            @test u[2:M_out] ≈ v[1:M_out-1]
        end

        count = 0
        for i in 1:100
            u = cu([1f3, 1f0, 1f0, 1f3])
            idx = Array(outer_indices!(u))
            for j in 1:4
                if idx[j] == 2 || idx[j] == 3
                    count += 1
                end
            end
        end
        @test count / 400 < 0.05
    end

    @testset "sortedness of idx" begin
        using BinomialSynapses: outer_indices!

        for M_out in [4, 8, 16, 32, 64, 128, 256, 512, 1024]
            u = rand(M_out)
            idx = outer_indices!(u)
            @test issorted(idx)
        end
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Compat.Test

using Distributions, PDMats, StatsBase

struct test_density <: AbstractDensity
end

BAT.nparams(td::test_density) = Int(3)
BAT.sampler(td::test_density) = BAT.sampler(MvNormal(ones(3), PDMat(eye(3))))

@testset "density" begin

    mvec = [-0.3, 0.3]
    cmat = [1.0 1.5; 1.5 4.0]
    Σ = @inferred PDMat(cmat)
    mvnorm = @inferred  MvNormal(mvec, Σ)

    density = @inferred GenericDensity(params -> logpdf(mvnorm, params), 2)

    econtext = @inferred ExecContext()

    params = [0.0 -0.3 0.5; 0.0 0.3 0.6]

    @testset "rand" begin
        td = test_density()
        @test rand(MersenneTwister(7002), td, Float64) ≈
            [-2.415270938, 0.7070171342, 1.0224848653]
        @test rand(MersenneTwister(7002), td, Float64, 2) ≈
            [-2.415270938 1.2951273090;
            0.7070171342  0.8896838714;
            1.0224848653  0.8824274590]
        x = ones(3,2)
        rand!(MersenneTwister(7002), td, x)
        @test x ≈
            [-2.415270938 1.2951273090;
            0.7070171342  0.8896838714;
            1.0224848653  0.8824274590]
    end

    @testset "param_bounds" begin
        pbounds = @inferred param_bounds(density)
        @test typeof(pbounds) == NoParamBounds
        @test pbounds.ndims == 2
    end

    @testset "convert" begin
        cdensity = @inferred convert(
            AbstractDensity, HyperRectBounds([-1., 0.5], [2.,1], BAT.hard_bounds))
        @test typeof(cdensity) <: ConstDensity
    end

    @testset "density_logval!" begin
        r = zeros(Float64, size(params, 2))

        density_logval!(r, density, params, econtext)
        @test r ≈ logpdf(mvnorm, params)
        BAT.unsafe_density_logval!(r, density, params, econtext)
        @test r ≈ logpdf(mvnorm, params)
    end

    @testset "density_logval" begin
        @test density_logval(density, params[:, 1], econtext) ≈ logpdf(mvnorm, params[:, 1])
        @test BAT.unsafe_density_logval(density, params[:, 1], econtext) ≈ logpdf(mvnorm, params[:, 1])
    end

    @testset "exec_capabilities" begin
        ecap = @inferred BAT.exec_capabilities(
            density_logval!, similar(params[1, :]), density, params)
        @test ecap.nthreads == 1
        @test ecap.threadsafe == false
        @test ecap.nprocs == 1
        @test ecap.remotesafe == true

        ecap = @inferred BAT.exec_capabilities(
            BAT.unsafe_density_logval!, similar(params[1, :]), density, params)
        @test ecap.nthreads == 1
        @test ecap.threadsafe == false
        @test ecap.nprocs == 1
        @test ecap.remotesafe == true

        ecap = @inferred BAT.exec_capabilities(density_logval, density, params[:, 1])
        @test ecap.nthreads == 1
        @test ecap.threadsafe == true
        @test ecap.nprocs == 1
        @test ecap.remotesafe == true

        ecap = @inferred BAT.exec_capabilities(
            BAT.unsafe_density_logval, density, params[:, 1])
        @test ecap.nthreads == 1
        @test ecap.threadsafe == true
        @test ecap.nprocs == 1
        @test ecap.remotesafe == true

        ecap = @inferred BAT.exec_capabilities(
            BAT.unsafe_density_logval, test_density())
        @test ecap.nthreads == 1
        @test ecap.threadsafe == false
        @test ecap.nprocs == 1
        @test ecap.remotesafe == true

    end

    @testset "parent" begin
        @test @inferred parent(density) == density.log_f
    end

    @testset "nparams" begin
        @test @inferred nparams(density) == 2
    end
end


text:
using WordTokenizers

"""
The function in the origin gpt code

fixes some issues the spacy tokenizer had on books corpus
also does some whitespace standardization
"""
function text_standardize(text)
    text = lowercase(text)
    text = replace(text, r"([a-z])(,|\.)"=>s"\1 \2")
    text = replace(text, "—"=>"-")
    text = replace(text, "–"=>"-")
    text = replace(text, "―"=>"-")
    text = replace(text, "…"=>"...")
    text = replace(text, "´"=>"'")
    text = replace(text, r"""(-+|~+|!+|"+|;+|\?+|\++|,+|\)+|\(+|\\+|\/+|\*+|\[+|\]+|}+|{+|\|+|_+)"""=>s" \1 ")
    text = replace(text, r"\s*\n\s*"=>" \n ")
    text = replace(text, r"[^\S\n]+"=>" ")
    strip(text)
end

"""
    gpt_tokenizer(x)

An alternative for origin tokenizer (spacy tokenizer) used in gpt model.
"""
gpt_tokenizer(x) = toktok_tokenize(text_standardize(x))


text:
module test1
name=1
end

module test2
name=2
println( Main.name )
end
println( test1.name, " ",test2.name)

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
# Lattice generation

module MyLattice

export Lattice, coord, flavor, nns, nnns, num_sites, num_uc, adj
export bravais_full

# This will not be accessible outside the module
function bravais_int(Xvecs)
    # Xvecs should be a d by d matrix of integers giving winding vectors as columns
    # This function returns a unique set of integer coordinates of points defining a cluster given the integer winding numbers

    const nmax = 40
    dim = size(Xvecs,1)
    if dim == 3
    	newxvecs = Xvecs
    elseif dim == 2
    	newxvecs = zeros(Int64,3,3)
    	newxvecs[1:2,1:2] = Xvecs
    	newxvecs[3,3] = 1
    else
        error("Lattice must have dimension 2 or 3")
    end
    tp = dot(newxvecs[:,1],cross(newxvecs[:,2],newxvecs[:,3]))

    numpts = abs(tp)
    function reduced_pt(pt)
        # gives a unique integer point inside a supercell defined by Xvecs
        # Reciprocal lattice for the supercell

        dim = length(pt)
        if dim == 2
            pt3d = [pt; 0]
        else
            pt3d = pt
        end
        recip = Array{Float64}(3,3)
        for a in 1:3
            recip[:,a] = cross(newxvecs[:,mod(a,3)+1],newxvecs[:,mod(a+1,3)+1])/tp
        end
        fracpt = mod(round(recip'*float(pt3d),3),1)
        newpt = round(Int64,newxvecs*fracpt)[1:dim]
        return newpt
    end

    myset=Set{Vector{Int64}}()
    skip = false
    if dim == 3
        for i in 0:nmax
            for j in 0:nmax
                for k in 0:nmax
                    thispt = [i,j,k]
                    modpt = reduced_pt(thispt)
                    push!(myset,modpt)
                    length(myset)==numpts && (skip = true; break)
                end
                skip && break
            end
            skip && break
        end
    else
    	for i in 0:nmax
    		for j in 0:nmax
    			thispt = [i,j]
    			modpt = reduced_pt(thispt)
    			push!(myset,modpt)
    			length(myset)==numpts && (skip = true; break)
    		end
    		skip && break
    	end
    end

    intpts = Array{Int64}(dim,numpts)
    for (i,vec) in enumerate(myset)
        intpts[:,i] = vec
    end

    return intpts
end

# Lattice type.  Stores all sorts of good stuff

type Lattice
    sites::Matrix{Float64}
    flavor::Vector{Int64}
    firstneighbors::Matrix{Int64}
    secondneighbors::Matrix{Int64}
    adjacency::Matrix{Int64}
    nsites::Int64
    nuc::Int64
    numnns::Int64
    numnnns::Int64
end

coord(lat::Lattice,i) = lat.sites[:,i]
flavor(lat::Lattice,i) = lat.flavor[i]
nns(lat::Lattice,i) = lat.firstneighbors[:,i]
nnns(lat::Lattice,i) = lat.secondneighbors[:,i]
num_sites(lat::Lattice) = lat.nsites
num_uc(lat::Lattice) = lat.nuc
adj(lat::Lattice,i,j) = lat.adjacency[i,j]

function bravais_full(avecs,basis,Xvecs)
    # Will construct a Lattice object describing a cluster of lattice defined by these parameters
    # avecs should be a d by d matrix of primitive lattice vectors arranged in columns
    # basis should be a d by n matrix of basis vectors for n sites in a unit cell.
    # Xvecs should be a d by d matrix of integers giving winding vectors as columns
    # This function should return a full Lattice object

    dim = size(avecs,1)
    int_pts = bravais_int(Xvecs) # integer coordinates of distinct lattice points
    nuc = size(int_pts,2) # number of unit cells
    nbasis = size(basis,2) # number of sites in the primitive unit cell
    nsites = nuc*nbasis # total number of sites in the cluster

    # form real space winding vectors
    wind = avecs*Xvecs

    # form reciprocal space vectors for supercell
    if dim == 3
    	newwind = wind
    elseif dim == 2
    	newwind = zeros(3,3)
    	newwind[1:2,1:2] = wind
    	newwind[3,3] = 1.0
    else
        error("Lattice must have dimension 2 or 3")
    end
    tp = dot(newwind[:,1],cross(newwind[:,2],newwind[:,3]))
    bigrecip = Array{Float64}(3,3)
    for a in 1:3
        bigrecip[:,a] = cross(newwind[:,mod(a,3)+1],newwind[:,mod(a+1,3)+1])/tp
    end
    if dim == 3
        recip = bigrecip
    elseif dim == 2
        recip = bigrecip[1:2,1:2]
    end

    function reduced_pt(pt)
        # gives a unique point inside a supercell defined by wind
        fracpt = mod(recip'*pt,1)
        newpt = wind*fracpt
        return newpt
    end

    function distance(pt1,pt2)
        # finds closest distance between two points taking into account periodicity
        # round to 3 digits
        diff = pt1-pt2
        fracdiff = recip'*diff
        normfrac = mod(fracdiff+0.5*ones(dim),1)-0.5*ones(dim)
        newdiff = wind*normfrac
        return round(norm(newdiff),3)
    end

    # Compute all the sites in the unit cell
    bravais_points = avecs*int_pts
    sites = Array{Float64}(dim,nuc*nbasis) # real space coordinate of site i
    flavors = Array{Int64}(nuc*nbasis) # basis type of site i
    for i in 1:nuc
        for a in 1:nbasis
            sites[:,nuc*(a-1)+i] = reduced_pt(bravais_points[:,i]+basis[:,a])
            flavors[nuc*(a-1)+i] = a
        end
    end

    # Compute neighbors
    # For the moment, make simple assumptions:
    # Find only first and second neighbors
    # Assume the same number of these for every site
    # Store them in nns[a,i] and nnns[a,i], here i is the site, and a runs over the nns and nnns of that site

    # First figure out how many there are and at what distance
    alldist = Array{Float64}(0)
    for i in 2:nsites
        push!(alldist,distance(sites[:,1],sites[:,i]))
    end
    dists = unique(alldist)
    sort!(dists)
    nndist = dists[1]
    nnndist = dists[2]
    numnns = length(find(x->x==nndist,alldist))
    numnnns = length(find(x->x==nnndist,alldist))

    # Now construct all of them and store in two 2d arrays
    thenns = Array{Int64}(numnns,nsites)
    thennns = Array{Int64}(numnnns,nsites)
    for i in 1:nsites
        dists = [distance(sites[:,i],sites[:,j]) for j in 1:nsites]
        thenns[:,i] = find(x->x==nndist,dists)
        thennns[:,i] = find(x->x==nnndist,dists)
    end

    # Also compute adjacency matrix: entries = 1 if NN, 2 if NNN, 0 otherwise
    adjacency = Array{Int64}(nsites,nsites)
    for i in 1:nsites
        for j in 1:nsites
            dist = distance(sites[:,i],sites[:,j])
            if dist == nndist
                adjacency[i,j]=1
            elseif dist == nnndist
                adjacency[i,j]=2
            else
                adjacency[i,j]=0
            end
        end
    end

    Lattice(sites,flavors,thenns,thennns,adjacency,nsites,nuc,numnns,numnnns)
end

end

max_stars_repo_name: balents/lattice

max_issues_repo_name: balents/lattice

max_forks_repo_name: balents/lattice


text:
"""
Farquhar–von Caemmerer–Berry (FvCB) model for C3 photosynthesis (Farquhar et al., 1980;
von Caemmerer and Farquhar, 1981). Direct implementation of the model.

See [`Fvcb`](@ref) for more details.

# See also

- [`Fvcb`](@ref) for the coupled assimilation / conductance model
- [`FvcbIter`](@ref) for the coupled assimilation / conductance model with an iterative resolution
- [`get_J`](@ref)
- [`photosynthesis`](@ref)

# References

Caemmerer, S. von, et G. D. Farquhar. 1981. « Some Relationships between the Biochemistry of
Photosynthesis and the Gas Exchange of Leaves ». Planta 153 (4): 376‑87.
https://doi.org/10.1007/BF00384257.

Farquhar, G. D., S. von von Caemmerer, et J. A. Berry. 1980. « A biochemical model of
photosynthetic CO2 assimilation in leaves of C3 species ». Planta 149 (1): 78‑90.

# Examples

```julia
Get the fieldnames:
fieldnames(FvcbRaw)
# Using default values for the model:
A = FvcbRaw()

A.Eₐᵥ
```
"""
struct FvcbRaw{T} <: AbstractAModel
    Tᵣ::T
    VcMaxRef::T
    JMaxRef::T
    RdRef::T
    TPURef::T
    Eₐᵣ::T
    O₂::T
    Eₐⱼ::T
    Hdⱼ::T
    Δₛⱼ::T
    Eₐᵥ::T
    Hdᵥ::T
    Δₛᵥ::T
    α::T
    θ::T
end

function FvcbRaw(;Tᵣ = 25.0, VcMaxRef = 200.0, JMaxRef = 250.0, RdRef = 0.6, TPURef = 9999., Eₐᵣ = 46390.0,
    O₂ = 210.0, Eₐⱼ = 29680.0, Hdⱼ = 200000.0, Δₛⱼ = 631.88, Eₐᵥ = 58550.0, Hdᵥ = 200000.0,
    Δₛᵥ = 629.26, α = 0.24, θ = 0.7)

    FvcbRaw(promote(Tᵣ, VcMaxRef, JMaxRef, RdRef, TPURef, Eₐᵣ, O₂, Eₐⱼ, Hdⱼ, Δₛⱼ, Eₐᵥ, Hdᵥ, Δₛᵥ, α, θ)...)
end

function inputs(::FvcbRaw)
    (:PPFD, :Tₗ, :Cᵢ)
end

function outputs(::FvcbRaw)
(:A,)
end

Base.eltype(x::FvcbRaw) = typeof(x).parameters[1]

"""
    assimilation!(leaf::LeafModels{I,E,<:FvcbRaw,<:AbstractGsModel,S},constants = Constants())

Direct implementation of the photosynthesis model for C3 photosynthesis from Farquhar–von
Caemmerer–Berry (Farquhar et al., 1980; von Caemmerer and Farquhar, 1981).

# Returns

Modify the first argument in place for A, the carbon assimilation (μmol[CO₂] m-2 s-1).

# Arguments

- `leaf::LeafModels{.,.,<:FvcbRaw,<:AbstractGsModel,.}`: A [`LeafModels`](@ref) struct holding the parameters for
the model with initialisations for:
    - `Tₗ` (°C): leaf temperature
    - `PPFD` (μmol m-2 s-1): absorbed Photosynthetic Photon Flux Density
    - `Cₛ` (ppm): surface CO₂ concentration.
    - `Dₗ` (kPa): vapour pressure difference between the surface and the air saturation
    vapour pressure in case you're using the stomatal conductance model of [`Medlyn`](@ref).
- `constants = Constants()`: physical constants. See [`Constants`](@ref) for more details

# Note

`Tₗ`, `PPFD`, `Cₛ` (and `Dₗ` if you use [`Medlyn`](@ref)) must be initialised by providing
them as keyword arguments (see examples). If in doubt, it is simpler to compute the energy
balance of the leaf with the photosynthesis to get those variables. See
[`energy_balance`](@ref) for more details.

# Examples

```julia
leaf = LeafModels(photosynthesis = FvcbRaw(), Tₗ = 25.0, PPFD = 1000.0, Cᵢ = 400.0)
# NB: we need Tₗ, PPFD and Cᵢ as inputs (see [`inputs`](@ref))

assimilation!(leaf)
leaf.status.A
leaf.status.Cᵢ

# using several time-steps:
leaf = LeafModels(photosynthesis = FvcbRaw(), Tₗ = [20., 25.0], PPFD = 1000.0, Cᵢ = [380.,400.0])
# NB: we need Tₗ, PPFD and Cᵢ as inputs (see [`inputs`](@ref))

assimilation!(leaf)
DataFrame(leaf) # fetch the leaf status as a DataFrame
```

# References

Baldocchi, Dennis. 1994. « An analytical solution for coupled leaf photosynthesis and
stomatal conductance models ». Tree Physiology 14 (7-8‑9): 1069‑79.
https://doi.org/10.1093/treephys/14.7-8-9.1069.

Duursma, R. A., et B. E. Medlyn. 2012. « MAESPA: a model to study interactions between water
limitation, environmental drivers and vegetation function at tree and stand levels, with an
example application to [CO2] × drought interactions ». Geoscientific Model Development 5
(4): 919‑40. https://doi.org/10.5194/gmd-5-919-2012.

Farquhar, G. D., S. von von Caemmerer, et J. A. Berry. 1980. « A biochemical model of
photosynthetic CO2 assimilation in leaves of C3 species ». Planta 149 (1): 78‑90.

Leuning, R., F. M. Kelliher, DGG de Pury, et E.D. Schulze. 1995. « Leaf nitrogen,
photosynthesis, conductance and transpiration: scaling from leaves to canopies ». Plant,
Cell & Environment 18 (10): 1183‑1200.

Lombardozzi, L. D. et al. 2018.« Triose phosphate limitation in photosynthesis models
reduces leaf photosynthesis and global terrestrial carbon storage ». Environmental Research
Letters 13.7: 1748-9326. https://doi.org/10.1088/1748-9326/aacf68.
"""
function assimilation!(leaf::LeafModels{I,E,<:FvcbRaw,Gs,<:MutableNamedTuples.MutableNamedTuple}, meteo = nothing,
    constants = Constants()) where {I,E,Gs}

    Tₖ =  leaf.status.Tₗ - constants.K₀
    Tᵣₖ = leaf.photosynthesis.Tᵣ - constants.K₀
    Γˢ = Γ_star(Tₖ, Tᵣₖ, constants.R) # Gamma star (CO2 compensation point) in μmol mol-1
    Km = get_km(Tₖ, Tᵣₖ, leaf.photosynthesis.O₂, constants.R) # effective Michaelis–Menten coefficient for CO2
    JMax = arrhenius(leaf.photosynthesis.JMaxRef, leaf.photosynthesis.Eₐⱼ, Tₖ, Tᵣₖ, leaf.photosynthesis.Hdⱼ, leaf.photosynthesis.Δₛⱼ, constants.R)
    VcMax = arrhenius(leaf.photosynthesis.VcMaxRef, leaf.photosynthesis.Eₐᵥ, Tₖ, Tᵣₖ, leaf.photosynthesis.Hdᵥ, leaf.photosynthesis.Δₛᵥ, constants.R)
    Rd = arrhenius(leaf.photosynthesis.RdRef, leaf.photosynthesis.Eₐᵣ, Tₖ, Tᵣₖ, constants.R)
    J = get_J(leaf.status.PPFD, JMax, leaf.photosynthesis.α, leaf.photosynthesis.θ) # in μmol m-2 s-1
    Vⱼ = J / 4
    Wⱼ = Vⱼ * (leaf.status.Cᵢ - Γˢ) / (leaf.status.Cᵢ + 2.0 * Γˢ) # also called Aⱼ
    Wᵥ = VcMax * (leaf.status.Cᵢ - Γˢ) / (leaf.status.Cᵢ + Km)
    ag = 0.
    Wₚ = (leaf.status.Cᵢ - Γˢ) * 3 * leaf.photosynthesis.TPURef / (leaf.status.Cᵢ - (1 .+ 3 * ag) * Γˢ)
    leaf.status.A = min(Wᵥ, Wⱼ, Wₚ) - Rd

    return leaf.status.A
end


text:
struct IntrinsicError <: Exception
    msg::String
end

macro intrinsic_stub(ctxname::Symbol, ex::Expr)
    esc(intrinsic_stub_m(__module__, ctxname, ex))
end

function intrinsic_stub_m(mod::Module, ctxname::Symbol, ex::Expr)
    @match ex begin
        Expr(:call, ::Symbol, Expr(:parameters, _...), _...) => error("syntax: stub cannot have kwargs")
        Expr(:call, ::Symbol, args...) => nothing
        _ => error("syntax: expect function call")
    end

    msg = "$(ex.args[1]) must be executed inside @$ctxname"
    jlfn = JLFunction(;
        name=ex.args[1],
        args=ex.args[2:end],
        body=quote
            throw($CompilerPluginTools.IntrinsicError($msg))
        end
    )

    isintrinsic_def = if isdefined(mod, :__INTRINSICS__)
        nothing
    else
        quote
            const __INTRINSICS__ = Set{Any}()
            isintrinsic(x) = x in __INTRINSICS__
        end
    end

    quote
        $isintrinsic_def
        @noinline $(codegen_ast(jlfn))
        push!(__INTRINSICS__, $(jlfn.name))
    end
end


text:
"""
    _A_inv_gain_method(parameters::Parameters1PL, latent::Latent1D)
"""
function _A_inv_gain_method(parameters::Parameters1PL, latent::Latent1D)
    return _item_expected_information(parameters, latent)
end

"""
    _A_inv_gain_method(parameters::Parameters2PL, latent::Latent1D)
"""
function _A_inv_gain_method(parameters::Parameters2PL, latent::Latent1D)
    old_exp_info = copy(parameters.expected_information)
    return LinearAlgebra.tr(
        LinearAlgebra.inv(old_exp_info + _item_expected_information(parameters, latent)) -
        LinearAlgebra.inv(old_exp_info),
    )
end

"""
    A_inv_gain_method(item::AbstractItem, examinee::AbstractExaminee)

# Description
Computes the gain in the trace of the inverse of the expected information matrix for an item.

# Arguments
- **`item::AbstractItem`**: The  item.
- **`examinee::AbstractExaminee`**: The examinee at which computing the information.

# Output
It returns a `Float64` scalar.
"""
function A_inv_gain_method(item::AbstractItem, examinee::AbstractExaminee)
    return _A_inv_gain_method(item.parameters, examinee.latent)
end


text:
using RoME
using DistributedFactorGraphs

# start with an empty factor graph object
fg = LightGraphsDFG{SolverParams}(    params=SolverParams())

# Add the first pose :x0
x0 = addVariable!(fg, :x0, Pose2)

# Add at a fixed location PriorPose2 to pin :x0 to a starting location (10,10, pi/4)
prior = addFactor!(fg, [:x0], PriorPose2( MvNormal([10; 10; 1.0/8.0], Matrix(Diagonal([0.1;0.1;0.05].^2))) ) )

# Drive around in a hexagon in the cloud
for i in 0:5
    psym = Symbol("x$i")
    nsym = Symbol("x$(i+1)")
    addVariable!(fg, nsym, Pose2)
    pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
    addFactor!(fg, [psym;nsym], pp )
end

# Alrighty! At this point, we should be able to solve locally...
# perform inference, and remember first runs are slower owing to Julia's just-in-time compiling
# Can do with graph too!
tree, smt, hist = solveTree!(fg)


text:
module MYASM1

using Modia

setLogMerge(false)

include("../src/ModelMechanic.jl")
include("../src/ModelEPSphasor.jl")

##################################

VarVolts = Model(   # Voltage dip
    Tstart = 0.0,
    dT = 0.15,
    Vmin = 0.0,
    y = output,
    equations = :[
        y = if after(Tstart) && !after(Tstart+dT); Vmin else; 1.0 end;
    ]
)

ASMtest = Model(
    LOAD = parameter | Map(value=0.0, fixed=false),
    wReference= parameter | 1.0 | Var(info = "ref. Frame"),
    Omegarated=100.0*pi,

    GridVolts   = EPSphasor.VoltageSource,
    GC          = EPSphasor.GroundVec,
    asmUnit     = EPSphasor.InductionMachine,
    transformer = EPSphasor.Transformer | Map( xt=0.10, rt=0.1/40.0, Strafo=119e6),
    meterGRID   = EPSphasor.Meter,
    meterGEN    = EPSphasor.Meter,
    varvolts    = VarVolts | Map(Vmin=0.7),

    FIXED = Mechanic.Fixed,

    J = Mechanic.Inertia | Map(
        J=3.0,
        w=Var(init=0.4),
        phi=Var(init=0.0),
        ),

    equations = :[
        J.flange_b.tau = LOAD
    ],

    connect = :[
        (GridVolts.V, varvolts.y)
        (GridVolts.p, meterGRID.n)
        (meterGRID.p, transformer.pin1)
        (transformer.pin2, meterGEN.n)
        (meterGEN.p, asmUnit.p)
        (asmUnit.Flange_a, J.flange_a)
        (asmUnit.Flange_b, FIXED.flange)

        (GridVolts.n, GC.p)
    ]
)

instModel = @instantiateModel(ASMtest,
    unitless=true,
    log=false,
    logTiming=false,
    logDetails=false,
    logCode=false,
    logStateSelection=false,
#     saveCodeOnFile="CODE.txt",
)

# include("../src/Utils.jl")
# TemplateInit(instModel)

end  # module MYGEN2

max_stars_repo_name: johhell/EPSphasor

max_issues_repo_name: johhell/EPSphasor

max_forks_repo_name: johhell/EPSphasor


text:
# ---
# title: 506. Relative Ranks
# id: problem506
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories:
# link: <https://leetcode.com/problems/relative-ranks/description/>
# hidden: true
# ---
#
# Given scores of **N** athletes, find their relative ranks and the people with
# the top three highest scores, who will be awarded medals: "Gold Medal",
# "Silver Medal" and "Bronze Medal".
#
# **Example 1:**
#
#
#
#     Input: [5, 4, 3, 2, 1]
#     Output: ["Gold Medal", "Silver Medal", "Bronze Medal", "4", "5"]
#     Explanation: The first three athletes got the top three highest scores, so they got "Gold Medal", "Silver Medal" and "Bronze Medal".
#     For the left two athletes, you just need to output their relative ranks according to their scores.
#
#
# **Note:**
#
#   1. N is a positive integer and won't exceed 10,000.
#   2. All the scores of athletes are guaranteed to be unique.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module DataModel

# this data model provide necessary definitions to work with JSON objects fetched from RestAPI
# most of these struct are one-way ready-only data pull
# means, OliveAI front-end is not used to store/update any value and is only used to query, read and display information
# last_updated field used in below structs are not client generated and reflect when RestAPI is last last_updated

# TODO: why structypes are mutable below, if these are read only structs, they need to be ordered structs

import Base: ==

using StructTypes, Dates

export SystemAlert, AnomalyDetection, AIAlert, ItemDetail, ElasticSearch, User

mutable struct SystemAlert
    id::Int64
    entity::String
    unit::String
    status::String
    last_updated::DateTime
    message::String
end
# define structType of JSON deserialization
StructTypes.StructType(::Type{SystemAlert}) = StructTypes.Mutable()
StructTypes.idproperty(::Type{SystemAlert}) = :id

mutable struct AnomalyDetection
    price::String
    purchaseRecommendation::String
    qtyOrdered::String
    qtyOnHand::String
    matchExceptionRaised::String
    last_updated::DateTime
end
# define structType of JSON deserialization
StructTypes.StructType(::Type{AnomalyDetection}) = StructTypes.Mutable()

mutable struct AIAlert
    anomalyDetection::AnomalyDetection
    purchaseRecommendation::String
end
# define structType of JSON deserialization
StructTypes.StructType(::Type{AIAlert}) = StructTypes.Mutable()

mutable struct ItemDetail
    transactionType::String
    itemID::String
    UNSPSC::String
    entity::String
    category::String
    resultType::String
    preferredItem::String
    preferredVendor::String
    onContract::String
    description::String
    message::String
end
# define structType of JSON deserialization
StructTypes.StructType(::Type{ItemDetail}) = StructTypes.Mutable()

mutable struct ElasticSearch
    txnType::String
    UID::String
    sdesc::String
    ldesc::String
    alert::String
    last_updated::DateTime
end
# define structType of JSON deserialization
StructTypes.StructType(::Type{ElasticSearch}) = StructTypes.Mutable()

mutable struct User
    id::Int64 # service-managed
    username::String
    password::String
end

==(x::User, y::User) = x.id == y.id
User() = User(0, "", "")
User(username::String, password::String) = User(0, username, password)
User(id::Int64, username::String) = User(id, username, "")
StructTypes.StructType(::Type{User}) = StructTypes.Mutable()
StructTypes.idproperty(::Type{User}) = :id

end # module


text:
const Model = Union{typeof(gaussian), typeof(normal), typeof(airydisk), typeof(moffat)}

"""
    PSFModels.fit(model, params, image, inds=axes(image); func_kwargs=(;), loss=abs2, maxfwhm=Inf, alg=LBFGS(), kwargs...)

Fit a PSF model (`model`) defined by the given `params` as a named tuple of the parameters to fit and their default values. This model is fit to the data in `image` at the specified `inds` (by default, the entire array). To pass extra keyword arguments to the `model` (i.e., to "freeze" a parameter), pass them in a named tuple to `func_kwargs`. The default loss function is the chi-squared loss, which uses the the square of the difference (i.e., the L2 norm). You can change this to the L1 norm, for example, by passing `loss=abs`. The maximum FWHM can be set with `maxfwhm` as a number or tuple.

Additional keyword arguments, as well as the fitting algorithm `alg`, are passed to `Optim.optimize`. By default we use forward-mode auto-differentiation (AD) to derive Jacobians for the LBFGS optimization algorithm. Refer to the [Optim.jl documentation](https://julianlsolvers.github.io/Optim.jl/stable/) for more information.

# Choosing parameters

The `fit` function is very powerful because it gives you a great amount of flexibility in the way you fit your models. To demonstrate this, let's start with a simple isotropic [`gaussian`](@ref).

```julia
model = gaussian
# match params to arguments of PSF
params = (x=20, y=20, fwhm=3, amp=1)
```

Note that `params` can follow any order

```julia
params = (amp=1, x=20, y=20, fwhm=3)
```

Now, to extend this interface to the bivariate PSF case, where `fwhm` is two values, all you need to do is use a tuple or vector

```julia
params = (x=20, y=20, fwhm=(3, 3))
```

and, again, the order does not matter

```julia
model = moffat
params = (alpha=1, x=20, y=20, fwhm=3, amp=10)
```

# Fitting a PSF

After selecting your model and parameters, fitting data is easy

```julia
P = (x=12, y=13, fwhm=13.2, amp=0.1)
psf = gaussian.(CartesianIndicies(1:25, 1:15); P...)

params, synthpsf = PSFModels.fit(gaussian, P, psf)
```

here `params` is a named tuple of the best fitting parameters. It will not include any fixed parameters.

`synthpsf` is the best-fitting model, for direct comparison with the input data.

```julia
psf_fit = synthpsf.(CartesianIndicies(psf))
```

To alter parameters without fitting them (i.e., "freeze" them) use `func_kwargs`

```julia
P = (x=12, y=13, fwhm=(12.4, 13.2), amp=0.1)
func_kwargs = (alpha=2)
params, synthpsf = PSFModels.fit(moffat, P, psf; func_kwargs)
```
"""
function fit(model::Model,
             params,
             image::AbstractMatrix{T},
             inds=axes(image);
             func_kwargs=(;),
             loss=abs2,
             alg=LBFGS(),
             maxfwhm=Inf,
             kwargs...) where T
    _keys = keys(params)
    cartinds = CartesianIndices(inds)
    minind = map(minimum, inds)
    maxind = map(maximum, inds)
    function _loss(X::AbstractVector{T}) where T
        P = generate_params(_keys, X)
        minind[1] - 0.5 ≤ P.x ≤ maxind[1] + 0.5 || return T(Inf)
        minind[2] - 0.5 ≤ P.y ≤ maxind[2] + 0.5 || return T(Inf)
        all(0 .< P.fwhm .< maxfwhm) || return T(Inf)
        if :ratio in _keys
            0 < P.ratio < 1 || return T(Inf)
        end
        if :theta in _keys
            -45 < P.theta < 45 || return T(Inf)
        end
        # mean of errors (by default, with L2 norm == chi-squared)
        mse = mean(cartinds) do idx
            resid = model(T, idx; P..., func_kwargs...) - image[idx]
            return loss(resid)
        end
        return mse
    end
    X0 = vector_from_params(T, params)
    result = optimize(_loss, X0, alg; autodiff=:forward, kwargs...)
    Optim.converged(result) || @warn "optimizer did not converge" result
    X = Optim.minimizer(result)
    P_best = generate_params(_keys, X)
    return P_best, model(T; P_best..., func_kwargs...)
end

function vector_from_params(T, params)
    _keys = keys(params)
    _vals = values(params)
    if :fwhm in _keys && params.fwhm isa BivariateLike
        _ind = findfirst(==(:fwhm), _keys)
        first_half = @views _vals[begin:_ind - 1]
        fwhmx, fwhmy = _vals[_ind]
        second_half = @views _vals[_ind + 1:end]
        return T[first_half..., fwhmx, fwhmy, second_half...]
    end
    if :theta in _keys && !(:fwhm isa BivariateLike)
        throw(ArgumentError("cannot fit theta for isotropic distribution"))
    end
    return collect(T, _vals)
end

function generate_params(names, values)
    if length(values) > length(names) && :fwhm in names
        _ind = findfirst(==(:fwhm), names)
        fwhm = values[_ind], values[_ind + 1]
        first_half = @views zip(names[begin:_ind - 1], values[begin:_ind - 1])
        second_half = @views zip(names[_ind + 1:end], values[_ind + 2:end])
        return (;first_half..., fwhm, second_half...)
    end
    return (;zip(names, values)...)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# this is used to test that relative include paths work on other processes

i27429 = 27429

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia