text:
#=##############################################################################
# DESCRIPTION
    Methods for VTK Legacy formatting
# AUTHORSHIP
  * Author    : Eduardo J Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Nov 2017
  * License   : MIT License
=###############################################################################

"""
  `lines2vtk(lines; values=nothing)`

  Receives an array of lines, with lines[i] the i-th line, and lines[i][j] the
  j-th point in the i-th line, and formats them in vtk format. Give it point
  data through `values` in the same format, and it will output the data formated
  as corresponding to the vtk lines.

  ** Output **
  * (points, vtk_lines, vtk_values) where `points` is an array with all points
    in all lines, `lines` is an array of arrays containing the indices of
    every point in each line, and `vtk_values` is an array with the point data
    of each point in `points`.
"""
function lines2vtk(lines; values=nothing)
  points = Array{Float64,1}[]
  vtk_lines = Array{Int64,1}[]

  if values!=nothing
    vtk_values = []
  else
    vtk_values = nothing
  end

  n_points = 0
  for (i,line) in enumerate(lines)
    vtk_line = []
    for (j,point) in enumerate(line)
      push!(points, point)
      push!(vtk_line, n_points)
      if values!=nothing
        push!(vtk_values, values[i][j])
      end
      n_points += 1
    end
    push!(vtk_lines, vtk_line)
  end

  return (points, vtk_lines, vtk_values)
end

"""
  `generateCells(line1::Array{Array{Float64,1}},line2::Array{Array{Float64,1}};
                    point_data1, point_data2)`

Given two lines with the same amount of points, it generate cells by matching
the points between lines. For instance:

  ```julia
  julia> line1 = [p11,p12,p13,p14,p15,p16,p17]
  julia> line2 = [p21,p22,p23,p24,p25,p26,p27]
  julia> generateCells(line1, line2)
  ([cell1, cell2, ...], point_data)
  ```
where `cell1=[p11, p12, p22, p21]`, `cell2=[p12,p13,p23,p22]`, etc.

Give it point data corresponding to points in each line through `point_data1`
and `point_data2`, and it return it through `point_data` formatted for
`lines2vtk`'s `values` input.
"""
function generateCells(line1, line2; point_data1=nothing, point_data2=nothing)
  # Error case
  if size(line1)!=size(line2)
    error("Cells can't be generated from lines of unequal divisions!")
  end

  ncells = size(line1)[1]-1   # Number of cells
  cells = Array{Float64,1}[]  # Cells
                              # Point data on each cell
  point_data = (nothing in [point_data1, point_data2]) ? nothing : []

  # Builds cells
  for i in 1:ncells
    push!(cells, [line1[i], line1[i+1], line2[i+1], line2[i]])
    if point_data!=nothing
      push!(point_data, [point_data1[i], point_data1[i+1], point_data2[i+1],
                        point_data2[i]])
    end
  end

  return cells, point_data
end

"""
  `lines2vtkcells(line1::Array{Array{Float64,1}},line2::Array{Array{Float64,1}};
                    point_data1, point_data2)`

Given two lines with the same amount of points, it generate cells in VTK format
by matching the points between lines. For instance:

  ```julia
  julia> line1 = [p11,p12,p13]
  julia> line2 = [p21,p22,p23]
  julia> lines2vtkcells(line1, line2)
  (points, vtk_cells, point_data)
  ```
where `points=[p11,p12,p13,p21,p22,p23]`, and `vtk_cells=[[0,1,4,3],[1,2,5,4]]`.

Give it point data corresponding to points in each line through `point_data1`
and `point_data2`, and it return it through `point_data` formatted for
`generateVTK`'s `point_data` input.

Prefer this method over `generateCells()` since it will store the data
efficiently when generating the VTK.
"""
function lines2vtkcells(line1, line2; point_data1=nothing, point_data2=nothing)
  # Error case
  if size(line1)!=size(line2)
    error("Cells can't be generated from lines of unequal divisions!")
  end

  npoints = size(line1)[1]        # Number of points per line
  ncells = npoints-1              # Number of cells
  points = vcat(line1, line2)     # Points
  vtk_cells = Array{Int64,1}[]    # Cells
  # Point data
  if !(nothing in [point_data1, point_data2])
    if size(point_data1)!=size(point_data2)
      error("Invalid point data! "*
            "$(size(point_data1))!=$(size(point_data2))"*
            "(size(point_data1)!=size(point_data2))")
    else
      point_data = vcat(point_data1, point_data2)
    end
  else
    point_data = nothing
  end

  # Builds cells
  for i in 1:ncells
    point_index = i-1
    push!(vtk_cells, [point_index, point_index+1,
                        npoints+point_index+1, npoints+point_index])
  end

  return points, vtk_cells, point_data
end

"""
  `lines2vtkmulticells(line1, line2,
                          sections::Array{Tuple{Float64,Int64,Float64,Bool},1},
                          point_data1=nothing, point_data2=nothing)`

Like `lines2vtkcells()` it generates cells between two lines, but allows to
define multiple rows of cells in between. The rows are given by `sections`,
which is the same variable `sections` defined in `VTKtools_geometry.jl`'s
`multidiscretize()` function, which divides the space in between the lines.
For more details see docstring of `lines2vtkcells()` and `multidiscretize()`.
"""
function lines2vtkmulticells(line1, line2,
                          sections::Array{Tuple{Float64,Int64,Float64,Bool},1};
                          point_data1=nothing, point_data2=nothing)
  # ERROR CASES
  if size(line1)!=size(line2)
    error("Cells can't be generated from lines of unequal divisions!")
  end
  if !(nothing in [point_data1, point_data2])
    if size(point_data1)!=size(point_data2)
      error("Invalid point data! "*
            "$(size(point_data1))!=$(size(point_data2))"*
            "(size(point_data1)!=size(point_data2))")
    else
      point_data = []
    end
  else
    point_data = nothing
  end

  npoints = size(line1)[1]        # Number of points per line
  points = []                     # Points
  vtk_cells = Array{Int64,1}[]    # Cells

  # Discretize the space between the lines
  f(x) = x              # Function to discretize
  xlow, xhigh = 0, 1    # The space is parameterized between 0 and 1
  row_pos = multidiscretize(f, xlow, xhigh, sections) # Position of each row
                                                      # between the lines

  # Iterates over each row
  prev_line2 = line1
  prev_point_data2 = point_data1
  for (i,this_pos) in enumerate(row_pos[2:end])

    # Creates the upper and lower lines of this row
    this_line1 = prev_line2
    this_line2 = [line1[j] + this_pos*(line2[j]-line1[j]) for j in 1:npoints]

    # Linearly interpolates point data
      this_point_data1 = prev_point_data2
    if point_data!=nothing
      this_point_data2 = [point_data1[j] + this_pos*(point_data2[j]-point_data1[j]) for j in 1:npoints]
    else
      this_point_data2 = nothing
    end

    # Discretizes this row
    row = lines2vtkcells(this_line1, this_line2; point_data1=this_point_data1,
                                                  point_data2=this_point_data2)
    row_points, row_vtk_cells, row_point_data = row

    # Adds this row to the overall mesh
    if i==1 # Case of the first row of cells (adds both upper and lower bound)
      # Adds points
      points = vcat(points, row_points)
      # Adds point data
      if point_data!=nothing
        point_data = vcat(point_data, row_point_data)
      end
    else # Case of other rows (only adds upper bound)
      # Adds points
      points = vcat(points, row_points[npoints+1:end])
      # Adds point data
      if point_data!=nothing
        point_data = vcat(point_data, row_point_data[npoints+1:end])
      end
    end
    # Adds point indexing of cells
    vtk_cells = vcat(vtk_cells,
                            [cells .+ (i-1)*npoints for cells in row_vtk_cells])

    prev_line2 = this_line2
    prev_point_data2 = this_point_data2
  end

  return points, vtk_cells, point_data
end

"""
  `multilines2vtkmulticells(lines::Array{Array{Array{Float64,1 },1},1},
                  sections::Array{Array{Tuple{Float64,Int64,Float64,Bool},1},1};
                  point_datas=nothing)`

Expanding `lines2vtkmulticells()` it generates cells in multiple rows of cells
in between multiple lines. The rows in between each line are given by `sections`
which is an array of `sections` defined in `lines2vtkmulticells()` which divides
the space in between the lines.

For more details see docstring of `lines2vtkmulticells()`.
"""
function multilines2vtkmulticells(lines,
                  sections::Array{Array{Tuple{Float64,Int64,Float64,Bool},1},1};
                  point_datas=nothing)
  nlines = size(lines)[1]         # Number of lines
  npoints = size(lines[1])[1]     # Number of points on each line
  pflag = point_datas!=nothing     # Flag for point data

  # ERROR CASES
  if nlines<=1
    error("Two or more lines are required. Received $nlines.")
  elseif pflag && size(point_datas)[1]!=nlines
    error("Invalid point data. "*
            "$(size(point_datas)[1])!=$(nlines) (size(point_datas)[1]!=nlines)")
  elseif size(sections)[1]!=(nlines-1)
    error("Invalid number of sections."*
            " Expected $(nlines-1), received $(size(sections)[1])")
  else
    for line in lines
      if size(line)[1]!=npoints
        error("All lines must have the same number of points."*
            " Expected $npoints, found $(size(line)[1]).")
      end
    end
  end

  # Outputs
  points = []
  vtk_cells = Array{Int64,1}[]
  point_data = pflag ? [] : nothing

  for (i,line) in enumerate(lines)
    if i==1
      nothing
    else

      cur_npoints = size(points)[1]   # Current number of points

      # Meshes this section
      out = lines2vtkmulticells(lines[i-1], line,
                              sections[i-1];
                              point_data1=pflag ? point_datas[i-1] : nothing,
                              point_data2=pflag ? point_datas[i] : nothing)

      this_points, this_vtk_cells, this_point_data = out

      # Adds this section to the overall mesh
      if i==2 # Case of the first section (adds both upper and lower bound)
        # Adds points
        points = vcat(points, this_points)
        # Adds point data
        if pflag
          point_data = vcat(point_data, this_point_data)
        end
      else # Case of any other section (only adds upper bound)
        # Adds points
        points = vcat(points, this_points[npoints+1:end])
        # Adds point data
        if point_data!=nothing
          point_data = vcat(point_data, this_point_data[npoints+1:end])
        end
      end
      # Adds point indexing of cells
      this_zero = cur_npoints - npoints*(i!=2) # Point index corresponding to
                                               # this cell's zero
      vtk_cells = vcat(vtk_cells, [cells.+this_zero for cells in this_vtk_cells])
    end
  end

  return points, vtk_cells, point_data
end

"""
    `read_vtk(filename; path="")`

Read the VTK legacy file `filename` in the directory `path`, and returns
`(points, cells, cell_types, data)`.
"""
function read_vtk(filename::String; path::String="", fielddata=Dict())

    f = open(joinpath(path, filename), "r")

    # -------------- HEADER ------------------------------------------
    version = readline(f)             # VTK version
    header = readline(f)              # Header
    format = readline(f)              # VTK format

    if format != "ASCII"
        error("Only ASCII files currently supported; found $format.")
    end

    dataset_type = readline(f)[9:end] # Dataset type

    if dataset_type != "UNSTRUCTURED_GRID"
        error("Only UNSTRUCTURED_GRID dataset type current supported;"*
                    " found $dataset_type.")
    end

    # -------------- POINTS -------------------------------------------
    ln = readline(f)

    # Read fields
    while ln[1:5]=="FIELD"      # Iterate over fields
        dataname = split(ln, " ")[2]    # Data name
        narr = Base.parse(Int, split(ln, " ")[3])  # Number of arrays

        fielddata[dataname] = Dict{String, Array{Float64}}()

        for ai in 1:narr    # Read arrays
            ln = readline(f)
            splt = split(ln, " ")

            arrayname = splt[1]         # Array name
            nc = Base.parse(Int, splt[2])    # Number of components
            nt = Base.parse(Int, splt[3])    # Number of tuples
            datatype = splt[4]          # Data type

            if datatype!="double"
                @warn("Reading data type $datatype as Float64!")
            end

            fielddata[dataname][arrayname] = Array{Float64}(undef, nt, nc) #pre-allocate

            for ti in 1:nt      # Read tuples
                comps = Base.parse.(Float64, split(readline(f), " ")) # Read components
                if length(comps)!=nc
                    error("LOGIC ERROR! Expected $nc components, got $(length(comps)).")
                end
                fielddata[dataname][arrayname][ti, :] = comps
            end
        end

        ln = readline(f)
    end

    if ln[1:6]!="POINTS"
        error("Expected to find POINTS, found $(ln[1:7]).")
    end

    np = Base.parse(Int, split(ln, " ")[2])  # Number of points
    preal = split(ln, " ")[3]           # Point data real type

    # Read points
    points = Array{Float64}(undef, 3, np) #pre-allocate
    for pi in 1:np
        points[:, pi] .= Base.parse.(Float64, split(readline(f), " "))
    end

    # -------------- CELLS --------------------------------------------

    # Skip empty lines
    ln = skip_empty_lines(f)

    if ln[1:5]!="CELLS"
        error("Expected to find CELLS, found $(ln[1:7]).")
    end

    nc = Base.parse(Int, split(ln, " ")[2])      # Number of cells
    csize = Base.parse(Int, split(ln, " ")[3])   # Cell list size

    # Read cells
    cells = Array{Int, 1}[]
    for ci in 1:nc
        aux = readline(f)
        while occursin("  ", aux); aux = replace(aux, "  "=>" "); end;
        ln = Base.parse.(Int, split(aux, " "))
        nn = ln[1]                     # Number of nodes
        push!(cells, ln[2:end])
    end

    # Skip empty lines
    ln = skip_empty_lines(f)

    if ln[1:10]!="CELL_TYPES"
        error("Expected to find CELL_TYPES, found $(ln[1:7]).")
    end

    nc2 = Base.parse(Int, split(ln, " ")[2])      # Number of cells again
    if nc != nc2
        error("Found $nc2 cell types but there's $nc cells.")
    end

    # Read cells
    cell_types = Vector{Int}(undef, nc2) #pre-allocate
    for ci in 1:nc2
        cell_types[ci] = Base.parse(Int, readline(f))    # Cell type
    end

    # -------------- DATA FIELDS ------------------------------------------
    data = Dict()

    # Skip empty lines
    ln = skip_empty_lines(f)

    while ln != nothing    # Iterate until end of file

        dataparent = split(ln, " ")[1]  # Parent of this data
        nd = Base.parse(Int, split(ln, " ")[2])   # Number of data entries

        if dataparent == "POINT_DATA"
            if nd != np
                error("Found $nd point data but there's $np points.")
            end

        elseif dataparent == "CELL_DATA"
            if nd != nc
                error("Found $nd cell data but there's $nc cells.")
            end

        else
            error("Found invalid data parent $dataparent."*
                   " Valid types are POINT_DATA and CELL_DATA.")
        end

        if dataparent in keys(data)
            error("LOGIC ERROR: Found data parent $dataparent more than once!")
        end

        data[dataparent] = Dict()

        ln = skip_empty_lines(f)

        # Read datasets until next data parent is encountered
        while ln != nothing && split(ln, " ")[1] in ["SCALARS", "VECTORS"]

            datatype, name, dreal = split(ln, " ")

            if name in keys(data[dataparent])
                error("LOGIC ERROR: Found data name $name more than once!")
            end

            if datatype == "SCALARS"

                entrytype, tag = split(readline(f), " ")

                if entrytype != "LOOKUP_TABLE"
                    error("Only LOOKUP_TABLE currently supported; found $entrytype.")
                end

                this_data = Vector{Float64}(undef, nd) #pre-allocate
                for i in 1:nd
                     this_data[i] = Base.parse(Float64, readline(f))
                end
            elseif datatype == "VECTORS"

                this_data = Array{Float64}(undef, 3, nd) #pre-allocate
                for i in 1:nd
                     this_data[:, i] = Base.parse.(Float64, split(readline(f), " "))
                end

            else
                error("Found invalid data type $datatype."*
                       " Valid types are SCALARS and VECTORS.")
            end

            data[dataparent][name] = this_data

            ln = skip_empty_lines(f)
        end

    end

    close(f)

    return points, cells, cell_types, data
end

"""
  `generateVTK(filename, points; lines, cells, point_data, path, num, time)`

Generates a vtk file with the given data.

  **Arguments**
  * `points::Array{Array{Float64,1},1}`     : Points to output.

  **Optional Arguments**
  * `lines::Array{Array{Int64,1},1}`  : line definitions. lines[i] contains the
                            indices of points in the i-th line.
  * `cells::Array{Array{Int64,1},1}`  : VTK polygons definiton. cells[i]
                            contains the indices of points in the i-th polygon.
  * `data::Array{Dict{String,Any},1}` : Collection of data point fields in the
                            following format:
                              data[i] = Dict(
                                "field_name" => field_name::String
                                "field_type" => "scalar" or "vector"
                                "field_data" => point_data
                              )
                            where point_data[i] is the data at the i-th point.

See `examples.jl` for an example on how to use this function.
"""
function generateVTK(filename::String, points;
                    lines::Array{Array{Int64,1},1}=Array{Int64,1}[],
                    cells::Array{Array{Int64,1},1}=Array{Int64,1}[],
                    point_data=nothing, cell_data=nothing,
                    num=nothing, time=nothing,
                    path="", comments="", _griddims::Int64=-1,
                    keep_points::Bool=false,
                    override_cell_type::Int64=-1,
                    rnd_d=32)

  aux = num!=nothing ? ".$num" : ""
  ext = aux*".vtk"
  if !isfile(path)
    mkpath(path)
  end
  f = open(joinpath(path, filename*ext), "w")

  # HEADER
  header = "# vtk DataFile Version 4.0" # File version and identifier
  header = string(header, "\n", " ", comments) # Title
  header = string(header, "\n", "ASCII") # File format
  header = string(header, "\n", "DATASET UNSTRUCTURED_GRID")
  write(f, header)

  # TIME
  if time!=nothing
    line0 = "\nFIELD FieldData 1"
    line1 = "\nSIM_TIME 1 1 double"
    line2 = "\n$(round.(time, digits=rnd_d))"
    write(f, line0*line1*line2)
  end

  np = size(points)[1]
  nl = size(lines)[1]
  nc = size(cells)[1]

  _keep_points = keep_points || (nl==0 && nc==0)

  # POINTS
  write(f, string("\n", "POINTS ", np, " float"))
  for i in 1:np
    print(f, "\n", round.(points[i][1], digits=rnd_d), " ",
            round.(points[i][2], digits=rnd_d), " ", round.(points[i][3], digits=rnd_d))
  end

  # We do this to avoid outputting points as cells if outputting a Grid
  # or if we simply want to ignore points
  if _griddims!=-1 || !_keep_points
    auxnp = np
    np = 0
  end

  # CELLS
  auxl = size(lines)[1]
  for line in lines
    auxl += size(line)[1]
  end
  auxc = size(cells)[1]
  for cell in cells
    auxc += size(cell)[1]
  end
  write(f, "\n\nCELLS $(np+nl+nc) $(2*np+auxl+auxc)")

  for i in 1:np+nl+nc
    if i<=np
      pts = [i-1]
    elseif i<=np+nl
      pts = lines[i-np]
    else
      pts = cells[i-(nl+np)]
    end
    print(f, "\n", size(pts,1))
    for pt in pts
      print(f, " ", pt)
    end
  end

  write(f, "\n\nCELL_TYPES $(np+nl+nc)")
  for i in 1:np+nl+nc
    if i<=np
      tpe = 1
    elseif i<=np+nl
      tpe = 4
    else
      if override_cell_type==-1
        if _griddims!=-1
          if _griddims==1
            tpe = 3
          elseif _griddims==2
            tpe = 9
          elseif _griddims==3
            tpe = 12
          else
            error("Generation of VTK cells of $_griddims dimensions not implemented")
          end
        else
          tpe = 7
        end
      else
        tpe = override_cell_type
      end
    end
    print(f, "\n", tpe)
  end

  if _griddims!=-1 || !_keep_points
    np = auxnp
  end

  # POINT DATA
  if point_data!=nothing
      write(f, "\n\nPOINT_DATA $np")
  end
  _p_data = point_data!=nothing ? point_data : []
  for field in _p_data
    field_name = field["field_name"]
    field_type = field["field_type"]
    data = field["field_data"]
    if size(data)[1]!=np
      @warn("Corrupted field $(field_name)! Field size != number of points.")
    end
    if field_type=="scalar"
      write(f, "\n\nSCALARS $field_name float\nLOOKUP_TABLE default")
      for entry in data
        print(f, "\n", round.(entry, digits=rnd_d))
      end
    elseif field_type=="vector"
      write(f, "\n\nVECTORS $field_name float")
      for entry in data
        print(f, "\n", round.(entry[1], digits=rnd_d), " ", round.(entry[2], digits=rnd_d), " ", round.(entry[3], digits=rnd_d))
      end
    else
      error("Unknown field type $(field_type).")
    end
  end

    # CELL DATA
    if cell_data!=nothing
        write(f, "\n\nCELL_DATA $nc")
    end
    _c_data = cell_data!=nothing ? cell_data : []
    for field in _c_data
      field_name = field["field_name"]
      field_type = field["field_type"]
      data = field["field_data"]
      if size(data)[1]!=nc
        @warn("Corrupted field $(field_name)! Field size != number of cells.")
      end
      if field_type=="scalar"
        write(f, "\n\nSCALARS $field_name float\nLOOKUP_TABLE default")
        for entry in data
          print(f, "\n", round.(entry, digits=rnd_d))
        end
      elseif field_type=="vector"
        write(f, "\n\nVECTORS $field_name float")
        for entry in data
          print(f, "\n", round.(entry[1], digits=rnd_d), " ", round.(entry[2], digits=rnd_d), " ", round.(entry[3], digits=rnd_d))
        end
      else
        error("Unknown field type $(field_type).")
      end
    end

  close(f)
  return filename*ext*";"
end


text:
include("velocity.jl")

mutable struct PSO <: AbstractParameters
    N::Int
    C1::Float64
    C2::Float64
    ω::Float64
    v::Array{Float64} # velocity
    flock::Array
end

"""
    PSO(;
        N  = 0,
        C1 = 2.0,
        C2 = 2.0,
        ω  = 0.8,
        information = Information(),
        options = Options()
    )

Parameters for Particle Swarm Optimization (PSO) algorithm: learning rates `C1` and `C2`,
`N` is the population size and `ω` controls the inertia weight.
# Example

```jldoctest
julia> f(x) = sum(x.^2)
f (generic function with 1 method)

julia> optimize(f, [-1 -1 -1; 1 1 1.0], PSO())
+=========== RESULT ==========+
  iteration: 1000
    minimum: 1.40522e-49
  minimizer: [3.0325415595139883e-25, 1.9862212295897505e-25, 9.543772256546461e-26]
    f calls: 30000
 total time: 0.1558 s
+============================+

julia> optimize(f, [-1 -1 -1; 1 1 1.0], PSO(N = 100, C1=1.5, C2=1.5, ω = 0.7))
+=========== RESULT ==========+
  iteration: 300
    minimum: 2.46164e-39
  minimizer: [-3.055334698085433e-20, -8.666986835846171e-21, -3.8118413472544027e-20]
    f calls: 30000
 total time: 0.1365 s
+============================+
```

"""
function PSO(;
    N::Int = 0,
    C1 = 2.0,
    C2 = 2.0,
    ω = 0.8,
    v = Float64[],
    flock = xf_indiv[],
    information = Information(),
    options = Options(),
)

parameters = PSO(N, promote(Float64(C1), C2, ω)..., v, flock)

Algorithm(
    parameters,
    information = information,
    options = options,
)
end

function update_state!(
        status,
        parameters::PSO,
        problem::AbstractProblem,
        information::Information,
        options::Options,
        args...;
        kargs...
    )
    xGBest = get_position(status.best_sol)

    X_new = zeros(parameters.N, size(problem.bounds, 2))

    # For each elements in population
    for i in 1:parameters.N
        x = get_position(parameters.flock[i])
        xPBest = get_position(status.population[i])
        parameters.v[i, :] = velocity(x, parameters.v[i, :], xPBest, xGBest, parameters)
        x += parameters.v[i, :]
        reset_to_violated_bounds!(x, problem.bounds)
        X_new[i,:] = x
    end

    for (i, sol) in enumerate(create_solutions(X_new, problem;ε = options.h_tol))
        if is_better(sol, status.population[i])
            status.population[i] = sol
            if is_better(sol, status.best_sol)
                status.best_sol = sol
                xGBest = get_position(status.best_sol)
            end
        end

        parameters.flock[i] = sol

        # stop condition
        stop_criteria!(status, parameters, problem, information, options)
        status.stop && break
    end

end

function initialize!(
    status,
    parameters::PSO,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
       )
    D = size(problem.bounds, 2)

    if parameters.N == 0
        parameters.N = 10 * D
    end

    if options.f_calls_limit == 0
        options.f_calls_limit = 10000D
        options.debug &&
        @warn( "f_calls_limit increased to $(options.f_calls_limit)")
    end

    if options.iterations == 0
        options.iterations = div(options.f_calls_limit, parameters.N) + 1
    end

    status = gen_initial_state(problem,parameters,information,options,status)

    parameters.v = zeros(parameters.N, D)

    parameters.flock = status.population

    status

end

function final_stage!(
    status,
    parameters::PSO,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
    )
    status.final_time = time()

end

ext: jl


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

n = 200 #3000 #1000 #20000
t = 0.35 #0.15 #0.4 #0.15

function randlines(n=1000,t=0.2)
	V = zeros(Float64,2,2*n)
	EV = [zeros(Int64,2) for k=1:n]
	for k=1:n
		v1 = rand(Float64,2)
		v2 = rand(Float64,2)
		vm = (v1+v2)/2
		transl = rand(Float64,2)
		V[:,k] = (v1-vm)*t + transl
		V[:,n+k] = (v2-vm)*t + transl
		EV[k] = [k,n+k]
	end
	V = GL.normalize2(V)
	model = (V,EV)
	return model
end

model = randlines(n,t)
Sigma = Lar.spaceindex(model)
GL.VIEW([ GL.GLLines(model...)] )

W,EW = Lar.fraglines(1.5,1.5,1.5)(model)
GL.VIEW([ GL.GLLines(W,EW,GL.COLORS[1])] );

W,EW = Lar.fragmentlines(model)
V,EVs = Lar.biconnectedComponent((W,EW::Lar.Cells)) # 2-connected components (H & T)

graphs = [ GL.GLLines(V,EVs[k],GL.COLORS[(k-1)%12+1]) for k=1:length(EVs) ]
GL.VIEW(graphs);


text:
module ZMQServer

using ..SceneTrees: SceneNode

using HttpServer
using WebSockets
using ZMQ

export ZMQWebSocketBridge, zmq_url, web_url

const VIEWER_ROOT = abspath(joinpath(@__DIR__, "..", "..", "viewer", "static"))
const VIEWER_HTML = "meshcat.html"
const DEFAULT_FILESERVER_PORT = 7000
const MAX_ATTEMPTS = 1000
const DEFAULT_ZMQ_METHOD = "tcp"
const DEFAULT_ZMQ_PORT = 6000

function handle_file_request(req, res)
    # This file handler is *extremely* simple, by design. I'm not currently
    # confident that I can write a secure file server that only serves files
    # inside the viewer root, so instead I am manually whitelisting the files
    # that might need to be served.
    if req.resource == "/static" || req.resource == "/static/" || req.resource == joinpath("/static", VIEWER_HTML)
        file = open(joinpath(VIEWER_ROOT, VIEWER_HTML))
    elseif req.resource in [
        "/static/js/LoaderSupport.js",
        "/static/js/OBJLoader.js",
        "/static/js/OBJLoader2.js",
        "/static/js/OrbitControls.js",
        "/static/js/dat.gui.js",
        "/static/js/meshcat.js",
        "/static/js/msgpack.min.js",
        "/static/js/split.min.js",
        "/static/js/three.js",
        "/static/dat.gui/build/dat.gui.js"
        ]
        file = open(joinpath(VIEWER_ROOT, req.resource[9:end]))
    else
        return Response(404)
    end
    return Response(read(file))
end

function find_available_port(f::Function, default=8000, max_attempts=1000)
    for i in 1:max_attempts
        port = default + i - 1
        try
            return f(port), port
        catch e
            if e isa Base.UVError || e isa ZMQ.StateError
                # info("Port $port in use, trying another")
            else
                rethrow(e)
            end
        end
    end
    error("Could not find an available port in the range [$default, $(default+max_attempts))")
end

mutable struct ZMQWebSocketBridge
    host::IPv4
    websockets::Set{WebSocket}
    tree::SceneNode
    web_port::Int
    web_server::Server
    zmq_context::ZMQ.Context
    zmq_socket::ZMQ.Socket
    zmq_url::String

    function ZMQWebSocketBridge(zmq_url=nothing, host::IPv4=IPv4(127,0,0,1), port=nothing)
        HttpServer.initcbs()
        bridge = new(host, Set{WebSocket}(), SceneNode())

        if port === nothing
            bridge.web_server, bridge.web_port = find_available_port(DEFAULT_FILESERVER_PORT, MAX_ATTEMPTS) do port
                start_server(bridge, port)
            end
        else
            bridge.web_port = port
            bridge.web_server = start_server(bridge, port)
        end
        yield()

        bridge.zmq_context = Context()

        if zmq_url === nothing
            (bridge.zmq_socket, bridge.zmq_url), _ = find_available_port(DEFAULT_ZMQ_PORT, MAX_ATTEMPTS) do port
                zmq_url = "$(DEFAULT_ZMQ_METHOD)://$(bridge.host):$port"
                start_zmq(bridge, zmq_url), zmq_url
            end
        else
            bridge.zmq_url = zmq_url
            bridge.zmq_socket = start_zmq(bridge, zmq_url)
        end
        println("zmq_url=", zmq_url)
        println("web_url=", web_url(bridge))
        bridge
    end
end

zmq_url(b::ZMQWebSocketBridge) = b.zmq_url
web_url(b::ZMQWebSocketBridge) = "http://$(b.host):$(b.web_port)/static/"

function WebSockets.WebSocketHandler(bridge::ZMQWebSocketBridge)
    WebSocketHandler() do req, client
        push!(bridge.websockets, client)
        send_scene(bridge, client)
        try
            while true
                read(client)
            end
        catch e
            if e isa WebSockets.WebSocketClosedError
                if client in bridge.websockets
                    delete!(bridge.websockets, client)
                end
            end
        end
    end
end

function start_server(bridge::ZMQWebSocketBridge, port)
    server = Server(HttpHandler(handle_file_request), WebSocketHandler(bridge))
    listen(server, bridge.host, port)
    @async HttpServer.handle_http_request(server)
    return server
end

function start_zmq(bridge::ZMQWebSocketBridge, url::AbstractString)
    socket = Socket(bridge.zmq_context, REP)
    ZMQ.bind(socket, url)
    socket
end

function wait_for_websockets(bridge::ZMQWebSocketBridge)
    while length(bridge.websockets) == 0
        sleep(0.1)
    end
end

function recv_multipart(sock::ZMQ.Socket)
    frames = [ZMQ.recv(sock)]
    while ZMQ.ismore(sock)
        push!(frames, ZMQ.recv(sock))
    end
    frames
end

function send_to_websockets(bridge::ZMQWebSocketBridge, msg)
    @sync begin
        for websocket in bridge.websockets
            @async begin
                if isopen(websocket)
                    write(websocket, msg)
                end
            end
        end
    end
end

function update_tree!(bridge::ZMQWebSocketBridge, command::String, path::AbstractVector, data)
    if command == "set_object"
        bridge.tree[path].object = data
    elseif command == "set_transform"
        bridge.tree[path].transform = data
    else
        @assert command == "delete"
        if length(path) == 0
            bridge.tree = SceneNode()
        else
            delete!(bridge.tree, path)
        end
    end
end

function send_scene(bridge::ZMQWebSocketBridge, websocket::WebSocket)
    @sync begin
        foreach(bridge.tree) do node
            if !isnull(node.object)
                @async write(websocket, get(node.object))
            end
            if !isnull(node.transform)
                @async write(websocket, get(node.transform))
            end
        end
    end
end

function Base.run(bridge::ZMQWebSocketBridge)
    if !isfile(joinpath(VIEWER_ROOT, VIEWER_HTML))
        error("Could not find `viewer/static/meshcat.html`. Please run `Pkg.build(\"MeshCat\")`.")
    end
    while true
        frames = recv_multipart(bridge.zmq_socket)
        command = unsafe_string(frames[1])
        if command == "url"
            ZMQ.send(bridge.zmq_socket, web_url(bridge))
        elseif command == "wait"
            wait_for_websockets(bridge)
            ZMQ.send(bridge.zmq_socket, "ok")
        elseif command in ["set_object", "set_transform", "delete"]
            path = filter(x -> length(x) > 0, split(unsafe_string(frames[2]), '/'))
            data = take!(convert(IOStream, frames[3]))
            send_to_websockets(bridge, data)
            update_tree!(bridge, command, path, data)
            ZMQ.send(bridge.zmq_socket, "ok")
        else
            ZMQ.send(bridge.zmq_socket, "error: unrecognized command")
        end
    end
end

end


text:
using Documenter
using NemoMod

makedocs(
    sitename = "NEMO",
    format = Documenter.HTML(),
    pages = [
        "Introduction" => "index.md"
        "Quick start" => "quick_start.md"
        "Installation" => "installation.md"
        "Model concept" => "model_concept.md"
        "Inputs" => [
            "Dimensions" => "dimensions.md"
            "Parameters" => "parameters.md"
            "Time slicing" => "time_slicing.md"
            "Units of measure" => "units.md"
        ]
        "Outputs" => [
            "Variables" => "variables.md"
        ]
        #"Mathematical model"
        "Scenario databases" => "scenario_db.md"
        "Calculating a scenario" => [
            "Julia syntax" => "calculatescenario.md"
            "Calculating selected years" => "selected_years.md"
            "Configuration files" => "configuration_file.md"
            "Performance tips" => "performance_tips.md"
        ]
        "Custom constraints" => "custom_constraints.md"
        "FAQs" => "faqs.md"
        "Function reference" => "functions.md"
        "Release notes" => "release_notes.md"
    ],
)

deploydocs(repo = "github.com/sei-international/NemoMod.jl.git")

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#

#= Markdown syntax for referring to one docstring

```@docs
#NemoMod.logmsg(msg::String, suppress=false, dtm=now()::DateTime)
```

=#


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

abstract type Permission end
struct Read <: Permission end
struct Write <: Permission end

struct Protected{T}
    data::T
    perms::NamedTuple
end

function Base.getproperty(protected::Protected, field::Symbol)
    return getfieldperm(
        getfield(protected, :data),
        field,
        get(getfield(protected, :perms), field, nothing),
    )
end

function Base.setproperty!(protected::Protected, field::Symbol, x)
    return setfieldperm!(
        getfield(protected, :data),
        field,
        x,
        get(getfield(protected, :perms), field, nothing),
    )
end

function Base.getindex(protected::Protected, idx)
    return getindexperm(
        getfield(protected, :data),
        idx,
        get(getfield(protected, :perms), field, nothing),
    )
end

function Base.setindex!(protected::Protected, idx, x)
    return setindexperm!(
        getfield(protected, :data),
        idx,
        x,
        get(getfield(protected, :perms), field, nothing),
    )
end

getindexperm(data::D, idx, perm::Nothing) where D =
    throw(ProtectedException("Read access to $(D)[$(string(idx))] disallowed."))
getindexperm(data, idx, perm::Union{Read, Write}) = getindex(data, idx)
getindexperm(data, idx, perm::NamedTuple) = protect(getindex(data, idx), perm)

setindexperm!(data::D, idx, x, perm::Union{Read, Nothing, NamedTuple}) where D =
    throw(ProtectedException("Write access to $(D)[$(string(idx))] disallowed."))
setindexperm!(data, idx, x, perm::Write) = setindex!(data, idx, x)

getfieldperm(data::D, field, perm::Nothing) where D =
    throw(ProtectedException("Read access to $(D).$(string(field)) disallowed."))
getfieldperm(data, field, perm::Union{Read, Write}) = getproperty(data, field)
getfieldperm(data, field, perm::NamedTuple) = protect(getproperty(data, field), perm)

setfieldperm!(data::D, field, x, perm::Union{Read, Nothing, NamedTuple}) where D =
    throw(ProtectedException("Write access to $(D).$(string(field)) disallowed."))
setfieldperm!(data, field, x, perm::Write) = setproperty!(data, field, x)

Base.fieldnames(protected::Protected) = fieldnames(typeof(getfield(protected, :data)))

protect(x, perms) = Protected(x, perms)

struct PropDict{V}
    d::Dict
    function PropDict(d::Dict{K, V}) where {K, V}
        return new{V}(d)
    end
end


text:
# The translation pass turns Julia expressions into SQL expressions
# Frustratingly, primary tables are translated differently from secondary tables, so translate must propagate the `primary` keyword
struct SQLExpression
    call::Symbol
    arguments::Tuple
    SQLExpression(call, arguments...) = new(call, arguments)
end

function translate(something; primary = true)
    something
end
function translate(source_row::SourceRow; primary = true)
    source_row.table_name
end
function translate(call::Expr; primary = true)
    arguments, keywords = split_call(call)
    # TODO: figure out a way for users to pass a keyword named primary
    translate_call(arguments...; primary = primary, keywords...)
end

# A 1-1 mapping between Julia functions and SQL functions
function translate_default(location, function_type, SQL_call)
    result = :(
        function translate_call($function_type, arguments...; primary = true)
            $SQLExpression($SQL_call, $map_unrolled(
                argument -> $translate(argument; primary = primary),
                arguments
            )...)
        end
    )
    result.args[2].args[1] = location
    result
end

macro translate_default(a_function, SQL_call)
    translate_default(__source__, a_function, SQL_call) |> esc
end

@translate_default ::typeof(==) :(=)

@translate_default ::typeof(!=) Symbol("<>")

@translate_default ::typeof(!) :NOT

@translate_default ::typeof(&) :AND

@translate_default ::typeof(|) :OR

@translate_default ::typeof(*) :*

@translate_default ::typeof(+) :+

@translate_default ::typeof(%) :%

@translate_default ::typeof(abs) :ABS

function as(pair; primary = true)
    SQLExpression(:AS,
        translate(pair.second.code; primary = primary),
        pair.first
    )
end

@translate_default ::Type{Char} :CHAR

@translate_default ::typeof(coalesce) :COALESCE

@translate_default ::typeof(QueryOperators.drop) :OFFSET

function translate_call(::typeof(QueryOperators.filter), iterator, call, call_expression; primary = true)
    SQLExpression(:WHERE,
        translate(iterator; primary = primary),
        translate(call(model_row(iterator)).code; primary = primary)
    )
end

function translate_call(::typeof(getproperty), source_tables::Database, table_name; primary = true)
    translated = translate(table_name; primary = primary)
    if primary
        SQLExpression(:FROM, translated)
    else
        translated
    end
end
function translate_call(::typeof(getproperty), source_row::SourceRow, column_name; primary = true)
    translated = translate(column_name; primary = primary)
    if primary
        translated
    else
        SQLExpression(:., source_row.table_name, translated)
    end
end

function translate_call(::typeof(QueryOperators.groupby), ungrouped, group_function, group_function_expression, map_selector, map_function_expression; primary = true)
    # TODO: map_selector
    model = model_row(ungrouped)
    SQLExpression(Symbol("GROUP BY"),
        SQLExpression(:AS, SQLExpression(:FROM, translate_call(
            QueryOperators.map,
            ungrouped,
            map_selector, map_function_expression,
            primary = primary
        )), :__TABLE__),
        translate(group_function(model).code; primary = primary)
    )
end

@translate_default ::typeof(if_else) :CASE

@translate_default ::typeof(in) :IN

@translate_default ::typeof(isequal) Symbol("IS NOT DISTINCT FROM")

@translate_default ::typeof(isless) :<

@translate_default ::typeof(ismissing) Symbol("IS NULL")

function translate_call(::typeof(QueryOperators.join), source1, source2, key1, key1_expression, key2, key2_expression, combine, combine_expression; primary = true)
    model_row_1 = model_row(source1)
    model_row_2 = model_row(source2)
    SQLExpression(:SELECT,
        SQLExpression(:ON,
            SQLExpression(Symbol("INNER JOIN"),
                translate(source1),
                # mark as not primary to suppress FROM
                translate(source2; primary = false)
            ),
            # mark both as not primary to always be explicit about table
            SQLExpression(:(=),
                translate(key1(model_row_1).code; primary = false),
                translate(key2(model_row_2).code; primary = false)
            )
        ),
        # mark both as not primary to always be explicit about table
        Generator(
            pair -> as(pair; primary = false),
            pairs(combine(model_row_1, model_row_2))
        )...
    )
end

@translate_default ::typeof(length) :COUNT

@translate_default ::typeof(lowercase) :LOWER

function translate_call(::typeof(QueryOperators.map), select_table, call, call_expression; primary = true)
    SQLExpression(
        Symbol("SELECT"), translate(select_table; primary = primary),
        Generator(
            pair -> as(pair; primary = primary),
            pairs(call(model_row(select_table)))
        )...
    )
end

@translate_default ::typeof(max) :max

@translate_default ::typeof(min) :min

translate_call(::typeof(occursin), needle::AbstractString, haystack; primary = true) =
    SQLExpression(
        :LIKE,
        translate(haystack; primary = primary),
        string('%', needle, '%')
    )
translate_call(::typeof(occursin), needle::Regex, haystack; primary = true) =
    SQLExpression(
        :LIKE,
        translate(haystack; primary = primary),
        # * => %, . => _
        replace(replace(needle.pattern, r"(?<!\\)\.\*" => "%"), r"(?<!\\)\." => "_")
    )
@translate_default ::typeof(occursin) :LIKE

function translate_call(::typeof(QueryOperators.orderby), unordered, key_function, key_function_expression; primary = true)
    SQLExpression(Symbol("ORDER BY"),
        translate(unordered; primary = primary),
        translate(key_function(model_row(unordered)).code; primary = primary)
    )
end

function translate_call(::typeof(QueryOperators.orderby_descending), unordered, key_function, key_function_expression; primary = true)
    SQLExpression(Symbol("ORDER BY"),
        translate(unordered; primary = primary),
        SQLExpression(:DESC,
            translate(key_function(model_row(unordered)).code; primary = primary)
        )
    )
end

translate_call(::typeof(startswith), full, prefix::AbstractString; primary = primary) =
    SQLExpression(
        :LIKE,
        translate(full),
        string(prefix, '%')
    )

@translate_default ::typeof(string) :||

@translate_default ::typeof(QueryOperators.take) :LIMIT

function translate_call(::typeof(QueryOperators.thenby), unordered, key_function, key_function_expression; primary = true)
    original = translate(unordered; primary = primary)
    SQLExpression(original.call, original.arguments...,
        translate(key_function(model_row(unordered)).code; primary = primary)
    )
end

function translate_call(::typeof(QueryOperators.thenby_descending), unordered, key_function, key_function_expression; primary = true)
    original = translate(unordered; primary = primary)
    SQLExpression(original.call, original.arguments...,
        SQLExpression(:DESC,
            translate(key_function(model_row(unordered)).code; primary = primary)
        )
    )
end

function translate_call(::typeof(QueryOperators.unique), repeated, key_function, key_function_expression; primary = true)
    result = translate(repeated; primary = primary)
    SQLExpression(Symbol(string(result.call, " DISTINCT")), result.arguments...)
end

@translate_default ::typeof(uppercase) :UPPER


text:
module Radau

using DocStringExtensions
using ForwardDiff
using StaticArrays
using LinearAlgebra
using DelimitedFiles

# using Polynomials: Poly, polyval, polyder, coeffs
# using PolynomialRoots: roots
# using GenericLinearAlgebra: eigvals
# using GenericSVD: svd

include("load_table_from_file.jl")
include("radau_struct.jl")
include("radau_functions.jl")
include("radau_utilities.jl")
include("radau_solve.jl")
include("interpolate.jl")
include("adaptive.jl")

export
    # load_table_from_file.jl

    # radau_struct.jl
    RadauIntegrator,
    RadauTable,
    RadauStep,
    RadauRule,

    # radau_functions.jl
    calcJacobian!,
    updateInvC!,
    zeroFill!,
    updateFX!,
    calcEw!,
    updateStageX!,

    # radau_utilities.jl
    makeRadauIntegrator,
    get_X_final,
    get_exponent,
    print_exit_flag,

    # radau_solve.jl
    solveRadau,

    # adaptive.jl
    calc_x̂_minus_x,
    calc_x_err_norm,
    update_x_err_norm!,
    update_h!,
    calc_h_new_estimate_1,

    # interpolate.jl
    interpolate_output!

end


text:
insts = OrderedDict()
insts["minimal"] = [
    ((2,),),
    ]
insts["fast"] = [
    ((5,),),
    ((10,),),
    ]
insts["various"] = [
    ((5,),),
    ((10,),),
    ((15,),),
    ]
return (NearestCorrelationJuMP, insts)


text:
@testset "Gain and phase mismatch and crosstalk" begin

    @testset "Static" begin
        gain_and_phase_mism_crosstalk = SMatrix{2,2}(1,2,3,4)
        @test GNSSSimulator.propagate(
            gain_and_phase_mism_crosstalk,
            1ms,
            Random.GLOBAL_RNG
        ) == gain_and_phase_mism_crosstalk

        @test get_gain_phase_mism_crosstalk(
            gain_and_phase_mism_crosstalk
        ) == SMatrix{2,2}(1,2,3,4)
    end

    @testset "Asymptotic" begin
        a = SVector(1/2, -1/2)
        b = SVector(3/4, 5/4)
        gain_phase_mism_crosstalk = @inferred AsymptoticGainPhaseMismCrosstalk(a, b, e = 2)
        @test gain_phase_mism_crosstalk.a == SVector(1/2, -1/2)
        @test gain_phase_mism_crosstalk.b == SVector(3/4, 5/4)
        @test gain_phase_mism_crosstalk.e == 2
        @test gain_phase_mism_crosstalk.C == SMatrix{2,2,ComplexF64}(I)
        @test gain_phase_mism_crosstalk.t == 0.0s

        gpmc = @inferred get_gain_phase_mism_crosstalk(gain_phase_mism_crosstalk)
        @test gpmc == Diagonal(cis.([3/4; 5/4]))

        next_gain_phase_mism_crosstalk = @inferred GNSSSimulator.propagate(
            gain_phase_mism_crosstalk,
            1000000s,
            Random.GLOBAL_RNG
        )
        next_gpmc = @inferred get_gain_phase_mism_crosstalk(next_gain_phase_mism_crosstalk)
        @test next_gpmc ≈ Diagonal(cis.([1/2, -1/2]))
    end
end


text:
using Documenter, RoadRunner

makedocs(
    modules = [RoadRunner],
    authors = "Jin Xu",
    repo="https://github.com/sunnyXu/RoadRunner.jl/blob/{commit}{path}#L{line}",
    sitename = "RoadRunner.jl",
    format=Documenter.HTML(;
    prettyurls=get(ENV, "CI", "false") == "true",
    canonical="https://sunnyXu.github.io/RoadRunner.jl",
    assets=String[],
    ),
    pages = [
        "Home" => "index.md",
    ],
    doctest = true
)

deploydocs(;
    deps = Deps.pip("mkdocs","python-markdown-math"),
    repo = "github.com/SunnyXu/RoadRunner.jl",
)


text:
using ITensors,
      Test
import Random

Random.seed!(12345)

include("2d_classical_ising.jl")

"""
trg(T::ITensor; χmax::Int, nsteps::Int) -> κ,T

Perform the TRG algorithm on the partition function composed of the ITensor T.
T is assumed to have Indices with tags "left", "right", "up", and "down".

The indices of T must obey:

`firstind(T,"left") = tags(firstind(T,"right"),"right->left")`

`firstind(T,"up") = tags(firstind(T,"down"),"down->up")`

χmax is the maximum renormalized bond dimension.

nsteps are the number of renormalization steps performed.

The outputs are κ, the partition function per site, and the final renormalized
ITensor T (also with Indices with tags "left","right","up", and "down").
"""
function trg(T::ITensor, horiz_inds, vert_inds;
             χmax::Int, nsteps::Int)

  l = horiz_inds[1]
  r = horiz_inds[2]
  u = vert_inds[1]
  d = vert_inds[2]

  @assert hassameinds((l,r,u,d),T)

  T = addtags(T,"orig")
  l = addtags(l,"orig")
  r = addtags(r,"orig")
  u = addtags(u,"orig")
  d = addtags(d,"orig")

  # Keep track of the partition function per site
  κ = 1.0

  for n = 1:nsteps
    Fr,Fl = factorize(T, (l,u); ortho="none",
                                maxdim=χmax,
                                tags="renorm")
    Fd,Fu = factorize(T, (r,u); ortho="none",
                                maxdim=χmax,
                                tags="renorm")

    Fl = addtags(Fl,"left","renorm")
    Fr = addtags(Fr,"right","renorm")
    Fu = addtags(Fu,"up","renorm")
    Fd = addtags(Fd,"down","renorm")

    Fl = replacetags(replacetags(Fl,"down","dnleft","orig"),
                                    "right","upleft","orig")
    Fu = replacetags(replacetags(Fu,"left","upleft","orig"),
                                    "down","upright","orig")
    Fr = replacetags(replacetags(Fr,"up","upright","orig"),
                                    "left","dnright","orig")
    Fd = replacetags(replacetags(Fd,"right","dnright","orig"),
                                    "up","dnleft","orig")

    T = Fl*Fu*Fr*Fd
    T = replacetags(T,"renorm","orig")

    l = firstind(T,"left")
    r = firstind(T,"right")
    u = firstind(T,"up")
    d = firstind(T,"down")

    trT = abs((T*δ(l,r)*δ(u,d))[])
    T = T/trT
    κ *= trT^(1.0/2^n)
  end
  T = removetags(T,"orig")
  l = firstind(T,"left")
  r = firstind(T,"right")
  u = firstind(T,"up")
  d = firstind(T,"down")
  return κ,T,(l,r),(u,d)
end

@testset "trg" begin
  # Make Ising model MPO
  β = 1.1*βc
  d = 2
  s = Index(d)
  l = addtags(s,"left")
  r = addtags(s,"right")
  u = addtags(s,"up")
  d = addtags(s,"down")
  T = ising_mpo((l,r),(u,d),β)

  χmax = 20
  nsteps = 20
  κ,T,(l,r),(u,d) = trg(T,(l,r),(u,d);χmax=χmax,nsteps=nsteps)

  @test κ≈exp(-β*ising_free_energy(β)) atol=1e-4
end

nothing


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ 3e0c149a-3350-11eb-0b17-7bea5d1a28d4
begin
	DAY = 6
	md"# 2020 Day $DAY"
end

# ╔═╡ 18080fe2-339b-11eb-34f0-f393638506fc
md"## Part 1"

# ╔═╡ 4e084a4e-339b-11eb-3ffc-41118215f226
md"*Count if any members in the group have said yes to a question*"

# ╔═╡ 24e1613c-339b-11eb-2976-59884053d1ee
input = read("../data/day_$DAY.txt", String)

# ╔═╡ df1d1be2-36c6-11eb-07c6-493a612dd182
test_input = """
abc

a
b
c

ab
ac

a
a
a
a

b
"""

# ╔═╡ 6cc1348c-3786-11eb-07c8-bb4f6ddfe769
part1(input) = sum(length.(Set.(replace.((split(input,"\n\n")),"\n"=>""))))

# ╔═╡ cdea65a0-3786-11eb-2444-b7d482414e20
@assert part1(test_input) == 11

# ╔═╡ d5b57b6c-3786-11eb-13da-33cee18cfa8b
@assert part1(input)  == 6590

# ╔═╡ 5f797e22-339d-11eb-1ff3-dda90f946838
md"## Part 2"

# ╔═╡ e3177eae-339e-11eb-3f9f-1f7230976832
md"*Only count if all members in the group have said yes to a question*"

# ╔═╡ 5218903e-3789-11eb-11c2-ebed36baf6fe
function part2(input::AbstractString)
	group_strings = split.((split(input,"\n\n")))
	groups = [Set.(ls) for ls in group_strings] # array with list of sets
	intersects = [intersect(group...) for group in groups]
	lengths = length.(intersects)
	return sum(lengths)
end

# ╔═╡ 6a2ab5da-3789-11eb-221f-7d829e2967b6
@assert part2(test_input) == 6

# ╔═╡ a8b55df0-3789-11eb-28c1-edc492b7060e
@assert part2(input) == 3288

# ╔═╡ 4afd939a-378d-11eb-1609-45e4f3fcd2f7
begin
	items = [(1,2,3),(4,5,6)]
	zitems = zip(items...)
	fun(a,b,c) =  a*b+c
	fun.(zitems...)
end

# ╔═╡ Cell order:
# ╠═3e0c149a-3350-11eb-0b17-7bea5d1a28d4
# ╟─18080fe2-339b-11eb-34f0-f393638506fc
# ╠═4e084a4e-339b-11eb-3ffc-41118215f226
# ╠═24e1613c-339b-11eb-2976-59884053d1ee
# ╠═df1d1be2-36c6-11eb-07c6-493a612dd182
# ╠═6cc1348c-3786-11eb-07c8-bb4f6ddfe769
# ╠═cdea65a0-3786-11eb-2444-b7d482414e20
# ╠═d5b57b6c-3786-11eb-13da-33cee18cfa8b
# ╟─5f797e22-339d-11eb-1ff3-dda90f946838
# ╠═e3177eae-339e-11eb-3f9f-1f7230976832
# ╠═5218903e-3789-11eb-11c2-ebed36baf6fe
# ╠═6a2ab5da-3789-11eb-221f-7d829e2967b6
# ╠═a8b55df0-3789-11eb-28c1-edc492b7060e
# ╠═4afd939a-378d-11eb-1609-45e4f3fcd2f7

max_stars_repo_name: RensDimmendaal/advent-of-code-julia

max_issues_repo_name: RensDimmendaal/advent-of-code-julia

max_forks_repo_name: RensDimmendaal/advent-of-code-julia


text:
abstract SVM{TSpec <: SVMSpec} <: RegressionModel

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

svmModel(spec::SVMSpec, solution::PrimalSolution, predmodel::Predictor, X::AbstractMatrix, Y::AbstractArray) =
  primalSVMModel(spec, solution, predmodel, X, Y)

svmModel(spec::SVMSpec, solution::DualSolution, predmodel::Predictor, X::AbstractMatrix, Y::AbstractArray) =
  dualSVMModel(spec, solution, predmodel, X, Y)

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

"""
`PrimalSVM <: SVM`

Description
============

The result to the primal problem of a support vector machine.
It is the return value of the `svm` function if the parameter `dual = false`

Fields
=======

- **`params(this)`** : The structure and parameters of the utilized support vector machine

- **`nsv(this)`** : The number of support vectors

- **`svindex(this)`** : The indicies (i) of the support vectors within the training set

- **`features(this)`** : The full training set observations

- **`targets(this)`** : The full training set targets

- **`predmodel(this)`** : The underlying prediction model. It defines if the intercept
is present, as well as the input and output dimensions (univariate vs multivariate prediction)

- **`iterations(this)`** : Number of iterations used to fit the model.

- **`isconverged(this)`** : true, if the algorithm converged during training.

- **`fval(this)`** : The final objective value achieved during training.

- **`coef(this)`** : The fitted coefficients w⃗.

- **`details(this)`** : The training information returned by the solver itself.
It stores the final objective value, the number of used iterations,
the fitted coefficients w⃗, and a boolean that states if the algorithm was able to converge.

Methods
========

- **`predict(this, X)`** : predicts the response for the given observations in `X`.
Note that `X` should be of the same type as the data used for training the SVM.

- **`classify(this, X)`** : predicts the response for the given observations in `X`
and performs the decision function on the result.
Note that `X` should be of the same type as the data used for training the SVM.

- **`accuracy(this, X, y)`** : computes the accuracy by calculating the predictions
of `X` and comparing the results to `y`.
Note that `X` and `y` should have the same number of observations.

See also
=========

`svm`, `DualSVM`, `CSVM`

"""
abstract PrimalSVM{TSpec<:SVMSpec} <: SVM{TSpec}

"""
`DensePrimalSVM <: PrimalSVM`

See `PrimalSVM`
"""
type DensePrimalSVM{TSpec<:SVMSpec, TDetails<:PrimalSolution, TPred<:Predictor, XT<:DenseMatrix, YT<:AbstractVector} <: PrimalSVM{TSpec}
  params::TSpec
  details::TDetails
  predmodel::TPred
  nsv::Int
  svindex::Vector{Int}
  Xtrain::XT
  Ytrain::YT
end

"""
`SparsePrimalSVM <: PrimalSVM`

See `PrimalSVM`
"""
type SparsePrimalSVM{TSpec<:SVMSpec, TDetails<:PrimalSolution, TPred<:Predictor, XT<:AbstractSparseMatrix, YT<:AbstractVector} <: PrimalSVM{TSpec}
  params::TSpec
  details::TDetails
  predmodel::TPred
  nsv::Int
  svindex::Vector{Int}
  Xtrain::XT
  Ytrain::YT
end

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

"""
`DualSVM <: SVM`

Description
============

The result to the dual problem of a support vector machine.
It is the return value of the `svm` function if the parameter `dual = true`

Fields
=======

- **`params(this)`** : The structure and parameters of the utilized support vector machine

- **`nsv(this)`** : The number of support vectors

- **`svindex(this)`** : The indicies (i) of the support vectors within the training set

- **`features(this)`** : The full training set observations to which the model was fit

- **`targets(this)`** : The full training set targets to which the model was fit

- **`predmodel(this)`** : The underlying prediction model. It defines if the intercept
is present, as well as the input and output dimensions (univariate vs multivariate prediction)

- **`iterations(this)`** : Number of iterations used to fit the model.

- **`isconverged(this)`** : true, if the algorithm converged during training.

- **`fval(this)`** : The final objective value achieved during training.

- **`coef(this)`** : The coefficients (α) of the support vectors

- **`xsv(this)`** : The support vectors (if dense), or the training set observations (if sparse)

- **`ysv(this)`** : The targets of the support vectors

- **`details(this)`** : The training information returned by the solver itself.
It includes the final objective value, the number of used iterations,
the full coefficient vector α, the bias if present, and a boolean that
states if the algorithm was able to converge.

Methods
========

- **`predict(this, X)`** : predicts the response for the given observations in `X`.
Note that `X` should be of the same form as the data used for training the SVM.

- **`classify(this, X)`** : predicts the response for the given observations in `X`
and performs the decision function on the result.
Note that `X` should be of the same type as the data used for training the SVM.

- **`accuracy(this, X, y)`** : computes the accuracy by calculating the predictions
of `X` and comparing the results to `y`.
Note that `X` and `y` should have the same number of observations.

See also
=========

`svm`, `PrimalSVM`, `CSVM`

"""
abstract DualSVM{TSpec<:SVMSpec} <: SVM{TSpec}

"""
`DenseDualSVM <: DualSVM`

See `DualSVM`
"""
type DenseDualSVM{TSpec<:SVMSpec, TPred<:Predictor, XT<:DenseMatrix, YT<:AbstractVector} <: DualSVM{TSpec}
  params::TSpec
  details::DualSolution
  predmodel::TPred
  nsv::Int
  alpha::Vector{Float64}
  svindex::Vector{Int}
  Xsv::Vector{ContiguousView{Float64,1,Array{Float64,2}}}
  Ysv::Vector{Float64}
  Xtrain::XT
  Ytrain::YT
end

"""
`SparseDualSVM <: DualSVM`

See `DualSVM`
"""
type SparseDualSVM{TSpec<:SVMSpec, TPred<:Predictor, XT<:AbstractSparseMatrix, YT<:AbstractVector} <: DualSVM{TSpec}
  params::TSpec
  details::DualSolution
  predmodel::TPred
  nsv::Int
  alpha::Vector{Float64}
  svindex::Vector{Int}
  Xsv::XT
  Ysv::Vector{Float64}
  Xtrain::XT
  Ytrain::YT
end

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

function primalSVMModel{TSpec<:SVMSpec, TPred<:Predictor, TReal<:Real}(
    params::TSpec,
    s::PrimalSolution,
    predmodel::TPred,
    Xtrain::StridedMatrix,
    Ytrain::AbstractVector{TReal})
  n = size(Xtrain, 2)
  p = value(predmodel, Xtrain, minimizer(s))
  l = zeros(n)
  @inbounds for i = 1:n
    lc = value(params.loss, Float64(Ytrain[i]), p[i] - sign(p[i]) * 0.0005)
    lr = value(params.loss, Float64(Ytrain[i]), p[i] + sign(p[i]) * 0.0005)
    l[i] = lc != 0 || lr != 0
  end
  svindex = find(l)
  nsv = length(svindex)
  DensePrimalSVM{TSpec, typeof(s), TPred, typeof(Xtrain), typeof(Ytrain)}(params, s, predmodel, nsv, svindex, Xtrain, Ytrain)
end

function primalSVMModel{TSpec<:SVMSpec, TPred<:Predictor, TReal<:Real}(
    params::TSpec,
    s::PrimalSolution,
    predmodel::TPred,
    Xtrain::AbstractSparseMatrix,
    Ytrain::AbstractVector{TReal})
  n = size(Xtrain, 2)
  w = minimizer(s)
  p = if typeof(predmodel) <: LinearPredictor{true}
    fill(w[end] * predmodel.bias, size(Ytrain))
  else
    zeros(size(Ytrain))
  end
  @inbounds for i = 1:n
    tstart = Xtrain.colptr[i]
    tstop  = Xtrain.colptr[i+1] - 1
    for j = tstart:tstop
      k = Xtrain.rowval[j]
      p[i] += Xtrain.nzval[j] * w[k]
    end
  end
  l = zeros(n)
  @inbounds for i = 1:n
    lc = value(params.loss, Float64(Ytrain[i]), p[i] - sign(p[i]) * 0.0005,)
    lr = value(params.loss, Float64(Ytrain[i]), p[i] + sign(p[i]) * 0.0005)
    l[i] = lc != 0 || lr != 0
  end
  svindex = find(l)
  nsv = length(svindex)
  SparsePrimalSVM{TSpec, typeof(s), TPred, typeof(Xtrain), typeof(Ytrain)}(params, s, predmodel, nsv, svindex, Xtrain, Ytrain)
end

function dualSVMModel{TSpec<:SVMSpec, TPred<:Predictor, TReal<:Real}(
    params::TSpec,
    s::DualSolution,
    p::TPred,
    Xtrain::StridedMatrix,
    Ytrain::AbstractVector{TReal})
  svindex = find(minimizer(s))
  alpha = minimizer(s)[svindex]
  nsv = length(alpha)
  Xsv = Array(ContiguousView{Float64,1,Array{Float64,2}}, nsv)
  @inbounds for i in 1:nsv
    Xsv[i] = view(Xtrain, :, svindex[i])
  end
  Ysv = Ytrain[svindex]
  DenseDualSVM{TSpec, TPred, typeof(Xtrain), typeof(Ytrain)}(params, s, p, nsv, alpha, svindex, Xsv, Ysv, Xtrain, Ytrain)
end

function dualSVMModel{TSpec<:SVMSpec, TPred<:Predictor, TReal<:Real}(
    params::TSpec,
    s::DualSolution,
    p::TPred,
    Xtrain::AbstractSparseMatrix,
    Ytrain::AbstractVector{TReal})
  svindex = find(minimizer(s))
  alpha = minimizer(s)[svindex]
  nsv = length(alpha)
  Xsv = Xtrain
  Ysv = Ytrain[svindex]
  SparseDualSVM{TSpec, TPred, typeof(Xtrain), typeof(Ytrain)}(params, s, p, nsv, alpha, svindex, Xsv, Ysv, Xtrain, Ytrain)
end

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

@inline labels{TSpec<:SVCSpec}(svm::SVM{TSpec}) = [-1., 1]
@inline nobs(fit::SVM) = length(fit.Ytrain)
@inline features(fit::SVM) = fit.Xtrain
@inline targets(fit::SVM) = fit.Ytrain
@inline model_response(fit::SVM) = fit.Ytrain

@inline details(fit::SVM) = fit.details
@inline isconverged(fit::SVM) = isconverged(details(fit))
@inline iterations(fit::SVM) = iterations(details(fit))
@inline params(fit::SVM) = fit.params
@inline minimum(fit::SVM) = minimum(details(fit))
@inline minimizer(fit::SVM) = minimizer(details(fit))
@inline nsv(fit::SVM) = fit.nsv
@inline svindex(fit::SVM) = fit.svindex
@inline intercept(fit::SVM) = typeof(predmodel(fit)) <: LinearPredictor{true}
@inline predmodel(fit::SVM) = fit.predmodel
@inline coef(fit::PrimalSVM) = coef(details(fit))
@inline coef(fit::DualSVM) = fit.alpha
@inline bias(fit::DualSVM) = bias(details(fit))
@inline xsv(fit::DualSVM) = fit.Xsv
@inline ysv(fit::DualSVM) = fit.Ysv

@inline predict(fit::SVM) = predict(fit, features(fit))
@inline classify{TSpec<:SVCSpec}(svm::SVM{TSpec}) = classify(svm, features(svm))
@inline accuracy{TSpec<:SVCSpec}(svm::SVM{TSpec}) = accuracy(svm, features(svm), targets(svm))

function classify{TSpec<:SVCSpec}(svm::SVM{TSpec}, X)
  ŷ = predict(svm, X)
  t = ndims(ŷ) == 1 ? sign(ŷ) : vec(mapslices(indmax, ŷ, 1))
  t
end

function accuracy{TSpec<:SVCSpec}(svm::SVM{TSpec}, X, y)
  n = size(X,2)
  n == length(y) || throw(DimensionMismatch("X and y have to have the same number of observations"))
  ȳ = classify(svm, X)
  countnz(ȳ .== y) / n
end

function predict(fit::PrimalSVM, X::DenseMatrix)
  p = value(predmodel(fit), X, minimizer(details(fit)))
  size(p, 1) == 1 ? vec(p) : p
end

function predict(fit::PrimalSVM, X::AbstractSparseMatrix)
  n = size(X,2)
  w = minimizer(details(fit))
  p = if typeof(predmodel(fit)) <: LinearPredictor{true}
    fill(w[end] * predmodel(fit).bias, size(targets(fit)))
  else
    zeros(size(targets(fit)))
  end
  @inbounds for i = 1:n
    tstart = X.colptr[i]
    tstop  = X.colptr[i+1] - 1
    for j = tstart:tstop
      k = X.rowval[j]
      p[i] += X.nzval[j] * w[k]
    end
  end
  p
end

function predict(fit::DenseDualSVM, X::AbstractMatrix)
  n = size(X,2)
  result = zeros(n)
  @inbounds for i in 1:n
    for j in 1:nsv(fit)
      result[i] += coef(fit)[j] * ysv(fit)[j] * dot(xsv(fit)[j], view(X,:,i))
    end
    result[i] += bias(fit)
  end
  result
end

function predict(fit::SparseDualSVM, X::AbstractMatrix)
  n = size(X,2)
  result = zeros(n)
  tmp = 0.
  @inbounds for i in 1:n
    for j in 1:nsv(fit)
      tstart = xsv(fit).colptr[svindex(fit)[j]]
      tstop  = xsv(fit).colptr[svindex(fit)[j]+1] - 1
      tmp = 0.
      for k = tstart:tstop
        tmp += X[xsv(fit).rowval[k],i] * xsv(fit).nzval[k]
      end
      tmp *= ysv(fit)[j]
      tmp *= coef(fit)[j]
      result[i] += tmp
    end
    result[i] += bias(fit)
  end
  result
end

# ==========================================================================
# convert

function convert{TSpec <: SVMSpec{ScalarProductKernel{Float64}}}(
    ::Type{PrimalSVM},
    dual::DualSVM{TSpec})
  l = length(xsv(dual))
  d = length(xsv(dual)[1])
  # w = ∑ yᵢαᵢxᵢ
  w = zeros(d+1)
  for j = 1:d
    for i = 1:l
      @inbounds w[j] += ysv(dual)[i] * coef(dual)[i] * xsv(dual)[i][j]
    end
  end
  sol = if typeof(predmodel(dual)) <: LinearPredictor{true}
    w[end] = bias(dual)
    PrimalSolution(w, minimum(dual), iterations(dual), isconverged(dual))
  else
    PrimalSolution(w[1:d], minimum(dual), iterations(dual), isconverged(dual))
  end
  svmModel(params(dual), sol, predmodel(dual), features(dual), targets(dual))
end

function convert{TSpec <: SVMSpec{ScalarProductKernel{Float64}}}(
    ::Type{DualSVM},
    primal::PrimalSVM{TSpec})
  k = size(features(primal), 1)
  n = size(features(primal), 2)
  nsv = nsv(primal)
  svindex = svindex(primal)
  nsv <= k || throw(DimensionMismatch("Converting to dual solution is only possible if there are less (or equal) support vectors than there are features"))
  Q = zeros(k, nsv)
  for i = 1:nsv
    for j = 1:k
      @inbounds Q[j, i] = targets(primal)[svindex[i]] * features(primal)[j, svindex[i]]
    end
  end
  # DOESN'T WORK FOR NO BIAS
  w = minimizer(details(primal))[1:k]
  α = Q \ w
  alpha = zeros(n)
  alpha[svindex] = α
  sol = DualSolution(alpha, minimizer(details(primal))[end], minimum(primal), iterations(primal), isconverged(primal))
  svmModel(params(primal), sol, predmodel(primal), features(primal), targets(primal))
end

# ==========================================================================
# Plotting

function scatterplot{TSpec<:SVCSpec}(
    fit::PrimalSVM{TSpec};
    title::AbstractString = "Primal SVM Classification Plot",
    xlim = [0.,0.],
    ylim = [0.,0.],
    lbl = map(string, labels(fit)),
    nargs...)
  size(features(fit),1) == 2 || throw(DimensionMismatch("Can only plot the SVM classification for a two-dimensional featurespace (i.e. size(X,1) == 2)"))
  intercept_fit = typeof(predmodel(fit)) <: LinearPredictor{true}
  offset = intercept_fit ? -(minimizer(details(fit))[3] * predmodel(fit).bias) / minimizer(details(fit))[2] : 0.
  slope = -minimizer(details(fit))[1] / minimizer(details(fit))[2]
  x1 = vec(view(features(fit), 1, :))
  x2 = vec(view(features(fit), 2, :))
  x1sv = x1[svindex(fit)]
  x2sv = x2[svindex(fit)]
  xmin = minimum(x1); xmax = maximum(x1)
  ymin = minimum(x2); ymax = maximum(x2)
  xlim = xlim == [0.,0.] ? [xmin, xmax] : xlim
  ylim = ylim == [0.,0.] ? [ymin, ymax] : ylim
  notalphaindex = setdiff(1:size(features(fit),2), svindex(fit))
  x1 = x1[notalphaindex]
  x2 = x2[notalphaindex]
  y = targets(fit)[notalphaindex]
  fig = scatterplot(x1[y.<0], x2[y.<0]; title = title, xlim = xlim, ylim = ylim, name = lbl[1], nargs...)
  scatterplot!(fig, x1[y.>0], x2[y.>0], name = lbl[2])
  scatterplot!(fig, x1sv, x2sv, color = :yellow, name = "support vectors")
  lineplot!(fig, offset, slope, color = :white)
  xlabel!(fig, "X₁")
  ylabel!(fig, "X₂")
  fig
end

function scatterplot{TSpec<:SVCSpec}(
    fit::DualSVM{TSpec};
    title::AbstractString = "Dual SVM Classification Plot",
    xlim = [0.,0.],
    ylim = [0.,0.],
    lbl = map(string, labels(fit)),
    nargs...)
  size(features(fit),1) == 2 || throw(DimensionMismatch("Can only plot the SVM classification for a two-dimensional featurespace (i.e. size(X,1) == 2)"))
  x1 = vec(view(features(fit), 1, :))
  x2 = vec(view(features(fit), 2, :))
  x1sv = x1[svindex(fit)]
  x2sv = x2[svindex(fit)]
  xmin = minimum(x1); xmax = maximum(x1)
  ymin = minimum(x2); ymax = maximum(x2)
  xlim = xlim == [0.,0.] ? [xmin, xmax] : xlim
  ylim = ylim == [0.,0.] ? [ymin, ymax] : ylim
  notalphaindex = setdiff(1:size(features(fit),2), svindex(fit))
  x1 = x1[notalphaindex]
  x2 = x2[notalphaindex]
  y = targets(fit)[notalphaindex]
  fig = scatterplot(x1[y.<0], x2[y.<0]; title = title, xlim = xlim, ylim = ylim, name = lbl[1], nargs...)
  scatterplot!(fig, x1[y.>0], x2[y.>0], name = lbl[2])
  scatterplot!(fig, x1sv, x2sv, color = :yellow, name = "support vectors")
  xlabel!(fig, "X₁")
  ylabel!(fig, "X₂")
  fig
end

# ==========================================================================
# Base.show

function _showprimal(io::IO, fit)
  _printconverged(io, isconverged(fit), iterations(fit))
  _printvariable(io, 19, "details()", typeof(details(fit)))
  _printvariable(io, 19, "isconverged()", isconverged(fit))
  _printvariable(io, 19, "iterations()", iterations(fit))
  println(io, "\n  ◦  support vector machine:")
  _printvariable(io, 14, "params()", params(fit))
  println(io, "\n  ◦  objective value (f):")
  _printvariable(io, 17, "minimum()", minimum(fit))
  println(io, "\n  ◦  fitted coefficients (w⃗):")
  _printvariable(io, 17, "coef()", coef(fit))
  _printvariable(io, 17, "intercept()", intercept(fit))
  _printvariable(io, 17, "predmodel()", predmodel(fit))
  println(io, "\n  ◦  support vectors (estimated):")
  _printvariable(io, 17, "nsv()", nsv(fit))
  _printvariable(io, 17, "svindex()", svindex(fit))
  if size(features(fit),1) == 2 && size(features(fit),2) < 500
    println(io, "\n  ◦  classification plot (UnicodePlots.scatterplot(..)):")
    fig = scatterplot(fit, margin = 5, width = 30, height = 10, title = "")
    print(io, fig)
  end
end

function _showdual(io::IO, fit)
  _printconverged(io, isconverged(fit), iterations(fit))
  _printvariable(io, 19, "details()", typeof(details(fit)))
  _printvariable(io, 19, "isconverged()", isconverged(fit))
  _printvariable(io, 19, "iterations()", iterations(fit))
  println(io, "\n  ◦  support vector machine:")
  _printvariable(io, 14, "params()", params(fit))
  println(io, "\n  ◦  objective value (f):")
  _printvariable(io, 17, "minimum()", minimum(fit))
  println(io, "\n  ◦  fitted coefficients (α):")
  _printvariable(io, 17, "coef()", coef(fit))
  intercept(fit) && _printvariable(io, 17, "bias()", bias(fit))
  _printvariable(io, 17, "intercept()", intercept(fit))
  _printvariable(io, 17, "predmodel()", predmodel(fit))
  println(io, "\n  ◦  support vectors:")
  _printvariable(io, 17, "nsv()", nsv(fit))
  _printvariable(io, 17, "svindex()", svindex(fit))
  _printvariable(io, 17, "xsv()", typeof(xsv(fit)))
  _printvariable(io, 17, "ysv()", ysv(fit))
  if size(features(fit),1) == 2 && size(features(fit),2) < 500
    println(io, "\n  ◦  classification plot (UnicodePlots.scatterplot(..)):")
    fig = scatterplot(fit, margin = 5, width = 30, height = 10, title = "")
    print(io, fig)
  end
end

function show(io::IO, fit::PrimalSVM)
  println(io, typeof(fit), "\n")
  _showprimal(io, fit)
end

function show(io::IO, fit::DualSVM)
  println(io, typeof(fit), "\n")
  _showdual(io, fit)
end


text:
function makeGridTest()
	@testset "Make grid" begin
		T = 3;
		nCourseV = [2, 3];
		gridV = make_ncompleted_grid(T, nCourseV);
		@test gridV[1] == [0]
		@test gridV[2] == nCourseV
		@test gridV[3] == [4, 5, 6]

		T = 6;
		nCourseV = [5, 8, 10];
		gridV = make_ncompleted_grid(T, nCourseV);
		for t = 1 : T
			g = gridV[t];
			@test isequal(g, sort(g))
			@test all(g .>= (t-1) * nCourseV[1])
			@test all(g .<= (t-1) * nCourseV[end])
		end
    end
end

function updateTest()
	@testset "Update grid" begin
		T = 3;
		nCourseV = [2, 3];
		cg = make_course_grid(T, nCourseV);

		gridV = make_ncompleted_grid(T, nCourseV);
		@test n_completed_grid(cg, 3) == gridV[3]

		t = 2;
		nTaken = gridV[t][1];
		idx_nTakenV = collect(1 : 2);
		idxOutV = update_ncourses(cg, t, nTaken, idx_nTakenV);
		@test idxOutV == [1, 2]
		@inferred update_ncourses(cg, t, nTaken, idx_nTakenV)

		nTaken = gridV[t][2]
		idxOutV = update_ncourses(cg, t, nTaken, idx_nTakenV);
		@test idxOutV == [2, 3]

		idxOut2V = update_ncourses(cg, t, nTaken);
		@test idxOut2V[1:2] == idxOutV
	end
end

function indices_test()
	@testset "indices" begin
		T = 4;
		cg = make_test_course_grid(T);
		t = 2;
		gridV = cg.nGridV[t];
		idxV = ncompleted_to_indices(cg, t, gridV);
		@test idxV == collect(1 : length(gridV))

		idx2V = ncompleted_to_indices(cg, t, gridV .+ 100);
		@test all(idx2V .== 0)

		isValid = true;
		nV = collect(0 : 20);
		idx3V = ncompleted_to_indices(cg, t, nV);
		for (j, n) in enumerate(nV)
			if idx3V[j] > 0
				isValid = isValid && (gridV[idx3V[j]] == n);
			else
				isValid = isValid && !any(gridV .== n);
			end
		end
		@test isValid
	end
end

function singleton_grid_test()
	@testset "Singleton grid" begin
		T = 4;
		cg = make_course_grid(T, [3]);
		t = 2;
		@test n_tried_grid(cg, t) == [3]
		@test length(n_completed_grid(cg, t)) == 1
	end
end

@testset "CourseGrid" begin
	makeGridTest()
	updateTest()
	indices_test()
	singleton_grid_test()
end

# ------------


text:
""" Dynamic """

function dynamic(kf::KalmanFilter, x::AbstractVector, u::AbstractVector)
    return kf.A*x + kf.B*u + cholesky(kf.Q).L*randn(size(kf.Q, 1))
end

function prediction(kf::KalmanFilter, s::State, u::AbstractVector)
    x_hat = kf.A*s.x + kf.B*u # predicted state prior
    P_hat = kf.A*s.P*kf.A' + kf.Q # a priori state covariance
    return State(x_hat, P_hat)
end

function param_dynamic(θ, param_kf::ParamKalmanFilter, x::AbstractVector, u::AbstractVector)
    return param_kf.A(θ)*x + param_kf.B(θ)*u + cholesky(param_kf.Q(θ)).L*randn(size(param_kf.Q(θ), 1))
end

function param_prediction(θ, param_kf::ParamKalmanFilter, s::State, u::AbstractVector)
    x_hat = param_kf.A(θ)*s.x + param_kf.B(θ)*u # predicted state prior
    P_hat = param_kf.A(θ)*s.P*param_kf.A(θ)' + param_kf.Q(θ) # a priori state covariance
    return State(x_hat, P_hat)
end

""" Observation """

function observation(kf::KalmanFilter, x::AbstractVector)
    return kf.H*x + cholesky(kf.R).L*randn(size(kf.R, 1))
end

function param_observation(θ, param_kf::ParamKalmanFilter, x::AbstractVector)
    return param_kf.H(θ)*x + cholesky(param_kf.R(θ)).L*randn(size(param_kf.R(θ), 1))
end

function correction(kf::KalmanFilter, s::State, y::AbstractVector)
    v = y - kf.H*s.x # measurement pre fit residual
    S = kf.H*s.P*kf.H' + kf.R # pre fit residual covariance
    K = s.P*kf.H'*inv(S) # Kalman gain
    x_post = s.x + K*v # a posteriori state estiamate
    P_post = (I - K*kf.H)*s.P # a posteriori covariance estiamate
    return State(x_post, P_post)
end

function param_correction(θ, param_kf::ParamKalmanFilter, s::State, y::AbstractVector)
    v = y - param_kf.H(θ)*s.x # measurement pre fit residual
    S = param_kf.H(θ)*s.P*param_kf.H(θ)' + param_kf.R(θ) # pre fit residual covariance
    K = s.P*param_kf.H(θ)'*inv(S) # Kalman gain
    x_post = s.x + K*v # a posteriori state estiamate
    P_post = (I - K*param_kf.H(θ))*s.P # a posteriori covariance estiamate
    return State(x_post, P_post)
end

""" k-state MSE """

function state_mse(filter::KalmanFilter, s::State, u::AbstractVector, y::AbstractVector)
    ϵx = norm(s.x - (filter.A*s.x + filter.B*u))
    ϵy = norm(y - filter.H*s.x)
    return ϵx + ϵy
end


text:
using CompScienceMeshes
using BEAST

Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))
X = raviartthomas(Γ)

d = BEAST.CurlSingleLayerDP3D(1.0im, 1.0)
s = Helmholtz3D.singlelayer(wavenumber=1.0)

X = raviartthomas(Γ,BEAST.Continuity{:none})
Y = lagrangec0d1(Γ)

x = refspace(X)
y = refspace(Y)
charts = [chart(Γ,c) for c in cells(Γ)]

qs = BEAST.defaultquadstrat(s, x, x)
qd = quaddata(s, x, x, charts, charts, qs)
m1,m2 = 10,11
ql = quadrule(s, x, x, m1, charts[m1], m2, charts[m2], qd, qs)
@show @which quadrule(s, x, x, m1, charts[m1], m2, charts[m2], qd, qs)
@show typeof(ql)
BEAST.momintegrals!(s, x, x, charts[10], charts[20],
    zeros(ComplexF64,3,3), ql)

# sop = BEAST.singularpart(s)
ctr = CompScienceMeshes.center(charts[m1])
# BEAST.momintegrals!(sop, ctr, x, x, charts[m1], charts[m1], zeros(3,3), ql, 1.0)

dvg = divergence

@hilbertspace j p
@hilbertspace k q

κ = 1.0
curlA = strace(((x,y,z),)->-exp(-im*κ*z)*ŷ, Γ)
ndotA = dot(n, ((x,y,z),)->-x*exp(-im*κ*z)*ẑ)
a = @varform s[k,j] + s[dvg(k),p] +
    s[q,dvg(j)] - κ^2*s[q,p] == curlA[k]+ndotA[q]

Eq = @discretise a k∈X q∈Y j∈X p∈Y

u = solve(Eq)

uj = u[1:numfunctions(X)]
up = u[numfunctions(X)+1:end]

xrange = range(-2.0, stop=2.0, length=40)
# xrange = [0.0]
# yrange = [0.0]
zrange = range(-2.0, stop=2.0, length=40)
pts = [point(x,0,z) for x ∈ xrange, z ∈ zrange]

A1 = BEAST.DPVectorialTerm(1.0, 1.0*im)
A2 = BEAST.DPScalarTerm(1.0, 1.0*im)

p1 = potential(A1,pts,uj,X)
p2 = potential(A2,pts,up,Y)

p = p1 + p2


text:
import LinearAlgebra.rank

function rank(A::AbstractMatrix{T}) where T <: Number
  rank(qrfact(A, Val{true}))
end

function rank(QR::LinearAlgebra.QRPivoted{T, S} ) where S <: AbstractMatrix{T} where T <: Number
  R = QR[:R]
  sv = sort(abs.(diag(R)), rev = true)
  tol = eps(sv[1]) * maximum(size(A))
  n = length(sv)
  while n > 0 && sv[n] <= tol
    n -= 1
  end
  n
end

"""
Produce orthogonal Matrix, whose vectors are basis of the range of Matrix A.
If the matrix is surjective, the identity matrix is returned.
"""
function span(A::AbstractMatrix{T}) where T <: Number

  QR = qrfact(A, Val{true})
  n, m = size(A)
  r = rank(QR)
  if r < n
    QR.Q[:,1:r]
  else
    Matrix{T}(I, n, n)
  end
end

"""
Produce othogonal Matrix, whose vectors are basis of the kernel of Matrix A.
If the matrix is injective, the identity matrix is returned.
"""
function kern(A::AbstractMatrix{T}) where T <: Number

  QR = qrfact(A, Val{true})
  r = rank(QR)
  n, m = size(A)
  K = zeros(T, m, m - r)
  for k = r+1:m
    K[k,k-r] = one(T)
  end
  if r < m
    r1r = 1:r
    rrm = r+1:m
    R = QR[:R]
    RV = view(R, r1r, r1r)
    K[r1r,:] = -RV \ R[r1r, rrm]
    p = sortperm(QR[:p]) # inverse of permutation
    K, R = qr(K[p,:])
  end
  K
end

max_stars_repo_name: KlausC/LinAlgBigFloat

max_issues_repo_name: KlausC/LinAlgBigFloat

max_forks_repo_name: KlausC/LinAlgBigFloat


text:
@testset "Element stiffness matrix" begin

    # 1x1 linear quadrilateral cell
    # Isoparameteric
    # 2nd order quadrature rule

    E = 1
    nu = 0.3
    k = [
        1 / 2 - nu / 6,
        1 / 8 + nu / 8,
        -1 / 4 - nu / 12,
        -1 / 8 + 3 * nu / 8,
        -1 / 4 + nu / 12,
        -1 / 8 - nu / 8,
        nu / 6,
        1 / 8 - 3 * nu / 8,
    ]
    Ke =
        E / (1 - nu^2) * [
            k[1] k[2] k[3] k[4] k[5] k[6] k[7] k[8]
            k[2] k[1] k[8] k[7] k[6] k[5] k[4] k[3]
            k[3] k[8] k[1] k[6] k[7] k[4] k[5] k[2]
            k[4] k[7] k[6] k[1] k[8] k[3] k[2] k[5]
            k[5] k[6] k[7] k[8] k[1] k[2] k[3] k[4]
            k[6] k[5] k[4] k[3] k[2] k[1] k[8] k[7]
            k[7] k[4] k[5] k[2] k[3] k[8] k[1] k[6]
            k[8] k[3] k[2] k[5] k[4] k[7] k[6] k[1]
        ]

    problem = HalfMBB((2, 2), (1.0, 1.0), 1.0, 0.3, 1.0)

    #@test ElementFEAInfo(problem).Kes[1] ≈ Ke
end


text:
@testset "bpe" begin
    infile = joinpath(dirname(@__FILE__), "data/corpus.en")
    bpefile = joinpath(dirname(@__FILE__), "data/bpe.out")
    refile = joinpath(dirname(@__FILE__), "data/corpus.bpe.ref.en")

    bpe = Bpe(bpefile; sepsym="@@", endsym="")

    open(refile) do fr
        open(infile) do fi
            for (li, lr) ∈ zip(eachline(fi), eachline(fr))
                @test process_line(bpe, li) == lr
            end
        end
    end

    @test process_line(bpe, "  iron cement  \n") == "  ir@@ on c@@ ement  \n"
    @test process_line(bpe, "iron" * "\ua0" * "cement\n") == "ir@@ on@@ "*"\ua0"*"@@ c@@ ement\n"
    @test process_line(bpe, "\n") == "\n"
end


text:
immutable SBM{T<:RS} <: AbstractDEA{Convex,T}
  Data::DEAData

  function SBM{T}(X,Y) where T<:RS
    new(DEAData(X,Y))
  end
end

function getdata(DMU::SBM)
  return DMU.Data
end

function (DMU::SBM{T})(Xk::Array,Yk::Array; sx::Array = [],sy::Array = [],lambda::Array = []) where T<:RS
  K = getnrdmu(DMU.Data)
  N,M = getiodim(DMU.Data)

  # Set the appropriate RTS constraint depending on T
  RSA,RSb,RSsense = RSconstraint(K,T())

  f = zeros(K+1+N+M)
  f[1] = 1
  f[K+2:K+N+1] = -1./(N.*reshape(Xk,:,1))
  l = zeros(K+1+N+M)
  l[1] = eps()
  u = [Inf for i=1:K+1+N+M]

  sense = Array{Char}(N+M+1)
  sense[1:N+M+1] = '='

  X,Y = DMU[1:end]

  A = [reshape(Xk',:,1) -X' -eye(N) zeros(N,M); reshape(Yk',:,1) -Y' zeros(M,N) eye(M); 1 zeros(1,K+N) reshape(1./(M.*Yk),1,:); RSA zeros(1,N+M)]
  b = [zeros(N); zeros(M); 1; RSb]
  b = b[:]

  # Solve linear program
  sol = linprog(f,A,[sense;RSsense],b,l,u,ClpSolver())

  if sol.status == :Optimal
    rho = sol.objval
    t = sol.sol[1]
    lambda = sol.sol[2:K+1]./t
    sx = sol.sol[K+2:K+1+N]./t
    sy = sol.sol[K+1+N+1:K+1+N+M]./t
    res = DEAResult(rho,sx,sy,lambda)
  else
    res = DEAResult(NaN)
    println("Error: solution status $(sol.status)")
  end

	return res
end

function (DMU::SBM{T})() where T<:RS
  	Data = getdata(DMU)
  	res = Array{DEAResult}(getnrdmu(Data))

  	@sync @parallel for k in eachindex(Data)
  		res[k] = DMU(Data[k]...)
  	end

  	return res
end

max_stars_repo_name: kepiej/NEPA


text:
using TypeDBORM
using Documenter

DocMeta.setdocmeta!(TypeDBORM, :DocTestSetup, :(using TypeDBORM); recursive=true)

makedocs(;
    modules=[TypeDBORM],
    authors="Frank Urbach",
    repo="https://github.com/Humans-of-Julia/TypeDBORM.jl/blob/{commit}{path}#{line}",
    sitename="TypeDBORM.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Humans-of-Julia.github.io/TypeDBORM.jl/dev",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Guide" => "guide.md",
        "Structural Infos" => "structural_infos.md",
    ],
)

deploydocs(;
    repo="GitHub.com/Humans-of-Julia/TypeDBORM.jl.git",
)

max_issues_repo_name: Humans-of-Julia/TypeDBORM

max_forks_repo_name: Humans-of-Julia/TypeDBORM


text:
@testset "`majority`" begin
    @test majority([1, 2, 1, 3, 2, 2, 3], 1:3) == 2
    @test majority([1, 2, 1, 3, 2, 2, 3, 4], 3:4) == 3
    rng = StableRNG(42)
    picked = [majority(rng, 1:2, 1:2) for _ in 1:1_000_000]
    @test isapprox(count(==(1), picked), 500_000; atol=1000)
end

@testset "`Lighthouse.area_under_curve`" begin
    @test_throws AssertionError Lighthouse.area_under_curve([0, 1, 2], [0, 1])
    @test isapprox(Lighthouse.area_under_curve(collect(0:0.01:1), collect(0:0.01:1)), 0.5;
                   atol=0.01)
    @test isapprox(Lighthouse.area_under_curve(collect(0:0.01:(2π)), sin.(0:0.01:(2π))),
                   0.0; atol=0.01)
end

@testset "`Lighthouse.area_under_curve_unit_square`" begin
    @test_throws AssertionError Lighthouse.area_under_curve_unit_square([0, 1, 2], [0, 1])
    @test isapprox(Lighthouse.area_under_curve_unit_square(collect(0:0.01:1),
                                                           collect(0:0.01:1)), 0.5;
                   atol=0.01)
    @test isapprox(Lighthouse.area_under_curve_unit_square(collect(0:0.01:(2π)),
                                                           sin.(0:0.01:(2π))), 0.459;
                   atol=0.01)
end


text:
#####
# !!! intentional type piracy
# https://github.com/oolong-dev/OpenTelemetry.jl/issues/32
#####

function Base.schedule(t::Task)
    ctx = current_context()
    if ctx !== OpenTelemetryAPI.Context()
        if isnothing(t.storage)
            t.storage = IdDict()
        end
        t.storage[OpenTelemetryAPI.CONTEXT_KEY] = ctx
    end
    Base.enq_work(t)
end

#####
# GarishPrint
#####

function GarishPrint.pprint_struct(
    io::GarishPrint.GarishIO,
    mime::MIME"text/plain",
    @nospecialize(x),
    fields,
)
    t = typeof(x)
    isstructtype(t) || throw(ArgumentError("expect a struct type, got $t"))
    GarishPrint.pprint_struct_type(io, mime, t)
    print(io.bland_io, "(")

    nf = nfields(x)::Int
    nf == 0 && return print(io.bland_io, ")")

    io.compact || println(io.bland_io)

    # findout fields to print
    fields_to_print = Int[]
    for i in 1:nf
        f = fieldname(t, i)
        value = getfield(x, i)
        if !io.include_defaults &&
           GarishPrint.is_option(x) &&
           value == GarishPrint.field_default(t, f)
        elseif f in fields  # !!! only change
            push!(fields_to_print, i)
        end
    end

    GarishPrint.within_nextlevel(io) do
        for i in fields_to_print
            f = fieldname(t, i)
            value = getfield(x, i)
            GarishPrint.print_indent(io)
            GarishPrint.print_token(io, :fieldname, f)
            if io.compact
                print(io.bland_io, "=")
            else
                print(io.bland_io, " = ")
                io.state.offset = 3 + length(string(f))
            end

            if !isdefined(x, f) # print undef as comment color
                pprint(io, undef)
            else
                new_io = GarishPrint.GarishIO(io; limit = true)
                GarishPrint.pprint_field(new_io, mime, value)
            end

            if !io.compact || i != last(fields_to_print)
                print(io.bland_io, ", ")
            end

            if i != last(fields_to_print)
                io.compact || println(io.bland_io)
            end
        end
    end
    io.compact || println(io.bland_io)
    GarishPrint.print_indent(io)
    print(io.bland_io, ")")
end

function GarishPrint.pprint_struct(
    io::GarishPrint.GarishIO,
    mime::MIME"text/plain",
    a::AbstractInstrument,
)
    GarishPrint.pprint_struct(io, mime, a, (:name, :description, :unit))
end

function GarishPrint.pprint_struct(
    io::GarishPrint.GarishIO,
    mime::MIME"text/plain",
    a::AggregationStore,
)
    GarishPrint.pprint_struct(io, mime, a, (:unique_points, :n_max_points, :n_max_attrs))
end


text:
using EzXML
using DataStructures
using LightGraphs
using Vulkan_Headers_jll:vk_xml
using DataFrames
using Query

xdoc = readxml(vk_xml)
xroot = xdoc.root

include("spec/utils.jl")
include("spec/aliases.jl")
include("spec/dataframes.jl")
include("spec/graph.jl")
include("spec/structure_types.jl")


text:
using Test, LinearMaps, LinearAlgebra

@testset "wrapped maps" begin
    A = rand(10, 20)
    B = rand(ComplexF64, 10, 20)
    SA = Hermitian(A'A + I)
    SB = Hermitian(B'B + I)
    L = @inferred LinearMap{Float64}(A)
    @test summary(L) == "10×20 LinearMaps.WrappedMap{Float64}"
    @test occursin("10×20 LinearMaps.WrappedMap{Float64}", sprint((t, s) -> show(t, "text/plain", s), L))
    MA = @inferred LinearMap(SA, isposdef=true)
    MB = @inferred LinearMap(SB, isposdef=true)
    @test eltype(Matrix{Complex{Float32}}(LinearMap(A))) <: Complex
    @test size(L) == size(A)
    @test @inferred !issymmetric(L)
    @test @inferred issymmetric(MA)
    @test @inferred !issymmetric(MB)
    @test @inferred isposdef(MA)
    @test @inferred isposdef(MB)

    A = 2 * rand(ComplexF64, (20, 10)) .- 1
    v = rand(ComplexF64, 10)
    w = rand(ComplexF64, 20)
    u = rand(ComplexF64, 20, 1)
    V = rand(ComplexF64, 10, 3)
    W = rand(ComplexF64, 20, 3)
    Av = A * v
    AV = A * V
    M = @inferred LinearMap(A)
    N = @inferred LinearMap(M)

    @test @inferred M' * w == A' * w
    @test @inferred mul!(copy(V), adjoint(M), W) ≈ A' * W
    @test @inferred transpose(M) * w == transpose(A) * w
    @test @inferred transpose(LinearMap(transpose(M))) * v == A * v
    @test @inferred LinearMap(M')' * v == A * v
    @test @inferred(transpose(transpose(M))) === M
    @test @inferred(adjoint(adjoint(M))) === M
    Mherm = @inferred LinearMap(Hermitian(A'A), isposdef=true)
    @test @inferred ishermitian(Mherm)
    @test @inferred !issymmetric(Mherm)
    @test @inferred !issymmetric(transpose(Mherm))
    @test @inferred ishermitian(transpose(Mherm))
    @test @inferred ishermitian(Mherm')
    @test @inferred isposdef(Mherm)
    @test @inferred isposdef(transpose(Mherm))
    @test @inferred isposdef(adjoint(Mherm))
    @test @inferred !(transpose(M) == adjoint(M))
    @test @inferred !(adjoint(M) == transpose(M))
    @test @inferred transpose(M') * v ≈ transpose(A') * v
    @test @inferred transpose(LinearMap(M')) * v ≈ transpose(A') * v
    @test @inferred LinearMap(transpose(M))' * v ≈ transpose(A)' * v
    @test @inferred transpose(LinearMap(transpose(M))) * v ≈ Av
    @test @inferred adjoint(LinearMap(adjoint(M))) * v ≈ Av

    for w in (vec(u), u)
        @test @inferred mul!(copy(w), transpose(LinearMap(M')), v) ≈ transpose(A') * v
        @test @inferred mul!(copy(w), LinearMap(transpose(M))', v) ≈ transpose(A)' * v
        @test @inferred mul!(copy(w), transpose(LinearMap(transpose(M))), v) ≈ Av
        @test @inferred mul!(copy(w), adjoint(LinearMap(adjoint(M))), v) ≈ Av
    end
    @test @inferred mul!(copy(W), transpose(LinearMap(M')), V) ≈ transpose(A') * V
    @test @inferred mul!(copy(W), LinearMap(transpose(M))', V) ≈ transpose(A)' * V
    @test @inferred mul!(copy(W), transpose(LinearMap(transpose(M))), V) ≈ AV
    @test @inferred mul!(copy(W), adjoint(LinearMap(adjoint(M))), V) ≈ AV
    @test @inferred mul!(copy(V), transpose(M), W) ≈ transpose(A) * W
    @test @inferred mul!(copy(V), adjoint(M), W) ≈ A' * W

    B = @inferred LinearMap(Symmetric(rand(10, 10)))
    @test transpose(B) == B
    @test B == transpose(B)

    B = @inferred LinearMap(Hermitian(rand(ComplexF64, 10, 10)))
    @test adjoint(B) == B
    @test B == B'

    N = 10
    Id = LinearMap(identity, N; issymmetric=true)
    A = rand(N, N)
    B = similar(A)
    @test mul!(copy(B), Id, A) == A
    @test mul!(B, A, Id) == B == A
    @test mul!(copy(B), A, LinearMap(Matrix(Id))) == A
    @test mul!(B, Id, A, true, true) ≈ 2A
    @test mul!(B, A, Id, true, true) == B == 3A
    @test mul!(copy(B), A, LinearMap(Matrix(Id)), true, true) == 4A

    # wrapped vectors viewed as column matrices
    u = ones(1); U = @inferred LinearMap(u)
    @test U isa LinearMap{Float64}
    @test issymmetric(U)
    @test ishermitian(U)
    @test isposdef(U)
    U = @inferred LinearMap(u.*im)
    @test U isa LinearMap{ComplexF64}
    @test issymmetric(U)
    @test !ishermitian(U)
    @test !isposdef(U)
    # symmetry
    O4 = ones(4)
    O3 = ones(3)
    Z3 = zeros(3)
    for M in (Bidiagonal(O4, O3, :U),
              Bidiagonal(O4, Z3, :L),
              Diagonal(O4),
              Diagonal(im*O4),
              SymTridiagonal(O4, O3),
              Tridiagonal(O3, O4, O3),
              Tridiagonal(O3, O4, Z3),
              sparse(Tridiagonal(O3, O4, O3)),
              sparse(Diagonal(im*O4)),
            )
        @test issymmetric(LinearMap(M)) == issymmetric(M)
        @test ishermitian(LinearMap(M)) == ishermitian(M)
    end
end


text:
struct GeneralizedPoisson <: Distributions.DiscreteUnivariateDistribution
    θ::Float64
    λ::Float64
end

function Distributions.pdf(G::GeneralizedPoisson, x::Int64)
    G.θ*(G.θ+G.λ*x)^(x-1)*exp(-G.θ - G.λ*x)/factorial(x)
end

function Distributions.logpdf(G::GeneralizedPoisson, x::Int64)
    a = log(G.θ)
    a += (x-1)*log(G.θ+G.λ*x)
    a -= G.θ
    a -= G.λ*x
    a -= lgamma(x+1)
    a
end

Distributions.mean(G::GeneralizedPoisson) = G.θ/(1-G.λ)
Distributions.var(G::GeneralizedPoisson) = G.θ/(1-G.λ)^3


text:
macro S(gen)
    :(SVector{length($(esc(gen)))}($(esc(gen))))
end

macro animated(pfunc, k_range, filename)
    @assert(pfunc.head == :call, "Animated must be used with a functioncall
             that returns a Plots.plot and has a kwarg for the time step `k` to
             render.")

    return :(animate_plot($((esc(a) for a in pfunc.args)...);
                          k_range=$(esc(k_range)),
                          filename=$(esc(filename))))
end


text:
using Printf
using LinearAlgebra
using SpecialFunctions

include("constants.jl")
include("alias.jl")
include("Atoms.jl")
include("utils.jl")
include("PGBF.jl")
include("CGBF.jl")
include("overlap.jl")
include("kinetic.jl")
include("nuclear.jl")
include("coulomb.jl")
include("sto3g.jl")
include("BasisSet.jl")
include("RHF.jl")

function test_main(filename)
    atoms = init_atoms_xyz(filename)
    energy, E, U = RHF(atoms)
end

#test_main("H2_v2.xyz")
#test_main("H2O.xyz")
#test_main("LiH.xyz")
test_main("I2.xyz")