text:
module MaxHelpingHandReskinnableSwapBlock

using ..Ahorn, Maple

@pardef ReskinnableSwapBlock(x1::Integer, y1::Integer, x2::Integer=x1+16, y2::Integer=y1, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight,
    spriteDirectory::String="objects/swapblock", particleColor1::String="fbf236", particleColor2::String="6abe30",
    moveSound::String="event:/game/05_mirror_temple/swapblock_move", returnSound::String="event:/game/05_mirror_temple/swapblock_return",
    moveEndSound::String="event:/game/05_mirror_temple/swapblock_move_end", returnEndSound::String="event:/game/05_mirror_temple/swapblock_return_end") =
    Entity("MaxHelpingHand/ReskinnableSwapBlock", x=x1, y=y1, nodes=Tuple{Int, Int}[(x2, y2)], width=width, height=height, spriteDirectory=spriteDirectory, particleColor1=particleColor1, particleColor2=particleColor2,
    moveSound=moveSound, returnSound=returnSound, moveEndSound=moveEndSound, returnEndSound=returnEndSound)

function swapFinalizer(entity)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    entity.data["nodes"] = [(x + width, y)]
end

const placements = Ahorn.PlacementDict(
    "Reskinnable Swap Block (max480's Helping Hand)" => Ahorn.EntityPlacement(
        ReskinnableSwapBlock,
        "rectangle",
        Dict{String, Any}(),
        swapFinalizer
    )
)

Ahorn.editingOrder(entity::ReskinnableSwapBlock) = String["x", "y", "width", "height", "moveSound", "moveEndSound", "returnSound", "returnEndSound"]

Ahorn.editingOptions(entity::ReskinnableSwapBlock) = Dict{String, Any}(
    "spriteDirectory" => String["objects/swapblock", "objects/swapblock/moon"]
)

Ahorn.nodeLimits(entity::ReskinnableSwapBlock) = 1, 1

Ahorn.minimumSize(entity::ReskinnableSwapBlock) = 16, 16
Ahorn.resizable(entity::ReskinnableSwapBlock) = true, true

function Ahorn.selection(entity::ReskinnableSwapBlock)
    x, y = Ahorn.position(entity)
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return [Ahorn.Rectangle(x, y, width, height), Ahorn.Rectangle(stopX, stopY, width, height)]
end

function getTextures(entity::ReskinnableSwapBlock)
    spriteDirectory = get(entity, "spriteDirectory", "objects/swapblock")
    return "$(spriteDirectory)/blockRed", "$(spriteDirectory)/target", "$(spriteDirectory)/midBlockRed00"
end

function renderTrail(ctx, x::Number, y::Number, width::Number, height::Number, trail::String)
    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    offset = (trail == "objects/swapblock/target" || trail == "objects/swapblock/moon/target" ? 0 : 2)

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, trail, x + (i - 1) * 8, y + 2, 6 + offset, 0 + offset, 8, 6)
        Ahorn.drawImage(ctx, trail, x + (i - 1) * 8, y + height - 8, 6 + offset, 14 + offset, 8, 6)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, trail, x + 2, y + (i - 1) * 8, 0 + offset, 6 + offset, 6, 8)
        Ahorn.drawImage(ctx, trail, x + width - 8, y + (i - 1) * 8, 14 + offset, 6 + offset, 6, 8)
    end

    for i in 2:tilesWidth - 1, j in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, trail, x + (i - 1) * 8, y + (j - 1) * 8, 6 + offset, 6 + offset, 8, 8)
    end

    Ahorn.drawImage(ctx, trail, x + width - 8, y + 2, 14 + offset, 0 + offset, 6, 6)
    Ahorn.drawImage(ctx, trail, x + width - 8, y + height - 8, 14 + offset, 14 + offset, 6, 6)
    Ahorn.drawImage(ctx, trail, x + 2, y + 2, 0 + offset, 0 + offset, 6, 6)
    Ahorn.drawImage(ctx, trail, x + 2, y + height - 8, 0 + offset, 14 + offset, 6, 6)
end

function renderSwapBlock(ctx::Ahorn.Cairo.CairoContext, x::Number, y::Number, width::Number, height::Number, midResource::String, frame::String)
    midSprite = Ahorn.getSprite(midResource, "Gameplay")

    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y, 8, 0, 8, 8)
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y + height - 8, 8, 16, 8, 8)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, x, y + (i - 1) * 8, 0, 8, 8, 8)
        Ahorn.drawImage(ctx, frame, x + width - 8, y + (i - 1) * 8, 16, 8, 8, 8)
    end

    for i in 2:tilesWidth - 1, j in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y + (j - 1) * 8, 8, 8, 8, 8)
    end

    Ahorn.drawImage(ctx, frame, x, y, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, x + width - 8, y, 16, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, x, y + height - 8, 0, 16, 8, 8)
    Ahorn.drawImage(ctx, frame, x + width - 8, y + height - 8, 16, 16, 8, 8)

    Ahorn.drawImage(ctx, midSprite, x + div(width - midSprite.width, 2), y + div(height - midSprite.height, 2))
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::ReskinnableSwapBlock, room::Maple.Room)
    sprite = get(entity.data, "sprite", "block")
    startX, startY = Int(entity.data["x"]), Int(entity.data["y"])
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    frame, trail, mid = getTextures(entity)

    renderSwapBlock(ctx, stopX, stopY, width, height, mid, frame)
    Ahorn.drawArrow(ctx, startX + width / 2, startY + height / 2, stopX + width / 2, stopY + height / 2, Ahorn.colors.selection_selected_fc, headLength=6)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::ReskinnableSwapBlock, room::Maple.Room)
    sprite = get(entity.data, "sprite", "block")

    startX, startY = Int(entity.data["x"]), Int(entity.data["y"])
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    frame, trail, mid = getTextures(entity)

    renderTrail(ctx, min(startX, stopX), min(startY, stopY), abs(startX - stopX) + width, abs(startY - stopY) + height, trail)
    renderSwapBlock(ctx, startX, startY, width, height, mid, frame)
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
# Automatically generated using Clang.jl

const CURAND_VER_MAJOR = 10
const CURAND_VER_MINOR = 2
const CURAND_VER_PATCH = 3
const CURAND_VER_BUILD = 152
const CURAND_VERSION = CURAND_VER_MAJOR * 1000 + CURAND_VER_MINOR * 100 + CURAND_VER_PATCH

@cenum curandStatus::UInt32 begin
    CURAND_STATUS_SUCCESS = 0
    CURAND_STATUS_VERSION_MISMATCH = 100
    CURAND_STATUS_NOT_INITIALIZED = 101
    CURAND_STATUS_ALLOCATION_FAILED = 102
    CURAND_STATUS_TYPE_ERROR = 103
    CURAND_STATUS_OUT_OF_RANGE = 104
    CURAND_STATUS_LENGTH_NOT_MULTIPLE = 105
    CURAND_STATUS_DOUBLE_PRECISION_REQUIRED = 106
    CURAND_STATUS_LAUNCH_FAILURE = 201
    CURAND_STATUS_PREEXISTING_FAILURE = 202
    CURAND_STATUS_INITIALIZATION_FAILED = 203
    CURAND_STATUS_ARCH_MISMATCH = 204
    CURAND_STATUS_INTERNAL_ERROR = 999
end

const curandStatus_t = curandStatus

@cenum curandRngType::UInt32 begin
    CURAND_RNG_TEST = 0
    CURAND_RNG_PSEUDO_DEFAULT = 100
    CURAND_RNG_PSEUDO_XORWOW = 101
    CURAND_RNG_PSEUDO_MRG32K3A = 121
    CURAND_RNG_PSEUDO_MTGP32 = 141
    CURAND_RNG_PSEUDO_MT19937 = 142
    CURAND_RNG_PSEUDO_PHILOX4_32_10 = 161
    CURAND_RNG_QUASI_DEFAULT = 200
    CURAND_RNG_QUASI_SOBOL32 = 201
    CURAND_RNG_QUASI_SCRAMBLED_SOBOL32 = 202
    CURAND_RNG_QUASI_SOBOL64 = 203
    CURAND_RNG_QUASI_SCRAMBLED_SOBOL64 = 204
end

const curandRngType_t = curandRngType

@cenum curandOrdering::UInt32 begin
    CURAND_ORDERING_PSEUDO_BEST = 100
    CURAND_ORDERING_PSEUDO_DEFAULT = 101
    CURAND_ORDERING_PSEUDO_SEEDED = 102
    CURAND_ORDERING_PSEUDO_LEGACY = 103
    CURAND_ORDERING_QUASI_DEFAULT = 201
end

const curandOrdering_t = curandOrdering

@cenum curandDirectionVectorSet::UInt32 begin
    CURAND_DIRECTION_VECTORS_32_JOEKUO6 = 101
    CURAND_SCRAMBLED_DIRECTION_VECTORS_32_JOEKUO6 = 102
    CURAND_DIRECTION_VECTORS_64_JOEKUO6 = 103
    CURAND_SCRAMBLED_DIRECTION_VECTORS_64_JOEKUO6 = 104
end

const curandDirectionVectorSet_t = curandDirectionVectorSet
const curandDirectionVectors32_t = NTuple{32, UInt32}
const curandDirectionVectors64_t = NTuple{64, Culonglong}
const curandGenerator_st = Cvoid
const curandGenerator_t = Ptr{curandGenerator_st}
const curandDistribution_st = Cdouble
const curandDistribution_t = Ptr{curandDistribution_st}
const curandDistributionShift_st = Cvoid
const curandDistributionShift_t = Ptr{curandDistributionShift_st}
const curandDistributionM2Shift_st = Cvoid
const curandDistributionM2Shift_t = Ptr{curandDistributionM2Shift_st}
const curandHistogramM2_st = Cvoid
const curandHistogramM2_t = Ptr{curandHistogramM2_st}
const curandHistogramM2K_st = UInt32
const curandHistogramM2K_t = Ptr{curandHistogramM2K_st}
const curandHistogramM2V_st = curandDistribution_st
const curandHistogramM2V_t = Ptr{curandHistogramM2V_st}
const curandDiscreteDistribution_st = Cvoid
const curandDiscreteDistribution_t = Ptr{curandDiscreteDistribution_st}

@cenum curandMethod::UInt32 begin
    CURAND_CHOOSE_BEST = 0
    CURAND_ITR = 1
    CURAND_KNUTH = 2
    CURAND_HITR = 3
    CURAND_M1 = 4
    CURAND_M2 = 5
    CURAND_BINARY_SEARCH = 6
    CURAND_DISCRETE_GAUSS = 7
    CURAND_REJECTION = 8
    CURAND_DEVICE_API = 9
    CURAND_FAST_REJECTION = 10
    CURAND_3RD = 11
    CURAND_DEFINITION = 12
    CURAND_POISSON = 13
end

const curandMethod_t = curandMethod

ext: jl


text:
using DataFrames, CSV
using CategoricalArrays

@testset "MLJ Integration Classification" begin
    # ============== Problem setup ============ #
    headers = [
        "id", "Clump Thickness",
        "Uniformity of Cell Size", "Uniformity of Cell Shape",
        "Marginal Adhesion", "Single Epithelial Cell Size",
        "Bare Nuclei", "Bland Chromatin",
        "Normal Nucleoli", "Mitoses", "class"
    ]

    path = "wbc.csv"
    # Replace the "?" to `missing`
    data = CSV.File(path; header=headers, missingstring="?") |> DataFrame

    # Don't include the id's
    select!(data, Not(:id))

    # Change the class tags
    replace!(data.class, 2 => -1)
    replace!(data.class, 4 => 1)

    # Transform to categorical
    transform!(data, :class => categorical, renamecols=false)

    # We don't need `missing`'s
    data = dropmissing(data)

    # Split the training and test data
    y, X = unpack(data, ==(:class), colname -> true)
    train, test = partition(eachindex(y), 2 / 3, shuffle=true, rng=15)

    # Define a good set of hyperparameters for this problem
    pipe = MLJ.@pipeline(Standardizer(count=true), LSSVClassifier(γ=80.0, σ=0.233333))
    mach = MLJ.machine(pipe, X, y)
    MLJ.fit!(mach, rows=train, verbosity=0)

    results = MLJ.predict(mach, rows=test)
    acc = MLJ.accuracy(results, y[test])

    # Test for correctness
    @test isreal(acc)
    # Test for accuracy, at least 95% for this problem
    @test acc >= 0.95
end

@testset "MLJ Integration Regressor" begin
    X, y = MLJ.make_regression(100, 5; noise=0.5, rng=18)
    df = DataFrame(X)
    df.y = y
    dfnew = coerce(df, autotype(df))

    y, X = unpack(dfnew, ==(:y), colname -> true)
    train, test = partition(eachindex(y), 0.75, shuffle=true, rng=20)

    # Define a good set of hyperparameters for this problem
    model = MLJ.@pipeline(Standardizer(), LSSVRegressor(γ=10.0, σ=0.5))
    mach = MLJ.machine(model, X, y)
    MLJ.fit!(mach, rows=train, verbosity=0)

    ŷ = MLJ.predict(mach, rows=test)
    result = round(MLJ.rms(ŷ, y[test]), sigdigits=4)

    # Test for correctness
    @test isreal(result)
end

@testset "Multiclass classification" begin
    X, y = MLJ.@load_iris
    train, test = partition(eachindex(y), 0.6, shuffle=true, rng=30)
    pipe = MLJ.@pipeline(Standardizer(), LSSVClassifier(γ=80.0, σ=0.1))
    mach = MLJ.machine(pipe, X, y)
    MLJ.fit!(mach, rows=train, verbosity=0)
    results = MLJ.predict(mach, rows=test)
    acc = MLJ.accuracy(results, y[test])

    # Check that it is not NaN, and never zero
    @test isreal(acc) && acc >= 0.9
end

max_stars_repo_name: edwinb-ai/LeastSquaresSVM

max_issues_repo_name: edwinb-ai/Elysivm

max_forks_repo_name: edwinb-ai/LeastSquaresSVM


text:
abstract type AbstractIntensityMap{T,S} <: AbstractMatrix{T} end

#abstract type AbstractPolarizedMap{I,Q,U,V} end

@inline function intensitymap(s::M,
                              fovx::Number, fovy::Number,
                              nx::Int, ny::Int;
                              pulse=ComradeBase.DeltaPulse(),
                              executor=SequentialEx()) where {M}
    return intensitymap(imanalytic(M), s, fovx, fovy, nx, ny; pulse, executor)
end

@inline function intensitymap!(img::AbstractIntensityMap, s::M, executor=SequentialEx()) where {M}
    return intensitymap!(imanalytic(M), img, s, executor)
end

function intensitymap(::IsAnalytic, s,
                      fovx::Number, fovy::Number,
                      nx::Int, ny::Int;
                      pulse=ComradeBase.DeltaPulse(),
                      executor=SequentialEx())
    T = typeof(intensity_point(s, 0.0, 0.0))
    img = IntensityMap(zeros(T, ny, nx), fovx, fovy, pulse)
    intensitymap!(IsAnalytic(), img, s, executor)
    return img
end

function intensitymap!(::IsAnalytic, img::AbstractIntensityMap, s, executor=SequentialEx())
    x,y = imagepixels(img)
    dx, dy = pixelsizes(img)
    @floop executor for I in CartesianIndices(img)
        iy,ix = Tuple(I)
        img[I] = intensity_point(s, x[ix], y[iy])*dx*dy
    end
    return img
end

# function intensitymap(::IsAnalytic, s, fovx::Number, fovy::Number, nx::Int, ny::Int; pulse=ComradeBase.DeltaPulse())
#     x,y = imagepixels(fovx, fovy, nx, ny)
#     pimg = map(CartesianIndices((1:ny,1:nx))) do I
#         iy,ix = Tuple(I)
#         @inbounds f = intensity_point(s, x[ix], y[iy])
#         return f
#     end
#     return IntensityMap(pimg, fovx, fovy, pulse)
# end

# function intensitymap!(::IsAnalytic, im::AbstractIntensityMap, m)
#     xitr, yitr = imagepixels(im)
#     @inbounds for (i,x) in pairs(xitr), (j,y) in pairs(yitr)
#         im[j, i] = intensity_point(m, x, y)
#     end
#     return im
# end

include("pulse.jl")
include("polarizedtypes.jl")
include("intensitymap.jl")
#include("polarizedmap.jl")


text:
module FOCs

export foc_mcost, foc_retail, foc_whsale!, foc_retail!, foc_retailV!
export foc_retailPartial
using TypeDef, Helper, Tools, Compute
using NLsolve, Printf

function foc_whsale!(
  fx,
  wprice0,
  price,
  mcost,
  fail,
  dim ::dimType,
  mat ::matType,
  par ::parType,
  sim ::simType,
  indx::indxType,
  mask::maskType)

  if any(fail)==true
    @printf("Aborting foc_whsale due to previous failure.\n")
    fx .= 0
    return
  end
  wprice = zeros(dim.product)
  wprice[indx.mvno] = wprice0
  p0 = view(price,indx.anyActive)

  # foc2(x) = foc_retail(x,mcost,wprice,dim,mat,par,sim,indx,mask)
  # Op = nlsolve(foc2,p0,show_trace=:true,autodiff=:finiteforward)

  foc3(f,x) = foc_retail!(f,x,mcost,wprice,dim,mat,par,sim,indx,mask)
  f0 = zeros(sum(indx.anyActive))
  df = OnceDifferentiable(foc3,p0,f0,autodiff=:finiteforward)
  Op = nlsolve(df,p0)
  iter = 1
  while (checkConvg(iter,Op)==false) & (iter<10)
    # @printf("Solving foc_retail failed %d times.\n",iter)
    Op = nlsolve(df,p0+2randn(length(p0)))
    iter += 1
  end
  while (checkConvg(iter,Op)==false) & (iter<15)
    @printf("Solving foc_retail failed %d times.\n",iter)
    Op = nlsolve(df,0.5p0 + 2randn(length(p0)))
    iter += 1
  end
  if checkConvg(iter,Op)==false
    @printf("Failed to solve foc_retail within foc_whsale %d times.\n",iter)
    fx .= 0
    fail .= true
    return
  end

  p = zeros(dim.product)
  p[indx.anyActive] = Op.zero

  markupR,markupW,markupV = getMarkup(p,wprice,mcost)
  demand!(sim,p,dim,mat,par,indx)
  x = passThrough(markupR,markupW,markupV,dim,mat,sim,indx,mask)

  f = zeros(dim.mvno)
  for n=1:3
    Ln = (mat.firm[:,2].==n)
    Lc = (mat.firm[:,2].!=n)
    L0 = mat.mvno[:,n]
    n0 = mat.mvnp[:,n]
    dΠdp = (markupR.*Ln)'*sim.∂D∂p + (markupW.*L0)'*sim.∂D∂p
    f[n0] = sim.D[L0] + (dΠdp[:,Lc]*x[Lc,n0])'
  end
  fx .= f[indx.anyActive[mask.mvno]]
end

function foc_mcost(
  mcostinc,
  mcostfree,
  p::Vector{Float64},
  wprice::Vector{Float64},
  dim::dimType,
  mat::matType,
  par::parType,
  sim::simType,
  indx::indxType,
  mask::maskType)

  mcost = mergeMcosts(mcostinc,mcostfree,dim,indx,mask)
  markupR,markupW,markupV = getMarkup(p,wprice,mcost)

  x  = passThrough(markupR,markupW,markupV,dim,mat,sim,indx,mask)
  f  = zeros(dim.product)
  p0 = view(p,indx.anyActive)
  f1 = foc_retail(p0,mcost,wprice,dim,mat,par,sim,indx,mask,"mcost")
  f[mask.mno .& indx.anyActive] = f1

  for n=1:3
    Ln = (mat.firm[:,2].==n)
    Lc = (mat.firm[:,2].!=n)
    L0 = mat.mvno[:,n]
    n0 = mat.mvnp[:,n]
    dΠdp = (markupR.*Ln)' * sim.∂D∂p + (markupW.*L0)' * sim.∂D∂p
    f[L0] = sim.D[L0] + (dΠdp[:,Lc]*x[Lc,n0])'
  end
  return f[indx.anyActive .& .!mask.free]

end

function foc_retail(
  p0,
  mc  ::Vector{Vector{Float64}},
  wp  ::Vector{Float64},
  dim ::dimType,
  mat ::matType,
  par ::parType,
  sim ::simType,
  indx::indxType,
  mask::maskType,
  opt ::String="retail")

  p = zeros(dim.product)
  p[indx.anyActive] = p0

  demandD1only!(sim,p,dim,mat,par,indx)

  markup1 = p - mc[2]
  markup1[mask.mvno] = wp[mask.mvno] - mc[1][mask.mvno]
  markup2 = p - wp - mc[2]

  Dm1 = (sim.∂D∂p.*mat.host)'*(markup1.*indx.anyActive)
  Dm2 = (sim.∂D∂p.*mat.host)'*(markup2.*indx.anyActive)

  f = zeros(dim.product)
  f[mask.mno] = sim.D[mask.mno] + Dm1[mask.mno]
  if opt=="retail"
    f[mask.mvno] = sim.D[mask.mvno] + Dm2[mask.mvno]
    return f[indx.anyActive]
  elseif opt=="mcost"
    return f[mask.mno .& indx.anyActive]
  end
end

function foc_retail!(
  f   ::Vector{Float64},
  p0  ::Vector{Float64},
  mc  ::Vector{Vector{Float64}},
  wp  ::Vector{Float64},
  dim ::dimType,
  mat ::matType,
  par ::parType,
  sim ::simType,
  indx::indxType,
  mask::maskType)

  p = zeros(dim.product)
  p[indx.anyActive] = p0

  demandD1only!(sim,p,dim,mat,par,indx)

  idx = findall(indx.anyActive)
  @inbounds for i=1:length(idx)
    j = idx[i]
    Dm1 = 0.0
    Dm2 = 0.0
    for k in idx
      markup1 = (mask.mno[k]) ? p[k]-mc[2][k] : wp[k]-mc[1][k]
      markup2 = p[k] - wp[k] - mc[2][k]
      Dm1 += sim.∂D∂p[k,j] * mat.host[k,j] * markup1
      Dm2 += sim.∂D∂p[k,j] * mat.host[k,j] * markup2
    end
    f[i] = sim.D[j] + Dm1*mask.mno[j] + Dm2*mask.mvno[j]
  end

end

function foc_retailV!(
  f   ::Vector{Float64},
  p0  ::Vector{Float64},
  mc  ::Vector{Vector{Float64}},
  dim ::dimType,
  mat ::matType,
  par ::parType,
  sim ::simType,
  indx::indxType,
  mask::maskType)

  p = zeros(dim.product)
  p[indx.anyActive] = p0
  demandD1only!(sim,p,dim,mat,par,indx)
  markup1 = p - mc[2]

  Dm1 = (sim.∂D∂p.*mat.host)'*(markup1.*indx.anyActive)
  idx = findall(indx.anyActive)
  for i=1:length(idx)
    j = idx[i]
    f[i] = sim.D[j] + Dm1[j]
  end

  # idx = findall(indx.anyActive)
  # for i=1:length(idx)
  #   j = idx[i]
  #   Dm1 = 0.0
  #   for k=1:dim.product
  #     markup1 = p[k]-mc[k]
  #     Dm1 += sim.∂D∂p[k,j] * mat.host[k,j] * markup1
  #   end
  #   f[i] = sim.D[j] + Dm1
  # end

end

function foc_retailPartial(
  p0  ::Vector{Float64},
  pj  ::Vector{Float64},
  wp  ::Vector{Float64},
  mc  ::Vector{Vector{Float64}},
  dim ::dimType,
  mat ::matType,
  par ::parType,
  sim ::simType,
  indx::indxType,
  mask::maskType)

  p = copy(pj)
  p[indx.anyActive] = p0

  demandD1only!(sim,p,dim,mat,par,indx)

  idx = findall(indx.anyActive)
  f = zeros(length(idx))
  @inbounds for i=1:length(idx)
    j = idx[i]
    Dm1 = 0.0
    Dm2 = 0.0
    for k in idx
      markup1 = (mask.mno[k]) ? p[k]-mc[2][k] : wp[k]-mc[1][k]
      markup2 = p[k] - wp[k] - mc[2][k]
      Dm1 += sim.∂D∂p[k,j] * mat.host[k,j] * markup1
      Dm2 += sim.∂D∂p[k,j] * mat.host[k,j] * markup2
    end
    f[i] = sim.D[j] + Dm1*mask.mno[j] + Dm2*mask.mvno[j]
  end
  return f
end

end

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand


text:
using Documenter, GalacticOptim

makedocs(
    sitename="GalacticOptim.jl",
    authors="Chris Rackauckas, Vaibhav Kumar Dixit et al.",
    clean=true,
    doctest=false,
    modules=[GalacticOptim],

    format=Documenter.HTML(assets=["assets/favicon.ico"],
                           canonical="https://galacticoptim.sciml.ai/stable/"),

    pages=[
        "GalacticOptim.jl: Unified Global Optimization Package" => "index.md",
        "Usage" => "usage.md",
        "Available optimizers" => "optimizers.md",
        "Examples" => "examples.md"
    ]
)

deploydocs(
    repo="github.com/SciML/GalacticOptim.jl";
    push_preview=true
)


text:
#=
This file contains ParamCtx, Param and KnetModule abstractions.
=#

"""Data structure used to store parameters. Currently, it is a 1d Any array"""
const ParamCtx = Array{Any, 1}

# The context abstraction

let
    global active_ctx, switch_ctx!, default_ctx, reset_ctx!

    #=default = ParamCtx()
    active = default=#
    active = ParamCtx()

    """`active_ctx()::ParamCtx` returns the context in use"""
    active_ctx()::ParamCtx = active

    """`switch_ctx!(new)::ParamCtx` switches to a new context and returns it"""
    switch_ctx!(new::ParamCtx)::ParamCtx = (active=new)

    #="""`default_ctx()::ParamCtx` returns the default context provided by the system"""
    default_ctx()::ParamCtx = default=#

    """`reset_ctx!()::ParamCtx` removes everything from the context"""
    reset_ctx!()::ParamCtx = (active=ParamCtx())
end

"""
Shorthand form of `active_ctx()`
"""
actx() = active_ctx()

"""`Param` stores an address to the parameters value.
Trainable arrays should be stored as `Param` inside modules.
User is responsible to keep track of `ParamCtx` object the `Param`
belongs to.

# Constructors
    `Param(w; ctx=active_ctx())` registers a parameter to the context ctx.
"""
type Param
    index::Integer
    take_grad::Bool
    Param(w; ctx::ParamCtx=active_ctx(), take_grad=true) =
        new(length(push!(ctx, w)), take_grad)
end

"""
`val(ctx, p::Union{Param, Void})` returns the value of parameter in
ctx.
This function should be used used to access parameters for AutoGrad to work.
"""
val(ctx, p::Union{Param, Void}) =
    p.take_grad ? ctx[p.index] : getval(ctx[p.index])

"""
`aval(p::Param)` returns value of p in the current `active_ctx()`.
It will stop differentiation when used inside forward funtions.
"""
aval(p::Param) = val(active_ctx(), p)

"""
`setval!(p::Param, w; ctx=active_ctx())` sets the value of `p` to `w` in `ctx`.
"""
setval!(p::Param, w; ctx=active_ctx()) = (ctx[p.index] = w)

"""`KnetModule` is the abstract type inherited by other modules."""
abstract type KnetModule end

import AutoGrad.grad

"""
`grad(m::KnetModule, loss)` returns a grad function of the form `(args..., y)->g`.
Recorded parameters are lived in `active_ctx`.

 # Arguments

    `m`: A KnetModule instance

    `loss`: A function of the form `(ypred, ygold)->scalar_loss`.
"""
function grad(m::KnetModule, loss::Function)
    _predict(w, x...; o...) = m(w, x...; o...)
    _loss(w, args...; kwargs...) = loss(_predict(w, args[1:end-1]...; kwargs...),
                                        args[end])
    lossgrad = grad(_loss)
    return (args...; kwargs...)->lossgrad(active_ctx(), args...; kwargs...)
end

import AutoGrad.gradloss
"""
`gradloss(m::Module, loss::Function)`
Same as `grad`, with only difference of returning the loss
"""
function gradloss(m::KnetModule, loss::Function)
    _predict(w, x...; o...) = m(w, x...; o...)
    _loss(w, args...; kwargs...) = loss(_predict(w, args[1:end-1]...; kwargs...),
                                        args[end])
    lossgrad = gradloss(_loss)
    return (args...; kwargs...)->lossgrad(active_ctx(), args...; kwargs...)
end

"""
`getgrad(p::Param, grads)` returns the corresponding gradient of `p` in
the grad array returned by any gradfn (generated by `grad` functions of Knet
or KnetModules).
"""
getgrad(p::Param, grads) = grads[p.index]

import Knet.optimizers

"""`optimizers(m::KnetModule, otype;o...)`
Creates a group of optimizers that correspond to parameters in module m.
otype specifies a Knet optimizer and `o...` is its options.
"""
optimizers(m::KnetModule, otype; o...) =
    Array{Any, 1}(map(_->otype(;o...), params(m)))

import Knet.update!
"""
`update!(m::KnetModule, grads, optims)`

Applies optimizers to parameters of `m`, using gradients stored in `grads`.
`optims` should have returned by `KnetModules.optimizers` for consistent behaviour.
"""
function update!(m::KnetModule, grads, optims)
    for (p, o) in zip(params(m), optims)
        update!(aval(p), getgrad(p, grads), o)
    end
end

function _populate_recursive(m, list, match_type)
    if isa(m, match_type)
        push!(list, m)
    end
    if isa(m, KnetModule)
        for fn in fieldnames(m)
            _populate_recursive(getfield(m, fn), list, match_type)
        end
    elseif isa(m, Array) || isa(m, Tuple)
        etype = eltype(m)
        # TODO: expand this filter
        if ~(etype <: Number || etype<:AbstractString)
            for l in m
                _populate_recursive(l, list, match_type)
            end
        end
    elseif isa(m, Associative)
        for k in keys(m)
            _populate_recursive(m[k], list, match_type)
        end
    else
        return
    end
end

"""
`params(m::KnetModule; kwargs)` returns parameters of the KnetModule m.
These include any value of type `Param` stored in fields, sub-modules,
dictionaries, arrays and tuples. Parameters are sorted based on their
locations in their respective context, so they are always returned in
the same order.
"""
function params(m::KnetModule)
    res = []
    _populate_recursive(m, res, Param)
    # For determinism (not solely for duplicate removal)
    sort!(res; lt=(r1, r2)->r1.index < r2.index)
    # Remove duplicates from the list
    #return res
    isempty(res) && return res
    res2 = Any[res[1]]
    for i = 2:length(res)
        (res[i-1].index != res[i].index) && push!(res2, res[i])
    end
    return res2
end

"""
`submodules(m::KnetModule)` returns list of modules contained in s.
These include any value of type `Param` stored in fields, sub-modules,
dictionaries, arrays and tuples.
"""
function modules(m::KnetModule)
    res = Set{KnetModule}()
    _populate_recursive(m, res, KnetModule)
    return Any[r for r in res]
end

"""
`convert_buffers!(m::KnetModule, atype)`
Creates a method of this function to convert non-parameter
state to `atype`.
You may create a method of this function to change cpu! and
gpu! behaviour of a particular module type.
"""
function convert_buffers!(m::KnetModule, atype)
    nothing
end

"""
`convert_params!(m::KnetModule, atype)`
Creates a method of this function to convert parameters
to `atype`.
You may create a method of this function to change `cpu!` and
`gpu!` behaviour of a particular module type.
"""
function convert_params!(m::KnetModule, atype)
    for p in params(m)
        setval!(p, atype(aval(p)))
    end
end

"""
`gpu!(m::KnetModule)` transfers all parameters in `m` to supported
gpu devices, which is identical to converting their valus to `KnetArray`.
If no gpu is avaliable, a warning is raised and transfer is ignored without
any error.
"""
function gpu!(m::KnetModule)
    if gpu() < 0
        warn("Gpu transfer is ignored, no available gpu")
        return m
    end
    for sm in modules(m)
        convert_params!(sm, ka)
        convert_buffers!(sm, ka)
    end
    return m
end

"""`cpu!(m::KnetModule)` transfers all parameters in `m` cpu,
which is identical to converting their values to `Array`"""
function cpu!(m::KnetModule)
    for sm in modules(m)
        convert_params!(sm, Array)
        convert_buffers!(sm, Array)
    end
    return m
end

"""
`training!(m::KnetModule)` makes the fields called `train` true if exists,
which may effect the execution mode in some modules. It also make `take_grad`
field of parameters true, which make them included in the training process.
"""
function training!(m::KnetModule)
    for sm in modules(m)
        if :train in fieldnames(sm)
            sm.train = true
        end
        for p in params(sm)
            p.take_grad = true
        end
    end
    return m
end

"""
`testing!(m::KnetModule)` makes the fields called `train` `false` if exists,
and freezes the parameters so that gradient is not taken w.r.t them. This
is done by unrecording the parameters during training, so that particular modules`
weights can be treated as constants by the AutoGrad.
"""
function testing!(m::KnetModule)
    for sm in modules(m)
        if :train in fieldnames(sm)
            sm.train = false
        end
        for p in params(sm)
            p.take_grad = false
        end
    end
    return m
end

# Macros for simple module operations
# TODO: support kwargs
#="""
`@mc expr` Converts the expression `m(a...)` to `forward(ctx, m, a...)`.
`ctx` should come from the upper scope.

# Example
    l = Linear(...)
    @mc l(x)
"""
macro mc(expr)
    if typeof(expr) == Expr && expr.head == :call
        return esc(:(forward(ctx, $(expr.args[1]),
                             $(expr.args[2:end]...))))
    end
    return expr
end

"""
`@run expr` Converts the expression m(a...) to forward(active_ctx(), m, a...).

# Example
    l = Linear(...)
    @run l(x)
"""
macro run(expr)
    if typeof(expr) == Expr && expr.head == :call
        return esc(:(forward(active_ctx(),
                             $(expr.args[1]),
                             $(expr.args[2:end]...))))
    end
    return expr
end

# Legacy
forward(ctx, m::KnetModule, args...) =
    m(ctx, args...)=#
    #=error("This forward signature is not implemented",
          "for abstract types and/or type ", typeof(m))=#

#=macro callable(Module, fn=forward)
    esc(:((m::$(Module))(ctx, args...; kwargs...) =
          $(fn)(ctx, m, args...; kwargs...)))
end=#

"""
`clean_ctx!(m::KnetModule)`
Switches to a new ctx where only parameters of `m` live.
"""
function clean_ctx!(m::KnetModule)
    is_ctx_clean(m) && return m
    new = ParamCtx()
    for p in params(m)
        p.index = length(push!(new, aval(p)))
    end
    switch_ctx!(new)
    return m
end

function is_ctx_clean(m::KnetModule)
    ps = Set{Int}(map(x->x.index, params(m)))
    as = Set{Int}(1:length(active_ctx()))
    return length(setdiff(as, ps)) == 0
end

function ctx_dict(m::KnetModule)
    cd = Dict()
    for p in params(m)
        cd[p.index] = aval(p)
    end
    return cd
end

# This thing is used with a loaded (not initialized) Knet
# module.
function from_ctx_dict!(m::KnetModule, cd::Dict; ctx=nothing)::ParamCtx
    if ctx == nothing; ctx = ParamCtx(); end
    for p in params(m)
        # FIXME: try not to break the abstraction
        p.index = length(push!(ctx, cd[p.index]))
    end
    return ctx
end

if Pkg.installed("JLD") !== nothing
    import JLD
else
    JLD = nothing
end

@inline assert_jld(fnname="this") =
    @assert (JLD !== nothing) "Install JLD for using $fname function"

function save_module(filename::String, m::KnetModule,
                     other_state::Dict=Dict())
    assert_jld("save_module")
    args = ["model", m, "ctx_dict", ctx_dict(m)]
    for k in keys(other_state)
        push!(args, k, other_state[k])
    end
    JLD.save(filename, args...)
end

function load_module(filename::String; ctx_switch=false)
    assert_jld("load_module")
    loaded = JLD.load(filename)
    model = loaded["model"]
    cdict = loaded["ctx_dict"]
    if ctx_switch
        switch_ctx!(from_ctx_dict!(model, cdict))
    else
        from_ctx_dict!(model, cdict; ctx=active_ctx())
    end
    if length(loaded) > 2
        other_state = loaded
        delete!(other_state, "model")
        delete!(other_state, "ctx_dict")
        return model, other_state
    end
    return model
end


text:
using Base: Int64
using BioCCP
@testset "BioCCP" begin

    n = 20

    # Equal probabilities
    p_uniform = ones(n)/n

    # Probabilities following Zipf's law
    ρ = 10
	α = exp(log(ρ)/(n-1))
	p_zipf = collect(α.^-(1:n))
    p_zipf = p_zipf ./ sum(p_zipf)

    @testset "Expectation" begin
        @test expectation_minsamplesize(n; p = p_uniform) isa Int64
        @test expectation_minsamplesize(n; p = p_uniform) == Int(ceil(n*sum(1 ./ (1:n))))
        @test expectation_minsamplesize(n; p = p_uniform) < expectation_minsamplesize(n; p = p_zipf)
        @test expectation_minsamplesize(n; p = p_uniform, m = 1) < expectation_minsamplesize(n; p = p_uniform, m = 2)
        @test Int(ceil(expectation_minsamplesize(n; p = p_uniform, m = 1, r = 1)/3)) == expectation_minsamplesize(n; p = p_uniform, m = 1, r = 3)
    end

    @testset "Standard deviation" begin
        @test std_minsamplesize(n; p = p_uniform) isa Int64
        @test std_minsamplesize(n; p = p_uniform) < std_minsamplesize(n; p = p_zipf)
        @test std_minsamplesize(n; p = p_uniform, m = 1) < std_minsamplesize(n; p = p_uniform, m = 2)
    end

    t = 500

    @testset "Success probability" begin
        @test success_probability(n, t; p = p_uniform) isa Float64
        @test success_probability(n, t; p = p_uniform) > success_probability(n, t; p = p_zipf)
        @test success_probability(n, t; p = p_uniform, m = 1) > success_probability(n, t; p = p_uniform, m = 2)
        @test all(success_probability(n, t; p = p_uniform) .< success_probability.(n, [t+10, t+100, t+1000]; p = p_uniform))
        @test all(success_probability.(n, [1, 5, 10, 100, 500]) .<= 1)
        @test all(success_probability.(n, [1, 5, 10, 100, 500]) .>= 0)

    end

    @testset "Expected fraction collected" begin
        @test expectation_fraction_collected(n, t; p = p_uniform) isa Float64
        @test expectation_fraction_collected(n, t; p = p_uniform) > expectation_fraction_collected(n, t; p = p_zipf)
        @test all(expectation_fraction_collected(n, t; p = p_uniform) .< expectation_fraction_collected.(n, [t+10, t+100, t+1000]; p = p_uniform))
        @test all(expectation_fraction_collected.(n, [1, 10, 100, 1000, 10000]) .<= 1)
        @test all(expectation_fraction_collected.(n, [1, 10, 100, 1000, 10000]) .>= 0)

    end

    k = 2
    pᵢ = 0.005
    r = 3

    @testset "Probability occurence" begin
        @test prob_occurrence_module(pᵢ, t, r, k) isa Float64
        @test all(prob_occurrence_module.(pᵢ, [0, 10, 100, 1000, 10000], r, k) .>= 0)
        @test all(prob_occurrence_module.(pᵢ, [0, 10, 100, 1000, 10000], r, k) .<= 1)
    end
end

max_issues_repo_name: kirstvh/BioCCP

max_forks_repo_name: kirstvh/BioCCP


text:
# Copyright 2020 Google LLC

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

#     https://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
raw_input = readlines("day05input.txt")

function find_seat_id(seat)
    lo = 0
    hi = 127
    mid = 0
    for half ∈ seat[1:7]
        mid = (lo + hi) ÷ 2
        if half == 'F'
            hi = mid
        elseif half == 'B'
            lo = mid
        end
    end
    if seat[7] == 'B'
        mid += 1
    end
    row = mid

    lo = 0
    hi = 7
    mid = 0
    for side ∈ seat[8:10]
        mid = (lo + hi) ÷ 2
        if side == 'L'
            hi = mid
        elseif side == 'R'
            lo = mid
        end
    end
    if seat[10] == 'R'
        mid += 1
    end
    col = mid

    return row * 8 + col
end

@assert find_seat_id("FBFBBFFRLR") == 357
@assert find_seat_id("BFFFBBFRRR") == 567
@assert find_seat_id("FFFBBBFRRR") == 119
@assert find_seat_id("BBFFBBFRLL") == 820

part1(inputs) = maximum(map(find_seat_id, inputs))

function part2(inputs)
    ids = sort(map(find_seat_id, inputs))
    for i in 1:(length(ids)-1)
        if ids[i] + 1 != ids[i+1]
            return ids[i] + 1
        end
    end
end

part1(raw_input)
part2(raw_input)


text:
export amplitudeDAC, frequencyDAC, phaseDAC, modeDAC,
       DCSignDAC, signalTypeDAC, offsetDAC, jumpSharpnessDAC, passPDMToFastDAC, DACConfig, configureFastDAC,
       waveforms, DACPerformanceData, computeRamping, rampUp, rampUpSteps, rampUpTotalSteps, sequenceRepetitions,
       prepareSlowDAC, slowDACStepsPerFrame, slowDACStepsPerSequence, samplesPerSlowDACStep,
       enableDACLUT, setArbitraryLUT, setConstantLUT, setPauseLUT, setRangeLUT, numSlowDACChan,
       appendSequence, popSequence, clearSequence, prepareSequence, resetAfterSequence, numLostStepsSlowADC,
       AbstractSequence, ArbitrarySequence, ConstantSequence, PauseSequence, RangeSequence, fastDACConfig

struct DACPerformanceData
  uDeltaControl::UInt8
  uDeltaSet::UInt8
  minDeltaControl::UInt8
  maxDeltaSet::UInt8
end

struct DACConfig
  amplitudes::Array{Union{Float64, Nothing}}
  offsets::Array{Union{Float64, Nothing}}
  frequencies::Array{Union{Float64, Nothing}}
  phases::Array{Union{Float64, Nothing}}
  signalTypes::Array{Union{String, Nothing}}
  jumpSharpness::Array{Union{Float64, Nothing}}
  function DACConfig()
    new(Array{Union{Float64, Nothing}}(nothing, 2,4), Array{Union{Float64, Nothing}}(nothing, 2),
    Array{Union{Float64, Nothing}}(nothing, 2,4), Array{Union{Float64, Nothing}}(nothing, 2,4),
    Array{Union{String, Nothing}}(nothing, 2), Array{Union{Float64, Nothing}}(nothing, 2))
  end
end

function getDACScpiPrefix(forSequence)
  prefix = "RP:DAC"
  if forSequence
    prefix = prefix * ":SEQ"
  end
  return prefix
end

function passPDMToFastDAC(rp::RedPitaya, val::Bool)
  valStr = val ? "ON" : "OFF"
  send(rp, string("RP:PassPDMToFastDAC ", valStr))
end
passPDMToFastDAC(rp::RedPitaya) = occursin("ON", query(rp,"RP:PassPDMToFastDAC?"))

function amplitudeDAC(rp::RedPitaya, channel, component)
  command = string("RP:DAC:CH", Int(channel)-1, ":COMP", Int(component)-1, ":AMP?")
  return query(rp, command, Float64)
end
function amplitudeDAC(rp::RedPitaya, channel, component, value; forSequence=false)
  if value > 1.0
    error("$value is larger than 1.0 V!")
  end
  command = string(getDACScpiPrefix(forSequence), ":CH", Int(channel)-1, ":COMP",
                   Int(component)-1, ":AMP ", Float64(value))
  return send(rp, command)
end
function amplitudeDAC(config::DACConfig, channel, component, value)
  if value > 1.0
    error("$value is larger than 1.0 V!")
  end
  config.amplitudes[channel, component] = value
end

function offsetDAC(rp::RedPitaya, channel)
  command = string("RP:DAC:CH", Int(channel)-1, ":OFF?")
  return query(rp, command, Float64)
end
function offsetDAC(rp::RedPitaya, channel, value; forSequence=false)
  if value > 1.0
    error("$value is larger than 1.0 V!")
  end
  command = string(getDACScpiPrefix(forSequence), ":CH", Int(channel)-1, ":OFF ", Float64(value))
  return send(rp, command)
end
function offsetDAC(config::DACConfig, channel, value)
  if value > 1.0
    error("$value is larger than 1.0 V!")
  end
  config.offsets[channel] = value
end

function frequencyDAC(rp::RedPitaya, channel, component)
  command = string("RP:DAC:CH", Int(channel)-1, ":COMP", Int(component)-1, ":FREQ?")
  return query(rp, command, Float64)
end
function frequencyDAC(rp::RedPitaya, channel, component, value; forSequence=false)
  command = string(getDACScpiPrefix(forSequence), ":CH", Int(channel)-1, ":COMP",
                   Int(component)-1, ":FREQ ", Float64(value))
  send(rp, command)
end
function frequencyDAC(config::DACConfig, channel, component, value)
  config.frequencies[channel, component] = value
end

function phaseDAC(rp::RedPitaya, channel, component)
  command = string("RP:DAC:CH", Int(channel)-1, ":COMP", Int(component)-1, ":PHA?")
  return query(rp, command, Float64)
end
function phaseDAC(rp::RedPitaya, channel, component, value; forSequence=false)
  command = string(getDACScpiPrefix(forSequence), ":CH", Int(channel)-1, ":COMP",
                   Int(component)-1, ":PHA ", Float64(value))
  send(rp, command)
end
function phaseDAC(config::DACConfig, channel, component, value)
  config.phases[channel, component] = value
end

function jumpSharpnessDAC(rp::RedPitaya, channel)
  command = string("RP:DAC:CH", Int(channel)-1, ":JumpSharpness?")
  return query(rp, command, Float64)
end
function jumpSharpnessDAC(rp::RedPitaya, channel, value; forSequence=false)
  command = string(getDACScpiPrefix(forSequence), ":CH", Int(channel)-1, ":JumpSharpness ", Float64(value))
  send(rp, command)
end
function jumpSharpnessDAC(config::DACConfig, channel, value)
  config.jumpSharpness[channel] = value
end

#"STANDARD" or "AWG" (not yet supported)
function modeDAC(rp::RedPitaya)
  modeDAC_(rp, 1)
  modeDAC_(rp, 2)
end
function modeDAC(rp::RedPitaya, mode::String)
  modeDAC_(rp, mode, 1)
  modeDAC_(rp, mode, 2)
end

function modeDAC_(rp::RedPitaya, channel=1)
  #return query(rp, "RP:DAC:CH", Int(channel)-1,":MODe?")[2:end-1]
  return query(rp, "RP:DAC:MODe?")[2:end-1]
end
function modeDAC_(rp::RedPitaya, mode::String, channel=1)
  #send(rp, string("RP:DAC:CH", Int(channel)-1, ":MODe ", mode))
  send(rp, string("RP:DAC:MODe ", mode))
end

function signalTypeDAC(rp::RedPitaya, channel)
  command = string("RP:DAC:CH", Int(channel)-1, ":SIGnaltype?")
  @show command
  return query(rp, command)
end

waveforms() = ["SINE","SQUARE","TRIANGLE","SAWTOOTH"]

function signalTypeDAC(rp::RedPitaya, channel, sigType::String; forSequence=false)
  if !(sigType in waveforms() )
    error("Signal type $sigType not supported!")
  end

  command = string(getDACScpiPrefix(forSequence), ":CH", Int(channel)-1, ":SIGnaltype ", sigType)
  return send(rp, command)
end
function signalTypeDAC(config::DACConfig, channel, sigType::String)
  if !(sigType in waveforms() )
    error("Signal type $sigType not supported!")
  end
  config.signalTypes[channel] = sigType
end

function DCSignDAC(rp::RedPitaya, channel)
  command = string("RP:DAC:CH", Int(channel)-1, ":SIGn?")
  return query(rp, command, Int64) == "POSITIVE" ? 1 : -1
end

function DCSignDAC(rp::RedPitaya, channel, sign::Integer)
  command = string("RP:DAC:CH", Int(channel)-1, ":SIGn ",
             sign > 0 ? "POSITIVE" : "NEGATIVE")
  return send(rp, command)
end

function configureFastDAC(rp::RedPitaya, config::DACConfig; forSequence = false)
  for ch = 1:2

    for cmp = 1:4
      amplitude = config.amplitudes[ch, cmp]
      isnothing(amplitude) || amplitudeDAC(rp, ch, cmp, amplitude, forSequence = forSequence)
      frequency = config.frequencies[ch, cmp]
      isnothing(frequency) || frequencyDAC(rp, ch, cmp, frequency, forSequence = forSequence)
      phase = config.phases[ch, cmp]
      isnothing(phase) || phaseDAC(rp, ch, cmp, phase, forSequence = forSequence)
    end

    offset = config.offsets[ch]
    isnothing(offset) || offsetDAC(rp, ch, offset, forSequence = forSequence)
    signalType = config.signalTypes[ch]
    isnothing(signalType) || signalTypeDAC(rp, ch, signalType, forSequence = forSequence)

  end
end

function readDACPerformanceData(rp::RedPitaya)
  perf = read!(rp.dataSocket, Array{UInt8}(undef, 4))
  return DACPerformanceData(perf[1], perf[2], perf[3], perf[4])
end

numLostStepsSlowADC(rp::RedPitaya) = query(rp,"RP:DAC:SEQ:LostSteps?", Int64) # TODO slowADC vs slowDAC in name and string

numSlowDACChan(rp::RedPitaya) = query(rp,"RP:DAC:SEQ:CHan?", Int64)
function numSlowDACChan(rp::RedPitaya, value)
  if value <= 0 || value > 4
    error("Num slow DAC channels needs to be between 1 and 4!")
  end
  send(rp, string("RP:DAC:SEQ:CHan ", Int64(value)))
end

function setValueLUT(rp::RedPitaya, lut::Union{Array, Nothing}, type::String="ARBITRARY")
  send(rp, string("RP:DAC:SEQ:LUT:", type))
  @debug "Writing arbitrary LUT"
  if !isnothing(lut)
    lutFloat32 = map(Float32, lut)
    write(rp.dataSocket, lutFloat32)
  end
end

function enableDACLUT(rp::RedPitaya, lut::Array)
  lutBool = map(Bool, lut)
  send(rp, string("RP:DAC:SEQ:LUT:ENaBle"))
  @debug "Writing enable DAC LUT"
  write(rp.dataSocket, lutBool)
end

samplesPerSlowDACStep(rp::RedPitaya) = query(rp,"RP:DAC:SEQ:SAMPlesPerStep?", Int64)
function samplesPerSlowDACStep(rp::RedPitaya, value)
  send(rp, string("RP:DAC:SEQ:SAMPlesPerStep ", value))
end

slowDACStepsPerSequence(rp::RedPitaya) = query(rp,"RP:DAC:SEQ:STEPsPerSequence?", Int64)
function slowDACStepsPerSequence(rp::RedPitaya, value)
  send(rp, string("RP:DAC:SEQ:STEPsPerSequence ", value))
end

function prepareSlowDAC(rp::RedPitaya, samplesPerStep, stepsPerSequence, numOfChan)
  numSlowDACChan(rp, numOfChan)
  samplesPerSlowDACStep(rp, samplesPerStep)
  slowDACStepsPerSequence(rp, stepsPerSequence)
end

function slowDACStepsPerFrame(rp::RedPitaya, stepsPerFrame)
  samplesPerFrame = rp.periodsPerFrame * rp.samplesPerPeriod
  samplesPerStep = div(samplesPerFrame, stepsPerFrame)
  samplesPerSlowDACStep(rp, samplesPerStep) # Sets PDMClockDivider
  #slowDACStepsPerSequence(rp, stepsPerFrame)
end

ramping(rp::RedPitaya, rampSteps::Int32, rampTotalSteps::Int32) = send(rp, string("RP:DAC:SEQ:RaMPing ", rampSteps, ",", rampTotalSteps))
rampUp(rp::RedPitaya, rampUpSteps::Int32, rampUpTotalSteps::Int32) = send(rp, string("RP:DAC:SEQ:RaMPing:UP ", rampUpSteps, ",", rampUpTotalSteps))
rampDown(rp::RedPitaya, rampDownSteps::Int32, rampDownTotalSteps::Int32) = send(rp, string("RP:DAC:SEQ:RaMPing:DOWn ", rampDownSteps, ",", rampDownTotalSteps))

rampingSteps(rp::RedPitaya) = query(rp, "RP:DAC:SEQ:RaMPing:STEPs?", Int32)
rampingSteps(rp::RedPitaya, value::Int32) = send(rp, string("RP:DAC:SEQ:RaMPing:STEPs ", value))
rampingTotalSteps(rp::RedPitaya) = query(rp, "RP:DAC:SEQ:RaMPing:TOTAL?", Int32)
rampingTotalSteps(rp::RedPitaya, value::Int32) = send(rp, string("RP:DAC:SEQ:RaMPing:TOTAL ", value))

rampUpSteps(rp::RedPitaya) = query(rp, "RP:DAC:SEQ:RaMPing:UP:STEPs?", Int32)
rampUpSteps(rp::RedPitaya, value::Int32) = send(rp, string("RP:DAC:SEQ:RaMPing:UP:STEPs ", value))
rampUpTotalSteps(rp::RedPitaya) = query(rp, "RP:DAC:SEQ:RaMPing:UP:TOTAL?", Int32)
rampUpTotalSteps(rp::RedPitaya, value::Int32) = send(rp, string("RP:DAC:SEQ:RaMPing:UP:TOTAL ", value))

rampDownSteps(rp::RedPitaya) = query(rp, "RP:DAC:SEQ:RaMPing:DOWN:STEPs?", Int32)
rampDownSteps(rp::RedPitaya, value::Int32) = send(rp, string("RP:DAC:SEQ:RaMPing:DOWN:STEPs ", value))
rampDownTotalSteps(rp::RedPitaya) = query(rp, "RP:DAC:SEQ:RaMPing:DOWN:TOTAL?", Int32)
rampDownTotalSteps(rp::RedPitaya, value::Int32) = send(rp, string("RP:DAC:SEQ:RaMPing:DOWN:TOTAL ", value))

function resetAfterSequence(rp::RedPitaya, val::Bool)
  valStr = val ? "ON" : "OFF"
  send(rp, string("RP:DAC:SEQ:RESETafter ", valStr))
end
resetAfterSequence(rp::RedPitaya) = occursin("ON", query(rp,"RP:DAC:SEQ:RESETafter?"))

sequenceRepetitions(rp::RedPitaya) = query(rp, "RP:DAC:SEQ:REPetitions?", Int32)
function sequenceRepetitions(rp::RedPitaya, value::Int)
  send(rp, string("RP:DAC:SEQ:REPetitions ", value))
end

appendSequence(rp::RedPitaya) = send(rp, "RP:DAC:SEQ:APPend")
popSequence(rp::RedPitaya) = send(rp, "RP:DAC:SEQ:POP")
clearSequence(rp::RedPitaya) = send(rp, "RP:DAC:SEQ:CLEAR")
prepareSequence(rp::RedPitaya) = query(rp, "RP:DAC:SEQ:PREPare?", Bool)

# Helper function for sequences
abstract type AbstractSequence end

mutable struct ArbitrarySequence <: AbstractSequence
  lut::Array{Float32}
  enable::Union{Array{Bool}, Nothing}
  stepsPerRepetition::Int
  repetitions::Int
  rampUpSteps::Int32
  rampUpTotalSteps::Int32
  rampDownSteps::Int32
  rampDownTotalSteps::Int32
  fastDAC::DACConfig
  resetAfter::Bool
end
ArbitrarySequence(lut, enable, stepsPerRepetition, repetitions, upSteps, upTotalSteps, downSteps, rampDownTotalSteps, reset=false) = ArbitrarySequence(lut, enable, stepsPerRepetition, repetitions, upSteps, upTotalSteps, downSteps, rampDownTotalSteps, DACConfig(), reset)
ArbitrarySequence(lut, enable, stepsPerRepetition, repetitions, steps, totalSteps, reset=false) = ArbitrarySequence(lut, enable, stepsPerRepetition, repetitions, steps, totalSteps, steps, totalSteps, DACConfig(), reset)
ArbitrarySequence(lut, enable, stepsPerRepetition, repetitions, (steps, totalSteps)::Tuple, reset=false) = ArbitrarySequence(lut, enable, stepsPerRepetition, repetitions, steps, totalSteps, steps, totalSteps, DACConfig(), reset)
ArbitrarySequence(lut, enable, stepsPerRepetition, repetitions, (upSteps, upTotalSteps)::Tuple, (downSteps, downTotalSteps)::Tuple, reset=false) = ArbitrarySequence(lut, enable, stepsPerRepetition, repetitions, upSteps, upTotalSteps, downSteps, downTotalSteps, DACConfig(), reset)

stepsPerRepetition(seq::ArbitrarySequence) = seq.stepsPerRepetition
rampUpSteps(seq::ArbitrarySequence) = seq.rampUpSteps
rampUpTotalSteps(seq::ArbitrarySequence) = seq.rampUpTotalSteps
rampDownSteps(seq::ArbitrarySequence) = seq.rampDownSteps
rampDownTotalSteps(seq::ArbitrarySequence) = seq.rampDownTotalSteps
repetitions(seq::ArbitrarySequence) = seq.repetitions
enableLUT(seq::ArbitrarySequence) = seq.enable
fastDACConfig(seq::ArbitrarySequence) = seq.fastDAC
resetAfterSequence(seq::ArbitrarySequence) = seq.resetAfter

mutable struct ConstantSequence <: AbstractSequence
  lut::Array{Float32}
  enable::Union{Array{Bool}, Nothing}
  stepsPerRepetition::Int
  repetitions::Int
  rampUpSteps::Int32
  rampUpTotalSteps::Int32
  rampDownSteps::Int32
  rampDownTotalSteps::Int32
  fastDAC::DACConfig
  resetAfter::Bool
end
ConstantSequence(lut, enable, stepsPerRepetition, repetitions, upSteps, upTotalSteps, downSteps, downTotalSteps, reset=false) = ConstantSequence(lut, enable, stepsPerRepetition, repetitions, upSteps, upTotalSteps, downSteps, downTotalSteps, DACConfig(), reset)
ConstantSequence(lut, enable, stepsPerRepetition, repetitions, steps, totalSteps, reset=false) = ConstantSequence(lut, enable, stepsPerRepetition, repetitions, steps, totalSteps, steps, totalSteps, DACConfig(), reset)
ConstantSequence(lut, enable, stepsPerRepetition, repetitions, (steps, totalSteps)::Tuple, reset=false) = ConstantSequence(lut, enable, stepsPerRepetition, repetitions, steps, totalSteps, steps, totalSteps, DACConfig(), reset)
ConstantSequence(lut, enable, stepsPerRepetition, repetitions, (upSteps, upTotalSteps)::Tuple, (downSteps, downTotalSteps)::Tuple, reset=false) = ConstantSequence(lut, enable, stepsPerRepetition, repetitions, upSteps, upTotalSteps, downSteps, downTotalSteps, DACConfig(), reset)

stepsPerRepetition(seq::ConstantSequence) = seq.stepsPerRepetition
rampUpSteps(seq::ConstantSequence) = seq.rampUpSteps
rampUpTotalSteps(seq::ConstantSequence) = seq.rampUpTotalSteps
rampDownSteps(seq::ConstantSequence) = seq.rampDownSteps
rampDownTotalSteps(seq::ConstantSequence) = seq.rampDownTotalSteps
repetitions(seq::ConstantSequence) = seq.repetitions
enableLUT(seq::ConstantSequence) = seq.enable
fastDACConfig(seq::ConstantSequence) = seq.fastDAC
resetAfterSequence(seq::ConstantSequence) = seq.resetAfter

mutable struct PauseSequence <: AbstractSequence
  enable::Union{Array{Bool}, Nothing}
  stepsPerRepetition::Int
  repetitions::Int
  fastDAC::DACConfig
  resetAfter::Bool
end
PauseSequence(enable, stepsPerRepetition, repetitions, reset=false) = PauseSequence(enable, stepsPerRepetition, repetitions, DACConfig(), reset)

stepsPerRepetition(seq::PauseSequence) = seq.stepsPerRepetition
rampUpSteps(seq::PauseSequence) = Int32(0)
rampUpTotalSteps(seq::PauseSequence) = Int32(0)
rampDownSteps(seq::PauseSequence) = Int32(0)
rampDownTotalSteps(seq::PauseSequence) = Int32(0)
repetitions(seq::PauseSequence) = seq.repetitions
enableLUT(seq::PauseSequence) = seq.enable
fastDACConfig(seq::PauseSequence) = seq.fastDAC
resetAfterSequence(seq::PauseSequence) = seq.resetAfter

mutable struct RangeSequence <: AbstractSequence
  lut::Array{Float32}
  enable::Union{Array{Bool}, Nothing}
  stepsPerRepetition::Int
  repetitions::Int
  rampUpSteps::Int32
  rampUpTotalSteps::Int32
  rampDownSteps::Int32
  rampDownTotalSteps::Int32
  fastDAC::DACConfig
  resetAfter::Bool
end
RangeSequence(lut, enable, stepsPerRepetition, repetitions, upSteps, upTotalSteps, downSteps, downTotalSteps, reset=false) = RangeSequence(lut, enable, stepsPerRepetition, repetitions, upSteps, upTotalSteps, downSteps, downTotalSteps, DACConfig(), reset)
RangeSequence(lut, enable, stepsPerRepetition, repetitions, steps, totalSteps, reset=false) = RangeSequence(lut, enable, stepsPerRepetition, repetitions, steps, totalSteps, steps, totalSteps, DACConfig(), reset)
RangeSequence(lut, enable, stepsPerRepetition, repetitions, (steps, totalSteps)::Tuple, reset=false) = RangeSequence(lut, enable, stepsPerRepetition, repetitions, steps, totalSteps, steps, totalSteps, DACConfig(), reset)
RangeSequence(lut, enable, stepsPerRepetition, repetitions, (upSteps, upTotalSteps)::Tuple, (downSteps, downTotalSteps)::Tuple, reset=false) = RangeSequence(lut, enable, stepsPerRepetition, repetitions, upSteps, upTotalSteps, downSteps, downTotalSteps, DACConfig(), reset)

stepsPerRepetition(seq::RangeSequence) = seq.stepsPerRepetition
rampUpSteps(seq::RangeSequence) = seq.rampUpSteps
rampUpTotalSteps(seq::RangeSequence) = seq.rampUpTotalSteps
rampDownSteps(seq::RangeSequence) = seq.rampDownSteps
rampDownTotalSteps(seq::RangeSequence) = seq.rampDownTotalSteps
repetitions(seq::RangeSequence) = seq.repetitions
enableLUT(seq::RangeSequence) = seq.enable
fastDACConfig(seq::RangeSequence) = seq.fastDAC
resetAfterSequence(seq::RangeSequence) = seq.resetAfter

setLUT(rp::RedPitaya, seq::AbstractSequence) = error("Sequence did not implement setLUT")

function setLUT(rp::RedPitaya, seq::ArbitrarySequence)
  setValueLUT(rp, seq.lut, "ARBITRARY")
end

function setLUT(rp::RedPitaya, seq::ConstantSequence)
  setValueLUT(rp, seq.lut, "CONSTANT")
end

function setLUT(rp::RedPitaya, seq::PauseSequence)
  setValueLUT(rp, nothing, "PAUSE")
end

function setLUT(rp::RedPitaya, seq::RangeSequence)
  setValueLUT(rp, seq.lut, "RANGE")
end

function computeRamping(dec, samplesPerStep, stepsPerSeq, rampTime, rampFraction)
  bandwidth = 125e6/dec
  samplesPerRotation = samplesPerStep * stepsPerSeq
  totalRotations = Int32(ceil(rampTime/(samplesPerRotation/bandwidth)))
  totalSteps = totalRotations * stepsPerSeq
  steps = Int32(ceil(rampTime * rampFraction/(samplesPerStep/bandwidth)))
  #@show samplesPerRotation totalRotations totalSteps steps rampTime rampFraction
  return (steps, totalSteps)
end
computeRamping(rp::RedPitaya, stepsPerSeq ,rampTime, rampFraction) = computeRamping(decimation(rp), samplesPerSlowDACStep(rp), stepsPerSeq, rampTime, rampFraction)

function appendSequence(rp::RedPitaya, seq::AbstractSequence)
  slowDACStepsPerSequence(rp, stepsPerRepetition(seq))
  rampUp(rp, rampUpSteps(seq), rampUpTotalSteps(seq))
  rampDown(rp, rampDownSteps(seq), rampDownTotalSteps(seq))
  sequenceRepetitions(rp, repetitions(seq))
  setLUT(rp, seq)
  enable = enableLUT(seq)
  if !isnothing(enable)
    enableDACLUT(rp, enable)
  end
  configureFastDAC(rp, fastDACConfig(seq), forSequence = true)
  resetAfterSequence(rp, resetAfterSequence(seq))
  appendSequence(rp)
end

function length(seq::AbstractSequence)
  result = stepsPerRepetition(seq) * repetitions(seq) + rampUpTotalSteps(seq) + rampDownTotalSteps(seq)
  result = resetAfterSequence(seq) ? result + 1 : result
  return result
end

function start(seq::AbstractSequence)
  return rampUpTotalSteps(seq)
end

max_stars_repo_name: tknopp/RedPitayaDAQServer

max_issues_repo_name: tknopp/RedPitayaDAQServer

max_forks_repo_name: tknopp/RedPitayaDAQServer


text:
verbose = true
 using Dates
 if verbose println(now()) end
 if verbose && !isdefined(Main,:RvSpectML)  println("# Loading RvSpecML")    end
 using PDMats  # SOme dependancy needs this explicit first
 using RvSpectMLBase
 using EchelleCCFs
 using EchelleInstruments#, EchelleInstruments.NEID
 using RvSpectML
 if verbose println("# Loading NeidSolarScripts")    end
 using SunAsAStar
 using NeidSolarScripts
 println("# Loading other packages 1/2")
 using ArgParse
 using Markdown

println("# Parsing arguments...")
   #lsf_width = 3.0e3
   function parse_commandline()
     s = ArgParseSettings( description = "Calculate RVs from JLD2 file with NEID's daily order CCFs.")
     #import_settings!(s, s_files_only, args_only=false)
     add_arg_group!(s, "Files", :argg_files)
     @add_arg_table! s begin
         "input"
             help = "Filename for input with daily order CCFs (jld2)"
             arg_type = String
             default = "daily_ccfs.jld2"
         "output"
             help = "Filename for output dailiy RVs (csv)"
             arg_type = String
             default = "daily_rvs.csv"
         "--overwrite"
            help = "Specify it's ok to overwrite the output file."
            #default = true
            action = :store_true
      end
      #=
      add_arg_group!(s, "CCF parameters", :argg_ccf_param)
      @add_arg_table! s begin
        "--orders_to_use"
           help = "First and last order _index_ to compute CCFs for"
           nargs = 2
           arg_type = Int64
           #default = [ min_order(NEID2D()), max_order(NEID2D()) ]
           #default = [ first(orders_to_use_default(NEID2D())), last(orders_to_use_default(NEID2D())) ]
        "--mask_scale_factor"
            help = "Specify CCF mask width scale as multiple of NEID default v width " * string(default_ccf_mask_v_width(NEID2D()))
            arg_type = Float64
            #default = round(lsf_width/default_ccf_mask_v_width(NEID2D()), sigdigits=3)
        "--line_width_50_default"
            help = "Specify default line full width half maximum"
            arg_type = Float64
            #default = 7.9e3 # m/2
      end
      =#
      add_arg_group!(s, "Algorithm for measuring RVs ", :argg_alg_param)
      @add_arg_table! s begin
         "--rv_alg"
            help = "Algorithm for measuring RVs."
            arg_type = String
            default = "template"
            #default = "gaussian"
         "--frac_of_width_to_fit"
            help = "Regions of CCF to fit, as fraction of measured width."
            arg_type = Float64
            default = 0.5
         "--measure_width_at_frac_depth"
            help = "Measure CCF width at this fractional depth."
            arg_type = Float64
            default = 0.5
         "--template_file"
            help = "Filename with CCF template stored as 'order_ccfs' with 'clean_obs_mask' for which observations to include (jld2)."
            arg_type = String
      end
      add_arg_group!(s, "Filter manifest for ", :argg_filter_param)
      @add_arg_table! s begin
         "--target"
            help = "Target field"
            arg_type = String
            default = "Sun"
         "--datestr"
            help = "Filenames that contain date string."
            arg_type = String
        "--max_solar_hour_angle"
           help = "Maximum absolute value of solar hour angle."
           arg_type = Float64
           default = 3.12
         "--max_airmass"
            help = "Maximum airmass."
            arg_type = Float64
            default = 2.0
         "--min_expmeter"
            help = "Minimum exposure meter flux relative to model flux."
            arg_type = Float64
            default = 6e4
#=
         "--min_expmeter_factor"
            help = "Minimum exposure meter flux relative to model flux."
            arg_type = Float64
            default = 0.9
         "--max_expmeter_factor"
            help = "Maximum exposure meter flux relative to model flux."
            arg_type = Float64
            default = 1.1
=#
         "--max_expmeter_frac_rms"
            help = "Maximum fractional RMS of exposure meter flux."
            arg_type = Float64
            default = 0.01
#=
         "--min_pyrhelio_factor"
            help = "Minimum pyrheliometer flux relative to model flux."
            arg_type = Float64
            default = 0.9
         "--max_pyrhelio_factor"
            help = "Maximum pyrheliometer flux relative to model flux."
            arg_type = Float64
            default = 1.1
=#
         "--max_pyrhelio_frac_rms"
            help = "Maximum fractional RMS of pyrheliometer meter flux."
            arg_type = Float64
            default = 0.0035
         #=
         "--min_drp_snr"
            help = "Minimum extracted SNR reported by NEID DRP."
            arg_type = Float64
            #default = 0.0
         =#
            #=
         "--start_time"
            help = "Specify daily start time for CCFs to be processed (HH MM)"
            nargs = 2
            arg_type = Int64
            default = [0, 0]
            #default = [18, 30]
        "--stop_time"
           help = "Specify daily stop time for CCFs to be processed (HH MM)"
           nargs = 2
           arg_type = Int64
           #default = [ min_order(NEID2D()), max_order(NEID2D()) ]
           #default = [22, 30]
           default = [23, 59]
           =#
         "--max_num_spectra"
            help = "Maximum number of spectra to process."
            arg_type = Int
            default = 300   # TODO: Increase  after finish testing
      end
      return parse_args(s)
 end
 args = parse_commandline()
 @assert any(map(alg->occursin(alg,args["rv_alg"]), ["template","gaussian","quadratic"] ))

if verbose   println("# Loading other packages 2/2")    end
 using CSV, DataFrames, Query, Dates
 using JLD2, FileIO, MD5 #, SHA
 using StatsBase, Statistics, NaNMath

 # Filename arguments
 @assert isfile(args["input"]) || islink(args["input"])
 @assert match(r"\.jld2$",args["input"]) != nothing
 daily_ccf_filename = args["input"]

 args["overwrite"] = true   # TODO: Comment after done testing
 if isfile(args["output"]) && !args["overwrite"]
    error("# Output file " * args["output"] * " already exists (size = " * string(filesize(args["output"])) * " ).")
 end
 @assert !isfile(args["output"]) || args["overwrite"]
 @assert match(r"\.csv$",args["output"]) != nothing
 daily_rvs_filename = args["output"]
 touch(daily_rvs_filename)

#=
 if isfile(args["report"]) && !args["overwrite"]
    error("# Report file " * args["report"] * " already exists (size = " * string(filesize(args["report"])) * " ).")
 end
 @assert !isfile(args["report"]) || args["overwrite"] == true
 @assert match(r"\.md$",args["report"]) != nothing
 daily_report_filename = args["report"]
 touch(daily_report_filename)
=#

 file_hashes = Dict{String,String}()
 start_processing_time = now()

if filesize(daily_ccf_filename) > 0
   println("# Reading daily CCFs from ", daily_ccf_filename)
   input_data = load(daily_ccf_filename)
else
   println("# Empy daily CCF file.  Creating empty ", daily_rvs_filename)
   touch(daily_rvs_filename)
   exit(0)
end

println("# Filtering for usable observations...")
manifest = input_data["manifest"]

manifest_use = manifest |>
 @filter( args["target"] == nothing || _.target == args["target"] ) |>
 @filter( args["datestr"] == nothing || occursin(args["datestr"],_.target) ) |>
 @filter( _.driftfun == "dailymodel0" ) |>
 @filter( args["max_airmass"] == nothing || _.airmass <= args["max_airmass"] ) |>
 @filter( args["max_solar_hour_angle"] == nothing || abs(_.hour_angle) <= args["max_solar_hour_angle"] ) |>
 #@filter( args["start_time"] == nothing || Time(julian2datetime(_.bjd)) >= start_time ) |>
 #@filter( args["stop_time"] == nothing || Time(julian2datetime(_.bjd)) <= stop_time ) |> # TODO for other instruments may need to deal wtih cross end of 24 UTC
 @filter( args["min_expmeter"] == nothing || _.expmeter_mean >= args["min_expmeter"] ) |>
 DataFrame

if hasproperty(manifest_use,:mean_pyroflux) && hasproperty(manifest_use,:rms_pyroflux)
   manifest_use = manifest_use |>
      @filter( args["max_pyrhelio_frac_rms"] == nothing || _.rms_pyroflux <= args["max_pyrhelio_frac_rms"] * _.mean_pyroflux ) |>
      DataFrame
else
   println("# Warning:  mean_pyroflux and rms_pyroflux not found.  Reverting to expmeter_mean and expmeter_rms.")
   manifest_use = manifest_use |>
      @filter( args["max_expmeter_frac_rms"] == nothing || _.expmeter_rms <= args["max_expmeter_frac_rms"] * _.expmeter_mean ) |>
      DataFrame
end

manifest_use = manifest_use |>
   @take(args["max_num_spectra"] ) |> @orderby(_.bjd) |>
   DataFrame

if !(size(manifest_use,1) >= 1)
   println("# No usable observations.  Creating empty ", daily_rvs_filename)
   touch(daily_rvs_filename)
   exit(0)
end

obs_to_use = map(i->searchsortedfirst(manifest.bjd,manifest_use.bjd[i]),1:size(manifest_use,1) )
(num_vels, num_orders_to_use, num_obs) = size(input_data["order_ccfs"])

order_ccfs = input_data["order_ccfs"]
order_ccf_vars = input_data["order_ccf_vars"]
non_nan_mask = map(i->(!any(isnan.(view(order_ccfs,:,:,i)))) && (!any(isnan.(view(order_ccf_vars,:,:,i)))), obs_to_use)
obs_to_use = obs_to_use[non_nan_mask]
   mask = falses(1:num_obs)
   mask[obs_to_use] .= true

if !(sum(mask) >= 1)
   println("# No usable observations without any NaNs in CCFs.  Creating empty ", daily_rvs_filename)
   touch(daily_rvs_filename)
   exit(0)
end

println("# Found ", length(obs_to_use), " files of ",  size(manifest,1), " to use for RVs.")
   #df_out = manifest_use |> @rename(:drp_ccfjdmod=>:jd_drp,:drp_ccfrvmod => :rv_drp,:drp_dvrmsmod => :σrv_drp,
   df_out = manifest |> @rename(:drp_ccfjdmod=>:jd_drp,:drp_ccfrvmod => :rv_drp,:drp_dvrmsmod => :σrv_drp,
                                    :mean_pyroflux=>:pyrflux_mean, :rms_pyroflux=>:pyrflux_rms ) |> DataFrame
   df_out[!,"mask"] = mask

ccf_sum = reshape(sum(order_ccfs,dims=2),num_vels,num_obs)
ccf_var_sum = reshape(sum(input_data["order_ccf_vars"],dims=2),num_vels,num_obs)

println("# Measuring RVs from RvSpectML's CCFs.")
v_grid = collect(input_data["v_grid"])
@assert num_vels == length(v_grid)
#rvs = measure_rv_from_ccf(v_grid,view(ccf_sum,:,1),view(ccf_var_sum,:,1))

if occursin("template",args["rv_alg"])
   if isnothing(args["template_file"])
      template_ccf = view(ccf_sum,:,obs_to_use)
      template_ccf_var = view(ccf_var_sum,:,obs_to_use)
   else
      template_filename = args["template_file"]
      @assert isfile(template_filename) || islink(template_filename)
      @assert match(r"\.jld2$",template_filename) != nothing
      @assert filesize(template_filename) > 0
      println("# Reading template CCFs from ", template_filename)
      template_data = load(template_filename)
      template_order_ccfs = template_data["order_ccfs"]
      template_order_ccf_vars = template_data["order_ccf_vars"]
      @assert size(template_order_ccfs,1) == size(order_ccfs,1)
      @assert size(template_order_ccfs,2) == size(order_ccfs,2) # Could relax and check which orders are in each
      template_obs_to_use = template_data["clean_obs_mask"]
      template_num_obs = sum(template_obs_to_use)
      @assert template_num_obs >= 1
      template_ccf = reshape(sum(view(template_order_ccfs,:,:,template_obs_to_use),dims=2),num_vels,template_num_obs)
      template_ccf_var = reshape(sum(view(template_order_ccf_vars,:,:,template_obs_to_use),dims=2),num_vels,template_num_obs)
   end
   #ccf_template = EchelleCCFs.calc_ccf_template(view(ccf_sum,:,obs_to_use), view(ccf_var_sum,:,obs_to_use) )
   ccf_template = EchelleCCFs.calc_ccf_template(template_ccf, template_ccf_var,frac_of_width_to_fit=args["frac_of_width_to_fit"], measure_width_at_frac_depth=args["measure_width_at_frac_depth"])
   alg_template = MeasureRvFromCCFTemplate(v_grid=collect(v_grid),template=ccf_template)
   rvs_template = DataFrame(map(i->any(isnan.(view(ccf_sum,:,i))) ? (rv=NaN, σ_rv=NaN) : measure_rv_from_ccf(v_grid,view(ccf_sum,:,i),view(ccf_var_sum,:,i),alg=alg_template),1:num_obs))
   df_out[!,Symbol("rv_template")] = rvs_template.rv
   df_out[!,Symbol("σrv_template")] = rvs_template.σ_rv
end
if occursin("gaussian",args["rv_alg"])
   alg_gauss = MeasureRvFromCCFGaussian(frac_of_width_to_fit=args["frac_of_width_to_fit"], measure_width_at_frac_depth=args["measure_width_at_frac_depth"])
   rvs_gauss = DataFrame(map(i->any(isnan.(view(ccf_sum,:,i))) ? (rv=NaN, σ_rv=NaN)  : measure_rv_from_ccf(v_grid,view(ccf_sum,:,i),view(ccf_var_sum,:,i),alg=alg_gauss),1:num_obs))
   df_out[!,Symbol("rv_template")]= rvs_gauss.rv
   df_out[!,Symbol("σrv_gaussian")] = rvs_gauss.σ_rv
end
if occursin("quadratic",args["rv_alg"])
   alg_quad = MeasureRvFromCCFQuadratic(frac_of_width_to_fit=args["frac_of_width_to_fit"], measure_width_at_frac_depth=args["measure_width_at_frac_depth"])
   rvs_quad = DataFrame(map(i->any(isnan.(view(ccf_sum,:,i))) ? (rv=NaN, σ_rv=NaN)  : measure_rv_from_ccf(v_grid,view(ccf_sum,:,i),view(ccf_var_sum,:,i),alg=alg_quad),1:num_obs))
   df_out[!,Symbol("rv_quadratic")] = rvs_quad.rv
   df_out[!,Symbol("σrv_quadratic")] = rvs_quad.σ_rv
end

orders_physical = collect(174 .-input_data["orders_to_use"])
for j in 1:num_orders_to_use
   if iszero(maximum(view(order_ccfs,:,j,:))) || any(isnan.(view(order_ccfs,:,j,:)))   continue    end
   #println(" j = ", j, "  order_fits_idx = ", input_data["orders_to_use"][j], "  order_physical = ", orders_physical[j])
   if occursin("template",args["rv_alg"])
   # Templated-based
   if isnothing(args["template_file"])
      global template_ccf = view(order_ccfs,:,j,obs_to_use)
      global template_ccf_var = view(order_ccf_vars,:,j,obs_to_use)
   else
      global template_ccf = view(template_order_ccfs,:,j,template_obs_to_use)
      global template_ccf_var = view(template_order_ccf_vars,:,j,template_obs_to_use)
   end
   #local ccf_template = EchelleCCFs.calc_ccf_template(view(order_ccfs,:,j,obs_to_use), view(order_ccf_vars,:,j,obs_to_use) )
   local ccf_template = EchelleCCFs.calc_ccf_template(template_ccf, template_ccf_var)
   local alg_template = MeasureRvFromCCFTemplate(v_grid=collect(v_grid),template=ccf_template, frac_of_width_to_fit=args["frac_of_width_to_fit"], measure_width_at_frac_depth=args["measure_width_at_frac_depth"])
   order_rvs_template = DataFrame(map(i->any(isnan.(view(order_ccfs,:,j,i))) ? (rv=NaN, σ_rv=NaN) : measure_rv_from_ccf(v_grid,view(order_ccfs,:,j,i),view(order_ccf_vars,:,j,i),alg=alg_template),1:num_obs))
   df_out[!,Symbol("rv_" * string(orders_physical[j]) * "_template")] = order_rvs_template.rv
   df_out[!,Symbol("σrv_" * string(orders_physical[j]) * "_template")] = order_rvs_template.σ_rv
   end
   if occursin("gaussian",args["rv_alg"])
   # Gaussian fitting
   order_rvs_gauss = DataFrame(map(i->any(isnan.(view(order_ccfs,:,j,i))) ? (rv=NaN, σ_rv=NaN) : measure_rv_from_ccf(v_grid,view(order_ccfs,:,j,i),view(order_ccf_vars,:,j,i),alg=alg_gauss) ,1:num_obs))
   df_out[!,Symbol("rv_" * string(orders_physical[j]) * "_gauss")] = order_rvs_gauss.rv
   df_out[!,Symbol("σrv_" * string(orders_physical[j]) * "_gauss")] = order_rvs_gauss.σ_rv
   end
   if occursin("quadratic",args["rv_alg"])
   # Quadratic fitting
   order_rvs_quad = DataFrame(map(i->any(isnan.(view(order_ccfs,:,j,i))) ? (rv=NaN, σ_rv=NaN) : measure_rv_from_ccf(v_grid,view(order_ccfs,:,j,i),view(order_ccf_vars,:,j,i),alg=alg_quad),1:num_obs))
   df_out[!,Symbol("rv_" * string(orders_physical[j]) * "_quad")] = order_rvs_quad.rv
   df_out[!,Symbol("σrv_" * string(orders_physical[j]) * "_quad")] = order_rvs_quad.σ_rv
   end
end

println("# Writing daily RVs.")
 CSV.write(daily_rvs_filename,df_out)


text:
using PkgBenchmark
using DataFrames
using CSV
using Statistics: median
import SymbolicRegression
config = BenchmarkConfig(;juliacmd=`julia -O3`,
                          env=Dict("JULIA_NUM_THREADS" => 4))
_results = benchmarkpkg(SymbolicRegression, config; script="benchmark/.benchmarks.jl")
results = vcat([[["evaluation_"*k, median(v.times)] for (k, v) in bigV]
      for (bigK, bigV) in _results.benchmarkgroup.data]...)
df = DataFrame(commit=_results.commit,
               name=[results[i][1] for i=1:length(results)],
               time=[results[i][2] for i=1:length(results)])

CSV.write("output.csv", df; append=true)


text:
#= ========== Consul KVS ==================
Consul KVS specific code.
For the API refer https://www.consul.io/api/kv.html
=#
# using Requests

struct KVSConsul <: AbstractKVS
    name::String     # Name of db
end

function kvs_close(hdl::KVSConsul)
    return true
end

struct KVSConsulException <: Exception
    status::Int
end

struct KVSConsulCASExpired <: Exception
    tries::Int
end

function kvs_put(hdl::KVSConsul, key, value)
    k = hdl.name * base64encode(byte_array(key))
    v = byte_array(value)
    r = Requests.put("http://127.0.0.1:8500/v1/kv/$(k)"; data = v)
    (r.status != 200 || Requests.json(r)[1] == false) && throw(KVSConsulException(r.status))
    r
end

function kvs_put_atomic(hdl::KVSConsul, key, value, index)
    k = hdl.name * base64encode(byte_array(key))
    v = byte_array(value)
    r = Requests.put("http://127.0.0.1:8500/v1/kv/$(k)?cas=$(index)"; data = v)
    r.status != 200 && throw(KVSConsulException(r.status))
    Requests.json(r)[1]
end

function kvs_put_sync(hdl::KVSConsul, key, value)
    # Nothing fancy
    kvs_put(hdl, key, value)
end

function kvs_delete(hdl::KVSConsul, key)
    k = hdl.name * base64encode(byte_array(key))
    r = Requests.delete("http://127.0.0.1:8500/v1/kv/$(k)")
    r.status != 200 && return false
    true
end

function kvs_get(hdl::KVSConsul, key)
    k = hdl.name * base64encode(byte_array(key))
    r = Requests.get("http://127.0.0.1:8500/v1/kv/$(k)?raw")
    r.status != 200 && return nothing
    array_to_type(r.data)
end

"""
    kvs_get_index(hdl::KVSConsul, key)
For the given key, returns (value, index), where index is an integer
that is used to atomically update value.
"""
function kvs_get_index(hdl::KVSConsul, key)
    k = hdl.name * base64encode(byte_array(key))
    r = Requests.get("http://127.0.0.1:8500/v1/kv/$(k)")
    r.status != 200 && return (nothing, 0)
    j = Requests.json(r)[1]
    (array_to_type(base64decode(j["Value"])), j["ModifyIndex"])
end

function kvs_get_many(hdl::KVSConsul, T::Type)
    prefix = hdl.name
    r = Requests.get("http://127.0.0.1:8500/v1/kv/$(prefix)?recurse")
    r.status != 200 && return nothing
    jsn = Requests.json(r)
    keys = Vector{Any}()
    for i in jsn
        #println("$i $(base64decode(i[length(hdl.name)+1:length(i)]))")
        k = array_to_type(base64decode(i["Value"]))[1]
        push!(keys, k)
    end
    keys
end

using Test
function runtests_consul()
    @testset "Consul tests:" begin
        db = KVSConsul("mydb")
        @test db != nothing
        k = rand(UInt128)
        v = rand(UInt8, 64)
        @test kvs_put(db, k, v) != nothing
        ret = kvs_get(db, k)
        @test ret != nothing
        @test ret == v
        (ret, idx) = kvs_get_index(db, k)
        @test ret != nothing
        @test ret == v
        v1 = rand(UInt8, 64)
        @test kvs_put_atomic(db, k, v1, idx-1) == false
        @test kvs_put_atomic(db, k, v1, idx+1) == false
        @test kvs_put_atomic(db, k, v1, idx) == true
        (ret, idx1) = kvs_get_index(db, k)
        @test ret == v1
        @test idx1 != idx
    end
end


text:
# Parses distribution expression and returns expression for a segment cost function
# Required for PELT and BS macros
function cost_function(data::Any, dist_expr::Expr)
    if !Meta.isexpr(dist_expr,:call)
        error("Syntax error: expected distribution as argument")
    end

    dist_type = dist_expr.args[1]

    if dist_type == :Normal
        if length(dist_expr.args) != 3
            error("Normal distribution has two parameters")
        end
        μ, σ = dist_expr.args[2:3]
        if μ == :(:?) && σ != :(:?)
            println("Changepoint distribution is Normal with changing mean and fixed variance")
            return :(NormalMeanSegment($data , $σ))
        elseif μ != :(:?) && σ == :(:?)
            println("Changepoint distribution is Normal with fixed mean and changing variance")
            #μ = eval(Main, μ)
            return :(NormalVarSegment($data, $μ))
        elseif μ == :(:?) && σ == :(:?)
            println("Changepoint distribution is Normal with changing mean and changing variance")
            return :(NormalMeanVarSegment($data))
        else
            error("Must mark at least one Normal parameter as changing with a :? symbol")
        end

    elseif dist_type == :Exponential
        println("Changepoint distribution is Exponential with changing mean")
        return :(ExponentialSegment($data))

    elseif dist_type == :Poisson
        println("Changepoint distribution is Poisson with changing mean")
        return :(PoissonSegment($data))

    elseif dist_type == :Gamma
        if length(dist_expr.args) != 3
            error("Gamma distribution has two parameters")
        end
        alpha , beta = dist_expr.args[2:3]
        if alpha == :(:?) && beta != :(:?)
            println("Changepoint distribution is Gamma with changing shape and fixed rate")
            return :(GammaShapeSegment($data, $beta))
        elseif alpha != :(:?) && beta == :(:?)
            println("Changepoint distribution is Gamma with fixed shape and changing rate")
            return :(GammaRateSegment($data, $alpha))
        else
            error("Must mark at least one Gamma parameter as changing with a :? symbol")
        end

    elseif dist_type == :Nonparametric
        K = dist_expr.args[2]
        println("Changepoint method is Nonparametric")
        return :(NonparametricSegment($data, $K))

    elseif dist_type == :OLS
        println("Changepoints in piecewise linear regressions")
        return :(OLSSegment($data))

    elseif dist_type == :CUSUM
        println("CUSUM cost function")
        sigma = dist_expr.args[2]
        return :(CUSUM($data , $sigma))

    else
        error("Distribution $(dist_type) has no implemented cost functions")
    end
end

function cost_function_CUSUM(data::Any, sigma::Any)#dist_expr::Expr)
    println("CUSUM cost function")    #sigma = dist_expr.args[2]
    return :(CUSUM($data , $sigma))
end

macro mad(data::Any)
    println("Estimate sigma by MAD")
    return   :(mad($data) )
end

"""
    @segment_cost data changepoint_model

Creates a segment cost function given data and changepoint model expression.

# Cost functions

`changepoint_model` is an expression which describes what segment cost function
should be constructed for use with PELT. For parametric segment cost functions,
this is represented by the name of a distribution (as in the Distributions package)
with some parameters replaced by ':?' to indicate that the parameters change.
The full list of available cost functions based on parametric distributions is
as follows:

* `Normal(:?, σ)`: Normal model with changing mean and fixed standard deviation `σ` (see also [`NormalMeanSegment`](@ref))
* `Normal(μ, :?)`: Normal model with fixed mean `μ` and changing standard deviation (see also [`NormalVarSegment`](@ref))
* `Normal(:?, :?)`: Normal model with changing mean and standard deviation (see also [`NormalMeanVarSegment`](@ref))
* `Exponential(:?)`: Exponential model with changing mean (see also [`ExponentialSegment`](@ref))
* `Poisson(:?)`: Poisson model with changing mean (see also [`PoissonSegment`](@ref))
* `Gamma(:?, β)`: Gamma model with fixed rate parameter `β` and changing shape parameter (see also [`GammaShapeSegment`](@ref))
* `Gamma(α, :?)`: Gamma model with fixed shape parameter `α` and changing rate parameter (see also [`GammaRateSegment`](@ref))

The following others models are also available:

* `Nonparametric(k)`: Nonparametric cost function with parameter `k` (see also [`NonparametricSegment`](@ref))
* `OLS` : Cost function for piecewise linear regressions (see also [`OLSSegment`](@ref))

# Example
```julia-repl
n = 1000
λ = 100
μ, σ = Normal(0.0, 10.0), 1.0
# Samples changepoints from Normal distribution with changing mean
sample, cps = @changepoint_sampler n λ Normal(μ, σ)
# Create cost function
seg_cost = @segment_cost sample Normal(:?, σ)
# Calculate changepoints using PELT and BS
pelt_cps, cost = PELT(seg_cost, n)
bs_cps = BS(seg_cost, n)
```
"""
macro segment_cost(data, dist)
    esc(cost_function(data, dist))
end

"""
    @PELT data changepoint_model [β₁ [β₂] ]

Runs the PELT algorithm on time series `data` using a specified `changepoint_model` and penalties.
If no penalty `β₁` provided, a default of value `log(length(data))` is used.
If two penalties `β₁` and `β₂` are provided then the CROPS algorithm is run which finds
all optimal segmentations for all penalties between `β₁` and `β₂`.

See also: [`PELT`](@ref), [`CROPS`](@ref)

# Example
```
n = 1000
λ = 100
μ, σ = Normal(0.0, 10.0), 1.0
# Samples changepoints from Normal distribution with changing mean
sample, cps = @changepoint_sampler n λ Normal(μ, σ)
# Run PELT on sample
pelt_cps, pelt_cost = @PELT sample Normal(:?, σ)
```
"""
macro PELT(data, dist, args...)
    cost_func = cost_function(data, dist)
    if length(args) == 0
        return esc(:(PELT($(cost_func), length($data))))
    elseif length(args) == 1
        return esc(:(PELT($(cost_func), length($data), pen=$(args[1]))))
    else
        return esc(:(CROPS($(cost_func), length($data), $(args[1]), $(args[2]))))
    end
end

"""
    @BS data changepoint_model [β₁]

Runs the Binary Segmentation algorithm using a specified `changepoint_model` for a given penalty `β₁`.
If no penalty specified, use penalty `log(length(data))` by default.

See also: [`BS`](@ref), [`@segment_cost`](@ref)

# Example
```
n = 1000   # Length of time series
λ = 100    # Frequency of changepoints
α, β = Uniform(0.0, 10.0), 1.0
# Samples changepoints from Gamma distribution with changing shape
sample, cps = @changepoint_sampler n λ Gamma(α, β)
# Run binary segmentation on sample
bs_cps = @BS sample Gamma(:?, β)
```
"""
macro BS(data, dist, args...)
    cost_func = cost_function(data, dist)
    if length(args) == 0
        return esc(:(BS($(cost_func), length($data))))
    else
        return esc(:(BS($(cost_func), length($data), pen=$(args[1]))))
    end
end


text:
# Linear fractional transformation "LFT"

module LinearFractionalTransformations

import Base.inv, Base.isequal, Base.show, Base.hash
import Base.==, Base.*, Base.getindex
export LFT, isequal, call

const complex_infinity = Inf + Inf * im

struct LFT
    M::Array{Complex{Float64},2}
    function LFT(a, b, c, d)
        if isinf(a) || isinf(b) || isinf(c) || isinf(d)
            error("Arguments must be finite: " * string((a, b, c, d)))
        end

        if a * d - b * c == 0
            error("Singularity detected: " * string((a, b, c, d)))
        end

        new([a b; c d])
    end
end

"""
    LFT
A linear fractional transformation `f(z)=(a*z+b)/(c*z+d)` where `a,b,c,d` are complex numbers.

Constructors:
* `LFT(a,b,c,d)`: create the function described above.
* `LFT()`: create the identity function `f(z) = z`. Equivalent to `LFT(1,0,0,1)`.
* `LFT(a,b,c)`: create the function with `a ↦ 0`, `b ↦ 1`, and `c ↦ ∞`.
* `LFT(a,aa,b,bb,c,cc)`: creates the function with `a ↦ aa`, `b ↦ bb`, and `c ↦ cc`.
"""
function LFT(M::Array)
    return LFT(M[1, 1], M[1, 2], M[2, 1], M[2, 2])
end

function LFT()
    return LFT(1, 0, 0, 1)
end

# a --> 0, b-->1, and c--> oo
function LFT(a::Number, b::Number, c::Number)
    if a == b || b == c || a == c
        error("Three arguments must be distinct: " * string((a, b, c)))
    end

    if isinf(a)
        return LFT(0, b - c, 1, -c)
    end

    if isinf(b)
        return LFT(1, -a, 1, -c)
    end

    if isinf(c)
        return LFT(1, -a, 0, b - a)
    end

    # if all args are finite
    aa = (b - c)
    bb = (-a) * (b - c)
    cc = (b - a)
    dd = (-c) * (b - a)
    return LFT(aa, bb, cc, dd)
end

# a-->aa, b-->bb, c-->cc
function LFT(a::Number, aa::Number, b::Number, bb::Number, c::Number, cc::Number)
    f = LFT(a, b, c)
    g = LFT(aa, bb, cc)
    return inv(g) * f
end

#### Equality checking ####

# These need to be rewritten

function isequal(f::LFT, g::LFT)
    return f[0] == g[0] && f[1] == g[1] && f[Inf] == g[Inf]
    # h = f * inv(g)
    # I = h.M
    # return I[1,1]==I[2,2] && I[1,2]==I[2,1]==0
end

==(f::LFT, g::LFT) = isequal(f, g)

# Inverse transformation
function inv(L::LFT)
    a = L.M[1, 1]
    b = L.M[1, 2]
    c = L.M[2, 1]
    d = L.M[2, 2]
    return LFT(d, -b, -c, a)
end

# Composition
*(A::LFT, B::LFT) = LFT(A.M * B.M)

# Function application
function getindex(A::LFT, x::Number)
    if isinf(x)
        a = A.M[1, 1]
        b = A.M[2, 1]
        if b == 0
            return complex_infinity
        end
        return a / b
    end
    w = A.M * [x + 0im; 1.0 + 0im]
    if w[2] == 0
        return complex_infinity
    end
    return w[1] / w[2]
end

# call(A::LFT, x::Number) = A[x]

(A::LFT)(x::Number) = A[x]

function show(io::IO, L::LFT)
    print(
        io,
        "LFT( ",
        L.M[1, 1],
        " , ",
        L.M[1, 2],
        " , ",
        L.M[2, 1],
        " , ",
        L.M[2, 2],
        " )",
    )
end

function hash(f::LFT, h::UInt64 = UInt64(0))
    z = 0.0 + 0.0 * im # kludge to make -0.0 and -0.0im into +versions
    a = f[0] + z
    b = f[1] + z
    c = f[Inf] + z
    return hash(a, hash(b, hash(c, h)))
end

include("projection.jl")

end # end of module "LFTs"


text:
const SVMLinearClassifier_ = sksv(:LinearSVC)
@sk_clf mutable struct SVMLinearClassifier <: MMI.Deterministic
    penalty::String = "l2"::(_ in ("l1","l2"))
    loss::String    = "squared_hinge"::(_ in ("hinge", "squared_hinge"))
    dual::Bool      = true
    tol::Float64    = 1e-4::(_ > 0)
    C::Float64      = 1.0::(_ > 0)
    multi_class::String        = "ovr"::(_ in ("crammer_singer", "ovr"))
    fit_intercept::Bool        = true
    intercept_scaling::Float64 = 1.0::(_ > 0)
    random_state::Any          = nothing
    max_iter::Int              = 1000::(_ > 0)
end
MMI.fitted_params(m::SVMLinearClassifier, (f, _, _)) = (
    coef      = f.coef_,
    intercept = f.intercept_,
    classes   = f.classes_
    )
meta(SVMLinearClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    descr   = "Linear Support Vector Classification."
    )

# ----------------------------------------------------------------------------
const SVMClassifier_ = sksv(:SVC)
@sk_clf mutable struct SVMClassifier <: MMI.Deterministic
    C::Float64          = 1.0::(_ > 0)
    kernel::Union{String,Function}  = "rbf"::(_ isa Function || _ in ("linear", "poly", "rbf", "sigmoid", "precomputed"))
    degree::Int         = 3::(_ > 0)
    gamma::Union{Float64,String}    = "auto"::(_ isa Float64 && _ > 0 || _ in ("scale", "auto"))
    coef0::Float64      = 0.0
    shrinking::Bool     = true
    tol::Float64        = 1e-3::(_ > 0)
    cache_size::Int     = 200::(_ > 0)
    max_iter::Int       = -1
    decision_function_shape::String = "ovr"::(_ in ("ovo", "ovr"))
    random_state        = nothing
end
MMI.fitted_params(m::SVMClassifier, (f, _, _)) = (
    support         = f.support_,
    support_vectors = f.support_vectors_,
    n_support       = f.n_support_,
    dual_coef       = f.dual_coef_,
    coef            = m.kernel == "linear" ? f.coef_ : nothing,
    intercept       = f.intercept_,
    fit_status      = f.fit_status_,
    classes         = f.classes_
    # probA = f.probA_,
    # probB = f.probB_,
    )
meta(SVMClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    descr   = "C-Support Vector Classification."
    )

# ============================================================================
const SVMLinearRegressor_ = sksv(:LinearSVR)
@sk_reg mutable struct SVMLinearRegressor <: MMI.Deterministic
    epsilon::Float64 = 0.0::(_ >= 0)
    tol::Float64     = 1e-4::(_ > 0)
    C::Float64       = 1.0::(_ > 0)
    loss::String      = "epsilon_insensitive"::(_ in ("epsilon_insensitive", "squared_epsilon_insensitive"))
    fit_intercept::Bool        = true
    intercept_scaling::Float64 = 1.0::(_ > 0)
    dual::Bool                 = true
    random_state::Any          = nothing
    max_iter::Int              = 1000::(_ > 0)
end
MMI.fitted_params(m::SVMLinearRegressor, (f, _, _)) = (
    coef = f.coef_,
    intercept = f.intercept_
    )
meta(SVMLinearRegressor,
    input   = Table(Continuous),
    target  = AbstractVector{Continuous},
    descr   = "Linear Support Vector Regression."
    )

# ----------------------------------------------------------------------------
const SVMRegressor_ = sksv(:SVR)
@sk_reg mutable struct SVMRegressor <: MMI.Deterministic
    kernel::Union{String,Function} = "rbf"::(_ isa Function || _ in ("linear", "poly", "rbf", "sigmoid", "precomputed"))
    degree::Int      = 3::(_ > 0)
    gamma::Union{Float64,String}   = "auto"::(_ isa Float64 && _ > 0 || _ in ("scale", "auto"))
    coef0::Float64   = 0.0
    tol::Float64     = 1e-3::(_ > 0)
    C::Float64       = 1.0::(_ > 0)
    epsilon::Float64 = 0.1::(_ >= 0)
    shrinking        = true
    cache_size::Int  = 200::(_ > 0)
    max_iter::Int    = -1
end
MMI.fitted_params(m::SVMRegressor, (f, _, _)) = (
    support         = f.support_,
    support_vectors = f.support_vectors_,
    dual_coef       = f.dual_coef_,
    coef            = m.kernel == "linear" ? f.coef_ : nothing,
    intercept       = f.intercept_,
    fit_status      = f.fit_status_
    )
meta(SVMRegressor,
    input   = Table(Continuous),
    target  = AbstractVector{Continuous},
    descr   = "Epsilon-Support Vector Regression."
    )

# ============================================================================
const SVMNuClassifier_ = sksv(:NuSVC)
@sk_clf mutable struct SVMNuClassifier <: MMI.Deterministic
    nu::Float64 = 0.5::(0 < _ <= 1)
    kernel::Union{String,Function}  = "rbf"::(_ isa Function || _ in ("linear", "poly", "rbf", "sigmoid", "precomputed"))
    degree::Int         = 3::(_ > 0)
    gamma::Union{Float64,String}    = "auto"::(_ isa Float64 && _ > 0 || _ in ("scale", "auto"))
    coef0::Float64      = 0.0
    shrinking::Bool     = true
    # probability
    tol::Float64        = 1e-3::(_ > 0)
    cache_size::Int     = 200::(_ > 0)
    max_iter::Int       = -1
    decision_function_shape::String = "ovr"::(_ in ("ovo", "ovr"))
    # break_ties::Bool    = false # XXX only >= 0.22
    random_state        = nothing
end
MMI.fitted_params(m::SVMNuClassifier, (f, _, _)) = (
    support         = f.support_,
    support_vectors = f.support_vectors_,
    n_support       = f.n_support_,
    dual_coef       = f.dual_coef_,
    coef            = m.kernel == "linear" ? f.coef_ : nothing,
    intercept       = f.intercept_,
    fit_status      = f.fit_status_,
    classes         = f.classes_
    # probA = f.probA_,
    # probB = f.probB_,
    )
meta(SVMNuClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    descr   = "Nu-Support Vector Classification."
    )

# ============================================================================
const SVMNuRegressor_ = sksv(:NuSVR)
@sk_reg mutable struct SVMNuRegressor <: MMI.Deterministic
    nu::Float64 = 0.5::(0 < _ <= 1)
    C::Float64       = 1.0::(_ > 0)
    kernel::Union{String,Function} = "rbf"::(_ isa Function || _ in ("linear", "poly", "rbf", "sigmoid", "precomputed"))
    degree::Int      = 3::(_ > 0)
    gamma::Union{Float64,String}   = "auto"::(_ isa Float64 && _ > 0 || _ in ("scale", "auto"))
    coef0::Float64   = 0.0
    shrinking        = true
    tol::Float64     = 1e-3::(_ > 0)
    cache_size::Int  = 200::(_ > 0)
    max_iter::Int    = -1
end
MMI.fitted_params(m::SVMNuRegressor, (f, _, _)) = (
    support         = f.support_,
    support_vectors = f.support_vectors_,
    dual_coef       = f.dual_coef_,
    coef            = m.kernel == "linear" ? f.coef_ : nothing,
    intercept       = f.intercept_,
    )
meta(SVMNuRegressor,
    input   = Table(Continuous),
    target  = AbstractVector{Continuous},
    descr   = "Nu-Support Vector Regression."
    )