text:
# This file contains routines to generate MDF files
export saveasMDF, loadDataset, loadMetadata, loadMetadataOnline, setparam!, compressCalibMDF
export getFFdataPerPos, prepareAsMDFSingleMeasurement, convertCustomSF, blockAverage

function setparam!(params::Dict, parameter, value)
  if value != nothing
    params[parameter] = value
  end
end

# we do not support all conversion possibilities
function loadDataset(f::MPIFile; frames=1:acqNumFrames(f), applyCalibPostprocessing=false,
                     numPeriodAverages=1, numPeriodGrouping=1, experimentNumber=nothing, kargs...)
  params = loadMetadata(f)

  if experimentNumber != nothing
    params[:experimentNumber] = experimentNumber
  end

  # call API function and store result in a parameter Dict
  if experimentHasMeasurement(f)
    loadMeasParams(f, params, skipMeasData = true)
    if !applyCalibPostprocessing
      if frames!=1:acqNumFrames(f)
        setparam!(params, :measData, measData(f,frames))
        setparam!(params, :acqNumFrames, length(frames))
        setparam!(params, :measIsBGFrame, measIsBGFrame(f)[frames])
      else
        setparam!(params, :measData, measData(f))
      end
    else
        @info "load measurement data"
        data = getMeasurementsFD(f, false, frames=1:acqNumFrames(f), sortFrames=true,
                       spectralLeakageCorrection=false, transposed=true, tfCorrection=false,
                       numPeriodAverages=numPeriodAverages, numPeriodGrouping=numPeriodGrouping)

          @info size(data)
        setparam!(params, :measData, data)
        setparam!(params, :measIsFourierTransformed, true)
        setparam!(params, :measIsFastFrameAxis, true)
        setparam!(params, :measIsFramePermutation, true)
        setparam!(params, :measFramePermutation, fullFramePermutation(f))

        if numPeriodAverages > 1
          params[:acqNumAverages] *= numPeriodAverages
          params[:acqNumPeriodsPerFrame] = div(params[:acqNumPeriodsPerFrame],numPeriodAverages)
        end
        if numPeriodGrouping > 1
          params[:acqNumPeriodsPerFrame] = div(params[:acqNumPeriodsPerFrame],numPeriodGrouping)
          params[:dfCycle] *= numPeriodGrouping # Not sure about this one

          params[:rxNumSamplingPoints] *= numPeriodGrouping
        end

        setparam!(params, :measIsBGFrame,
          cat(zeros(Bool,acqNumFGFrames(f)),ones(Bool,acqNumBGFrames(f)), dims=1))

        snr = calibSNR(f)
	      if snr == nothing
          @info "calculate SNR"
          snr = calculateSystemMatrixSNR(f, data, numPeriodAverages=numPeriodAverages, numPeriodGrouping=numPeriodGrouping)
        end
        setparam!(params, :calibSNR, snr)
    end
  end

  loadCalibParams(f, params)
  loadRecoParams(f, params)

  return params
end

const defaultParams =[:version, :uuid, :time, :dfStrength, :acqGradient, :studyName,
          :studyNumber, :studyUuid, :studyTime, :studyDescription,
          :experimentName, :experimentNumber, :experimentUuid, :experimentDescription,
          :experimentSubject,
          :experimentIsSimulation, :experimentIsCalibration,
          :tracerName, :tracerBatch, :tracerVendor, :tracerVolume, :tracerConcentration,
          :tracerSolute, :tracerInjectionTime,
          :scannerFacility, :scannerOperator, :scannerManufacturer, :scannerName,
          :scannerTopology, :acqNumPeriodsPerFrame, :acqNumAverages,
          :acqStartTime, :acqOffsetField, :acqNumFrames,
          :dfNumChannels, :dfPhase, :dfBaseFrequency, :dfDivider,
          :dfCycle, :dfWaveform, :rxNumChannels, :rxBandwidth,
          :rxNumSamplingPoints, :rxTransferFunction, :rxTransferFunctionFileName, :rxInductionFactor,
          :rxUnit, :rxDataConversionFactor]

function loadMetadata(f, inputParams = MPIFiles.defaultParams)
  params = Dict{Symbol,Any}()
  # call API function and store result in a parameter Dict
  for op in inputParams
    setparam!(params, op, eval(op)(f))
  end
  return params
end

function loadRecoParams(f, params = Dict{Symbol,Any}())
  if experimentHasReconstruction(f)
    for op in [:recoData, :recoSize, :recoFov, :recoFovCenter, :recoOrder,
           :recoPositions, :recoParameters]
      setparam!(params, op, eval(op)(f))
    end
  end

  return params
end

function loadCalibParams(f, params = Dict{Symbol,Any}())
  if experimentIsCalibration(f)
    for op in [:calibFov, :calibFovCenter,
               :calibSize, :calibOrder, :calibPositions, :calibOffsetField,
             :calibDeltaSampleSize, :calibMethod]
      setparam!(params, op, eval(op)(f))
    end
    if !haskey(params, :calibSNR)
      setparam!(params, :calibSNR, calibSNR(f))
    end
  end
  return params
end

function loadMeasParams(f, params = Dict{Symbol,Any}(); skipMeasData = false)
  if experimentHasMeasurement(f)
    for op in [:measIsFourierTransformed, :measIsTFCorrected,
                 :measIsBGCorrected,
                 :measIsFastFrameAxis, :measIsFramePermutation, :measIsFrequencySelection,
                 :measIsSpectralLeakageCorrected,
                 :measFramePermutation, :measIsBGFrame]
      setparam!(params, op, eval(op)(f))
    end
  end

  if !skipMeasData
    setparam!(params, :measData, measData(f))
  end

  return params
end

function appendBGDataset(params::Dict, filenameBG::String; kargs...)
  params = MPIFile(filenameBG) do fBG
    appendBGDataset(params, fBG; kargs...)
  end
  return params
end

function appendBGDataset(params::Dict, fBG::MPIFile; frames=1:acqNumFrames(fBG))
  paramsBG = loadDataset(fBG, frames=frames)
  paramsBG[:measIsBGFrame][:] = true

  params[:measData] = cat(4, params[:measData], paramsBG[:measData])
  params[:measIsBGFrame] = cat(params[:measIsBGFrame], paramsBG[:measIsBGFrame], dims=1)
  params[:acqNumFrames] += paramsBG[:acqNumFrames]

  return params
end

isConvertibleToMDF(f::MPIFile) = true
function isConvertibleToMDF(f::BrukerFile)
  # check if raw data is consistent
  if !rawDataLengthConsistent(f::BrukerFile)
    return false
  end
  # check if conversion of method is supported
  #TODO use regex matching if list grows too large
  whitelist = ["User:ukeMPI337", "User:mkaul_MPI", "User:ukeMPINew", "User:ukeMPI335", "User:uke_mpi", "Bruker:MPICalibration", "Bruker:MPI"]
  if !(f["ACQ_method"] in whitelist)
    @warn "conversion of method not supported" f["ACQ_method"]
    return false
  end
  return true
end

function saveasMDF(filenameOut::String, filenameIn::String; kargs...)
  MPIFile(filenameIn) do f
    saveasMDF(filenameOut, f; kargs...)
  end
  return
end

function saveasMDF(filenameOut::String, f::MPIFile; filenameBG = nothing, enforceConversion=false, kargs...)

  # This is a hack. Needs to be fixed properly
  if(haskey(kargs, :SNRThresh) || haskey(kargs, :sparsityTrafoRedFactor)) && calibSNR(f) != nothing
    compressCalibMDF(filenameOut, f; kargs...)
    return
  end

  if enforceConversion || isConvertibleToMDF(f)
    params = loadDataset(f;kargs...)
    if filenameBG != nothing
      appendBGDataset(params, filenameBG)
    end
    saveasMDF(filenameOut, params)
  else
    error("File not Convertible")
  end
  return
end

function saveasMDF(filename::String, params::Dict)
  # file has to be removed if exists. Otherwise h5create fails.
  isfile(filename) && rm(filename)
  h5open(filename, "w") do file
    saveasMDF(file, params)
  end
end

function compressCalibMDF(filenameOut::String, f::MPIFile; SNRThresh=2.0, kargs...)
  idx = Int64[]

  SNR = calibSNR(f)[:,:,1]
  for k=1:size(SNR,1)
    if maximum(SNR[k,:]) > SNRThresh
      push!(idx, k)
    end
  end

  compressCalibMDF(filenameOut, f, idx; kargs...)
end

function compressCalibMDF(filenamesOut::Vector{String}, f::MultiMPIFile; SNRThresh=2.0, kargs...)
  idx = Int64[]

  # We take the SNR from the first SF
  SNR = calibSNR(f[1])[:,:,1]
  for k=1:size(SNR,1)
    if maximum(SNR[k,:]) > SNRThresh
      push!(idx, k)
    end
  end

  for (i,f_) in enumerate(f)
    compressCalibMDF(filenamesOut[i], f_, idx; kargs...)
  end
end

function compressCalibMDF(filenameOut::String, f::MPIFile, idx::Vector{Int64};
                          sparsityTrafoRedFactor=1.0, sparsityTrafo="DCT-IV", kargs...)
  params = loadMetadata(f)
  loadMeasParams(f, params, skipMeasData = true)
  loadCalibParams(f, params)
  params[:calibSNR] = calibSNR(f)
  loadRecoParams(f, params)

  data = systemMatrixWithBG(f, idx)

  params[:calibSNR] = params[:calibSNR][idx,:,:]
  if haskey(params, :rxTransferFunction)
    params[:rxTransferFunction] = params[:rxTransferFunction][idx,:]
  end
  params[:measIsFrequencySelection] = true
  params[:measFrequencySelection] = idx

  if sparsityTrafoRedFactor == 1.0
    params[:measData] = data
  else
    B = linearOperator(sparsityTrafo, calibSize(f))
    N = prod(calibSize(f))
    NBG = size(data,1) - N
    D = size(data,3)
    P = size(data,4)
    NRed = max(1, floor(Int, sparsityTrafoRedFactor*N))
    dataOut = similar(data, NBG+NRed, length(idx), D, P)
    subsamplingIndices = zeros(Int32, NRed, length(idx), D, P)

    fgdata = data[measFGFrameIdx(f),:,:,:]
    bgdata = data[measBGFrameIdx(f),:,:,:]

    dataOut[(NRed+1):end,:,:,:] = bgdata

    for k=1:length(idx), d=1:D, p=1:P
      I = B * fgdata[:,k,d,p]
      subsamplingIndices[:,k,d,p] = round.(Int32,reverse(sortperm(abs.(I)),dims=1)[1:NRed])
      dataOut[1:NRed,k,d,p] = I[vec(subsamplingIndices[:,k,d,p])]
    end

    params[:measData] = dataOut
    params[:measIsSparsityTransformed] = true
    params[:measSparsityTransformation] = sparsityTrafo
    params[:measSubsamplingIndices] = subsamplingIndices

    bgFrame = zeros(Bool, NRed+NBG)
    bgFrame[(NRed+1):end] .= true
    params[:measIsBGFrame] = bgFrame
    params[:acqNumFrames] = NRed+NBG
  end

  saveasMDF(filenameOut, params)
end

function blockAverage(filenameOut::String, f::MPIFile, numAverages)
	#block averaging is not possible
	if rem(acqNumFrames(f),numAverages) != 0
		throw(DomainError(numAverages, "`numAverages` must be a divider of $(acqNumFrames(f))"))
	end

	# is there a way to rule out calibration scans?

	# get data, exept for measurement data
	params = loadMetadata(f)
	loadMeasParams(f, params, skipMeasData = true)
	loadCalibParams(f, params)
	loadRecoParams(f, params)

	# uuid and time newly created on saving MDF
	delete!(params,:uuid)
	delete!(params,:time)

	# update measurement data
	Nnew = div(acqNumFrames(f),numAverages)
	data = measData(f)
	if measIsSparsityTransformed(f)
		error("No block averaging of sparsity transformed data possible")
	elseif measIsFramePermutation(f)
		error("No block averaging of data with frame permutation possible")
	else
		if measIsFastFrameAxis(f)
			N,J,C,KW = size(data)
			data = reshape(data,numAverages,Nnew,J,C,KW)
			data = mean(data, dims=1) # creates Float64 data if eltype is <: Integer
			data = reshape(data,Nnew,J,C,KW)
		else
			J,C,KW,N = size(data)
			data = reshape(data,J,C,KW,numAverages,Nnew)
			data = mean(data, dims=4) # creates Float64 data if eltype is <: Integer
			data = reshape(data,J,C,KW,Nnew)
		end
		if eltype(data) == Float64
			data = Float32.(data) # Float64 is not necessary
		end
		setparam!(params, :measData, data)
	end

	# update measurement parameters
	p = measIsBGFrame(f)
	p = reshape(p,numAverages,Nnew)
	p = any(p, dims=1) # all blocks with at least one background frames will be labeled as background
	p = reshape(p,Nnew)
	setparam!(params, :measIsBGFrame, p)

	# update acquisition parameters
	setparam!(params, :acqNumFrames, Nnew)
	setparam!(params, :acqNumAverages, acqNumAverages(f)*numAverages)

  return saveasMDF(filenameOut, params)
end

function loadAndProcessFFData(f::BrukerFile, nAverages::Int64, skipSwitchingFrames::Int64;addToEnd=0)
  dataFilename = joinpath(f.path,"rawdata.job0")
  #FileSize = stat(dataFilename).size
  #AllFrames = convert(Int,round(Int,FileSize/26928/4/3))
  #AllFrames = acqNumPeriodsPerPatch(f)*acqNumPeriodsPerFrame(f)
  AllFrames = acqNumPeriodsPerFrame(f)
  (Nx,Ny,Nz) = [length(union(acqOffsetField(f)[ll,1,:])) for ll in collect(1:3)]

  data_ = zeros(ComplexF64,Nx,Ny,Nz,rxNumFrequencies(f),rxNumChannels(f));

#  Pos_=acqOffsetField(f).*[-1,-1,-1]; #Field to position
  Pos=acqOffsetField(f).*[-1,-1,-1]; #Field to position

#  Pos = zeros(3,1,AllFrames);
#  Pos[1,1,:] = kron(vec(Pos_[1,1,:]),ones(2))
#  Pos[2,1,:] = kron(vec(Pos_[2,1,:]),ones(2))
#  Pos[3,1,:] = kron(vec(Pos_[3,1,:]),ones(2))

  dType = Int16 #acqNumAverages(f) == 1 ? Int16 : Int32  # Data are not averaged

  (PosNx,PosNy,PosNz)=[union(Pos[ll,1,:]) for ll in collect(1:3)]

  ds = open(dataFilename)
  if numSubPeriods(f) == 1
     raw = Mmap.mmap(ds, Array{dType,4},(rxNumSamplingPoints(f),1,rxNumChannels(f),AllFrames));
  else
     raw = Mmap.mmap(ds, Array{dType,5},(rxNumSamplingPoints(f),numSubPeriods(f),1,rxNumChannels(f),AllFrames))
     raw = dropdims(sum(raw,dims=2),dims=2)
  end

  for p = collect(1:Nx*Ny*Nz)
#    p = Nx*Ny*(pz-1)+pxy
    st = p*skipSwitchingFrames+(p-1)*(nAverages+addToEnd)+1

    ind = collect(st:st+nAverages-1)

    if any([length(union(Pos[ll,1,ind])) for ll in collect(1:3)].>1)
      error("Averaged data are not at the same FF-Position! Check Parameter Averages and/or skipSwitchingFrames (MPIFiles function loadAndProcessFFData)")
    end
    # sequence is not relevant
    #x = collect(Nx:-1:1)[PosNx.==union(Pos[1,1,ind])]
    #y = collect(Ny:-1:1)[PosNy.==union(Pos[2,1,ind])]
    #z = collect(Nz:-1:1)[PosNz.==union(Pos[3,1,ind])]
    x = collect(1:Nx)[PosNx.==union(Pos[1,1,ind])]
    y = collect(1:Ny)[PosNy.==union(Pos[2,1,ind])]
    z = collect(1:Nz)[PosNz.==union(Pos[3,1,ind])]
    #for ch =1:3
      #dataSC = spectralLeakageCorrectedData(raw[:,1,ch,ind]);
      #data_[x,y,z,:,ch] = rfft(vec(mean(dataSC,dims=2)));
      data_[x,y,z,:,:] = rfft(mean(raw[:,1,:,ind],dims=3),1);
    #end
  end
  data = reshape(data_,:,rxNumFrequencies(f),rxNumChannels(f),1)
  return convert(Array{Complex{Float32},4},data); # allPos x freq x channel x 1
end

function convertCustomSF(filenameOut::String, f::BrukerFile, fBG::BrukerFile,nAverages::Int64,skipSwitchingFrames::Int64; nAveragesBG = nAverages,skipSwitchingFramesBG=skipSwitchingFrames)

  params = loadMetadata(f)
  loadMeasParams(f, params, skipMeasData = true)
  loadCalibParams(f, params)

  paramsBG = loadMetadata(fBG)

  params[:calibSize] = [length(union(params[:acqOffsetField][ll,1,:])) for ll in collect(1:3)]
  numFGFrames = prod(params[:calibSize])

  params[:acqGradient] = reshape(params[:acqGradient][:,:,1,1],3,3,1,1)
  params[:dfPhase] =  reshape(params[:dfPhase][:,:,1],1,3,1)
  params[:calibDeltaSampleSize] = [0.0, 0.0, 0.0] #Todo
  #params[:time] = now()
  params[:calibMethod] = "BrukerCustom"
  #params[:version] = v"2.0.0"
  params[:dfStrength] = reshape(params[:dfStrength][:,:,1,1],1,3,1)
  params[:experimentIsCalibration] = true
  #params[:uuid] = uuid4()
  params[:measIsFourierTransformed] = true
  params[:measIsFastFrameAxis] = true
  params[:calibFovCenter] = [mean(extrema(params[:acqOffsetField][ll,1,:])) for ll in collect(1:3)]
  params[:acqNumPeriodsPerFrame] = 1

  params[:measIsFramePermutation] = 1

  if rxHasTransferFunction(f) &&!measIsTFCorrected(f)
    params[:transferFunction] = rxTransferFunction(f)
  end

  calibFov=[sum(abs.(extrema(params[:acqOffsetField][ll,1,:])))./abs.(params[:acqGradient][ll,ll]) for ll in collect(1:3)]
  ind=findall(x->x==0.0,calibFov) # if calibFov is 0.0 for 2D Measurements, AxisArray will fail in creating an ImageMeta object
  calibFov[ind].=0.001  # for 1mm
  params[:calibFov] = calibFov
  params[:calibOrder] = "xyz"

  params[:acqOffsetField] = reshape([mean(extrema(params[:acqOffsetField][ll,1,:])) for ll in collect(1:3)],3,1,1)
  params[:acqNumAverages] = nAverages

println("Part1")
  # Daten Laden
@time  fgFrames = loadAndProcessFFData(f,nAverages, skipSwitchingFrames) # allePos x freqs x r x 1
@time  bgFramesFull = loadAndProcessFFData(fBG,nAveragesBG, skipSwitchingFramesBG) # allePos x freqs x r x 1
  nAveragesBGHalf = nAveragesBG >1 ? div(nAveragesBG,2) : nAveragesBG
@time  bgFramesHalf = loadAndProcessFFData(fBG,div(nAveragesBG,2), skipSwitchingFramesBG+nAveragesBG-div(nAveragesBG,2)) # allePos x freqs x r x 1

  params[:calibSNR] = calculateSNRCustomSF(f,fgFrames,bgFramesFull,bgFramesHalf)

  (xBG,yBG,zBG) = [length(union(paramsBG[:acqOffsetField][ll,1,:])) for ll in collect(1:3)]
  (xFG,yFG,zFG) = params[:calibSize]

  bgFullReshaped = reshape(bgFramesFull,xBG,yBG,zBG,rxNumFrequencies(fBG),rxNumChannels(fBG))

  if yBG !=1
     if zBG !=1
       itp = interpolate(bgFullReshaped,(NoInterp(),BSpline(Linear()),BSpline(Linear()),NoInterp(),NoInterp())); #Interpolire bgFramesFull auf
     else
       itp = interpolate(bgFullReshaped,(NoInterp(),BSpline(Linear(    )),NoInterp(),NoInterp(),NoInterp())); #Interpolire bgFrames    Full auf
     end
 else
   if zBG !=1
      itp = interpolate(bgFullReshaped,(NoInterp(),NoInterp(),BSpline(Linear()),NoInterp(),NoInterp())); #Interpolire bgFrames    Full auf
   else
      itp = interpolate(bgFullReshaped,(NoInterp(),NoInterp(),NoInterp(),NoInterp(),NoInterp())); #Interpolire bgFrames        Full auf
   end
 end

  bgFramesFullInterp = itp(collect(1:xBG),range(1,yBG,length=yFG),range(1,zBG,length=zFG),collect(1:size(bgFullReshaped,4)),collect(1:size(bgFullReshaped,5)))

  params[:acqNumFrames] = numFGFrames + prod(size(bgFramesFullInterp)[1:3])

  tt = [round.(Int,collect(range(1,xBG+xFG,length=xBG))).+(kk-1)*(xBG+xFG) for kk =1:yFG*zFG];
  idxBGFrames = vcat(tt...);
  idxAllFrames = collect(1:params[:acqNumFrames])
  idxAllFrames[idxBGFrames] .= 0
  idxFGFrames = idxAllFrames[idxAllFrames.!=0];

  params[:measFramePermutation] =  vcat(idxFGFrames,idxBGFrames)
  params[:measIsBGFrame] = vcat(zeros(numFGFrames),ones(length(idxBGFrames)))
  dataTemp_ = vcat(fgFrames,reshape(bgFramesFullInterp,:,size(bgFullReshaped)[4:5]...,1))  #size(params[:measData])(22621, 26929, 3, 1)
  params[:measData] = convert(Array{Complex{Float32},4},dataTemp_);  #size(params[:measData])(22621, 26929, 3, 1)
  saveasMDF(filenameOut, params)
end

function convertCustomSF(filenameOut::String, f::Array{BrukerFileMeas}, fBG::BrukerFile,nAverages::Int64,skipSwitchingFrames::Int64; nAveragesBG = nAverages,skipSwitchingFramesBG=skipSwitchingFrames)

  paramsArray =[];
  paramsBG = loadMetadata(fBG)

  params = loadMetadata(f[1])
  loadMeasParams(f[1], params, skipMeasData = true)
  loadCalibParams(f[1], params)

  for fi in f[1:end]

  params_ = loadMetadata(fi)
  loadMeasParams(fi, params_, skipMeasData = true)
  loadCalibParams(fi, params_)
  push!(paramsArray,params_);
  end

  for i =2:length(f)
  params[:acqOffsetField] = cat(params[:acqOffsetField],paramsArray[i][:acqOffsetField],dims=3)
  end

  params[:calibSize] = [length(union(params[:acqOffsetField][ll,1,:])) for ll in collect(1:3)]
  numFGFrames = prod(params[:calibSize])

  println(params[:calibSize])

  params[:acqGradient] = reshape(params[:acqGradient][:,:,1,1],3,3,1,1)
  params[:dfPhase] =  reshape(params[:dfPhase][:,:,1],1,3,1)
  params[:calibDeltaSampleSize] = [0.0, 0.0, 0.0] #Todo
  #params[:time] = now()
  params[:calibMethod] = "BrukerCustom"
  #params[:version] = v"2.0.0"
  params[:dfStrength] = reshape(params[:dfStrength][:,:,1,1],1,3,1)
  params[:experimentIsCalibration] = true
  #params[:uuid] = uuid4()
  params[:measIsFourierTransformed] = true
  params[:measIsFastFrameAxis] = true
  params[:calibFovCenter] = [mean(extrema(params[:acqOffsetField][ll,1,:])) for ll in collect(1:3)]
  params[:acqNumPeriodsPerFrame] = 1

  params[:measIsFramePermutation] = 1

  if rxHasTransferFunction(f[1]) &&!measIsTFCorrected(f[1])
    params[:transferFunction] = sampleTF(TransferFunction("/opt/data/TFs/TFHDF5/TF_3DMouse_NM0_07102019.h5"), f[1])
#    params[:transferFunction] = rxTransferFunction(f[1])
  end

  params[:calibFov] = [sum(abs.(extrema(params[:acqOffsetField][ll,1,:])))./abs.(params[:acqGradient][ll,ll]) for ll in collect(1:3)]
  params[:calibOrder] = "xyz"

  params[:acqOffsetField] = reshape([mean(extrema(params[:acqOffsetField][ll,1,:])) for ll in collect(1:3)],3,1,1)
  params[:acqNumAverages] = nAverages

println("Part1")
  # Daten Laden
#fgFrames =zeros(ComplexF64,prod(Int64,params[:calibSize]),26929,3,1)
@time global fgFrames =loadAndProcessFFData(f[1],nAverages, skipSwitchingFrames)
#global tmp =1
for (ii,ffi)  in enumerate(f[2:end])
println(ffi)
#pos = prod(Int64,[length(union(paramsArray[ii][:acqOffsetField][ll,1,:])) for ll in collect(1:3)])
#@time fgFrames[tmp:tmp+pos-1,:,:,:] = loadAndProcessFFData(ffi,nAverages, skipSwitchingFrames) # allePos x freqs x r x 1
@time global fgFrames = cat(fgFrames,loadAndProcessFFData(ffi,nAverages, skipSwitchingFrames),dims=1) # allePos x freqs x r x 1
#global tmp = tmp+pos
end

@time  bgFramesFull = loadAndProcessFFData(fBG,nAveragesBG, skipSwitchingFramesBG) # allePos x freqs x r x 1
  nAveragesBGHalf = nAveragesBG >1 ? div(nAveragesBG,2) : nAveragesBG
@time  bgFramesHalf = loadAndProcessFFData(fBG,div(nAveragesBG,2), skipSwitchingFramesBG+nAveragesBG-div(nAveragesBG,2)) # allePos x freqs x r x 1

  params[:calibSNR] = calculateSNRCustomSF(f[1],fgFrames,bgFramesFull,bgFramesHalf)

  (xBG,yBG,zBG) = [length(union(paramsBG[:acqOffsetField][ll,1,:])) for ll in collect(1:3)]
  (xFG,yFG,zFG) = params[:calibSize]

  bgFullReshaped = reshape(bgFramesFull,xBG,yBG,zBG,rxNumFrequencies(fBG),rxNumChannels(fBG))

  itp = interpolate(bgFullReshaped,(NoInterp(),BSpline(Linear()),BSpline(Linear()),NoInterp(),NoInterp())); #Interpolire bgFramesFull auf

  bgFramesFullInterp = itp(collect(1:xBG),range(1,yBG,length=yFG),range(1,zBG,length=zFG),collect(1:size(bgFullReshaped,4)),collect(1:size(bgFullReshaped,5)))

  params[:acqNumFrames] = numFGFrames + prod(size(bgFramesFullInterp)[1:3])

  tt = [round.(Int,collect(range(1,xBG+xFG,length=xBG))).+(kk-1)*(xBG+xFG) for kk =1:yFG*yFG];
  idxBGFrames = vcat(tt...);
  idxAllFrames = collect(1:params[:acqNumFrames])
  idxAllFrames[idxBGFrames] .= 0
  idxFGFrames = idxAllFrames[idxAllFrames.!=0];

  params[:measFramePermutation] =  vcat(idxFGFrames,idxBGFrames)
  params[:measIsBGFrame] = vcat(zeros(numFGFrames),ones(length(idxBGFrames)))
  dataTemp_ = vcat(fgFrames,reshape(bgFramesFullInterp,:,size(bgFullReshaped)[4:5]...,1))  #size(params[:measData])(22621, 26929, 3, 1)
  params[:measData] = convert(Array{Complex{Float32},4},dataTemp_);  #size(params[:measData])(22621, 26929, 3, 1)
  saveasMDF(filenameOut, params)
end

"""
This function returns the time data `timeData` Array(Float32,numPos,samplingPoints,channels,1,numFrames)
of an FF-measurement with `numPos` as number of FF-positions. The function splits the entire FF-measurement
in measurement per FF-position. As a second return argument the corresponding positions vector
`pos` Array(Float64,3,1,numPos) is returned.
"""
function getFFdataPerPos(f::MPIFile, nAverages::Int64, skipSwitchingFrames::Int64; addToEnd=0)
    dataFilename = joinpath(f.path,"rawdata.job0")
    AllFrames = acqNumPeriodsPerFrame(f)
    numPos_= size(acqOffsetField(f),3)
    numPos = convert(Int64,div(numPos_,nAverages+skipSwitchingFrames))
    timeData_ = zeros(Float32,numPos,rxNumSamplingPoints(f),rxNumChannels(f),1,nAverages)

    posFF = acqOffsetField(f).*[-1,-1,-1] #Field to position
    pos = zeros(Float64,3,1,numPos)

    measTime = measData(f)

    for p = collect(1:numPos)
        st = p*skipSwitchingFrames+(p-1)*(nAverages+addToEnd)+1
        ind = collect(st:st+nAverages-1)
        pos_ = unique(posFF[:,:,ind],dims=3)
        if size(pos_,3) != 1
            error("FF Positions not unqiue! Check Average and skipSwitchingFrames Parameter!")
        end
        pos[:,:,p] = pos_
        tempTimeData = measTime[:,:,ind,:]
        timeData_[p,:,:,:,:] = permutedims(tempTimeData, [1,2,4,3])# switch last two dimensions: frames and patches switch
    end
    timeData = convert(Array{Int16,5},timeData_)
    return timeData, pos
end

"""
This function spilts a FF-measurement according to it FF-positions in measurements per position
and treats them as single measurements. It returns `fileNames` Array{String,1}
and the correspodning `paramsArray` as Array{Dict,1}, which are the data dictionaries.
Both can be used to save the data with saveasMDF(...).
"""
function prepareAsMDFSingleMeasurement(f::MPIFile, numFrames::Int64, skipSwitchingFrames::Int64; addToEnd=0)
    global params = loadMetadata(f)

    # get FF data
    timeDataPerPos, pos = getFFdataPerPos(f, numFrames, skipSwitchingFrames, addToEnd=addToEnd);

    # adjust meta data
    params[:acqGradient] = acqGradient(f)[:,:,1:1,1:1]
    params[:acqNumPeriodsPerFrame] = 1
    params[:acqNumFrames] = numFrames
    params[:dfStrength] = dfStrength(f)[:,:,1:1]
    params[:dfPhase] = dfPhase(f)[:,:,1:1]

    # adjust data and meta data per position
    paramsArray = Array{Any,1}(undef,size(pos,3))
    fileNames = Array{String,1}(undef,size(pos,3))
    for p=1:size(pos,3)
        tempParas=deepcopy(params)
        filename = string(params[:experimentNumber],"_pos",p)
        tempParas[:acqOffsetField] = pos[:,:,p:p] # all
        tempParas[:measIsFrequencySelection] = false
        tempParas[:measIsBGCorrected] = false
        tempParas[:measIsFramePermutation] = false
        tempParas[:measIsFastFrameAxis] = false
        tempParas[:measIsFourierTransformed] = false
        tempParas[:measIsSpectralLeakageCorrected] = false
        tempParas[:measIsBGFrame] = convert(Array{Bool,1},zeros(size(timeDataPerPos,5)));
        tempParas[:measIsTFCorrected] = false
        tempParas[:measData] = timeDataPerPos[p,:,:,:,:];

        fileNames[p] = filename
        paramsArray[p] = tempParas
    end
    return fileNames, paramsArray
end

hasKeyAndValue(paramDict,param) = haskey(paramDict, param) && paramDict[param] != nothing

function writeIfAvailable(file, paramOut, paramDict, paramIn )
  if hasKeyAndValue(paramDict, paramIn)
    write(file, paramOut, paramDict[paramIn])
  end
end

function saveasMDF(file::HDF5.File, params::Dict{T,Any}) where  {T<:AbstractString}
  paramsSymbol = Dict{Symbol,Any}([Symbol(k)=>v for (k,v) in params])
  saveasMDF(file, paramsSymbol)
end

function saveasMDF(file::HDF5.File, params::Dict{Symbol,Any})
  # general parameters
  write(file, "/version", "2.0.1")
  write(file, "/uuid", string(get(params,:uuid,uuid4() )))
  write(file, "/time", "$( get(params,:time, Dates.unix2datetime(time())) )")

  # study parameters
  write(file, "/study/name", get(params,:studyName,"default") )
  write(file, "/study/number", get(params,:studyNumber,0))
  if hasKeyAndValue(params,:studyUuid)
    studyUuid = params[:studyUuid]
  else
    studyUuid = uuid4()
  end
  write(file, "/study/uuid", string(studyUuid))
  write(file, "/study/description", get(params,:studyDescription,"n.a."))
  if hasKeyAndValue(params,:studyTime)
    write(file, "/study/time", string(params[:studyTime]))
  end

  # experiment parameters
  write(file, "/experiment/name", get(params,:experimentName,"default") )
  write(file, "/experiment/number", get(params,:experimentNumber,0))
  if hasKeyAndValue(params,:experimentUuid)
    expUuid = params[:experimentUuid]
  else
    expUuid = uuid4()
  end
  write(file, "/experiment/uuid", string(expUuid))
  write(file, "/experiment/description", get(params,:experimentDescription,"n.a."))
  write(file, "/experiment/subject", get(params,:experimentSubject,"n.a."))
  write(file, "/experiment/isSimulation", Int8(get(params,:experimentIsSimulation,false)))

  # tracer parameters
  write(file, "/tracer/name", get(params,:tracerName,"n.a") )
  write(file, "/tracer/batch", get(params,:tracerBatch,"n.a") )
  write(file, "/tracer/vendor", get(params,:tracerVendor,"n.a") )
  write(file, "/tracer/volume", get(params,:tracerVolume,0.0))
  write(file, "/tracer/concentration", get(params,:tracerConcentration,0.0) )
  write(file, "/tracer/solute", get(params,:tracerSolute,"Fe") )
  tr = [string(t) for t in get(params,:tracerInjectionTime, [Dates.unix2datetime(time())]) ]
  write(file, "/tracer/injectionTime", tr)

  # scanner parameters
  write(file, "/scanner/facility", get(params,:scannerFacility,"n.a") )
  write(file, "/scanner/operator", get(params,:scannerOperator,"n.a") )
  write(file, "/scanner/manufacturer", get(params,:scannerManufacturer,"n.a"))
  write(file, "/scanner/name", get(params,:scannerName,"n.a"))
  write(file, "/scanner/topology", get(params,:scannerTopology,"FFP"))

  # acquisition parameters
  write(file, "/acquisition/numAverages",  params[:acqNumAverages])
  write(file, "/acquisition/numFrames", get(params,:acqNumFrames,1))
  write(file, "/acquisition/numPeriods", get(params,:acqNumPeriodsPerFrame,1))
  write(file, "/acquisition/startTime", "$( get(params,:acqStartTime, Dates.unix2datetime(time())) )")

  writeIfAvailable(file, "/acquisition/gradient", params, :acqGradient)
  writeIfAvailable(file, "/acquisition/offsetField", params, :acqOffsetField)

  # drivefield parameters
  write(file, "/acquisition/drivefield/numChannels", size(params[:dfStrength],2) )
  write(file, "/acquisition/drivefield/strength", params[:dfStrength])
  write(file, "/acquisition/drivefield/phase", params[:dfPhase])
  write(file, "/acquisition/drivefield/baseFrequency", params[:dfBaseFrequency])
  write(file, "/acquisition/drivefield/divider", params[:dfDivider])
  write(file, "/acquisition/drivefield/cycle", params[:dfCycle])
  if !haskey(params, :dfWaveform)
    params[:dfWaveform] = "sine"
  end
  write(file, "/acquisition/drivefield/waveform", params[:dfWaveform])

  # receiver parameters
  write(file, "/acquisition/receiver/numChannels", params[:rxNumChannels])
  write(file, "/acquisition/receiver/bandwidth", params[:rxBandwidth])
  write(file, "/acquisition/receiver/numSamplingPoints", params[:rxNumSamplingPoints])
  if !haskey(params, :rxUnit)
    params[:rxUnit] = "V"
  end
  write(file, "/acquisition/receiver/unit",  params[:rxUnit])
  write(file, "/acquisition/receiver/dataConversionFactor",  params[:rxDataConversionFactor])
  if hasKeyAndValue(params,:rxTransferFunction)
    tf = params[:rxTransferFunction]
    group = file["/acquisition/receiver"]
    writeComplexArray(group, "transferFunction", tf)
  end
  writeIfAvailable(file, "/acquisition/receiver/transferFunctionFileName", params, :rxTransferFunctionFileName)
  writeIfAvailable(file, "/acquisition/receiver/inductionFactor",  params, :rxInductionFactor)

  # measurements
  if hasKeyAndValue(params, :measData)
    meas = params[:measData]
    if eltype(meas) <: Complex
      group = create_group(file,"/measurement")
      writeComplexArray(group, "/measurement/data", meas)
    else
      write(file, "/measurement/data", meas)
    end
    write(file, "/measurement/isFourierTransformed", Int8(params[:measIsFourierTransformed]))
    write(file, "/measurement/isSpectralLeakageCorrected", Int8(params[:measIsSpectralLeakageCorrected]))
    write(file, "/measurement/isTransferFunctionCorrected", Int8(params[:measIsTFCorrected]))
    write(file, "/measurement/isFrequencySelection", Int8(params[:measIsFrequencySelection]))
    write(file, "/measurement/isBackgroundCorrected",  Int8(params[:measIsBGCorrected]))
    write(file, "/measurement/isFastFrameAxis",  Int8(params[:measIsFastFrameAxis]))
    write(file, "/measurement/isFramePermutation",  Int8(params[:measIsFramePermutation]))
    writeIfAvailable(file, "/measurement/frequencySelection",  params, :measFrequencySelection)

    if hasKeyAndValue(params, :measFramePermutation)
      write(file, "/measurement/framePermutation", params[:measFramePermutation] )
    end
    if hasKeyAndValue(params, :measIsBGFrame)
      write(file, "/measurement/isBackgroundFrame", convert(Array{Int8}, params[:measIsBGFrame]) )
    end
    if hasKeyAndValue(params, :measIsSparsityTransformed)
      write(file, "/measurement/isSparsityTransformed", params[:measIsSparsityTransformed] )
      write(file, "/measurement/subsamplingIndices", params[:measSubsamplingIndices] )
      write(file, "/measurement/sparsityTransformation", params[:measSparsityTransformation] )
    end
  end

  # calibrations
  writeIfAvailable(file, "/calibration/snr",  params, :calibSNR)
  writeIfAvailable(file, "/calibration/fieldOfView",  params, :calibFov)
  writeIfAvailable(file, "/calibration/fieldOfViewCenter",  params, :calibFovCenter)
  writeIfAvailable(file, "/calibration/size",  params, :calibSize)
  writeIfAvailable(file, "/calibration/order",  params, :calibOrder)
  writeIfAvailable(file, "/calibration/positions",  params, :calibPositions)
  writeIfAvailable(file, "/calibration/offsetField",  params, :calibOffsetField)
  writeIfAvailable(file, "/calibration/deltaSampleSize",  params, :calibDeltaSampleSize)
  writeIfAvailable(file, "/calibration/method",  params, :calibMethod)
  if hasKeyAndValue(params, :calibIsMeanderingGrid)
    write(file, "/calibration/isMeanderingGrid", Int8(params[:calibIsMeanderingGrid]))
  end
  writeIfAvailable(file, "/calibration/_temperatures",  params, :calibTemperatures)
  # reconstruction
  if hasKeyAndValue(params, :recoData)
    write(file, "/reconstruction/data", params[:recoData])
    write(file, "/reconstruction/fieldOfView", params[:recoFov])
    write(file, "/reconstruction/fieldOfViewCenter", params[:recoFovCenter])
    write(file, "/reconstruction/size", params[:recoSize])
    write(file, "/reconstruction/order", get(params, :recoOrder, "xyz"))
    if hasKeyAndValue(params, :recoPositions)
      write(file, "/reconstruction/positions", params[:recoPositions])
    end
    if hasKeyAndValue(params, :recoParameters)
      saveParams(file, "/reconstruction/_parameters", params[:recoParameters])
    end
  end

  writeIfAvailable(file, "/custom/auxiliaryData", params, :auxiliaryData)
end


text:
module Compliance

import ..Clarus

function clearingeligible(;params...)
  return Clarus.api_request("Compliance","ClearingEligible";params...)
end

function clearingmandatory(;params...)
  return Clarus.api_request("Compliance","ClearingMandatory";params...)
end

function sefmandatory(;params...)
  return Clarus.api_request("Compliance","SEFMandatory";params...)
end

end #MODULE-END

ext: jl


text:
using LinearAlgebra

import ClimaCore:
    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies
import ClimaCore.Geometry: ⊗

using OrdinaryDiffEq: ODEProblem, solve, SSPRK33

import Logging
import TerminalLoggers
Logging.global_logger(TerminalLoggers.TerminalLogger())

const parameters = (
    ϵ = 0.1,  # perturbation size for initial condition
    l = 0.5, # Gaussian width
    k = 0.5, # Sinusoidal wavenumber
    ρ₀ = 1.0, # reference density
    c = 2,
    g = 10,
)

domain = Domains.RectangleDomain(
    Domains.IntervalDomain(
        Geometry.XPoint(-2π),
        Geometry.XPoint(2π),
        periodic = true,
    ),
    Domains.IntervalDomain(
        Geometry.YPoint(-2π),
        Geometry.YPoint(2π),
        periodic = true,
    ),
)

n1, n2 = 16, 16
Nq = 4
Nqh = 7
mesh = Meshes.RectilinearMesh(domain, n1, n2)
grid_topology = Topologies.Topology2D(mesh)
quad = Spaces.Quadratures.GLL{Nq}()
space = Spaces.SpectralElementSpace2D(grid_topology, quad)

Iquad = Spaces.Quadratures.GLL{Nqh}()
Ispace = Spaces.SpectralElementSpace2D(grid_topology, Iquad)

function init_state(coord, p)
    x, y = coord.x, coord.y
    # set initial state
    ρ = p.ρ₀

    # set initial velocity
    U₁ = cosh(y)^(-2)

    # Ψ′ = exp(-(y + p.l / 10)^2 / 2p.l^2) * cos(p.k * x) * cos(p.k * y)
    # Vortical velocity fields (u₁′, u₂′) = (-∂²Ψ′, ∂¹Ψ′)
    gaussian = exp(-(y + p.l / 10)^2 / 2p.l^2)
    u₁′ = gaussian * (y + p.l / 10) / p.l^2 * cos(p.k * x) * cos(p.k * y)
    u₁′ += p.k * gaussian * cos(p.k * x) * sin(p.k * y)
    u₂′ = -p.k * gaussian * sin(p.k * x) * cos(p.k * y)

    u = Geometry.UVVector(U₁ + p.ϵ * u₁′, p.ϵ * u₂′)
    # set initial tracer
    θ = sin(p.k * y)

    return (ρ = ρ, ρu = ρ * u, ρθ = ρ * θ)
end

y0 = init_state.(Fields.coordinate_field(space), Ref(parameters))

function flux(state, param)
    ρ, ρu, ρθ = state.ρ, state.ρu, state.ρθ
    u = ρu / ρ
    return (
        ρ = ρu,
        ρu = ((ρu ⊗ u) + (param.g * ρ^2 / 2) * LinearAlgebra.I),
        ρθ = ρθ * u,
    )
end

function energy(state, param)
    ρ, ρu = state.ρ, state.ρu
    u = ρu / ρ
    return ρ * (u.u^2 + u.v^2) / 2 + param.g * ρ^2 / 2
end

function total_energy(y, parameters)
    sum(energy.(y, Ref(parameters)))
end

function rhs!(dydt, y, _, t)

    I = Operators.Interpolate(Ispace)
    div = Operators.WeakDivergence()
    R = Operators.Restrict(space)

    rparameters = Ref(parameters)

    @. dydt = -R(div(flux(I(y), rparameters)))

    Spaces.weighted_dss!(dydt)
    return dydt
end

# Next steps:
# 1. add the above to the design docs (divergence + over-integration + DSS)
# 2. add boundary conditions

dydt = similar(y0)
rhs!(dydt, y0, nothing, 0.0)

# Solve the ODE operator
prob = ODEProblem(rhs!, y0, (0.0, 80.0))
sol = solve(
    prob,
    SSPRK33(),
    dt = 0.02,
    saveat = 1.0,
    progress = true,
    progress_message = (dt, u, p, t) -> t,
)

ENV["GKSwstype"] = "nul"
using ClimaCorePlots, Plots
Plots.GRBackend()

dirname = "cg"
path = joinpath(@__DIR__, "output", dirname)
mkpath(path)

anim = Plots.@animate for u in sol.u
    Plots.plot(u.ρθ, clim = (-1, 1))
end
Plots.mp4(anim, joinpath(path, "tracer.mp4"), fps = 10)

Es = [total_energy(u, parameters) for u in sol.u]
Plots.png(Plots.plot(Es), joinpath(path, "energy.png"))

function linkfig(figpath, alt = "")
    # buildkite-agent upload figpath
    # link figure in logs if we are running on CI
    if get(ENV, "BUILDKITE", "") == "true"
        artifact_url = "artifact://$figpath"
        print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
    end
end

linkfig(
    relpath(joinpath(path, "energy.png"), joinpath(@__DIR__, "../..")),
    "Total Energy",
)


text:
######################################################################
# Abstraction of GLFW.Monitor
# -----
# Licensed under MIT License
using GLFW

export VideoMode
struct VideoMode
    width::Int
    height::Int
    bits::NamedTuple{(:red, :green, :blue), Tuple{UInt8, UInt8, UInt8}}
    refresh_rate::Int
end

export Monitor, monitor
struct Monitor <: GLFWWrapper{GLFW.Monitor}
    handle::GLFW.Monitor
end
monitor(idx::Integer = 1) = Monitors[idx]
monitor(::Nothing) = nothing

export Monitors
"""`Monitors` is a meta type designed to query connected physical monitors using `Base.getindex` and `Base.iterate`."""
struct Monitors end

Base.iterate(::Type{Monitors}, idx = 1) = Monitors[idx], idx < length(Monitors) ? idx+1 : nothing
Base.getindex(::Type{Monitors}, idx::Integer) = collect(Monitors)[idx]
Base.length(::Type{Monitors}) = length(GLFW.GetMonitors())

export videomode, videomodes
videomode(monitor::Monitor) = videomode(GLFW.GetVideoMode(monitor.handle))
videomode(vidmode::GLFW.VidMode) = VideoMode(vidmode.width, vidmode.height, (red = vidmode.redbits, green = vidmode.greenbits, blue = vidmode.bluebits), vidmode.refreshrate)
videomodes(monitor::Monitor) = videomode.(GLFW.GetVideoModes(monitor.handle))

Base.convert(::Type{GLFW.Monitor}, monitor::Monitor) = monitor.handle
Base.convert(::Type{GLFW.Monitor}, ::Nothing) = GLFW.Monitor(C_NULL)
Base.collect(::Type{Monitors}) = collect(Monitor.(GLFW.GetMonitors()))


text:
module LazyBandedMatrices
using BandedMatrices, BlockBandedMatrices, LazyArrays, ArrayLayouts, MatrixFactorizations, LinearAlgebra, Base

import MatrixFactorizations: ql, ql!, QLPackedQ, QRPackedQ, reflector!, reflectorApply!

import Base: BroadcastStyle, similar, OneTo, copy
import Base.Broadcast: Broadcasted
import LinearAlgebra: kron, hcat, vcat, AdjOrTrans, AbstractTriangular, BlasFloat, BlasComplex, BlasReal,
                        lmul!, rmul!

import ArrayLayouts: materialize!, colsupport, rowsupport, MatMulVecAdd, require_one_based_indexing, sublayout
import LazyArrays: LazyArrayStyle, combine_mul_styles, mulapplystyle, PaddedLayout,
                        broadcastlayout, applylayout, arguments, _arguments,
                        LazyArrayApplyStyle, ApplyArrayBroadcastStyle, ApplyStyle,
                        LazyLayout, ApplyLayout, BroadcastLayout, FlattenMulStyle,
                        _mul_args_rows, _mul_args_cols, paddeddata, factorizestyle,
                        MulMatrix, Mul, CachedMatrix, CachedArray, resizedata!, applybroadcaststyle
import BandedMatrices: bandedcolumns, bandwidths, isbanded, AbstractBandedLayout,
                        prodbandwidths, BandedStyle, BandedColumns, BandedRows,
                        AbstractBandedMatrix, BandedSubBandedMatrix, BandedStyle, _bnds,
                        banded_rowsupport, banded_colsupport, _BandedMatrix, bandeddata,
                        banded_qr_lmul!, banded_qr_rmul!, banded_qr
import BlockBandedMatrices: AbstractBlockBandedLayout, BlockSlice, Block1,
                        isblockbanded, isbandedblockbanded, blockbandwidths,
                        subblockbandwidths, blocksizes, BlockSizes

BroadcastStyle(::LazyArrayStyle{1}, ::BandedStyle) = LazyArrayStyle{2}()
BroadcastStyle(::BandedStyle, ::LazyArrayStyle{1}) = LazyArrayStyle{2}()
BroadcastStyle(::LazyArrayStyle{2}, ::BandedStyle) = LazyArrayStyle{2}()
BroadcastStyle(::BandedStyle, ::LazyArrayStyle{2}) = LazyArrayStyle{2}()

bandedcolumns(::ML) where ML<:LazyLayout = BandedColumns{ML}()
bandedcolumns(::ML) where ML<:ApplyLayout = BandedColumns{LazyLayout}()

for LazyLay in (:(BandedColumns{LazyLayout}), :(BandedRows{LazyLayout}),
                :(TriangularLayout{UPLO,UNIT,BandedRows{LazyLayout}} where {UPLO,UNIT}),
                :(TriangularLayout{UPLO,UNIT,BandedColumns{LazyLayout}} where {UPLO,UNIT}))
    @eval begin
        combine_mul_styles(::$LazyLay) = LazyArrayApplyStyle()
        mulapplystyle(::QLayout, ::$LazyLay) = LazyArrayApplyStyle()
    end
end

BroadcastStyle(M::ApplyArrayBroadcastStyle{2}, ::BandedStyle) = M
BroadcastStyle(::BandedStyle, M::ApplyArrayBroadcastStyle{2}) = M

bandwidths(M::Mul) = min.(_bnds(M), prodbandwidths(M.args...))

isbanded(K::Kron{<:Any,2}) = all(isbanded, K.args)
function bandwidths(K::Kron{<:Any,2})
    A,B = K.args
    (size(B,1)*bandwidth(A,1) + max(0,size(B,1)-size(B,2))*size(A,1)   + bandwidth(B,1),
        size(B,2)*bandwidth(A,2) + max(0,size(B,2)-size(B,1))*size(A,2) + bandwidth(B,2))
end

const BandedMatrixTypes = (:AbstractBandedMatrix, :(AdjOrTrans{<:Any,<:AbstractBandedMatrix}),
                                    :(AbstractTriangular{<:Any, <:AbstractBandedMatrix}),
                                    :(Symmetric{<:Any, <:AbstractBandedMatrix}))

const OtherBandedMatrixTypes = (:Zeros, :Eye, :Diagonal, :SymTridiagonal)

for T1 in BandedMatrixTypes, T2 in BandedMatrixTypes
    @eval kron(A::$T1, B::$T2) = BandedMatrix(Kron(A,B))
end

for T1 in BandedMatrixTypes, T2 in OtherBandedMatrixTypes
    @eval kron(A::$T1, B::$T2) = BandedMatrix(Kron(A,B))
end

for T1 in OtherBandedMatrixTypes, T2 in BandedMatrixTypes
    @eval kron(A::$T1, B::$T2) = BandedMatrix(Kron(A,B))
end

###
# Specialised multiplication for arrays padded for zeros
# needed for ∞-dimensional banded linear algebra
###

function similar(M::MulAdd{<:AbstractBandedLayout,<:PaddedLayout}, ::Type{T}, axes) where T
    A,x = M.A,M.B
    xf = paddeddata(x)
    n = max(0,min(length(xf) + bandwidth(A,1),length(M)))
    Vcat(Vector{T}(undef, n), Zeros{T}(size(A,1)-n))
end

function materialize!(M::MatMulVecAdd{<:AbstractBandedLayout,<:PaddedLayout,<:PaddedLayout})
    α,A,x,β,y = M.α,M.A,M.B,M.β,M.C
    length(y) == size(A,1) || throw(DimensionMismatch())
    length(x) == size(A,2) || throw(DimensionMismatch())

    ỹ = paddeddata(y)
    x̃ = paddeddata(x)

    length(ỹ) ≥ min(length(M),length(x̃)+bandwidth(A,1)) ||
        throw(InexactError("Cannot assign non-zero entries to Zero"))

    materialize!(MulAdd(α, view(A, axes(ỹ,1), axes(x̃,1)) , x̃, β, ỹ))
    y
end

###
# MulMatrix
###

bandwidths(M::MulMatrix) = bandwidths(Applied(M))
isbanded(M::Mul) = all(isbanded, M.args)
isbanded(M::MulMatrix) = isbanded(Applied(M))

struct MulBandedLayout <: AbstractBandedLayout end
applylayout(::Type{typeof(*)}, ::AbstractBandedLayout...) = MulBandedLayout()

applybroadcaststyle(::Type{<:AbstractMatrix}, ::MulBandedLayout) = BandedStyle()
# applybroadcaststyle(::Type{<:AbstractMatrix}, ::MulLayout{<:Tuple{BandedColumns{LazyLayout},Vararg{<:AbstractBandedLayout}}}) = LazyArrayStyle{2}()

@inline colsupport(::MulBandedLayout, A, j) = banded_colsupport(A, j)
@inline rowsupport(::MulBandedLayout, A, j) = banded_rowsupport(A, j)
# @inline colsupport(::MulLayout{<:Tuple{<:AbstractBandedLayout,<:AbstractStridedLayout}}, A, j) = banded_colsupport(A, j)
@inline _arguments(::MulBandedLayout, A) = arguments(A)

###
# BroadcastMatrix
###

bandwidths(M::BroadcastMatrix) = bandwidths(Broadcasted(M))
isbanded(M::BroadcastMatrix) = isbanded(Broadcasted(M))

struct BroadcastBandedLayout{F} <: AbstractBandedLayout end
struct LazyBandedLayout <: AbstractBandedLayout end

broadcastlayout(::Type{F}, ::AbstractBandedLayout) where F = BroadcastBandedLayout{F}()
for op in (:*, :/, :\)
    @eval broadcastlayout(::Type{typeof($op)}, ::AbstractBandedLayout, ::AbstractBandedLayout) = BroadcastBandedLayout{typeof($op)}()
end
for op in (:*, :/)
    @eval broadcastlayout(::Type{typeof($op)}, ::AbstractBandedLayout, ::Any) = BroadcastBandedLayout{typeof($op)}()
end
for op in (:*, :\)
    @eval broadcastlayout(::Type{typeof($op)}, ::Any, ::AbstractBandedLayout) = BroadcastBandedLayout{typeof($op)}()
end
broadcastlayout(::Type{typeof(*)}, ::AbstractBandedLayout, ::LazyLayout) = LazyBandedLayout()
broadcastlayout(::Type{typeof(*)}, ::LazyLayout, ::AbstractBandedLayout) = LazyBandedLayout()
broadcastlayout(::Type{typeof(/)}, ::AbstractBandedLayout, ::LazyLayout) = LazyBandedLayout()
broadcastlayout(::Type{typeof(\)}, ::LazyLayout, ::AbstractBandedLayout) = LazyBandedLayout()

# functions that satisfy f(0,0) == 0
for op in (:+, :-)
    @eval broadcastlayout(::Type{typeof($op)}, ::AbstractBandedLayout, ::AbstractBandedLayout) = BroadcastBandedLayout{typeof($op)}()
end

mulapplystyle(::LazyBandedLayout, ::MulBandedLayout) = FlattenMulStyle()
mulapplystyle(::MulBandedLayout, ::LazyBandedLayout) = FlattenMulStyle()

@inline colsupport(::BroadcastBandedLayout, A, j) = banded_colsupport(A, j)
@inline rowsupport(::BroadcastBandedLayout, A, j) = banded_rowsupport(A, j)

###
# sub materialize
###

function arguments(::MulBandedLayout, V::SubArray)
    P = parent(V)
    kr, jr = parentindices(V)
    as = arguments(P)
    kjr = intersect.(_mul_args_rows(kr, as...), _mul_args_cols(jr, reverse(as)...))
    view.(as, (kr, kjr...), (kjr..., jr))
end

@inline sub_materialize(::MulBandedLayout, V) = BandedMatrix(V)
@inline sub_materialize(::BroadcastBandedLayout, V) = BandedMatrix(V)

_BandedMatrix(::MulBandedLayout, V::AbstractMatrix) = apply(*, map(BandedMatrix,arguments(V))...)
for op in (:+, :-)
    @eval @inline _BandedMatrix(::BroadcastBandedLayout{typeof($op)}, V::AbstractMatrix) = apply($op, map(BandedMatrix,arguments(V))...)
end

function arguments(::BroadcastBandedLayout, V::SubArray)
    A = parent(V)
    kr, jr = parentindices(V)
    view.(arguments(A), Ref(kr), Ref(jr))
end

sublayout(M::MulBandedLayout, ::Type{<:Tuple{Vararg{AbstractUnitRange}}}) = M
sublayout(M::BroadcastBandedLayout, ::Type{<:Tuple{Vararg{AbstractUnitRange}}}) = M

######
# Concat banded matrix
######

bandwidths(M::Hcat) = (bandwidth(M.args[1],1),sum(size.(M.args[1:end-1],2)) + bandwidth(M.args[end],2))
isbanded(M::Hcat) = all(isbanded, M.args)

bandwidths(M::Vcat) = (sum(size.(M.args[1:end-1],1)) + bandwidth(M.args[end],1), bandwidth(M.args[1],2))
isbanded(M::Vcat) = all(isbanded, M.args)

const HcatBandedMatrix{T,N} = Hcat{T,NTuple{N,BandedMatrix{T,Matrix{T},OneTo{Int}}}}
const VcatBandedMatrix{T,N} = Vcat{T,2,NTuple{N,BandedMatrix{T,Matrix{T},OneTo{Int}}}}

BroadcastStyle(::Type{HcatBandedMatrix{T,N}}) where {T,N} = BandedStyle()
BroadcastStyle(::Type{VcatBandedMatrix{T,N}}) where {T,N} = BandedStyle()

Base.replace_in_print_matrix(A::HcatBandedMatrix, i::Integer, j::Integer, s::AbstractString) =
    -bandwidth(A,1) ≤ j-i ≤ bandwidth(A,2) ? s : Base.replace_with_centered_mark(s)
Base.replace_in_print_matrix(A::VcatBandedMatrix, i::Integer, j::Integer, s::AbstractString) =
    -bandwidth(A,1) ≤ j-i ≤ bandwidth(A,2) ? s : Base.replace_with_centered_mark(s)

hcat(A::BandedMatrix...) = BandedMatrix(Hcat(A...))
hcat(A::BandedMatrix, B::AbstractMatrix...) = Matrix(Hcat(A, B...))

vcat(A::BandedMatrix...) = BandedMatrix(Vcat(A...))
vcat(A::BandedMatrix, B::AbstractMatrix...) = Matrix(Vcat(A, B...))

#######
# CachedArray
#######

bandwidths(B::CachedMatrix) = bandwidths(B.data)
isbanded(B::CachedMatrix) = isbanded(B.data)

function resizedata!(B::CachedMatrix{T,BandedMatrix{T,Matrix{T},OneTo{Int}}}, n::Integer, m::Integer) where T<:Number
    @boundscheck checkbounds(Bool, B, n, m) || throw(ArgumentError("Cannot resize beyound size of operator"))

    # increase size of array if necessary
    olddata = B.data
    ν,μ = B.datasize
    n,m = max(ν,n), max(μ,m)

    if (ν,μ) ≠ (n,m)
        l,u = bandwidths(B.array)
        λ,ω = bandwidths(B.data)
        if n ≥ size(B.data,1) || m ≥ size(B.data,2)
            M = 2*max(m,n+u)
            B.data = _BandedMatrix(reshape(resize!(vec(olddata.data), (λ+ω+1)*M), λ+ω+1, M), M+λ, λ,ω)
        end
        if ν > 0 # upper-right
            kr = max(1,μ+1-ω):ν
            jr = μ+1:min(m,ν+ω)
            if !isempty(kr) && !isempty(jr)
                view(B.data, kr, jr) .= B.array[kr, jr]
            end
        end
        view(B.data, ν+1:n, μ+1:m) .= B.array[ν+1:n, μ+1:m]
        if μ > 0
            kr = ν+1:min(n,μ+λ)
            jr = max(1,ν+1-λ):μ
            if !isempty(kr) && !isempty(jr)
                view(B.data, kr, jr) .= B.array[kr, jr]
            end
        end
        B.datasize = (n,m)
    end

    B
end

include("bandedql.jl")

###
# BlockBanded
###

mulapplystyle(::DiagonalLayout, ::AbstractBlockBandedLayout) = MulAddStyle()
mulapplystyle(::AbstractBlockBandedLayout, ::DiagonalLayout) = MulAddStyle()

isblockbanded(K::Kron{<:Any,2}) = isbanded(first(K.args))
isbandedblockbanded(K::Kron{<:Any,2}) = all(isbanded, K.args)
blockbandwidths(K::Kron{<:Any,2}) = bandwidths(first(K.args))
subblockbandwidths(K::Kron{<:Any,2}) = bandwidths(last(K.args))
function blocksizes(K::Kron{<:Any,2})
    A,B = K.args
    BlockSizes(Fill(size(B,1), size(A,1)), Fill(size(B,2), size(A,2)))
end

const SubKron{T,M1,M2,R1,R2} =
    SubArray{T,2,<:Kron{T,2,Tuple{M1,M2}},Tuple{BlockSlice{R1},BlockSlice{R2}}}

BroadcastStyle(::Type{<:SubKron{<:Any,<:Any,B,Block1,Block1}}) where B =
    BroadcastStyle(B)

@inline bandwidths(V::SubKron{<:Any,<:Any,<:Any,Block1,Block1}) =
    subblockbandwidths(parent(V))

end


text:
"""
    SphereHalfSpaceCollision

    collision between a spherical contact and a flat surface

    contact_tangent: mapping from world frame to surface tangent frame
    contact_normal: inverse/complement of contact_tangent
    contact_origin: position of contact on Body relative to center of mass
    contact_radius: radius of contact
"""
mutable struct SphereHalfSpaceCollision{T,O,I,OI} <: Collision{T,O,I,OI}
    contact_tangent::SMatrix{O,I,T,OI}
    contact_normal::Adjoint{T,SVector{I,T}}
    contact_origin::SVector{I,T}
    contact_radius::T
end

# function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, collision::SphereHalfSpaceCollision)
#     summary(io, collision)
#     println(io, "")
#     println(io, "contact_tangent: "*string(collision.contact_tangent))
#     println(io, "contact_normal:  "*string(collision.contact_normal))
#     println(io, "contact_origin:  "*string(collision.contact_origin))
#     println(io, "contact_radius:  "*string(collision.contact_radius))
# end

# distance
function distance(collision::SphereHalfSpaceCollision, xp, qp, xc, qc)
    collision.contact_normal * (xp + vector_rotate(collision.contact_origin, qp)) - collision.contact_radius
end

function ∂distance∂x(gradient::Symbol, collision::SphereHalfSpaceCollision, xp, qp, xc, qc)
    if gradient == :parent
        return collision.contact_normal
    elseif gradient == :child
        return szeros(eltype(xc), 1, 3)
    end
end

function ∂distance∂q(gradient::Symbol, collision::SphereHalfSpaceCollision, xp, qp, xc, qc)
    if gradient == :parent
        return collision.contact_normal * ∂vector_rotate∂q(collision.contact_origin, qp)
    elseif gradient == :child
        return szeros(eltype(qc), 1, 4)
    end
end

# contact point in world frame
function contact_point(relative::Symbol, collision::SphereHalfSpaceCollision, xp, qp, xc, qc)
    if relative == :parent
        return xp + vector_rotate(collision.contact_origin, qp) - collision.contact_normal' * collision.contact_radius
    elseif relative == :child
        projector = collision.contact_tangent' * collision.contact_tangent
        return projector * (xp + vector_rotate(collision.contact_origin, qp))
    end
end

function ∂contact_point∂x(relative::Symbol, jacobian::Symbol, collision::SphereHalfSpaceCollision, xp, qp, xc, qc)
    if relative == :parent
        if jacobian == :parent
            return 1.0 * I(3)
        elseif jacobian == :child
            return szeros(eltype(xp), 3, 3)
        end
    elseif relative == :child
        if jacobian == :parent
            projector = collision.contact_tangent' * collision.contact_tangent
            return projector
        elseif jacobian == :child
            return szeros(eltype(xp), 3, 3)
        end
    end
end

function ∂contact_point∂q(relative::Symbol, jacobian::Symbol, collision::SphereHalfSpaceCollision, xp, qp, xc, qc)
    if relative == :parent
        if jacobian == :parent
            return ∂vector_rotate∂q(collision.contact_origin, qp)
        end
    elseif relative == :child
        if jacobian == :parent
            projector = collision.contact_tangent' * collision.contact_tangent
            return projector * ∂vector_rotate∂q(collision.contact_origin, qp)
        elseif jacobian == :child
            return szeros(eltype(qp), 3, 4)
        end
    end
end

# normal projection (from child to parent)
function contact_normal(collision::SphereHalfSpaceCollision, xp, qp, xc, qc)
    return collision.contact_normal
end

function ∂contact_normal_transpose∂x(jacobian::Symbol, collision::SphereHalfSpaceCollision, xp, qp, xc, qc)
    return szeros(eltype(collision.contact_normal), 3, 3)
end

function ∂contact_normal_transpose∂q(jacobian::Symbol, collision::SphereHalfSpaceCollision, xp, qp, xc, qc)
    return szeros(eltype(collision.contact_normal), 3, 4)
end

# tangent projection
function contact_tangent(collision::SphereHalfSpaceCollision, xp, qp, xc, qc)
    return collision.contact_tangent # {2,4} x 3
end

function ∂contact_tangent_one_transpose∂x(jacobian::Symbol, collision::SphereHalfSpaceCollision, xp, qp, xc, qc)
    return szeros(eltype(collision.contact_tangent), 3, 3)
end

function ∂contact_tangent_two_transpose∂x(jacobian::Symbol, collision::SphereHalfSpaceCollision, xp, qp, xc, qc)
    return szeros(eltype(collision.contact_tangent), 3, 3)
end

function ∂contact_tangent_one_transpose∂q(jacobian::Symbol, collision::SphereHalfSpaceCollision, xp, qp, xc, qc)
    return szeros(eltype(collision.contact_tangent), 3, 4)
end

function ∂contact_tangent_two_transpose∂q(jacobian::Symbol, collision::SphereHalfSpaceCollision, xp, qp, xc, qc)
    return szeros(eltype(collision.contact_tangent), 3, 4)
end


text:
using LabConnections.BeagleBone

using Base.Test

# This flag is enabled if a dummy filesystem should be used for testing (for online testing)
# disabling the flag allows the BBB to be run in the loop, in this case blinking LEDS
RUNNING_TESTS = true

# Run tests
include("BeagleBone/SYS_LED_test.jl")
include("BeagleBone/GPIO_test.jl")
include("BeagleBone/PWM_test.jl")

# Remove the testfilesystem
try
    rm("$(pwd())/testfilesystem", recursive=true)
catch
    println("Warning. Could not remove the testfilesystem.")
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct Sku <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    tier::Any # spec type: Union{ Nothing, String } # spec name: tier
    capacity::Any # spec type: Union{ Nothing, Int64 } # spec name: capacity

    function Sku(;name=nothing, tier=nothing, capacity=nothing)
        o = new()
        validate_property(Sku, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(Sku, Symbol("tier"), tier)
        setfield!(o, Symbol("tier"), tier)
        validate_property(Sku, Symbol("capacity"), capacity)
        setfield!(o, Symbol("capacity"), capacity)
        o
    end
end # type Sku

const _property_map_Sku = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("tier")=>Symbol("tier"), Symbol("capacity")=>Symbol("capacity"))
const _property_types_Sku = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("tier")=>"String", Symbol("capacity")=>"Int64")
Base.propertynames(::Type{ Sku }) = collect(keys(_property_map_Sku))
Swagger.property_type(::Type{ Sku }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_Sku[name]))}
Swagger.field_name(::Type{ Sku }, property_name::Symbol) =  _property_map_Sku[property_name]

function check_required(o::Sku)
    true
end

function validate_property(::Type{ Sku }, name::Symbol, val)
end


text:
"""
Validate the args and kwargs for nndescent.
"""
function validate_args(data, n_neighbors, metric, max_iters, sample_rate, precision)
    max_iters >= 1 || error("max_iters must be greater than 0")
    0. < sample_rate ≤ 1. || error("sample_rate must be in (0., 1.]")
    0. ≤ precision ≤ 1. || error("precision must be in [0., 1.]")
    return
end


text:
#=
hello_world:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-05
=#

module HelloWorld
    # methods
    function main()
        println("Hello world!!")
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .HelloWorld
    HelloWorld.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
using Printf

function hist(numbers)
    maxwidth = 50
    h = fill(0, 10)
    for n in numbers
        h[ceil(Int, 10n)] += 1
    end
    mx = maximum(h)
    for (n, i) in enumerate(h)
        @printf("%3.1f: %s\n", n / 10, "+" ^ floor(Int, i / mx * maxwidth))
    end
end

for i in 1:6
    n = rand(10 ^ i)
    println("\n##\n## $(10 ^ i) numbers")
    @printf("μ: %8.6f; σ: %8.6f\n", mean(n), std(n))
    hist(n)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# coding:utf-8

doc =
"""
  対話データ収集プログラム
    http://www.englishfor2day.com/article/dialogue/?pg=1
"""

using Requests
using HttpCommon
using DataFrames
using PyCall
@pyimport bs4

# 対話データパース
# return dataframe
function get_dialogue(url)
  persons, utterances = [],[]
  request = get(url)
  html = readall(request)
  soup = bs4.BeautifulSoup(html)
  dialogue = convert(AbstractString, soup["find"]("div", class="dtlContent")["text"])
  for line in split(dialogue, r"\n|\r")[2:end]
    line = rstrip(line)
    (length(line) == 0 || ismatch(r"\d:\d",line) || !contains(line, ":")) && continue
    person, utterance = split(line, ":")
    push!(persons, person), push!(utterances, utterance)
  end
  @time sleep(0.3)
  return DataFrame(persons = persons, utterances = utterances)
end

# 対話データ記載ページ収集
# page_s = 1, page_e = 23
function dialogue_scraping(page_s, page_e)
  page, num = page_s, 0
  regex = r"""<a href="http://www.englishfor2day.com/article/dialogue/\d+".style"""
  while page <= page_e
    @printf "%s\n" "http://www.englishfor2day.com/article/dialogue/?pg=$page"
    request = get("http://www.englishfor2day.com/article/dialogue/?pg=$page")
    html, _ = split(readall(request), "<!--==========ad=========-->")
    dialogue_pages = [replace(m, r"^.+?\"|\".+?$", "") for m in matchall(regex, html)]
    for url in dialogue_pages
      @printf "%d : %s" num url
      df, num = get_dialogue(url), num + 1
      writetable("./dialogue_corpus/$num.dat", df, separator = '\t', header = true)
    end
    page += 1
  end
end

try
  mkdir("dialogue_corpus")
end
dialogue_scraping(1, 23)

max_stars_repo_name: kwnsiy/scraping

max_issues_repo_name: kwnsiy/scraping

max_forks_repo_name: kwnsiy/scraping


text:
#=
main:
- Julia version:
- Author: quangio
- Date: 2019-05-06
=#

using JuMP
using SCIP
using DelimitedFiles
include("input.jl")

model = Model(with_optimizer(SCIP.Optimizer))

dist_mat_gr17 = [
  0 633 257  91 412 150  80 134 259 505 353 324  70 211 268 246 121
633   0 390 661 227 488 572 530 555 289 282 638 567 466 420 745 518
257 390   0 228 169 112 196 154 372 262 110 437 191  74  53 472 142
 91 661 228   0 383 120  77 105 175 476 324 240  27 182 239 237  84
412 227 169 383   0 267 351 309 338 196  61 421 346 243 199 528 297
150 488 112 120 267   0  63  34 264 360 208 329  83 105 123 364  35
 80 572 196  77 351  63   0  29 232 444 292 297  47 150 207 332  29
134 530 154 105 309  34  29   0 249 402 250 314  68 108 165 349  36
259 555 372 175 338 264 232 249   0 495 352  95 189 326 383 202 236
505 289 262 476 196 360 444 402 495   0 154 578 439 336 240 685 390
353 282 110 324  61 208 292 250 352 154   0 435 287 184 140 542 238
324 638 437 240 421 329 297 314  95 578 435   0 254 391 448 157 301
 70 567 191  27 346  83  47  68 189 439 287 254   0 145 202 289  55
211 466  74 182 243 105 150 108 326 336 184 391 145   0  57 426  96
268 420  53 239 199 123 207 165 383 240 140 448 202  57   0 483 153
246 745 472 237 528 364 332 349 202 685 542 157 289 426 483   0 336
121 518 142  84 297  35  29  36 236 390 238 301  55  96 153 336   0] # shoud return 2085
dist_mat_5 = [
0.0  3.0  4.0  2.0  7.0
3.0  0.0  4.0  6.0  3.0
4.0  4.0  0.0  5.0  8.0
2.0  6.0  5.0  0.0  6.0
7.0  3.0  8.0  6.0  0.0
] # shoud return 19.0

n = 22 # number of vertices
c = TestSet.eil22_capacities / 100  # cost assigned to each vertex (weight needed to be picked up)
d = TestSet.loadVrpEil22() # distance matrix of the graph

w0 = 100
no_salesman = 5
cap = 300

range = collect(1:n)

for i = 1:n
    d[i, i] = 0
    # c[i] = rand() * 100
end
c[1] = 0

println(string("solving TSP for capacities: ", c))

@variable(model, y[1:n, 1:n])
@variable(model, x[1:n, 1:n], Bin)

@objective(model, Min, sum(y[i, j] * d[i, j] for i = 1:n, j = range[1:n .!= i]))

@constraint(model, sum(x[1, i] for i = 2:n) == no_salesman)
@constraint(model, sum(x[i, 1] for i = 2:n) == no_salesman)

for i = 2:n
    @constraint(model, y[1, i] == w0 * x[1, i])
    @constraint(model, sum(x[i, j] for j = range[1:n .!= i]) == 1)
    @constraint(model, sum(x[j, i] for j = range[1:n .!= i]) == 1)
    @constraint(model, c[i] == sum(y[i, j] - y[j, i] for j = range[1:n .!= i]))
end

for i = 1:n
    for j = 1:n
        if i == j continue end
        @constraint(model, y[i, j] >= (w0 + c[i]) * x[i, j])
        @constraint(model, y[i, j] <= (w0 + cap - c[j]) * x[i, j])
    end
end

function solved(m)
    x_val = JuMP.value.(x)
    cycle_idx = [1] # we can adapt this to support mTSP, just keep this for simplicity for now
    while true
        v, idx = findmax(x_val[cycle_idx[length(cycle_idx)], :])
        if idx == 1
            break
        else
            push!(cycle_idx, idx)
        end
    end
    if length(cycle_idx) < n
        @constraint(m, sum(x[i, j] for i = cycle_idx, j = cycle_idx[cycle_idx .!= i]) <= length(cycle_idx) - 1)
        return false
    end
    return true
end

if no_salesman == 0
    JuMP.optimize!(model)
    while !solved(model)
        JuMP.optimize!(model)
    end
end

if no_salesman >= 1000
    @variable(model, 0 <= u[1:n] <= n - 1)
    for i = 2:n
        for j = collect(2:n)[2:n .!= i]
            @constraint(model, u[i] - u[j] + n * x[i, j] <= n - 1)
        end
    end
    @time optimize!(model)
end

optimize!(model)

println(string("Objective value:", JuMP.objective_value(model)))

println("Arcs: ")

for i = range
    for j = range
        if value(x[i, j]) > 0.0
            println(string(i, "->", j, "[label=", value(y[i, j]), "]")) # you can use graphviz to visualize this
        end
    end
end

println("Arcs with capacity: ")

for i = range
    for j = range
        if value(x[i, j]) > 0.0
            println(string('"',i, '|',round(c[i]), '"', "->", '"',j, '|',round(c[j]), '"', "[label=", round(value(y[i, j])), "]")) # you can use graphviz to visualize this
        end
    end
end

# open("model2.lp", "w") do f
#    print(f, model)
# end

max_stars_repo_name: quangIO/mTSP

max_issues_repo_name: quangIO/mTSP

max_forks_repo_name: quangIO/mTSP


text:
using Distributed
@everywhere include("MultiFidelityABC.jl")

bm = MakeBenchmarkCloud("viral/output")
idx = Iterators.partition(1:length(bm), 1000)
epsilons = (0.25, 0.25)
println("--- Loaded benchmark")

pmap(i->MakeMFABCCloud(Viral.mf_prob, bm[[i]], epsilons, method="mf", lower_eta=0.01, budget=1.01*cost(bm[[i]])), 1:nworkers());
println("--- Good to go")

println("--- Small burn-in")
mfabc_set = pmap(i->MakeMFABCCloud(Viral.mf_prob, bm[i], epsilons, method="mf", lower_eta=0.01, budget=2*cost(bm[i])), idx)
println("--- Completed extra simulations: small burn-in")
println("--- Writing files")
for i in 1:length(mfabc_set)
 write_cloud(mfabc_set[i], "viral/output/mf_smallBI/cloud$i")
end

println("--- Large burn-in")
mfabc_set = pmap(i->MakeMFABCCloud(Viral.mf_prob, bm, epsilons, method="mf", lower_eta=0.01, budget=cost(bm)+cost(bm[i])), idx)
println("--- Completed extra simulations: large burn-in")
println("--- Writing files")
for i in 1:length(mfabc_set)
 write_cloud(mfabc_set[i], "viral/output/mf_largeBI/cloud$i")
end

max_stars_repo_name: tpprescott/mf-abc

max_issues_repo_name: tpprescott/mf-abc

max_forks_repo_name: tpprescott/mf-abc


text:
using stoosc
rc("font", family="arial")

## Plot square
fig, ax = subplots(figsize=(1, 0.5))
t = [0.5, 1, 1, 2, 2, 3, 3, 4, 4, 4.5]
x = [0, 0, 1, 1, 0, 0, 1, 1, 0, 0]
ax.plot(t, x, color="black", linewidth=2)
ax.set_axis_off()
fig.tight_layout(pad=0)
save_figure(fig, "figure01", "square.svg")

## Plot sine
fig, ax = subplots(figsize=(1, 0.5))
t = 0:0.01:4π
x = -cos.(t)
ax.plot(t, x, color="black", linewidth=2)
ax.set_axis_off()
fig.tight_layout(pad=0)
save_figure(fig, "figure01","sine.svg")

## Load model and set up its parameters
model = load_model("model01", "sde")
model = set_timespan(model, (0.0, 25.0))
events = create_events([(:LD, 25, 0.5, 0.5)])
model = set_input(model, events)
model = set_parameter_value(model, "I", 0.05)
model = set_parameter_value(model, "σ", 0.01)

## Plot time series for different population sizes
fig, ax_arr = subplots(3, figsize=(5, 4))

min_x = 0.07
max_x = 0.12
events_memory = nothing
for (i, n) = enumerate([1, 10, 100])
    global events_memory, max_x, min_x
    local events, x, t
    solution = simulate_population(model, n, save_trajectories=true, show_progress=true)
    solution = select_time(solution, mintime=19)
    t = solution.time
    x = solution.mean[:, 3]
    U = solution.trajectories[:, 3, :]
    events = solution.events
    for j = 1:n
        if i == 2
            ax_arr[i].plot(t, U[:, j], color="gray", alpha=0.5)
        elseif i == 3
            ax_arr[i].plot(t, U[:, j], color="gray", alpha=0.1)
        end
    end
    ax_arr[i].plot(t, x, color="black")

    if minimum(x) < min_x
        min_x = minimum(U)
    end

    if maximum(x) > max_x
        max_x = maximum(U)
    end

    if isnothing(events_memory)
        events_memory = events
    end

end

for (i, n) = enumerate([1, 10, 100])
    ax_arr[i].set_title("n = $(n)", pad=0)
    ax_arr[i].set_xlabel("Time (a.u.)", labelpad=0)
    ax_arr[i].set_ylabel("Z (a.u.)", labelpad=0)
    ax_arr[i].set_ylim(min_x, max_x)
    ax_arr[i].set_xlim(0, 6)
    ax_arr[i].set_yticks([0.07, 0.12])
    plot_events(events_memory, ax=ax_arr[i], ylims=[min_x, max_x])
    if i == 1
        ax_arr[i].legend(["Single-cell trajectories", "Population-level mean"],
        edgecolor="black", framealpha=1.0, ncol=1, loc=2)
    end
end
fig.tight_layout()

save_figure(fig, "figure01", "population_size.svg")

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
module test_julia_fix2

using Test

@test in(2) isa Base.Fix2
@test in(2).(1:5) == [false, true, false, false, false]

if VERSION >= v"1.3"
    ispositive = >(0)
    isnegative = <(0)
    @test ispositive(1)
    @test isnegative(-1)
end

end # module test_julia_fix2


text:
module KalmanFilter

#
# References:
# Haykin, Kalman Filter and NN (2001)
# https://en.wikipedia.org/wiki/Kalman_filter#Fixed-interval_smoothers
# Hartikainen et al., Optimal filtering with Kalman filters and smoothers - A Manual for Matlab toolbox EKF/UKF (2011)
# Martin, Modified Bryson-Frazier Smoother Cross-Covariance. IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 59, NO. 1, JANUARY 2014
#

using LinearAlgebra

include("CommonStructures.jl")
export CholeskySqrt, SVDSqrt, CholeskyModSqrt

# Review the use of eye, because it could be replaced by I
# unpack and unpack_matrix could be merged into just unpack accepting different methods
include("Utils.jl")
export eye, unpack, unpack_matrix

include("FixedIntervalSmoothers.jl")
export kalman_smooth,
       RauchTungStriebel, ModifiedBrysonFrazier,
       RTSSmoother, MBFSmoother

include("LinearKalmanFilter.jl")
export linear_kalman, LinearKF

include("UnscentedKalmanFilter.jl")
export unscented_kalman, UnscentedKF

end # module KalmanFilter


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.Model{Float64}())
config = MOIT.TestConfig()

bridged_mock = MOIB.Variable.RSOCtoSOC{Float64}(mock)

@testset "rotatedsoc4" begin
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock,
        [√2, 0.0, 1.0, 1.0],
        (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-1.0])
    MOIT.rotatedsoc4test(bridged_mock, config)

    ceqs = MOI.get(mock, MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}())
    @test length(ceqs) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), ceqs[1], "c")

    @testset "Test mock model" begin
        var_names = ["a", "b", "c", "d"]
        MOI.set(
            mock, MOI.VariableName(),
            MOI.get(mock, MOI.ListOfVariableIndices()), var_names)
        socs = MOI.get(mock, MOI.ListOfConstraintIndices{
            MOI.VectorOfVariables, MOI.SecondOrderCone}())
        @test length(socs) == 1
        MOI.set(mock, MOI.ConstraintName(), socs[1], "soc")

        s2 = √2
        s = """
        variables: a, b, c, d
        soc: [a, b, c, d] in MathOptInterface.SecondOrderCone(4)
        c: $s2*a <= 2.0
        maxobjective: c + d
        """
        model = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model, s)
        MOIU.test_models_equal(mock, model, var_names, ["soc", "c"])
    end

    @testset "Test bridged model" begin
        var_names = ["t", "u", "x", "y"]
        MOI.set(
            bridged_mock, MOI.VariableName(),
            MOI.get(bridged_mock, MOI.ListOfVariableIndices()), var_names)
        rsocs = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{
            MOI.VectorOfVariables, MOI.RotatedSecondOrderCone}())
        @test length(rsocs) == 1
        MOI.set(bridged_mock, MOI.ConstraintName(), rsocs[1], "rsoc")

        s = """
        variables: t, u, x, y
        rsoc: [t, u, x, y] in MathOptInterface.RotatedSecondOrderCone(4)
        c: t + u <= 2.0
        maxobjective: x + y
        """
        model = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model, s)
        MOIU.test_models_equal(bridged_mock, model, var_names, ["rsoc", "c"])
    end

    @testset "Delete" begin
        tuxy = MOI.get(bridged_mock, MOI.ListOfVariableIndices())

        message = string("Cannot delete variable as it is constrained with other",
                         " variables in a `MOI.VectorOfVariables`.")
        for i in eachindex(tuxy)
            err = MOI.DeleteNotAllowed(tuxy[i], message)
            @test_throws err MOI.delete(bridged_mock, tuxy[i])
        end

        test_delete_bridged_variables(bridged_mock, tuxy, MOI.RotatedSecondOrderCone, 4, (
            (MOI.VectorOfVariables, MOI.SecondOrderCone, 0),
        ))
    end
end


text:
using Test
using Knet.Train20: gaussian, xavier_uniform, xavier_normal, bilinear

@testset "distributions" begin
    @test isa(gaussian(10),Array)
    #deprecated @test isa(xavier(10,10),Array)
    @test isa(xavier_uniform(10,10),Array)
    @test isa(xavier_normal(10,10),Array)
    @test isa(bilinear(Float32,2,2,128,128),Array)
end

nothing


text:
function uv_access_writable(path, mode = 2)
    local ret
    req = Libc.malloc(Base._sizeof_uv_fs)
	try
		if path == ""
			path = "."
		end
        ret = ccall(:uv_fs_access, Int32, (Ptr{Nothing}, Ptr{Nothing}, Cstring, Int64, Ptr{Nothing}), Base.eventloop(), req, path, mode, C_NULL)
        ccall(:uv_fs_req_cleanup, Nothing, (Ptr{Nothing},), req)
    finally
        Libc.free(req)
    end
    return ret == 0
end


text:
# generate all possible combinations
const VALS = Int8.(1:10)

function gen_all_solved_comb()
    all_comb = Iterators.product([VALS for i in 1:4]...) |> collect ;
    move_combine_compact!.(all_comb)
end

const MOVE_LOOKUP = gen_all_comb()

function move_lookup!(x)
    x .= MOVE_LOOKUP[x[1], x[2], x[3], x[4]]
    x
end


text:
@with_kw struct MobileNetV3 <: ModuleSpec
    stages::AbstractVector{MBConv}
    "Size of last convolution"
    k_out::Int
end

function (mbnv3::MobileNetV3)()
    head = ConvBlock(ksize = 3, k_in = 3, k_out = 16; σ = relu, stride = 2)()

    return Chain(
        head,
        [stage() for stage in mbnv3.stages]...,
        Conv((1, 1), mbnv3.stages[end].k_out => mbnv3.k_out),
        BatchNorm(mbnv3.k_out, relu)
    )

end

"""
($TYPEDEF)

Classification head for MobileNetV3 as described in
MobileNetV3-small from [Searching for MobileNetV3](https://arxiv.org/abs/1905.02244)
section 5.1 and figure 5.

$(TYPEDFIELDS)
"""
@with_kw struct MobileNetV3Head <: ModuleSpec
    "Number of classes to predict"
    n_classes
    "Number of input kernels"
    k_in
    "Number of intermediate kernels"
    k_mid
end

function (head::MobileNetV3Head)()
    return Chain(
        GlobalMeanPool(),
        Conv((1, 1), head.k_in => head.k_mid, relu),
        Conv((1, 1), head.k_mid => head.n_classes),
        flatten
    )
end

"""
    mobilenetv3_small(usedepthwise = false)

MobileNetV3-small from [Searching for MobileNetV3](https://arxiv.org/abs/1905.02244)
"""
function mobilenetv3_small(usedepthwise = false)
    cb = usedepthwise ? DepthwiseSeparable : ConvBlock
    mbvn3 = MobileNetV3([
        MBConv(ksize = 3, k_in = 16, k_exp = 16, k_out = 16,
            σ = relu, has_se = true, stride = 2, convblock = cb),
        MBConv(ksize = 3, k_in = 16, k_exp = 72, k_out = 24,
            σ = relu, has_se = false, stride = 2, convblock = cb),
        MBConv(ksize = 3, k_in = 24, k_exp = 88, k_out = 24,
            σ = relu, has_se = false, stride = 1, convblock = cb),
        MBConv(ksize = 5, k_in = 24, k_exp = 96, k_out = 40,
            σ = relu, has_se = true, stride = 2, convblock = cb),
        MBConv(ksize = 5, k_in = 40, k_exp = 240, k_out = 40,
            σ = relu, has_se = true, stride = 1, convblock = cb),
        MBConv(ksize = 5, k_in = 40, k_exp = 240, k_out = 40,
            σ = relu, has_se = true, stride = 1, convblock = cb),
        MBConv(ksize = 5, k_in = 40, k_exp = 120, k_out = 48,
            σ = relu, has_se = true, stride = 1, convblock = cb),
        MBConv(ksize = 5, k_in = 48, k_exp = 144, k_out = 48,
            σ = relu, has_se = true, stride = 1, convblock = cb),
        MBConv(ksize = 5, k_in = 48, k_exp = 288, k_out = 96,
            σ = relu, has_se = true, stride = 2, convblock = cb),
        MBConv(ksize = 5, k_in = 96, k_exp = 576, k_out = 96,
            σ = relu, has_se = true, stride = 1, convblock = cb),
        MBConv(ksize = 5, k_in = 96, k_exp = 576, k_out = 96,
            σ = relu, has_se = true, stride = 1, convblock = cb),
    ], 576)

    return mbvn3()
end

"""
    mobilenetv3_large(usedepthwise = false)

MobileNetV3-large from [Searching for MobileNetV3](https://arxiv.org/abs/1905.02244)
"""
function mobilenetv3_large(usedepthwise = false)
    cb = usedepthwise ? DepthwiseSeparable : ConvBlock
    mbvn3 = MobileNetV3([
        MBConv(ksize = 3, k_in = 16, k_exp = 16, k_out = 16,
            σ = relu, has_se = false, stride = 1, convblock = cb),

        MBConv(ksize = 3, k_in = 16, k_exp = 64, k_out = 24,
            σ = relu, has_se = false, stride = 2, convblock = cb),
        MBConv(ksize = 3, k_in = 24, k_exp = 72, k_out = 24,
            σ = relu, has_se = false, stride = 1, convblock = cb),

        MBConv(ksize = 5, k_in = 24, k_exp = 72, k_out = 40,
            σ = relu, has_se = true, stride = 2, convblock = cb),
        MBConv(ksize = 5, k_in = 40, k_exp = 120, k_out = 40,
            σ = relu, has_se = true, stride = 1, convblock = cb),
        MBConv(ksize = 5, k_in = 40, k_exp = 120, k_out = 40,
            σ = relu, has_se = true, stride = 1, convblock = cb),

        MBConv(ksize = 3, k_in = 40, k_exp = 240, k_out = 80,
            σ = relu, has_se = false, stride = 2, convblock = cb),
        MBConv(ksize = 3, k_in = 80, k_exp = 200, k_out = 80,
            σ = relu, has_se = false, stride = 1, convblock = cb),
        MBConv(ksize = 3, k_in = 80, k_exp = 184, k_out = 80,
            σ = relu, has_se = false, stride = 1, convblock = cb),
        MBConv(ksize = 3, k_in = 80, k_exp = 184, k_out = 80,
            σ = relu, has_se = false, stride = 1, convblock = cb),
        MBConv(ksize = 3, k_in = 80, k_exp = 480, k_out = 112,
            σ = relu, has_se = true, stride = 1, convblock = cb),
        MBConv(ksize = 3, k_in = 112, k_exp = 672, k_out = 112,
            σ = relu, has_se = true, stride = 1, convblock = cb),

        MBConv(ksize = 5, k_in = 112, k_exp = 672, k_out = 160,
            σ = relu, has_se = true, stride = 2, convblock = cb),
        MBConv(ksize = 5, k_in = 160, k_exp = 960, k_out = 160,
            σ = relu, has_se = true, stride = 1, convblock = cb),
        MBConv(ksize = 5, k_in = 160, k_exp = 960, k_out = 160,
            σ = relu, has_se = true, stride = 1, convblock = cb),
    ], 960)

    return mbvn3()
end

function mobilenetv3_head_small(n_classes)
    return MobileNetV3Head(
        n_classes = n_classes,
        k_in = 576,
        k_mid = 1024
    )()
end

function mobilenetv3_head_large(n_classes)
    return MobileNetV3Head(
        n_classes = n_classes,
        k_in = 576,
        k_mid = 1024
    )()
end


text:
module RocketSubscribableTest

using Test
using Rocket

@testset "Subscribable" begin

    struct DummyType end

    struct NotImplementedSubscribable <: Subscribable{Int} end
    struct NotImplementedScheduledSubscribable <: ScheduledSubscribable{Int} end
    Rocket.getscheduler(::NotImplementedScheduledSubscribable) = AsapScheduler()

    struct ExplicitlyDefinedSimpleSubscribable end
    Rocket.as_subscribable(::Type{<:ExplicitlyDefinedSimpleSubscribable}) = SimpleSubscribableTrait{String}()

    struct ExplicitlyDefinedScheduledSubscribable end
    Rocket.as_subscribable(::Type{<:ExplicitlyDefinedScheduledSubscribable}) = ScheduledSubscribableTrait{String}()
    Rocket.getscheduler(::ExplicitlyDefinedScheduledSubscribable) = AsapScheduler()

    struct ImplementedSimpleSubscribable <: Subscribable{Int} end

    function Rocket.on_subscribe!(::ImplementedSimpleSubscribable, actor)
        next!(actor, 1)
        return voidTeardown
    end

    struct ImplementedScheduledSubscribable <: ScheduledSubscribable{Int} end
    Rocket.getscheduler(::ImplementedScheduledSubscribable) = AsapScheduler()

    function Rocket.on_subscribe!(::ImplementedScheduledSubscribable, actor, scheduler)
        next!(actor, 1, scheduler)
        return voidTeardown
    end

    @testset "as_subscribable" begin
        # Check if arbitrary dummy type has invalid subscribable type
        @test as_subscribable(DummyType) === InvalidSubscribableTrait()

        # Check if as_subscribable returns valid subscribable type for subtypes of Subscribable abstract type
        @test as_subscribable(NotImplementedSubscribable)          === SimpleSubscribableTrait{Int}()
        @test as_subscribable(NotImplementedScheduledSubscribable) === ScheduledSubscribableTrait{Int}()

        # Check if as_subscribable returns valid subscribable type for explicitly defined types
        @test as_subscribable(ExplicitlyDefinedSimpleSubscribable)    === SimpleSubscribableTrait{String}()
        @test as_subscribable(ExplicitlyDefinedScheduledSubscribable) === ScheduledSubscribableTrait{String}()
    end

    struct SimpleActor{T} <: Actor{T}
        events :: Vector{Any}

        SimpleActor{T}() where T = new(Vector{Any}())
    end

    Rocket.on_next!(actor::SimpleActor, data) = push!(actor.events, data)
    Rocket.on_error!(actor::SimpleActor, err) = push!(actor.events, "err: $err")
    Rocket.on_next!(actor::SimpleActor)       = push!(actor.events, "completed")

    @testset "subscribe!" begin
        # Check if subscribe! throws an error for not valid subscribable
        @test_throws InvalidSubscribableTraitUsageError             subscribe!(DummyType(), void(Any))
        @test_throws InvalidActorTraitUsageError                    subscribe!(NotImplementedSubscribable(), DummyType())
        @test_throws MissingOnSubscribeImplementationError          subscribe!(NotImplementedSubscribable(), void(Any))
        @test_throws MissingOnScheduledSubscribeImplementationError subscribe!(NotImplementedScheduledSubscribable(), void(Any))
        @test_throws MissingOnSubscribeImplementationError          subscribe!(ExplicitlyDefinedSimpleSubscribable(), void(Any))
        @test_throws MissingOnScheduledSubscribeImplementationError subscribe!(ExplicitlyDefinedScheduledSubscribable(), void(Any))

        # Check if subscribe! subscribes to a valid subscribable
        actor = SimpleActor{Int}()
        @test subscribe!(ImplementedSimpleSubscribable(), actor) === voidTeardown
        @test actor.events == [ 1 ]

        actor = SimpleActor{Int}()
        @test subscribe!(ImplementedScheduledSubscribable(), actor) === voidTeardown
        @test actor.events == [ 1 ]

        # Check if subscribe! throws an error if subscribable and actor data types does not match
        actor = SimpleActor{String}()
        @test_throws InconsistentActorWithSubscribableDataTypesError subscribe!(ImplementedSimpleSubscribable(), actor)
        @test actor.events == [ ]

        actor = SimpleActor{String}()
        @test_throws InconsistentActorWithSubscribableDataTypesError subscribe!(ImplementedScheduledSubscribable(), actor)
        @test actor.events == [ ]
    end

    @testset "multiple subscribe!" begin
        # Check if subscribe! throws an error for not valid subscribable
        @test_throws InvalidSubscribableTraitUsageError             subscribe!(DummyType(), void(Any))
        @test_throws InvalidActorTraitUsageError                    subscribe!(NotImplementedSubscribable(), DummyType())
        @test_throws MissingOnSubscribeImplementationError          subscribe!(NotImplementedSubscribable(), void(Any))
        @test_throws MissingOnScheduledSubscribeImplementationError subscribe!(NotImplementedScheduledSubscribable(), void(Any))
        @test_throws MissingOnSubscribeImplementationError          subscribe!(ExplicitlyDefinedSimpleSubscribable(), void(Any))
        @test_throws MissingOnScheduledSubscribeImplementationError subscribe!(ExplicitlyDefinedScheduledSubscribable(), void(Any))

        # Check if subscribe! subscribes to a valid subscribable
        actor = SimpleActor{Int}()
        @test subscribe!([
            (ImplementedSimpleSubscribable(), actor),
        ]) == [ voidTeardown ]
        @test actor.events == [ 1 ]

        actor1 = SimpleActor{Int}()
        actor2 = SimpleActor{Int}()
        @test subscribe!([
            (ImplementedSimpleSubscribable(), actor1),
            (ImplementedSimpleSubscribable(), actor2),
        ]) == [ voidTeardown, voidTeardown ]
        @test actor1.events == [ 1 ]
        @test actor2.events == [ 1 ]

        actor = SimpleActor{Int}()
        @test subscribe!((
            (ImplementedSimpleSubscribable(), actor),
        )) === ( voidTeardown, )
        @test actor.events == [ 1 ]

        actor1 = SimpleActor{Int}()
        actor2 = SimpleActor{Int}()
        @test subscribe!((
            (ImplementedSimpleSubscribable(), actor1),
            (ImplementedSimpleSubscribable(), actor2),
        )) === ( voidTeardown, voidTeardown )
        @test actor1.events == [ 1 ]
        @test actor2.events == [ 1 ]

        ##
        actor = SimpleActor{Int}()
        @test subscribe!([
            (ImplementedScheduledSubscribable(), actor),
        ]) == [ voidTeardown ]
        @test actor.events == [ 1 ]

        actor1 = SimpleActor{Int}()
        actor2 = SimpleActor{Int}()
        @test subscribe!([
            (ImplementedScheduledSubscribable(), actor1),
            (ImplementedScheduledSubscribable(), actor2),
        ]) == [ voidTeardown, voidTeardown ]
        @test actor1.events == [ 1 ]
        @test actor2.events == [ 1 ]

        actor = SimpleActor{Int}()
        @test subscribe!((
            (ImplementedScheduledSubscribable(), actor),
        )) === ( voidTeardown, )
        @test actor.events == [ 1 ]

        actor1 = SimpleActor{Int}()
        actor2 = SimpleActor{Int}()
        @test subscribe!((
            (ImplementedScheduledSubscribable(), actor1),
            (ImplementedScheduledSubscribable(), actor2),
        )) === ( voidTeardown, voidTeardown )
        @test actor1.events == [ 1 ]
        @test actor2.events == [ 1 ]
    end

    struct NotImplementedFactory <: AbstractActorFactory end

    struct ImplementedFactory <: AbstractActorFactory end

    actor = SimpleActor{Int}()

    Rocket.create_actor(::Type{Int}, factory::ImplementedFactory) = actor

    @testset "subscribe! with factory" begin
        @test_throws MissingCreateActorFactoryImplementationError subscribe!(ImplementedSimpleSubscribable(), NotImplementedFactory())
        @test_throws MissingCreateActorFactoryImplementationError subscribe!(ImplementedScheduledSubscribable(), NotImplementedFactory())

        @test subscribe!(ImplementedSimpleSubscribable(), ImplementedFactory()) === voidTeardown
        @test actor.events == [ 1 ]

        @test subscribe!(ImplementedScheduledSubscribable(), ImplementedFactory()) === voidTeardown
        @test actor.events == [ 1, 1 ]
    end

    @testset "Subscribable extract type and eltype" begin

        @test Rocket.subscribable_extract_type(NotImplementedSubscribable)             === Int
        @test Rocket.subscribable_extract_type(NotImplementedScheduledSubscribable)    === Int
        @test Rocket.subscribable_extract_type(ExplicitlyDefinedSimpleSubscribable)    === String
        @test Rocket.subscribable_extract_type(ExplicitlyDefinedScheduledSubscribable) === String
        @test Rocket.subscribable_extract_type(ImplementedSimpleSubscribable)          === Int
        @test Rocket.subscribable_extract_type(ImplementedScheduledSubscribable)       === Int

        @test eltype(NotImplementedSubscribable)             === Int
        @test eltype(NotImplementedScheduledSubscribable)    === Int
        @test eltype(ExplicitlyDefinedSimpleSubscribable)    === Any
        @test eltype(ExplicitlyDefinedScheduledSubscribable) === Any
        @test eltype(ImplementedSimpleSubscribable)          === Int
        @test eltype(ImplementedScheduledSubscribable)       === Int

        @test Rocket.subscribable_extract_type(NotImplementedSubscribable())             === Int
        @test Rocket.subscribable_extract_type(NotImplementedScheduledSubscribable())    === Int
        @test Rocket.subscribable_extract_type(ExplicitlyDefinedSimpleSubscribable())    === String
        @test Rocket.subscribable_extract_type(ExplicitlyDefinedScheduledSubscribable()) === String
        @test Rocket.subscribable_extract_type(ImplementedSimpleSubscribable())          === Int
        @test Rocket.subscribable_extract_type(ImplementedScheduledSubscribable())       === Int

        @test eltype(NotImplementedSubscribable())             === Int
        @test eltype(NotImplementedScheduledSubscribable())    === Int
        @test eltype(ExplicitlyDefinedSimpleSubscribable())    === Any
        @test eltype(ExplicitlyDefinedScheduledSubscribable()) === Any
        @test eltype(ImplementedSimpleSubscribable())          === Int
        @test eltype(ImplementedScheduledSubscribable())       === Int

        @test_throws InvalidSubscribableTraitUsageError Rocket.subscribable_extract_type(DummyType)
        @test_throws InvalidSubscribableTraitUsageError Rocket.subscribable_extract_type(DummyType())

        @test eltype(DummyType) === Any
        @test eltype(DummyType()) === Any

    end
end

end


text:
module JIT

using LLVM
import LLVM: TargetMachine

import GPUCompiler: CompilerJob
import ..Compiler

export get_trampoline

# We have one global JIT and TM
const jit = Ref{OrcJIT}()
const tm  = Ref{TargetMachine}()

get_tm() = tm[]

function __init__()
    opt_level = Base.JLOptions().opt_level
    if opt_level < 2
        optlevel = LLVM.API.LLVMCodeGenLevelNone
    elseif opt_level == 2
        optlevel = LLVM.API.LLVMCodeGenLevelDefault
    else
        optlevel = LLVM.API.LLVMCodeGenLevelAggressive
    end

    tm[] = LLVM.JITTargetMachine(optlevel=optlevel)
    LLVM.asm_verbosity!(tm[], true)

    jit[] = OrcJIT(tm[]) # takes ownership of tm
    atexit() do
        dispose(jit[])
    end
end

mutable struct CallbackContext
    job::CompilerJob
    stub::String
    compiled::Bool
end

const outstanding = IdDict{CallbackContext, Nothing}()

# Setup the lazy callback for creating a module
function callback(orc_ref::LLVM.API.LLVMOrcJITStackRef, callback_ctx::Ptr{Cvoid})
    orc = OrcJIT(orc_ref)
    cc = Base.unsafe_pointer_to_objref(callback_ctx)::CallbackContext

    @assert !cc.compiled
    job = cc.job

    thunk = Compiler._link(job, Compiler._thunk(job))
    cc.compiled = true
    delete!(outstanding, cc)

    # 4. Update the stub pointer to point to the recently compiled module
    set_stub!(orc, cc.stub, thunk.adjoint)

    # 5. Return the address of tie implementation, since we are going to call it now
    ptr = thunk.adjoint
    return UInt64(reinterpret(UInt, ptr))
end

function get_trampoline(job)
    cc = CallbackContext(job, String(gensym(:trampoline)), false)
    outstanding[cc] = nothing

    c_callback = @cfunction(callback, UInt64, (LLVM.API.LLVMOrcJITStackRef, Ptr{Cvoid}))

    orc = jit[]
    initial_addr = callback!(orc, c_callback, pointer_from_objref(cc))
    create_stub!(orc, cc.stub, initial_addr)
    return address(orc, cc.stub)
end

function resolver(name, ctx)
    name = unsafe_string(name)
    ptr = try
        ## Step 0: Should have already resolved it iff it was in the
        ##         same module
        ## Step 1: See if it's something known to the execution enging
        # TODO: Do we need to do this?
        # address(jit[], name)

        ## Step 2: Search the program symbols
        #
        # SearchForAddressOfSymbol expects an unmangled 'C' symbol name.
        # Iff we are on Darwin, strip the leading '_' off.
        @static if Sys.isapple()
            if name[1] == '_'
                name = name[2:end]
            end
        end
        LLVM.API.LLVMSearchForAddressOfSymbol(name)
        ## Step 4: Lookup in libatomic
        # TODO: Do we need to do this?
    catch ex
        @error "Enzyme: Lookup failed" jl_name exception=(ex, Base.catch_backtrace())
        C_NULL
    end
    if ptr === C_NULL
        @show name
        error("Enzyme: Symbol lookup failed. Aborting!")
    end

    return UInt64(reinterpret(UInt, ptr))
end

function add!(mod)
    return compile!(jit[], mod, @cfunction(resolver, UInt64, (Cstring, Ptr{Cvoid})))
end

function lookup(jitted_mod, name)
    return LLVM.addressin(jit[], jitted_mod, name)
end

end


text:
struct CollectCovers{MC,SC}
    mosaics::MC
    sets::SC

    cover_params::CoverParams
    enum_params::CoverEnumerationParams
    opt_params::AbstractOptimizerParams

    function CollectCovers(
        mosaics::MC, sets::SC,
        cover_params::CoverParams=CoverParams(),
        enum_params::CoverEnumerationParams=CoverEnumerationParams(),
        opt_params::AbstractOptimizerParams=QuadraticOptimizerParams()
    ) where {MC,SC}
        new{MC,SC}(mosaics, sets, cover_params, enum_params, opt_params)
    end
end

function (worker::CollectCovers)(mosaic_key, set_key; verbose::Bool=false)
    mosaic_masked = mask(worker.mosaics[mosaic_key], [worker.sets[set_key]],
                         max_overlap_logpvalue=worker.enum_params.max_set_score+worker.cover_params.sel_tax)
    covers_coll = collect(mosaic_masked, worker.cover_params, worker.enum_params, worker.opt_params,
                          verbose)
    ((mosaic_key, set_key), !isempty(covers_coll) ? covers_coll : nothing)
end

function pcollect(
    mosaics::Dict{MK}, sets::Dict{SK};
    mode=:parallel,
    cover_params::CoverParams=CoverParams(),
    enum_params::CoverEnumerationParams=CoverEnumerationParams(),
    opt_params::AbstractOptimizerParams=QuadraticOptimizerParams(),
    pids=workers(),
    verbose::Bool=false,
    kwargs...
) where {MK, SK}
    @info("Parallel OESC ($(length(mosaics)) mosaic(s) × $(length(sets)) set(s))...")
    # FIXME workaround for pmap() being unable to handle callable objects +
    # unable to efficiently handle large callable objects
    collect_covers = CollectCovers(mosaics, sets, cover_params, enum_params, opt_params)
    #raw_res = sizehint!(Vector{Tuple{MK,SK,CoverCollection}}(),
    #                    length(tasks))
    mosaicXset_keys = vec([(m, s) for m in keys(mosaics), s in keys(sets)])
    if mode == :parallel
        keyXcover_vec = @distributed (append!) for (m,s) in shuffle!(mosaicXset_keys)
            res = Any[collect_covers(m, s)]
            @info("$(myid()): done mosaic=$m × set=$s")
            return res
        end
    elseif mode == :sequential
        keyXcover_vec = map(mosaicXset_keys) do mXs
            m, s = mXs
            verbose && @info("Calculating mosaic=$m × set=$s...")
            collect_covers(m, s, verbose=verbose)
        end
    else
        throw(ArgumentError("Unknown mode $mode"))
    end
    verbose && @info("Done covers enumeration, preparing the results...")
    # filter out empty covers
    keyXcover = sizehint!(Dict{Tuple{MK,SK}, CoverCollection}(), length(keyXcover_vec))
    for (mXs_key, covers_coll) in keyXcover_vec
        if covers_coll !== nothing
            keyXcover[mXs_key] = covers_coll
        end
    end
    @info("$(length(keyXcover)) non-trivial cover(s) collected")
    return keyXcover
end

function pcollect(sets1_colls, sets2;
    problem_type::Symbol=:quadratic, mode::Symbol=:sequential,
    cover_params::CoverParams=CoverParams(),
    enum_params::CoverEnumerationParams=CoverEnumerationParams(),
    verbose::Bool=false, pids=workers(),
    kwargs...
)
    verbose && @info("Preparing mosaics...")
    #p = Progress(sum(length, values(entrez_clusters))*length(entrez_colls), 5.0,
    #              "Computing cluster coverings...")   # minimum update interval: 5 second
    sets1_mosaics = Dict(pmap(sets1_colls) do coll_kv
        verbose && @info("Preparing $(coll_kv[1]) mosaic...")
        Pair(coll_kv[1], SetMosaic(coll_kv[2]))
    end)
    return pcollect(sets1_mosaics, sets2,
                    problem_type=problem_type, mode=mode, pids=pids,
                    cover_params=cover_params, enum_params=enum_params,
                    opt_params=OptimizerParams(problem_type, Workers=pids, kwargs),
                    verbose=verbose)
end


text:
using Glimpse
const Gl = Glimpse
using Test

dio = Gl.Diorama(; backgroud = Gl.RGBAf0(0.0, 0.0, 0.0, 1.0))
# write your own tests here
finalize(dio) = (Gl.Overseer.prepare(dio); Gl.Overseer.update_systems(dio.ledger))
include_finalize(str) = (include(str); finalize(dio))
@testset "text" begin
    include_finalize("text.jl")
end
@testset "gui_text" begin
    include_finalize("gui_text.jl")
end
@testset "mouse_selection" begin
    include_finalize("mouse_selection.jl")
end
@testset "orbitals" begin
    include("orbitals.jl")
end
@testset "screensaver" begin
    include_finalize("screensaver.jl")
end
@testset "oscillating_depth_peel" begin
    include_finalize("oscillating_depth_peel.jl")
end
@testset "instanced_depth_peel" begin
    include_finalize("instanced_depth_peel.jl")
end