text:
export RobotBasedSystemMatrixProtocol, RobotBasedSystemMatrixProtocolParams
"""
Parameters for the RobotBasedSystemMatrixProtocol
"""
Base.@kwdef mutable struct RobotBasedSystemMatrixProtocolParams <: RobotBasedProtocolParams
  "Minimum wait time between robot movements"
  waitTime::Float64
  "Number of background frames to measure for a background position"
  bgFrames::Int64
  "Number of frames that are averaged for a foreground position"
  fgFrameAverages::Int64
  "Flag if the calibration should be saved as a system matrix or not"
  saveAsSystemMatrix::Bool = true
  "Number of background measurements to take"
  bgMeas::Int64 = 0
  "If set the tx amplitude and phase will be set with control steps"
  controlTx::Bool = false
  "Sequence used for the calibration at each position"
  sequence::Union{Sequence, Nothing} = nothing
  positions::Union{GridPositions, Nothing} = nothing
end
function RobotBasedSystemMatrixProtocolParams(dict::Dict, scanner::MPIScanner)
  if haskey(dict, "Positions")
    posDict = dict["Positions"]

    positions = Positions(posDict)
    delete!(dict, "Positions")
  else
    positions = nothing
  end

  sequence = nothing
  if haskey(dict, "sequence")
    sequence = Sequence(scanner, dict["sequence"])
    dict["sequence"] = sequence
    delete!(dict, "sequence")
  end

  params = params_from_dict(RobotBasedSystemMatrixProtocolParams, dict)
  params.positions = positions
  params.sequence = sequence
  return params
end

# Based on https://github.com/MagneticParticleImaging/MPIMeasurements.jl/tree/cde1c72b820a72b3c3dfa4235b2b37bd506b0109
mutable struct SystemMatrixRobotMeas
  task::Union{Task,Nothing}
  consumer::Union{Task, Nothing}
  producerFinalizer::Union{Task, Nothing}
  #store::DatasetStore # Do we need this?
  positions::GridPositions
  currPos::Int
  signals::Array{Float32,4}
  currentSignal::Array{Float32,4}
  measIsBGPos::Vector{Bool}
  posToIdx::Vector{Int64}
  measIsBGFrame::Vector{Bool}
  temperatures::Matrix{Float64}
end

Base.@kwdef mutable struct RobotBasedSystemMatrixProtocol <: Protocol
  name::AbstractString
  description::AbstractString
  scanner::MPIScanner
  params::RobotBasedSystemMatrixProtocolParams
  biChannel::Union{BidirectionalChannel{ProtocolEvent}, Nothing} = nothing

  executeTask::Union{Task, Nothing} = nothing
  systemMeasState::Union{SystemMatrixRobotMeas, Nothing} = nothing
  txCont::Union{TxDAQController, Nothing} = nothing
  stopped::Bool = false
  cancelled::Bool = false
  restored::Bool = false
  finishAcknowledged::Bool = false
end

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

function requiredDevices(protocol::RobotBasedSystemMatrixProtocol)
  result = [AbstractDAQ, Robot, SurveillanceUnit]
  if protocol.params.controlTx
    push!(result, TxDAQController)
  end
  return result
end

function _init(protocol::RobotBasedSystemMatrixProtocol)
  if isnothing(protocol.params.sequence)
    throw(IllegalStateException("Protocol requires a sequence"))
  end
  protocol.systemMeasState = SystemMatrixRobotMeas()

  # Prepare Positions
  # Extend Positions to include background measurements, TODO behaviour if positions already includes background pos
  cartGrid = protocol.params.positions
  if protocol.params.bgMeas == 0
    positions = cartGrid
  else
    bgIdx = round.(Int64, range(1, stop=length(cartGrid)+protocol.params.bgMeas, length=protocol.params.bgMeas ) )
    bgPos = namedPosition(getRobot(protocol.scanner),"park")
    positions = BreakpointGridPositions(cartGrid, bgIdx, bgPos)
  end
  protocol.params.positions = positions
  measIsBGPos = isa(protocol.params.positions,BreakpointGridPositions) ? MPIFiles.getmask(protocol.params.positions) : zeros(Bool,length(protocol.params.positions))
  numBGPos = sum(measIsBGPos)
  numFGPos = length(measIsBGPos) - numBGPos
  numTotalFrames = numFGPos + protocol.params.bgFrames*numBGPos
  # The following looks like a secrete line but it makes sense
  posToIdx = cumsum(vcat([false],measIsBGPos)[1:end-1] .* (protocol.params.bgFrames - 1) .+ 1)
  measIsBGFrame = zeros(Bool, numTotalFrames)

  protocol.systemMeasState.measIsBGPos = measIsBGPos
  protocol.systemMeasState.posToIdx = posToIdx
  protocol.systemMeasState.measIsBGFrame = measIsBGFrame
  protocol.systemMeasState.currPos = 1
  protocol.systemMeasState.positions = protocol.params.positions # TODO remove redundancy

  if !checkPositions(protocol)
    throw(IllegalStateException("Protocol has illegal positions"))
  end

  #Prepare Signals
  numRxChannels = length(rxChannels(protocol.params.sequence)) # kind of hacky, but actual rxChannels for RedPitaya are only set when setupRx is called
  rxNumSamplingPoints = rxNumSamplesPerPeriod(protocol.params.sequence)
  numPeriods = acqNumPeriodsPerFrame(protocol.params.sequence)
  filenameSignals = "/tmp/sysObj.bin"
  io = open(filenameSignals, "w+");
  signals = Mmap.mmap(io, Array{Float32,4}, (rxNumSamplingPoints,numRxChannels,numPeriods,numTotalFrames));
  signals[:] .= 0.0 # Does this not erase our stored .bin in recovery case?
  protocol.systemMeasState.signals = signals

  protocol.systemMeasState.currentSignal = zeros(Float32,rxNumSamplingPoints,numRxChannels,numPeriods,1)

  # TODO implement properly
  protocol.systemMeasState.temperatures = zeros(0, 0)

  # Init TxDAQController
  if protocol.params.controlTx
    controllers = getDevices(protocol.scanner, TxDAQController)
    if length(controllers) > 1
      throw(IllegalStateException("Cannot unambiguously find a TxDAQController as the scanner has $(length(controllers)) of them"))
    end
    protocol.txCont = controllers[1]
    protocol.txCont.currTx = nothing
  else
    protocol.txCont = nothing
  end
end

function checkPositions(protocol::RobotBasedSystemMatrixProtocol)
  rob = getRobot(protocol.scanner)
  valid = true
  if hasDependency(rob, AbstractCollisionModule)
    cms = dependencies(rob, AbstractCollisionModule)
    for cm in cms
      valid &= all(checkCoords(cm, protocol.params.positions))
    end
  end
  for pos in protocol.params.positions
    valid &= checkAxisRange(rob, toRobotCoords(rob, ScannerCoords(pos)))
  end
  return valid
end

function enterExecute(protocol::RobotBasedSystemMatrixProtocol)
  protocol.stopped = false
  protocol.cancelled = false
  protocol.finishAcknowledged = false
  protocol.restored = false
end

function _execute(protocol::RobotBasedSystemMatrixProtocol)
  @info "Start System Matrix Protocol"
  if !isReferenced(getRobot(protocol.scanner))
    throw(IllegalStateException("Robot not referenced! Cannot proceed!"))
  end

  finished = false
  started = false
  notifiedStop = false
  while !finished
    if !started
      if isfile("/tmp/sysObj.toml")
        message = """Found existing calibration file! \n
        Should it be resumed?"""
        if askConfirmation(protocol, message)
          restore(protocol)
        end
      end
    end
    finished = performCalibration(protocol)
    started = true

    # Stopped
    notifiedStop = false
    while protocol.stopped
      handleEvents(protocol)
      protocol.cancelled && throw(CancelException())
      if !notifiedStop
        put!(protocol.biChannel, OperationSuccessfulEvent(StopEvent()))
        notifiedStop = true
      end
      if !protocol.stopped
        put!(protocol.biChannel, OperationSuccessfulEvent(ResumeEvent()))
      end
      sleep(0.05)
    end
  end

  put!(protocol.biChannel, FinishedNotificationEvent())
  while !protocol.finishAcknowledged
    handleEvents(protocol)
    protocol.cancelled && throw(CancelException())
    sleep(0.01)
  end
  @info "Protocol finished."
  close(protocol.biChannel)
end

function cleanup(protocol::RobotBasedSystemMatrixProtocol)
  # TODO should cleanup remove temp files? Would require a handler to differentiate between successful and unsuccesful "end"
  removeTempFiles(protocol)
end

function removeTempFiles(protocol::RobotBasedSystemMatrixProtocol)
  filename = "/tmp/sysObj.toml"
  filenameSignals = "/tmp/sysObj.bin"
  rm(filename, force=true)
  rm(filenameSignals, force=true)
end

function performCalibration(protocol::RobotBasedSystemMatrixProtocol)
  @info "Enter calibration loop"
  finished = false
  calib = protocol.systemMeasState
  su = getSurveillanceUnit(protocol.scanner)
  daq = getDAQ(protocol.scanner)
  robot = getRobot(protocol.scanner)

  positions = calib.positions
  numPos = length(calib.positions)
  @info "Store SF"

  stopTx(daq)
  while true
    @info "Curr Pos in performCalibrationInner $(calib.currPos)"
    handleEvents(protocol)

    if protocol.stopped
      wait(calib.consumer)
      wait(calib.producerFinalizer)
      @info "Stop calibration loop"
      finished = false
      break
    end

    if calib.currPos <= numPos
      pos = ScannerCoords(uconvert.(Unitful.mm, positions[calib.currPos]))
      performCalibration(protocol, pos)
      calib.currPos +=1
    end

    if calib.currPos > numPos
      wait(calib.consumer)
      wait(calib.producerFinalizer)
      stopTx(daq)
      disableACPower(su)
      enable(robot)
      movePark(robot)
      disable(robot)

      finished = true
      break
    end

  end
    @info "Exit calibration loop"
  return finished
end

function performCalibration(protocol::RobotBasedSystemMatrixProtocol, pos)
  calib = protocol.systemMeasState
  robot = getRobot(protocol.scanner)

  enable(robot)
  timePreparing = 0
  try
    timePreparing = @elapsed prepareMeasurement(protocol, pos) # TODO params
  catch ex
    if ex isa CompositeException
      @error "CompositeException while preparing measurement:"
      for e in ex
        @error e
      end
    end
    rethrow(ex)
  end

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

  disable(robot)

  timeMeasuring = @elapsed measurement(protocol)

  @info "Preptime $timePreparing, meas time: $(timeMeasuring)"

end

function prepareMeasurement(protocol::RobotBasedSystemMatrixProtocol, pos)
  calib = protocol.systemMeasState
  robot = getRobot(protocol.scanner)
  daq = getDAQ(protocol.scanner)
  su = getSurveillanceUnit(protocol.scanner)
  timeMove = 0
  timePrepDAQ = 0
  timeFinalizer = 0
  timeFrameChange = 0
  timeSeq = 0
  timeTx = 0
  timeConsumer = 0
  timeWaitSU = 0

  @sync begin
    # Prepare Robot/Sample
    timeMove = @elapsed moveAbs(robot, pos)

    # Prepare DAQ
    @async timePrepDAQ = @elapsed begin
      allowControlLoop = mod1(calib.currPos, 11) == 1
      timeFinalizer = @elapsed wait(calib.producerFinalizer)
      # The following tasks can only be started after the finalizer and mostly only in this order
      suTask = @async enableACPower(su)

      timeFrameChange = @elapsed begin
        if protocol.restored || (calib.currPos == 1) || (calib.measIsBGPos[calib.currPos] != calib.measIsBGPos[calib.currPos-1])
          acqNumFrames(protocol.params.sequence, calib.measIsBGPos[calib.currPos] ? protocol.params.bgFrames : 1)
          acqNumFrameAverages(protocol.params.sequence, calib.measIsBGPos[calib.currPos] ? 1 : protocol.params.fgFrameAverages)
          setup(daq, protocol.params.sequence) #TODO setupTx might be fine once while setupRx needs to be done for each new sequence
          setSequenceParams(daq, protocol.params.sequence)
          protocol.restored = false
        end
      end

      timeSeq = @elapsed prepareSequence(daq, protocol.params.sequence)
      timeWaitSU = @elapsed wait(suTask)
      # TODO check again if controlLoop can be run while robot is active
      timeTx = @elapsed begin
        if protocol.params.controlTx
          if allowControlLoop
            controlTx(protocol.txCont, protocol.params.sequence, protocol.txCont.currTx)
          else
            setTxParams(daq, txFromMatrix(protocol.txCont, protocol.txCont.currTx)...)
          end
        else
          prepareTx(daq, protocol.params.sequence)
        end
      end
    end

    @async timeConsumer = @elapsed wait(calib.consumer)
  end
  @info "############### Preparing: Move $timeMove Prep DAQ time $timePrepDAQ, Finalizer $timeFinalizer, Frame $timeFrameChange, Seq $timeSeq, SU $timeWaitSU, Tx $timeTx, Consumer $timeConsumer"
end

function measurement(protocol::RobotBasedSystemMatrixProtocol)
  index = protocol.systemMeasState.currPos
  calib = protocol.systemMeasState
  @info "Measurement" index length(calib.positions)

  timeGetThings = @elapsed begin
    # TODO getSafety and getTempSensor if necessary
    #safety = getSafety(protocol.scanner)
    daq = getDAQ(protocol.scanner)
    su = getSurveillanceUnit(protocol.scanner)
    channel = Channel{channelType(daq)}(32)
  end

  @info "Starting Measurement"
  timeEnableSlowDAC = @elapsed begin
    calib.consumer = @tspawnat protocol.scanner.generalParams.consumerThreadID asyncConsumer(channel, protocol)
    producer = @tspawnat protocol.scanner.generalParams.producerThreadID asyncProducer(channel, daq, protocol.params.sequence, prepTx = false, prepSeq = false, endSeq = false)
    while !istaskdone(producer)
      handleEvents(protocol)
      # Dont want to throw cancel here
      sleep(0.05)
    end
  end
  close(channel)

  @show timeEnableSlowDAC
  start, endFrame = getFrameTiming(daq)
  calib.producerFinalizer = @async begin
    endSequence(daq, endFrame)
    disableACPower(su)
  end
end

function asyncConsumer(channel::Channel, protocol::RobotBasedSystemMatrixProtocol)
  index = protocol.systemMeasState.currPos
  calib = protocol.systemMeasState
  @info "readData"
  daq = getDAQ(protocol.scanner)
  #tempSensor = getTemperatureSensor(protocol.scanner)
  tempSensor = nothing
  asyncBuffer = AsyncBuffer(daq)
  numFrames = acqNumFrames(protocol.params.sequence)
  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(protocol.scanner))

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

  if calib.measIsBGPos[index]
    calib.signals[:,:,:,startIdx:stopIdx] = uMeas
  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(protocol)
  @info "done after $timeStore"
end

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

  sysObj = protocol.systemMeasState
  params = MPIFiles.toDict(sysObj.positions)
  params["currPos"] = sysObj.currPos
  #params["stopped"] = protocol.stopped
  #params["currentSignal"] = sysObj.currentSignal
  params["waitTime"] = protocol.params.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(protocol::RobotBasedSystemMatrixProtocol)
  filename = "/tmp/sysObj.toml"
  filenameSignals = "/tmp/sysObj.bin"

  sysObj = protocol.systemMeasState
  params = MPIFiles.toDict(sysObj.positions)
  if isfile(filename)
    params = TOML.parsefile(filename)
    sysObj.currPos = params["currPos"]
    protocol.stopped = false
    #params["stopped"]
    #sysObj.currentSignal = params["currentSignal"]
    protocol.params.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)

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

    numTotalFrames = numFGPos + protocol.params.bgFrames*numBGPos
    seq = protocol.params.sequence
    numRxChannels = length(rxChannels(seq)) # kind of hacky, but actual rxChannels for RedPitaya are only set when setupRx is called
    rxNumSamplingPoints = rxNumSamplesPerPeriod(seq)
    numPeriods = acqNumPeriodsPerFrame(seq)

    io = open(filenameSignals, "r+");
    sysObj.signals = Mmap.mmap(io, Array{Float32,4},
          (rxNumSamplingPoints,numRxChannels, numPeriods, numTotalFrames));
    protocol.restored = true
    @info "Restored system matrix measurement"
  end
end

function stop(protocol::RobotBasedSystemMatrixProtocol)
  calib = protocol.systemMeasState
  if calib.currPos <= length(calib.positions)
    # OperationSuccessfulEvent is put when it actually is in the stop loop
    protocol.stopped = true
  else
    # Stopped has no concept once all measurements are done
    put!(protocol.biChannel, OperationUnsuccessfulEvent(StopEvent()))
  end
end

function resume(protocol::RobotBasedSystemMatrixProtocol)
  protocol.stopped = false
  # OperationSuccessfulEvent is put when it actually leaves the stop loop
end

function cancel(protocol::RobotBasedSystemMatrixProtocol)
  protocol.cancelled = true # Set cancel s.t. exception can be thrown when appropiate
  protocol.stopped = true # Set stop to reach a known/save state
end

function handleEvent(protocl::RobotBasedSystemMatrixProtocol, event::ProgressQueryEvent)
  put!(protocl.biChannel, ProgressEvent(protocl.systemMeasState.currPos, length(protocl.systemMeasState.positions), "Position", event))
end

function handleEvent(protocol::RobotBasedSystemMatrixProtocol, event::DataQueryEvent)
  put!(protocol.biChannel, DataAnswerEvent(protocol.systemMeasState.currentSignal, event))
end

function handleEvent(protocol::RobotBasedSystemMatrixProtocol, event::DatasetStoreStorageRequestEvent)
  if false
    # TODO this should be some sort of storage failure event
    put!(protocol.biChannel, IllegaleStateEvent("Calibration measurement is not done yet. Cannot save!"))
  else
    store = event.datastore
    scanner = protocol.scanner
    params = event.params
    data = protocol.systemMeasState.signals
    positions = protocol.systemMeasState.positions
    isBackgroundFrame = protocol.systemMeasState.measIsBGFrame
    params["storeAsSystemMatrix"] = protocol.params.saveAsSystemMatrix
    filename = saveasMDF(store, scanner, protocol.params.sequence, data, positions, isBackgroundFrame, params)
    @show filename
    put!(protocol.biChannel, StorageSuccessEvent(filename))
  end
end

handleEvent(protocol::RobotBasedSystemMatrixProtocol, event::FinishedAckEvent) = protocol.finishAcknowledged = true

protocolInteractivity(protocol::RobotBasedSystemMatrixProtocol) = Interactive()


text:
@testset "(In)Finite" begin
    cv = categorical(['x', 'y'])
    c  = cv[1]
    uv = categorical(['x', 'y'], ordered=true)
    u  = uv[1]
    @test scitype((4, 4.5, c, u, "X")) == Tuple{
                       Count,Continuous,Multiclass{2},OrderedFactor{2},Textual}
end

A = Any[2 4.5;
        4 4.5;
        6 4.5]

@testset "scitype_union" begin
    @test scitype_union(A) == Union{Count,Continuous}
    @test scitype_union(randn(1000000)) == Continuous
    @test scitype_union(1) == Count
    @test scitype_union([1]) == Count
    @test scitype_union(Any[1]) == Count
    @test scitype_union([1, 2.0, "3"]) == Union{Continuous, Count, Textual}
end

@testset "elscitype" begin
    X = randn(5, 5)
    @test scitype(X) == Arr{Continuous, 2}
    @test elscitype(X) == Continuous
    X = [1, 2, missing, 5]
    @test elscitype(X) == Union{Missing,Count}
    X = categorical([1, 2, 3, missing, 2, 1, missing])
    @test elscitype(X) == Union{Missing,Multiclass{3}}
    X = categorical("lksfjalksjdflkjsdlkfjasldkfj" |> collect)
    @test elscitype(X) <: Multiclass
    X = coerce(categorical([1,2,3,1,2,3]), OrderedFactor)
    @test elscitype(X) <: Union{Missing,OrderedFactor}
end

@testset "Arrays" begin
    # COUNT
    @test scitype(A)               == Arr{Union{Count,Continuous}, 2}
    @test scitype([1, 2, 3])       == Vec{Count}
    @test scitype([1, missing, 3]) == Vec{Union{Missing,Count}}
    @test scitype(Any[1, 2, 3])    == Vec{Count}
    @test scitype(Any[1, missing]) == Vec{Union{Missing,Count}}

    # CONTINUOUS
    @test scitype([1.0, 2.0, 3.0])    == Vec{Continuous}
    @test scitype(Any[1.0, missing])  == Vec{Union{Missing,Continuous}}
    @test scitype(Any[1.0, 2.0, 3.0]) == Vec{Continuous}
    @test scitype(Any[1.0, missing])  == Vec{Union{Missing,Continuous}}

    # MULTICLASS
    @test scitype(categorical(1:4))         == Vec{Multiclass{4}}
    @test scitype(view(categorical(1:4), 1:3)) == Vec{Multiclass{4}}
    @test scitype(Any[categorical(1:4)...]) == Vec{Multiclass{4}}
    @test scitype(categorical([1, missing, 3])) ==
        Vec{Union{Multiclass{2},Missing}}
    # NOTE: the slice here does not contain missings but the machine type
    # still contains a missing so the scitype remains with a missing
    @test scitype(categorical([1, missing, 3])[1:1]) ==
        Vec{Union{Multiclass{2},Missing}}

    # ORDEREDFACTOR
    @test scitype(categorical(1:4, ordered=true)) ==
        Vec{OrderedFactor{4}}
    @test scitype(Any[categorical(1:4, ordered=true)...]) ==
        Vec{OrderedFactor{4}}
    @test scitype(categorical([1, missing, 3], ordered=true)) ==
        Vec{Union{OrderedFactor{2},Missing}}
    # NOTE: see note for multiclass
    @test scitype(categorical([1, missing, 3], ordered=true)[1:1]) ==
        Vec{Union{OrderedFactor{2},Missing}}
end

@testset "Images" begin
    black = RGB(0, 0, 0)
    color_image = fill(black, (10, 20))
    @test scitype(color_image) == ColorImage{10,20}

    color_image2 = fill(black, (5, 3))
    v = [color_image, color_image2, color_image2]
    @test scitype(v) ==
        Vec{Union{ColorImage{10,20},ColorImage{5,3}}}

    white = Gray(1.0)
    gray_image = fill(white, (10, 20))
    @test scitype(gray_image) == GrayImage{10,20}
end

@testset "PersistenceDiagrams" begin
    diagram = PersistenceDiagram([(1, Inf), (2, 3)], dim=0)
    @test scitype(diagram) == PersistenceDiagram

    diagrams = [diagram, diagram, diagram]
    @test scitype(diagrams) == Vec{PersistenceDiagram}
end

@testset "temporal types" begin
    d = Date(2020, 4, 21)
    t = Time(8, 15, 42)
    dt = now()
    @test scitype(d) == ScientificDate
    @test scitype(t) == ScientificTime
    @test scitype(dt) == ScientificDateTime
    @test scitype(vcat([missing, ], fill(d, 3))) ==
        AbstractVector{Union{Missing,ScientificDate}}
    @test scitype(vcat([missing, ], fill(t, 3))) ==
        AbstractVector{Union{Missing,ScientificTime}}
    @test scitype(vcat([missing, ], fill(dt, 3))) ==
        AbstractVector{Union{Missing,ScientificDateTime}}
    @test scitype([missing, d, t]) ==
        AbstractVector{Union{Missing,ScientificTimeType}}
end

@testset "Type coercion" begin
    X = (x=10:10:44, y=1:4, z=collect("abcd"))

    @test_throws ArgumentError coerce(X, Count())
    @test_throws ArgumentError coerce(X, Dict(:x=>Count(), :z=>Multiclass))
    @test_throws ArgumentError coerce(X, :x=>Count(), :z=>Multiclass)

    # test fix for issue 39
    y = collect(Int64, 1:5)
    @test_throws MLJScientificTypes.CoercionError coerce(y, Float64)
    @test_throws MLJScientificTypes.CoercionError coerce(y, Textual)
    @test_throws MLJScientificTypes.CoercionError coerce(X, :x=>Float64)
    @test_throws MLJScientificTypes.CoercionError coerce(X, :x=>Textual)

    types = Dict(:x => Continuous, :z => Multiclass)
    X_coerced = coerce(X, types)
    @test X_coerced ==  coerce(X, :x => Continuous, :z => Multiclass)
    @test scitype_union(X_coerced.x) === Continuous
    @test scitype_union(X_coerced.z) <: Multiclass
    @test !X_coerced.z.pool.ordered
    @test_throws MethodError coerce(["a", "b", "c"], Count)

    y = collect(Float64, 1:5)
    y_coerced = coerce(y, Count)
    @test scitype_union(y_coerced) === Count
    @test y_coerced == y
    y = [1//2, 3//4, 6//5]
    y_coerced = coerce(y, Continuous)
    @test scitype_union(y_coerced) === Continuous
    @test y_coerced ≈ y
    X_coerced = coerce(X, Dict(:z => OrderedFactor))
    @test X_coerced == coerce(X, :z => OrderedFactor)
    @test X_coerced.x === X.x
    @test scitype_union(X_coerced.z) <: OrderedFactor
    @test X_coerced.z.pool.ordered

    # Check no-op coercion
    y = rand(Float64, 5)
    @test coerce(y, Continuous) === y
    y = rand(Float32, 5)
    @test coerce(y, Continuous) === y
    y = rand(BigFloat, 5)
    @test coerce(y, Continuous) === y
    y = rand(Int, 5)
    @test coerce(y, Count) === y
    y = big.(y)
    @test coerce(y, Count) === y
    y = rand(UInt32, 5)
    @test coerce(y, Count) === y
    X_coerced = coerce(X, Dict{Symbol, Type}())
    @test X_coerced.x === X.x
    @test X_coerced.z === X.z
    z = categorical(X.z)
    @test coerce(z, Multiclass) === z
    z = categorical(X.z, compress=true, ordered=false)
    @test coerce(z, Multiclass) === z
    z = categorical(X.z, compress=true, ordered=true)
    @test coerce(z, OrderedFactor) === z

    # missing values
    y_coerced = @test_logs(
        (:info, r"Trying to coerce from `Union{Missing,"),
        coerce([4, 7, missing], Continuous))
    @test ismissing(y_coerced == [4.0, 7.0, missing])
    @test scitype_union(y_coerced) === Union{Missing,Continuous}
    y_coerced = @test_logs(
        (:info, r"Trying to coerce from `Any"),
        coerce(Any[4, 7.0, missing], Continuous))
    @test ismissing(y_coerced == [4.0, 7.0, missing])
    @test scitype_union(y_coerced) === Union{Missing,Continuous}
    y_coerced = @test_logs(
        (:info, r"Trying to coerce from `Union{Missing,"),
        coerce([4.0, 7.0, missing], Count))
    @test ismissing(y_coerced == [4, 7, missing])
    @test scitype_union(y_coerced) === Union{Missing,Count}
    y_coerced = @test_logs(
        (:info, r"Trying to coerce from `Any"),
        coerce(Any[4, 7.0, missing], Count))
    @test ismissing(y_coerced == [4, 7, missing])
    @test scitype_union(y_coerced) === Union{Missing,Count}
    @test scitype_union(@test_logs(
        (:info, r"Trying to coerce from `Union{Missing,"),
        coerce(['x', 'y', missing], Multiclass))) ===
            Union{Missing, Multiclass{2}}
    @test scitype_union(@test_logs(
        (:info, r"Trying to coerce from `Union{Missing,"),
        coerce(['x', 'y', missing], OrderedFactor))) ===
            Union{Missing, OrderedFactor{2}}
    # non-missing Any vectors
    @test coerce(Any[4, 7],     Continuous) == [4.0, 7.0]
    @test coerce(Any[4.0, 7.0], Continuous) == [4, 7]

    # Finite conversions:
    @test scitype_union(coerce(['x', 'y'], Finite)) === Multiclass{2}
    @test scitype_union(@test_logs(
        (:info, r"Trying to coerce from `Union{Missing,"),
        coerce(['x', 'y', missing], Finite))) === Union{Missing, Multiclass{2}}

    # More finite conversions (to check resolution of #48):
    y = categorical([1, 2, 3, missing]) # unordered
    yc = coerce(y, Union{Missing,OrderedFactor})
    @test isordered(yc)
    @test yc[1].pool.ordered
    @test scitype(yc) == Vec{Union{Missing,OrderedFactor{3}}}
    @test scitype_union(yc) == Union{Missing,OrderedFactor{3}}
    @test scitype_union(y) == Union{Missing,Multiclass{3}}

end

@testset "coerce arrays" begin
    A = rand(Int, 2, 3)
    z = rand(Char, 2, 3)
    y = Any[1.0 2; 3 4]
    @test scitype_union(coerce(A, Continuous)) == Continuous
    @test scitype_union(coerce(A, OrderedFactor)) <: OrderedFactor
    @test scitype_union(coerce(z, Multiclass)) <: Multiclass
    @test scitype_union(coerce(y, Count)) === Count
end

@testset "Real->OF" begin
    v = [0.1, 0.2, 0.2, 0.3, missing, 0.1]
    w = [0.1, 0.2, 0.2, 0.3, 0.1]
    @test_logs((:info, r"Trying to coerce from `Union{Missing,"),
                   global cv = coerce(v, OrderedFactor))
    cw = coerce(w, OrderedFactor)
    @test all(skipmissing(unique(cv)) .== [0.1, 0.2, 0.3])
    @test all(unique(cw) .== [0.1, 0.2, 0.3])
end

@testset "Any->MC" begin
    v2 = Symbol.(collect("asdfghj"))
    global v2c = @test_logs((:warn, r"Converting array"),
                            coerce(v2, Multiclass))
    @test scitype_union(v2c) == Multiclass{7}
    @test eltype(v2c) <: Union{Missing,CategoricalValue{String}}

    v3 = Any[1, 2, 3]
    coerce(v3, Multiclass) # should be no warning here

    # normal behaviour
    v1 = categorical([1,2,1,2,1,2,missing])
    v2 = collect("aksldjfalsdjkfslkjdfalksjdf")
    @test_logs((:info, r"Trying to coerce from `Union"),
               global v1c = coerce(v1, Multiclass))
    global v2c # otherwise julia complains...
    v2c = coerce(v2, Multiclass)
    @test scitype_union(v1c) == Union{Missing,Multiclass{2}}
    @test scitype_union(v2c) == Multiclass{7}
    @test eltype(v1c) <: Union{Missing,CategoricalValue{Int64}}
    @test eltype(v2c) <: CategoricalValue{Char}
end

@testset "Cat->C" begin
    a = categorical(["a","b","a","b",missing])
    a1 = coerce(a, Union{Count,Missing})
    @test scitype_union(a1) == Union{Missing,Count}
    @test all(skipmissing(a1 .== [1, 2, 1, 2, missing]))
    a1 = coerce(a, Union{Continuous,Missing})
    @test scitype_union(a1) == Union{Missing,Continuous}
    @test all(skipmissing(a1 .== [1., 2., 1., 2., missing]))

    y = categorical(1:10, ordered=true)
    new_order = [4, 10, 9, 7, 6, 2, 8, 3, 1, 5]
    levels!(y, new_order)
    @test all(coerce(y, Count) .== sortperm(new_order))
    @test all(coerce(y, Count) .== [9, 6, 8, 1, 10, 5, 4, 7, 3, 2])

    y = categorical([1:10..., missing, 11], ordered=true)
    @test all(skipmissing(coerce(y, Union{Continuous, Missing}) .==
                          float([1:10...,missing,11])))
end

# issue #62 (ScientficTypes)
@testset "Type=>Type" begin
    X = (x=[1,2,1,2,5,1,0,7],
         y=[0,1,0,1,0,1,0,1],
         z=['a','b','a','b','a','a',missing,missing])
    Xc = coerce(X, :y=>OrderedFactor)
    Xc = coerce(Xc, Count=>Continuous)
    @test elscitype(Xc.x) == Continuous
    @test elscitype(Xc.y) == OrderedFactor{2}
    Xc = coerce(Xc, OrderedFactor=>Count)
    @test elscitype(Xc.y) == Count
    Xc = coerce(Xc, :z=>Multiclass, verbosity=0)
    Xc = coerce(Xc, Multiclass=>OrderedFactor, verbosity=0)
    @test elscitype(Xc.z) == Union{Missing,OrderedFactor{2}}
    Xc = coerce(X, Count=>Continuous, Unknown=>Multiclass, verbosity=0)
    @test elscitype(Xc.x) == Continuous
    @test elscitype(Xc.y) == Continuous
    @test elscitype(Xc.z) == Union{Missing,Multiclass{2}}
end

# issue #13

@testset "mixture of Type=>Type and Symbol=>Type" begin
    X = (x=10:10:44, y=1:4, z=collect("abcd"), w=["a", "b", "c", missing])
    @test MLJScientificTypes.feature_scitype_pairs(:x => Continuous, X) ==
        [:x => Continuous, ]
    @test MLJScientificTypes.feature_scitype_pairs(Count => Continuous, X) ==
        [:x => Union{Continuous}, :y=> Union{Continuous}]
    X1 = coerce(X, :z => Multiclass, Count=>Continuous)
    @test schema(X1).scitypes ==
        (Continuous, Continuous, Multiclass{4}, Union{Missing,Textual})
    X2 = coerce(X, :z => Multiclass, :w=>Multiclass, Count=>Continuous,
                verbosity=0)
    @test schema(X2).scitypes ==
        (Continuous, Continuous, Multiclass{4}, Union{Missing,Multiclass{3}})
    @test_throws ArgumentError coerce(X, Count => Continuous, :x=>Multiclass)
end

@testset "Image" begin
    #1. Single GrayImage:
    img = rand(10, 10)
    @test coerce(img, GrayImage) |> size == size(img)
    @test scitype(coerce(img, GrayImage)) == GrayImage{10, 10}

    #2. Single colorimage
    img = rand(10, 10, 3)
    @test coerce(img, ColorImage) |> size == (10, 10)
    @test scitype(coerce(img, ColorImage)) == ColorImage{10, 10}

    #3. Collection: 3dcollection -> GrayImage
    imgs = rand(10, 10, 3)
    @test coerce(imgs, GrayImage) |> size == (3,)
    @test scitype(coerce(imgs, GrayImage)) == AbstractArray{GrayImage{10, 10},1}

    #4. Collection: 4dcollection -> GrayImage
    imgs = rand(10, 10, 1, 3)
    @test coerce(imgs, GrayImage) |> size == (3,)
    @test scitype(coerce(imgs, GrayImage)) == AbstractArray{GrayImage{10, 10},1}

    #5. Collection : 4dcollection -> ColorImage
    imgs = rand(10, 10, 3, 5)
    @test coerce(imgs, ColorImage) |> size == (5,)
    @test scitype(coerce(imgs, ColorImage)) == AbstractArray{ColorImage{10,10},1}
end

ext: jl


text:
struct Marginalization{D} <: AbstractVector{D}
    samples::DensitySampleVector
    vsel::Any
end

struct MarginalDist
    dist::ReshapedDist
end

function _get_edges(data::Tuple, nbins::Tuple{Vararg{<:Integer}}, closed::Symbol)
    return StatsBase.histrange(data, StatsBase._nbins_tuple(data, nbins), closed)
end

function _get_edges(data::Any, nbins::Integer, closed::Symbol)
    return StatsBase.histrange((data,), StatsBase._nbins_tuple((data,), (nbins,)), closed)[1]
end

function _get_edges(data::Any, nbins::Union{AbstractRange, Tuple{AbstractRange}}, closed::Symbol)
    return nbins
end

function MarginalDist(
    samples::Union{DensitySampleVector, StructArrays.StructVector},
    vsel;
    bins = 200, #::Union{I, Vector{I}, Tuple{I}} where I<:Integer
    closed::Symbol = :left,
    filter::Bool = false
)

    marg_samples = bat_marginalize(samples, vsel)
    vs = varshape(marg_samples)

    if filter
        marg_samples = BAT.drop_low_weight_samples(marg_samples)
    end

    marg_samples = flatview(unshaped.(marg_samples).v)
    cols = Tuple(Vector.(eachrow(marg_samples)))

    edges = if isa(bins, Integer)
        _get_edges(cols, (bins,), closed)
    else
        Tuple(_get_edges(cols[i], bins[i], closed) for i in 1:length(bins))
    end

    hist = fit(Histogram, cols, edges, closed = closed)

    binned_dist = EmpiricalDistributions.MvBinnedDist(hist)
    binned_dist = vs(binned_dist) isa ReshapedDist ? vs(binned_dist) : ReshapedDist(vs(binned_dist), vs)

    return MarginalDist(binned_dist)
end

#for prior or MarginalDist
#TODO: think about implementing parsing of already marginalized names for remarginalization
function MarginalDist(
    dist::Union{NamedTupleDist, MarginalDist},
    vsel;
    bins = 200, #::Union{I, Vector{I}, Tuple{I}} where I<:Integer
    closed::Symbol = :left,
    nsamples::Integer = 10^6,
    filter::Bool = false
)
    dist = dist isa MarginalDist ? dist.dist : dist

    vs = varshape(dist)
    samples = vs.(DensitySampleVector(unshaped.(rand(dist, nsamples)), fill(NaN, nsamples)))

    return MarginalDist(samples, vsel; bins, closed, filter)
end


text:
using BenchmarkTools
using DataFrames
using PrettyTables
using CUDA

using Oceananigans
using Oceananigans.Architectures
using JULES

using BenchmarkTools: prettytime, prettymemory
using JULES: intermediate_thermodynamic_field, compute_temperature!

Archs = [CPU]
@hascuda Archs = [CPU, GPU]

Ns = [32, 64]
@hascuda Ns = [32, 192]

Tvars = [Energy, Entropy]
Gases = [DryEarth, DryEarth3]

tags = ["arch", "N", "gases", "tvar"]

suite = BenchmarkGroup(
    "temperature" => BenchmarkGroup(),
    "pressure" => BenchmarkGroup()
)

_compute_temperature!(model) = compute_temperature!(model)
_compute_temperature!(model::CompressibleModel{GPU}) = CUDA.@sync compute_temperature!(model)

for Arch in Archs, N in Ns, Gases in Gases, Tvar in Tvars
    @info "Running temperature computation benchmark [$Arch, N=$N, $Tvar, $Gases]..."

    grid = RegularCartesianGrid(size=(N, N, N), extent=(1, 1, 1))
    model = CompressibleModel(architecture=Arch(), grid=grid, thermodynamic_variable=Tvar(),
                              gases=Gases())

    _compute_temperature!(model) # warmup

    b = @benchmark _compute_temperature!($model) samples=10

    key = (Arch, N, Gases, Tvar)
    suite["temperature"][key] = b
end

function benchmarks_to_dataframe(suite)
    df = DataFrame(architecture=[], size=[], gases=[],
                   thermodynamic_variable=[], min=[], median=[],
                   mean=[], max=[], memory=[], allocs=[])

    for (key, b) in suite
        Arch, N, Gases, Tvar = key

        d = Dict(
            "architecture" => Arch,
            "size" => "$(N)³",
            "gases" => Gases,
            "thermodynamic_variable" => Tvar,
            "min" => minimum(b.times) |> prettytime,
            "median" => median(b.times) |> prettytime,
            "mean" => mean(b.times) |> prettytime,
            "max" => maximum(b.times) |> prettytime,
            "memory" => prettymemory(b.memory),
            "allocs" => b.allocs
        )

        push!(df, d)
    end

    return df
end

header = ["Arch" "Size" "Gases" "ThermoVar" "min" "median" "mean" "max" "memory" "allocs"]

df = benchmarks_to_dataframe(suite["temperature"])
pretty_table(df, header, title="Temperature computation benchmarks", nosubheader=true)

function gpu_speedups(suite)
    df = DataFrame(size=[], gases=[], thermodynamic_variable=[], speedup=[])

    for N in Ns, Gas in Gases, Tvar in Tvars
        b_cpu = suite[CPU, N, Gas, Tvar] |> median
        b_gpu = suite[GPU, N, Gas, Tvar] |> median
        b_ratio = ratio(b_cpu, b_gpu)

        d = Dict(
            "size" => "$(N)³",
            "gases" => Gas,
            "thermodynamic_variable" => Tvar,
            "speedup" => @sprintf("%.3fx", b_ratio.time)
        )

        push!(df, d)
    end

    return df
end

header = ["Size" "Gases" "ThermoVar" "speedup"]

df = gpu_speedups(suite["temperature"])
pretty_table(df, header, title="Temperature computation speedups", nosubheader=true)


text:
using Clockwork
using Base.Test

@test 🕟+🕝== 🕖
@test 🕟-🕝== 🕐
@test 🕑^50 == 🕓
@test Clock("3:30") == 🕞
let cp = Dates.Hour(3) + Dates.Minute(30)
  @test Clock(cp) == 🕞
  @test Dates.CompoundPeriod(🕞) == cp
end
@test Dates.CompoundPeriod(🕒) == Dates.Hour(3)


text:
include("struct.jl")

_NORMAL_SELECT = true

"""
Return the number of cells in a group.
"""
function n_group_cells(cells::Array{Cell}, g::Int)
    return size(findall(x -> x.group == g, cells), 1)
end

"""
Return all coordinates that are stars
"""
function get_stars(cells)
    return findall(x -> x.is_star, cells)
end

function fill_cells!(cells, puzzle)
	d, n = size(puzzle)
	for i in 1:d
		for j in 1:n
			cell = Cell(i, j, puzzle[i, j]+1, false, false)
			cells[i, j] = cell
		end
	end
end

function fill_regions!(regions::Dict, cells::Array{Cell})
	d, n = size(cells)

	for i in 1:d
		regions["groups"][i] = Region(n_group_cells(cells, i), 0)
		regions["rows"][i] = Region(n, 0)
		regions["cols"][i] = Region(d, 0)
	end
end

function update_regions!(regions::Dict, cells::Array{Cell, 2})
	fill_regions!(regions, cells)
	d, n = size(cells)

	for cell in cells
		if cell.is_blocked || cell.is_star
			regions["rows"][cell.row].candidates -= 1
			regions["cols"][cell.col].candidates -= 1
			regions["groups"][cell.group].candidates -= 1
			if cell.is_star
				regions["rows"][cell.row].stars += 1
				regions["cols"][cell.col].stars += 1
				regions["groups"][cell.group].stars += 1
			end
		end
	end
end

function illegal(regions, k)
	to_check = [regions["rows"] regions["cols"] regions["groups"]]
	for region in to_check
		if (region.candidates < (k - region.stars))
			return true
		end
	end
	return false
end

function block!(cells, regions, k::Int)
	d, n = size(cells)
	stars = get_stars(cells)
	blocked = []

	for star in stars
		b = cells[star]

		check_row = regions["rows"][b.row].stars ≥ k
		check_col = regions["cols"][b.col].stars ≥ k
		check_groups = regions["groups"][b.group].stars ≥ k

		# adjacent
		check_points = findall(
			a -> !a.is_star
				&& ((check_row && a.row == b.row)			  	     # row
					|| (check_col && a.col == b.col)     		 	 # col
					|| (check_groups && a.group == b.group) 	 	 # group
					|| sqrt((a.row - b.row)^2+(a.col - b.col)^2) < 2 # adjacent
			), cells)
		blocked = cat(blocked, check_points, dims=1)
	end

	for coor in unique(blocked)
		cells[coor].is_blocked = true
	end

	update_regions!(regions, cells)
end

function smallest_region(regions)
	d = size(regions["rows"], 1)

	region_type = nothing
	min_region = nothing
	min_size = d*d

	for i in 1:d
		row_candidates = regions["rows"][i].candidates
		if 0 < row_candidates < min_size
			region_type = "rows"
			min_region = i
			min_size = row_candidates
		end

		col_candidates = regions["cols"][i].candidates
		if 0 < col_candidates < min_size
			region_type = "cols"
			min_region = i
			min_size = col_candidates
		end

		group_candidates = regions["groups"][i].candidates
		if 0 < group_candidates < min_size
			region_type = "groups"
			min_region = i
			min_size = group_candidates
		end
	end

	return min_region, region_type
end

function place_star!(cells::Array{Cell}, regions, normal_select)
	d, n = size(cells)

	idx = nothing

	if normal_select
		idx = findfirst(x -> !(x.is_blocked || x.is_star), cells)
	else
		region, type = smallest_region(regions)
		if type == "rows"
			idx = findfirst(x -> !(x.is_blocked || x.is_star)
				&& x.row == region, cells)
		elseif type == "cols"
			idx = findfirst(x -> !(x.is_blocked || x.is_star)
				&& x.col == region, cells)
		elseif type == "groups"
			idx = findfirst(x -> !(x.is_blocked || x.is_star)
				&& x.group == region, cells)
		end
	end

	if !isnothing(idx)
		cells[idx].is_star = true
		update_regions!(regions, cells)
		return cells[idx]
	end

	# throw(ArgumentError("could not find an available cell"))
	return
end

struct StarBattleResult
	cells::Union{Nothing, Array{Cell, 2}}
	steps::Int
end

function solve(puzzle::Array{Int, 2}, k::Int; normal_select::Bool=_NORMAL_SELECT, logging::Bool=false)
	d, n = size(puzzle)

	if logging
		println("$k-star $n x $n puzzle")
	end

	regions = Dict{String,Array}(
		"rows" => Vector{Region}(undef, d),
		"cols" => Vector{Region}(undef, d),
		"groups" => Vector{Region}(undef, d),
	)
	cells = Array{Cell}(undef, d, d)
	fill_cells!(cells, puzzle)
	fill_regions!(regions, cells)

	n = Node(data=deepcopy(cells))
	steps = 0
	while true
		steps += 1
		if size(get_stars(cells), 1) == k*d
			return StarBattleResult(cells, steps)
			break
		elseif !illegal(regions, k)
			cell = place_star!(cells, regions, normal_select)
			if !isnothing(cell)
				block!(cells, regions, k)
				insert!(n, cells, cell)
				n = n.children[end]
				cells = deepcopy(n.data)
			else
				return StarBattleResult(nothing, steps)
				break
			end
		elseif !isnothing(n.parent)
			bad_cell = n.cell
			n = n.parent
			delete_children!(n)
			n.data[bad_cell.row, bad_cell.col].is_blocked = true
			n.data[bad_cell.row, bad_cell.col].is_star = false
			cells = deepcopy(n.data)
			update_regions!(regions, cells)
		else
			return StarBattleResult(nothing, steps)
			break
		end
	end

	return n.cell
end


text:
using heuristicopt
using Test

@testset "heuristicopt.jl" begin
    @test "hoge" == "hoge"
    # Write your own tests here.
end

max_stars_repo_name: ariririri/heuristicopt

max_issues_repo_name: ariririri/heuristicopt


text:
using Bertweet
using Test

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


text:
# =============================================================
#		MODULE: plots
# =============================================================
module plots
	using PGFPlots

	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	#		FUNCTION : TERNARY
	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	function TERNARY(Path,  Sand, Silt, Clay)
		# Work in progress labels not sure
		Plot_SandSiltClay=TernaryAxis(Plots.Linear3(Sand,Silt,Clay,style="solid, blue, thick", onlyMarks=true),
		style="ternary limits relative=false",
		xlabel=L"Silt", ylabel=L"Clay", zlabel=L"Sand")

		save(Path, Plot_SandSiltClay)

		# label=L"Sand %", ylabel=L"Silt %", zlabel=L"Clay%"

	end  # function TERNARY

end  # module plots
# ............................................................

max_stars_repo_name: joseph-pollacco/TextureTriangleNZ

max_issues_repo_name: joseph-pollacco/TextureTriangleNZ

max_forks_repo_name: joseph-pollacco/TextureTriangleNZ


text:
using Documenter, REoptLite, JuMP

makedocs(
    sitename="REoptLite.jl Documentation",
    pages = [
        "Home" => "index.md",
        "REopt Lite" => Any[
            "reopt/examples.md",
            "reopt/inputs.md",
            "reopt/outputs.md",
            "reopt/methods.md"
        ],
        "Model Predictive Control" => Any[
            "mpc/examples.md",
            "mpc/inputs.md",
            "mpc/methods.md",
        ],
        "Developer" => Any[
            "developer/concept.md",
            "developer/organization.md",
            "developer/inputs.md",
            "developer/adding_tech.md"
        ]
    ],
    workdir = joinpath(@__DIR__, "..")
)

deploydocs(
    repo = "github.com/NREL/REoptLite.git",
)

max_stars_repo_name: NREL/REoptLite

max_issues_repo_name: NREL/REoptLite

max_forks_repo_name: NREL/REoptLite


text:
import StateSpaceReconstruction:
    Embeddings,
    assign_bin_labels,
    assign_coordinate_labels

import ..TransferOperators:
    BinVisits,
    TransferOperatorRectangularBinning,
    get_binvisits,
    estimate_transferoperator_from_binvisits

import DelayEmbeddings:
    Dataset

import StaticArrays:
    SVector, MVector

import CausalityToolsBase: get_minima_and_edgelengths, encode, RectangularBinning, RectangularBinningScheme, CustomReconstruction

export RectangularInvariantMeasure, rectangularinvariantmeasure

"""
    RectangularInvariantMeasure{T} <: AbstractRectangularInvariantMeasure where {T}

A RectangularInvariantMeasure created from a set of points representing a
state space by discretizing the state space into rectangular bins with
edge lengths dictated by the binning scheme `ϵ`. The invariant measure is
then computed from an approximation to the transfer operator over the
discretized state space.

## Fields
- **`points::AbstractArray{T, 2}`**: The points for which to estimate the invariant measure.
    Each column is a point.

- **`binning_scheme`**: The binning scheme. See docs for `RectangularBinning` for details.

- **`axisminima`**: The minima along each axis.

- **`edgelengths`**: `ϵ` converted to absolute edge lengths.

- **`visited_bins_inds::AbstractArray{Int,2}`**: Counting from the start of each
    coordinate axis in steps dictated by `ϵ`, which bins do each
    of the points lie in? One column vector of indices for each point in `points`
    (one index per coordinate axis).

- **`visited_bins_coordinates::AbstractArray{Float64,2}`**: The coordinates of the origin
    of each bin that is visited. Each column vector is a bin origin. The bins (rectangular
    boxes) have edge lengths given by `ϵ`.

- **`binvisits::BinVisits`**: A `BinVisits` instance, indicating which points
    visits which bin. The indices in this `BinVisits` instance refers to the unique columns
    of `visited_bins_inds` (one column for each point in `points`, of which some may
    be identical), or, equivalently,
    `visited_bins_coordinates` (also one column for each point in `points`).

- **`transfermatrix::TransferOperatorRectangularBinning`**: The transfer matrix from which
    the invariant measure is obtained.

- **`measure::InvariantDistribution`**: The invariant measure over the visited bins.
"""
struct RectangularInvariantMeasure{T} <: AbstractRectangularInvariantMeasure where {T}
    points
    binning_scheme::RectangularBinning
    axisminima::Vector{Float64}
    edgelengths::Vector{Float64}
    encoded_points::Vector{T}
    visited_bins_coordinates::Vector{Vector{Float64}}
    binvisits::BinVisits
    transfermatrix::TransferOperatorRectangularBinning
    measure::InvariantDistribution
end

"""
    rectangularinvariantmeasure(data,
        binning_scheme::RectangularBinning,
        estimator = TransferOperatorEstimatorRectangularBinning;
        kwargs...)

Estimate the invariant measure from a rectangular partition of `data`.

This is done by discretizing the state space into rectangular bins with
edge lengths dictated by the `binning_scheme`. We then approximate the
transfer operator over the discretized state space, and compute the
invariant measure over the bins from the transfer operator.

Returns a `RectangularInvariantMeasure` instance.

## Arguments
- **`data`**: The data from which to estimate the invariant measure. The following data
    types are currently accepted:
    - `Dataset` instances from `DynamicalSystems.jl`/`DelayEmbeddings.jl`.
    - `AbstractEmbedding` subtypes from `StateSpaceReconstruction.jl`.
    - `AbstractArray{T, 2}` instances where each column represents a point.
    - `Vector{Vector{T}}`
    - `Vector{SVector{D, T}}`
    - `Vector{MVector{D, T}}`

- **`binning_scheme`**: A valid `RectangularBinning` instance. See docs for
    `RectangularBinning` for details.

- **`estimator`**: A transfer operator estimator yielding a
    `TransferOperatorRectangularBinning`. Defaults to
    `:TransferOperatorEstimatorRectangularBinning`.

- **`kwargs`**: Keyword arguments when calling `invariantmeasure` on the
    transfer operator.
"""
function rectangularinvariantmeasure end

function rectangularinvariantmeasure(data::Vector{T},
        binning_scheme::RectangularBinning,
        estimator::Symbol = :TransferOperatorEstimatorRectangularBinning;
        kwargs...) where {T <: Union{Vector, SVector, MVector}}

    if estimator == :TransferOperatorEstimatorRectangularBinning
        # Identify which bins of the partition resulting from using ϵ each
        # point of the embedding visits.

        # Get absolute bin sizes
        mini, edgelengths = get_minima_and_edgelengths(data, binning_scheme)
        encoded_pts = encode(data, mini, edgelengths)

        # The coordinate of the bin origins
        visited_bins_coordinates = [edgelengths .* pt .+ mini for pt in encoded_pts] #for pt #assign_coordinate_labels(data, ϵ)

        # Which are the visited bins, which points
        # visits which bin, repetitions, etc...
        binvisits = get_binvisits(encoded_pts)

        # Use that information to estimate transfer operator
        TO = estimate_transferoperator_from_binvisits(binvisits)

        # Compute invariant measure
        ivm = invariantmeasure(TO; kwargs...)

        RectangularInvariantMeasure(
            data,
            (binning_scheme isa RectangularBinning) ? binning_scheme : RectangularBinning(binning_scheme),
            mini,
            edgelengths,
            encoded_pts,
            visited_bins_coordinates,
            binvisits,
            TO,
            ivm
        )
    end
end

function rectangularinvariantmeasure(data::Dataset,
    binning_scheme::RectangularBinning,
    estimator = :TransferOperatorEstimatorRectangularBinning;
    kwargs...)

    rectangularinvariantmeasure(data.data, binning_scheme, estimator, kwargs...)
end

function rectangularinvariantmeasure(data::CustomReconstruction,
    binning_scheme::RectangularBinning,
    estimator = :TransferOperatorEstimatorRectangularBinning;
    kwargs...)

    rectangularinvariantmeasure(data.reconstructed_pts, binning_scheme, estimator, kwargs...)
end

function rectangularinvariantmeasure(data::AbstractArray{T, 2},
    binning_scheme::RectangularBinning,
    estimator::Symbol = :TransferOperatorEstimatorRectangularBinning;
    kwargs...) where {T <: Real}

    if size(data, 1) > size(data, 2)
        rectangularinvariantmeasure(Dataset(data), binning_scheme, estimator; kwargs...)
    else
        rectangularinvariantmeasure(Dataset(transpose(data)), binning_scheme, estimator;
            kwargs...)
    end
end

function summarise(invm::RectangularInvariantMeasure)
    if invm.points isa AbstractArray{T, 2} where T
        if size(invm.points, 1) > size(invm.points, 2)
            npoints = size(invm.points, 1)
        else
            npoints = size(invm.points, 2)
        end
    else
        npoints = length(invm.points)
    end

    D = length(invm.encoded_points[1])
    unique_states_visited = length(unique(invm.encoded_points))
    ϵ = invm.binning_scheme

    return join([typeof(invm), " from $npoints $D-dimensional points visiting $unique_states_visited unique states in the partition formed by the binning scheme $ϵ"])
end

Base.show(io::IO, invm::RectangularInvariantMeasure) = println(io, summarise(invm))


text:
#=
#
# Reads a mesh file and builds a MeshData struct
=#
export read_mesh

# numbers of CORNER nodes for first and second order elements as defined by GMSH
# TODO add higher order types
etypetonv = [2, 3, 4, 4, 8, 6, 5, 2, 3, 4, 4, 8, 6, 5, 1, 4, 8, 6, 5]; # number of vertices for each type
etypetodim= [1, 2, 2, 3, 3, 3, 3, 1, 2, 2, 3, 3, 3, 3, 1, 2, 2, 3, 3]; # dimension of each type

# Reads from the file stream
# Returns a MeshData struct
function read_mesh(file)
    # Determine the MSH format version
    mesh_file_type = 0; # 1=msh_v2, 2=msh_v4, 3=medit
    while !eof(file)
        line = readline(file);
        if occursin("\$MeshFormat", line)
            # Check if the MSH version is old(2) or new(4)
            vals = split(readline(file), " ", keepempty=false);
            if parse(Float64, vals[1]) >= 4
                mesh_file_type = 2;
            else
                mesh_file_type = 1;
            end
            break;
        elseif occursin("MeshVersionFormatted", line)
            mesh_file_type = 3;
        end
    end

    # use the appropriate reader
    seekstart(file);
    if mesh_file_type == 1
        return read_msh_v2(file);
    elseif mesh_file_type == 2
        return read_msh_v4(file);
    elseif mesh_file_type == 3
        return read_medit(file);
    else
        printerr("Couldn't recognize mesh file type. Use GMSH .msh version 2 or 4, or Medit .mesh")
    end
end

function read_msh_v2(file)
    # Find the beginning of the Nodes and Elements sections.
    # Then read in the numbers.
    nodes_done = false;
    elements_done = false;
    nx = 0;
    nel = 0;
    nodes = [];
    indices = [];
    elements = [];
    etypes = [];
    nv = [];
    while((!nodes_done || !elements_done) && !eof(file))
        line = readline(file);
        if occursin("\$Nodes", line)
            # The Nodes section
            line = readline(file);
            nx = parse(Int, split(line, " ", keepempty=false)[1]);
            if nx > 0
                # parse each node's coordinates
                nodes = zeros(3, nx);
                indices = zeros(Int, nx);
                i = 1;
                line = readline(file);
                while !occursin("\$EndNodes", line) && !eof(file)
                    vals = split(line, " ", keepempty=false);
                    indices[i] = parse(Int, vals[1]);
                    nodes[1,i] = parse(Float64, vals[2]);
                    nodes[2,i] = parse(Float64, vals[3]);
                    nodes[3,i] = parse(Float64, vals[4]);
                    i += 1;
                    line = readline(file);
                end

                nodes_done = true;
            end
        elseif occursin("\$Elements", line)
            # The Elements section
            line = readline(file);
            nel = parse(Int, split(line, " ", keepempty=false)[1]);
            if nel > 0
                # parse each element's numbers
                nv = zeros(Int, nel);
                etypes = zeros(Int, nel);
                elements = zeros(Int, 8, nel);
                i = 1;
                line = readline(file);
                while !occursin("\$EndElements", line) && !eof(file)
                    vals = split(line, " ", keepempty=false);
                    # Skip lower dimensional elements
                    tmptype = parse(Int, vals[2]);
                    if etypetodim[tmptype] == config.dimension
                        etypes[i] = tmptype;
                        nv[i] = etypetonv[tmptype];
                        offset = parse(Int, vals[3]) + 3;
                        for j=1:nv[i]
                            elements[j,i] = parse(Int, vals[offset + j]);
                        end
                        i += 1;
                    end

                    line = readline(file);
                end
                # adjust to correct number of elements
                nel = i-1;
                nv = nv[1:nel];
                etypes = etypes[1:nel];
                elements = elements[:, 1:nel];

                elements_done = true;
            end
        end
    end
    nodes = nodes[1:config.dimension, :];

    return MeshData(nx, nodes, indices, nel, elements, etypes, nv);
end

function read_msh_v4(file)
    # Find the beginning of the Nodes and Elements sections.
    # Then read in the numbers.
    nodes_done = false;
    elements_done = false;
    nx = 0;
    nel = 0;
    nodes = [];
    indices = [];
    elements = [];
    etypes = [];
    nv = [];
    while((!nodes_done || !elements_done) && !eof(file))
        line = readline(file);
        if occursin("\$Nodes", line)
            # The Nodes section
            #=
            $Nodes
            numEntityBlocks(size_t) numNodes(size_t) minNodeTag(size_t) maxNodeTag(size_t)
            entityDim(int) entityTag(int) parametric(int; 0 or 1) numNodesInBlock(size_t)
                nodeTag(size_t)
                ...
                x(double) y(double) z(double)
                < u(double; if parametric and entityDim >= 1) >
                < v(double; if parametric and entityDim >= 2) >
                < w(double; if parametric and entityDim == 3) >
                ...
            ...
            $EndNodes
            =#
            line = readline(file);
            nx = parse(Int, split(line, " ", keepempty=false)[2]);
            if nx > 0
                # parse node info
                nodes = zeros(3,nx);
                indices = zeros(Int, nx);
                i = 1;
                line = readline(file);
                while !occursin("\$EndNodes", line) && !eof(file)
                    # process entity blocks
                    # The first line has entity info.
                    entinfo = split(line, " ", keepempty=false);
                    entdim = parse(Int, entinfo[1]);
                    parametric = parse(Int, entinfo[3]);
                    entnx = parse(Int, entinfo[4]);
                    # Node indices
                    for ni=1:entnx
                        line = readline(file);
                        vals = split(line, " ", keepempty=false);
                        indices[i-1 + ni] = parse(Int, vals[1]);
                    end
                    # Node coordinates
                    for ni=1:entnx
                        line = readline(file);
                        vals = split(line, " ", keepempty=false);
                        nodes[1,i-1 + ni] = parse(Float64, vals[1]);
                        nodes[2,i-1 + ni] = parse(Float64, vals[2]);
                        nodes[3,i-1 + ni] = parse(Float64, vals[3]);
                    end
                    # Parametric (ignore)
                    if parametric > 0
                        printerr("Paremetric entity found in mesh file. Ignoring parameters.")
                        for ni=1:entdim
                            line = readline(file);
                        end
                    end
                    i += entnx;
                    line = readline(file);
                end

                nodes_done = true;
            end
        elseif occursin("\$Elements", line)
            # The Elements section
            #=
            $Elements
            1 2 1 2          1 entity bloc, 2 elements total, min/max element tags: 1 and 2
            2 1 3 2          2D entity (surface) 1, element type 3 (4-node quad), 2 elements
            1 1 2 3 4          quad tag #1, nodes 1 2 3 4
            2 2 5 6 3          quad tag #2, nodes 2 5 6 3
            $EndElements
            =#
            line = readline(file);
            nel = parse(Int, split(line, " ", keepempty=false)[2]);
            if nel > 0
                # parse elements
                nv = zeros(Int, nel);
                etypes = zeros(Int, nel);
                elements = zeros(Int, 8, nel);
                i = 1;
                line = readline(file);
                while !occursin("\$EndElements", line) && !eof(file)
                    # Process entity blocks
                    # The first line has entity info.
                    vals = split(line, " ", keepempty=false);
                    entdim = parse(Int, vals[3]);
                    enttype = parse(Int, vals[3]);
                    entnel = parse(Int, vals[4]);

                    # Ignore elements with lower dimension (faces, points etc.)
                    if entdim == config.dimension
                        # read element info
                        for ei=1:entnel
                            line = readline(file);
                            vals = split(line, " ", keepempty=false);
                            etypes[i-1 + ei] = enttype;
                            entnv = etypetonv[enttype];
                            nv[i-1 + ei] = entnv;
                            for j=1:entnv
                                elements[j, i-1 + ei] = parse(Int, vals[1 + j]);
                            end
                        end

                        i += entnel;
                        line = readline(file);
                    else
                        for ei=1:entnel
                            line = readline(file); # skip lines
                        end
                        line = readline(file);
                    end

                end

                # adjust nel to remove unwanted elements
                nel = i-1;
                elements = elements[:,1:nel];
                etypes = etypes[1:nel];
                nv = nv[1:nel];

                elements_done = true;
            end
        end
    end
    nodes = nodes[1:config.dimension, :];

    return MeshData(nx, nodes, indices, nel, elements, etypes, nv);
end

function read_medit(file)
    # Whitespace is inconsistent, but essentially it will look like
    #
    # MeshVersionFormatted 1
    # Dimension 3
    # Vertices 12345
    # 0.0509745 -0.268407 -0.198911 0      <-- Note: always 4 values (4th is a tag?)
    # ...
    # Tetrahedra 1234
    # 1 2 3 4 0     <-- again an extra tag value
    # ...
    # End
    #
    # There may be lower dimensional elements for surfaces, but it's inconsistent
    nodes_done = false;
    elements_done = false;
    nx = 0;
    nel = 0;
    nodes = [];
    indices = [];
    elements = [];
    etypes = [];
    nv = [];
    while((!nodes_done || !elements_done) && !eof(file))
        line = readline(file);
        if occursin("Dimension", line)
            tokens = split(line, " ", keepempty=false);
            if length(tokens) > 1
                dim = parse(Int, tokens[2]);
            else
                line = readline(file);
                dim = parse(Int, line);
            end
            if config.dimension != dim
                # Maybe someone forgot to set dimension?
                printerr("Dimension in mesh file doesn't match config. Updating config.")
                config.dimension = dim;
            end
        elseif occursin("Vertices", line)
            println(line)
            tokens = split(line, " ", keepempty=false);
            if length(tokens) > 1
                nx = parse(Int, tokens[2]);
            else
                line = readline(file);
                nx = parse(Int, line);
            end

            # parse vertex info
            nodes = zeros(3,nx);
            flags = zeros(Int, nx);
            for i=1:nx
                line = readline(file);
                vals = split(line, " ", keepempty=false);
                nodes[1,i] = parse(Float64, vals[1]);
                nodes[2,i] = parse(Float64, vals[2]);
                nodes[3,i] = parse(Float64, vals[3]);
                flags[i] = parse(Float64, vals[4]);
            end

            nodes_done = true;

        elseif config.dimension == 1 && (occursin("Lines", line))
            # Is this even an option for medit?
        elseif config.dimension == 2 && (occursin("Triangles", line) || occursin("Quadrilaterals", line))
            tokens = split(line, " ", keepempty=false);
            if length(tokens) > 1
                nel = parse(Int, tokens[2]);
            else
                line = readline(file);
                nel = parse(Int, line);
            end
            if tokens[1] == "Triangles"
                nv = fill(3, nel);
                etypes = fill(2, nel);
                elements = zeros(Int, 8, nel);
                eflags = zeros(Int, nel);
                for i=1:nel
                    line = readline(file);
                    vals = split(line, " ", keepempty=false);
                    elements[1,i] = parse(Float64, vals[1]);
                    elements[2,i] = parse(Float64, vals[2]);
                    elements[3,i] = parse(Float64, vals[3]);
                    eflags[i] = parse(Float64, vals[4]);
                end

            elseif tokens[1] == "Quadrilaterals"
                nv = fill(4, nel);
                etypes = fill(3, nel);
                elements = zeros(Int, 8, nel);
                eflags = zeros(Int, nel);
                for i=1:nel
                    line = readline(file);
                    vals = split(line, " ", keepempty=false);
                    elements[1,i] = parse(Float64, vals[1]);
                    elements[2,i] = parse(Float64, vals[2]);
                    elements[3,i] = parse(Float64, vals[3]);
                    elements[4,i] = parse(Float64, vals[4]);
                    eflags[i] = parse(Float64, vals[5]);
                end
            end
            elements_done = true;

        elseif config.dimension == 3 && (occursin("Tetrahedra", line) || occursin("Hexahedra", line))
            println(line)
            tokens = split(line, " ", keepempty=false);
            if length(tokens) > 1
                nel = parse(Int, tokens[2]);
            else
                line = readline(file);
                nel = parse(Int, line);
            end

            if tokens[1] == "Tetrahedra"
                nv = fill(4, nel);
                etypes = fill(4, nel);
                elements = zeros(Int, 8, nel);
                eflags = zeros(Int, nel);
                for i=1:nel
                    line = readline(file);
                    vals = split(line, " ", keepempty=false);
                    elements[1,i] = parse(Float64, vals[1]);
                    elements[2,i] = parse(Float64, vals[2]);
                    elements[3,i] = parse(Float64, vals[3]);
                    elements[4,i] = parse(Float64, vals[4]);
                    eflags[i] = parse(Float64, vals[5]);
                end

            elseif tokens[1] == "Hexahedra"
                nv = fill(8, nel);
                etypes = fill(5, nel);
                elements = zeros(Int, 8, nel);
                eflags = zeros(Int, nel);
                for i=1:nel
                    line = readline(file);
                    vals = split(line, " ", keepempty=false);
                    elements[1,i] = parse(Float64, vals[1]);
                    elements[2,i] = parse(Float64, vals[2]);
                    elements[3,i] = parse(Float64, vals[3]);
                    elements[4,i] = parse(Float64, vals[4]);
                    elements[5,i] = parse(Float64, vals[5]);
                    elements[6,i] = parse(Float64, vals[6]);
                    elements[7,i] = parse(Float64, vals[7]);
                    elements[8,i] = parse(Float64, vals[8]);
                    eflags[i] = parse(Float64, vals[9]);
                end
            end
            elements_done = true;
        end
    end
    if !nodes_done || !elements_done
        printerr("Did not successfully read the mesh file.")
        return nothing;
    end
    nodes = nodes[1:config.dimension, :];
    indices = 1:nx;

    return MeshData(nx, nodes, indices, nel, elements, etypes, nv);
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 00145d19-9a62-4343-89c8-c623688c5af6
using PlutoUI

# ╔═╡ a2cd8e7c-a730-44a5-a111-a18376eeddad
using HypertextLiteral

# ╔═╡ 585c2b10-793a-44f2-8c98-c102ec64bc9b
md"""
# Basic
"""

# ╔═╡ 99c57e5a-4f2f-4f57-bbe5-b9b435cd7619
md"""
This should say `123`:
"""

# ╔═╡ 538c90e0-76c9-4030-83fb-a2f76871b6e2
embed_display(123)

# ╔═╡ 24155278-a621-49d4-831c-e2adc78719d2
md"""
These two should look exactly the same:
"""

# ╔═╡ 640ae34f-86ec-4f21-8afb-202aa1b4c70d
begin
	demo_img_data = let
		url = "https://user-images.githubusercontent.com/6933510/116753174-fa40ab80-aa06-11eb-94d7-88f4171970b2.jpeg"
		data = read(download(url))
	end

	demo_img() = PlutoUI.Show(MIME"image/jpg"(), demo_img_data)

	demo_img()
end

# ╔═╡ 041886d2-c38f-4a9e-83da-c4dd5e0bce27
const zzez = md"""
The next array should appear

not here 👇, but here 👉 $(embed_display([1,2,3])) 👈 here

And the next image should appear

not here 👇, but here 👉 $(embed_display(demo_img())) 👈 here
"""

# ╔═╡ 62b8be72-3bf8-429a-ba22-52802dfe8639
embed_display(embed_display(demo_img()))

# ╔═╡ 986cf446-dca9-4468-ba0b-73351d02f991
md"""
These two should look exactly the same:
"""

# ╔═╡ 4f24288e-3f19-442b-a440-8427773e840d
html"
Hello
"

# ╔═╡ aa560488-d50d-4b90-80f3-40e691b5a65d
embed_display(html"
Hello
")

# ╔═╡ 4d17f528-1480-44c0-b492-8ea651831e57
md"""
These three should look exactly the same:
"""

# ╔═╡ 41a3fdf0-9b8b-4444-b3e3-f926f075ffba
["asdfasdf a" for x in 1:50]

# ╔═╡ ba5d08b3-7b20-4958-bf0d-838b7b400434
embed_display(["asdfasdf a" for x in 1:50])

# ╔═╡ 4e8fc9f8-671a-465f-ba7a-2d41b69f7c49
md"""
These four should look exactly the same:
"""

# ╔═╡ 0a94496e-23eb-4ba6-b748-f0eb28eb1534
begin
	long_array = rand(50)
end

# ╔═╡ 3464d183-5c85-48f4-b754-e0ef07533e70
embed_display(long_array)

# ╔═╡ 0dd0857a-82e5-4144-9eb8-df54470bfdcf
md"""
$(embed_display(long_array))
"""

# ╔═╡ 0d1fd955-a0dc-4230-9fd8-5db13d1c6d0b

# ╔═╡ d55644f6-9f3f-424a-89d4-2cf053d1988a
nested(n, x) = if n == 0
	x
else
	embed_display(nested(n - 1, x))
end

# ╔═╡ d7d87f41-5417-412f-aab8-03ad773c519d
nested(n) = x -> nested(n, x)

# ╔═╡ e93308d3-d1c9-4df1-b796-1bf6bc0ad0c2
nested(4, long_array)

# ╔═╡ ba397445-db5f-457b-bedb-f3dcf77e5b5f
md"""
These three should look exactly the same:
"""

# ╔═╡ 0652dfdc-cb6b-4e79-997f-6e64b8f0e7b2
Any[Any[1,2],Any[3,md"## 4"]]

# ╔═╡ 9df6da34-df5d-47ff-a799-af6ff3312b2e
to_any(x) = Any[x...]

# ╔═╡ 5070715c-9e3d-4de3-a535-733dd8f1f2bd
to_any(embed_display.(["asdfasdf a" for x in 1:50]))

# ╔═╡ 3984969b-7af9-44d9-bfe3-b620ca9f4a3d
to_any(embed_display.(long_array))

# ╔═╡ f0cd5783-f5a9-4ab6-9c22-6df8d17b7c43
to_any(embed_display.([[1,2],[3,md"## 4"]]))

# ╔═╡ 0e5fee32-0cfb-434f-8c6b-533150641c06
Any[to_any(embed_display.([1,2])),to_any(embed_display.([3,md"## 4"]))]

# ╔═╡ 5dbecdca-f91b-48fa-ab8f-7510afd2055e
md"""
# Running scripts
"""

# ╔═╡ 39ac42fb-6221-40b1-a811-0db673329d77
md"""
Both checkboxes should each trigger exactly one alert:
"""

# ╔═╡ baf99bdf-21ac-49d9-aac0-666a3e2f2927
@bind test_script_runs_1 html"<input type=checkbox>"

# ╔═╡ 1843f036-f003-4da8-8c91-4a8d25fab6ad
test_script_runs_1 === true && embed_display(@htl("""

		$(embed_display([1,2]))

		"""))

# ╔═╡ ce78e7b8-6d74-43b5-8559-1688b933666d
@bind test_script_runs_2 html"<input type=checkbox>"

# ╔═╡ 7921e510-f2cd-433e-8675-a1d3e698e3fe
test_script_runs_2 === true && nested(3,@htl("""

		$(embed_display([1,2]))

		"""))

# ╔═╡ 1a9d606d-54fc-4948-be9e-52cbdcc79627
md"""
# More items
"""

# ╔═╡ 90f2103d-7cc1-4a82-a056-f141519239af
md"""
Expand the three arrays below and click on "More". You should:
- load more items
- remain in the expanded view after more items are loaded
"""

# ╔═╡ 28e30787-950a-4044-a16b-651e6ce5c649
@bind reset_more html"<input type=button value=Reset>"

# ╔═╡ be19ac0a-0849-4dfe-a786-cf5cf25ff385
reset_more; embed_display(embed_display(embed_display(rand(50))))

# ╔═╡ 43d5a69b-9fdd-43cf-bfa6-0c2706c0c96a
reset_more; embed_display(rand(50))

# ╔═╡ 3a8fa3d0-02f3-4f65-9634-95cb1e533e3a
md"""
# Bonds
"""

# ╔═╡ 777d670e-cdb7-48f0-ba3e-59b0044c6880
md"""
These four sliders:
- should be synchronised
- should all control `test_b`
"""

# ╔═╡ 22ab869f-8126-4f2a-a05e-12fdded20625
begin
	test_bond = @bind test_b html"<input type=range>"
end

# ╔═╡ 293579ef-421c-428e-96a7-aa66f21de87b
test_b

# ╔═╡ 413e9286-94c5-4c08-a4a4-9ed9e6f50584
test_bond |> embed_display

# ╔═╡ 6999e86a-cd68-4e97-a875-38a9c1ac8939
md"""
# Persistance
"""

# ╔═╡ e579aafa-60ce-4175-9348-02e6b5687323
md"""
## Same objects, same displays
"""

# ╔═╡ 5636ff9c-b887-4c4c-a02e-b9fe03124a29
md"""
When you move the slider:

- There should be no flickering
- If the array is expanded, it remains expanded
- _(Failing, but that's fine)_ If the array is showing more items, it will keep those items
"""

# ╔═╡ 51ba2785-2931-4c82-9d1e-8ac162a765af
@bind pers_1_trigger html"<input type=range>"

# ╔═╡ 95f6e0d0-7a2f-4956-b907-9e0e358056f4
md"""
## _(Failing)_ Same objects, new displays
"""

# ╔═╡ c276ef71-913b-476c-b2ff-510b9d5b1d51
md"""
When you move the slider:

- _(Failing)_ There should be no flickering
- _(Failing)_ If the array is expanded, it remains expanded
- _(Failing, but that's fine)_ If the array is showing more items, it will keep those items
"""

# ╔═╡ 41b5dd76-0e7a-40c8-9e33-a9286739bf9a
@bind pers_2_trigger html"<input type=range>"

# ╔═╡ 7738ab88-0717-4007-830e-e2a50ff0aa1d
pers_2_items = (rand(50), demo_img());

# ╔═╡ e77670a0-0574-4bed-9a60-956f0436caa3
md"""
## _(Failing, but that's fine)_ New objects, new displays
"""

# ╔═╡ 3f684dfe-605c-4b80-a551-341eb6e816f7
md"""
When you move the slider:

- _(Failing)_ There should be no flickering
- _(Failing)_ If the array is expanded, it remains expanded
- _(Failing, but that's fine)_ If the array is showing more items, it will keep those items
"""

# ╔═╡ 70895e40-0192-4f82-b891-5fa43c699a18
@bind pers_3_trigger html"<input type=range>"

# ╔═╡ 89f19fd0-d190-40cb-ba3e-093ec8acbc6e
pers_3_items() = (rand(50), demo_img());

# ╔═╡ a6560e5e-05bb-4966-ae00-05532ba091de
md"""
# Stuff
"""

# ╔═╡ 0fe7980e-0942-4619-aee2-db8f0f7d8f19
showhtml(r) = repr(MIME"text/html"(), r) |> Text

# ╔═╡ b7944093-6eb1-41b3-8f6b-096bd83a9315
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ 75c5b6c7-2705-4151-a28b-92522db7b592
Layout = ingredients(download("https://fonsp-disorganised-mess.netlify.app/Layout.jl"))

# ╔═╡ abd620ef-fd30-4adb-8681-a26267d21f6b
equal_container(x) = Layout.Div(x, Dict(
			"flex" => "1 1 0px",
			"overflow-x" => "auto",
		))

# ╔═╡ 8139f5ef-1d87-406d-b26b-e0baac01d104
hbox_even(elements...) = Layout.flex(equal_container.(elements)...)

# ╔═╡ e418c6d5-0eb2-4b38-ae8d-74e115ce6fe8
reset_more; hbox_even(embed_display(rand(50)), embed_display(demo_img()))

# ╔═╡ e97c5574-e808-4e35-a1fe-ffba661132de
hbox_even(test_bond |> nested(4), test_bond |> nested(2))

# ╔═╡ 518d06da-c828-47bb-b3c2-92b36a95a52c
pers_1 = hbox_even(embed_display(rand(50)), embed_display(demo_img()));

# ╔═╡ ac9b8afe-8d54-41ce-8ed3-fefd6c9ba3ae
pers_1_trigger; pers_1

# ╔═╡ bc03bc0f-0d62-4661-b7d3-f7f44c364bff
pers_2() = hbox_even(embed_display.(pers_2_items)...);

# ╔═╡ fb98383e-5e92-487b-964d-986f6c12b70f
pers_2_trigger; pers_2()

# ╔═╡ 17020b64-04bc-4124-aa6a-eaf4210cf6d2
pers_3() = hbox_even(embed_display.(pers_3_items())...);

# ╔═╡ 1c72b7e1-5fb3-4f8b-86a4-c508348ffdf0
pers_3_trigger; pers_3()

# ╔═╡ Cell order:
# ╟─585c2b10-793a-44f2-8c98-c102ec64bc9b
# ╟─99c57e5a-4f2f-4f57-bbe5-b9b435cd7619
# ╠═538c90e0-76c9-4030-83fb-a2f76871b6e2
# ╟─041886d2-c38f-4a9e-83da-c4dd5e0bce27
# ╟─24155278-a621-49d4-831c-e2adc78719d2
# ╟─640ae34f-86ec-4f21-8afb-202aa1b4c70d
# ╟─62b8be72-3bf8-429a-ba22-52802dfe8639
# ╟─986cf446-dca9-4468-ba0b-73351d02f991
# ╟─4f24288e-3f19-442b-a440-8427773e840d
# ╟─aa560488-d50d-4b90-80f3-40e691b5a65d
# ╟─4d17f528-1480-44c0-b492-8ea651831e57
# ╟─41a3fdf0-9b8b-4444-b3e3-f926f075ffba
# ╟─ba5d08b3-7b20-4958-bf0d-838b7b400434
# ╟─5070715c-9e3d-4de3-a535-733dd8f1f2bd
# ╟─4e8fc9f8-671a-465f-ba7a-2d41b69f7c49
# ╟─0a94496e-23eb-4ba6-b748-f0eb28eb1534
# ╟─3464d183-5c85-48f4-b754-e0ef07533e70
# ╟─3984969b-7af9-44d9-bfe3-b620ca9f4a3d
# ╟─e93308d3-d1c9-4df1-b796-1bf6bc0ad0c2
# ╟─0dd0857a-82e5-4144-9eb8-df54470bfdcf
# ╟─0d1fd955-a0dc-4230-9fd8-5db13d1c6d0b
# ╟─d55644f6-9f3f-424a-89d4-2cf053d1988a
# ╟─d7d87f41-5417-412f-aab8-03ad773c519d
# ╟─ba397445-db5f-457b-bedb-f3dcf77e5b5f
# ╟─0652dfdc-cb6b-4e79-997f-6e64b8f0e7b2
# ╟─f0cd5783-f5a9-4ab6-9c22-6df8d17b7c43
# ╟─0e5fee32-0cfb-434f-8c6b-533150641c06
# ╟─9df6da34-df5d-47ff-a799-af6ff3312b2e
# ╟─5dbecdca-f91b-48fa-ab8f-7510afd2055e
# ╟─39ac42fb-6221-40b1-a811-0db673329d77
# ╟─baf99bdf-21ac-49d9-aac0-666a3e2f2927
# ╟─1843f036-f003-4da8-8c91-4a8d25fab6ad
# ╟─ce78e7b8-6d74-43b5-8559-1688b933666d
# ╟─7921e510-f2cd-433e-8675-a1d3e698e3fe
# ╟─1a9d606d-54fc-4948-be9e-52cbdcc79627
# ╟─90f2103d-7cc1-4a82-a056-f141519239af
# ╟─be19ac0a-0849-4dfe-a786-cf5cf25ff385
# ╟─43d5a69b-9fdd-43cf-bfa6-0c2706c0c96a
# ╟─e418c6d5-0eb2-4b38-ae8d-74e115ce6fe8
# ╟─28e30787-950a-4044-a16b-651e6ce5c649
# ╟─3a8fa3d0-02f3-4f65-9634-95cb1e533e3a
# ╟─777d670e-cdb7-48f0-ba3e-59b0044c6880
# ╠═293579ef-421c-428e-96a7-aa66f21de87b
# ╟─22ab869f-8126-4f2a-a05e-12fdded20625
# ╟─413e9286-94c5-4c08-a4a4-9ed9e6f50584
# ╟─e97c5574-e808-4e35-a1fe-ffba661132de
# ╟─6999e86a-cd68-4e97-a875-38a9c1ac8939
# ╟─e579aafa-60ce-4175-9348-02e6b5687323
# ╟─5636ff9c-b887-4c4c-a02e-b9fe03124a29
# ╟─51ba2785-2931-4c82-9d1e-8ac162a765af
# ╟─ac9b8afe-8d54-41ce-8ed3-fefd6c9ba3ae
# ╠═518d06da-c828-47bb-b3c2-92b36a95a52c
# ╟─95f6e0d0-7a2f-4956-b907-9e0e358056f4
# ╟─c276ef71-913b-476c-b2ff-510b9d5b1d51
# ╟─41b5dd76-0e7a-40c8-9e33-a9286739bf9a
# ╟─fb98383e-5e92-487b-964d-986f6c12b70f
# ╠═bc03bc0f-0d62-4661-b7d3-f7f44c364bff
# ╠═7738ab88-0717-4007-830e-e2a50ff0aa1d
# ╟─e77670a0-0574-4bed-9a60-956f0436caa3
# ╟─3f684dfe-605c-4b80-a551-341eb6e816f7
# ╟─70895e40-0192-4f82-b891-5fa43c699a18
# ╟─1c72b7e1-5fb3-4f8b-86a4-c508348ffdf0
# ╠═17020b64-04bc-4124-aa6a-eaf4210cf6d2
# ╠═89f19fd0-d190-40cb-ba3e-093ec8acbc6e
# ╟─a6560e5e-05bb-4966-ae00-05532ba091de
# ╠═0fe7980e-0942-4619-aee2-db8f0f7d8f19
# ╠═00145d19-9a62-4343-89c8-c623688c5af6
# ╠═a2cd8e7c-a730-44a5-a111-a18376eeddad
# ╠═abd620ef-fd30-4adb-8681-a26267d21f6b
# ╠═8139f5ef-1d87-406d-b26b-e0baac01d104
# ╠═75c5b6c7-2705-4151-a28b-92522db7b592
# ╟─b7944093-6eb1-41b3-8f6b-096bd83a9315


text:
using PFVG_AA2021a
using Test

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


text:
"""
    accola_chars_r(e, g1, g2, g3)

Compute all characteristics of the form e+g, where g is in the group generated by g1, g2, g3.
"""
function accola_chars_r(e, g1, g2, g3)
    return [e, remainder_char(e+g1), remainder_char(e+g2), remainder_char(e+g3), remainder_char(e+g1+g2), remainder_char(e+g1+g3), remainder_char(e+g2+g3), remainder_char(e+g1+g2+g3)];
end

"""
    accola_chars_sr(chars, indices)

Compute the characteristics used in Accola's special theta relations.
"""
function accola_chars_sr(chars, indices)
    g1 = remainder_char(sum([chars[i] for i in indices[2]]));
    g2 = remainder_char(sum([chars[i] for i in indices[3]]));
    g3 = remainder_char(sum([chars[i] for i in indices[4]]));
    return [accola_chars_r(chars[1], g1, g2, g3), accola_chars_r(chars[2], g1, g2, g3), accola_chars_r(chars[3], g1, g2, g3), accola_chars_r(chars[indices[1][4]], g1, g2, g3)];
end

"""
    accola_chars()

Compute all characteristics in the 8 Accola special theta relations, for a fixed fundamental system.

# Examples
```julia
julia> accola_chars()
```
"""
function accola_chars()
    chars = [[[1,0,0,0,0],[0,0,0,0,0]],
             [[0,1,0,0,0],[1,0,0,0,0]],
             [[0,0,1,0,0],[1,1,0,0,0]],
             [[0,0,0,1,0],[1,1,1,0,0]],
             [[0,0,0,0,1],[1,1,1,1,0]],
             [[1,0,0,0,0],[0,1,1,1,1]],
             [[0,1,0,0,0],[0,0,1,1,1]],
             [[0,0,1,0,0],[0,0,0,1,1]],
             [[0,0,0,1,0],[0,0,0,0,1]],
             [[1,1,1,1,0],[0,1,0,1,0]],
             [[1,1,1,1,0],[1,0,1,0,1]]];
    chars_sr_indices = [[[1,2,3,4], [5,6,7,8], [5,6,9,10], [5,7,9,11]], #SR1234
                        [[1,2,3,5], [4,6,7,8], [4,6,9,10], [4,7,9,11]], #SR1235
                        [[1,2,3,6], [5,4,7,8], [5,4,9,10], [5,7,9,11]], #SR1236
                        [[1,2,3,7], [5,6,4,8], [5,6,9,10], [5,4,9,11]], #SR1237
                        [[1,2,3,8], [5,6,7,4], [5,6,9,10], [5,7,9,11]], #SR1238
                        [[1,2,3,9], [5,6,7,8], [5,6,4,10], [5,7,4,11]], #SR1239
                        [[1,2,3,10], [5,6,7,8], [5,6,9,4], [5,7,9,11]], #SR12310
                        [[1,2,3,11], [5,6,7,8], [5,6,9,10], [5,7,9,4]]]; #SR12311
    return [accola_chars_sr(chars, i) for i in chars_sr_indices];
end

"""
    accola(τ, chars=accola_chars())

Compute the 8 Accola special theta relations, without taking the product over the signs of the square roots. Return the largest absolute value of the special theta relation.

# Arguments
- `τ::Array{<:Number}`: 2-dimensional array of complex numbers
- `chars=accola_chars()`: characteristics in the 8 Accola special theta relations

# Examples
```julia
julia> τ = random_siegel(5);
julia> accola(τ)
```
"""
function accola(τ::Array{<:Number}, chars=accola_chars())
    R = RiemannMatrix(τ, nderivs=0);
    return accola(R, chars);
end

"""
    accola(R, chars=accola_chars())

Compute the 8 Accola special theta relations, without taking the product over the signs of the square roots. Return the largest absolute value of the special theta relation.

# Arguments
- `R::RiemannMatrix`: RiemannMatrix type containing matrix.
- `chars=accola_chars()`: characteristics in the 8 Accola special theta relations

# Examples
```julia
julia> R = RiemannMatrix(random_siegel(5));
julia> accola(R)
```
"""
function accola(R::RiemannMatrix, chars=accola_chars())
    z = zeros(5);
    theta_constants = [[[theta(z, R, char=c) for c in chars_r] for chars_r in chars_sr] for chars_sr in chars];
    accola_r = [[sqrt(prod(theta_constant)) for theta_constant in theta_r] for theta_r in theta_constants];
    accola_sr = [[r[1] + r[2] + r[3] + r[4],
                  r[1] + r[2] + r[3] - r[4],
                  r[1] + r[2] - r[3] + r[4],
                  r[1] + r[2] - r[3] - r[4],
                  r[1] - r[2] + r[3] + r[4],
                  r[1] - r[2] + r[3] - r[4],
                  r[1] - r[2] - r[3] + r[4],
                  r[1] - r[2] - r[3] - r[4]]
                 for r in accola_r];
    max_accola_sr = maximum([minimum(abs.(sr)) for sr in accola_sr]);
    return max_accola_sr;
end

"""
    random_nonaccola(tol=0.1, trials=100)

Find a random genus 5 matrix in the Siegel upper half space which is not in the Accola locus, such that the largest Accola relation has absolute value at least tol, using input number of trials.

# Arguments
- `tol::Real=0.1`: tolerance for deciding whether absolute value of largest Accola relation is nonzero.
- `trials::Integer=100`: number of trials

# Examples
```julia
julia> random_nonaccola(tol=0.1, trials=10)
```
"""
function random_nonaccola(tol::Real=0.1, trials::Integer=100)
    chars = accola_chars();
    t = 0; # stores value of largest Accola relation encountered
    i = 0; # counter for number of trials
    max_matrix = rand(5,5); # stores candidate matrix
    while t < tol && i < trials
        τ = convert(Array, random_siegel(5));
        s = accola(τ, chars);
        if s > t
            t = s;
            max_matrix = τ;
        end
        i += 1;
    end
    return [max_matrix, t];
end


text:
# Run most Stan.jl examples, remove all tmp dirs.

println("Running Stan.jl examples:")

function rmtmpdirs()
  dirs = [
    "Examples/Mamba/Bernoulli",
    "Examples/Mamba/BernoulliDiagnose",
    "Examples/Mamba/BernoulliOptimize",
    "Examples/Mamba/BernoulliVariational",
    "Examples/Mamba/BernoulliInitTheta",
    "Examples/Mamba/BernoulliScalar",
    "Examples/Mamba/Binomial",
    "Examples/Mamba/Binormal",
    "Examples/Mamba/EightSchools",
    "Examples/Mamba/Dyes",
    "Examples/Mamba/ARM/Ch03/Kid"
  ]

  for dir in dirs
    cd(joinpath(Pkg.dir("Stan"), dir)) do
      isdir("tmp") && rm("tmp", recursive=true)
    end
  end
end

examples = [
  "../Examples/Mamba/Bernoulli/bernoulli.jl",
  "../Examples/Mamba/BernoulliOptimize/bernoulli_optimize.jl",
  "../Examples/Mamba/BernoulliDiagnose/bernoulli_diagnose.jl",
  "../Examples/Mamba/BernoulliVariational/bernoulli_variational.jl",
  "../Examples/Mamba/BernoulliInitTheta/bernoulliinittheta.jl",
  "../Examples/Mamba/BernoulliScalar/bernoulliscalar.jl",
  "../Examples/Mamba/Binomial/binomial.jl",
  "../Examples/Mamba/Binormal/binormal.jl",
  "../Examples/Mamba/EightSchools/schools8.jl",
  "../Examples/Mamba/Dyes/dyes.jl",
  "../Examples/Mamba/ARM/Ch03/Kid/kidscore.jl"
]

for example in examples
    println("\n\n\n  * $(example) *")
    include(example)
    rmtmpdirs()
    println("\n$(example) done!\n")
end


text:
using Test

include("../coverage.jl")

### total_combinations
tc_trials = [
    [[2, 2, 2], 2, 12],
    [[2, 3, 2], 2, 16],
    [[2, 3, 2, 2], 2, 6+4+4+6+6+4],
    [[2, 3, 2, 2], 3, 12+12+8+12]
]
for trial in tc_trials
    res = total_combinations(trial[1], trial[2])
    @test res == trial[3]
end

### next_multiplicative
arity0 = [3, 7, 9, 4]
vv0 = copy(arity0)
next_multiplicative!(vv0, arity0)
@test vv0 == [1, 1, 1, 1]
vv = [1,1,1]
nmarity = [2,3,2]
for i in 1:prod(nmarity)
    next_multiplicative!(vv, nmarity)
end
@test vv == [1, 1, 1]

### all_combinations
rng = MersenneTwister(9237425)
for ac_trial_idx in 1:5
    n_way = [2, 3, 4][rand(rng, 1:3)]
    param_cnt = rand(rng, (n_way + 1):(n_way + 3))
    arity = rand(rng, 2:4, param_cnt)
    coverage = all_combinations(arity, n_way)
    # It has the right column dimnsion.
    @test size(coverage, 2) == length(arity)
    # Every combination is nonzero.
    @test sum(sum(coverage, dims = 2) == 0) == 0
    # The total number of combinations agrees with expectations.
    @test total_combinations(arity, n_way) == size(coverage, 1)
    # Generate some random combinations and check that they are in there.
    for comb_idx in 1:100
        comb = [rand(rng, 1:arity[cj]) for cj in 1:param_cnt]
        comb[randperm(rng, param_cnt)[1:(param_cnt - n_way)]] .= 0
        @test sum(comb .!= 0) == n_way
        found = false
        for sidx in 1:size(coverage, 1)
            if coverage[sidx, :] == comb
                found = true
            end
        end
        @test found
    end
end

## multi_way_coverage
mwc = multi_way_coverage([2,3,4,2,2,3], Dict(3 => [[1,4,4,5]]), 2)
@test maximum(sum(mwc .!= 0, dims = 2)) == 3

### coverage_by_paramter
cbp_cases = [
    [[0 0 0; 0 0 0; 0 1 0], 3, [0, 1, 0]],
    [[0 0 0; 0 0 0; 0 1 0], 2, [0, 0, 0]],
    [[0 0 7; 0 0 4; 0 1 0], 3, [0, 1, 2]]
]
for cbp_case in cbp_cases
    res = coverage_by_parameter(cbp_case[1], cbp_case[2])
    @test res == cbp_case[3]
end

### coverage_by_value

cbv_cases = [
    [[0 0 0; 0 0 0; 0 0 0], 3, [2,2,2], 1, [0, 0]],
    [[0 0 0; 0 0 0; 0 0 0], 3, [2,2,2], 2, [0, 0]],
    [[0 0 0; 0 0 0; 0 0 0], 3, [2,3,2], 2, [0, 0, 0]], # sees arity of choice
    [[0 1 0; 0 1 0; 0 0 0], 3, [2,3,2], 2, [2, 0, 0]], # does multiple rows
    [[0 1 0; 0 1 0; 0 3 0], 3, [2,3,2], 2, [2, 0, 1]], # does multiple columns
    [[0 1 0; 0 1 0; 0 3 0], 2, [2,3,2], 2, [2, 0, 0]]  # excludes extra rows
]
for cbv_case in cbv_cases
    res = coverage_by_value(cbv_case[1], cbv_case[2], cbv_case[3], cbv_case[4])
    @test res == cbv_case[5]
end

### most_matches_existing(allc, row_cnt, arity, existing, param_idx)
mmm_allc = [1 1 1 0 0; 1 2 1 0 0; 0 0 1 2 2]
mmm_arity = [2, 2, 2, 2, 2]
mmm_cases = [
    [[1, 0, 0, 0, 0], 2, [1, 1]],
    [[1, 0, 0, 0, 0], 4, [0, 0]],
    [[0, 0, 0, 0, 1], 2, [0, 0]],
    [[0, 0, 0, 0, 1], 3, [0, 0]],
    [[0, 0, 0, 0, 2], 3, [1, 0]],
    [[1, 0, 0, 0, 0], 3, [2, 0]],
    [[0, 0, 1, 0, 0], 1, [2, 0]],
    [[1, 0, 1, 0, 0], 2, [1, 1]],
    [[1, 1, 0, 0, 0], 3, [1, 0]],
    [[1, 2, 0, 0, 0], 3, [1, 0]],
    [[0, 0, 1, 0, 2], 4, [0, 1]]
]
for mmm_case in mmm_cases
    res = most_matches_existing(mmm_allc, 3, mmm_arity, mmm_case[1], mmm_case[2])
    @test res == mmm_case[3]
end

### combination_number

for cn_a in 1:6
    for cn_b in 1:cn_a
        @test combination_number(cn_a, cn_b) == factorial(cn_a) ÷ (factorial(cn_a - cn_b) * factorial(cn_b))
    end
end

### add_coverage!(allc, row_cnt, entry)
ac_cases = [
    [[1 1 0; 1 2 0; 1 0 1; 1 0 2], 4, 2, [1, 1, 1], 2, [1 0 2; 1 2 0; 1 1 0; 1 0 1]],
    [[1 1 0; 1 2 0; 1 0 1; 1 0 2], 3, 2, [1, 1, 1], 1, [1 2 0; 1 1 0; 1 0 1; 1 0 2]]
]
for ac_case in ac_cases
    cover0 = copy(ac_case[1])
    res = add_coverage!(cover0, ac_case[2], ac_case[4])
    @test res == ac_case[5]
    println(cover0)
    println(ac_case[6])
    @test cover0 == ac_case[6]
end

### match_score
ms_cases = [
    [[1 1 0; 1 2 0; 0 1 3], 3, 2, [4,4,4], 0],
    [[1 1 0; 1 2 0; 0 2 1; 0 1 3], 4, 2, [1,1,1], 1],
    [[1 1 0; 1 2 0; 0 2 1; 0 1 3], 4, 2, [1,2,1], 2],
    [[1 1 0; 1 2 0; 0 2 1; 0 1 3], 4, 2, [1,1,3], 2],
    [[1 1 0; 1 2 0; 0 2 1; 0 1 3], 3, 2, [1,1,3], 1]
]
for ms_case in ms_cases
    cnt = match_score(ms_case[1], ms_case[2], ms_case[4])
    @test cnt == ms_case[5]
end

### argmin_rand
ar_cases = [
    [[1,2,3,1,4], [1,4]],
    [[-5, -7, -6, -3], [2]],
    [[0, -1, -2, -2], [3, 4]],
    [[1,1,1,1], [1,2,3,4]]
]
ar_rng = MersenneTwister(342234)
for ar_case in ar_cases
    values = Set{Int}()
    for i in 1:10
        idx = argmin_rand(ar_rng, ar_case[1])
        push!(values, idx)
    end
    res = sort([x for x in values])
    @test res == ar_case[2]
end

## remove_combinations
rc_cases = [
    [[1 1 1; 1 1 2; 1 2 1; 1 2 2], (x -> x[2] == 1 && x[3] == 1), [1 1 2; 1 2 1; 1 2 2]],
    [[1 1 1; 1 1 2; 1 2 1; 1 2 2], (x -> x[2] == 1), [1 2 1; 1 2 2]],
    [[1 1 1; 1 1 2; 1 2 1; 1 2 2], (x -> x[1] == 1 && x[3] == 1), [1 1 2; 1 2 2]],
    [[1 1 1; 1 1 2; 1 2 1; 1 2 2], (x -> x[2] < x[3]), [1 1 1; 1 2 1; 1 2 2]]
]
for rc_case in rc_cases
    res = remove_combinations(rc_case[1], rc_case[2])
    @test res == rc_case[3]
end

### tuples_in_trials
tt_cases = [
    [[[1,2,3], [1,2,1], [2,3,4]], 2, 8],
    [[[1,1,1], [1,1,1], [1,1,1]], 2, 3],
    [[[1,1,1], [1,1,1], [1,1,1]], 3, 1],
]
for tt_case in tt_cases
    res = tuples_in_trials(tt_case[1], tt_case[2])
    total_covered = sum([length(aset) for aset in values(res)])
    @test total_covered == tt_case[3]
end

### n_way_coverage
rng = Random.MersenneTwister(9234724)
arity = [2,3,2,3]
n_way = 2

M = 50
cover = n_way_coverage(arity, n_way, M, rng)
tuple_cnt = coverage_by_tuple(cover, n_way)
@test tuple_cnt == total_combinations(arity, n_way)

n_way = 3
cover = n_way_coverage(arity, n_way, M, rng)
tuple_cnt = coverage_by_tuple(cover, n_way)
@test tuple_cnt == total_combinations(arity, n_way)

cover = n_way_coverage_init(arity, n_way, [], M, rng)
tuple_cnt = coverage_by_tuple(cover, n_way)
@test tuple_cnt == total_combinations(arity, n_way)

seed = [
    1 1 1 1;
    1 3 1 1;
    2 2 2 2;
    2 1 2 1
]
cover = n_way_coverage_init(arity, n_way, seed, M, rng)
tuple_cnt = coverage_by_tuple(cover, n_way)
# It is possible for this to fail randomly.
@test tuple_cnt < total_combinations(arity, n_way)

cover = n_way_coverage_filter(arity, n_way, x -> x[3] == 1 && x[4] != 1, [], M, rng)
n_way_coverage([4,4,4,4,4,4,4,4,4], 2, M, rng)

### n_way_coverage_multi
rng = MersenneTwister(789607)
M = 50
arity = [2,3,4,2,2,3]
n_way = 2
high_indices = [1,3,4,5]
mwc = multi_way_coverage(arity, Dict(3 => [high_indices]), n_way)
@test maximum(sum(mwc .!= 0, dims = 2)) == 3
cov = n_way_coverage_multi(arity, mwc, x->false, [], M, rng)
cov_mat = vcat([cv' for cv in cov]...)
cov_high = cov_mat[:, high_indices]
singles = unique(sort(cov_mat[:, high_indices], dims = 1), dims = 1)

cov_high_list = [cov_high[idx, :] for idx in 1:size(cov_high, 1)]
accumulated = tuples_in_trials(cov_high_list, 3)
tuple_cnt = coverage_by_tuple(cov_high_list, 3)
# This is it. The total number of 3-way tuples should match straight
# count of combinations of the four items for 3-way.
@test tuple_cnt == total_combinations(arity[high_indices], 3)

max_stars_repo_name: adolgert/test-scientific-code

max_issues_repo_name: adolgert/test-scientific-code

max_forks_repo_name: adolgert/test-scientific-code


text:
using LuaCall
using Test

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


text:
##include( "../src/lossFunctions.jl" )
using Statistics
using Test

println( "running loss functions' unit tests" );

#== testing for 0,1 labels error functions ==#
a = zeros( Float64, 10000 )
b = ones( Float64, 10000 )

@test mean( lossFunctions.MSE( a, b ) ) == 1.0
@test mean( lossFunctions.MSE( a, a ) ) == 0.0
@test mean( lossFunctions.MSE( b, b ) ) == 0.0
@test mean( lossFunctions.zeroOneLoss( a, b ) ) == 1.0
@test mean( lossFunctions.zeroOneLoss( a, a ) ) == 0.0
@test mean( lossFunctions.zeroOneLoss( b, b ) ) == 0.0

#== testing for 1,-1 labels error functions ==#
a       = Bool.( rand( 0:1, 10000 ) );
b       = .!a;
index_a = ( a .== false );
a       = Float64.( a );
a[index_a] .= -1.0;
index_b     = ( b .== false );
b           = Float64.( b );
b[index_b] .= -1.0;
@test mean( lossFunctions.hingeLoss( a, b ) ) == 2.0
@test mean( lossFunctions.hingeLoss( a, a ) ) == 0.0
@test mean( lossFunctions.hingeLoss( b, b ) ) == 0.0

@test lossFunctions.MSE( 1.0, 1.0 ) == 0.0
@test lossFunctions.MSE( 0.0, 0.0 ) == 0.0
@test lossFunctions.MSE( 1.0, 0.0 ) == 1.0
@test lossFunctions.MSE( 0.0, 1.0 ) == 1.0
@test lossFunctions.MSE( 0.0, 10.0 ) == 100.0
@test lossFunctions.MSE( 10.0, 0.0 ) == 100.0


text:
# Add MDP to Julia search path
push!(LOAD_PATH, MDPpath * "src/Julia");
push!(LOAD_PATH, MDPpath * "src/Julia/Gencode");
push!(LOAD_PATH, MDPpath * "src/Julia/Preprocessing");
push!(LOAD_PATH, MDPpath * "src/Julia/Postprocessing");
push!(LOAD_PATH, MDPpath * "src/Julia/Potential");
push!(LOAD_PATH, MDPpath * "src/Julia/Optimization");

# Set Julia's PATH enviroment variable so that MDP can call external programs
ENV["PATH"] = "/usr/local/bin:/usr/bin:/opt/local/bin:/bin:/usr/sbin:/sbin";

using Preprocessing, Potential, Optimization

# load Snap library
Potential.loadsnap(MDPpath)

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using  LinearAlgebra
using Plots

function secant_method(f′, x0, x1, ϵ)
    g0 = f′(x0)
    Δ = Inf
    while abs(Δ) > ϵ
        g1 = f′(x1)
        Δ = (x1 - x0)/(g1 - g0)*g1
        x0, x1, g0 = x1, x1 - Δ, g1
    end
    return x1
end

f = x->sin(10x)+cos(3x)
df(x) = 10*cos(10*x)-3*sin(3*x)

#dois pontos iniciais
x0 = 3.4
x1 = 3.8

#convergencia
convergence = 1e-8

xs = secant_method(df,x0,x1,convergence)
xs,f(xs),df(xs)

plot(f,0.,8.,lw=2,draw_arrow="true")
#scatter!([a,c],[f(a),f(c)])
scatter!([x0,x1],[f(x0),f(x1)])
scatter!([xs],[f(xs)])


text:
module problem121

description = """
Disc game prize fund

A bag contains one red disc and one blue disc. In a game of chance a player takes a disc at random and its colour is noted. After each turn the disc is returned to the bag, an extra red disc is added, and another disc is taken at random.
The player pays £1 to play and wins if they have taken more blue discs than red discs at the end of the game.
If the game is played for four turns, the probability of a player winning is exactly 11/120, and so the maximum prize fund the banker should allocate for winning in this game would be £10 before they would expect to incur a loss. Note that any payout will be a whole number of pounds and also includes the original £1 paid to play the game, so in the example given the player actually wins £9.
Find the maximum prize fund that should be allocated to a single game in which fifteen turns are played.
"""

# A brute force calculation of the probability of winning with given number of turns
function pwin(turns :: Int)
  p = zero(Rational)
  for numblues = (div(turns, 2) + 1) : turns		# Win when num blues is over half the turns
    for blueturns  in combinations(1:turns, numblues)   # Find all combinations of turns that could produce the blues
      ps = [n//(n+1) for n = 1:turns]			# Calculate probability of getting a red at every turn
      for t in blueturns				# Update for getting a blue at some turns
        ps[t] = 1 // (t+1)
      end
      p += prod(ps)					# take the product for overall probability of winning exactly this way
    end
  end
  return p
end

# The maximum prize fund that should be allocated given the chance of winning
prizefund(pwin :: Rational) = int(floor(inv(pwin)))

# Solver
solve(n = 15) = println("prize fund for $n turns should be: $(prizefund(pwin(n)))")

using Base.Test

@test pwin(4) == 11//120
@test prizefund(pwin(4)) == 10

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
using Drawing

d = 0.05
k = 2

with(
     PNG("text-align.png", 300, 150),
#     TK(300, 150),
     Axes(centred=true, scale=0.5)) do
    for align in 1:9
        x = k * ((align + 2) % 3 - 1) / 2
        y = (1 - round((align - 2) / 3)) / 2
        draw(Pen(0.03), Ink("orange")) do
            move(x-d, y)
            line(x+d, y)
            move(x, y-d)
            line(x, y+d)
            move(x, y)
        end
        paint(Font(size=0.2)) do
            text("align=$(align)"; align=align)
        end
    end
end