text:
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Active forces
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
export AbstractActiveForce
export ActiveBrownian
export RunAndTumble

abstract type AbstractActiveForce end

"""
Stores properties of active Brownian force

    ActiveBrownian(; γv, D_rot, align)

Initializes an active Brownian force with propulsion `γv`, rotational diffusion
`D_rot`.  If `align == true`, the active Brownian force will align with the
particle's orientation and allow us to introduce orientational interactions.
"""
mutable struct ActiveBrownian <: AbstractActiveForce
    γv_x::Float64
    γv_y::Float64
    θ::Float64

    γv::Float64
    amplitude::Float64

    align::Bool

    function ActiveBrownian(; γv::Float64, D_rot::Float64, align::Bool = false)
        θ = 2 * pi * rand()
        amplitude = sqrt(2 * D_rot)
        new(γv * cos(θ), γv * sin(θ), θ, γv, amplitude, align)
    end
end

"""
Stores properties of run-and-tumble force

    RunAndTumble(; γv, α, align)

Initializes a run-and-tumble force with propulsion `γv` and tumble rate `α`. If
`align == true`, the active Brownian force will align with the particle's
orientation and allow us to introduce orientational interactions.
"""
mutable struct RunAndTumble <: AbstractActiveForce
    γv_x::Float64
    γv_y::Float64
    θ::Float64
    τ_run::Float64

    γv::Float64
    α::Float64

    align::Bool

    function RunAndTumble(; γv::Float64, α::Float64, align::Bool = false)
        θ = 2 * pi * rand()
        τ_run = -log(rand()) / α

        new(γv * cos(θ), γv * sin(θ), θ, τ_run, γv, α, align)
    end
end

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Particle
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
export Particle

"""
Stores properties of a particle

    Particle(; ptype, x, y, θ, R, γ_trans, γ_rot, D_trans, D_rot, active_force)

# Arguments
- `ptype::Symbol = :particle`: particle type
- `x::Float64`: x position
- `y::Float64`: y position
- `θ::Float64 = 2 * pi * rand()`: in-plane orientation
- `R::Float64`: radius
- `γ_trans::Float64 = -1.0`: translational friction.  If `γ_trans < 0`, use `γ_trans = 6 * pi * 8.9e-4 * R`.
- `γ_rot::Float64 = -1.0`: rotational friction. If `γ_rot < 0`, use `γ_rot = 8 * pi * 8.9e-4 * R^3`.
- `D_trans::Float64`: translational diffusion
- `D_rot::Float64`: rotational diffusion
- `active_force::Union{AbstractActiveForce, Nothing} = nothing`: active force.
"""
mutable struct Particle
    ptype::Symbol # particle type

    x::Float64 # x position
    y::Float64 # y position
    θ::Float64 # orientation

    f_x::Float64 # x force
    f_y::Float64 # y force
    t_θ::Float64 # torque in z

    R::Float64 # particle radius
    γ_trans::Float64 # translational friction
    γ_rot::Float64 # rotational friction
    D_trans::Float64 # translational diffusivity
    D_rot::Float64 # rotational diffusivity

    active_force::Union{AbstractActiveForce, Nothing} # active force/self-propulsion

    function Particle(; ptype::Symbol = :particle,
                        x::Float64, y::Float64, θ::Float64 = 2 * pi * rand(),
                        R::Float64, γ_trans::Float64 = -1.0, γ_rot::Float64 = -1.0, D_trans::Float64, D_rot::Float64 = 0.0,
                        active_force::Union{AbstractActiveForce, Nothing} = nothing)
        if γ_trans < 0.0
            γ_trans = 6 * pi * 8.9e-4 * R
        end
        if γ_rot < 0.0
            γ_rot = 8 * pi * 8.9e-4 * R^3
        end

        new(ptype, x, y, θ, 0.0, 0.0, 0.0, R, γ_trans, γ_rot, D_trans, D_rot, active_force)
    end
end

"""
    Molecule(particles, pairs, triplets)

...
"""
struct Molecule
    particles::Array{Particle, 1}
    pairs::Array{NTuple{2, Particle}, 1}
    triplets::Array{NTuple{3, Particle}, 1}

    function Molecule(; particles::Array{Particle, 1} = Array{Particle, 1}(),
                        pairs::Array{NTuple{2, Particle}, 1} = Array{NTuple{2, Particle}, 1}(),
                        triplets::Array{NTuple{3, Particle}, 1} = Array{NTuple{3, Particle}, 1}())
        new(particles, pairs, triplets)
    end
end

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Cell list
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
export CellList

"""
Stores properties of a cell list

    CellList(; particles, L_x, L_y, cutoff)

Creates a cell list for `particles` in region/simulation with dimensions
`L_x`, `L_y`.  `cutoff` gives the approximate size of the cells, the actual size
of which is chosen so that an integer number of cells fit in each dimension.
"""
struct CellList
    start_pid::Array{Int64, 2}
    next_pid::Array{Int64, 1}

    particles::Array{Particle, 1}

    num_cells_x::Int64
    num_cells_y::Int64
    cell_spacing_x::Float64
    cell_spacing_y::Float64

    function CellList(; particles::Array{Particle, 1}, L_x::Float64, L_y::Float64, cutoff::Float64)
        num_cells_x = trunc(Int64, L_x / cutoff)
        num_cells_y = trunc(Int64, L_y / cutoff)
        cell_spacing_x = L_x / num_cells_x
        cell_spacing_y = L_y / num_cells_y

        start_pid = -ones(Int64, num_cells_x, num_cells_y)
        next_pid = -ones(Int64, length(particles))

        for (n, particle) in enumerate(particles)
            i = trunc(Int64, particle.x / cell_spacing_x) + 1
            j = trunc(Int64, particle.y / cell_spacing_y) + 1

            if start_pid[i, j] > 0
                next_pid[n] = start_pid[i, j]
            end
            start_pid[i, j] = n
        end

        new(start_pid, next_pid, particles, num_cells_x, num_cells_y, cell_spacing_x, cell_spacing_y)
    end
end

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Interactions
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
export AbstractInteraction
export LennardJones
export HarmonicBond
export HarmonicCosineAngle

abstract type AbstractInteraction end

"""
Stores properties of a Lennard-Jones interaction

    LennardJones(; particles, cell_list, ϵ, σ, cutoff, multithreaded, use_newton_3rd)
...
# Arguments
- `particles::Array{Particle, 1}`: particles to compute Lennard-Jones for
- `cell_list::CellList`: cell list for neighbors
- `ϵ::Float64`: energy scale of Lennard-Jones potential
- `σ::Float64 = -1.0`: length scale of Lennard-Jones potential.  If `σ < 0`, use diameters of interacting particles.
- `cutoff::Float64 = -1.0`: distance to truncate interaction.  If `cutoff < 0`, use `cutoff = 2^(1.0 / 6.0) * σ`, which corresponds to WCA interaction.
- `multithreaded::Bool = false`: if `true`, split particles between threads.
- `use_newton_3rd::Bool = false`: if `true`, use Newton's 3rd law to also compute force on neighbors.  Be CAREFUL with overcounting and race conditions with threads.
...
"""
struct LennardJones <: AbstractInteraction
    ϵ::Float64
    σ::Float64
    cutoff::Float64

    particles::Array{Particle, 1}
    cell_list::CellList

    multithreaded::Bool
    use_newton_3rd::Bool

    function LennardJones(; particles::Array{Particle, 1}, cell_list::CellList,
                            ϵ::Float64, σ::Float64 = -1.0, cutoff::Float64 = -1.0,
                            multithreaded::Bool = false, use_newton_3rd::Bool = false)
        new(ϵ, σ, cutoff, particles, cell_list, multithreaded, use_newton_3rd)
    end
end

"""
Stores properties of a harmonic bond between pair of particles

    HarmonicBond(; pairs, k_bond, l_rest, multithreaded)

Initialize a harmonic bond between `pairs` of particles with stiffness `k_bond`
and rest length `l_rest`.  For now `multithreaded` should only really be used
for dimers to avoid race conditions.  This will be changed with the
introduction of Molecules.
"""
struct HarmonicBond <: AbstractInteraction
    k_bond::Float64
    l_rest::Float64

    particle_pairs::Array{NTuple{2, Particle}, 1}

    multithreaded::Bool

    function HarmonicBond(; pairs::Array{NTuple{2, Particle}, 1}, k_bond::Float64, l_rest::Float64 = 0.0, multithreaded::Bool = false)
        new(k_bond, l_rest, pairs, multithreaded)
    end
end

"""
    etc

etc
"""
struct HarmonicAngle <: AbstractInteraction
    k_θ::Float64
    θ_rest::Float64

    particle_triplets::Array{NTuple{3, Particle}, 1}

    multithreaded::Bool

    function HarmonicAngle(; triplets::Array{NTuple{3, Particle}, 1}, k_θ::Float64, θ_rest::Float64, multithreaded::Bool = false)
        new(k_θ, θ_rest, triplets, multithreaded)
    end
end

"""
    ...

...
"""
struct CosineAngle <: AbstractInteraction

end

"""
Stores properties of bending force for triplets of particles

    HarmonicCosineAngle(triplets, k_cosθ, cosθ_rest, multithreaded)

Initialize a bending force for `triplets` of particles with
with bending stiffness `k_cosθ` with rest cosine angle
`cosθ_rest`.  For now `multithreaded` should only really be used
for trimers to avoid race conditions.  This will be changed with
the introduction of Molecules.
"""
struct HarmonicCosineAngle <: AbstractInteraction
    k_cosθ::Float64
    cosθ_rest::Float64

    particle_triplets::Array{NTuple{3, Particle}, 1}

    multithreaded::Bool

    function HarmonicCosineAngle(; triplets::Array{NTuple{3, Particle}, 1}, k_cosθ::Float64, cosθ_rest::Float64, multithreaded::Bool = false)
        new(k_cosθ, cosθ_rest, triplets, multithreaded)
    end
end

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# External forces
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
export AbstractExternalForce
export ConstantForce
export HarmonicTrap

abstract type AbstractExternalForce end

"""
Stores properties of constant force

    ConstantForce(particles; f_x, f_y)

Initialize a constant force `f_x`, `f_y` for 'particles'.
"""
struct ConstantForce <: AbstractExternalForce
    f_x::Float64
    f_y::Float64

    particles::Array{Particle, 1}

    multithreaded::Bool

    function ConstantForce(particles::Array{Particle, 1}; f_x::Float64, f_y::Float64, multithreaded::Bool = false)
        new(f_x, f_y, particles, multithreaded)
    end
end

"""
Stores properties of harmonic trap

    HarmonicTrap(particles; x_center, y_center, k_trap)

Initialize a harmonic trap with stiffness `k_trap` centered at `x_center`,
`y_center` for `particles`.
"""
struct HarmonicTrap <: AbstractExternalForce
    x_center::Float64
    y_center::Float64
    k_trap::Float64

    particles::Array{Particle, 1}

    multithreaded::Bool

    function HarmonicTrap(particles::Array{Particle, 1}; x_center::Float64, y_center::Float64, k_trap::Float64, multithreaded::Bool = false)
        new(x_center, y_center, k_trap, particles, multithreaded)
    end
end

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Integrators
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
export AbstractIntegrator
export Brownian

abstract type AbstractIntegrator end

"""
Stores properties of a Brownian integrator

    Brownian(; particles, dt, rotations, multithreaded)

Initialize a Brownian integrator for `particles` with timestep `dt`.  If
`rotations == true`, integrate the orientational degree of freedom.  If
`multithreaded == true`, split particles between threads.
"""
struct Brownian <: AbstractIntegrator
    dt::Float64
    particles::Array{Particle, 1}

    rotations::Bool

    multithreaded::Bool

    function Brownian(; particles::Array{Particle, 1}, dt::Float64, rotations::Bool = false, multithreaded::Bool = false)
        new(dt, particles, rotations, multithreaded)
    end
end

struct Langevin <: AbstractIntegrator
    # to be added
end

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Simulation data
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
export Simulation

"""
Stores all information needed to run simulation

    Simulation(...)
...
# Parameters
- `descriptor::String`: a description of the simulation for easy referencing
- `L_x::Float64`: width of simulation region
- `L_y::Float64`: height of simulation region
- `periodic_in_x::Bool = true`: if `true`, apply periodicity along x
- `periodic_in_y::Bool = true`: if `true`, apply periodicity along y
- `particles::Array{Particle, 1}`: all particles in simulation
- `cell_lists::Array{CellList, 1}`: all cell lists used
- `interactions::Array{AbstractInteraction, 1}`: all interactions used
- `external_forces::Array{AbstractExternalForce, 1}`: all external forces used
- `dt::Float64`: timestep
- `integrators::Array{AbstractIntegrator}`: all integrators used
- `num_steps::Int64`: total steps to run
- `save_interval::Int64`: store data periodically
- `save_particles::Array{Particle, 1}`: particles to store
- `history::Array{Array{Particle, 1}}`: history of particles
...
"""
mutable struct Simulation
    descriptor::String

    L_x::Float64
    L_y::Float64
    periodic_in_x::Bool
    periodic_in_y::Bool

    particles::Array{Particle, 1}
    bonds::Array{NTuple{2, Particle}, 1}
    angles::Array{NTuple{3, Particle}, 1}

    cell_lists::Array{CellList, 1}
    interactions::Array{AbstractInteraction, 1}
    external_forces::Array{AbstractExternalForce, 1}

    dt::Float64
    integrators::Array{AbstractIntegrator, 1}
    num_steps::Int64

    save_interval::Int64
    things_to_save::NamedTuple{(:particles, :bonds), Tuple{Array{Particle, 1}, Array{NTuple{2, Particle}, 1}}}
    history::Array{NamedTuple{(:particles, :bonds), Tuple{Array{Particle, 1}, Array{NTuple{2, Particle}, 1}}}, 1}

    function Simulation(; descriptor::String = "No description given...",
                          L_x::Float64 = 0.0,
                          L_y::Float64 = 0.0,
                          periodic_in_x::Bool = true,
                          periodic_in_y::Bool = true,
                          particles::Array{Particle, 1} = Particle[],
                          bonds::Array{NTuple{2, Particle}, 1} = NTuple{2, Particle}[],
                          angles::Array{NTuple{3, Particle}, 1} = NTuple{3, Particle}[],
                          cell_lists::Array{CellList, 1} = CellList[],
                          interactions::Array{AbstractInteraction, 1} = AbstractInteraction[],
                          external_forces::Array{AbstractExternalForce, 1} = AbstractExternalForce[],
                          dt::Float64 = 0.0,
                          integrators::Array{AbstractIntegrator, 1} = AbstractIntegrator[],
                          num_steps::Int64 = 0,
                          save_interval::Int64 = 0,
                          things_to_save::NamedTuple{(:particles, :bonds), Tuple{Array{Particle, 1}, Array{NTuple{2, Particle}, 1}}} = (particles = Particle[], bonds = NTuple{2, Particle}[]),
                          history::Array{NamedTuple{(:particles, :bonds), Tuple{Array{Particle, 1}, Array{NTuple{2, Particle}, 1}}}, 1} = NamedTuple{(:particles, :bonds), Tuple{Array{Particle, 1}, Array{NTuple{2, Particle}, 1}}}[])
        new(descriptor, L_x, L_y, periodic_in_x, periodic_in_y, particles, bonds, angles, cell_lists, interactions, external_forces, dt, integrators, num_steps, save_interval, things_to_save, history)
    end
end

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Useful functions and macros
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
export use_threads
export wrap_displacement
export wrap_positions
export print_message

"""
    @use_threads multithreaded ...

Applies `Threads.@threads` if `multithreaded == true`.  This is mostly used to
shorten code.
"""
macro use_threads(multithreaded::Union{Expr, Symbol}, expr::Expr)
    esc(quote
        if $multithreaded
            Threads.@threads $expr
        else
            $expr
        end
    end)
end

"""
    wrap_displacement(displacement; period)

Returns a new displacement after applying periodic boundary conditions.  The
periodicity is given by `period`.  If period < 0, then no periodicity is
applied.
"""
@inline function wrap_displacement(displacement::Float64; period::Float64)
    if period > 0.0 && abs(displacement) > period / 2
        return displacement - sign(displacement) * period
    end
    return displacement
end

"""
    wrap_position(position; period)

Returns a new position after applying periodic boundary conditions.  The
periodicity is given by `period`.  If `period < 0`, then no periodic boundary
condition is applied.
"""
@inline function wrap_position(position::Float64; period::Float64)
    if period > 0.0
        return mod(position, period)
    end
    return position
end

"""

"""
function print_message(message::String)
    message = string("++++++ ", message, " ++++++")
    N_space = max(div(62 - length(message), 2), 0)
    for _ = 1 : N_space
        message = string(" ", message)
    end

    println("")
    println(message)
    println("")
end


text:
module DBaseReader

# File format references used:
# [1] http://www.dbf2002.com/dbf-file-format.html
# [2] http://www.clicketyclick.dk/databases/xbase/format/dbf.html#DBF_NOTE_1_TARGET

using DataArrays

export readdbf, readdbfheader, DBFHeader

function bytestostring(ary)
  # findfirst plus sub-indexing is faster than rstrip(..., ['\0']])
  iend = (findfirst(ary, 0x00)-1)
  if iend > 0
    return rstrip(String(ary[1:iend]))
  else
    return rstrip(String(ary))
  end
end

mutable struct DBFFieldInfo
  name::Symbol
  kind::Char
  displacement::Int64
  bytesize::Int64
  decimalplaces::Int64
  rest::Array{UInt8}
end

# References [1] and [2] disagree about the lengths and meanings of bytes
# after offset 17. We just stick it all into the "rest" field.
function readfieldinfo(f)
  info = DBFFieldInfo(                        # Byte offset
    Symbol(bytestostring(read(f, 11))),  # 0-10
    Char(read(f, 1)[1]),                 # 11
    read(f, Int32),                           # 12-15
    read(f, Int8),                            # 16
    read(f, Int8),                            # 17
    read(f, 14))                         # 18-31
  # FoxPro 'F' is same as 'N'
  if info.kind == 'F'; info.kind = 'N' end
  return info
end

function store!(col, datum, typ, decimalplaces, name)
  try
    if typ == 'C'
      push!(col, bytestostring(datum))
    elseif typ == 'N'
      if decimalplaces == 0
        push!(col, parse(Int64, bytestostring(datum)))
      else
        push!(col, parse(Float64, bytestostring(datum)))
      end
    elseif typ == 'D' # Date in YYYYMMDD character format
      yyyy = parse(Int64, bytestostring(datum[1:4]))
      mm = parse(Int64, bytestostring(datum[5:6]))
      dd = parse(Int64, bytestostring(datum[7:8]))
      if yyyy == 0 || mm == 0 || dd == 0
        push!(col, NA)
      else
        push!(col, Date(yyyy, mm, dd))
      end
    elseif typ == 'L' # Character boolean
      if datum in "TtYy"
        push!(col, true)
      elseif datum in "FfNn"
        push!(col, false)
      else
        push!(col, NA)
      end
    else
      push!(col, NA)
    end
  catch # If a conversion failed, call it an NA and move on.
    warn("Failed to convert '$typ' value in column $name")
    push!(col, NA)
  end
  return nothing
end

function parserecords!(cols, chunk, info, recordlength)
  coloffset = 0
  # Iterate column-wise. Yields better performance than row-wise.
  for inf in info
    name = inf.name
    col = cols[name]
    istart = coloffset + 2
    records = length(chunk) / recordlength
    iend = 0
    for x in 1:records
      iend = istart + inf.bytesize - 1
      store!(col, chunk[istart:iend], inf.kind, inf.decimalplaces, name)
      istart += recordlength
    end
    coloffset += inf.bytesize
  end
  return nothing
end

function stubcolumns(info)
  cols = Dict()
  for inf in info
    if inf.kind == 'C'
      cols[inf.name] = DataArray(String[], Bool[])
    elseif inf.kind == 'D'
      cols[inf.name] = DataArray(Date[], Bool[])
    elseif inf.kind == 'L'
      cols[inf.name] = DataArray(Bool[], Bool[])
    elseif inf.kind == 'N' && inf.decimalplaces == 0
      cols[inf.name] = DataArray(Int64[], Bool[])
    elseif inf.kind == 'N' && inf.decimalplaces != 0
      cols[inf.name] = DataArray(Float64[], Bool[])
    else
      warn("Unsupported column type $(inf.kind); will fill with NA's")
      cols[inf.name] = DataArray(Bool[], Bool[])
    end
  end
  return cols
end

"""
`DBFHeader`

All the information you ever expected to get out of a dBase header.
Fields are: `dbasetype`, `yearmodified`, `monthmodified`, `daymodified`,
`numrecords`, `firstrecordoffset`, `recordlength`, `reservedspace`,
`tableflags`, `codepagemark`, `reservedzeroes`.
"""
mutable struct DBFHeader
  dbasetype::UInt8
  yearmodified::Int64
  monthmodified::Int64
  daymodified::Int64
  numrecords::Int64
  firstrecordoffset::Int64
  recordlength::Int64
  reservedspace::Array{UInt8}
  tableflags::UInt8
  codepagemark::UInt8
  reservedzeroes::Array{UInt8}
end
DBFHeader() = DBFHeader(0x00,0,0,0,0,0,0,UInt8[],0x00,0x00,UInt8[])

"""
`readdbfheader(path::AbstractString)`

Reads the header from a dBase file and returns a `DBFHeader` object.
`path` specifies the desired .dbf file, and must include the file extension.
"""
function readdbfheader(path::AbstractString)
  f = open(path)

  h = DBFHeader()
                                                      # Byte offset
  h.dbasetype = read(f, 1)[1]                         # 0
  h.yearmodified = 1900 + read(f, Int8)               # 1
  h.monthmodified = read(f, Int8)                     # 2
  h.daymodified = read(f, Int8)                       # 3
  h.numrecords = read(f, Int32)                       # 4-7
  h.firstrecordoffset = read(f, Int16)                # 8-9
  h.recordlength = read(f, Int16)                     # 10-11
  h.reservedspace = read(f, 16)                       # 12-27
  h.tableflags = read(f, 1)[1]                        # 28
  h.codepagemark = read(f, 1)[1]                      # 29
  h.reservedzeroes = read(f, 2)                       # 30-31

  close(f)

  return h
end

"""
`readdbf(path::AbstractString, maxbytes::Int64=10000000)`

Reads a dBase file into a `Dict` containing `Symbol`s as keys, and
`DataArray`s as values. `path` specifies the desired file, and must
include the file extension. `maxbytes` is optional, and specifies
the maximum number of bytes to read from the file at once.
"""
function readdbf(path::AbstractString, maxbytes::Int64=10000000)
  header = readdbfheader(path)

  f = open(path)

  # Skip over header
  skip(f, 32)

  info = []
  for x in 1:((header.firstrecordoffset - 33) // 32)
    push!(info, readfieldinfo(f))
  end
  read(f, 1) # eat the terminator

  cols = stubcolumns(info)

  # Add data column-wise rather than row-wise because the looping is
  # significantly faster. To avoid issues with streams that can't
  # (or shouldn't) be read all at once, read the stream in contiguous chunks,
  # stepping column-wise through each chunk.
  chunking = div(maxbytes, header.recordlength)
  chunks, extra = divrem(header.numrecords, chunking)

  for x in 1:chunks
    parserecords!(cols, read(f, chunking * header.recordlength), info,
                  header.recordlength)
  end
  extra != 0 && parserecords!(cols, read(f, extra * header.recordlength),
                              info, header.recordlength)

  d = Dict()

  for (name, col) in cols
    d[name] = col
  end

  close(f)

  return d
end

end # module

ext: jl


text:
module Tree

  include("node.jl")

end


text:
using LMCLUS
using Test
using Serialization
using Statistics

@testset "LMCLUS Types" begin

    s = LMCLUS.Separation(10., 10., 10., 10, 10.0, 10.0, 10)
    io = IOBuffer()
    serialize(io, s)
    seek(io, 0)
    d = deserialize(io)
    @test s.depth == d.depth
    @test s.discriminability == d.discriminability
    @test s.threshold == d.threshold
    @test s.globalmin == d.globalmin
    @test s.threshold == threshold(s)
    @test s.maxdist == d.maxdist
    @test s.mindist == d.mindist
    @test s.bins == d.bins

    m = Manifold(2, [2., 2.], [2. 2.; 2. 2.], [1, 2, 3], 10., 10.)
    @test outdim(m) == 2
    @test size(m) == 3
    @test points(m) == [1, 2, 3]
    @test mean(m) == [2., 2.]
    @test projection(m) == [2. 2.; 2. 2.]
    @test threshold(m) == (10., 10.)

    io = IOBuffer()
    serialize(io, m)
    seek(io, 0)
    n = deserialize(io)
    @test outdim(m) == outdim(n)
    @test size(m) == size(n)
    @test points(m) == points(n)
    @test mean(m) == mean(n)
    @test projection(m) == projection(n)
    @test threshold(m) == threshold(n)

end


text:
module TestPhylo

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

using Bio.Phylo
using LightGraphs

@testset "Phylogenies" begin

    @testset "Creation" begin
        # Create an unresolved star phylogeny.
        tips = [:A, :B, :C]
        stringtips = ["A", "B", "C"]
        tree = Phylogeny(tips)
        tree2 = Phylogeny(stringtips)

        names = [:A, :B, :C, :Root]

        @testset "Roots" begin
            @test isrooted(tree) == false
            @test isrerootable(tree)
            @test root(tree) == 4
            for i in 1:5
                r = i == 4
                @test isroot(tree, i) == r
            end
            for i in names
                r = i == :Root
                @test isroot(tree, i) == r
            end
        end

        @testset "Leaves" begin
            @test leaves(tree) == 1:3
            for i in 1:5
                l = i < 4
                @test isleaf(tree, i) == l
            end
            for i in names
                l = i != :Root
                @test isleaf(tree, i) == l
            end
        end

        @testset "Clades" begin
            @test clades(tree) == 4:5
            for i in 1:5
                c = i > 3
                p = i == 4
                @test isclade(tree, i) == c
                @test ispreterminal(tree, i) == p
                @test issemipreterminal(tree, i) == false
            end
            for i in names
                c = i == :Root
                @test isclade(tree, i) == c
            end
        end

        @testset "Parent - Child" begin
            for i in 1:5
                hasp = i < 4
                nop = i >= 4
                hasp2 = i < 4
                nchild = i == 4 ? 3 : 0
                hasc = i == 4
                hasc2 = i < 4
                noc = i != 4
                @test hasparent(tree, i) == hasp
                @test !hasparent(tree, i) == nop
                @test hasparent(tree, i, 4) == hasp2
                @test hasparent(tree, i, 5) == false
                @test nchildren(tree, i) == nchild
                @test haschildren(tree, i) == hasc
                @test !haschildren(tree, i) == noc
                @test haschild(tree, 4, i) == hasc2
                @test isredundant(tree, i) == false
                if i < 4
                    @test Phylo.parent(tree, i) == 4
                    @test children(tree, i) == Int[]
                else
                    @test_throws ErrorException Phylo.parent(tree, i)
                end
            end
            for i in names
                p = i != :Root
                nchild = i != :Root ? 0 : 3
                @test hasparent(tree, i) == p
                @test nchildren(tree, i) == nchild
            end
            @test hasparent(tree, collect(1:5)) == [true, true, true, false, false]
            @test hasparent(tree) == [true, true, true, false, false]
            @test !hasparent(tree, collect(1:5)) == [false, false, false, true, true]
            @test !hasparent(tree) == [false, false, false, true, true]
            @test nchildren(tree, collect(1:5)) == [0, 0, 0, 3, 0]
            @test nchildren(tree) == [0, 0, 0, 3, 0]
            @test haschildren(tree, collect(1:5)) == [false, false, false, true, false]
            @test haschildren(tree) == [false, false, false, true, false]
            @test !haschildren(tree, collect(1:5)) == [true, true, true, false, true]
            @test !haschildren(tree) == [true, true, true, false, true]
            @test isredundant(tree, collect(1:5)) == [false, false, false, false, false]
            @test isredundant(tree) == [false, false, false, false, false]
        end
        g = DiGraph(5)
        for i in 1:3
            add_edge!(g, Edge(4, i))
        end
        @test tree.graph == g
        @test tree.ntaxa == 3
        @test nv(tree.graph) == 5

        @testset "Copying" begin
            @test copy(tree) == tree
            @test copy(tree) !== tree
        end

        @testset "Branches" begin
            @test branchlength(branchlength!(tree, Edge(4, 1), 0.237), Edge(4, 1)) == 0.237
            @test parent_branch(tree, 1) == Edge(4, 1)
            @test child_branches(tree, 4) == [Edge(4, 1), Edge(4, 2), Edge(4, 3)]
            @test child_branches(rem_branch!(tree, Edge(4, 1)), 4) == [Edge(4, 2), Edge(4, 3)]
            @test child_branches(add_branch!(tree, Edge(4, 1), 0.5), 4) == [Edge(4, 1), Edge(4, 2), Edge(4, 3)]
        end

        @testset "Manipulation" begin
            @test Phylo.unconnected_clades(tree) == [5]
            @test Phylo.subtree_roots(tree) == [4]
            @test child_branches(Phylo.delete!(tree, 1), 4) == [Edge(4, 2), Edge(4, 3)]
            @test child_branches(Phylo.disconnect_root!(tree), 4) == []
        end

    end

end

@testset "Misc Functions" begin
    @test n_possible_rooted(3) == 3
    @test n_possible_rooted(4) == 15
    @test n_possible_rooted(5) == 105
    @test n_possible_unrooted(3) == 1
    @test n_possible_unrooted(4) == 3
    @test n_possible_unrooted(5) == 15
end

end # TestPhylo


text:
module AWSMetadata

using DataStructures: OrderedDict
using HTTP
using JSON

"""
    parse_aws_metadata()

Update the AWS service metadata definitions in `metadata.json` by parsing the `AWS SDK JS`.

Check to see if the AWS service exists in the metadata, if not then a new service has been
created and we should generate low and high level wrappers for it. Otherwise, check to see
if the service SHA hash matches what we keep on file, if they differ then the API has been
updated and we need to re-generate low and high level wrappers for the service.
"""
function parse_aws_metadata()
    # TODO:
    # - Only support the latest API version, need to add another filter!(files)
    # - Generate high-level wrapper for each service
    # - Only regenerate API definitions for services which have changed
    metadata_path = joinpath(@__DIR__, "metadata.json")
    metadata = JSON.parsefile(metadata_path, dicttype=OrderedDict)
    headers = ["User-Agent" => "JuliaCloud/AWSCore.jl"]
    url = "https://api.github.com/repos/aws/aws-sdk-js/contents/apis"
    req = HTTP.get(url, headers)
    files = JSON.parse(String(req.body), dicttype=OrderedDict)
    filter!(f -> occursin(r".normal.json$", f["name"]), files)  # Only get ${Service}.normal.json files
    data_changed = false
    services_modified = String[]

    function _process_service(file, version)
        data_changed = true
        push!(metadata, file["name"] => Dict("version" => version, "sha" => file["sha"]))
        push!(services_modified, file["name"])
    end

    for file in files
        filename = join(split(file["name"], '.')[1:end-2],'.')
        filename = split(filename, '-')
        service_name = join(filename[1:end-3], '-')
        version = join(filename[end-2:end], '-')

        filename = file["name"]

        if !haskey(metadata, filename)
            println(service_name, " does not exist in metadata.")
            _process_service(file, version)
        else
            if metadata[filename]["sha"] != file["sha"]
                println(service_name, " sha hashes do not match, updating.")
                _process_service(file, version)
            end
        end
    end

    if data_changed
        # Should move this generation to _process_service() so we don't regenerate for every
        # service, just the ones that change
        _generate_low_level_wrapper(files)
        _generate_high_level_wrapper(files)
        open(metadata_path, "w") do f
            print(f, json(OrderedDict(metadata), 2))
        end
    end
end

function _generate_low_level_wrapper(services)
    # TODO:
    # - How do we deal with multiple versions of the same service? e.g. cloudfront
    # - We now generate the same const over and over again
    service_definitions = _generate_service_definitions(services)

    template = """
    module AWSCorePrototypeServices

    include("../AWSCorePrototype.jl")

    $(join(service_definitions, "\n"))

    end
    """

    services_path = joinpath(@__DIR__, "AWSCorePrototypeServices.jl")
    open(services_path, "w") do f
        print(f, template)
    end
end

function _generate_service_definitions(services)
    service_definitions = String[]

    for service in services
        service_name = service["name"]
        println("Generating low level wrapper for $service_name")
        request = HTTP.get(service["download_url"])
        service_metadata = JSON.parse(String(request.body))["metadata"]

        definition = _generate_service_definition(service_metadata)
        if definition != nothing
            push!(service_definitions, definition)
        end
    end

    return service_definitions
end

function _generate_service_definition(service)
    println("Generating Metadata definitions for ", service["serviceId"])
    request_protocol = service["protocol"]
    service_name = service["endpointPrefix"]
    service_id = replace(lowercase(service["serviceId"]), ' ' => '_')
    api_version = service["apiVersion"]

    if request_protocol in ["rest-xml", "ec2", "query", "rest-json"]
        return "const $service_id = AWSCorePrototype.RestXMLService(\"$service_name\", \"$api_version\")"
    elseif request_protocol == "json"
        json_version = service["jsonVersion"]
        target = service["targetPrefix"]
        return "const $service_id = AWSCorePrototype.JSONService(\"$service_name\", \"$api_version\", \"$json_version\", \"$target\")"
    else
        # This should ping an alarm, we need to update this source code as its something we haven't seen before
        println(service_name, " uses a new protocol ", request_protocol)
    end
end

function _generate_rest_xml_high_level_wrapper(service_name, operations, shapes)
    # TODO:
    # - When generating S3 the documentation ends up with $ScanRange set in it which breaks things
    # - Pull down documentation for each input variable and write to the docstr
    function_definitions = String[]

    for operation in operations
        operation = operation[2]
        name = operation["name"]
        method = operation["http"]["method"]
        request_uri = operation["http"]["requestUri"]

        # Replace curly braces around parameters in the request_uri with a $, so Julia can pass the parameters along
        request_uri = replace(request_uri, '{' => "\$")
        request_uri = replace(request_uri, '}' => "")

        # Removes everything inbetween <> characters
        documentation = ""
        if haskey(operation, "documentation")
            documentation = operation["documentation"]
            documentation = replace(documentation, r"\<.*?\>" => "")
        end

        required_parameters = ""

        if haskey(operation, "input")
            input_shape = shapes[operation["input"]["shape"]]

            if haskey(input_shape, "required")
                required_parameters = input_shape["required"]
            end
        end

        definition = """
        \"\"\"
        $documentation
        \"\"\"
        $name($(join(required_parameters, ", "))) = $service_name(\"$method\", \"$request_uri\")
        """

        push!(function_definitions, definition)
    end

    service_path = joinpath(@__DIR__, "services/$service_name.jl")
    open(service_path, "w") do f
        println(f, "module aws_$service_name")
        println(f, "include(\"AWSCorePrototypeServices.jl\")")
        println(f, "using .AWSCorePrototypeServices: $service_name\n")
        print(f, join(function_definitions, "\n"))
        println(f, "end")
    end
end

function _generate_high_level_wrapper(services)
    # TODO:
    # - Create functions for query, rest-json, and json protocols
    #   These only seem to differ from rest-xml with how their shapes are defined
    #   We only need to pull the required uri params down for function definitions
    for service in services
        service_name = service["name"]
        println("Generating high level wrapper for $service_name")
        service = JSON.parse(String(HTTP.get(service["download_url"]).body))
        service_name = lowercase(service["metadata"]["serviceId"])
        service_name = replace(service_name, ' ' => '_')
        operations = service["operations"]
        shapes = service["shapes"]

        protocol = service["metadata"]["protocol"]

        if protocol in ["rest-xml"]
            _generate_rest_xml_high_level_wrapper(service_name, operations, shapes)
        end
    end
end
end

max_stars_repo_name: nicoleepp/JuliaCloud-AWS-State

max_issues_repo_name: nicoleepp/JuliaCloud-AWS-State

max_forks_repo_name: nicoleepp/JuliaCloud-AWS-State


text:
"""
    solve(problem::WaveguideProblem; [ρ₀, Nouts], kwargs...)

Solves the `WaveguideProblem`, using `QuantumOptics.timeevolution.master_nh[_dynamic]`.

# Arguments
- `problem`: the waveguide problem
- `Nouts::Array{Int}` (optional): specifies the dimension of the output Fock spaces
- `ρ₀` (optional): Alternative initial state for the simulation.
	Overrides the initial state provided by the `WavePacket` and the quantum system.
- Additional keyword arguments are passed to `QuantumOptics.jl`'s master equation solver,
  see [it's documentation](https://docs.qojulia.org/api/#QuantumOptics.timeevolution.master)

# Returns
- ts:    The times on which the output is determined.
- ρ(t):  The density matrix at each point in time. The Hilbert space depends on the specific
         problem; if unsure check out `basis(ρ[1])`.
"""
function solve(problem::_DrivenProblem{O, Coherent, ContinuousWave};
		ρ₀ = nothing, kwargs...) where {O}
	ψ₀, Ls, H_nh = _generateLindbladian(problem)
	ρ₀ = isnothing(ρ₀) ? ψ₀ : ρ₀

	return timeevolution.master_nh(problem.ts, ρ₀, H_nh, Ls; kwargs...)
end

function solve(problem::WaveguideProblem; Nouts::Union{Nothing,Array} = nothing,
		ρ₀ = nothing,  kwargs...)
	if problem isa _ScatterProblem
		if !isnothing(ρ₀)
			!isnothing(Nouts) && @warn "`Nouts` and `ρ₀` simultaneously set. Ignoring `Nouts`"
			Nouts = [basis.N for basis ∈ basis(ρ₀).bases[3:end]]
		end
		isnothing(Nouts) && ( Nouts = _expectedOutputPhotons(problem) )
		@assert length(Nouts) == length(problem.ψₒ)
	end

	ψ₀, Ls, H_nh = _generateLindbladian(problem, Nouts)
	ρ₀ = isnothing(ρ₀) ? ψ₀ : ρ₀

	function Operators(t, ψ)
		Ls_t = [Ls[1](t), Ls[2:end]...]
		return H_nh(t), H_nh(t)', Ls_t, dagger.(Ls_t)
	end

	return timeevolution.master_nh_dynamic(problem.ts, ρ₀, Operators; kwargs...)
end

function _expectedOutputPhotons(problem::_ScatterProblem{O, Coherent, ContinuousWave}) where {O}
	if problem.displace_output
		return fill(ceil(Int, 3 * sqrt(problem.ts[end] - problem.ts[1]) ), length(problem.ψₒ) )
	end
	fill(coherent_cutoff( sqrt( (problem.ts[end] - problem.ts[1]) * abs2(problem.ψᵢ.α) ) ),
		length(problem.ψₒ)
	)
end

function _expectedOutputPhotons(problem::_ScatterProblem{O, Coherent, WavePacket{Coherent}}) where {O}
	if problem.displace_output
		return fill(ceil(Int, 3 * sqrt(problem.ts[end] - problem.ts[1]) ), length(problem.ψₒ) )
	end
	N = problem.ψᵢ.state.N_cutoff
	fill(N, length(problem.ψₒ))
end

function _expectedOutputPhotons(problem::_ScatterProblem{O, S, WavePacket{S}}) where {O, S <: Displaced}
	if problem.displace_output
		return fill(
			roblem.ψᵢ.state.N_cutoff + ceil(Int, 3 * sqrt(problem.ts[end] - problem.ts[1]) ),
			length(problem.ψₒ)
		)
	end
	N = problem.ψᵢ.state.N_cutoff + coherent_cutoff(problem.ψᵢ.state.α)
	fill(N, length(problem.ψₒ))
end

function _expectedOutputPhotons(problem::_ScatterProblem{O, S, WavePacket{S}}) where {O, S}
	N = problem.ψᵢ.state.N_cutoff
	fill(N, length(problem.ψₒ))
end

# ------------------------------------------------------------
# _DrivenProblems

function _generateLindbladian(problem::_DrivenProblem{O, Coherent, ContinuousWave}) where {O}
	# ----------------------------------------
	# Scatter System
	H_sys, σ⁻, σ⁺ = problem.H, problem.σ, problem.σ'
	σ⁺σ⁻ = σ⁺ * σ⁻

	# ----------------------------------------
	# Constructing Disspator and Hamiltonian
	α  = problem.ψᵢ.α
	atomic_dissipation = isempty(problem.Ls) ?
		-0.5im * σ⁺σ⁻ :
		-0.5im * σ⁺σ⁻ - 0.5im * sum([L'*L for L ∈ problem.Ls])

	H_nh = atomic_dissipation - 1.0im * (σ⁺ * α - σ⁻ * α) + H_sys

	ψ₀ = problem.system_state
	return ψ₀, [σ⁻, problem.Ls...], H_nh
end

function _generateLindbladian(problem::_DrivenProblem{O, Coherent, WavePacket{Coherent}}, _) where {O}
	# ----------------------------------------
	# Scatter System
	H_sys, σ⁻, σ⁺ = problem.H, problem.σ, problem.σ'
	σ⁺σ⁻ = σ⁺ * σ⁻

	# ----------------------------------------
	# Constructing Disspator and Hamiltonian
	α  = problem.ψᵢ.state.α
	mf = problem.ψᵢ.mode.modeFunction
	atomic_dissipation = isempty(problem.Ls) ?
		-0.5im * σ⁺σ⁻ :
		-0.5im * σ⁺σ⁻ - 0.5im * sum([L'*L for L ∈ problem.Ls])

	H_nh(t) = atomic_dissipation - 1.0im * (σ⁺ * α * mf(t) - σ⁻ * α * mf(t)) + H_sys

	ψ₀ = problem.system_state
	return ψ₀, [t -> σ⁻, problem.Ls...], H_nh
end

function _generateLindbladian(problem::_DrivenProblem{O, T, WavePacket{T}}, _) where {O, T}
	# ----------------------------------------
	# Scatter System
	H_sys, σ⁻, σ⁺ = problem.H, problem.σ, problem.σ'
	σ⁺σ⁻ = σ⁺ * σ⁻
	idₛ = one(basis(H_sys))

	# ----------------------------------------
	# Input System
	inBasis = FockBasis(problem.ψᵢ.state.N_cutoff)
	gᵢ      = problem.ψᵢ.mode.gᵢ
	idᵢ     = identityoperator(inBasis)
	aᵢ, aᵢᵀ = destroy(inBasis) ⊗ idₛ, create(inBasis) ⊗ idₛ
	aᵢᵀaᵢ   = aᵢᵀ * aᵢ

	# ----------------------------------------
	# Embed Operators and abbreviations
	H_sys, σ⁻, σ⁺, σ⁺σ⁻ = idᵢ ⊗ H_sys, idᵢ ⊗ σ⁻, idᵢ ⊗ σ⁺, idᵢ ⊗ σ⁺σ⁻
	aᵢᵀσ⁻,  aᵢσ⁺ = aᵢᵀ * σ⁻, aᵢ * σ⁺

	# ----------------------------------------
	# Constructing Disspator and Hamiltonian
	L  = t -> σ⁻ + gᵢ(t)aᵢ
	Lᵀ = t -> σ⁺ + gᵢ(t)aᵢᵀ

	atomic_dissipation = isempty(problem.Ls) ?
		-0.5im * σ⁺σ⁻ :
		-0.5im * σ⁺σ⁻ - (0.5im * idᵢ) ⊗ sum([L'*L for L ∈ problem.Ls])

	H_nh = _createHamiltonian(problem, atomic_dissipation, aᵢᵀaᵢ, aᵢσ⁺, σ⁺, σ⁻, H_sys)

	ψ₀ = Ket(inBasis, createState(problem.ψᵢ.state)) ⊗ problem.system_state
	return ψ₀, [L, [idᵢ ⊗ L for L ∈ problem.Ls]...], H_nh
end

function _createHamiltonian(problem::_DrivenProblem{O, S, WavePacket{S}},
		atomic_dissipation, aᵢᵀaᵢ, aᵢσ⁺, σ⁺, σ⁻, H_sys) where {O, S <: Displaced}
	gᵢ, mf, α = problem.ψᵢ.mode.gᵢ, problem.ψᵢ.mode.modeFunction, problem.ψᵢ.state.α

	return t -> atomic_dissipation - 0.5im * gᵢ(t)^2 * aᵢᵀaᵢ -
		1.0im * (σ⁺ * α * mf(t) - σ⁻ * α' * mf(t)') - 1.0im * gᵢ(t)aᵢσ⁺ + H_sys
end

function _createHamiltonian(problem::_DrivenProblem{O, S, WavePacket{S}},
		atomic_dissipation, aᵢᵀaᵢ, aᵢσ⁺, σ⁺, σ⁻, H_sys) where {O, S <: NonDisplaced}
	gᵢ = problem.ψᵢ.mode.gᵢ

	return t -> atomic_dissipation - 0.5im * gᵢ(t)^2 * aᵢᵀaᵢ - 1.0im * gᵢ(t)aᵢσ⁺ + H_sys
end

# ------------------------------------------------------------
# _ScatterProblems

function _generateLindbladian(problem::_ScatterProblem{O, Coherent, ContinuousWave}, Nouts) where {O}
	_generateLindbladian(
		WaveguideProblem(
			[problem.H, problem.Ls, problem.σ, problem.system_state],
			ContinuousWave(t -> problem.ψᵢ.α),
			problem.ψₒ,
			problem.ts,
			problem.displace_output
		),
		Nouts
	)
end

function _generateLindbladian(problem::_ScatterProblem{O, Coherent, WavePacket{Coherent}}, Nouts) where {O}
	# ----------------------------------------
	# Output System
	outBasis, groundstateOutput, aₒs, aₒᵀs = _generateOutputCavities(Nouts)
	idₒ = one(outBasis)

	# ----------------------------------------
	# Scatter System
	H_sys, σ⁻, σ⁺ = problem.H ⊗ idₒ, problem.σ ⊗ idₒ, problem.σ' ⊗ idₒ
	σ⁺σ⁻ = σ⁺ * σ⁻
	idₛ = one(basis(problem.H))

	# ----------------------------------------
	# Embed Operators and abbreviations
	aₒs, aₒᵀs = [idₛ ⊗ a for a ∈ aₒs], [idₛ ⊗ aᵀ for aᵀ ∈ aₒᵀs]

	gₒs = [mode.gₒ for mode in problem.ψₒ]

	gₒaₒs(t)  = sum([gₒ(t) for gₒ ∈ gₒs] .* aₒs)
	gₒaₒᵀs(t) = sum([gₒ(t) for gₒ ∈ gₒs] .* aₒᵀs)

	aₒᵀσ⁻s = [aₒᵀ * σ⁻ for aₒᵀ in aₒᵀs]
	aₒᵀaₒs = aₒᵀs .* aₒs
	# TODO: maybe it is faster to do `gₒaₒs(t)  = idₛ ⊗ sum(gₒs(t) .* aₒs)` and leave aₒs unchanged
	# due to the sum

	# ----------------------------------------
	# Constructing Disspator and Hamiltonian
	L  = t -> (σ⁻ + gₒaₒs(t) )
	Lᵀ = t -> (σ⁺ + gₒaₒᵀs(t) )

	α  = problem.ψᵢ.state.α
	mf = problem.ψᵢ.mode.modeFunction
	atomic_dissipation = isempty(problem.Ls) ?
		-0.5im * σ⁺σ⁻ :
		-0.5im * σ⁺σ⁻ - 0.5im * sum([L'*L for L ∈ problem.Ls]) ⊗ idₒ

	if problem.displace_output
		H_nh = function(t)
			return atomic_dissipation - 0.5im * sum([gₒ(t)^2 for gₒ ∈ gₒs] .* aₒᵀaₒs) -
				1.0im * sum([gₒ(t) for gₒ ∈ gₒs] .* aₒᵀσ⁻s) -
				1.0im * (σ⁺ * α * mf(t) - σ⁻ * α' * mf(t)) +
				H_sys
		end
	else
		H_nh = function(t)
			return atomic_dissipation - 0.5im * sum([gₒ(t)^2 for gₒ ∈ gₒs] .* aₒᵀaₒs) -
				1.0im * sum([gₒ(t) for gₒ ∈ gₒs] .* aₒᵀσ⁻s) -
				1.0im * (Lᵀ(t) * α * mf(t) - L(t) * α' * mf(t)) +
				H_sys
		end
	end

	ψ₀ = problem.system_state ⊗ groundstateOutput
	return ψ₀, [t -> L(t), [L ⊗ idₒ for L ∈ problem.Ls]...], H_nh
end

function _generateLindbladian(problem::_ScatterProblem{O, T, WavePacket{T}}, Nouts) where {O, T}
	# ----------------------------------------
	# Output System
	outBasis, groundstateOutput, aₒs, aₒᵀs = _generateOutputCavities(Nouts)
	idₒ = one(outBasis)

	# ----------------------------------------
	# Scatter System
	H_sys, σ⁻, σ⁺  = problem.H, problem.σ, problem.σ'
	σ⁺σ⁻ = σ⁺ * σ⁻
	idₛ = one(basis(problem.H))

	# ----------------------------------------
	# Input System
	inBasis = FockBasis(problem.ψᵢ.state.N_cutoff)
	gᵢ      = problem.ψᵢ.mode.gᵢ
	idᵢ     = identityoperator(inBasis)
	aᵢ, aᵢᵀ = destroy(inBasis) ⊗ idₛ ⊗ idₒ, create(inBasis) ⊗ idₛ ⊗ idₒ
	aᵢᵀaᵢ = aᵢᵀ * aᵢ

	# ----------------------------------------
	# Embed Operators and abbreviations
	H_sys, σ⁻, σ⁺, σ⁺σ⁻ = idᵢ ⊗ H_sys ⊗ idₒ, idᵢ ⊗ σ⁻ ⊗ idₒ, idᵢ ⊗ σ⁺ ⊗ idₒ, idᵢ ⊗ σ⁺σ⁻ ⊗ idₒ
	aₒs, aₒᵀs = [idᵢ ⊗ idₛ ⊗ a for a ∈ aₒs], [idᵢ ⊗ idₛ ⊗ aᵀ for aᵀ ∈ aₒᵀs]

	gₒs = [mode.gₒ for mode in problem.ψₒ]

	gₒaₒs(t)  = sum([gₒ(t)' for gₒ ∈ gₒs] .* aₒs)
	gₒaₒᵀs(t) = sum([gₒ(t)  for gₒ ∈ gₒs] .* aₒᵀs)

	aᵢσ⁺   = aᵢ * σ⁺
	aₒᵀaᵢs = [aₒᵀ * aᵢ for aₒᵀ in aₒᵀs]
	aₒᵀσ⁻s = [aₒᵀ * σ⁻ for aₒᵀ in aₒᵀs]
	aₒᵀaₒs = aₒᵀs .* aₒs
	# TODO: maybe it is faster to do `gₒaₒs(t)  = idₛ ⊗ sum(gₒs(t) .* aₒs)` and leave aₒs unchanged
	# due to the sum

	# ----------------------------------------
	# Constructing Disspator and Hamiltonian
	L  = t -> σ⁻ + gᵢ(t)aᵢ  + gₒaₒs(t)

	atomic_dissipation = isempty(problem.Ls) ?
		-0.5im * σ⁺σ⁻ :
		-0.5im * σ⁺σ⁻ - (0.5im * idᵢ) ⊗ sum([L'*L for L ∈ problem.Ls]) ⊗ idₒ

	H_nh = _createHamiltonian(problem, atomic_dissipation,
		aₒᵀaₒs, aᵢᵀaᵢ, aᵢσ⁺, aₒᵀaᵢs, aₒᵀσ⁻s, σ⁺, σ⁻, gₒaₒᵀs, gₒaₒs, H_sys
	)

	ψ₀ = Ket(inBasis, createState(problem.ψᵢ.state)) ⊗ problem.system_state ⊗ groundstateOutput
	return ψ₀, [t -> L(t), [idᵢ ⊗ L ⊗ idₒ for L ∈ problem.Ls]...], H_nh
end

function _generateOutputCavities(Nouts)
	# Preconverting into `CompositeBasis` allows us to `embed` into it, even for a single system
	outBasis  = reduce(⊗, FockBasis.(Nouts) .|> CompositeBasis)

	groundstateOutput = reduce(⊗, [fockstate(FockBasis(N), 0) for N in Nouts])

	aₒs  = [embed(outBasis, i, destroy(FockBasis(N))) for (i, N) in enumerate(Nouts)]
	aₒᵀs = [embed(outBasis, i, create(FockBasis(N)))  for (i, N) in enumerate(Nouts)]

	return outBasis, groundstateOutput, aₒs, aₒᵀs
end

function _createHamiltonian(problem::_ScatterProblem{O, S, WavePacket{S}},
		atomic_dissipation, aₒᵀaₒs, aᵢᵀaᵢ, aᵢσ⁺, aₒᵀaᵢs, aₒᵀσ⁻s, σ⁺, σ⁻, gₒaₒᵀs, gₒaₒs,
		H_sys) where {O, S<: Displaced}
	gᵢ, mf, α = problem.ψᵢ.mode.gᵢ, problem.ψᵢ.mode.modeFunction, problem.ψᵢ.state.α
	gₒs   = [mode.gₒ for mode in problem.ψₒ]

	return t -> atomic_dissipation -
		0.5im * (gᵢ(t)^2 * aᵢᵀaᵢ + sum([gₒ(t)^2 for gₒ ∈ gₒs] .* aₒᵀaₒs) ) -
		1.0im * (gᵢ(t)aᵢσ⁺ + sum([gₒ(t) for gₒ ∈ gₒs] .* aₒᵀσ⁻s) +
			sum(gᵢ(t) .* [gₒ(t) for gₒ ∈ gₒs] .* aₒᵀaᵢs) ) -
		1.0im * (( σ⁺ + gₒaₒᵀs(t) ) * α * mf(t) - ( σ⁻ + gₒaₒs(t) )  * α' * mf(t)') + H_sys
end

function _createHamiltonian(problem::_ScatterProblem{O, S, WavePacket{S}},
		atomic_dissipation, aₒᵀaₒs, aᵢᵀaᵢ, aᵢσ⁺, aₒᵀaᵢs, aₒᵀσ⁻s, σ⁺, σ⁻, gₒaₒᵀs, gₒaₒs,
		H_sys) where {O, S<: NonDisplaced}
	gᵢ    = problem.ψᵢ.mode.gᵢ
	gₒs   = [mode.gₒ for mode in problem.ψₒ]

	return t -> atomic_dissipation -
		0.5im * (gᵢ(t)^2 * aᵢᵀaᵢ + sum([gₒ(t)^2 for gₒ ∈ gₒs] .* aₒᵀaₒs) ) -
		1.0im * (gᵢ(t)aᵢσ⁺ + sum([gₒ(t) for gₒ ∈ gₒs] .* aₒᵀσ⁻s) +
			sum(gᵢ(t) .* [gₒ(t) for gₒ ∈ gₒs] .* aₒᵀaᵢs) ) +
		H_sys
end


text:
@debug "BedgraphFiles loading DataFrames integration."

function Base.convert(::Type{Bedgraph.Record}, row::DataFrames.DataFrameRow) :: Bedgraph.Record
    return Bedgraph.Record(row[1], row[2], row[3], row[4]) # Note: using index to allow flexible column names.
end

function Base.convert(::Type{Vector{Bedgraph.Record}}, df::DataFrames.DataFrame) :: Vector{Bedgraph.Record}

    records = Vector{Bedgraph.Record}(undef, size(df)[1])

    for (i, row) in enumerate(eachrow(df))
        records[i] = convert(Bedgraph.Record, row)
    end

    return records
end


text:
"""
    CodeBook{U<:Unsigned,T}

The codebook structure. It holds the codes corresponding
to the vector prototypes and the mapping bethween the codes
and prototypes.

# Fields
  * `codes::Vector{U}` the codes
  * `vectors::Matrix{T}` the prototypes
  * `codemap::Dict{U,Int}` mapping from code to column in `vectors`
"""
struct CodeBook{U<:Unsigned,T}
    codes::Vector{U}
    vectors::Matrix{T}
    codemap::Dict{U,Int}
end

# Constructors
CodeBook(codes::Vector{U}, vectors::Matrix{T}) where {U,T} = begin
    length(codes) != size(vectors, 2) &&
        throw(DimensionMismatch("Dimension of codes and vectors are inconsistent"))
    codemap = Dict{U, Int}(c => i for (i,c) in enumerate(codes))
    return CodeBook(codes, vectors, codemap)
end

CodeBook{U}(vectors::Matrix{T}) where {U,T} = begin
    k = size(vectors,2)
    nb = TYPE_TO_BITS[U]
    nbk = TYPE_TO_BITS[quantized_eltype(k)]
    nbk > nb && throw(ErrorException("$k vectors cannot be coded in $nb bits."))
    codes = Vector{U}(0:k-1)
    codemap = Dict{U, Int}(c => i for (i,c) in enumerate(codes))
    return CodeBook(codes, vectors, codemap)
end

CodeBook{U,T}(vectors::Matrix{T2}) where {U,T<:Integer,T2} =
    CodeBook{U}(round.(T, vectors))

CodeBook{U,T}(vectors::Matrix{T2}) where {U,T<:AbstractFloat,T2} =
    CodeBook{U}(convert.(T, vectors))

CodeBook(vectors::Matrix{T}) where {T} =
    CodeBook{quantized_eltype(size(vectors, 2)), T}(vectors)

# Show method
Base.show(io::IO, cb::CodeBook{U,T}) where {U,T} = begin
    len_vecs, num_codes = size(cb.vectors)
    print(io, "CodeBook{$U,$T} $(num_codes) codes, $(len_vecs)-element vectors")
end

# Basic indexing: for a key, get the vector
Base.getindex(cb::CodeBook{U,T}, key::U) where {U,T} =
    cb.vectors[:, cb.codemap[key]]

"""
    encode(codebook, aa [;distance=DEFAULT_DISTANCE])

Encodes the input array aa using `distance` to calculate the
closest vector prototype from the `codebook`.
"""
function encode(cb::CodeBook{U,T},
                aa::AbstractMatrix{T};
                distance::Distances.PreMetric=DEFAULT_DISTANCE
               ) where {U,T}
    dists = Distances.pairwise(distance, cb.vectors, aa, dims=2)  # distances between codebook vectors and data
    best_idx = getindex(findmin(dists, dims=1), 2)      # get positions vector of all minima
    best_idx_cols = vec(getindex.(best_idx, 1))         # get row values from positions
    return cb.codes[best_idx_cols]                      # get codes corresponding to minima
end

"""
    decode(codebook, codes)

Returns a partial array reconstruction for using the `codes` and `codebook`.
"""
function decode(cb::CodeBook{U,T}, codes::Vector{U}) where {U,T}
    m, n = size(cb.vectors, 1), length(codes)
    aa = Matrix{T}(undef, m, n)
    @inbounds @simd for i in 1:n
        aa[:,i] = cb[codes[i]]
    end
    return aa
end

"""
    build_codebooks(X, k, m, U [;method=DEFAULT_METHOD, distance=DEFAULT_DISTANCE, kwargs])

Generates `m` codebooks of `k` prototypes each for the input matrix `X`
using the algorithm specified my `method` and distance `distance`.
Specific codebook aglorithm keyword arguments can be specified as well.

# Arguments
  * `X::AbstractMatrix{T}` input matrix of type `T`
  * `k::Int` number of prototypes/codebook
  * `m::Int` number of codebooks
  * `U::Type{<:Unsigned}` type for codebook codes

# Keyword arguments
  * `method::Symbol` the algorithm to be employed for codebook
generation; possible values are `:sample` (default), `:pq` for
classical k-means clustering codebooks and `:opq` for 'cartesian'
k-means clustering codebooks
  * `distance::PreMetric` the distance to be used in the
codebook generation methods and data encoding
  * `kwargs...` other codebook generation algorithm specific
keyword arguments such as `maxiter::Int`.
"""
function build_codebooks(X::AbstractMatrix{T},
                         k::Int,
                         m::Int,
                         ::Type{U};
                         method::Symbol=DEFAULT_METHOD,
                         distance::Distances.PreMetric=DEFAULT_DISTANCE,
                         kwargs...) where {U,T}
    nrows, ncols = size(X)
    k = min(k, ncols)              # number of codes for a quantizer
    m = ifelse(method == :rvq, m, min(m, nrows))  # number of quantizers
    R = diagm(0=>ones(T, nrows))
    if method == :sample
        vectors = sampling_codebooks(X, k, m)  # does not use distances
    elseif method == :pq
        vectors = pq_codebooks(X, k, m, distance=distance; kwargs...)
    elseif method == :opq
        vectors, R = opq_codebooks(X, k, m, distance=distance; kwargs...)
    elseif method == :rvq
        vectors = rvq_codebooks(X, k, m, distance=distance; kwargs...)
    else
        @error "Unknown codebook generation method '$(method)'"
    end
    return map(CodeBook{U,T}, vectors), R
end


text:
import Pkg

path_to_env = joinpath(@__DIR__, "../../")
Pkg.activate(path_to_env)

using CytofResearch
using BSON
using Distributions
using LaTeXStrings
using DataFrames
using CSV
using DelimitedFiles
include(joinpath(@__DIR__, "salso.jl"))

import PyCall, PyPlot
const plt = PyPlot.plt
PyPlot.matplotlib.use("Agg")

#= Interactive plot
PyPlot.matplotlib.use("TkAgg")
=#
#= Non-interactive plot
PyPlot.matplotlib.use("Agg")
=#

# TODO
# Steal things from: Cytof5/sims/repfam_fs/test

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
abstract type InferenceOptimizer{T} end

## Abstract type for analytical optimiser (closed forms are known)
abstract type AOptimizer{T} <: InferenceOptimizer{T} end
## Abstract type for numerical optimiser (need for numerical integration)
abstract type NOptimizer{T} <: InferenceOptimizer{T} end

mutable struct AVIOptimizer{T<:Real,O} <: AOptimizer{T}
    optimiser::O #Learning rate for stochastic updates
    ∇η₁::Vector{T}
    ∇η₂::Matrix{T}
    function AVIOptimizer{T}(n::Int, opt::O) where {T,O}
        return new{T,O}(opt, zeros(T, n), zeros(T, n, n))
    end
end

mutable struct NVIOptimizer{T<:Real,O} <: NOptimizer{T}
    optimiser::O #Learning rate for stochastic updates
    ∇η₁::Vector{T}
    ∇η₂::Matrix{T}
    ν::Vector{T} #Derivative -<dv/dx>_qn
    λ::Vector{T} #Derivative  <d²V/dx²>_qm
    L::LowerTriangular{T,Matrix{T}}
    function NVIOptimizer{T}(n::Int, b::Int, opt::O) where {T,O}
        return new{T,O}(
            opt,
            zeros(T, n),
            zeros(T, n, n),
            zeros(T, b),
            zeros(T, b),
            LowerTriangular(diagm(sqrt(0.5) * ones(T, n))),
        )
    end
end

mutable struct SOptimizer{T<:Real,O} <: AOptimizer{T}
    optimiser::O
end

function SOptimizer{T}(opt::O) where {T,O}
    return SOptimizer{T,O}(opt)
end


text:
function bar(d, kw_args)
    x, y = d[:x], d[:y]
    nx, ny = length(x), length(y)
    axis = d[:subplot][isvertical(d) ? :xaxis : :yaxis]
    cv = [discrete_value!(axis, xi)[1] for xi=x]
    x = if nx == ny
        cv
    elseif nx == ny + 1
        0.5diff(cv) + cv[1:end-1]
    else
        error("bar recipe: x must be same length as y (centers), or one more than y (edges).\n\t\tlength(x)=$(length(x)), length(y)=$(length(y))")
    end
    if haskey(kw_args, :stroke_width) # stroke is inside for bars
        #kw_args[:stroke_width] = -kw_args[:stroke_width]
    end
    # compute half-width of bars
    bw = nothing
    hw = if bw == nothing
        ignorenan_mean(diff(x))
    else
        Float64[_cycle(bw,i)*0.5 for i=1:length(x)]
    end

    # make fillto a vector... default fills to 0
    fillto = d[:fillrange]
    if fillto == nothing
        fillto = 0
    end
    # create the bar shapes by adding x/y segments
    positions, scales = Array{Point2f0}(ny), Array{Vec2f0}(ny)
    m = Reactive.value(kw_args[:model])
    sx, sy = m[1,1], m[2,2]
    for i=1:ny
        center = x[i]
        hwi = abs(_cycle(hw,i)); yi = y[i]; fi = _cycle(fillto,i)
        if Plots.isvertical(d)
            sz = (hwi*sx, yi*sy)
        else
            sz = (yi*sx, hwi*2*sy)
        end
        positions[i] = (center-hwi*0.5, fi)
        scales[i] = sz
    end

    kw_args[:scale] = scales
    kw_args[:offset] = Vec2f0(0)
    visualize((GLVisualize.RECTANGLE, positions), Style(:default), kw_args)
    #[]
end

const _box_halfwidth = 0.4

notch_width(q2, q4, N) = 1.58 * (q4-q2)/sqrt(N)

function boxplot(d, kw_args)
    kwbox = copy(kw_args)
    range = 1.5; notch = false
    x, y = d[:x], d[:y]
    glabels = sort(collect(unique(x)))
    warning = false
    outliers_x, outliers_y = zeros(0), zeros(0)

    box_pos = Point2f0[]
    box_scale = Vec2f0[]
    outliers = Point2f0[]
    t_segments = Point2f0[]
    m = Reactive.value(kw_args[:model])
    sx, sy = m[1,1], m[2,2]
    for (i,glabel) in enumerate(glabels)
        # filter y
        values = y[filter(i -> _cycle(x,i) == glabel, 1:length(y))]
        # compute quantiles
        q1,q2,q3,q4,q5 = quantile(values, linspace(0,1,5))
        # notch
        n = Plots.notch_width(q2, q4, length(values))
        # warn on inverted notches?
        if notch && !warning && ( (q2>(q3-n)) || (q4<(q3+n)) )
            warn("Boxplot's notch went outside hinges. Set notch to false.")
            warning = true # Show the warning only one time
        end

        # make the shape
        center = Plots.discrete_value!(d[:subplot][:xaxis], glabel)[1]
        hw = d[:bar_width] == nothing ? Plots._box_halfwidth*2 : _cycle(d[:bar_width], i)
        l, m, r = center - hw/2, center, center + hw/2

        # internal nodes for notches
        L, R = center - 0.5 * hw, center + 0.5 * hw
        # outliers
        if Float64(range) != 0.0  # if the range is 0.0, the whiskers will extend to the data
            limit = range*(q4-q2)
            inside = Float64[]
            for value in values
                if (value < (q2 - limit)) || (value > (q4 + limit))
                    push!(outliers, (center, value))
                else
                    push!(inside, value)
                end
            end
            # change q1 and q5 to show outliers
            # using maximum and minimum values inside the limits
            q1, q5 = ignorenan_extrema(inside)
        end
        # Box
        if notch
            push!(t_segments, (m, q1), (l, q1), (r, q1), (m, q1), (m, q2))# lower T
            push!(box_pos, (l, q2));push!(box_scale, (hw*sx, n*sy)) # lower box
            push!(box_pos, (l, q4));push!(box_scale, (hw*sx, n*sy)) # upper box
            push!(t_segments, (m, q5), (l, q5), (r, q5), (m, q5), (m, q4))# upper T

        else
            push!(t_segments, (m, q2), (m, q1), (l, q1), (r, q1))# lower T
            push!(box_pos, (l, q2)); push!(box_scale, (hw*sx, (q3-q2)*sy)) # lower box
            push!(box_pos, (l, q4)); push!(box_scale, (hw*sx, (q3-q4)*sy)) # upper box
            push!(t_segments, (m, q4), (m, q5), (r, q5), (l, q5))# upper T
        end
    end
    kwbox = Dict{Symbol, Any}(
        :scale => box_scale,
        :model => kw_args[:model],
        :offset => Vec2f0(0),
    )
    extract_marker(d, kw_args)
    outlier_kw = Dict(
        :model => kw_args[:model],
        :color =>  scalar_color(d, :fill),
        :stroke_width => Float32(d[:markerstrokewidth]),
        :stroke_color => scalar_color(d, :markerstroke),
    )
    lines_kw = Dict(
        :model => kw_args[:model],
        :stroke_width =>  d[:linewidth],
        :stroke_color =>  scalar_color(d, :fill),
    )
    vis1 = GLVisualize.visualize((GLVisualize.RECTANGLE, box_pos), Style(:default), kwbox)
    vis2 = GLVisualize.visualize((GLVisualize.CIRCLE, outliers), Style(:default), outlier_kw)
    vis3 = GLVisualize.visualize(t_segments, Style(:linesegment), lines_kw)
    [vis1, vis2, vis3]
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
function add_tech_size_constraints(m, p; _n="")

    # PV techs can be constrained by space available based on location at site (roof, ground, both)
    @constraint(m, [loc in p.pvlocations],
        sum(m[Symbol("dvSize"*_n)][t] * p.pv_to_location[t][loc] for t in p.techs.pv) <= p.maxsize_pv_locations[loc]
    )

    # max size limit
    @constraint(m, [t in p.techs.all],
        m[Symbol("dvSize"*_n)][t] <= p.max_sizes[t]
    )

    ##Constraint (7c): Minimum size for each tech
    @constraint(m, [t in p.techs.all],
        m[Symbol("dvSize"*_n)][t] >= p.min_sizes[t]
    )

    @constraint(m, [t in p.techs.all],
        m[Symbol("dvPurchaseSize"*_n)][t] >= m[Symbol("dvSize"*_n)][t] - p.existing_sizes[t]
    )

    ## Constraint (7d): Non-turndown technologies are always at rated production
    @constraint(m, [t in p.techs.no_turndown, ts in p.time_steps],
        m[Symbol("dvRatedProduction"*_n)][t,ts] == m[Symbol("dvSize"*_n)][t]
    )
end

function add_no_curtail_constraints(m, p; _n="")
    for t in p.techs.no_curtail
        for ts in p.time_steps
            fix(m[Symbol("dvCurtail"*_n)][t, ts] , 0.0, force=true)
        end
    end
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
module Utils

export launch_config, work_layout, launch!
export cell_advection_timescale
export TimeStepWizard, update_Δt!
export prettytime, pretty_filesize
export tupleit, parenttuple, datatuple, datatuples
export validate_intervals, time_to_run
export ordered_dict_show
export with_tracers
export versioninfo_with_gpu, oceananigans_versioninfo
export instantiate
export TimeInterval, IterationInterval, WallTimeInterval, SpecifiedTimes, AndSchedule, OrSchedule

import CUDA  # To avoid name conflicts
# import AMDGPU

#####
##### Misc. small utils
#####

instantiate(x) = x
instantiate(X::DataType) = X()

#####
##### Include utils
#####

include("prettysummary.jl")
include("kernel_launching.jl")
include("cell_advection_timescale.jl")
include("prettytime.jl")
include("pretty_filesize.jl")
include("tuple_utils.jl")
include("output_writer_diagnostic_utils.jl")
include("ordered_dict_show.jl")
include("with_tracers.jl")
include("versioninfo.jl")
include("schedules.jl")
include("user_function_arguments.jl")

end # module


text:
module PatchMixtureKriging

    # using FFTW
    #
    # #using Distributed
    # #using SharedArrays
    #
    # #import Random
    # #import Printf
    #
    # import Interpolations
    #
    # import Convex
    # import SCS
    # #import SDPA
    #
    # import RWUtilities
    # import SignalTools
    #
    # import Optim
    # import VisualizationTools
    # import Colors
    # import RWRiemannianOptim

    using LinearAlgebra
    import Statistics
    using AbstractTrees

    include("./misc/declarations.jl")
    include("./misc/front_end.jl")

    include("./RKHS/RKHS.jl")
    include("./RKHS/kernel.jl")

    include("./RKHS/patchwork_RKHS.jl")

    include("./RKHS/querying.jl")

    include("./patchwork/partition.jl")
    include("./RKHS/mixtureGP.jl")

    include("./misc/utilities.jl")

    include("./patchwork/visualize_2D.jl")

    ### legacy and future clean-up/expansion.
    #include("./RKHS/ODEkernels.jl")
    ##include("./RKHS/interpolators.jl")
    #include("./RKHS/derivatives/SqExp_derivatives.jl")
    #include("./warp_map/Rieszanalysis.jl")
    #include("./optimization/sdp.jl")
    #include("./misc/utilities.jl")

    export  RKHSProblemType,
            fitRKHS!,
            query!,
            #fitnDdensity,
            #fitpdfviaSDP,
            #fitnDdensityRiemannian,
            constructkernelmatrix,
            evalkernel,
            evalquery,

            # mixGP.
            setuppartition,
            getpartitionlines!,
            organizetrainingsets,
            fetchhyperplanes,
            MixtureGPType,
            MixtureGPDebugType,
            fitmixtureGP!

end # module

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
module PlantingRunner
import Planting
import Plots

function __init__()
    Plots.pyplot(leg = false, grid = true, xticks = nothing, yticks = nothing)
end

type PlantingPlot
    plt::Nullable{Plots.Plot}
    active_list::Vector{Planting.AbstractPoint}
    final_list::Vector{Planting.AbstractPoint}
    radius::Float64
    poly::Nullable{Planting.AbstractPolygon}
    anim::Plots.Animation
end

PlantingPlot() = PlantingPlot(Nullable{Plots.Plot}(),[],[],0, Nullable{Planting.AbstractPolygon}(), Plots.Animation())

function drawCircle(plt::Plots.Plot,p::Planting.AbstractPoint, poly::Planting.AbstractPolygon, radius::Float64, args...; kwargs...)
    circlepoint = Plots.partialcircle(0,2π,1000,radius)
    (xs,ys) = Plots.unzip(circlepoint)

    xs += Planting.x(p)
    ys += Planting.y(p)

    ps = map(zip(xs,ys)) do x
        p = Planting.Point(x...)
        if ! (p in poly)
            p2 = [p...]
            closestLine = nothing
            minDist = Inf
            for idx in 1:(length(poly) - 1)

                distance = begin
                    l = norm(Planting.Line(poly[idx],poly[idx + 1]))
                    if l == 0
                        norm(Planting.Line(poly[idx],p))
                    else
                        v = [poly[idx]...]
                        w = [poly[idx + 1]...]

                        t = max(0,min(1,dot(p2 - v, w - v)/(l*l)))

                        proj = v + t * (w - v)
                        norm(Planting.Line(p,Planting.Point(proj...)))
                    end
                end

                if distance < minDist
                    minDist = distance
                    closestLine = idx
                end
            end

            p = begin
                l = norm(Planting.Line(poly[closestLine],poly[closestLine + 1]))
                if l == 0
                    poly[closestLine]
                else
                    v = [poly[closestLine]...]
                    w = [poly[closestLine + 1]...]

                    t = max(0,min(1,dot(p2 - v, w - v)/(l*l)))

                    proj = v + t * (w - v)
                    Planting.Point(proj...)
                end
            end
        end
        return (p...)
    end

    Plots.plot!(plt,Plots.Shape(Plots.unzip(ps)...),args...; kwargs...)
end

function init(plt::PlantingPlot, poly::Planting.AbstractPolygon, o::Planting.AbstractObject, active_list, final_list, grid)
    (lb,ub) = Planting.boundingbox(poly)

    width = Planting.x(ub) - Planting.x(lb)
    height = Planting.y(ub) - Planting.y(lb)

    midx = (Planting.x(lb) + Planting.x(ub)) / 2
    midy = (Planting.y(lb) + Planting.y(ub)) / 2

    xlim = (midx - ((width/2)*1.1),midx + ((width/2)*1.1))
    ylim = (midy - ((height/2)*1.1),midy + ((height/2)*1.1))

    plt.plt = Plots.plot(xlim = xlim, ylim = ylim, aspect_ratio = :equal, xticks = Planting.x(lb):grid.gridSize:Planting.x(ub),yticks = Planting.y(lb):grid.gridSize:Planting.y(ub))
    plt.active_list = active_list
    plt.final_list = final_list
    plt.radius = Planting.radius(o)
    plt.poly = poly

    xs = map(Planting.x,poly)
    ys = map(Planting.y,poly)
    Plots.plot!(get(plt.plt),Plots.Shape(xs[1:end],ys[1:end]), c = Plots.RGBA(1,1,1,0))
    Plots.gui()
    sleep(.01)
end

function test(plt::PlantingPlot, p::Planting.AbstractPoint, np::Planting.AbstractPoint, points::Vector, active_list, final_list, grid)
    px = Planting.x(p)
    py = Planting.y(p)
    poly = get(plt.poly)
    (lb,ub) = Planting.boundingbox(get(plt.poly))

    width = Planting.x(ub) - Planting.x(lb)
    height = Planting.y(ub) - Planting.y(lb)

    midx = (Planting.x(lb) + Planting.x(ub)) / 2
    midy = (Planting.y(lb) + Planting.y(ub)) / 2

    xlim = (midx - ((width/2)*1.1),midx + ((width/2)*1.1))
    ylim = (midy - ((height/2)*1.1),midy + ((height/2)*1.1))
    xticks = Planting.x(lb):grid.gridSize:Planting.x(ub)
    yticks = Planting.y(lb):grid.gridSize:Planting.y(ub)
    plt.plt = Plots.plot(xlim = xlim, ylim = ylim, aspect_ratio = :equal, xticks = (xticks,repmat([""],length(xticks))), yticks = (yticks,repmat([""],length(yticks))))

    xs = map(Planting.x,poly)
    ys = map(Planting.y,poly)
    Plots.plot!(get(plt.plt),Plots.Shape(xs[1:end],ys[1:end]), c = Plots.RGBA(1,1,1,0))

    # active_idx = findin(final_list,active_list)
    for point in points
        norm(Planting.Line(np,point)) <= plt.radius && Plots.plot!(get(plt.plt),map(Planting.x,[np,point]),map(Planting.y,[np,point]), linecolor = :black, marker = :none)
    end

    for plant in final_list
        circlepoint = Plots.partialcircle(0,2π,1000,plt.radius)
        (xs,ys) = Plots.unzip(circlepoint)

        xs += Planting.x(plant)
        ys += Planting.y(plant)

        drawCircle(get(plt.plt),plant,get(plt.poly), plt.radius; c = Plots.RGBA(.5,.5,.5,.1), linewidth = 0)

        if plant in active_list
            Plots.scatter!(get(plt.plt),[Planting.x(plant)],[Planting.y(plant)], markercolor = :red)
        else
            Plots.scatter!(get(plt.plt),[Planting.x(plant)],[Planting.y(plant)], markercolor = :black)
        end
    end

    drawCircle(get(plt.plt),p,get(plt.poly), plt.radius; c = Plots.RGBA(.5,.5,.5,0))
    drawCircle(get(plt.plt),p,get(plt.poly), 2 * plt.radius; c = Plots.RGBA(.5,.5,.5,.2))

    Plots.scatter!(get(plt.plt),[Planting.x(np)],[Planting.y(np)], markercolor = :white)

    Plots.gui()
    sleep(.01)
    # sleep(.1)

end

function step(plt::PlantingPlot, active_list, final_list, grid)
    px = Planting.x(p)
    py = Planting.y(p)
    poly = get(plt.poly)
    (lb,ub) = Planting.boundingbox(get(plt.poly))

    width = Planting.x(ub) - Planting.x(lb)
    height = Planting.y(ub) - Planting.y(lb)

    midx = (Planting.x(lb) + Planting.x(ub)) / 2
    midy = (Planting.y(lb) + Planting.y(ub)) / 2

    xlim = (midx - ((width/2)*1.1),midx + ((width/2)*1.1))
    ylim = (midy - ((height/2)*1.1),midy + ((height/2)*1.1))
    xticks = Planting.x(lb):grid.gridSize:Planting.x(ub)
    yticks = Planting.y(lb):grid.gridSize:Planting.y(ub)
    plt.plt = Plots.plot(xlim = xlim, ylim = ylim, aspect_ratio = :equal, xticks = (xticks,repmat([""],length(xticks))), yticks = (yticks,repmat([""],length(yticks))))

    xs = map(Planting.x,poly)
    ys = map(Planting.y,poly)
    Plots.plot!(get(plt.plt),Plots.Shape(xs[1:end],ys[1:end]), c = Plots.RGBA(1,1,1,0))

    for plant in final_list
        circlepoint = Plots.partialcircle(0,2π,1000,plt.radius)
        (xs,ys) = Plots.unzip(circlepoint)

        xs += Planting.x(plant)
        ys += Planting.y(plant)

        drawCircle(get(plt.plt),plant,get(plt.poly), plt.radius; c = Plots.RGBA(.5,.5,.5,.1), linewidth = 0)

        if plant in active_list
            Plots.scatter!(get(plt.plt),[Planting.x(plant)],[Planting.y(plant)], markercolor = :red)
        else
            Plots.scatter!(get(plt.plt),[Planting.x(plant)],[Planting.y(plant)], markercolor = :black)
        end
    end

    Plots.gui()
    sleep(.01)
end

function main()
    plt = PlantingPlot()

    callbacks = Dict{String,Function}(
        "init" => (poly,o,active_list, final_list, grid) -> init(plt,poly,o,active_list,final_list,grid),
        "test" => (p, np, points, active_list, final_list, grid) -> test(plt,p,np,points, active_list, final_list, grid)
    )

    poly = Planting.Polygon([ Planting.Point(0.0,0.0),Planting.Point(0.0,1.0),Planting.Point(1.0,1.0), Planting.Point(1.0,0.0), Planting.Point(.5,-.5)])

    obj = Planting.SimpleObject(.2)

    plants = Planting.plant(Planting.FastPoissonDisk(30), poly, obj; seed = 0, callback = callbacks)

end

end


text:
using Qiskit
using Test

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


text:
module ChirpsController

using App, SearchLight, Chirps, Validation
@dependencies

# before_hook() = println("I run before the action")
# after_hook() = println("I run after the action")
# rescue_hook(ex::Exception) = redirect_to("/omg")

const CHIRPS_PER_PAGE = 20

function index()
  total_chirps = SearchLight.count(Chirp)
  chirps = SearchLight.find(Chirp, SQLQuery(limit = CHIRPS_PER_PAGE, offset = Int(@params(:page, 0)) * CHIRPS_PER_PAGE, order = "created_at DESC"))
  respond_with_html(:chirps, :index, chirps = chirps, chirps_per_page = CHIRPS_PER_PAGE, total_chirps = total_chirps)
end

function new(chirp = Chirp(content = ""))
  respond_with_html(:chirps, :new, chirp = chirp)
end

function create()
  chirp = Chirp(content = @params(:content))
  if Validation.validate!(chirp) && SearchLight.save(chirp)
    flash("Your chirp was saved")
    redirect_to(:get_chirps)
  else
    new(chirp)
  end
end

end

max_stars_repo_name: essenciary/chirper

max_issues_repo_name: essenciary/chirper

max_forks_repo_name: essenciary/chirper


text:
const ALL_COLS = [:name, :charitytype, :ein, :fisyr, :fileyear, :invrealized,
  :invrealizeddirect, :totrev, :totexp,
  :invrealizedindirect, :name, :netassets,
  :pfexpbook, :pfexpcharity, :totinv, :category,
  :returnamt, :returncheckamt,:netincome, :naics]

const CATEGORY_CONSOLIDATED_THRESHOLD = 100_000

#=const :Arts, :Education, :Environmental, :AnimalWelfare, :Health, :MentalHealth,
:MedicalDiseases, :MedicalResearch, :CrimeAndLegal, :Employment, :FoodAndAgriculture,
:Housing, :PublicSafety, :Recreation, :YouthDevelopment, :HumanServices, :ForeignAffairs,
:SocialAction, :CommunityImprovement, :GrantmakingFoundations, :ScienceandTechResearch,
:SocialScienceResearch, :SocietyBenefitAndMultipurpose, :ReligionRelated,
:MutualBenefitOrg, :Unknown, :Other=#

function processReturns!(data::NCCSData)

  local N::Int = size(data.df,1) #length of dataframe
  local assetCutoff::Float64 = ASSET_CUTOFF #minimum allowed maximum assets
  local largeDifferenceCutoff::Float64 = LARGE_DIFFERENCE_CUTOFF #max discrepenacy between methods

  println("data.df N1: $N")
  println("data.df pf N1: $(sum(data.df[:charitytype] .== :pf))")
  #filter out organizations which are too small
  data.df[:tokeep] = trues(size(data.df,1))
  dfByein::GroupedDataFrame = groupby(data.df, :ein)
  for subdf::SubDataFrame ∈ dfByein
    subN::Int = size(subdf,1)
    maxAssets::Float64 = maximum(collect(skipmissing(subdf[:netassets])))
    if (maxAssets < assetCutoff)
      subdf[:tokeep] = falses(subN)
    end

    if subN ≠ length(unique(subdf[:fisyr])) # check if we have duplicate filings
      deleteOrg::Bool = false
      subdfByfisyr = groupby(subdf, :fisyr)
      for ssubdf::SubDataFrame ∈ subdfByfisyr
        nrows::Int = size(ssubdf, 1)
        if nrows > 1 # if we are at the org-year with the dup
          fileYears::Vector{Int} = unique(ssubdf[:fileyear])
          if length(fileYears) == nrows #if we can reconcile by fileYear
            maxFileYear::Int = maximum(fileYears) # keep the newest version
            ssubdf[ssubdf[:fileyear] .≠ maxFileYear, :tokeep] = false
          else
            deleteOrg = true #If we cannot correct, throw out the org
          end
        end
      end
      if deleteOrg #delete the org if we need to
        subdf[:tokeep] = false
      end
    end
  end

  #delete the invalid rows
  data.df = data.df[data.df[:tokeep] .== true, :]
  N = size(data.df,1)

  println("data.df N2: $N")

  #println("N2-pf: $(sum(data.df[:charitytype] .== :pf)) \nN2-pf-Return Amt: $(
  #  sum((data.df[:charitytype] .== :pf) .& ((!ismissing).(data.df[:returnamt]))))")

  #data.df[:returnamt] = Vector{Union{Float64,Missing}}(missings(N)) #don't need this since its done during mapping
  data.df[:return] = Vector{Union{Float64,Missing}}(missings(N))

  #data.df[:returncheckamt] = Vector{Union{Float64,Missing}}(missings(N))
  data.df[:returncheck] = Vector{Union{Float64,Missing}}(missings(N))

  data.df[:largeDifference] = Vector{Union{Bool,Missing}}(missings(N))

  sort!(data.df, [:ein, :fisyr])
  #ctr::Int = 0
  dfByein = groupby(data.df, :ein)
  for subdf::SubDataFrame ∈ dfByein
    subN = size(subdf,1)

    if (subdf[1,:charitytype] == :co) || (subdf[1,:charitytype] == :pc)

      for i ∈ 2:subN
        if subdf[i-1, :fisyr] == subdf[i, :fisyr] - 1

          #NOTE: returns = unrealized + inv income s.t. unrealized = Δassets - netincome
          returnCandidateAmt::MFloat64 = (subdf[i, :netassets] .- subdf[i-1, :netassets]
            ) .+ subdf[i, :invrealizeddirect] .- subdf[i, :netincome]
          returnCandidate::MFloat64 = returnCandidateAmt / subdf[i-1, :netassets]

          returnCheckCandidateAmt::MFloat64 = (subdf[i, :netassets] .- subdf[i-1, :netassets]
            ) .+ subdf[i, :invrealizedindirect] .- subdf[i, :netincome]
          returnCheckCandidate::MFloat64 = returnCheckCandidateAmt / subdf[i-1, :netassets]

          if ((!ismissing(returnCandidate)) && returnCandidate < 2. && returnCandidate > -0.9 &&
              abs(returnCandidate - returnCheckCandidate) ≤ largeDifferenceCutoff)# && abs(subdf[i, :totinv]) ≥ 1. #make sure we have valid numbers here
            subdf[i, :returnamt] = returnCandidateAmt
            subdf[i, :return] = returnCandidate

            subdf[i, :returncheckamt] = returnCheckCandidateAmt
            subdf[i, :returncheck] = returnCheckCandidate
          end
        end
      end
    elseif subdf[1,:charitytype] == :pf
      for i ∈ 2:subN
        if subdf[i-1, :fisyr] == subdf[i, :fisyr] - 1

          #NOTE: returns = unrealized + inv income s.t. unrealized = Δassets - netincome
          returnCandidateAmt::MFloat64 = subdf[i, :returnamt]
          returnCandidate::MFloat64 = returnCandidateAmt / subdf[i-1, :netassets]

          returnCheckCandidateAmt::MFloat64 = subdf[i, :returncheckamt]
          returnCheckCandidate::MFloat64 = returnCheckCandidateAmt / subdf[i-1, :netassets]

          if ( (!ismissing(returnCandidate)) && (returnCandidate < 2.) && (returnCandidate > -0.9) &&
              (!ismissing(returnCheckCandidate)) &&
              abs(returnCandidate - returnCheckCandidate) ≤ largeDifferenceCutoff)# && abs(subdf[i, :totinv]) ≥ 1. #make sure we have valid numbers here
            subdf[i, :returnamt] = returnCandidateAmt
            subdf[i, :return] = returnCandidate

            subdf[i, :returncheckamt] = returnCheckCandidateAmt
            subdf[i, :returncheck] = returnCheckCandidate
          #=elseif ctr < 50
            ctr+=1
            println("###\nreturn: $(returnCandidate)")
            println("returnchk: $(returnCheckCandidate)")
            println("returnamt: $(returnCandidateAmt)")
            println("returnchk: $((returnCheckCandidateAmt))")=#
          end
        end
      end
    end
  end

  sort!(data.df, [:ein, :fisyr])

  return data
end

function processWealth!(data::NCCSData)

  data.df[:lnetassets] = (log).(data.df[:netassets])
  data.df[:ladjnetassets] = (log).(data.df[:adjnetassets])

  return data
end

#data processing related to the categorization
function processCategories!(data::NCCSData;
  minPointsPerCategory::Int = MIN_POINTS_PER_CATEGORY,
  categoryConsolidatedThreshold::Int = CATEGORY_CONSOLIDATED_THRESHOLD)::Nothing

  data.df[:categoryname] = (s::MSymbol->((ismissing(s) || (!haskey(category, s))) ?
    missing : category[s])).(Vector{MSymbol}(data.df[:category]))

  data.df[:tokeep] .= true
  data.df[:categoryfiltered] = deepcopy(data.df[:category])
  data.df[:categorynamefiltered] = deepcopy(data.df[:categoryname])

  dfBycategoryfisyr::GroupedDataFrame = groupby(data.df, [:category, :fisyr])
  for subdf::SubDataFrame ∈ dfBycategoryfisyr
    if size(subdf,1) < minPointsPerCategory
      subdf[:categoryfiltered] = :ZZ
      subdf[:categorynamefiltered] = category[:ZZ]
    end
  end

  data.df[:categoryconsolidated] = deepcopy(data.df[:categoryfiltered])
  data.df[:categorynameconsolidated] = deepcopy(data.df[:categorynamefiltered])

  dfBycategorynamefiltered::GroupedDataFrame = groupby(data.df, :categorynamefiltered)
  for subdf::SubDataFrame ∈ dfBycategorynamefiltered
    if size(subdf,1) < categoryConsolidatedThreshold
      subdf[:categoryconsolidated] = :ZZ
      subdf[:categorynameconsolidated] = category[:ZZ]
    end
  end

  data.df[:naics2] = (s::MString->(!ismissing(s) ? Symbol(s[1:2]) : missing)).(data.df[:naics])
  data.df[:naics2name] = (s::MSymbol->(naicsCodes[s])).(data.df[:naics2])

  return nothing
end

#adjusts for inflation
#currently just does wealth, eventually can use this for return fields
function adjustInflation!(data::NCCSData; fredPath::String = FRED_PATH,
  fredInflationFile::String = FRED_INFLATION_FILE, inflationSeries::Symbol = :CPI,
  baseYear::Int = BASE_YEAR, fredDateFormat::String = FRED_DATE_FORMAT,
  adjustWealthCutoff::Bool = true)::Nothing

  #read in the inflation info
  freddf::DataFrame = CSV.read("$fredPath\\$fredInflationFile.csv")

  #make a dictionary
  cpiIndex::Dict =
    Dict(Dates.year(Date(freddf[i,:DATE], fredDateFormat))=>freddf[i, inflationSeries]
      for i ∈ 1:size(freddf,1))

  #get vectors for performance purpuses
  nrows::Int = size(data.df, 1)
  baseValue::Float64 = cpiIndex[baseYear]
  netassets::Vector{Int} = data.df[:netassets]
  adjnetassets = Vector{MFloat64}(undef, nrows)
  fisyr::Vector{Int} = data.df[:fisyr]
  @inbounds @simd for i ∈ 1:nrows
    if haskey(cpiIndex, fisyr[i])
      adjnetassets[i] = netassets[i] / (cpiIndex[fisyr[i]] / baseValue)
    end
  end

  data.df[:adjnetassets] = adjnetassets
  if adjustWealthCutoff #this is to account for the cutoff inflating forward
    inflatedCutoff::Float64 = ADJ_NET_ASSETS_CUTOFF * baseValue/cpiIndex[minimum(data.df[:fisyr])]
    println("Inflated cutoff of $inflatedCutoff")
    data.df = data.df[
      (f::MFloat64-> (!ismissing(f)) && (f ≥ inflatedCutoff)).(data.df[:adjnetassets]), :]
  end

  return nothing
end

#does an initial filter to get down the file size
function processAll!(data::NCCSData; allcols::Vector{Symbol} = ALL_COLS,
  adjNetAssetsCutoff::Float64 = ADJ_NET_ASSETS_CUTOFF)

  data = removeFields!(data, setdiff(names(data.df), allcols))
  data.df[:rows] = collect(1:size(data.df,1))
  println("N0: $(size(data.df,1))")
  deleterows!(data.df, data.df[(!).(completecases(data.df[[:fisyr, :netassets, :totrev, :totexp]])),:rows])

  adjustInflation!(data)
  processReturns!(data)
  processWealth!(data)
  processCategories!(data)

  println("N3: $(size(data.df,1)) N3-Return: $(sum((!ismissing).(data.df[:return])))")
  println("N3-pf: $(sum(data.df[:charitytype] .== :pf)) \nN3-pf-Return: $(
    sum((data.df[:charitytype] .== :pf) .& ((!ismissing).(data.df[:return]))))\nN3 check return: $(
    sum((data.df[:charitytype] .== :pf) .& ((!ismissing).(data.df[:returncheckamt]))))\nN3 net income: $(
    sum((data.df[:charitytype] .== :pf) .& ((!ismissing).(data.df[:netincome]))))")

  data.df = data.df[((!ismissing).(data.df[:adjnetassets])) .&
    (data.df[:adjnetassets] .≥ adjNetAssetsCutoff), :]

  println("N4: $(size(data.df,1)) N4-Return: $(sum((!ismissing).(data.df[:return])))")

  return data
end

function constructPrimaryOutput(; refreshFilter::Bool = true,
    primaryOutputName::String = PRIMARY_OUTPUT_NAME,
    writePrimaryCSV::Bool = false,
    dataPath::String = DATA_PATH)::NCCSData

  local data::NCCSData
  local outPath = "$dataPath\\$(primaryOutputName)_long.csv"
  if refreshFilter
    data = deserializeNCCSData()
    data = processAll!(data)
    serializeNCCSData(data, dataName = primaryOutputName)
  else
    data = deserializeNCCSData(dataName = primaryOutputName)
  end

  #for debugging
  #uCSV.write(outPath, data.df[400_000:450_000,:])

  writePrimaryCSV && uCSV.write(outPath, data.df) #write to a csv if desired

  return data

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using AbstractAlgebra
using Nemo
using SymbolicIntegration
SI = SymbolicIntegration

using Test

@testset "Chapter 5" begin
    @info "HermiteReduce, example 5.3.1, p. 140"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    D0 = BasicDerivation(QQx)
    kt, t = PolynomialRing(k, :t)
    D = SI.ExtensionDerivation(kt, D0, 1+t^2)

    f = (x-t)//t^2
    g, h, r = SI.HermiteReduce(f, D)
    expected = (-x*1//t, 0, -x)
    @test (g, h, r) == expected

    @info "PolynomialReduce, example 5.4.1, p. 141"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    D0 = BasicDerivation(QQx)
    kt, t = PolynomialRing(k, :t)
    D = SI.ExtensionDerivation(kt, D0, 1+t^2)

    p = 1 + x*t + t^2
    q, r = SI.PolynomialReduce(p, D)
    expected = (t, x*t)
    @test (q, r) == expected

#=
    @info "ResidueReduce, ConstantPart, example 5.6.3, p. 151"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    D0 = BasicDerivation(QQx)
    kt, t = PolynomialRing(k, :t)
    D = ExtensionDerivation(kt, D0, 1//x+0*t)

    f = (2*t^2-t-x^2)//(t^3-x^2*t)
    ss, Ss, ρ = SI.ResidueReduce(f, D)
    @test ρ == 0
    αs, gs, ss1, Ss1 = SI.ConstantPart(ss, Ss, D)
    p = sortperm(αs)
    αs = αs[p]
    gs = gs[p]
    @test αs == [-1//2, 1//2]
    @test gs == [-(2*x^2+3*x+1)*(t-x)//2, -(2*x^2-3*x+1)*(t+x)//2]
    Dg = sum([αs[i]*D(gs[i])//gs[i] for i=1:length(αs) ])
    @test f-Dg == 1//t - (6*x^2-3)//(4*x^4-5*x^2+1)
=#

    @info "IntegratePrimitivePolynomial, example 5.8.1, p. 158"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k0 = FractionField(QQx)
    D0 = BasicDerivation(QQx)
    k0t0, t0 = PolynomialRing(k0, :t0)
    k = FractionField(k0t0)
    D1 = ExtensionDerivation(k0t0, D0, 1//x+0*t0)
    kt, t = PolynomialRing(k, :t)
    D = ExtensionDerivation(kt, D1, 1//t0+0*t)

    p = (t0+1//t0)*t - x*1//t0
    q, ρ = SI.IntegratePrimitivePolynomial(p, D)
    @test ρ == 1
    @test q == t^2//2 + (x*t0-x)*t
    @test p - D(q) == -x

    @info "IntegrateHyperexponentialPolynomial, example 5.9.1, p. 162"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k0 = FractionField(QQx)
    D0 = BasicDerivation(QQx)
    k0t0, t0 = PolynomialRing(k0, :t0)
    k = FractionField(k0t0)
    D1 = ExtensionDerivation(k0t0, D0, 1+t0^2)
    kt, t = PolynomialRing(k, :t)
    D = ExtensionDerivation(kt, D1, (1+t0^2)*t)

    p = (t0^3+(x+1)*t0^2+t0+x+2)*t + 1//(x^2+1)
    # p must be in k[t, t⁻¹] => must be passed as a rational function
    q, ρ = SI.IntegrateHyperexponentialPolynomial(p//1, D)
    @test ρ == 1
    @test q == (t0+x)*t
    @test p - D(q) == 1//(x^2+1)

    @info "IntegrateHypertangentPolynomial, example 5.10.1, p. 167"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    kt, t = PolynomialRing(k, :t)
    D = ExtensionDerivation(kt,  BasicDerivation(QQx), 1+t^2)

    p = t^2 + x*t + 1
    q, c = SI.IntegrateHypertangentPolynomial(p, D)
    expected = (t, x//2)
    @test (q, c) == expected

    @info "IntegrateHypertangentReduced, example 5.10.2, p. 169"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    kt, t = PolynomialRing(k, :t)
    D = ExtensionDerivation(kt,  BasicDerivation(QQx), 1//2*(1+t^2))

    p = (2*t*1//x)//(t^2+1)
    q, ρ = SI.IntegrateHypertangentReduced(p, D)
    @test ρ == 0

    @info "IntegrateHypertangentReduced, example 5.10.3, p. 170"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    kt, t = PolynomialRing(k, :t)
    D = ExtensionDerivation(kt,  BasicDerivation(QQx), 1+t^2)

    p = (t^5 + t^3 + x^2*t + 1)//(t^2 + 1)^3
    q, ρ = SI.IntegrateHypertangentReduced(p, D)
    q_ex = ((1 + 1//3*x)*t - (x^2 - 1//18))//(6*(t^2 + 1)^3) +
           (5*(1 + 1//3*x)*t + 77//12)//(24*(t^2 + 1)^2) + (5*(1 + 1//3*x)*t -43//6)//(16*(t^2 + 1))
    @test ρ == 1
    @test q == q_ex
    @test p - D(q) == 5//16*(1 + x//3)
end

@testset "Chapter 6" begin
    @info "RdeNormalDenominator, example 6.1.1, p. 186"
    QQt, t = PolynomialRing(Nemo.QQ, :t)
    k = FractionField(QQt)
    D = BasicDerivation(QQt)

    f = 1 + 0//t
    g = 1//t
    (a, b, c, h, ρ) = SI.RdeNormalDenominator(f, g, D)
    @test ρ == 0

    @info "RdeNormalDenominator, example 6.1.2, p. 186"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    kt, t = PolynomialRing(k, :t)
    D = ExtensionDerivation(kt, BasicDerivation(k), 1+t^2)

    f = t^2+1+0//t
    g = 1//t^2
    (a, b, c, h, ρ) = SI.RdeNormalDenominator(f, g, D)
    expected = (t, (t-1)*(t^2+1), 1, t, 1)
    @test (a, b, c, h, ρ) == expected

    @info "RdeSpeciaDenominatorExp, example 6.2.1, p. 190"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    kt, t = PolynomialRing(k, :t)
    D = ExtensionDerivation(kt, BasicDerivation(k), t)

    a = t^2+2*x*t+x^2
    b = (1+1//x^2)*t^2 + (2*x-1+2//x)*t + x^2 +0//t
    c = t*1//x^2-1+2*1//x +0//t
    a, b, c, h = SI.RdeSpecialDenomExp(a, b, c, D)
    expected = (t^2+2*x*t+x^2, t^2*1//x^2+(2//x-1)*t, t^2*1//x^2+(2//x-1)*t, t)
    @test (a, b, c, h) == expected

    @info "RdeSpeciaDenominatorTan, example 6.2.2, p. 192"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    kt, t = PolynomialRing(k, :t)
    D = ExtensionDerivation(kt, BasicDerivation(k), 1 + t^2)

    a = t
    b = (t - 1)*(t^2 + 1) + 0//t
    c = 1 + 0//t
    a1, b1, c1, h = SI.RdeSpecialDenomTan(a, b, c, D)
    @test (a1, b1, c1, h) == (a, b, c, 1)

    @info "RdeBoundDegreePrim, example 6.3.1, p. 198"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k0 = FractionField(QQx)
    k0t0, t0 = PolynomialRing(k0, :t0)
    D0 = ExtensionDerivation(k0t0, BasicDerivation(k0), 1//x^2 + 0*t0)
    k = FractionField(k0t0)
    kt, t = PolynomialRing(k, :t)
    D = ExtensionDerivation(kt, D0, 1//x + 0*t)

    a = t^2
    b = -(1//x^2*t^2+1//(x+0*t0))
    c = (2*x-1)*t^4 + ((t0+x)//(x+0*t0))*t^3 - ((t0+4*x^2)//(2*x+0*t0))*t^2 + x*t
    n = n = SI.RdeBoundDegreePrim(a, b, c, D)
    @test n==3

    @info "RdeBoundDegreeBase, example 6.3.2, p. 199"
    QQt, t = PolynomialRing(Nemo.QQ, :t)
    D = BasicDerivation(QQt)

    a = 1 +0*t
    b = -2*t
    c = 1 +0*t
    n = SI.RdeBoundDegreeBase(a, b, c)
    @test n==0

    @info "RdeBoundDegreeExp, example 6.3.3, p. 201"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    kt, t = PolynomialRing(k, :t)
    D = ExtensionDerivation(kt, BasicDerivation(k), t)

    # a,b,c ... result of Ex 6.2.1
    a = t^2+2*x*t+x^2
    b = t^2*1//x^2+(2//x-1)*t
    c = b
    n = SI.RdeBoundDegreeExp(a, b, c, D)
    @test n==0

    @info "RdeBoundDegreeNonLinear, example 6.3.4, p. 202"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    kt, t = PolynomialRing(k, :t)
    D = ExtensionDerivation(kt, BasicDerivation(k), 1+t^2)

    a = t
    b = (t-1)*(t^2+1)
    c = 1 + 0*t
    n = SI.RdeBoundDegreeNonLinear(a, b, c, D)
    @test n==0

    @info "SPDE, example 6.4.1, p. 203"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    kt, t = PolynomialRing(k, :t)
    D = ExtensionDerivation(kt, BasicDerivation(k), 1+t^2)

    a = t
    b = (t-1)*(t^2+1)
    c = 1 + 0*t
    n = 0
    (b, c, m, α, β, ρ) =  SI.SPDE(a, b, c, D, n)
    @test ρ==0

    @info "SPDE, example 6.4.2, p. 203"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    kt, t = PolynomialRing(k, :t)

    D = ExtensionDerivation(kt, BasicDerivation(k), t)
    # a,b,c ... result of example 6.2.1
    a = t^2+2*x*t+x^2
    b = t^2*1//x^2+(2//x-1)*t
    c = b
    n = 0
    (b, c, m, α, β, ρ) =  SI.SPDE(a, b, c, D, n)
    expected = (0, 0, 0, 0, 1, 1)
    @test (b, c, m, α, β, ρ) == expected

    @info "SPDE, example 6.4.3, p. 204"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k0 = FractionField(QQx)
    k0t0, t0 = PolynomialRing(k0, :t0)
    k = FractionField(k0t0)
    kt, t = PolynomialRing(k, :t)
    D0 = ExtensionDerivation(k0t0, BasicDerivation(k0), t0*1//x^2)
    D = ExtensionDerivation(kt, D0, 1//x+zero(t))
    a = t^2
    b = -t^2*1//x^2-1//x
    c = (2*x-1)*t^4 + (t0+x)*1//x*t^3 - (t0+4*x^2)*1//(2*x)*t^2 + x*t
    n = 3
    (b, c, m, α, β, ρ) =  SI.SPDE(a, b, c, D, n)
    expected = (0, 0, 0, 0, (x^2+t0//2)*t^2-x^2*t, 1)
    @test (b, c, m, α, β, ρ) == expected

    @info "SPDE, example 6.4.4, p. 205"
    QQt, t = PolynomialRing(Nemo.QQ, :t)
    D = BasicDerivation(QQt)

    a = t^2+t+1
    b = -2*t-1
    c = 1//2*t^5+3//4*t^4+t^3-t^2+1
    n = 7 # arbitrary !
    (b, c, m, α, β, ρ) =  SI.SPDE(a, b, c, D, n)
    expected = (0, 1//2*t-1//4, n-2, t^2+t+1, 5//4*t, 1)
    @test (b, c, m, α, β, ρ) == expected

    @info "PolyRischDENoCancel1, example 6.5.1, p. 208"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    kt, t = PolynomialRing(k, :t)
    D = ExtensionDerivation(kt, BasicDerivation(k), 1+t^2)

    b = t^2+1
    c = t^3+(x+1)*t^2+t+x+2
    q, ρ = SI.PolyRischDENoCancel1(b, c, D)
    expected = (t+x, 1)
    @test (q, ρ) == expected

    @info "PolyRischDENoCancel2, example 6.5.2, p. 209"
    kt, t = PolynomialRing(Nemo.QQ, :t)
    D = BasicDerivation(kt)

    b = zero(kt)
    c = 1//2*t-1//4
    h, b, c, ρ  = SI.PolyRischDENoCancel2(b, c, D)
    @test ρ == 1
    @test h == 1//4*t^2-1//4*t

    @info "PolyRischDENoCancel3, example 6.5.3, p. 211"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    kt, t = PolynomialRing(k, :t)
    D = ExtensionDerivation(kt, BasicDerivation(k), 1 + t^2)

    b = 1-t
    c = t^3+t^2-2*x*t-2*x
    h, m, c, ρ   = SI.PolyRischDENoCancel3(b, c, D)
    expected = (t^2,  1, -2*(x + 1)*t - 2*x, 2)
    @test (h, m, c, ρ) == expected

    @info "PolyRischDECancelTan, example 6.6.1, p. 216"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    kt, t = PolynomialRing(k, :t)
    D = ExtensionDerivation(kt, BasicDerivation(k), 1 + t^2)

    b0 = 1 + 0//x
    c = -2*(x + 1)*t - 2*x
    n = 1
    q, ρ = SI.PolyRischDECancelTan(b0, c, D, n)
    expected = (-2*x*t, 1)
    @test  (q, ρ) == expected
end

@testset "Chapter 7" begin
    @info "ParametricLogarithmicDerivative, example 7.3.1, p. 254"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    kt, t = PolynomialRing(k, :t)
    D = SI.ExtensionDerivation(kt, SI.BasicDerivation(k), 1//x+0*t)

    w = -1//(x*t^2)
    f = (5*t^2+t-6)//(2*x*t^2)
    n, m, u, ρ = SI.ParametricLogarithmicDerivative(f, w, D)
    expected = (2, 6, x^5*t, 1)
    @test (n, m, u, ρ) == expected
end

@testset "Chapter 8" begin
    @info "CoupledDECancelTan, example 8.4.1, p. 265"
    QQx, x = PolynomialRing(Nemo.QQ, :x)
    k = FractionField(QQx)
    kt, t = PolynomialRing(k, :t)
    D = SI.ExtensionDerivation(kt, SI.BasicDerivation(k), 1 + t^2)

    b0 = 0//x
    b2 = 4*x + 0//x
    c1 = -t^2 + 2*t - 8*x^2 + 1
    c2 = 2*(1 - 2*x) + 0*t
    n = 2
    q1, q2, ρ = SI.CoupledDECancelTan(b0, b2,  c1, c2, D, n)
    expected = (t - 1, 2*x, 1)
    @test (q1, q2, ρ) == expected
end


text:
push!(LOAD_PATH, "../src/")

using Documenter
using Arrhenius

makedocs(
    format=Documenter.HTML(),
    sitename="Arrhenius.jl",
    modules=[Arrhenius],
    pages=[
        ##############################################
        ## MAKE SURE TO SYNC WITH docs/src/index.md ##
        ##############################################
        "Basics" => [
            "index.md",
            "install.md",
            "get_started.md",
            "concepts.md"
           ],
        "tutorial.md",
        "faq.md",
        "api.md"
    ],
)

deploydocs(;
    repo="github.com/DENG-MIT/Arrhenius.jl",
    devbranch="main"
)