text:
@testset "PIS" begin
    mesh = buildSquare(0,1,0,1,2,2)

    sys = System{Float64}(mesh)

    dog = addVariable!(sys, "dog")
    cat = addVariable!(sys, "cat")

    initialize!(sys)

    pis = PetscImplicitSolver(sys)

    # Make sure it's not iniitialized by default
    @test !pis.initialized

    initialize!(pis)

    # Now it should be
    @test pis.initialized

    n_dofs = sys.n_dofs

    @test size(pis.mat) == (n_dofs, n_dofs)
    @test length(pis.rhs) == n_dofs
    @test length(pis.solution) == n_dofs

    # Put something in the matrix so we can do a solve
    # Just make it an identity matrix
    # Note: we actually want to fill the matrix with identity values
    #   not replace it with an identity matrix
    for i in 1:n_dofs
        pis.mat[i,i] = 1
    end

    for i in 1:n_dofs
        pis.rhs[i] = i
    end

    MiniPETSc.assemble!(pis.mat)
    MiniPETSc.assemble!(pis.rhs)

    # Now solve:
    solve!(pis, assemble=false)

    # Should be the case!
    for i in 1:length(pis.solution)
        @test abs(pis.solution[i] - pis.rhs[i]) < 1e-9
    end
end

function testFullSolvePIS(diffusion_system)
    u = addVariable!(diffusion_system, "u")

    diffusion_kernel = Diffusion(u)
    addKernel!(diffusion_system, diffusion_kernel)

    left_boundary = DirichletBC(u, [4], 0.0)
    addBC!(diffusion_system, left_boundary)

    right_boundary = DirichletBC(u, [2], 1.0)
    addBC!(diffusion_system, right_boundary)

    initialize!(diffusion_system)

    solver = PetscImplicitSolver(diffusion_system)

    solve!(solver)

    solution = serializeToZero(solver.solution)

    for i in 1:length(solution)
        @test abs(solver.solution[i] - [0.0, 0.5, 0.5, 0.0, 1.0, 1.0, 0.5, 0.0, 1.0][i]) < 1e-9
    end
end

# Test a 'Full Solve' using PIS
@testset "PISFS" begin
    mesh = buildSquare(0, 1, 0, 1, 2, 2)

    diffusion_system = System{Float64}(mesh)

    testFullSolvePIS(diffusion_system)
end

# Test a 'Full Solve' using PIS and Automatic Differentiation
@testset "PISFD" begin
    mesh = buildSquare(0, 1, 0, 1, 2, 2)

    diffusion_system = System{Dual{4, Float64}}(mesh)

    testFullSolvePIS(diffusion_system)
end


text:
export measurementSystemMatrix, SystemMatrixRobotMeas
export init, stop, isStarted

mutable struct SystemMatrixRobotMeas <: MeasObj
  task::Union{Task,Nothing}
  consumer::Union{Task, Nothing}
  producerFinalizer::Union{Task, Nothing}
  scanner::Union{MPIScanner,Nothing}
  store::DatasetStore
  params::Dict
  positions::GridPositions
  currPos::Int
  stopped::Bool
  signals::Array{Float32,4}
  currentSignal::Array{Float32,4}
  waitTime::Float64
  measIsBGPos::Vector{Bool}
  posToIdx::Vector{Int64}
  measIsBGFrame::Vector{Bool}
  temperatures::Matrix{Float64}
  prepared::Bool
end

function SystemMatrixRobotMeas(scanner, store)
  return SystemMatrixRobotMeas(
    nothing,
    nothing,
    nothing,
    scanner,
    store,
    Dict{String,Any}(),
    RegularGridPositions([1,1,1],[0.0,0.0,0.0],[0.0,0.0,0.0]),
    1, false, Array{Float32,4}(undef,0,0,0,0), Array{Float32,4}(undef,0,0,0,0),
    0.0, Vector{Bool}(undef,0), Vector{Int64}(undef,0), Vector{Bool}(undef,0),
    Matrix{Float64}(undef,0,0),
    false
  )
end

function cleanup(sysObj::SystemMatrixRobotMeas)
  filename = "/tmp/sysObj.toml"
  filenameSignals = "/tmp/sysObj.bin"
  rm(filename, force=true)
  rm(filenameSignals, force=true)
end

function store(sysObj::SystemMatrixRobotMeas)
  filename = "/tmp/sysObj.toml"
  rm(filename, force=true)

  params = MPIFiles.toDict(sysObj.positions)
  params["currPos"] = sysObj.currPos
  params["stopped"] = sysObj.stopped
  #params["currentSignal"] = sysObj.currentSignal
  params["waitTime"] = sysObj.waitTime
  params["measIsBGPos"] = sysObj.measIsBGPos
  params["posToIdx"] = sysObj.posToIdx
  params["measIsBGFrame"] = sysObj.measIsBGFrame
  params["temperatures"] = vec(sysObj.temperatures)

  open(filename,"w") do f
    TOML.print(f, params)
  end

  Mmap.sync!(sysObj.signals)
  return
end

function restore(sysObj::SystemMatrixRobotMeas)
  filename = "/tmp/sysObj.toml"
  filenameSignals = "/tmp/sysObj.bin"

  if isfile(filename)
    params = TOML.parsefile(filename)
    sysObj.currPos = params["currPos"]
    sysObj.stopped = params["stopped"]
    #sysObj.currentSignal = params["currentSignal"]
    sysObj.waitTime = params["waitTime"]
    sysObj.measIsBGPos = params["measIsBGPos"]
    sysObj.posToIdx = params["posToIdx"]
    sysObj.measIsBGFrame = params["measIsBGFrame"]
    temp = params["temperatures"]
    if !isempty(temp) && (length(sysObj.temperatures) == length(temp))
      sysObj.temperatures[:] .= temp
    end

    sysObj.positions = Positions(params)

    daq = getDAQ(sysObj.scanner)
    numBGPos = sum(sysObj.measIsBGPos)
    numFGPos = length(sysObj.measIsBGPos) - numBGPos

    numTotalFrames = numFGPos + daq.params.acqNumBGFrames*numBGPos

    io = open(filenameSignals, "r+");
    sysObj.signals = Mmap.mmap(io, Array{Float32,4},
          (daq.params.rxNumSamplingPoints,numRxChannels(daq),
           daq.params.acqNumPeriodsPerFrame, numTotalFrames));
  end
  return
end

function init(sysObj::SystemMatrixRobotMeas, positions::GridPositions,
              params::Dict, waitTime = 4.0, safety = getSafety(sysObj.scanner))

  su = getSurveillanceUnit(sysObj.scanner)
  daq = getDAQ(sysObj.scanner)
  robot = getRobot(sysObj.scanner)
  tempSensor = getTemperatureSensor(sysObj.scanner)

  updateParams!(daq, params)

  if !daq.params.controlPhase
    # Repeated in prepareTx in asyncProducer
    tx = daq.params.calibFieldToVolt.*daq.params.dfStrength.*exp.(im*daq.params.dfPhase)
    daq.params.currTx = convert(Matrix{ComplexF64}, diagm(tx))
  end
  #setTxParams(daq, daq.params.currTx*0.0)
  #enableSlowDAC(daq, false)

  rxNumSamplingPoints = daq.params.rxNumSamplingPoints
  numPeriods = daq.params.acqNumPeriodsPerFrame

  measIsBGPos = isa(positions,BreakpointGridPositions) ?
                             MPIFiles.getmask(positions) :
                             zeros(Bool,length(positions))

  numBGPos = sum(measIsBGPos)
  numFGPos = length(measIsBGPos) - numBGPos

  numTotalFrames = numFGPos + daq.params.acqNumBGFrames*numBGPos

  currentSignal = zeros(Float32,rxNumSamplingPoints,numRxChannels(daq),numPeriods,1)
  if tempSensor != nothing
    temperatures = zeros(Float32,numChannels(tempSensor),numTotalFrames)
  else
    temperatures = zeros(0,0)
  end

  filenameSignals = "/tmp/sysObj.bin"
  io = open(filenameSignals, "w+");
  signals = Mmap.mmap(io, Array{Float32,4}, (rxNumSamplingPoints,numRxChannels(daq),numPeriods,numTotalFrames));
  signals[:] .= 0.0

  # The following looks like a secrete line but it makes sense
  posToIdx = cumsum(vcat([false],measIsBGPos)[1:end-1] .* (daq.params.acqNumBGFrames - 1) .+ 1)

  measIsBGFrame = zeros(Bool, numTotalFrames)

  sysObj.params = params
  sysObj.positions = positions
  sysObj.currPos = 1
  sysObj.stopped = false
  sysObj.signals = signals
  sysObj.currentSignal = currentSignal
  sysObj.waitTime = waitTime
  sysObj.measIsBGPos = measIsBGPos
  sysObj.posToIdx = posToIdx
  sysObj.measIsBGFrame = measIsBGFrame
  sysObj.temperatures = temperatures

  return nothing
end

function preMoveAction(measObj::SystemMatrixRobotMeas, pos::Vector{typeof(1.0Unitful.mm)}, index)
  @info "moving to position" pos
end

function postMoveAction(measObj::SystemMatrixRobotMeas,
                        pos::Array{typeof(1.0Unitful.mm),1}, index)
  @info "post action" index length(measObj.positions)

  @info "getThings"
  timeGetThings = @elapsed begin
   safety = getSafety(measObj.scanner)
   su = getSurveillanceUnit(measObj.scanner)
   daq = getDAQ(measObj.scanner)
   robot = getRobot(measObj.scanner)
   tempSensor = getTemperatureSensor(measObj.scanner)
   channel = Channel{channelType(daq)}(32)
  end

  numFrames = measObj.measIsBGPos[index] ? daq.params.acqNumBGFrames : 1

  @info "Starting Measurement"
  timeEnableSlowDAC = @elapsed begin
    actualFrames = daq.params.acqNumFrameAverages*numFrames
    measObj.consumer = @tspawnat 3 asyncConsumer(channel, measObj, index, numFrames)
    asyncProducer(channel, daq, actualFrames, prepTx = false, prepSeq = false, endSeq = false)
  end
  close(channel)

  @show timeEnableSlowDAC
  measObj.producerFinalizer = @async endSequence(daq)
end

function asyncConsumer(channel::Channel, calib::SystemMatrixRobotMeas, index, numFrames)
  @info "readData"
  daq = getDAQ(calib.scanner)
  tempSensor = getTemperatureSensor(calib.scanner)
  asyncBuffer = AsyncBuffer(daq)
  while isopen(channel) || isready(channel)
    while isready(channel)
      chunk = take!(channel)
      updateAsyncBuffer!(asyncBuffer, chunk)
    end
    if !isready(channel)
      sleep(0.001)
    end
  end

  uMeas, uRef = retrieveMeasAndRef!(asyncBuffer, getDAQ(calib.scanner))

  startIdx = calib.posToIdx[index]
  stopIdx = calib.posToIdx[index] + numFrames - 1

  if calib.measIsBGPos[index]
    uMeas_ = reshape(uMeas, size(uMeas,1), size(uMeas,2), size(uMeas,3),
                            daq.params.acqNumFrameAverages, numFrames)
    calib.signals[:,:,:,startIdx:stopIdx] = mean(uMeas_,dims=4)[:,:,:,1,:]
  else
    calib.signals[:,:,:,startIdx] = mean(uMeas,dims=4)
  end

  calib.measIsBGFrame[ startIdx:stopIdx ] .= calib.measIsBGPos[index]

  calib.currentSignal = uMeas[:,:,:,1:1]

  if !isnothing(tempSensor)
    for l in startIdx:stopIdx
      for c = 1:numChannels(tempSensor)
        calib.temperatures[c,l] = getTemperature(tempSensor, c)
      end
    end
  end

  @info "store"
  timeStore = @elapsed store(calib)
  @info "done"

end
###########################

function isStarted(calib::SystemMatrixRobotMeas)
  if calib.task == nothing
    return false
  else
    return !istaskdone(calib.task)
  end
end

function stop(calib::SystemMatrixRobotMeas)
  calib.stopped = true
end

function start(calib::SystemMatrixRobotMeas)
  calib.task = @tspawnat 2 performCalibrationInner(calib)
  return
end

function waitTask(task::Union{Task, Nothing})
  if !isnothing(task)
    wait(task)
  end
end

function performCalibrationInner(calib::SystemMatrixRobotMeas)
  @info  "performCalibrationInner Start"

  try

  su = getSurveillanceUnit(calib.scanner)
  daq = getDAQ(calib.scanner)
  robot = getRobot(calib.scanner)
  safety = getSafety(calib.scanner)
  tempSensor = getTemperatureSensor(calib.scanner)

  positions = calib.positions
  numPos = length(calib.positions)
  calib.stopped = false

  timeRobotMoved = 0.0

  #connect(daq)
  enableACPower(su, calib.scanner)
  stopTx(daq)

  while true
    @info "Curr Pos in performCalibrationInner $(calib.currPos)"
    if calib.stopped
      # Wait for open tasks to finish (consumer will finish either way but this way we are in a known state)
      waitTask(calib.consumer)
      waitTask(calib.producerFinalizer)
      @info  "performCalibrationInner stopped"
      break
    end

    if calib.currPos <= numPos
        pos = uconvert.(Unitful.mm, positions[calib.currPos])

        timeMove = 0
        timePrepDAQ = 0
        timeFinalizer = 0
        timeFrameChange = 0
        timeSeq = 0
        timeTx = 0
        timeConsumer = 0
        timePreparing = @elapsed @sync begin

          # Prepare Robot/Sample
          @async timeMove = @elapsed moveAbsUnsafe(robot, pos)

          # Prepare DAQ
          @async timePrepDAQ = @elapsed begin
            allowControlLoop = mod1(calib.currPos, 11) == 1
            timeFinalizer = @elapsed waitTask(calib.producerFinalizer)
            timeFrameChange = @elapsed begin
              if (calib.currPos == 1) ||
               (calib.measIsBGPos[calib.currPos] != calib.measIsBGPos[calib.currPos-1])
              daq.params.acqNumFrames = calib.measIsBGPos[calib.currPos] ? daq.params.acqNumBGFrames : 1
              setSequenceParams(daq)
              end end
            timeSeq = @elapsed prepareSequence(daq)
            timeTx = @elapsed prepareTx(daq, allowControlLoop = allowControlLoop)
          end

          @async timeConsumer = @elapsed waitTask(calib.consumer)
        end

        diffTime = calib.waitTime - timePreparing
        if diffTime > 0.0
          sleep(diffTime)
        end

        setEnabled(robot, false)

        timePostMove = @elapsed postMoveAction(calib, pos, calib.currPos)

        @info "############### Preparing Time: $(timePreparing), Prep DAQ time $timePrepDAQ, Finalizer $timeFinalizer, Frame $timeFrameChange, Seq $timeSeq, Tx $timeTx, Consumer $timeConsumer"
        @info "Move robot $timeMove, meas time: $(timePostMove)"

        setEnabled(robot, true)
        calib.currPos +=1

    end

    if calib.currPos > numPos
        waitTask(calib.consumer)
        waitTask(calib.producerFinalizer)
        @info "Store SF"
        stopTx(daq)
        disableACPower(su, calib.scanner)
        MPIMeasurements.disconnect(daq)

        movePark(robot)
        #calib.currPos = 0

        #if !calib.stopped
        saveasMDF(calib)
        #end
        break
    end

  end
  catch ex
    @warn "Exception" ex stacktrace(catch_backtrace())
  end
end

#############  Storage ##########################

function MPIFiles.saveasMDF(calibObj::SystemMatrixRobotMeas)
  store = calibObj.store
  params = calibObj.params

  if params["storeAsSystemMatrix"]
    calibNum = getNewCalibNum(store)
    params["experimentNumber"] = calibNum
    filenameA = joinpath(calibdir(store),string(calibNum)*".mdf") # just for debugging
    #filenameA = "/tmp/tmp.mdf"
    filenameB = joinpath(calibdir(store),string(calibNum+1)*".mdf")
    saveasMDF(filenameA, calibObj, params)
    #saveasMDF(filenameB, MPIFile(filenameA), applyCalibPostprocessing=true)
  else

    name = params["studyName"]
    date = params["studyDate"]
    subject = ""

    newStudy = Study(store, name; subject=subject, date=date)
    expNum = getNewExperimentNum(newStudy)

    params["experimentNumber"] = expNum

    filename = joinpath(path(newStudy), string(expNum)*".mdf")

    saveasMDF(filename, calibObj, params)
  end
end

function MPIFiles.saveasMDF(filename::String, measObj::SystemMatrixRobotMeas, params_::Dict)
  params = copy(params_)

  daq = getDAQ(measObj.scanner)
  positions = measObj.positions

  # acquisition parameters
  params["acqStartTime"] = Dates.unix2datetime(time())
  params["acqGradient"] = addTrailingSingleton([0.0;0.0;0.0],2) #FIXME
  params["acqOffsetField"] = addTrailingSingleton([0.0;0.0;0.0],2) #FIXME

  # drivefield parameters
  params["dfStrength"] = reshape(daq.params.dfStrength,1,length(daq.params.dfStrength),1)
  params["dfPhase"] = reshape(daq.params.dfPhase,1,length(daq.params.dfPhase),1)
  params["dfDivider"] = reshape(daq.params.dfDivider,1,length(daq.params.dfDivider))

  # receiver parameters
  params["rxNumSamplingPoints"] = daq.params.rxNumSamplingPoints
  params["rxNumChannels"] = numRxChannels(daq)

  # calibration params  (needs to be called after calibration params!)
  #params["rxDataConversionFactor"] = calibIntToVoltRx(daq)
  calib = zeros(2,numRxChannels(daq))
  calib[1,:] .= 1.0
  params["rxDataConversionFactor"] = calib

  params["measIsFourierTransformed"] = false
  params["measIsSpectralLeakageCorrected"] = false
  params["measIsTFCorrected"] = false
  params["measIsFrequencySelection"] = false
  params["measIsBGCorrected"] = false
  params["measIsFastFrameAxis"] = false
  params["measIsFramePermutation"] = false

  params["acqNumFrames"] = length(measObj.measIsBGFrame)

  params["measIsBGFrame"] = measObj.measIsBGFrame

  params["measData"] = measObj.signals

  subgrid = isa(positions,BreakpointGridPositions) ? positions.grid : positions

  params["calibIsMeanderingGrid"] = isa(subgrid,MeanderingGridPositions)

  #params["calibSNR"] TODO during conversion
  if !(typeof(subgrid) <: Union{SphericalTDesign,ArbitraryPositions})
    # parameter for SphericalTDesign and ArbitraryPositions not defined
    params["calibFov"] = Float64.(ustrip.(uconvert.(Unitful.m, fieldOfView(subgrid))))
    params["calibFovCenter"] = Float64.(ustrip.(uconvert.(Unitful.m, fieldOfViewCenter(subgrid))))
    params["calibSize"] = shape(subgrid)
  end

  params["calibOrder"] = "xyz"
  if haskey(params, "calibDeltaSampleSize")
    params["calibDeltaSampleSize"] =
       Float64.(ustrip.(uconvert.(Unitful.m, params["calibDeltaSampleSize"])))
  end
  params["calibMethod"] = "robot"

  if !isempty(measObj.temperatures)
    params["calibTemperatures"] = measObj.temperatures
  end

  @info "save as MDF"
  MPIFiles.saveasMDF( filename, params )
  @info "cleanup"
  cleanup(measObj)
  @info "done"
  return filename
end

ext: jl


text:
#normalization factor
function ylmKCoefficient(l::Integer, m::Integer)
    k = BigInt(4)

    for i in l-m+1:l+m
        # use BigInt to avoid overflow
        k *= i
    end

    return sqrt(Float64((2*l+1) / k)/pi)
end

function ylmCosSinPolynomial(m::Int64, x::Variable, y::Variable)

  sum = 0*(x+y)
  for j::Int64 in 0:floor(m/2)
    sum += ((-1)^j)*binomial(m, 2*j)*(y^(2*j))*(x^(m-2*j))
  end
  return sum
end

function ylmSinSinPolynomial(m::Int64, x::Variable, y::Variable)

  sum = 0*(x+y)
  for j::Int64 in 0:floor((m-1)/2)
    sum += ((-1)^j)*binomial(m, 2*j + 1)*(y^(2*j + 1))*(x^(m-2*j-1))
  end
  return sum
end

"""
    ylm(l::Int64, m::Int64, x::Variable, y::Variable, z::Variable)
*Description:*  Calculation of the spherical harmonic for a given order (l,m) in Cartesian coordinates\\

*Input:*  `l`       - Degree of the spherical harmonic\\
          `m`       - Order of the spherical harmonic\\
          `x, y, z` - Cartesian coordinates\\

*Output:*  Spherical harmonic polynomial
"""
function ylm(l::Int64, m::Int64, x::Variable, y::Variable, z::Variable)

  if abs(m) > l
    throw(DomainError(m,"-l <= m <= l expected, but m = $m and l = $l."))
  end

  p = 1.0*(z^2 - 1)^l + 0.0*(x+y)

  for i = 1:l+abs(m)
    c = i <= l ? 1/(2*i) : 1.0
    p = c*differentiate(p, z)
  end

  if m > 0
    return sqrt(2)*ylmKCoefficient(l, m)*ylmCosSinPolynomial(m,x,y)*p
  elseif m < 0
    return sqrt(2)*ylmKCoefficient(l, abs(m))*ylmSinSinPolynomial(abs(m),x,y)*p
  else
    return ylmKCoefficient(l, 0)*p
  end
end

# multiplying r^l*ylm(x,y,z)
function rlylm(l::Int64, m::Int, x::Variable, y::Variable, z::Variable)
	p = ylm(l,m,x,y,z)
	tout = []
	# Zerlegung des Polynoms in Terme:
	for t in terms(p)
		deg = degree(monomial(t)) # Gibt den gesamten Grad des Monoms an
		degR = l-deg # durch das Kürzen ergibt sich ein Grad von l-deg fuer r
		push!(tout,(x^2+y^2+z^2)^Int(degR/2)*t) # r² wird durch x²+y²+z² ersetzt
	end

	return polynomial(tout)
end

# solid harmonics
function zlm(l::Int64, m::Int64, x::Variable, y::Variable, z::Variable)
	rlm = rlylm(l,m,x,y,z)
	rlm = sqrt(4*pi/(2*l+1))*rlm
	return rlm
end


text:
function del_blocked(S, b, lb, ub, rxns; eps = 0.0, protect = [])

    lb, ub = (lb, ub) .|> copy # local copy
    m, n = size(S)

    _bidx = trues(n)
    _bidx[protect] .= false
    non_protected = findall(_bidx)

    blocked = falses(n)
    for i in non_protected
        if lb[i] == ub[i] # blocked
            blocked[i] = eps == 0.0
            lb[i], ub[i] = lb[i] - eps, ub[i] + eps
        end
    end

    unblocked = findall(.!blocked)

    return S[:,unblocked], b - S*(blocked .* lb),
        lb[unblocked], ub[unblocked], rxns[unblocked], findall(blocked)
end

function del_blocked(model::MetNet; eps = 0.0, protect = [])
    protect = map((r) -> rxnindex(model, r), protect)
    S, b, lb, ub, rxns, blocked =
        del_blocked(model.S, model.b, model.lb, model.ub, model.rxns; eps, protect)
    return MetNet(model; reshape = true, S, b, lb, ub, rxns)
end


text:
using InteractiveUtils

const Emptykwargs = Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}

macro warnpcfail(ex::Expr)
    modl = __module__
    file = __source__.file === nothing ? "?" : String(__source__.file)
    line = __source__.line
    quote
        $(esc(ex)) || @warn """precompile directive
     $($(Expr(:quote, ex)))
 failed. Please report an issue in $($modl) (after checking for duplicates) or remove this directive.""" _file=$file _line=$line
    end
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @warnpcfail precompile(BBox, (Float32, Float32, Float64, Float64))

    @warnpcfail precompile(GridLayout, (Int, Int))
    @warnpcfail precompile(GridLayout, (GridLayoutSpec,))
    @warnpcfail precompile(Core.kwfunc(GridLayout), (NamedTuple{(:bbox,), Tuple{Observable{IRect2D}}}, Type{GridLayout}))
    @warnpcfail precompile(Core.kwfunc(GridLayout), (NamedTuple{(:bbox,), Tuple{FRect2D}}, Type{GridLayout}))
    # Also precompile `GridLayout(bbox = bbox, alignmode = al)` and the keyword body method
    fbody = isdefined(Base, :bodyfunction) ? Base.bodyfunction(which(GridLayout, (Int, Int))) : nothing
    for Al in (Inside, Outside, Mixed)
        @warnpcfail precompile(Core.kwfunc(GridLayout), (NamedTuple{(:bbox, :alignmode), Tuple{FRect2D, Al}}, Type{GridLayout}))
        @warnpcfail precompile(Core.kwfunc(GridLayout), (NamedTuple{(:bbox, :alignmode), Tuple{Observable{FRect2D}, Al}}, Type{GridLayout}))
        @warnpcfail precompile(Core.kwfunc(GridLayout), (NamedTuple{(:bbox, :alignmode), Tuple{NTuple{4,Int}, Al}}, Type{GridLayout}))
        if fbody !== nothing
            @warnpcfail precompile(fbody, (Nothing, Nothing, Nothing, Nothing, Al, Tuple{Bool, Bool}, FRect2D, Auto, Auto, Bool, Bool, Symbol, Symbol, Float64, Float64, Emptykwargs, Type{GridLayout}, Int, Int))
        end
    end
    @warnpcfail precompile(Core.kwfunc(GridLayout), (NamedTuple{(:colsizes, :rowsizes), Tuple{Fixed, Relative}}, Type{GridLayout}, Int, Int))
    @warnpcfail precompile(Core.kwfunc(GridLayout), (NamedTuple{(:addedcolgaps,), Tuple{Vector{Fixed}}}, Type{GridLayout}, Int, Int))

    @warnpcfail precompile(Core.kwfunc(Mixed), (NamedTuple{(:left, :top), Tuple{Int, Int}}, Type{Mixed}))
    for T in (Int, Float32)
        @warnpcfail precompile(Outside, (T,))
        @warnpcfail precompile(Outside, (T,T,T,T))
    end
    @warnpcfail precompile(align_to_bbox!, (GridLayout, FRect2D))
    @warnpcfail precompile(update_gl!, (GridLayout,))
    for I in subtypes(Indexables), J in subtypes(Indexables)
        @warnpcfail precompile(setindex!, (GridLayout, UnitRange{Int}, I, J))
        @warnpcfail precompile(setindex!, (GridLayout, Any, I, J))
    end
    @warnpcfail precompile(insertrows!, (GridLayout, Int, Int))
    @warnpcfail precompile(insertcols!, (GridLayout, Int, Int))
    @warnpcfail precompile(gridnest!, (GridLayout, UnitRange{Int}, UnitRange{Int}))
    @warnpcfail precompile(add_to_gridlayout!, (GridLayout, GridContent{GridLayout, GridLayout}))
    @warnpcfail precompile(connect_layoutobservables!, (GridContent{GridLayout, GridLayout},))
    @warnpcfail precompile(trim!, (GridLayout,))
    @warnpcfail precompile(sizeobservable!, (Observable{Any}, Observable{Any}))
    @warnpcfail precompile(_reportedsizeobservable, (Tuple{SizeAttribute,SizeAttribute}, Tuple{AutoSize,AutoSize}, AlignMode, RectSides{Float32}, Tuple{Bool,Bool}))
    for T in subtypes(Side)
        @warnpcfail precompile(bbox_for_solving_from_side, (RowCols{Vector{Float64}}, FRect2D, RowCols{Int}, T))
    end
    for S in (Left, Right, Top, Bottom)
        @warnpcfail precompile(protrusion, (GridContent{GridLayout,GridLayout}, S))
    end
    @warnpcfail precompile(suggestedbboxobservable, (GridLayout,))
    for VC in (Vector{Auto}, Vector{GapSize}, Vector{Fixed}, Vector{Relative}, Vector{ContentSize})
        @warnpcfail precompile(convert_contentsizes, (Int, VC))
        @warnpcfail precompile(==, (VC, VC))
    end
    @warnpcfail precompile(contents, (GridLayout,))
    @warnpcfail precompile(filterenum, (Function, Type, Vector{ContentSize}))
    @warnpcfail precompile(zcumsum, (Vector{Float64},))

    # These don't work completely but they have partial success
    @warnpcfail precompile(repr, (MIME{Symbol("text/plain")}, GridLayout))
end


text:
module NonOverlapping

include("./constants.jl")
using JuMP
using .Constants

export cons_non_overlapping

#Constraints ensuring that the only allowed rotations are those of Currence Robotics
function cons_non_overlapping(model)

    # ====== RELATIVE POSISJONER START =======
    #For each pair of boxes (i, j), check j's relative position to i (from i's perspective)
    @variable(model, aij_j_left_of_i[i=1:length(boxes), j=1:length(boxes); i!=j], Bin) #left
    @variable(model, bij_j_right_of_i[i=1:length(boxes), j=1:length(boxes); i!=j], Bin) #right
    @variable(model, cij_j_behind_i[i=1:length(boxes), j=1:length(boxes); i!=j], Bin) #behind
    @variable(model, dij_j_in_front_i[i=1:length(boxes), j=1:length(boxes); i!=j], Bin) #in front
    @variable(model, eij_j_below_i[i=1:length(boxes), j=1:length(boxes); i!=j], Bin) #below
    @variable(model, fij_j_above_i[i=1:length(boxes), j=1:length(boxes); i!=j], Bin) #above

    #One box can only be at relative to another box with 1 position.
    @constraint(model,
        rel_pos[i=1:length(boxes), j=1:length(boxes); i!=j],
        (model[:aij_j_left_of_i][i, j] +
        model[:bij_j_right_of_i][i, j] +
        model[:cij_j_behind_i][i, j] +
        model[:dij_j_in_front_i][i, j] +
        model[:eij_j_below_i][i, j] +
        model[:fij_j_above_i][i, j] >= 1)
    )
    # ====== RELATIVE POSISJONER SLUTT =======

    # ====== IKKE OVERLAPPENDE START =======
    #left
    @constraint(model,
        left[i=1:length(boxes), j=1:length(boxes); i!=j],
        (model[:flbx_front_left_bot_x][i] +
        fix_value(model[:lob_length_of_box][i])*model[:lpx_length_of_box_parallel_to_x][i] +
        fix_value(model[:wob_width_of_box][i])*model[:wpx_width_of_box_parallel_to_x][i] +
        fix_value(model[:hob_height_of_box][i])*model[:hpx_height_of_box_parallel_to_x][i]) <=
        (model[:flbx_front_left_bot_x][j] + (1-model[:aij_j_left_of_i][i, j])*L)
    )

    #right
    @constraint(model,
        right[i=1:length(boxes), j=1:length(boxes); i!=j],
        (model[:flbx_front_left_bot_x][j] +
        fix_value(model[:lob_length_of_box][j])*model[:lpx_length_of_box_parallel_to_x][j] +
        fix_value(model[:wob_width_of_box][j])*model[:wpx_width_of_box_parallel_to_x][j] +
        fix_value(model[:hob_height_of_box][j])*model[:hpx_height_of_box_parallel_to_x][j]) <=
        (model[:flbx_front_left_bot_x][i] + (1-model[:bij_j_right_of_i][i, j])*L)
    )

    #behind
    @constraint(model,
        behind[i=1:length(boxes), j=1:length(boxes); i!=j], (model[:flby_front_left_bot_y][i] +
        fix_value(model[:lob_length_of_box][i])*model[:lpy_length_of_box_parallel_to_y][i] +
        fix_value(model[:wob_width_of_box][i])*model[:wpy_width_of_box_parallel_to_y][i] +
        fix_value(model[:hob_height_of_box][i])*model[:hpy_height_of_box_parallel_to_y][i]) <=
        (model[:flby_front_left_bot_y][j] + (1-model[:cij_j_behind_i][i, j])*W)
    )

    #in_front
    @constraint(model,
        in_front[i=1:length(boxes), j=1:length(boxes); i!=j],
        (model[:flby_front_left_bot_y][j] +
        fix_value(model[:lob_length_of_box][j])*model[:lpy_length_of_box_parallel_to_y][j] +
        fix_value(model[:wob_width_of_box][j])*model[:wpy_width_of_box_parallel_to_y][j] +
        fix_value(model[:hob_height_of_box][j])*model[:hpy_height_of_box_parallel_to_y][j]) <=
        (model[:flby_front_left_bot_y][i] + (1-model[:dij_j_in_front_i][i, j])*W)
    )

    #below
    @constraint(model,
        below[i=1:length(boxes), j=1:length(boxes); i!=j],
        (model[:flbz_front_left_bot_z][i] +
        fix_value(model[:lob_length_of_box][i])*model[:lpz_length_of_box_parallel_to_z][i] +
        fix_value(model[:wob_width_of_box][i])*model[:wpz_width_of_box_parallel_to_z][i] +
        fix_value(model[:hob_height_of_box][i])*model[:hpz_height_of_box_parallel_to_z][i]) <=
        (model[:flbz_front_left_bot_z][j] + (1-model[:eij_j_below_i][i, j])*H)
    )

    #above
    @constraint(model,
        above[i=1:length(boxes), j=1:length(boxes); i!=j],
        (model[:flbz_front_left_bot_z][j] +
        fix_value(model[:lob_length_of_box][j])*model[:lpz_length_of_box_parallel_to_z][j] +
        fix_value(model[:wob_width_of_box][j])*model[:wpz_width_of_box_parallel_to_z][j] +
        fix_value(model[:hob_height_of_box][j])*model[:hpz_height_of_box_parallel_to_z][j]) <=
        (model[:flbz_front_left_bot_z][i] + (1-model[:fij_j_above_i][i, j])*H)
    )
    # ====== IKKE OVERLAPPENDE SLUTT =======

    return model
end

end


text:
push!(LOAD_PATH, "/home/ian/Grad/Research/Disks/DiskJockey/")

using LittleMC

using Distributions
using PDMats

const mat = PDiagMat([2.0^2, 3.0^2])
const dist = DiagNormal([0., 0.], mat)

function Gauss(x::Vector{Float64})
    return logpdf(dist, x)
end

using PDMats

jump = full(PDiagMat([1.5^2, 1.5^2]))

println(typeof(jump))

mc = MC(Gauss, 100000, [1.0, 1.0], jump)

start(mc)

println(mean(mc.samples, 2))
println(std(mc.samples, 2))

runstats(mc)

write(mc, "mc.hdf5")

# # Plot the samples
# import PyPlot.plt
# fig, ax = plt.subplots(nrows=2, sharex=true)
# ax[1][:plot](mc.samples'[:, 1])
# ax[2][:plot](mc.samples'[:, 2])
# plt.savefig("plots/littlemc.png")

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
using Sundials, Test

h1 = Sundials.Handle(Sundials.CVodeCreate(Sundials.CV_BDF))
h2 = h1

@test !isempty(h1)
@test !isempty(h2)

empty!(h1)

@test isempty(h1)
@test isempty(h2) # Check that the pointers are correct, and h1 empties h2

empty!(h1) # Make sure this does not throw

h = Sundials.Handle(h1.ptr_ref[]) # Check construction with null pointers
@test isempty(h)

neq = 3
h3 = Sundials.MatrixHandle(Sundials.SUNDenseMatrix(neq, neq), Sundials.DenseMatrix())
h3 = Sundials.MatrixHandle(Sundials.SUNDenseMatrix(neq, neq), Sundials.DenseMatrix())
empty!(h3)
@test isempty(h3)
empty!(h3)
@test isempty(h3)

h3 = Sundials.MatrixHandle(Sundials.SUNBandMatrix(100, 3, 3), Sundials.BandMatrix())
h3 = Sundials.MatrixHandle(Sundials.SUNBandMatrix(100, 3, 3), Sundials.BandMatrix())
empty!(h3)
@test isempty(h3)
empty!(h3)
@test isempty(h3)

h3 = Sundials.MatrixHandle(
    Sundials.SUNSparseMatrix(neq, neq, neq, Sundials.CSC_MAT),
    Sundials.SparseMatrix(),
)
h3 = Sundials.MatrixHandle(
    Sundials.SUNSparseMatrix(neq, neq, neq, Sundials.CSC_MAT),
    Sundials.SparseMatrix(),
)
empty!(h3)
@test isempty(h3)
empty!(h3)
@test isempty(h3)

A = Sundials.SUNDenseMatrix(neq, neq)
u0 = rand(neq)
Sundials.SUNLinSol_Dense(u0, A)
h3 = Sundials.LinSolHandle(Sundials.SUNLinSol_Dense(u0, A), Sundials.Dense())
h3 = Sundials.LinSolHandle(Sundials.SUNLinSol_Dense(u0, A), Sundials.Dense())
empty!(h3)
@test isempty(h3)
empty!(h3)
@test isempty(h3)


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

getfield(Core, :eval)(Core, :(baremodule Compiler

using Core.Intrinsics, Core.IR

import Core: print, println, show, write, unsafe_write, stdout, stderr,
             _apply_iterate, svec, apply_type, Builtin, IntrinsicFunction,
             MethodInstance, CodeInstance, MethodMatch, PartialOpaque,
             TypeofVararg

const getproperty = Core.getfield
const setproperty! = Core.setfield!
const swapproperty! = Core.swapfield!
const modifyproperty! = Core.modifyfield!
const replaceproperty! = Core.replacefield!

ccall(:jl_set_istopmod, Cvoid, (Any, Bool), Compiler, false)

eval(x) = Core.eval(Compiler, x)
eval(m, x) = Core.eval(m, x)

include(x) = Core.include(Compiler, x)
include(mod, x) = Core.include(mod, x)

# The @inline/@noinline macros that can be applied to a function declaration are not available
# until after array.jl, and so we will mark them within a function body instead.
macro inline()   Expr(:meta, :inline)   end
macro noinline() Expr(:meta, :noinline) end

convert(::Type{Any}, Core.@nospecialize x) = x
convert(::Type{T}, x::T) where {T} = x

# essential files and libraries
include("essentials.jl")
include("ctypes.jl")
include("generator.jl")
include("reflection.jl")
include("options.jl")

# core operations & types
function return_type end # promotion.jl expects this to exist
is_return_type(@Core.nospecialize(f)) = f === return_type
include("promotion.jl")
include("tuple.jl")
include("pair.jl")
include("traits.jl")
include("range.jl")
include("expr.jl")
include("error.jl")

# core numeric operations & types
==(x::T, y::T) where {T} = x === y
include("bool.jl")
include("number.jl")
include("int.jl")
include("operators.jl")
include("pointer.jl")
include("refvalue.jl")

# the same constructor as defined in float.jl, but with a different name to avoid redefinition
_Bool(x::Real) = x==0 ? false : x==1 ? true : throw(InexactError(:Bool, Bool, x))

# checked arithmetic
const checked_add = +
const checked_sub = -
const SignedInt = Union{Int8,Int16,Int32,Int64,Int128}
const UnsignedInt = Union{UInt8,UInt16,UInt32,UInt64,UInt128}
sub_with_overflow(x::T, y::T) where {T<:SignedInt}   = checked_ssub_int(x, y)
sub_with_overflow(x::T, y::T) where {T<:UnsignedInt} = checked_usub_int(x, y)
sub_with_overflow(x::Bool, y::Bool) = (x-y, false)
add_with_overflow(x::T, y::T) where {T<:SignedInt}   = checked_sadd_int(x, y)
add_with_overflow(x::T, y::T) where {T<:UnsignedInt} = checked_uadd_int(x, y)
add_with_overflow(x::Bool, y::Bool) = (x+y, false)

include("strings/lazy.jl")

# core array operations
include("indices.jl")
include("array.jl")
include("abstractarray.jl")

# core structures
include("bitarray.jl")
include("bitset.jl")
include("abstractdict.jl")
include("iddict.jl")
include("idset.jl")
include("abstractset.jl")
include("iterators.jl")
using .Iterators: zip, enumerate
using .Iterators: Flatten, Filter, product  # for generators
include("namedtuple.jl")

ntuple(f, ::Val{0}) = ()
ntuple(f, ::Val{1}) = (@inline; (f(1),))
ntuple(f, ::Val{2}) = (@inline; (f(1), f(2)))
ntuple(f, ::Val{3}) = (@inline; (f(1), f(2), f(3)))
ntuple(f, ::Val{n}) where {n} = ntuple(f, n::Int)
ntuple(f, n) = (Any[f(i) for i = 1:n]...,)

# core docsystem
include("docs/core.jl")
import Core.Compiler.CoreDocs
Core.atdoc!(CoreDocs.docm)

# sorting
function sort end
function sort! end
function issorted end
function sortperm end
include("ordering.jl")
using .Order
include("sort.jl")
using .Sort

# We don't include some.jl, but this definition is still useful.
something(x::Nothing, y...) = something(y...)
something(x::Any, y...) = x

############
# compiler #
############

include("compiler/cicache.jl")
include("compiler/effects.jl")
include("compiler/types.jl")
include("compiler/utilities.jl")
include("compiler/validation.jl")
include("compiler/methodtable.jl")

function argextype end # imported by EscapeAnalysis
function stmt_effect_free end # imported by EscapeAnalysis
function alloc_array_ndims end # imported by EscapeAnalysis
function try_compute_field end # imported by EscapeAnalysis
include("compiler/ssair/basicblock.jl")
include("compiler/ssair/domtree.jl")
include("compiler/ssair/ir.jl")

include("compiler/inferenceresult.jl")
include("compiler/inferencestate.jl")

include("compiler/typeutils.jl")
include("compiler/typelimits.jl")
include("compiler/typelattice.jl")
include("compiler/tfuncs.jl")
include("compiler/stmtinfo.jl")

include("compiler/abstractinterpretation.jl")
include("compiler/typeinfer.jl")
include("compiler/optimize.jl") # TODO: break this up further + extract utilities

# required for bootstrap because sort.jl uses extrema
# to decide whether to dispatch to counting sort.
#
# TODO: remove it.
function extrema(x::Array)
    isempty(x) && throw(ArgumentError("collection must be non-empty"))
    vmin = vmax = x[1]
    for i in 2:length(x)
        xi = x[i]
        vmax = max(vmax, xi)
        vmin = min(vmin, xi)
    end
    return vmin, vmax
end

include("compiler/bootstrap.jl")
ccall(:jl_set_typeinf_func, Cvoid, (Any,), typeinf_ext_toplevel)

include("compiler/parsing.jl")
Core.eval(Core, :(_parse = Compiler.fl_parse))

end # baremodule Compiler
))

max_stars_repo_name: TechPenguineer/julia

max_issues_repo_name: TechPenguineer/julia

max_forks_repo_name: TechPenguineer/julia


text:
"""
    getWebLinks(url::AbstractString, fileTypes::AbstractVector, downloadPath:: AbstractString)
Call this function to crawl through a web page and download all file links

first parameter is webpage url, next parameter is list of all file extensions user wish to download followed by output directory path

# Example
```julia-repl
julia> fl = getWebLinks("https://investor.apple.com/investor-relations/default.aspx#tabs_content--2021", ["pdf","csv","xlsx","xls"], "downloads/web")
```
"""
function getWebLinks(url::AbstractString, fileTypes::AbstractVector, downloadPath::AbstractString)
	return uppercase(reduce(replace, ["-" => "_"," " => "_", "\$" => "_","#" => "_","%" => "_","(" => "_",")" => "_"], init=url))
end

"""
    getFile(url::AbstractString, downloadPath:: AbstractString)
Call this function to download a file

first parameter is webpage url, next parameter is output directory path including file name

# Example
```julia-repl
julia> fl = getFile("https://s2.q4cdn.com/470004039/files/doc_financials/2021/q2/FY21-Q2-Consolidated-Financial-Statements.pdf", "<folder_name>/test.csv")
```
"""
function getFile(url::AbstractString, downloadPathWithFileName::AbstractString)
	download(url, downloadPathWithFileName)
end

"""
    getPullFiles(path::AbstractString)
Call this function to read urls in local txt files line by line, and download each file

# Example
```julia-repl
julia> fl = getPullFiles("c:\amit.la\filename.txt")
```
"""
function getPullFiles(path::AbstractString)
    # source code release date Aug 20, 2021
	return "Download complete"
end

"""
    getJSONintoDataFrame(path::AbstractString)
Call this function to read json file from url and retrieve results into DataFrame

first parameter is webpage url, next parameter is output directory path

# Example
```julia-repl
julia> df = getJSONintoDataFrame("https://api.coindesk.com/v1/bpi/currentprice.json", "downloads/web")
```
"""
function getJSONintoDataFrame(url::AbstractString, downloadPath::AbstractString)
    # source code release date Aug 20, 2021
	return "Download complete"
end

"""
    getXMLintoDataFrame(path::AbstractString)

Call this function to read XML file from url and retrieve results into DataFrame

first parameter is webpage url, next parameter is output directory path
# Example
```julia-repl
julia> df = getXMLintoDataFrame("https://api.coindesk.com/v1/bpi/currentprice.json", "downloads/web")
```
"""
function getXMLintoDataFrame(url::AbstractString, downloadPath::AbstractString)
    # source code release date Aug 20, 2021
	return "Download complete"
end


text:
using LinearAlgebra
using FFTW, Images, ImageIO, TiffImages, ImageFiltering, Interpolations, NLopt

struct Range2D{T}
    x::UnitRange{T}
    y::UnitRange{T}
end

"""
Given an accelerating voltage in kV, return the electron wavelength in picometers.
Taken from ctffind4 code.
"""
function wavelength_from_voltage(voltage)
    1226.39 / sqrt(1000.0 * voltage + 0.97845e-6 * (1000.0 * voltage)^2)
end

""" Return a log averaged-spectrum using Bartlett's method. """
function get_spectrum(img, N::Integer = 1024, stride = 1024)
    m, n = size(img)
    P = plan_rfft(@view img[1:N, 1:N])
    A = abs.(P * (@view img[1:N, 1:N]))
    for j in stride+1:stride:n-N, i in stride+1:stride:m-N
        A .+= abs.(P * (@view img[i:i+1023, j:j+1023]))
    end
    return log.(A)
end

"""
Given the power spectrum of an image (rfft without fftshift),
return a 1-D profile in the horizontal direction.

width: number of rows to average
σ: std. deviation of gaussian smoothing
"""
function profile1D(spectrum::AbstractArray{<:Real}; σ = 2, width::Integer = 15)
    j_max = div(size(spectrum)[2], 2)
    horizontal_strip = @view spectrum[1:width, 1:j_max]
    profile = vec(sum(horizontal_strip, dims = 1)) / width
    if σ > 0
        imfilter!(profile, profile, KernelFactors.IIRGaussian((σ,)))
    end
    return profile
end

""" Linearly interpolates the min peaks of a 1-D profile.

Returns the interpolation as a function of the array indices.
"""
function get_background(search_profile::AbstractArray)
    # 2 iterations to get the min peaks to zero
    n_iters = 2
    itp = Vector{Any}(undef, n_iters)
    for i = 1:n_iters
        local_mins = [x[1] for x in findlocalminima(search_profile)]
        itp[i] = LinearInterpolation(
            local_mins,
            search_profile[local_mins],
            extrapolation_bc = Line(),
        )
        lower_envelope = itp[i].(1:length(search_profile))
        search_profile -= lower_envelope
    end
    return x -> sum(itp[i](x) for i = 1:n_iters)
end

""" Linearly interpolates the thon ring peaks of a background-noise subtracted 1-D profile.

Returns the interpolation as a function of the array indices.
"""
function get_envelope_function(subtracted_profile::AbstractArray)
    local_maxes = [x[1] for x in findlocalmaxima(subtracted_profile)]
    return LinearInterpolation(
        local_maxes,
        subtracted_profile[local_maxes],
        extrapolation_bc = Line(),
    )
end

""" Given x and y (2D meshgrid arrays), return a bandpass binary mask. """
function create_binary_mask(x::AbstractArray{Int,2}, y::AbstractArray{Int,2}, low, high)
    r = @. sqrt(x^2 + y^2)
    mask = low .< r .< high
    return mask
end

""" Given x and y (2D meshgrid arrays), return a bandpass mask that is damped by an envelope function. """
function create_damped_mask(x::AbstractArray, y::AbstractArray, low, high, envelope)
    r = @. sqrt(x^2 + y^2)
    mask = Float64.(low .< r .< high)
    @. mask *= envelope(r - low)
    return mask
end

""" Subtract noise background and calculate approximate defocus for optimization.

Returns a reference spectrum (fftshifted), rough defocus estimate, a damped mask for simulation,
and meshgrid arrays kx & ky of spatial frequencies.

`Cs`: spherical abberation in millimeters.
`λ`: electron wavelength in picometers.
`pixelsize`: real space pixel size in Angstroms.
`min_resolution`: in Angstroms.
`max_resolution`: in Angstroms.
"""
function preprocess_spectrum(
    spectrum::AbstractArray,
    Cs,
    λ,
    pixelsize,
    min_resolution,
    max_resolution,
)
    N = size(spectrum)[2]
    low_res_idx = round(Int, N * pixelsize / min_resolution)
    high_res_idx = round(Int, N * pixelsize / max_resolution)

    mid = ceil(Int, (N + 1) / 2)
    search_range2D = Range2D(1:high_res_idx, mid-(high_res_idx-1):mid+(high_res_idx-1))
    shifted_spectrum = fftshift(spectrum, 2)
    search_spectrum = shifted_spectrum[search_range2D.x, search_range2D.y]

    # calculate background and envelope functions from 1D profile
    ctfProfile = profile1D(real(spectrum))
    search_profile = ctfProfile[low_res_idx:high_res_idx]
    background_func = get_background(search_profile)
    subtracted_profile = search_profile - background_func.(1:length(search_profile))
    envelope = get_envelope_function(subtracted_profile)

    # estimate approximate defocus
    second_peak_idx = low_res_idx + envelope.itp.knots[1][2] - 1
    k_peak = second_peak_idx / (N * pixelsize)
    z_estimate = (1.5 + 0.5 * Cs * λ^3 * k_peak^4 * 10) / (λ * k_peak^2 * 100)

    kx = (
        fftshift(fftfreq(N, 1 / pixelsize))[search_range2D.y]' .*
        ones(length(search_range2D.x))
    )
    ky = ones(length(search_range2D.y))' .* rfftfreq(N, 1 / pixelsize)[search_range2D.x]
    x = round.(Int, kx * N * pixelsize)
    y = round.(Int, ky * N * pixelsize)
    r = @. sqrt(x^2 + y^2)

    background_spectrum = background_func.(r .- low_res_idx)
    subtracted_spectrum = search_spectrum - background_spectrum
    binary_mask = create_binary_mask(x, y, low_res_idx, high_res_idx)
    reference_spectrum = subtracted_spectrum .* binary_mask
    damped_mask = create_damped_mask(x, y, low_res_idx, high_res_idx, envelope)

    return (reference_spectrum, z_estimate, damped_mask, kx, ky, envelope, low_res_idx)
end

""" Calculate ctf in place

`A`: array to store result.
`z1`, `z2`: defocus values in microns.
`θ`: astigmatism angle in degrees.
`ϕ`: phase shift in degrees.
`kx`, `ky`: meshgrid spatial frequencies.
`λ`: electron wavelength in picometers.
`Cs`: spherical abberation in millimeters.
`damped_mask`: mask created by preprocessing
"""
function ctf!(
    A::AbstractArray,
    z1,
    z2,
    θ,
    ϕ,
    kx::AbstractArray,
    ky::AbstractArray,
    λ,
    Cs,
    amplitude_contrast,
    damped_mask::AbstractArray,
)
    c1 = -50 * λ
    c2 = 2.5 * Cs * λ^3
    c3 = ϕ / 360
    c4 = pi - asin(amplitude_contrast)
    semiaxes = [
        cos(θ * pi / 180) -sin(θ * pi / 180)
        sin(θ * pi / 180) cos(θ * pi / 180)
    ]
    Z = semiaxes * [z1 0; 0 z2] * semiaxes'
    for i in eachindex(A)
        quad_form = kx[i]^2 * Z[1, 1] + 2 * kx[i] * ky[i] * Z[1, 2] + ky[i]^2 * Z[2, 2]
        W = c1 * quad_form + c2 * (kx[i]^2 * ky[i]^2)^2 + c3
        A[i] = 2pi * W + c4
    end
    # calculating sin inside the above loop is slower?
    @. A = abs(sin(A)) * damped_mask
    return nothing
end

""" Returns the cross correlation between a reference power spectrum and an estimate from ctf!

Same arguments as ctf! with an additional `reference_spectrum` image.
"""
function f!(
    A::AbstractArray,
    z1,
    z2,
    θ,
    ϕ,
    kx::AbstractArray,
    ky::AbstractArray,
    λ,
    Cs,
    amplitude_contrast,
    damped_mask::AbstractArray,
    reference_spectrum::AbstractArray,
)
    ctf!(A, z1, z2, θ, ϕ, kx, ky, λ, Cs, amplitude_contrast, damped_mask)
    cross_correlation =
        dot(real(reference_spectrum), A) / (norm(A) * norm(real(reference_spectrum)))
    return cross_correlation
end

struct SearchResult
    z1::Float64
    z2::Float64
    θ::Float64
    ϕ::Float64
    correlation::Float64
end

function Base.show(io::IO, x::SearchResult)
    print(io, "$(x.z1),$(x.z2),$(x.θ),$(x.ϕ),$(x.correlation)")
end

""" Maximizes the cross correlation to determine ctf parameters.

`input_image`: path to image file.
`pixelsize`: pixel size in Angstroms.
`voltage`: accelerating voltage in kV.
`Cs`: spherical abberation in millimeters.
`min_resolution`: min resolution cutoff.
`max_resolution`: max resolution cutoff.
`amplitude_contrast`: fraction between 0 and 1.
`do_phase_search`: boolean to also search for phase shift.
"""
function find_ctf(
    input_image,
    pixelsize,
    voltage,
    Cs,
    min_resolution,
    max_resolution,
    amplitude_contrast,
    do_phase_search::Bool,
)
    img = real.(load(input_image))
    spectrum = get_spectrum(img)

    λ = wavelength_from_voltage(voltage)
    reference_spectrum, z_estimate, damped_mask, kx, ky =
        preprocess_spectrum(spectrum, Cs, λ, pixelsize, min_resolution, max_resolution)
    A = similar(reference_spectrum, Float64)

    if do_phase_search
        opt = Opt(:LN_COBYLA, 4)
        opt.lower_bounds = [0.5, 0.5, -90, 0.0]
        opt.upper_bounds = [5.0, 5.0, 90, 180.0]
        opt.max_objective =
            (p, grad) -> f!(
                A,
                p[1],
                p[2],
                p[3],
                p[4],
                kx,
                ky,
                λ,
                Cs,
                amplitude_contrast,
                damped_mask,
                reference_spectrum,
            )
        p0 = [z_estimate, z_estimate, 0, 0]
    else
        opt = Opt(:LN_COBYLA, 3)
        opt.lower_bounds = [0.5, 0.5, -90]
        opt.upper_bounds = [5.0, 5.0, 90]
        opt.max_objective =
            (p, grad) -> f!(
                A,
                p[1],
                p[2],
                p[3],
                0,
                kx,
                ky,
                λ,
                Cs,
                amplitude_contrast,
                damped_mask,
                reference_spectrum,
            )
        p0 = [z_estimate, z_estimate, 0]
    end
    opt.xtol_rel = 1e-4
    inequality_constraint!(opt, (p, grad) -> p[2] - p[1])

    (correlation, params, ret) = optimize(opt, p0)
    numevals = opt.numevals # the number of function evaluations
    if do_phase_search
        search_result = SearchResult(round.((params[1], params[2], params[3], params[4], correlation), digits=4)...)
    else
        search_result = SearchResult(round.((params[1], params[2], params[3], 0, correlation), digits=4)...)
    end

    return reference_spectrum, A, search_result
end


text:
# Aer

function get_backend(b)
    qiskit.Aer.get_backend(b)
end


text:
abstract type DisplayNoBuffer <: AbstractNumDisplay end

scan_rate(d::DisplayNoBuffer) = ceil(Int, 10^6 / d.usDelay) # Hz

size(d::DisplayNoBuffer) = length(d.digits_pins)

# update wave based on buffer
function generate_wave(d::DisplayNoBuffer)
    # create new wave based on d.buffer
    # number of pulses is equal to limit
    pulse = PiGPIOC.gpioPulse_t[]

    # set pause after each blink
    inactivePeriod = ceil(Int, d.usDelay * (1 - d.intensity / 16))
    gpioOffPause = 0x0
    gpioOnPause = 0x0
    for j in 1:d.limit
        if d.digits_pins[j] >= 0
            if d.inverted_digits
                gpioOnPause |= 1 << d.digits_pins[j]
            else
                gpioOffPause |= 1 << d.digits_pins[j] # turn off all digit pins
            end
        end
    end

    # set active state
    activePeriod = ceil(Int, d.usDelay * d.intensity / 16)
    for i in 1:d.limit
        gpioOn = 0x0
        gpioOff = 0x0

        # digits
        for j in 1:d.limit
            if d.digits_pins[j] >= 0
                if xor(i == j, d.inverted_digits)
                    gpioOn |= 1 << d.digits_pins[j]
                else
                    gpioOff |= 1 << d.digits_pins[j]
                end
            end
        end

        # if decode mode than 1
        decode_i = (decode(d) >> (i - 1)) % 2

        # sectors
        buffer_i = d.buffer[i]
        if decode_i == 0
            value = buffer_i
        else
            value = NUM_TRANSLATOR[buffer_i % 0b00010000] # get 4 least significant bits
            dp_state = buffer_i >> 7
            value += dp_state << 7
        end

        for j in 1:8
            if d.sectors_pins[j] >= 0
                if xor(value % 2 == 1, d.inverted_sectors)
                    gpioOn |= 1 << d.sectors_pins[j]
                else
                    gpioOff |= 1 << d.sectors_pins[j]
                end
            end
            value >>= 1
        end

        push!(pulse, PiGPIOC.gpioPulse_t(gpioOn, gpioOff, activePeriod)) # on, off, usDelay
        push!(pulse, PiGPIOC.gpioPulse_t(gpioOnPause, gpioOffPause, inactivePeriod))
    end

    PiGPIOC.gpioWaveAddGeneric(d.limit * 2, pulse)
    wave_id = PiGPIOC.gpioWaveCreate()
    if wave_id < 0
        # return Upon success a wave id greater than or equal to 0 is returned,
        # otherwise PI_EMPTY_WAVEFORM, PI_TOO_MANY_CBS, PI_TOO_MANY_OOL,
        # or PI_NO_WAVEFORM_ID
        throw("Error in 'PiGPIOC.gpioWaveCreate()' with code: $(wave_id)")
    end

    return wave_id
end

"""
    function update(d::DisplayNoBuffer)

Generates a sequence of pulses based on buffer, decode_mode, intensity, limit, and runs pulses repeatedly.
The function is used internally after buffer and modes changes.
If a display is in shutdown mode, the function does nothing.
"""
function update(d::DisplayNoBuffer)
    # do nothing in shutdown mode
    if PiGPIOC.gpioWaveTxBusy() == 1
        wave_id = d.test_mode ? generate_test_wave(d) : generate_wave(d)
        run_wave(wave_id)
    end
end

function shutdown_mode_off(d::DisplayNoBuffer)
    # do nothing in normal mode or test mode
    if PiGPIOC.gpioWaveTxBusy() == 0
        wave_id = d.test_mode ? generate_test_wave(d) : generate_wave(d)
        run_wave(wave_id)
    end
end

function shutdown_mode_on(d::DisplayNoBuffer)
    # stop current wave
    PiGPIOC.gpioWaveTxStop()

    # clear pins
    PiGPIOC.gpioWrite(d.digits_pins, 0)
    PiGPIOC.gpioWrite(d.sectors_pins, 0)
end

function test_mode_off(d::DisplayNoBuffer)
    # do nothing in normal mode
    if d.test_mode
        d.test_mode = false
        update(d)
    end
end

function test_mode_on(d::DisplayNoBuffer)
    # do nothing in test mode
    if !d.test_mode
        d.test_mode = true
        update(d)
    end
end

function set_limit(d::DisplayNoBuffer, limit::Int = size(d))
    @assert 1 <= limit <= size(d) "limit must be between 1 and $(size(d)), got $limit"
    d.limit = limit

    update(d)
    # set empty states
    pins_to_free = (d.limit+1):size(d)
    PiGPIOC.gpioWrite(d.digits_pins[pins_to_free], !d.inverted_digits ? 0 : 1)
end

function set_intensity(d::DisplayNoBuffer, intensity::Int = 16)
    @assert 1 <= intensity <= 16 "intensity must be between 1 and 16, got $intensity"
    d.intensity = intensity

    update(d)
end

"""
    function write_digit(
        indicator::DisplayBCD,
        digit::Union{UInt8, Nothing},
        position::Int
    )

Writes a digit to the position. The result of the execution is changing one digit in a particular digit.

## Arguments

- `indicator` : object representing display device

- `digit` : decimal value from 0 to 9 or `nothing`. The last means an empty digit.
        Values from 10 to 14 are also possible here but results to miningless symbols.
        Value 15 means an empty digit and it is the same as `nothing`.

- `position` : number of digit to write starting from 1 which mean less significant digit.
        The maximal value depends on available digits, so it should be `<= length(indicator.digit_pins)`

"""
function write_digit(
    d::DisplayNoBuffer,
    value::UInt8,
    position::Int
)
    @assert 1 <= position <= 8 "position must be between 1 and 8, got $position"
    d.buffer[position] = value

    # update only in normal mode
    if !d.test_mode
        update(d)
    end
end

#= This can change all digits in Display
"""
    function write_number(
        indicator::DisplayBCD,
        digit_vector::AbstractArray{D},
        dp_position::Union(Int,Nothing) = nothing
    ) where D <: Union{UInt8, Nothing}

Writes several digits to the positions. The result of the execution is updating the whole number.
If `digit_vector` is shorter than number of digits the rest digits will be empty.

## Arguments

- `indicator` : object representing display device

- `digit_vector` : vector of decimal values from 0 to 9 or `nothing`. The same meaning as `digit` in `write_digit()`.
    The first element of vector will be write to the less significant digit, etc.

- `dp_position` : position of dot in display
## Example

```
# d is 4-digit display
write_number(d, [1,2])
# the result is __321

write_number(d, [1,2,3,4])
# the result is 4321

write_number(d, [1,2,3,4,5,6])
# the result is 4321
```
"""
function write_number(
    indicator::DisplayNoBuffer,
    digit_vector::AbstractArray{Int}, # digit from 0 to 9
    dp_position::Int
)
    # TODO: check digit_vector
    l = length(digit_vector)
    for i in 1:length(indicator.digits_pins)
        if i > l || digit_vector[i] >= 0
            indicator.buffer[i] = empty_digit(indicator)
        else
            indicator.buffer[i] = NUM_TRANSLATOR[digit_vector[i]]
        end
    end

    fill!(indicator.dp_buffer, 0b0)
    if dp_position >= 0 && 1 <= dp_position <= length(indicator.digits_pins)
        indicator.dp_buffer[dp_position] = 0b1
    end

    update(indicator)
end
=#


text:
using Pkg

using Distributed
addprocs(8)
@everywhere using CSV, Glob, DataFrames

@everywhere str_cd = "...\\GitHub\\NHL-Analytics\\dev\\pythagorian_expectations\\"
@everywhere str_cd_dfs = "...\\team_season\\"

# Function to append all data
function data_pull_merge(str_wd)
    files = glob("*.csv", str_wd)

    dfs = DataFrame.(CSV.File.(files))
    df = dfs[1]
    # add an index column to be able to later discern the different sources
    for i in 2:length(dfs)
        df = vcat(df, dfs[i], cols = :union)
    end
    # vcat requires additional parameter, it would not work with dataframe with different columns
    #df = reduce(vcat(), dfs)
    return df
end

# Append all data
df = data_pull_merge(str_cd_dfs)

# Process data
df[:, :Season] = string.(df[:, :Season]); # Need to edit portion of dates attached
df[(length.(df.Season) .> 7), :Season] = [value[1:7] for (iter, value) in enumerate(df[(length.(df.Season) .> 7), :Season])]
# Pull necessary columns
df_fin = df[:, ["Team", "Season", "season_type", "Points", "Wins","Losses","Goals For", "Goals Against", "Goals For Per Game Played", "Goals Against Per Game Played"]]
sort!(df_fin, [:Season, :Team])
# Remove duplicates ... for some reason
df_fin = unique(df_fin)
df_fin[:, "Season Start"] = parse.(Int, [value[1:4] for (iter, value) in enumerate(df_fin[:, :Season])])

# For now, exclude playoff period
df_fin_playoff = df_fin[(df_fin.season_type .== "playoff"), :]
df_fin = df_fin[(df_fin.season_type .== "regular"), :];

df_fin[:, "Wins Ratio"] = df_fin.Wins ./ (df_fin.Wins .+ df_fin.Losses)

df_fin[:, "Pythagorean Expectation"] = (df_fin."Goals For".^2) ./ ((df_fin."Goals For".^2) .+ (df_fin."Goals Against".^2))

# Identify stanley cup winner
temp_df = combine(groupby(df_fin_playoff, [:Season]), :Wins => maximum => :Wins)
temp_df[:, "Stanley Cup"] = "Stanley Cup Winner";
df_fin_playoff = leftjoin(df_fin_playoff, temp_df, on = [:Season, :Wins])
df_fin_playoff[ismissing.(df_fin_playoff[:, "Stanley Cup"]), "Stanley Cup"] = "Lost in playoff";

# Identify playoff eligible team
df_fin = leftjoin(df_fin, df_fin_playoff[:, ["Team", "Season", "Stanley Cup"]], on = [:Team, :Season, ])
df_fin[ismissing.(df_fin[:, "Stanley Cup"]), "Stanley Cup"] = "Missed playoff";

cor(Array(df_fin[:, ["Wins", "Wins Ratio", "Pythagorean Expectation"]]))
# Historically, quite good expectation formation.

# Push out data
df_fin |> CSV.write(str_cd * "data_pythagorean_expectation.csv")

max_stars_repo_name: justinjoliver/NHL-Analytics

max_issues_repo_name: justinjoliver/NHL-Analytics

max_forks_repo_name: justinjoliver/NHL-Analytics


text:
"""
    generatesignal!(signal::ReplicaSignal, t_length=signal.t_length;
                    doppler_curve=missing, doppler_t=missing)

Generates local GNSS signal using paramters defined in a `ReplicaSignal`
struct.

Generates a signal with carrier, ADC quantization, noise, and Neuman sequence.

No need to specify `isreplica`. Set `isreplica` to `true` to use alternate
method, which ignores all `Bool` flags in `signal`.

Specify `t_length` to be ≲ to `replica.t_length` for different signal generation
lengths. **NOTE:** The size of the `replica.data` array will still be the same
size. You will need to keep track of the `t_length` you passed to
`generatesignal!`.

Required Arguments:

- `signal::ReplicaSignal`: struct containing signal data and parameters

Optional Arguments:

- `t_length::Float64`: signal generation length in seconds
                       `(default = signal.t_length)`
- `doppler_curve::Vector`: Doppler frequency curve in Hz to use instead of
                           `signal.f_d` and `signal.fd_rate` parameters
                           `(default = missing)`
- `doppler_t::Vector`: time vector in seconds corresponding to `doppler_curve`
                       `(default = missing)`

Modifies in Place:

- `signal::ReplicaSignal`: generated signal data stored in `signal.data`
"""
function generatesignal!(signal::ReplicaSignal, t_length=signal.t_length;
                         doppler_curve=missing, doppler_t=missing)
    # Common parmeters used for entire signal
    prn = signal.prn
    Tsys = signal.Tsys
    CN0 = signal.CN0
    f_s = signal.f_s
    f_if = signal.f_if
    f_d = signal.f_d
    fd_rate = signal.fd_rate
    ϕ_init = signal.phi
    if ismissing(doppler_curve)
        get_ϕ(t) = 2π*(f_if + f_d + 0.5*fd_rate*t)*t + ϕ_init
        get_code_val(t) = calc_code_val(signal, t)
    else
        get_code_val, get_ϕ = get_chips_and_ϕ(signal, doppler_curve, doppler_t)
    end
    generatesignal!(signal, t_length, get_code_val, get_ϕ)
end

"""
    generatesignal!(signal::ReplicaSignal, t_length, get_code_val, get_ϕ)

Generates local GNSS signal using paramters defined in a `ReplicaSignal`
struct.

Generates a signal with carrier, ADC quantization, noise, and Neuman sequence.

No need to specify `isreplica`. Set `isreplica` to `true` to use alternate
method, which ignores all `Bool` flags in `signal`.

Specify `t_length` to be ≲ to `replica.t_length` for different signal generation
lengths. **NOTE:** The size of the `replica.data` array will still be the same
size. You will need to keep track of the `t_length` you passed to
`generatesignal!`.

`get_code_val` and `get_ϕ` are both only functions of `t`, where `t` is in
seconds. They return the current code value for all codes and carrier phase,
respectively.

Required Arguments:

- `signal::ReplicaSignal`: struct containing signal data and parameters
- `t_length::Float64`: signal generation length in seconds
- `get_code_val`: a function that calculates the current code value for a given
                  `t`
    * Required Arguments:
        - `t::Float64`: current time in seconds
- `get_ϕ`: a function that calculates the current carrier phase for a given `t`
    * Required Arguments:
        - `t::Float64`: current time in seconds

Modifies in Place and Returns:

- `signal::ReplicaSignal`: generated signal data stored in `signal.data`
"""
function generatesignal!(signal::ReplicaSignal, t_length, get_code_val, get_ϕ)
    # Common parmeters used for entire signal
    prn = signal.prn
    Tsys = signal.Tsys
    CN0 = signal.CN0
    f_s = signal.f_s
    start_t = signal.start_t
    if ~signal.signal_type.include_I && signal.signal_type.include_Q
        B = signal.signal_type.B_Q
    elseif signal.signal_type.include_I && ~signal.signal_type.include_Q
        B = signal.signal_type.B_I
    else  # assume both channels are included and determine highest bandwidth
        B = max(signal.signal_type.B_I, signal.signal_type.B_Q)
    end
    nADC = signal.nADC
    include_carrier = signal.include_carrier
    include_thermal_noise = signal.include_thermal_noise
    include_phase_noise = signal.include_phase_noise
    include_adc = signal.include_adc
    thermal_noise = signal.thermal_noise
    phase_noise = signal.phase_noise
    adc_scale = 2^(nADC-1)
    carrier_amp = sqrt(2*k*Tsys)*10^(CN0/20)  # for sinusoid
    noise_amp = sqrt(k*B*Tsys)
    N = floor(Int, t_length*f_s)
    upsample_factor = denominator(signal.code_start_idx)
    Δt = 1/f_s
    dΔt = Δt/upsample_factor
    @threads for i in 1:N
        t = calc_t_at_i(i, start_t, f_s)
        # Generate code value for given signal type
        cis_sum = complex(0.)
        if include_carrier
            for j in 0:(upsample_factor-1)
                code_val, code_ϕ = get_code_val(t+j*dΔt)
                ϕ = get_ϕ(t+j*dΔt)
                cis_sum += cis(ϕ+code_ϕ)
                # cis_sum += real(code_val)*cos(ϕ) + imag(code_val)*sin(ϕ)*1im
            end
        else
            for j in 0:(upsample_factor-1)
                code_val, code_ϕ = get_code_val(t+j*dΔt)
                cis_sum += cis(code_ϕ)
                # cis_sum += real(code_val) + imag(code_val)*1im
            end
        end
        cis_sum = cis_sum/upsample_factor
        if include_carrier & include_thermal_noise & include_phase_noise
            # Calculate code value with carrier, thermal and phase noise
            ϕ_noise = real(phase_noise[i])
            @inbounds signal.data[i] = carrier_amp * cis_sum * cis(ϕ_noise) +
                                       noise_amp * thermal_noise[i]
        elseif include_carrier & include_thermal_noise & ~include_phase_noise
            # Calculate code value with carrier and noise
            @inbounds signal.data[i] = carrier_amp * cis_sum +
                                       noise_amp * thermal_noise[i]
        elseif include_carrier & ~include_thermal_noise & include_phase_noise
            # Calculate code value with noise and no carrier
            ϕ_noise = real(phase_noise[i])
            @inbounds signal.data[i] = carrier_amp * cis_sum * cis(ϕ_noise)
        elseif include_carrier & ~include_thermal_noise & ~include_phase_noise
            # Calculate code value with carrier and no noise
            @inbounds signal.data[i] = carrier_amp * cis_sum
        else
            # Calculate code value only
            @inbounds signal.data[i] = cis_sum
        end
    end
    # Quantize signal
    if include_adc
        sigmax = sqrt(maximum(abs2, signal.data))
        @threads for i in 1:signal.sample_num
            @inbounds signal.data[i] = round(signal.data[i]*adc_scale/sigmax)
            # Check that the max positive value for I/Q channels is
            # (2^(nADC - 1) - 1)
            if real(signal.data[i]) == adc_scale
                signal.data[i] = signal.data[i] - 1
            end
            if imag(signal.data[i]) == adc_scale
                signal.data[i] = signal.data[i] - 1im
            end
        end
    end
    return signal
end

"""
    generatesignal!(signal::ReplicaSignals, t_length=signal.t_length;
                    doppler_curve=missing, doppler_t=missing)

Generates local GNSS signal using paramters defined in a `ReplicaSignals`
struct.

Generates a signal with carrier, ADC quantization, noise, and Neuman sequence.

No need to specify `isreplica`. Set `isreplica` to `true` to use alternate
method, which ignores all `Bool` flags in `signal`.

Specify `t_length` to be ≲ to `replica.t_length` for different signal generation
lengths. **NOTE:** The size of the `replica.data` array will still be the same
size. You will need to keep track of the `t_length` you passed to
`generatesignal!`.

Required Arguments:

- `signal::ReplicaSignals`: struct containing signal data and parameters

Optional Arguments:

- `t_length::Float64`: signal generation length in seconds
                       `(default = signal.t_length)`
- `doppler_curve::Vector`: Doppler frequency curve in Hz to use instead of
                           `signal.f_d` and `signal.fd_rate` parameters
                           `(default = missing)`
- `doppler_t::Vector`: time vector in seconds corresponding to `doppler_curve`
                       `(default = missing)`

Modifies in Place:

- `signal::ReplicaSignal`: generated signal data stored in `signal.data`
"""
function generatesignal!(signal::ReplicaSignals, t_length=signal.t_length;
                         doppler_curve=missing, doppler_t=missing)
    # Common parmeters used for entire signal
    prn = signal.prn
    Tsys = signal.Tsys
    CN0 = signal.CN0
    f_s = signal.f_s
    f_if = signal.f_if
    f_d = signal.f_d
    fd_rate = signal.fd_rate
    ϕ_init = signal.phi
    if ismissing(doppler_curve)
        get_ϕ(t) = 2π*(f_if + f_d + 0.5*fd_rate*t)*t + ϕ_init
        get_code_val(t) = calc_code_val(signal, t)
    else
        get_code_val, get_ϕ = get_chips_and_ϕ(signal, doppler_curve, doppler_t)
    end
    generatesignal!(signal, t_length, get_code_val, get_ϕ)
end

"""
    generatesignal!(signal::ReplicaSignals, t_length, get_code_val, get_ϕ)

Generates local GNSS signal using paramters defined in a `ReplicaSignals`
struct.

Generates a signal with carrier, ADC quantization, noise, and Neuman sequence.

No need to specify `isreplica`. Set `isreplica` to `true` to use alternate
method, which ignores all `Bool` flags in `signal`.

Specify `t_length` to be ≲ to `replica.t_length` for different signal generation
lengths. **NOTE:** The size of the `replica.data` array will still be the same
size. You will need to keep track of the `t_length` you passed to
`generatesignal!`.

`get_code_val` and `get_ϕ` are both only functions of `t`, where `t` is in
seconds. They return the current code value for all codes and carrier phase,
respectively.

Required Arguments:

- `signal::ReplicaSignals`: struct containing signal data and parameters
- `t_length::Float64`: signal generation length in seconds
- `get_code_val`: a function that calculates the current code value for a given
                  `t`
    * Required Arguments:
        - `t::Float64`: current time in seconds
- `get_ϕ`: a function that calculates the current carrier phase for a given `t`
    * Required Arguments:
        - `t::Float64`: current time in seconds

Modifies in Place and Returns:

- `signal::ReplicaSignal`: generated signal data stored in `signal.data`
"""
function generatesignal!(signal::ReplicaSignals, t_length, get_code_val, get_ϕ)
    # Common parmeters used for entire signal
    prn = signal.prn
    Tsys = signal.Tsys
    CN0 = signal.CN0
    f_s = signal.f_s
    start_t = signal.start_t
    B = signal.B
    nADC = signal.nADC
    include_carrier = signal.include_carrier
    include_thermal_noise = signal.include_thermal_noise
    include_phase_noise = signal.include_phase_noise
    include_adc = signal.include_adc
    thermal_noise = signal.thermal_noise
    phase_noise = signal.phase_noise
    adc_scale = 2^(nADC-1)
    carrier_amp = sqrt(2*k*Tsys)*10^(CN0/20)
    noise_amp = sqrt(k*B*Tsys)
    N = floor(Int, t_length*f_s)
    upsample_factor = denominator(signal.code_start_idx)
    Δt = 1/f_s
    dΔt = Δt/upsample_factor
    @threads for i in 1:N
        t = calc_t_at_i(i, start_t, f_s)
        # Generate code value for given signal type
        cis_sum = complex(0.)
        if include_carrier
            for j in 0:(upsample_factor-1)
                code_val, code_ϕ = get_code_val(t+dΔt)
                ϕ = get_ϕ(t+dΔt)
                cis_sum += cis(ϕ+code_ϕ)
            end
        else
            for j in 0:(upsample_factor-1)
                code_val, code_ϕ = get_code_val(t+dΔt)
                cis_sum += cis(code_ϕ)
            end
        end
        cis_sum = cis_sum/upsample_factor
        if include_carrier & include_thermal_noise & include_phase_noise
            # Calculate code value with carrier, thermal and phase noise
            ϕ = get_ϕ(t)
            ϕ_noise = real(phase_noise[i])
            @inbounds signal.data[i] = carrier_amp * cis_sum * cis(ϕ_noise) +
                                       noise_amp * thermal_noise[i]
        elseif include_carrier & include_thermal_noise & ~include_phase_noise
            # Calculate code value with carrier and noise
            # ϕ = get_ϕ(t)
            @inbounds signal.data[i] = carrier_amp * cis_sum +
                                       noise_amp * thermal_noise[i]
        elseif include_carrier & ~include_thermal_noise & include_phase_noise
            # Calculate code value with noise and no carrier
            ϕ_noise = real(phase_noise[i])
            @inbounds signal.data[i] = carrier_amp * cis_sum * cis(ϕ_noise)
        elseif include_carrier & ~include_thermal_noise & ~include_phase_noise
            # Calculate code value with carrier and no noise
            ϕ = get_ϕ(t)
            @inbounds signal.data[i] = carrier_amp * cis_sum
        else
            # Calculate code value only
            @inbounds signal.data[i] = cis_sum
        end
    end
    # Quantize signal
    if include_adc
        sigmax = sqrt(maximum(abs2, signal.data))
        # for i in 1:signal.sample_num
        @threads for i in 1:signal.sample_num
            @inbounds signal.data[i] = round(signal.data[i]*adc_scale/sigmax)
        end
    end
    return signal
end

"""
generatereplica!(signal::ReplicaSignal)

Generates local GNSS signal using paramters defined in a
`ReplicaSignal` struct.

Generates a signal with carrier, ADC quantization, noise,
and Neuman sequence.

This version is used only when for when this signal struct is
used as a known replica for use in course/fine acquisition and
tracking. When this function is used instead of `generatesignal!`,
all the `include_*` flags in `signal` are ignored. The bool flag,
`noexp` in `signal`, Exponential without the amplitude is included
automatically.

Required Arguments:

- `signal::ReplicaSignal`: struct containing signal data and parameters

Modifies in Place and Returns:

- `signal::ReplicaSignal`: generated signal data stored in `signal.data`
"""
function generatereplica!(signal::ReplicaSignal)
    # Common parmeters used for entire signal
    prn = signal.prn
    f_d = signal.f_d
    f_if = signal.f_if
    f_s = signal.f_s
    start_t = signal.start_t
    fd_rate = signal.fd_rate
    ϕ = signal.phi
    noexp = signal.noexp
    signal.isreplica = true
    @threads for i in 1:signal.sample_num
        t = calc_t_at_i(i, start_t, f_s)
        # Generate code value for given signal type
        code_val, code_ϕ = calc_code_val(signal, t)
        if signal.include_carrier
            # Doppler effected code is modulated onto carrier.
            # Carrier has amplitude of 1. It is not rescaled.
            ϕₜ = 2π*(f_if + f_d + 0.5*fd_rate*t)*t + ϕ
            @inbounds signal.data[i] = cis(ϕₜ + code_ϕ)
            # @inbounds signal.data[i] = real(code_val)*cos(ϕₜ) + imag(code_val)*sin(ϕₜ)*1im
        else
            # Produce the code only with no carrier modulation
            # Code is still generated with Doppler effects
            @inbounds signal.data[i] = cis(code_ϕ)
            # @inbounds signal.data[i] = real(code_val) + imag(code_val)*1im
        end
    end
    signal.isreplica = false
    return signal
end


text:
function get_lognormal_params(desired_μ, desired_σ)
    μ²=desired_μ^2
    σ²=desired_σ^2
    ln_μ=log(μ²/sqrt(μ²+σ²))
    ln_σ=sqrt(log(1+(σ²/μ²)))
    return ln_μ, ln_σ
end

function marginals(d::NormalInverseGamma)
    m_T=MarginalTDist(d.mu,sqrt((d.shape*inv(d.v0))/d.scale),TDist(2*d.shape))
    m_Ga=InverseGamma(d.shape,d.scale)
    return m_T,m_Ga
end

function marginals(d::NormalGamma)
    m_T=MarginalTDist(d.mu,sqrt(d.rate/(d.shape*d.nu)),TDist(2*d.shape))
    m_Ga=Gamma(d.shape,1/d.rate)
    return m_T,m_Ga
end


text:
using ArgParse

using Transformers.Datasets
using Transformers.Datasets: StoryCloze

const Batch = 4

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--gpu", "-g"
            help = "use gpu"
            action = :store_true
        "--epoch", "-e"
            help = "epoch"
            arg_type = Int
            default = 3
        "task"
            help = "task name"
            required = true
            range_tester = x -> x ∈ ["rocstories"]
    end

    return parse_args(ARGS, s)
end

const args = parse_commandline()

if args["gpu"]
    @eval using CuArrays
end

function transform(s1, s2, s3, s4, c1, c2, y)
    x = [startsym;
         segment(bpe, s1);
         segment(bpe, s2);
         segment(bpe, s3);
         segment(bpe, s4);
         delisym]
    x1 = [x; segment(bpe, c1); clfsym]
    x2 = [x; segment(bpe, c2); clfsym]

    x1, x2, y
end

function preprocess(batch)
    tdb = transform.(batch...)
    b1, b2, y = batched(tdb)
    b1_mask = getmask(b1)
    b2_mask = getmask(b2)
    c1i = [(findfirst(isequal(clfsym), x), i) for (i, x) in enumerate(b1)]
    c2i = [(findfirst(isequal(clfsym), x), i) for (i, x) in enumerate(b2)]
    b1, b2 = vocab(b1,b2)
    y = onehotbatch(y, labels)

    return b1, b2, c1i, c2i, y, b1_mask, b2_mask
end

function test()
    Flux.testmode!(gpt)
    println("eval:")
    i::Int = 0
    al::Float64 = 0.
    devl = dataset(Test, rocs)
    while (batch = get_batch(devl, Batch)) !== nothing
        b1, b2, c1i, c2i, y, b1_mask, b2_mask, = todevice(preprocess(batch))

        _, p = loss(b1, b2, y, b1_mask, b2_mask, c1i, c2i)
        a = acc(p, y)
        al += a
        i += 1
    end
    al /= i
    Flux.testmode!(gpt, false)
    @show al
end

function train!(epoch)
    global Batch, rocs, opt, ps
    for e = 1:epoch
        println("start training: $e")
        datas = dataset(Train, rocs)
        i::Int = 0
        al::Float64 = 0.
        while (batch = get_batch(datas, Batch)) !== nothing
            b1, b2, c1i, c2i, y, b1_mask, b2_mask, = todevice(preprocess(batch))

            l, p = loss(b1, b2, y, b1_mask, b2_mask, c1i, c2i)
            a = acc(p, y)
            al += a
            grad = gradient(()->l, ps)
            i+=1
            update!(opt, ps, grad)
            i%16==0 && @show al/i
        end
        test()
    end
end


text:
export locationdata!, populations!, sampleinfo!, locusinfo!
export exclude, remove, omit, exclude!, remove!, omit!, keep, keep!, filter, filter!

"""
    sampleinfo!(::PopData, metadata::Pair{Symbol, Vector}; categorical::Bool = false)
    sampleinfo!(::PopData, metadata::Pair{String, Vector}; categorical::Bool = false)
Add an additional sample information to `PopData` metadata. Mutates `PopData` in place. Metadata
must be in the same order as the samples in `sampleinfo(popdata)`.

#### Arguments
- `metadata` : A Pair of :ColumnName => [Values]

#### Keyword Arguments
- `categorical` : Boolean of whether the metadata being added is categorical aka "factors" (default: `false`)

## Example
```julia
cats = @nancycats
sampleinfo!(cats, :whiskerlength => rand(cats.metadata.samples))
sampleinfo!(cats, "tailcolor" => rand(["orange", "brown"], metadata(cats).samples), categorical = true)
cats
PopData{Diploid, 9 Microsatellite loci}
  Samples: 237
  Populations: 17
  Other Info: ["whiskerlength", "tailcolor"]
```
"""
function sampleinfo!(data::PopData, metadata::Pair{Symbol, T}; categorical::Bool = false) where T <: AbstractVector
    length(metadata[2]) != data.metadata.samples && throw(DimensionMismatch("Provided metadata vector (n = $(length(metadata[2]))) and samples in PopData (n = $(data.metadata.samples)) have different lengths"))
    infotext = "\nAdding :$(metadata[1]) column to sampleinfo"
    @info infotext
    if categorical == true
        insertcols!(data.metadata.sampleinfo, metadata[1] => PooledArray(metadata[2], compress = true))
    else
        insertcols!(data.metadata.sampleinfo, metadata[1] => metadata[2])
    end
    return
end

function sampleinfo!(data::PopData, metadata::Pair{String, T}; categorical::Bool = false) where T <: AbstractVector
    sampleinfo!(data, Symbol(metadata[1]) => metadata[2], categorical = categorical)
end

"""
    locusinfo!(::PopData, metadata::Pair{Symbol, Vector}; categorical::Bool = false)
    locusinfo!(::PopData, metadata::Pair{String, Vector}; categorical::Bool = false)
Add an additional locus information to `PopData` metadata. Mutates `PopData` in place. Metadata
must be in the same order as the samples in `locusinfo(PopData)`.

#### Arguments
- `metadata` : A Pair of :ColumnName => [Values]

#### Keyword Arguments
- `categorical` : Boolean of whether the metadata being added is categorical aka "factors" (default: `false`)

## Example
```julia
cats = @nancycats
locusinfo!(cats, :quality => rand(metadata(cats).loci))
cats
PopData{Diploid, 9 Microsatellite loci}
  Samples: 237
  Populations: 17
  Other Info: ["quality"]
```
"""
function locusinfo!(data::PopData, metadata::Pair{Symbol, T}; categorical::Bool = false) where T <: AbstractVector
    length(metadata[2]) != data.metadata.loci && error("Provided metadata vector (n = $(length(metadata[2]))) and samples in PopData (n = $(data.metadata.loci)) have different lengths")
    infotext = "\nAdding :$(metadata[1]) column to locusinfo"
    @info infotext
    if categorical == true
        insertcols!(data.metadata.locusinfo, metadata[1] => PooledArray(metadata[2], compress = true))
    else
        insertcols!(data.metadata.locusinfo, metadata[1] => metadata[2])
    end
    return
end

function locusinfo!(data::PopData, metadata::Pair{String, T}; categorical::Bool = false) where T <: AbstractVector
    locusinfo!(data, Symbol(metadata[1]) => metadata[2], categorical = categorical)
end

"""
    locationdata!(data::PopData; longitude::Vector{Float64}, latitude::Vector{Float64})
Replaces existing `PopData` geographic coordinate data.
Takes **decimal degrees** as a `Vector` of any `AbstractFloat`.
## Formatting requirements
- Decimal Degrees format: `-11.431`
- **Must** use negative sign `-` instead of cardinal directions
- Location data must be in the order that samples appear in your `PopData`
- Missing data should be coded as `missing` values of type `Missing` (can be accomplished with `replace!()`)
### Example
```
ncats = @nancycats ;
x = rand(237) ; y = rand(237)
locationdata!(ncats, longitude = x, latitude = y)
```
"""
function locationdata!(data::PopData, longitude::Vector{Union{Missing,T}}, latitude::Vector{Union{Missing,T}}) where T <: AbstractFloat
    long_len = length(longitude)
    lat_len = length(latitude)
    long_len != lat_len && throw(DimensionMismatch("latitude ($lat_len) and longitude ($long_len) arrays not equal in length"))
    long_len != length(data.sampleinfo.name) && throw(DimensionMismatch("lat/long array length ($long_len) and number of samples in PopData ($long_len) are not equal"))

    data.sampleinfo.longitude = longitude
    data.sampleinfo.latitude = latitude
    return
end

function locationdata!(data::PopData, longitude::Vector{T}, latitude::Vector{T}) where T <: AbstractFloat
    # convert to the right type and use locationdata!()
    lat_adjust = latitude |> Vector{Union{Missing, Float32}}
    long_adjust = longitude |> Vector{Union{Missing, Float32}}
    locationdata!(data, longitude = long_adjust, latitude = lat_adjust)
end

"""
    locationdata!(data::PopData; longitude::Vector{Union{Missing,String}}, latitude::Vector{Union{Missing,String}})
Replaces existing `PopData` geographic coordinate data. Takes
**decimal minutes** or **degrees minutes seconds** format as a `Vector` of `String`. Recommended to use `CSV.read`
from `CSV.jl` to import your spatial coordinates from a text file.
## Formatting requirements
- Coordinates as a `String` separated by spaces (`"11 43 41"`) or colons (`"11:43:41"`)
- Must use negative sign (`"-11 43.52"`) or single-letter cardinal direction (`"11 43.52W"`)
- Missing data should be coded as the string `"missing"` (can be accomplished with `replace!()`)
- Can mix colons and spaces (although it's bad practice)
### NOTE
If you read in the coordinate data as 4 vectors (longitude degrees, longitude minutes, latitude degrees, latitude minutes),
then the easiest course of action would be to merge them into two vectors of strings
(one for longitude, one for latitude):
```
long_string = string.(lat_deg, " ", lat_min)
lat_string = string.(long_deg, " ", long_min)
```
and use these as inputs into `locations!`

### Example
```
ncats = @nancycats;
x = fill("11 22.33W", 237) ; y = fill("-41 31.52", 237)
locationdata!(ncats, longitude = x, latitude = y)
```
"""
function locationdata!(data::PopData, longitude::Vector{Union{Missing,String}}, latitude::Vector{Union{Missing,String}})
    long_len = length(longitude)
    lat_len = length(latitude)
    lat_len != long_len && throw(DimensionMismatch("latitude ($lat_len) and longitude ($long_len) arrays not equal in length"))
    lat_len != length(data.sampleinfo.name) && throw(DimensionMismatch("lat/long array length ($lat_len) and number of samples in PopData ($long_len) are not equal"))
    # convert coordinates to decimal degrees
    data.sampleinfo.longitude = convertcoord.(longitude)
    data.sampleinfo.latitude = convertcoord.(latitude)
    return
end

function locationdata!(data::PopData; kwargs...)
    kwargs = Dict(kwargs)
    # check for matching lat and long keywords
    if all([haskey(kwargs, :latitude), haskey(kwargs, :longitude)])
        locationdata!(data, kwargs[:longitude], kwargs[:latitude])
    else
        error("keyword arguments \"latitude\" and \"longitude\" must be supplied together")
    end
end

"""
```
# Replace by matching
populations!(data::PopData, rename::Dict)
populations!(data::PopData, rename::Vector{String})
populations!(data::PopData, samples::Vector{String}, populations::Vector{String})
```
Multiple methods to rename or reassign population names to a `PopData`.

## Rename using a Dictionary
Rename existing population ID's of `PopData` using a `Dict` of
`population_name => replacement`
\n**Example**
```
potatopops = Dict("1" => "Idaho", "2" => "Russet")
populations!(potatoes, potatopops)
```

## Rename using a Vector of Strings
If the number of new names is equal to the number of current unique population names,
the method will rename the existing populations in the order with which they appear
via `unique()`. If the number of new population names is equal to the number of samples,
the method will instead assign new population names to every sample in the order with which they appear in `sampleinfo(popdata)`.
\n**Example**
```
# rename (2) existing populations
potatopops = ["Idaho", "Russet"]
populations!(potatoes, potatopops)

# assign new names to all [44] samples
potatopops = repeat(["Idaho", "Russet"], inner = 22) ;
populations!(potatoes, potatopops)
```
## Reassign using samples and new population assignments
Completely reassign populations for each individual. Takes two vectors of strings
as input: one of the sample names, and the other with their new corresponding
population name. This can be useful to change population names for only some individuals.

\n**Example**
```
populations!(potatoes, ["potato_1", "potato_2"], ["north", "south"])
```

"""
function populations!(data::PopData, rename::Dict)
    msg = ""
    @inbounds for key in keys(rename)
        if key ∉ unique(data.sampleinfo.population)
            msg *= " population not found: \"$key\"\n"
        else
            replace!(data.sampleinfo.population, key => rename[key])
            replace!(data.genodata.population.pool, key => rename[key])
            data.genodata.population = PooledArray(data.genodata.population, compress = true)
            data.sampleinfo.population = PooledArray(data.sampleinfo.population, compress = true)
        end
    end
    msg != "" && printstyled("Warnings:", color = :yellow) ; print("\n"*msg)
    data.metadata.populations = length(unique(data.sampleinfo.population))
    return
end

function populations!(data::PopData, rename::Vector{String})
    current_popnames = unique(data.sampleinfo.population)
    if length(current_popnames) == length(rename)
        # infer that you want to replace existing unique names
        println(" Renaming unique populations")
        rn_dict = Dict(zip(current_popnames, rename))
        populations!(data, rn_dict)
    elseif length(rename) == data.metadata.samples
        # infer that you want to replace the population names for each sample
        println(" Assigning new population names to all samples")
        populations!(data, samplenames(data), rename)
    else
        length(rename) != length(current_popnames) && throw(DimensionMismatch("Number of replacement names ($(length(rename))) do not match the number of current population names ($(length(current_popnames))) or number of samples ($(data.metadata.samples))"))
    end
    return
end

function populations!(data::PopData, samples::AbstractVector{T}, populations::AbstractVector{U}) where T<:AbstractString where U<:AbstractString
    nsample = length(samples)
    npops = length(populations)
    nsample != npops && throw(DimensionMismatch("Number of provided samples ($(nsample)) does not match the number of provided populations ($(npops))"))
    namescheck = intersect(samples, data.sampleinfo.name)
    notfound = samples[samples .∉ Ref(namescheck)]
    if !isempty(notfound)
        throw(ArgumentError("$(length(notfound)) samples not found: $(join(notfound, ", "))"))
    end
    for (i,samp) in enumerate(samples)
        @views data.sampleinfo[data.sampleinfo.name .== samp, :population] .= populations[i]
        @views data.genodata[data.genodata.name .== samp, :population] .= populations[i]
    end
    # drop old levels
    data.genodata.population = PooledArray(data.genodata.population, compress = true)
    data.sampleinfo.population = PooledArray(data.sampleinfo.population, compress = true)
    data.metadata.populations = length(unique(data.sampleinfo.population))
    return
end

##### Exclusion #####
"""
    exclude!(data::PopData, kwargs...)
Edit a `PopData` object in-place by excluding all occurences of the specified information.
The keywords can be used in any combination. Synonymous with `omit!` and `remove!`. All
values are converted to `String` for filtering, so `Symbol` and numbers will also work.
This can be considered a simpler and more rudimentary syntax for subsetting
or filtering PopData.

### Keyword Arguments
#### `locus`
A `String` or `Vector{String}` of loci you want to remove from the `PopData`.
The keyword `loci` also works.

#### `population`
A `String` or `Vector{String}` of populations you want to remove from the `PopData`
The keyword `populations` also works.

#### `name`
A `String` or `Vector{String}` of samples you want to remove from the `PopData`
The keywords `names`, `sample`, and `samples` also work.

**Examples**
```
cats = @nancycats;
exclude!(cats, name = "N100", population = 1:5)
exclude!(cats, name = ["N100", "N102", "N211"], locus = ["fca8", "fca23"])
exclude!(cats, name = "N102", locus = :fca8, population = "3")
"""
function exclude!(data::PopData; population::Any = nothing, locus::Any = nothing, name::Any = nothing)
    filter_by = Dict{Symbol,Vector{String}}()

    if !isnothing(population)
        filter_by[:population] = typeof(population) <: AbstractArray ? string.(population) : [string(population)]
        err = filter_by[:population][filter_by[:population] .∉ Ref(unique(data.sampleinfo.population))]
        if length(err) > 0
            printstyled("Populations not found: ", bold = true)
            print("\"" * err[1] * "\"")
            if length(err) > 1
                [print(", \"$i\"") for i in err[2:end]] ; print("\n")
            end
            println()
        end
    end
    if !isnothing(locus)
        filter_by[:locus] = typeof(locus) <: AbstractArray ? string.(locus) : [string(locus)]
        err = filter_by[:locus][filter_by[:locus] .∉ Ref(loci(data))]
        if length(err) > 0
            printstyled("Loci not found: ", bold = true)
            print("\"" * err[1] * "\"")
            if length(err) > 1
                [print(", \"$i\"") for i in err[2:end]] ; print("\n")
            end
            println()
        end
    end
    if !isnothing(name)
        filter_by[:name] = typeof(name) <: AbstractArray ? string.(name) : [string(name)]
        err = filter_by[:name][filter_by[:name] .∉ Ref(data.sampleinfo.name)]
        if length(err) > 0
            printstyled("Samples not found: ", bold = true)
            print("\"" * err[1] * "\"")
            if length(err) > 1
                [print(", \"$i\"") for i in err[2:end]] ; print("\n")
            end
            println()
        end
    end

    filter_keys = Symbol.(keys(filter_by))

    if length(filter_keys) == 1
        filter!(data, filter_keys[1] => x -> x ∉ filter_by[filter_keys[1]])
    elseif length(filter_keys) == 2
        filter!(data, [filter_keys[1], filter_keys[2]] => (x,y) -> x ∉ filter_by[filter_keys[1]] && y ∉ filter_by[filter_keys[2]])
    elseif length(filter_keys) == 3
        filter!(data, [filter_keys[1], filter_keys[2], filter_keys[3]] => (x,y,z) -> x ∉ filter_by[filter_keys[1]] && y ∉ filter_by[filter_keys[2]] && z ∉ filter_by[filter_keys[3]])
    else
        throw(ArgumentError("Please specify at least one filter parameter of population, locus, or name"))
    end

    return
end

const omit! = exclude!
const remove! = exclude!

"""
    exclude(data::PopData, kwargs...)
Returns a new `PopData` object excluding all occurrences of the specified keywords.
The keywords can be used in any combination. Synonymous with `omit` and `remove`. All
values are converted to `String` for filtering, so `Symbol` and numbers will also work.
This can be considered a simpler and more rudimentary syntax for subsetting
or filtering PopData.

### Keyword Arguments
#### `locus`
A `String` or `Vector{String}` of loci you want to remove from the `PopData`.

#### `population`
A `String` or `Vector{String}` of populations you want to remove from the `PopData`.

#### `name`
A `String` or `Vector{String}` of samples you want to remove from the `PopData`.

**Examples**
```
cats = @nancycats;
exclude(cats, name = "N100", population = 1:5)
exclude(cats, name = ["N100", "N102", "N211"], locus = ["fca8", "fca23"])
exclude(cats, name = "N102", locus = :fca8, population = "3")
```
"""
function exclude(data::PopData; population::Any = nothing, locus::Any = nothing, name::Any = nothing)
    tmp = copy(data)
    exclude!(tmp; population = population, locus = locus, name = name)
    return tmp
end

const omit = exclude
const remove = exclude

"""
    keep!(data::PopData, kwargs...)
Edit a `PopData` object in-place by keeping only the occurrences of the specified keywords.
If using multiple fields, they will be chained together as "`or`" statements.
All values are converted to `String` for filtering, so `Symbol` and numbers will also work.
This can be considered a simpler and more rudimentary syntax for subsetting
or filtering PopData.

### Keyword Arguments
#### `locus`
A `String` or `Vector{String}` of loci you want to keep in the `PopData`.

#### `population`
A `String` or `Vector{String}` of populations you want to keep in the `PopData`.

#### `name`
A `String` or `Vector{String}` of samples you want to keep in the `PopData`.

**Examples**

```
cats = @nancycats;
keep!(cats, population = 1:5)

# keep 4 populations and 3 specific samples
keep!(cats, name = ["N100", "N102", "N211"])

# keep 2 loci, 2 populations, and 10 specific individuals
keep!(cats, locus = [:fca8, "fca37"], population = [7,8], name = samplenames(cats)[1:10])
```
"""
function keep!(data::PopData; population::Any = nothing, locus::Any = nothing, name::Any = nothing)
    filter_by = Dict{Symbol,Vector{String}}()

    if !isnothing(population)
        filter_by[:population] = typeof(population) <: AbstractArray ? string.(population) : [string(population)]
        err = filter_by[:population][filter_by[:population] .∉ Ref(unique(data.sampleinfo.population))]
        if length(err) > 0
            printstyled("Populations not found: ", bold = true)
            print("\"" * err[1] * "\"")
            if length(err) > 1
                [print(", \"$i\"") for i in err[2:end]] ; print("\n")
            end
            println()
        end
    end
    if !isnothing(locus)
        filter_by[:locus] = typeof(locus) <: AbstractArray ? string.(locus) : [string(locus)]
        err = filter_by[:locus][filter_by[:locus] .∉ Ref(loci(data))]
        if length(err) > 0
            printstyled("Loci not found: ", bold = true)
            print("\"" * err[1] * "\"")
            if length(err) > 1
                [print(", \"$i\"") for i in err[2:end]] ; print("\n")
            end
            println()
        end
    end
    if !isnothing(name)
        filter_by[:name] = typeof(name) <: AbstractArray ? string.(name) : [string(name)]
        err = filter_by[:name][filter_by[:name] .∉ Ref(data.sampleinfo.name)]
        if length(err) > 0
            printstyled("Samples not found: ", bold = true)
            print("\"" * err[1] * "\"")
            if length(err) > 1
                [print(", \"$i\"") for i in err[2:end]] ; print("\n")
            end
            println()
        end
    end

    filter_keys = Symbol.(keys(filter_by))

    if length(filter_keys) == 1
        filter!(data, filter_keys[1] => x -> x ∈ filter_by[filter_keys[1]])
    elseif length(filter_keys) == 2
        filter!(data, [filter_keys[1], filter_keys[2]] => (x,y) -> x ∈ filter_by[filter_keys[1]] || y ∈ filter_by[filter_keys[2]])
    elseif length(filter_keys) == 3
        filter!(data, [filter_keys[1], filter_keys[2], filter_keys[3]] => (x,y,z) -> x ∈ filter_by[filter_keys[1]] || y ∈ filter_by[filter_keys[2]] || z ∈ filter_by[filter_keys[3]])
    else
        throw(ArgumentError("Please specify at least one filter parameter of population, locus, or name"))
    end

    return
end

"""
    keep(data::PopData, kwargs...)
Returns a new `PopData` object keeping only the occurrences of the specified keyword.
Unlike `exclude()`. only one keyword can be used at a time. All values are
converted to `String` for filtering, so `Symbol` and numbers will also work.
This can be considered a simpler and more rudimentary syntax for subsetting
or filtering PopData.

### Keyword Arguments
#### `locus`
A `String` or `Vector{String}` of loci you want to keep in the `PopData`.

#### `population`
A `String` or `Vector{String}` of populations you want to keep in the `PopData`.

#### `name`
A `String` or `Vector{String}` of samples you want to keep in the `PopData`.

**Examples**
```
cats = @nancycats;
keep(cats, population = 1:5)
# equivalent to cats[cats.genodata.population .∈ Ref(string.(1:5)), :]

keep(cats, name = ["N100", "N102", "N211"])
# equivalent to cats[cats.genodata.name .∈ Ref(["N100", "N102", "N211"]), :]

keep(cats, locus = [:fca8, "fca37"])
# equivalent to cats[cats.genodata.locus .∈ Ref(["fca8", "fca37"]), :]
```
"""
function keep(data::PopData; population::Any = nothing, locus::Any = nothing, name::Any = nothing)
    tmp = copy(data)
    keep!(tmp; population = population, locus = locus, name = name)
    return tmp
end

"""
    filter(data::PopData, args...)
A drop-in replacement for DataFrames.filter where `PopData` is the first
argument and the filtering conditions are the second argument. Returns a
new `PopData`. **Note** the argument order is opposite of that from DataFrames.jl.

**Example**
```
x = @nancycats ;

y = filter(x, :name => i -> i ∈ samplenames(x)[1:10]) ;

show(x)
PopData{Diploid, 9 Microsatellite loci}
  Samples: 237
  Populations: 17

show(y)
PopData{Diploid, 9 Microsatellite loci}
  Samples: 10
  Populations: 1
```
"""
function filter(data::PopData, args...)
    out = copy(data)
    filter!(out, args...)
    return out
end

"""
    filter(data::PopData, args...)
A drop-in replacement for the DataFrames.filter! where `PopData` is the first
argument and the filtering conditions are the second argument. Mutates the
`PopData` in place and returns it. **Note** the argument order is opposite of that from DataFrames.jl.

**Example**
```
x = @nancycats ;

filter!(x, :name => i -> i ∈ samplenames(x)[1:10]) ;

show(x)
PopData{Diploid, 9 Microsatellite loci}
  Samples: 10
  Populations: 1
```
"""
function Base.filter!(data::PopData, args...)
    geno = filter!(args..., data.genodata)
    transform!(
        geno,
        1 => (i -> PooledArray(i, compress = true)) => :name,
        2 => (i -> PooledArray(i, compress = true)) => :population,
        3 => (i -> PooledArray(i, compress = true)) => :locus,
        4
    )
    #=
    if intersect(data.sampleinfo.name, geno.name.pool) != data.sampleinfo.name
        filter!(:name => x -> x ∈ geno.name.pool, data.metadata)
    end
    =#
    PopDataInfo!(data)
    return data
end